Line data Source code
1 : MODULE m_lhcal
2 : use m_juDFT
3 :
4 : !*********************************************************************
5 : ! determines the lattice harmonics for the given local
6 : ! symmetry defined in terms of the rotation matrices.
7 : !
8 : ! input: lmax max l for lattice harmonics
9 : ! nrot number of operations
10 : ! orth rotation matrices (cartesian coordinates)
11 : ! ( 1,1 1,2 1,3 )( x ) ( x' )
12 : ! ( 2,1 2,2 2,3 )( y ) = ( y' )
13 : ! ( 3,1 3,2 3,3 )( z ) ( z' )
14 : ! memd max. members/lattice harmonic dimensioned for
15 : ! nlhd max. number of lattice harmonics dimensioned for
16 : !
17 : ! output: nlh number of harmonic from l=0 to l=lmax
18 : ! lnu l value of harmonic
19 : ! mem number of members in lattice harmonic
20 : ! lmnu lm index of each member of lattice harmonic
21 : ! c coefficients of lattice harmonic
22 : !
23 : ! m. weinert 1989/98
24 : !*********************************************************************
25 :
26 : CONTAINS
27 372 : SUBROUTINE lhcal(
28 372 : > memd,nlhd,lmax,nrot,orth,
29 372 : < nlh,lnu,mem,lmnu,c)
30 : !DEC$ NOOPTIMIZE
31 : USE m_constants
32 : USE m_gaussp
33 : USE m_gtest
34 : USE m_ylm
35 : IMPLICIT NONE
36 :
37 : !---> Arguments
38 : INTEGER, INTENT (IN) :: memd,nlhd
39 : INTEGER, INTENT (IN) :: lmax !---> max. l to consider
40 : INTEGER, INTENT (IN) :: nrot !---> number of
41 : REAL, INTENT (IN) :: orth(3,3,nrot) !---> rotation matrices
42 :
43 : INTEGER, INTENT(OUT) :: nlh
44 : INTEGER, INTENT(OUT) :: lnu((lmax+1)**2),mem((lmax+1)**2)
45 : INTEGER, INTENT(OUT) :: lmnu(memd,(lmax+1)**2)
46 : COMPLEX, INTENT(OUT) :: c(memd,(lmax+1)**2)
47 :
48 : !---> Locals
49 : INTEGER :: j,l,lh,lm,lm0,lmmax,l2,m,mems,mp,n,nn
50 : REAL :: s
51 372 : REAL :: v(3),vr(3),ovlp(0:2*lmax)
52 372 : COMPLEX :: a(0:lmax,0:2*lmax,0:lmax)
53 372 : COMPLEX, DIMENSION((LMAX+1)**2) :: ylm,ylmr,ylms
54 :
55 : INTEGER :: ngpts ! gaussian
56 372 : REAL :: vgauss(3,(2*lmax+1)*(lmax+1 + mod(lmax+1,2))) ! integration
57 372 : REAL :: wt((2*lmax+1)*( lmax+1 + mod(lmax+1,2) )) ! points
58 :
59 : REAL, PARAMETER :: del = 1.e-5 ! parameter for machine roundoff, etc
60 : !
61 : !---> generate gaussian points and test
62 : !
63 372 : !$OMP MASTER
64 372 : ngpts = (2*lmax+1)*( lmax+1 + mod(lmax+1,2) )
65 : CALL gaussp(
66 : > lmax,
67 372 : < vgauss,wt)
68 : CALL gtest(
69 372 : > lmax,ngpts,vgauss,wt)
70 :
71 : !---> initialize
72 372 : lmmax = (lmax+1)**2
73 739152 : a = cmplx(0.0,0.0)
74 : !
75 : !===> loop over gaussian integration points
76 : !
77 73356 : DO nn = 1, ngpts
78 291936 : v(:) = vgauss(:,nn)
79 : CALL ylm4(
80 : > lmax,v,
81 72984 : < ylm)
82 7686128 : ylms(1:lmmax) = ylm(1:lmmax)
83 : !---> apply rotations
84 957368 : DO n =2, nrot
85 11496992 : vr = matmul( orth(:,:,n), v )
86 : CALL ylm4(
87 : > lmax,vr,
88 884384 : < ylmr)
89 100030616 : ylms(:) = ylms(:) + ylmr(:)
90 : ENDDO
91 7686128 : ylms = ylms/nrot
92 : !---> obtain coefficients
93 808580 : DO l = 0, lmax
94 735224 : lm0 = l*(l+1)+1
95 4909408 : DO mp = 0,l
96 : a(mp,0,l) = a(mp,0,l) +
97 4909408 : + wt(nn)*conjg(ylm(lm0+mp))*real(ylms(lm0))
98 : ENDDO
99 4247168 : DO m = 1, l
100 3438960 : lm = lm0+m
101 30935464 : DO mp = 0, l
102 : a(mp,m,l) = a(mp,m,l) +
103 26761280 : + wt(nn)*conjg(ylm(lm0+mp))* real(ylms(lm))
104 : a(mp,m+l,l) = a(mp,m+l,l) +
105 30200240 : + wt(nn)*conjg(ylm(lm0+mp))*aimag(ylms(lm))
106 : ENDDO
107 : ENDDO
108 : ENDDO
109 : ENDDO
110 : !
111 : !===> orthonormalize the projections for each l (many maybe zero)
112 : !
113 372 : nlh = 0
114 3928 : DO l = 0, lmax
115 3556 : lm0 = l*(l+1)+1
116 3556 : l2 = 2*l
117 38828 : DO m = 0, l2
118 : !---> calculate overlaps with previous harmonics for this l
119 251484 : DO j = 0, m-1
120 216584 : s = real( conjg(a(0,j,l)) * a(0,m,l))
121 1806576 : DO mp=1,l
122 1806576 : s = s + 2*real( conjg(a(mp,j,l)) * a(mp,m,l) )
123 : ENDDO
124 251484 : ovlp(j) = s
125 : ENDDO
126 : !---> Gram-Schmidt
127 251484 : DO j = 0,m-1
128 2058060 : a(0:l,m,l) = a(0:l,m,l) - ovlp(j)*a(0:l,j,l)
129 : ENDDO
130 : !---> normalize
131 34900 : s = real( conjg(a(0,m,l)) * a(0,m,l))
132 251484 : DO mp = 1, l
133 251484 : s = s + 2*real( conjg(a(mp,m,l)) * a(mp,m,l))
134 : ENDDO
135 38456 : IF (s.GT.del) THEN
136 14368 : s = sqrt(s)
137 113016 : a(0:l,m,l) = a(0:l,m,l)/s
138 : !---> store lattice harmonic
139 14368 : mems = 0
140 14368 : nlh = nlh + 1
141 14368 : IF (nlh>nlhd+1) CALL juDFT_error("nlhd too small",
142 0 : + calledby="lhcal")
143 14368 : lnu(nlh) = l
144 14368 : IF (abs(a(0,m,l)).GT.del) THEN
145 2608 : mems = mems+1
146 2608 : IF (mems>memd) CALL juDFT_error("memd too small"
147 0 : + ,calledby ="lhcal")
148 2608 : c(mems,nlh) = a(0,m,l)
149 2608 : lmnu(mems,nlh) = lm0
150 : ENDIF
151 98648 : DO mp=1,l
152 98648 : IF( abs(a(mp,m,l)).GT.del) THEN
153 12596 : mems = mems + 1
154 12596 : IF(mems>memd) CALL juDFT_error("memd too small"
155 0 : + ,calledby ="lhcal")
156 12596 : c(mems,nlh) = a(mp,m,l)
157 12596 : lmnu(mems,nlh) = lm0 + mp
158 12596 : mems = mems + 1
159 12596 : IF (mems>memd) CALL juDFT_error("memd too small"
160 0 : + ,calledby ="lhcal")
161 12596 : c(mems,nlh) = ((-1.)**mp)*conjg(c(mems-1,nlh))
162 12596 : lmnu(mems,nlh) = lm0 - mp
163 : ENDIF
164 : ENDDO
165 14368 : mem(nlh) = mems
166 : ELSE
167 173368 : a(0:l,m,l) = cmplx(0.0,0.0)
168 : ENDIF
169 : ENDDO ! m = 0, l2
170 : ENDDO ! l = 0, lmax
171 : !
172 : !===> test of lattice harmonics using an arbitary point
173 : !
174 372 : v(1) = sqrt(2.0)
175 372 : v(2) = sqrt(5.0)
176 372 : v(3) = sqrt(17.0)
177 : !---> generate lattice harmonic for this point
178 : CALL ylm4(
179 : > lmax,v,
180 372 : < ylm)
181 14740 : DO lh = 1, nlh
182 14368 : ylms(lh) = cmplx(0.0,0.0)
183 42540 : DO m = 1, mem(lh)
184 42168 : ylms(lh) = ylms(lh) + c(m,lh)*ylm(lmnu(m,lh))
185 : ENDDO
186 : ENDDO
187 : !---> rotate point and generate lattice harmonic
188 4664 : DO n = 2, nrot
189 55796 : vr = matmul( orth(:,:,n), v )
190 : CALL ylm4(
191 : > lmax,vr,
192 4292 : < ylm)
193 39640 : DO lh = 1, nlh
194 34976 : ylmr(lh) = cmplx(0.0,0.0)
195 117464 : DO m = 1, mem(lh)
196 117464 : ylmr(lh) = ylmr(lh) + c(m,lh)*ylm(lmnu(m,lh))
197 : ENDDO
198 39268 : IF ( abs(ylms(lh)-ylmr(lh)).GT.del ) THEN
199 0 : WRITE (oUnit,'(/," error for operation",i3)') n
200 0 : WRITE (oUnit,'( " lattice harmonic ",i3)') lh
201 0 : WRITE (oUnit,'(4f12.6)') ylms(lh),ylmr(lh)
202 0 : CALL juDFT_error("k_lv(Rr)",calledby="lhcal")
203 : ENDIF
204 : ENDDO
205 : ENDDO
206 : !$OMP END MASTER
207 372 : RETURN
208 : END SUBROUTINE lhcal
209 : END MODULE m_lhcal
|
