-
Notifications
You must be signed in to change notification settings - Fork 20
/
basis_mod.f90
2174 lines (1907 loc) · 73.4 KB
/
basis_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! fortran_dialect=elf
!-----------------------------------------------------------------------
! PSCF - Polymer Self-Consistent Field Theory
! Copyright (2002-2016) Regents of the University of Minnesota
! contact: David Morse, [email protected]
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation. A copy of this license is included in
! the LICENSE file in the top-level PSCF directory.
!-----------------------------------------------------------------------
!****m* scf/basis_mod
! MODULE
! basis_mod - create symmetry-adapted basis functions
!
! PURPOSE
! Define data structures and procedures involving reciprocal
! lattice vectors, stars, and symmetry-adapted basis functions
!
! AUTHOR
!
! David Morse (2002)
! Chris Tyler (2002-2003) various corrections and additions
!
! SOURCE
!----------------------------------------------------------------------
module basis_mod
use const_mod ! Defines constant long, and variable dim=1,2, or 3
use io_mod, only : output
use string_mod, only : int_string
use version_mod, only : version_type, output_version
use unit_cell_mod, only : G_basis, make_G_basis, output_unit_cell
use group_mod ! Defines types and operations for crystal groups
use space_groups_mod, only : space_groups
use grid_mod ! Utilities for fft grids
implicit none
private
! Public procedures
public:: make_basis, make_dGsq ! needed by scf
public:: make_waves, make_stars, release_basis
public:: f_basis
public:: reorder_basis
public:: scattering_intensity
public:: output_waves
! Public variables
public:: group
public:: N_star
public:: N_wave, wave, star_of_wave, coeff
public:: G_max, which_wave
public:: star_begin, star_end, star_count
public:: star_cancel, star_invert, basis_sign
public:: wave_of_star
public:: sort
public:: valid_wave
! Declaration of public variables
type(group_type) :: group ! space group
integer :: N_wave ! # of reciprocal lattice G vectors
integer :: N_star ! # of stars
integer :: G_max(3) ! maximum of values G(1), G(2), G(3)
integer :: wave(:,:) ! (dim,N_wave) list of G vectors
real(long) :: Gsq(:) ! (N_wave) list of values of |G|^{2}
complex(long):: coeff(:) ! (N_wave) coefficient of wave in star
integer :: star_of_wave(:) ! (N_wave) index of star containing wave
integer :: which_wave(:,:,:)! (-G_max(1):G_max(1), ..., ... )
! indices of waves listed by components
integer :: star_begin(:) ! (N_star) index of first wave in star
integer :: star_end(:) ! (N_star) index of last wave in star
integer :: star_count(:) ! (N_star) # of G vectors in star
integer :: wave_of_star(:,:)! (dim,N_star) components of one wave from
! the star, serves as a label
logical :: star_cancel(:) ! (N_star) see wave_format
integer :: star_invert(:) ! (N_star) see wave_format
integer :: basis_sign(:) ! (N_star) see wave_format
!***
allocatable :: wave, Gsq, coeff, star_of_wave, which_wave
allocatable :: star_begin, star_end, star_count, wave_of_star
allocatable :: star_invert, star_cancel, basis_sign
! private parameters
integer, parameter :: max_star = 48 ! Max. # of G vectors in star
integer, parameter :: max_list =110592! Max # of G vecs in a list
! private variables
logical :: grid ! true if PS method is used
! Documentation headers for public global variables
!----------------------------------------------------------------
!****v basis_mod/N_wave
! VARIABLE
! integer N_wave = # of G vectors (reciprocal lattice vectors)
!*** ------------------------------------------------------------
!****v basis_mod/N_star
! VARIABLE
! integer N_star = # of stars
!*** ------------------------------------------------------------
!****v basis_mod/G_max
! VARIABLE
! integer G_max(3) = Max values of components of integer waves
!*** ------------------------------------------------------------
!****v basis_mod/wave
! VARIABLE
! integer wave(:,:) - allocatable, dimension(dim,N_wave)
! wave(1:dim,i) = integer wavevectors # i, Bravais basis
!*** ------------------------------------------------------------
!****v basis_mod/Gsq
! VARIABLE
! real(long) Gsq(:) - allocatable, dimension(N_wave)
! Gsq(i) = value of |G|^{2} for wavevector i
!*** ------------------------------------------------------------
!****v basis_mod/coeff
! VARIABLE
! complex(long) coeff(:) - allocatable, dimension(N_wave)
! coeff(i) = complex coefficients of wavevector i
! in symmetry adapted function for star
!*** ------------------------------------------------------------
!****v basis_mod/star_of_wave
! VARIABLE
! integer star_of_wave(:) - allocatable, dimension(N_wave)
! star_of_wave(i) = index of star containing wave i
!*** ------------------------------------------------------------
!****v basis_mod/which_wave
! VARIABLE
! integer which_wave(G1,G2,G3) = index of wavevector G=(G1,G2,G3)
!
! If wave(i) = (G1,G2,G3), then which_wave(G1,G2,G3) = i
! allocatable, dimension( -G_max(1):G_max(1) , ... , ... )
!*** ------------------------------------------------------------
!****v basis_mod/star_begin
! VARIABLE
! integer star_begin(:) - allocatable, dimension(N_star)
! star_begin(i) = index of first G in star i
!*** ------------------------------------------------------------
!****v basis_mod/star_end
! VARIABLE
! integer star_end(:) - allocatable, dimension(N_star)
! star_end(i) = index of last G in star i
!*** ------------------------------------------------------------
!****v basis_mod/star_count
! VARIABLE
! integer star_count(:) - allocatable, dimension(N_star)
! star_count(i) = # of G vectors in star i
!*** ------------------------------------------------------------
!****v basis_mod/star_invert
! VARIABLE
! integer star_invert(:) - allocatable, dimension(N_star)
! star_invert(i) = invert flag for star i, values (0,1,-1)
! See discussion of wave_format
!*** ------------------------------------------------------------
!****v basis_mod/basis_sign
! VARIABLE
! integer basis_sign(:) - allocatable, dimension(N_star)
! basis_sign(i) = sign of basis function i under inversion
! +1 -> even under inversion
! -1 -> odd under inversion
! See discussion of wave_format
!*** ------------------------------------------------------------
!****v basis_mod/wave_of_star
! VARIABLE
! integer wave_of_star(:,:) - dimension(dim, N_star)
! wave_of_star(:,i) = G vector from star i, serves as a label
! See discussion of wave_format
!*** ------------------------------------------------------------
!****v basis_mod/star_cancel
! VARIABLE
! logical star_cancel(:) - allocatable, dimension(N_star)
! star_cancel(i) = true if star i is cancelled
!*** ------------------------------------------------------------
!-----------------------------------------------------------------------
!****c* basis_mod/wave_format
! COMMENT
!
! i) G vectors are listed in the array wave in non-decreasing order
! of eigen (for some G_basis passed to make_waves), and (after
! output from make_stars) are grouped by stars, each of which
! forms contiguous block of vectors.
!
! ii) Cancelled waves and stars are included in the lists of waves
! and stars only if logical variable keep_cancel is passed to
! make_waves with a value .true. If keep_cancel = true,
! cancelled waves are asigned values coeff()=0.0 and cancelled
! starts are assigned star_cancel = true. If keep_cancel = false,
! then star_cancel = false for all stars.
!
! iii) Within each star, vectors are listed in "descending" order,
! if the vectors in wave are read as dim=2, or 3 digit numbers,
! with more significant digits on the left. For example the
! [111] star of a cubic structure is listed in the order:
!
! 1 1 1 (first)
! 1 1 -1
! 1 -1 1
! 1 -1 -1
! -1 1 1
! -1 1 -1
! -1 -1 1
! -1 -1 -1 (last)
!
! iv) If star number i is closed under inversion, so that the
! negative -G of each wave G in the star is also an element
! of that star, then we assign star_invert(i) = 0. For
! space groups that contain inversion symmetry (i.e., for
! centro-symmetric groups), all stars must be closed under
! inversion.
!
! v) For non-centrosymmetric space groups, pairs of stars may be
! related by inversion, so that the -G of each vector in star
! number i is found not in star i but in a second star. Such
! pairs of stars that are related by inversion are always listed
! consecutively. The first star in such a pair is assigned
! star_invert(i) = 1, and the second has star_invert(i+1) = -1.
!
! vi) The combination of conventions (iii) and (v) implies that the
! inverse of the first vector in a star with star_invert = 1 is
! the last vector in the next star, which is the partner star
! with star_invert = -1. More generally, the inverse of the ith
! member of the first star in such a pair is the ith to last
! vector of the second star in the pair.
!
! vi) wave_of_star is first wave in star for star_invert=0 or 1, and
! last wave in star for star_invert=-1
!
! vii) The absolute magnitudes of the coefficients in array coeff()
! for (un-cancelled) stars are chosen so that each stars yields
! an orthonormal but generally complex function, i.e., values
! of coeff() have absolute magnitude 1.0/sqrt(N_in_star).
!
! viii) The phases of the coefficients in the array coeff() are chosen
! such that the coefficient of wave_of_star is positive and real.
! This is sufficient to guarantee that a vector -G has value of
! coeff that is the complex conjugate of the coeff for G, whether
! G and -G are in the same star star (for star_invert = 0) or in
! two consecutive stars (for star_invert = +-1). It also
! guarantees that the star basis function generated by the first
! star in an inversion-related pair is the complex conjugate of
! that generated by the second star in the pair. As a result, if
! the first star in a pair with star_invert = +-1 is multiplied
! by a coefficient and the second by the complex conjugate of
! this coefficient, the resulting sum of plane waves will yield
! a real basis function.
!
! A star with star_invert = 0 will yield a normalized real basis
! function, while two linearly independent real basis functions
! may be constructed out of a consecutive pair of stars with
! star_invert=1 and -1, by multiplying both by 1/sqrt(2) to
! obtain a cosine-like function, and by multiplying the first
! star by i/sqrt(2) and the second by -i/sqrt(2), to obtain a
! sine-like real function. The number of real basis functions
! will thus be equal to the number of stars N_star.
!
! ix) The stars are used to construct N_star real basis functions
! f_1(r), ... , f_{N_star}(r), which are defined and indexed
! as follows:
!
! For stars that are closed under inversion, with star_invert=0,
! and a real star basis function phi_j(r), we create a single
! real basis function that is is equal to the normalized star
! basis function
!
! f_j(r) = phi_j(r) ,
!
! and that is assigned a basis function index equal to star
! index of the star.
!
! For pairs of stars that are related by inversion, with star
! indices j and j+1, with associated star basis functions
! phi_{j}(r) and phi_{j+1}(r) we construct one cosine-like
! basis function, with basis i, given by the sum
!
! f_j(r) = [ phi_{j}(r) + phi_{j+1}(r) ]/sqrt(2) ,
!
! which is even under inversion, and a second sine-like basis,
! with a basis function index j+1, given by
!
! f_{j+1}(r) = [ phi_{j}(r) - phi_{j+1}(r) ]/sqrt(-2)
!
! which is odd under inversion. Note that sqrt(-2) is a pure
! imaginary number.
!
! xi) For non-centrosymmetric groups, basis functions may be
! either either even or odd under inversion. These two cases
! are distinguished by the value of variable basis_sign:
!
! basis_sign(i) = 1 implies f_{i}(r) = + f_{i}(-r)
! basis_sign(i) = -1 implies f_{i}(r) = - f_{i}(-r)
!
! In a centro-symmetric group basis_sign(i) = 1 for all
! basis functions.
!
! In a non-centrosymmetric group, a basis function that is
! closed under inversion, so that star_invert(i) = 0, may
! be either even or odd or under inversion, and may thus
! have either value of basis_sign(i). This is why we need
! to distinguish basis_sign from star_invert.
!
! In a non-centrosymmetric group, pairs of sine-like and
! cosine-like basis functions that are made constructed
! from a pairs of inversion-related stars are indexed such
! that:
!
! If star_invert(i) = 1, then basis_sign(i) = 1
! If star_invert(i) = -1, then basis_sign(i) = 1
!
!***
!-------------------------------------------------------------------
contains
!-------------------------------------------------------------------
!****p basis_mod/release_basis
!
! SUBROUTINE
! release_basis()
! PURPOSE
! release the memory allocated by constructing basis information
! SOURCE
!--------------------------------------------------------------------
subroutine release_basis
!***
integer :: error ! error index
if(allocated(wave)) deallocate(wave,stat=error)
if (error /= 0) stop "wave deallocation error!"
if(allocated(Gsq)) deallocate(Gsq,stat=error)
if (error /= 0) stop "Gsq deallocation error!"
if(allocated(coeff)) deallocate(coeff,stat=error)
if (error /= 0) stop "coeff deallocation error!"
if(allocated(star_of_wave)) deallocate(star_of_wave,stat=error)
if (error /= 0) stop "star_of_wave deallocation error!"
if(allocated(which_wave)) deallocate(which_wave,stat=error)
if (error /= 0) stop "which_wave deallocation error!"
if(allocated(star_begin)) deallocate(star_begin,stat=error)
if (error /= 0) stop "star_begin deallocation error!"
if(allocated(star_end)) deallocate(star_end,stat=error)
if (error /= 0) stop "star_end deallocation error!"
if(allocated(star_count)) deallocate(star_count,stat=error)
if (error /= 0) stop "star_count deallocation error!"
if(allocated(wave_of_star)) deallocate(wave_of_star,stat=error)
if (error /= 0) stop "wave_of_star deallocation error!"
if(allocated(star_invert)) deallocate(star_invert,stat=error)
if (error /= 0) stop "star_invert deallocation error!"
if(allocated(star_cancel)) deallocate(star_cancel,stat=error)
if (error /= 0) stop "star_cancel deallocation error!"
if(allocated(basis_sign)) deallocate(basis_sign,stat=error)
if (error /= 0) stop "basis_sign deallocation error!"
end subroutine release_basis
!====================================================================
!-------------------------------------------------------------------
!****p basis_mod/make_basis
! SUBROUTINE
!
! make_basis(R_basis,G_basis,group_name,Gabs_max)
!
! PURPOSE
!
! Generates information about basis functions needed by scf code:
! i) Generates reciprocal lattice basis vectors G_basis, using
! input R_basis (Bravais basis vectors)
! ii) Reads group generators from file named by the character
! string group_name, and generates remainder of the group
! iii) Calls make_waves and make_stars to generate all waves
! and stars with |G| < Gabs_max, excluding cancelled stars
!
! SOURCE
!--------------------------------------------------------------------
subroutine make_basis( &
R_basis, &! real lattice basis vectors
G_basis, &! reciprocal lattice basis vectors
group_name, &! file for generators of space group
N_grids, &! maximum value of |G|
grid_flag &! true if PS method used
)
logical, intent(IN) :: grid_flag
real(long) , intent(IN) :: R_basis(:,:) ! real lattice basis
real(long) , intent(OUT) :: G_basis(:,:) ! Reciprocal basis
character(*), intent(INOUT) :: group_name ! file for group
integer , intent(IN) :: N_grids(3) ! # of grid points
!***
! Local variables
integer :: i,j,k,l
logical :: keep_cancel
!logical, SAVE :: first_time = .true.
character(len = 3) :: Nchar
keep_cancel = .false.
grid=grid_flag
! Make reciprocal lattice basis
call make_G_basis(R_basis,G_basis)
! Select space group by name in space_groups
call space_groups(group_name,group)
call make_group(group,R_basis,G_basis)
! Make waves
call make_waves( &
G_basis, &! (dim,dim), basis for reciprocal lattice
R_basis, &! (dim,dim), basis for Bravais lattice
N_grids, &! (3), # of grid points
keep_cancel &! if true, keep cancelled waves
)
!group waves into stars related by space group symmetries
call make_stars
if (grid) call make_ksq(G_basis) ! true if PS method used
end subroutine make_basis
!====================================================================
!--------------------------------------------------------------------
!****ip basis_mod/make_waves
! SUBROUTINE
! make_waves(G_basis,R_basis,N_grids,keep_cancel)
! PURPOSE
! Subroutine fills arrays wave, Gsq with lists of all reciprocal
! lattice vectors in FFT grid, sort in ascending order of Gsq.
! ARGUMENTS
! Inputs
! real(long) G_basis(:,:) - basis for reciprocal lattice
! real(long) R_basis(:,:) - basis for Bravais lattice
! integer N_grids(3) - # of grid points in x, y, z
! logical keep_cancel - if true, keep cancelled waves
! SOURCE
!--------------------------------------------------------------------
subroutine make_waves( &
G_basis, &! (dim,dim), basis for reciprocal lattice
R_basis, &! (dim,dim), basis for Bravais lattice
N_grids, &! (3), # of grid points
keep_cancel &! if true, keep cancelled waves
)
real(long), intent(IN) :: G_basis(:,:), R_basis(:,:)
integer, intent(IN) :: N_grids(3)
logical :: keep_cancel
!***
! local variables
real(long) :: Gabs_max
real(long) :: Gsq_new, Gsq_max, G_vec(3), twopi, a_abs
integer :: i1, i2, i3, i_vec(3), i, j, k
logical :: cancel
! Temporary arrays to store index and values of wave, and Gsq
integer, allocatable, dimension(:) :: index
integer, allocatable, dimension(:,:) :: wave_temp
real(long), allocatable, dimension(:) :: Gsq_temp
Gabs_max = max_Gabs(G_basis) ! grid_mod/max_Gabs
! Calculate G_max
Gsq_max = Gabs_max*Gabs_max
twopi = 4.0_long*acos(0.0_long)
G_max = 0
do i=1, dim
a_abs = sqrt( R_basis(i,:) .dot. R_basis(i,:) )
G_max(i) = int( Gabs_max * a_abs / twopi ) + 1
enddo
! Calculate number of plane waves
N_wave = 0
if ( .not. grid ) then ! use spectral method
do i1= -G_max(1), G_max(1)
i_vec(1) = i1
if ( dim == 1) then
G_vec= i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
if ( Gsq_new <= Gsq_max ) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
N_wave = N_wave + 1
endif
endif
else
do i2= -G_max(2), G_max(2)
i_vec(2) = i2
if (dim == 2) then
G_vec= i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
if ( Gsq_new <= Gsq_max ) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
N_wave = N_wave + 1
endif
endif
else
do i3= -G_max(3), G_max(3)
i_vec(3) = i3
G_vec = i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
! write(6,*) i_vec(1),i_vec(2),i_vec(3), Gsq_new
if (Gsq_new <= Gsq_max) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
N_wave = N_wave + 1
endif
endif
enddo
endif
enddo
endif
enddo
else ! use shifted FFT grid to generate wave
do i1=0,N_grids(1)-1
i_vec(1)=i1
do i2=0,N_grids(2)-1
i_vec(2)=i2
do i3=0,N_grids(3)-1
i_vec(3)=i3
i_vec = G_to_bz(i_vec)
cancel= vector_cancel(i_vec)
if( (.not.cancel) .or. keep_cancel ) N_wave=N_wave+1
enddo
enddo
enddo
end if ! basis if-then
! Allocate module arrays
if (allocated(wave)) deallocate(wave)
allocate(wave(dim,N_wave), stat = j)
if (j.ne.0) stop 'make_basis allocate wave error'
if (allocated(Gsq)) deallocate(Gsq)
allocate(Gsq(N_wave), stat = j)
if (j.ne.0) stop 'make_basis allocate N_wave error'
if (allocated(star_of_wave)) deallocate(star_of_wave)
allocate(star_of_wave(N_wave), stat = j)
if (j.ne.0) stop 'Make_basis allocate star_of_wave error'
if (allocated(which_wave)) deallocate(which_wave)
allocate(which_wave&
(-G_max(1):G_max(1),-G_max(2):G_max(2),-G_max(3):G_max(3)),stat=j)
if (j.ne.0) stop 'Make_basis allocate which_wave error'
if (allocated(coeff)) deallocate(coeff)
allocate(coeff(N_wave),stat=j)
if (j.ne.0) stop 'Make_basis coef allocate error)'
! Allocate temporary arrays
allocate(wave_temp(dim,N_wave))
allocate(Gsq_temp(N_wave))
allocate(index(N_wave))
! Make wave, and Gsq
if ( .not. grid ) then ! use spectral method
j = 0
i_vec = 0
do i1= -G_max(1), G_max(1)
i_vec(1) = i1
if (dim == 1) then
G_vec= i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
if (Gsq_new <= Gsq_max) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
j = j + 1
Gsq(j) = Gsq_new
!do k=1, dim
! wave(k,j) = i_vec(k)
!enddo
call assign_ivec(wave(:,j),i_vec)
endif
endif
else
do i2= -G_max(2), G_max(2)
i_vec(2) = i2
if (dim == 2) then
G_vec= i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
if (Gsq_new <= Gsq_max) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
j = j + 1
Gsq(j) = Gsq_new
!do k=1, dim
! wave(k,j) = i_vec(k)
!enddo
call assign_ivec(wave(:,j),i_vec)
endif
endif
else
do i3= -G_max(3), G_max(3)
i_vec(3) = i3
G_vec= i_vec .dot. G_basis
Gsq_new = G_vec .dot. G_vec
if (Gsq_new <= Gsq_max) then
cancel = vector_cancel(i_vec)
if ((.not.cancel).or.keep_cancel) then
j = j + 1
Gsq(j) = Gsq_new
!do k=1, dim
! wave(k,j) = i_vec(k)
!enddo
call assign_ivec(wave(:,j),i_vec)
endif
endif
enddo
endif
enddo
endif
enddo
else ! use shifted fft grid to generate wave
j = 0
do i1=0,N_grids(1)-1
i_vec(1)=i1
do i2=0,N_grids(2)-1
i_vec(2)=i2
do i3=0,N_grids(3)-1
i_vec(3)=i3
i_vec = G_to_bz(i_vec)
cancel= vector_cancel(i_vec)
if( (.not.cancel) .or. keep_cancel ) then
j = j + 1
Gsq(j) = norm(i_vec, G_basis)
call assign_ivec(wave(:,j),i_vec)
endif
enddo
enddo
enddo
end if ! grid if--then
! Create ordered index
do i=1, N_wave
index(i)=i
enddo
! Sort Gsq in asending order, and reorder index accordingly
call sort(N_wave,Gsq,index)
! Use reordered index to re-order wave
do i=1, N_wave
wave_temp(:,i) = wave(:,index(i))
enddo
wave = wave_temp
!Deallocate temporary arrays
deallocate(wave_temp)
deallocate(Gsq_temp)
deallocate(index)
end subroutine make_waves
!===================================================================
!-------------------------------------------------------------------
!****ip basis_mod/make_stars
! SUBROUTINE
! make_stars
! PURPOSE
! Given output of make_waves, the subroutine divides the list of
! G vectors first into "lists" of G vectors with equal values of
! Gsq and then into "stars" of vectors that are related by the
! symmetry operations of the group.
!
! On input:
! The arrays wave(dim,N_wave) and Gsq(N_wave) contain lists
! of all waves and corresponding values of |G|^{2}, ordered by
! nondecreasing values of Gsq, but not yet grouped into stars.
!
! On output:
! 1) Gsq, wave are reordered so that members of each star are
! listed consecutively
! 2) star_of_wave(i) contains index of star for wave i
! 3) coeff(i) contains complex coefficient for wave i
! 4) wave(which_wave(G(1),G(2),G(3))) = G for any wave in list
! 5) star_begin(i) and star_end(i) contain wave indices for first
! and last wave in star i, and star_count(i) = # waves in star
! 6) star_invert(i) contains an integer flag giving information
! about effect of inversion upon star i (see below)
! 7) star_cancel(i) = .true. if star i is cancelled.
! (Relevant only if make_waves left cancelled waves in lists)
!-------------------------------------------------------------------
! Stars related by inversion, and the value of star_invert:
!
! 1) Stars that are closed under inversion (i.e., that contain
! the inverse -G of every vector G in the star) have
!
! star_invert(i) = 0
!
! For centrosymmetric groups (i.e., groups that include
! inversion symmetry) all stars will be closed under inversion
!
! 2) Pairs of stars that are related by inversion are list
! consecutively, with:
!
! star_invert(i) = 1 for for the 1st star in the pair
! star_invert(i+1) = -1 for the 2nd
!
! Such inversion-related pairs of stars can exist only for
! non-centrosymmetric groups
!
! SOURCE
!-------------------------------------------------------------------
subroutine make_stars
!***
! Parameter
real(long), parameter :: epsilon = 1.0E-7_long
! Local variables
integer :: list(3,max_list), star(3,max_star), G1(3), G2(3)
real(long) :: star_phase(max_star)
real(long) :: Gsq_max, twopi
integer :: list_index(max_list), star_index(max_star)
integer :: first_list, last_list, list_N, star_N
integer :: i, j, k, l, i_index, i_star, root, invert_flag
integer :: i_first, i_last
complex(long) :: c_norm
logical :: new_list, cancel
! Arrays to store index and temporary values of wave, and Gsq
!integer :: index(N_wave)
!integer :: wave_temp(dim,N_wave)
!real(long) :: Gsq_temp(N_wave)
integer :: index(:)
integer :: wave_temp(:,:)
real(long) :: Gsq_temp(:)
allocatable :: index, wave_temp, Gsq_temp
integer :: info
real(long) :: c1, c2
allocate(index(N_wave), stat=info)
if ( info /= 0 ) then
write(6,*) "index(N_wave) allocation error: N_wave = ", N_wave
stop
end if
allocate(wave_temp(dim,N_wave), stat=info)
if ( info /= 0 ) then
write(6,*) "wave_temp(dim,N_wave) allocation error: N_wave = ", N_wave
stop
end if
allocate(Gsq_temp(N_wave), stat=info)
if ( info /= 0 ) then
write(6,*) "Gsq_temp(N_wave) allocation error: N_wave = ", N_wave
stop
end if
twopi = 4.0_long*acos(0.0_long)
which_wave = 0
! Loop over Gsq identifying places where value changes
i_index = 0
i_star = 0
first_list = 1
invert_flag = 0
Gsq_max = 0.0_long
!-----------------------------------------------------------------!
! This loop creates star_of_wave, which_wave, and index
! The array index is subsequently used to reorder wave
!-----------------------------------------------------------------!
! Begin loop over wavevectors
do i=2, N_wave+1
!--------------------------------------------------------------!
! Check for end of a list, due to change in Gsq or if i=N_wave !
! A "list" is set of vectors of equal magnitude |G| !
!--------------------------------------------------------------!
if ( i > N_wave ) then
new_list = .false.
if( first_list == N_wave ) then
new_list = .true.
last_list = N_wave
endif
else
new_list = .false.
if (Gsq(i) > Gsq_max + epsilon) then
Gsq_max = Gsq(i)
last_list = i - 1
new_list = .true.
else if (i == N_wave) then
last_list = N_wave
new_list = .true.
else if (Gsq(i) < Gsq_max - epsilon) then
write(6,*) 'Error in make_stars: unordered Gsq array'
stop
endif
end if
! Begin processing completed list
if (new_list) then
list_N = last_list - first_list + 1
!Extract list
if ( list_N > max_list) then
write(6,*) 'Error: list_N > max_list in make_stars'
write(6,*) last_list, first_list, max_list, N_wave
write(6,'(3i5)') list(:,first_list-1)
write(6,'(3i5)') list(:,first_list:last_list)
stop
endif
list = 0
do j=first_list, last_list
k = j - first_list + 1
!do l=1, dim
! list(l,k) = wave(l,j)
!enddo
call assign_ivec(list(:,k),wave(:,j))
list_index(k) = j
enddo
! Sort vectors in list
call G_sort(list_N,list,list_index)
!-------------------------------------------------------
! Begin loop over stars within list:
! Pass list to routine get_star, which outputs a "star"
! constructed from vector # root in input list, and
! remaining "list" containing vectors from input list
! that are not in this star
! Add the resulting star to global arrays.
! Repeat if the remaining "list" is not empty.
!-------------------------------------------------------
root = 1
10 call get_star(root, &
list,list_index,list_N, &
star,star_index,star_N,star_phase)
! Check that star is not empty
if (star_N == 0) then
write(6,*) ' star_N <= 0 in make_star'
stop
endif
! Update i_star (index of new star)
i_star = i_star + 1
! Sort vectors in star and in remaining list
call G_sort(star_N,star,star_index,reorder_real = star_phase )
call G_sort(list_N,list,list_index)
! Output contents of star and remaining list (diagnostic)
!write(6,*)
!write(6,*) 'star ='
!do j=1, star_N
! write(6,FMT='(2I5,3X,3I5)') j, star_index(j), star(:,j)
!enddo
!write(6,*)
!write(6,*) 'remaining list='
!do j=1, list_N
! write(6,FMT='(2I5,3X,3I5)') j, list_index(j), list(:,j)
!enddo
! Check for cancellation of first vector in star
G1 = star(:,1)
cancel = vector_cancel(G1)
! Check that entire star is cancelled or not cancelled
do j=1, star_N
if (cancel .NEQV. vector_cancel(star(:,j))) then
write(6,*) 'Inconsistent results for cancellation:'
write(6,*) 'star(:,1), cancel=',G1,cancel
cancel = vector_cancel(star(:,j))
write(6,*) 'j, star(:,j), cancel=',star(:,j),cancel
stop
endif
enddo
!---------------------------------------------------------
! Calculate normalizing denominator for un-cancelled stars:
! Coefficient of last wave in star is positive real if
! invert_flag == 1 for previous star, so that the
! current star is the 2nd star in an inversion pair
! Otherwise, star_invert = 0 or 1 for current star, and
! the coefficient of first wave is positive real
! In either case abs(c_norm) = sqrt(star_N)
!---------------------------------------------------------
if (.not.cancel) then
if (invert_flag == 1) then
c_norm = exp((0.0_long,1.0_long)*star_phase(star_N))
else
c_norm = exp((0.0_long,1.0_long)*star_phase(1))
endif
c_norm = c_norm * sqrt(real(star_N))
endif
!------------------------------------------------------
! Add star to global arrays that are indexed by wave:
! Append star_index to index,
! Add values to star_of_wave and which_wave
! Calculate coeff for star
!------------------------------------------------------
do j=1, star_N
! Update i_index (global index of new wave)
i_index = i_index + 1
index(i_index) = star_index(j)
star_of_wave(i_index) = i_star
G1 = star(:,j)
which_wave(G1(1),G1(2),G1(3)) = i_index
if (cancel) then ! coeff = zero
coeff(i_index) = (0.0_long,0.0_long)
else
coeff(i_index) = &
exp( (0.0_long,1.0_long)*star_phase(j) )/c_norm
endif
end do
!------------------------------------------------------
! Set invert_flag for this star & root for next star:
!
! If star is closed under inversion (invert_flag = 0),
! or is second star of pair related by inversion
! (invert_flag = -1), set root of next star to the
! first vector in remaining "list" : root = 1
!
! If star is not closed under inversion, and is 1st star
! in pair related by inversion (invert_flag=1), set
! root of next star to be the inverse of first vector
! in the current star, so the next star will be related
! to the current one by inversion.
!------------------------------------------------------
G1 = -1.0_long*star(:,1)
if (grid) G1 = G_to_bz(G1)
!write(6,*) '-G1=',G1
if (in_list(G1,star,star_N) /= 0) then
if (invert_flag == 1) then
! Expect second star in pair, shouldn't be closed'
write(6,*) &
'Error: Expect invert_flag = -1, but -G1 is in star'
stop
else
invert_flag = 0
root = 1
!write(6,*) 'Parity flag =0, -G1 in star'
endif
else if (invert_flag == 1) then
invert_flag = -1
root = 1
!write(6,*) 'Parity flag = -1'
else
invert_flag = 1
root = in_list(G1,list,list_N)
if (root == 0) then
write(6,*) 'Error: -G1 in neither in star nor list'
stop
else
!write(6,*) 'invert_flag = 1, -G1 in remaining list'
endif
endif
!---------------------------------------------------------
! Upon exit of above section:
! invert_flag = 0 -> this star closed under inversion
! invert_flag = 1 -> this star is 1st in pair related
! by inversion
! invert_flag = -1 -> this star is 2nd in pair related
! by inversion
!---------------------------------------------------------