Line data Source code
1 : MODULE m_gradYlm
2 : IMPLICIT NONE
3 :
4 : CONTAINS
5 :
6 : ! Derivative routine from Spex.
7 : ! Calculates derivative of a function with scalar argument lying on a muffin-tin mesh
8 196620 : subroutine Derivative(f, itype, atoms, df)
9 :
10 : use m_types
11 :
12 : type(t_atoms), intent(in) :: atoms
13 :
14 : integer, intent(in) :: itype
15 : real, intent(in) :: f(atoms%jri(itype))
16 : real, intent(out) :: df(atoms%jri(itype))
17 : real :: h, r, d21, d32, d43, d31, d42, d41, df1, df2, s
18 : real :: y0, y1, y2
19 : integer :: i, n
20 :
21 196620 : n = atoms%jri(itype)
22 196620 : h = atoms%dx(itype)
23 196620 : r = atoms%rmsh(1, itype)
24 :
25 : ! use Lagrange interpolation of 3rd order (and averaging) for points 3 to n
26 196620 : d21 = r * (exp(h)-1) ; d32 = d21 * exp(h) ; d43 = d32 * exp(h)
27 196620 : d31 = d21 + d32 ; d42 = d32 + d43
28 196620 : d41 = d31 + d43
29 : df(1) = d31*d41 / (d21*d32*d42) * f(2) + ( -1.0/d21 - 1.0/d31 - 1.0/d41) * f(1)&
30 196620 : & - d21*d41 / (d31*d32*d43) * f(3) + d21*d31 / (d41*d42*d43) * f(4)
31 : df(2) = - d32*d42 / (d21*d31*d41) * f(1) + ( 1.0/d21 - 1.0/d32 - 1.0/d42) * f(2)&
32 196620 : & + d21*d42 / (d31*d32*d43) * f(3) - d21*d32 / (d41*d42*d43) * f(4)
33 : df1 = d32*d43 / (d21*d31*d41) * f(1) - d31*d43 / (d21*d32*d42) * f(2) +&
34 196620 : & ( 1.0/d31 + 1.0/d32 - 1.0/d43 ) * f(3) + d31*d32 / (d41*d42*d43) * f(4)
35 138338700 : do i = 3, n - 2
36 138142080 : d21 = d32 ; d32 = d43 ; d43 = d43 * exp(h)
37 138142080 : d31 = d42 ; d42 = d42 * exp(h)
38 138142080 : d41 = d41 * exp(h)
39 : df2 = - d32*d42 / (d21*d31*d41) * f(i-1) + ( 1.0/d21 - 1.0/d32 - 1.0/d42) * f(i) + &
40 138142080 : & d21*d42 / (d31*d32*d43) * f(i+1) - d21*d32 / (d41*d42*d43) * f(i+2)
41 138142080 : df(i) = ( df1 + df2 ) / 2
42 : df1 = d32*d43 / (d21*d31*d41) * f(i-1) - d31*d43 / (d21*d32*d42) * f(i) +&
43 138338700 : & ( 1.0/d31 + 1.0/d32 - 1.0/d43 ) * f(i+1) + d31*d32 / (d41*d42*d43) * f(i+2)
44 : enddo
45 196620 : df(n-1) = df1
46 : df(n) = - d42*d43 / (d21*d31*d41) * f(n-3) + d41*d43 / (d21*d32*d42) * f(n-2) -&
47 196620 : & d41*d42 / (d31*d32*d43) * f(n-1) + ( 1.0/d41 + 1.0/d42 + 1.0/d43 ) * f(n)
48 : ! for first two points use Lagrange interpolation of second order for log(f(i))
49 : ! or, as a fall-back, Lagrange interpolation with the conditions f(1), f(2), f(3), f'(3).
50 196620 : s = sign(1.0,f(1))
51 578913 : if(sign(1.0,f(2)) /= s .or. sign(1.0,f(3)) /= s .or. any(f(:3) == 0)) then
52 70897 : d21 = r * (exp(h)-1)
53 70897 : d32 = d21 * exp(h)
54 70897 : d31 = d21 + d32
55 70897 : s = df(3) / (d31*d32) - f(1) / (d21*d31**2) + f(2) / (d21*d32**2) - f(3) / (d31**2*d32) - f(3) / (d31*d32**2)
56 70897 : df(1) = - (d21+d31) / (d21*d31) * f(1) + d31 / (d21*d32) * f(2) - d21 / (d31*d32) * f(3) + d21*d31 * s
57 :
58 70897 : df(2) = - (d21-d32) / (d21*d32) * f(2) - d32 / (d21*d31) * f(1) + d21 / (d31*d32) * f(3) - d21*d32 * s
59 : else
60 125723 : y0 = log(abs(f(1)))
61 125723 : y1 = log(abs(f(2)))
62 125723 : y2 = log(abs(f(3)))
63 125723 : df(1) = ( - 3*y0/2 + 2*y1 - y2/2 ) * f(1) / (h*r)
64 125723 : df(2) = (y2-y0)/2 * f(2) / (h*r*exp(h))
65 : endif
66 196620 : end subroutine Derivative
67 :
68 194 : SUBROUTINE gradYlm(atoms, r2FshMt, r2GrFshMt)
69 : ! Based on work for juphon by C. Gerhorst.
70 : !---------------------------------------------------------------------------------------------------------------------------------
71 : !> @author
72 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
73 : !>
74 : !> @brief
75 : !> Calculates the spherical harmonic expansion coefficients of the muffin-tin gradient applied to an arbitrary function multiplied
76 : !> by $r^2$. The resulting gradient expansion coefficients are multiplied by a factor $r^2$.
77 : !>
78 : !> @note
79 : !> The ingoing function is assumed to be multiplied with $r^2$. The outgoing resulting function is also multiplied by $r^2$.
80 : !>
81 : !> @param[in] atoms : Contains atoms-related quantities; definition of its members in types.F90 file.
82 : !> @param[in] r2FshMt : Lattice harmonic coefficients of muffin-tin quantity multiplied by a factor of r**2.
83 : !> @param[out] r2GrFshMt : Spherical harmonic coefficients of muffin-tin quantity's gradient multiplied by a factor of r**2
84 : !---------------------------------------------------------------------------------------------------------------------------------
85 :
86 : use m_constants, only : fpi_const, ImagUnit
87 : use m_gaunt, only : gaunt1
88 : use m_types
89 :
90 :
91 : ! Type parameter
92 : ! ***************
93 : type(t_atoms), intent(in) :: atoms
94 :
95 : ! Array parameters
96 : ! ****************
97 : complex, allocatable, intent(in) :: r2FshMt(:, :, :)
98 : complex, allocatable, intent(out) :: r2GrFshMt(:, :, :, :)
99 :
100 :
101 : ! Local Scalar Variables
102 : ! **********************
103 : ! pfac : Prefactor
104 : ! tGaunt : Gaunt coefficient
105 : ! itype : Loop index for atom types
106 : ! ieqat : Loop index for equivalent atoms
107 : ! iatom : Loop index for all atoms
108 : ! imesh : Loop index for radial mesh point
109 : ! mqn_m : Magnetic quantum number m
110 : ! oqn_l : Orbital quantum number l
111 : ! mqn_mpp : Magnetic quantum number double primed to index the natural coordinates
112 : ! lmInput : Collective index for orbital and magnetic quantum number used for input function
113 : ! lmOutput : Collective index for orbital and magnetic quantum number used for output function
114 : real :: pfac
115 : real :: tGaunt
116 : integer :: itype
117 : integer :: ieqat
118 : integer :: iatom
119 : integer :: imesh
120 : integer :: mqn_m
121 : integer :: oqn_l
122 : integer :: mqn_mpp
123 : integer :: lmOutput
124 : integer :: lmInput
125 :
126 : ! Local Array Variables
127 : ! *********************
128 : ! rDerFlhMtre : Real part of the radial derrivative applied to the incoming fuction coefficients
129 : ! rDerFlhMtim : Imaginary part of the radial derrivative applied to the incoming fuction coefficients
130 : ! rDerFlhMt : Radial derrivative of the incoming fuction
131 : ! r2GrFshMtNat : Expansion coefficients of the muffin-tin gradient applied to the incoming function. The coefficients are given
132 : ! in natural coordinates and multiplied by $r^2$
133 194 : real, allocatable :: r2FshMtre(:)
134 194 : real, allocatable :: r2FshMtim(:)
135 194 : real, allocatable :: rDerFshMtre(:)
136 194 : real, allocatable :: rDerFshMtim(:)
137 194 : real, allocatable :: rDerJunk(:)
138 194 : complex, allocatable :: rDerFshMt(:)
139 194 : complex, allocatable :: r2GrFshMtNat(:, :, :, :)
140 : !Matrix syntax idea from http://stackoverflow.com/questions/3708307/how-to-initialize-two-dimensional-arrays-in-fortran
141 : complex, dimension(3, 3) :: Tmatrix !no parameter anymore since the init below is not OK for this
142 : Tmatrix = transpose(reshape([ &
143 : & cmplx(1 / sqrt(2.), 0), cmplx(0, 0), cmplx(-1 / sqrt(2.), 0), &
144 : & -ImagUnit / sqrt(2.), cmplx(0, 0), -ImagUnit / sqrt(2.), &
145 : & cmplx(0, 0), cmplx(1, 0), cmplx(0, 0) &
146 194 : & ], [3, 3] )) ! see my notes
147 :
148 :
149 : ! Initialization of additionaly required arrays.
150 1164 : allocate( r2GrFshMt(atoms%jmtd, ( atoms%lmaxd + 2 )**2, atoms%nat, 3) )
151 776 : allocate( r2GrFshMtNat(atoms%jmtd, ( atoms%lmaxd + 2 )**2, atoms%nat, 3) )
152 1746 : allocate( r2FshMtre(atoms%jmtd), r2FshMtim(atoms%jmtd), rDerFshMtre(atoms%jmtd), rDerFshMtim(atoms%jmtd), rDerJunk(atoms%jmtd), rDerFshMt(atoms%jmtd) )
153 :
154 135658 : r2FshMtre(:) = 0.
155 135658 : r2FshMtim(:) = 0.
156 135658 : rDerFshMtre(:) = 0.
157 135658 : rDerFshMtim(:) = 0.
158 135658 : rDerFshMt(:) = 0.
159 47975654 : r2GrFshMt = cmplx(0., 0.)
160 47975654 : r2GrFshMtNat = cmplx(0., 0.)
161 :
162 : pfac = sqrt( fpi_const / 3 )
163 776 : do mqn_mpp = -1, 1
164 582 : iatom = 0
165 1454 : do itype = 1, atoms%ntype
166 1938 : do ieqat = 1, atoms%neq(itype)
167 678 : iatom = iatom + 1
168 7458 : do oqn_l = 0, atoms%lmax(itype)
169 61698 : do mqn_m = -oqn_l, oqn_l
170 :
171 : ! l + 1 block
172 54918 : if ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 ) then
173 54918 : lmOutput = ( oqn_l + 1 ) * ( oqn_l + 2 ) + 1 + mqn_m - mqn_mpp
174 54918 : lmInput = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
175 38859102 : r2FshMtre(:) = 0.
176 38859102 : r2FshMtim(:) = 0.
177 38859102 : rDerFshMtre(:) = 0.
178 38859102 : rDerFshMtim(:) = 0.
179 : ! This is also a trade of between storage and performance, because derivative is called redundantly, maybe store it?
180 38859102 : r2FshMtre(:)=real(r2FshMt(:, lmInput, itype))
181 38859102 : r2FshMtim(:)=aimag(r2FshMt(:, lmInput, itype))
182 54918 : call Derivative( r2FshMtre, itype, atoms, rDerFshMtre )
183 54918 : call Derivative( r2FshMtim, itype, atoms, rDerFshMtim )
184 :
185 54918 : tGaunt = Gaunt1( oqn_l + 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd + 1 )
186 38859102 : do imesh = 1, atoms%jri(itype)
187 38804184 : rDerFshMt(imesh) = cmplx(rDerFshMtre(imesh), rDerFshMtim(imesh))
188 : r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) = r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) + pfac * &
189 : & (-1)**mqn_mpp * tGaunt &
190 38859102 : & * (rDerFshMt(imesh) - ((oqn_l + 2)* r2FshMt(imesh, lmInput, iatom) / atoms%rmsh(imesh, itype)))
191 : end do ! imesh
192 : end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 )
193 :
194 : ! l - 1 block
195 61020 : if ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 ) then
196 43392 : lmInput = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
197 43392 : lmOutput = (oqn_l - 1) * oqn_l + 1 + mqn_m - mqn_mpp
198 30703488 : r2FshMtre(:) = 0.
199 30703488 : r2FshMtim(:) = 0.
200 30703488 : rDerFshMtre(:) = 0.
201 30703488 : rDerFshMtim(:) = 0.
202 : ! This is also a trade of between storage and performance, because derivative is called redundantly, maybe store it?
203 30703488 : r2FshMtre(:)=real(r2FshMt(:, lmInput, itype))
204 30703488 : r2FshMtim(:)=aimag(r2FshMt(:, lmInput, itype))
205 43392 : call Derivative( r2FshMtre, itype, atoms, rDerFshMtre )
206 43392 : call Derivative( r2FshMtim, itype, atoms, rDerFshMtim )
207 :
208 43392 : tGaunt = Gaunt1( oqn_l - 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd + 1 )
209 30703488 : do imesh = 1, atoms%jri(itype)
210 30660096 : rDerFshMt(imesh) = cmplx(rDerFshMtre(imesh), rDerFshMtim(imesh))
211 : r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) = r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) + pfac * &
212 : & (-1)**mqn_mpp * tGaunt * ( rDerFshMt(imesh) + ( (oqn_l - 1) * r2FshMt(imesh, lmInput, iatom)&
213 30703488 : & / atoms%rmsh(imesh, itype) ) )
214 : enddo ! imesh
215 : end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 )
216 : end do ! mqn_m
217 : end do ! oqn_l
218 : end do ! ieqat
219 : end do ! itype
220 : end do ! mqn_mpp
221 :
222 : ! Conversion from natural to cartesian coordinates
223 194 : iatom = 0
224 420 : do itype = 1, atoms%ntype
225 646 : do ieqat = 1, atoms%neq(itype)
226 226 : iatom = iatom + 1
227 2712 : do oqn_l = 0, atoms%lmax(itype) + 1
228 25086 : do mqn_m = -oqn_l, oqn_l
229 22600 : lmOutput = oqn_l * (oqn_l + 1) + 1 + mqn_m
230 15993660 : do imesh = 1, atoms%jri(itype)
231 271492200 : r2GrFshMt(imesh, lmOutput, iatom, 1:3) = matmul( Tmatrix(1:3, 1:3), r2GrFshMtNat(imesh, lmOutput, iatom, 1:3) )
232 : end do ! imesh
233 : end do ! mqn_m
234 : end do ! oqn_l
235 : end do ! ieqat
236 : end do ! itype
237 :
238 194 : END SUBROUTINE gradYlm
239 :
240 64 : SUBROUTINE divYlm(gradMtx, gradMty, gradMtz, divMt)
241 : COMPLEX, INTENT(IN) :: gradMtx(:,:,:,:), gradMty(:,:,:,:), gradMtz(:,:,:,:) !r,lm,n,x
242 : COMPLEX, INTENT(OUT) :: divMt(:,:,:)
243 :
244 4235952 : divMt(:,:,:)=gradMtx(:,:,:,1)+gradMty(:,:,:,2)+gradMtz(:,:,:,3)
245 :
246 64 : END SUBROUTINE divYlm
247 :
248 0 : SUBROUTINE rotYlm(gradMtx, gradMty, gradMtz, rotMt)
249 : COMPLEX, INTENT(IN) :: gradMtx(:,:,:,:), gradMty(:,:,:,:), gradMtz(:,:,:,:) !r,lm,n,x
250 : COMPLEX, INTENT(OUT) :: rotMt(:,:,:,:)
251 :
252 0 : rotMt(:,:,:,1)=gradMtz(:,:,:,2)-gradMty(:,:,:,3)
253 0 : rotMt(:,:,:,2)=gradMtx(:,:,:,3)-gradMtz(:,:,:,1)
254 0 : rotMt(:,:,:,3)=gradMty(:,:,:,1)-gradMtx(:,:,:,2)
255 :
256 0 : END SUBROUTINE rotYlm
257 : END MODULE m_gradYlm
|
