Line data Source code
1 : module m_dfpt_eii2
2 :
3 : use m_types
4 :
5 : implicit none
6 :
7 : contains
8 :
9 :
10 : ! should be correct, has been reviewed
11 0 : subroutine GenPsDens2ndOrd(atoms, cell, ngpqdp, G0index, gpqdp, qpt, psDens2ndOrd, testMode)
12 :
13 : use m_sphbes
14 :
15 : implicit none
16 :
17 : ! Type parameter
18 : type(t_atoms), intent(in) :: atoms
19 : type(t_cell), intent(in) :: cell
20 :
21 : ! Scalar parameter
22 : integer, intent(in) :: ngpqdp
23 : logical, intent(in) :: testMode
24 : integer, intent(out):: G0index
25 :
26 : ! Array parameter
27 : integer, intent(in) :: gpqdp(:, :)
28 : real, intent(in) :: qpt(:)
29 : complex, allocatable, intent(out):: psDens2ndOrd(:, :, :, :)
30 :
31 : ! Scalar variable
32 : integer :: itype
33 : integer :: ii
34 : integer :: jj
35 : integer :: iG
36 : integer :: iatom
37 : integer :: iatomTemp
38 : integer :: ieqat
39 :
40 : ! Array variable
41 0 : real, allocatable :: sbes(:, :, :)
42 0 : real, allocatable :: psDensMat(:, :, :) ! Matrix part of pseudo density
43 0 : real, allocatable :: Gpqext(:, :)
44 0 : real, allocatable :: Gpq(:, :)
45 0 : complex, allocatable :: phaseFactor(:, :)
46 0 : real, allocatable :: GpqextRmt(:, :)
47 0 : real, allocatable :: prefactor(:)
48 :
49 : ! Initializiation of local arrays
50 0 : allocate(psDensMat(ngpqdp, 3, 3))
51 0 : allocate(psDens2ndOrd(ngpqdp, 3, atoms%nat, 3))
52 0 : allocate(Gpqext(3, ngpqdp))
53 0 : allocate(Gpq(3, ngpqdp))
54 0 : allocate( sbes( 0 : MAXVAL(atoms%ncv) + 1, ngpqdp, atoms%ntype) )
55 0 : allocate( phaseFactor(ngpqdp, atoms%nat) )
56 0 : allocate( GpqextRmt(ngpqdp, atoms%ntype) )
57 0 : allocate( prefactor(atoms%ntype) )
58 :
59 0 : sbes(:, :, :) = 0.
60 0 : GpqextRmt(:, :) = 0.
61 0 : psDensMat(:, :, :) = 0.
62 0 : psDens2ndOrd(:, :, :, :) = cmplx(0.0, 0.0)
63 0 : Gpqext(:, :) = 0.
64 0 : Gpq(:, :) = 0.
65 0 : G0index = -1
66 0 : prefactor(:) = 1.
67 0 : phaseFactor(:, :) = cmplx(0., 0.)
68 :
69 : ! Precalculate the matrix-like part of the pseudodensity
70 : ! todo If we optimize for the dynamical matrix in the end, we need one non-linear run through memory. As we have to evaluate every matrix
71 : ! element of the dynamical matrix seperately, it makes sense to shift the indices of the 3x3 matrix after iG to avoid redundant
72 : ! operations later.
73 0 : do iG = 1, ngpqdp
74 : ! If denominator gets zero skip G-vector, that happens only for q = 0, basically.
75 0 : if ( norm2( gpqdp(1:3, iG) + qpt(1:3) ) <= 1e-9 ) then
76 0 : G0index = iG
77 0 : cycle
78 : end if
79 :
80 : ! (G+q)(G+q)^T and ((G+q)(G+q)^T - (G+q)^2/3)
81 0 : Gpq(1:3, iG) = real(gpqdp(1:3, iG) + qpt(1:3))
82 0 : Gpqext(1:3, iG) = matmul(Gpq(1:3, iG), cell%bmat(1:3, 1:3))
83 0 : if ( testMode ) then
84 0 : psDensMat(iG, 1:3, 1:3) = outerProduct(Gpqext(1:3, iG), Gpqext(1:3, iG))
85 : else
86 0 : psDensMat(iG, 1:3, 1:3) = outerProduct(Gpqext(1:3, iG), Gpqext(1:3, iG)) - (norm2(Gpqext(1:3, iG))**2 * id3x3(1:3, 1:3) / 3.)
87 : end if
88 : end do ! iG
89 :
90 : iatom = 0
91 0 : do itype = 1, atoms%ntype
92 :
93 : ! Calculate double factorial (2 N + 7)!! in 7.78 (PhD thesis Klüppelberg). Note, atoms%ncv(itype) is already the value taken from
94 : ! Table 1 in the Weinert paper. For the orbital quantum number l = 2, this leads to
95 : ! (2 * atoms%ncv(itype) + 3)!! = (2 * N + 2 * l + 3)!! = (2 * N + 7)!!.
96 : ! The right hand side of the recent equation is consistent with 7.78 (PhD thesis Klüppelberg), where N is the Weinert parameter.
97 0 : do ii = 1, 2 * atoms%ncv(itype) + 3, 2
98 0 : prefactor(itype) = prefactor(itype) * ii
99 : end do ! ii
100 :
101 : ! Complete prefactor.
102 0 : prefactor(itype) = prefactor(itype) * atoms%zatom(itype) / cell%omtil
103 :
104 0 : iatomTemp = iatom
105 0 : do iG = 1, ngpqdp
106 :
107 : ! No pseudo density contribution for G + q = 0
108 0 : if (iG == G0index) cycle
109 :
110 0 : GpqextRmt(iG, itype) = norm2(Gpqext(1:3, iG)) * atoms%rmt(itype)
111 : ! sbes is initialized within sphbes with first parameter of sphbes
112 0 : call sphbes(atoms%ncv(itype) + 1, GpqextRmt(iG, itype), sbes(0:atoms%ncv(itype) + 1, iG, itype))
113 0 : iatom = iatomTemp
114 0 : do ieqat = 1, atoms%neq(itype)
115 0 : iatom = iatom + 1
116 0 : phaseFactor(iG, iatom) = exp(-tpi_const * ImagUnit * dot_product(Gpq(1:3, iG), atoms%taual(1:3, iatom)))
117 : end do ! ieqat
118 : end do ! iG
119 : end do ! itype
120 :
121 0 : do jj = 1, 3
122 : iatom = 0
123 0 : do itype = 1, atoms%ntype
124 0 : do ieqat = 1, atoms%neq(itype)
125 0 : iatom = iatom + 1
126 0 : do ii = 1, 3
127 0 : do iG = 1, ngpqdp
128 0 : if (iG == G0index) cycle
129 : ! Note that atoms%ncv = N + l, so for l = 2, we need only an exponent or the spherical Bessel function, respectively,
130 : ! of ncv + 1.
131 : psDens2ndOrd(iG, ii, iatom, jj) = prefactor(itype) / GpqextRmt(iG, itype)**(atoms%ncv(itype) + 1) &
132 0 : & * sbes(atoms%ncv(itype) + 1, iG, itype) * psDensMat(iG, ii, jj) * phaseFactor(iG, iatom)
133 : end do ! iG
134 : end do ! jj
135 : end do ! ieqat
136 : end do ! itype
137 : end do ! jj
138 :
139 0 : end subroutine GenPsDens2ndOrd
140 :
141 0 : subroutine CalcIIEnerg2MatElem( atoms, cell, qpt, ngpqdp, gpqdp, E2ndOrdII )
142 :
143 : use m_sphbes
144 :
145 : implicit none
146 :
147 : ! Type parameter
148 : type(t_atoms), intent(in) :: atoms
149 : type(t_cell), intent(in) :: cell
150 :
151 : ! Scalar parameter
152 : integer, intent(in) :: ngpqdp
153 :
154 : ! Array parameter
155 : integer, intent(in) :: gpqdp(:, :)
156 : real, intent(in) :: qpt(:)
157 : complex, intent(out):: E2ndOrdII(:, :)
158 :
159 : ! Scalar variable
160 : integer :: G0index
161 : integer :: iG
162 : integer :: iAdir
163 : integer :: iBdir
164 : integer :: oqn_l
165 : integer :: mqn_m
166 : integer :: lm
167 : integer :: t
168 : logical :: testMode
169 : integer :: iAtype
170 : integer :: iBtype
171 : integer :: iAatom
172 : integer :: iBatom
173 : integer :: iAeqat
174 : integer :: iBeqat
175 :
176 :
177 : ! Array variables
178 0 : complex, allocatable :: psDens2ndOrd(:, :, :, :)
179 0 : complex, allocatable :: expAlpha(:)
180 0 : real, allocatable :: Gpqext(:, :)
181 0 : real, allocatable :: sbes(:, :)
182 :
183 0 : allocate( expAlpha(ngpqdp) )
184 0 : allocate( Gpqext(3, ngpqdp))
185 0 : allocate( sbes(0:atoms%lmaxd, ngpqdp) )
186 :
187 0 : E2ndOrdII(:, :) = cmplx(0.0, 0.0)
188 0 : expAlpha(:) = cmplx(0.0, 0.0)
189 0 : Gpqext(:, :) = 0.0
190 0 : sbes(:, :) = 0.0
191 :
192 : ! Generates pseudodensity as it is given in 7.95 Aaron Phd thesis, or for q = 0 in 7.78 Aaron Phd thesis
193 : ! We leave the test mode on, leading to the fact that the trace is not subtracted. Therefore, we have a non-vanishing diagonal.
194 0 : testMode = .true.
195 0 : call GenPsDens2ndOrd(atoms, cell, ngpqdp, G0index, gpqdp, qpt, psDens2ndOrd, testMode)
196 :
197 : ! Leave it here so it needs not be calculated 3N x 3N times.
198 0 : do iG = 1, ngpqdp
199 0 : Gpqext(1:3, iG) = matmul(real(gpqdp(1:3, iG) + qpt(1:3)), cell%bmat(1:3, 1:3))
200 : end do ! iG
201 :
202 0 : iAatom = 0
203 0 : do iAtype = 1, atoms%ntype
204 0 : sbes(:, :) = 0.0
205 0 : do iG = 1, ngpqdp
206 : ! If we precalcalculate the scalar factor we have 7 multiplications less
207 : !todo we can also only calculate it to 0 not to lmaxd for performance reasons
208 0 : call Sphbes( atoms%lmaxd, norm2(Gpqext(1:3, iG)) * atoms%rmt(iAtype), sbes(0:atoms%lmaxd, iG) )
209 : end do
210 0 : do iAeqat = 1, atoms%neq(iAtype)
211 0 : iAatom = iAatom + 1
212 0 : expAlpha(:) = cmplx(0.0, 0.0)
213 0 : do iG = 1, ngpqdp
214 0 : expAlpha(iG) = exp(tpi_const * ImagUnit * dot_product(gpqdp(1:3, iG) + qpt(1:3), atoms%taual(1:3, iAatom)))
215 : end do
216 0 : do iAdir = 1, 3
217 : iBatom = 0
218 0 : do iBtype = 1, atoms%ntype
219 0 : do iBeqat = 1, atoms%neq(iBtype)
220 0 : iBatom = iBatom + 1
221 0 : if (iBatom /= iAatom) then
222 0 : do iBdir = 1, 3
223 : ! Add contribution 7.99
224 0 : do iG = 1, ngpqdp
225 0 : if (iG == G0index) cycle
226 : E2ndOrdII( iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1) ) = &
227 : & E2ndOrdII( iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1) ) &
228 : + atoms%zatom(iAtype) * fpi_const / &
229 0 : & norm2(Gpqext(1:3, iG))**2 * psDens2ndOrd(iG, iBdir, iBatom, iAdir) * expAlpha(iG)
230 : end do ! iG
231 : end do ! iBdir
232 : else ! iBatom = iAatom
233 : oqn_l = 0
234 : mqn_m = 0
235 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
236 0 : do iBdir = 1, 3
237 0 : do iG = 1, ngpqdp
238 0 : if (iG == G0index) cycle
239 : E2ndOrdII( iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1) ) = &
240 : E2ndOrdII( iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1) ) + atoms%zatom(iAtype) * fpi_const * &
241 0 : & psDens2ndOrd(iG, iBdir, iBatom, iAdir) * sbes(0, iG) / norm2(Gpqext(:, iG))**2 * expAlpha(iG)
242 : end do
243 : ! Constant term is switched off because it is subtracted away anyway in Equation 7.89, as q-independent
244 0 : if (.false.) then
245 : do t = -1, 1
246 : E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1)) = &
247 : & E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom -1)) &
248 : - atoms%zatom(iAtype) * atoms%zatom(iBtype) / atoms%rmt(iBtype)**3 * &
249 : & ( 3 / fpi_const * c_im(iBdir, t + 2) * c_im(iAdir, 2 - t) * (-1)**t )
250 : end do ! t
251 : end if ! Constant term switched off
252 : end do ! iBdir
253 : end if ! iBatom = iAatom ?
254 : end do ! iBeqat
255 : end do ! iBtype
256 :
257 : ! Constant term is switched off because it is subtracted away anyway in Equation 7.89 because it is not dependent on q
258 0 : if (.false.) then
259 : ! Factor 3 because it is within the t-sum
260 : E2ndOrdII(iAdir + 3 * (iAatom - 1), iAdir + 3 * (iAatom -1)) = E2ndOrdII(iAdir, iAdir) + 3 * atoms%zatom(iAtype) * atoms%zatom(iAtype) / atoms%rmt(iAtype)**3
261 : end if ! Constant term switched off
262 :
263 : end do ! iAdir
264 : end do ! iAeqat
265 : end do ! iAtype
266 :
267 0 : end subroutine CalcIIEnerg2MatElem
268 :
269 0 : SUBROUTINE getConstTerm(atoms, cell, constTerm)
270 : TYPE(t_atoms), INTENT(IN) :: atoms
271 : TYPE(t_cell), INTENT(IN) :: cell
272 :
273 : REAL, INTENT(OUT) :: constTerm(:, :)
274 :
275 : INTEGER :: iAlpha, iBeta
276 :
277 : REAL :: tauExtAlpha(3), tauExtBeta(3), vecR(3)
278 :
279 0 : constTerm = 0.0
280 :
281 0 : DO iBeta = 1, atoms%nat
282 0 : tauExtBeta = MATMUL(cell%amat, atoms%taual(:, iBeta))
283 0 : DO iAlpha = 1, atoms%nat
284 0 : IF (iBeta==iAlpha) CYCLE
285 0 : tauExtAlpha = MATMUL(cell%amat, atoms%taual(:, iAlpha))
286 0 : vecR = tauExtAlpha - tauExtBeta
287 : constTerm(3*(iBeta-1) + 1:3*(iBeta-1) + 3, 3*(iAlpha-1) + 1:3*(iAlpha-1) + 3) &
288 0 : & = atoms%zatom(iBeta) * atoms%zatom(iBeta) * (3*outerProduct(vecR,vecR)-id3x3*norm2(vecR)**2) / norm2(vecR)**5
289 : END DO
290 : END DO
291 0 : END SUBROUTINE
292 :
293 0 : subroutine CalcIIEnerg2(atoms, cell, qpts, stars, input, iqpt, ngdp, gdp, E2ndOrdII)
294 :
295 : use m_sphbes
296 :
297 : implicit none
298 :
299 :
300 : ! Type parameter
301 : type(t_atoms), intent(in) :: atoms
302 : type(t_cell), intent(in) :: cell
303 : type(t_kpts), intent(in) :: qpts
304 : type(t_stars), intent(in) :: stars
305 : type(t_input), intent(in) :: input
306 :
307 : ! Scalar parameter
308 : integer, intent(in) :: iqpt
309 : integer, intent(in) :: ngdp
310 :
311 : ! Array parameter
312 : integer, intent(in) :: gdp(:, :)
313 : complex, allocatable, intent(out):: E2ndOrdII(:, :)
314 :
315 : ! Scalar variables
316 : integer :: iAtype
317 : integer :: iBtype, iCtype
318 : integer :: iAeqat
319 : integer :: iBeqat
320 : integer :: iAatom
321 : integer :: iBatom
322 : integer :: iAdir
323 : integer :: iBdir, iCdir
324 : integer :: ngpqdp2km
325 : integer :: ngpqdp
326 :
327 : LOGICAL :: oldStuff
328 :
329 : ! Array variables
330 : complex, allocatable :: E2ndOrdIIatFinQ(:, :)
331 : complex, allocatable :: E2ndOrdIIatQ0(:, :)
332 : real, allocatable :: constTerm(:, :)
333 0 : integer, allocatable :: gpqdp(:, :)
334 :
335 0 : oldStuff = .FALSE.
336 :
337 : ! We get the same results for -q and q, probably because Eii2 is a real quantity
338 : ! todo only generate G-vectors once in the beginning #56, leave the -q version here so that we can test it to be the same
339 0 : call genPertPotDensGvecs( stars, cell, input, ngpqdp, ngpqdp2km, -qpts%bk(1:3, iqpt), gpqdp )
340 :
341 : #ifdef DEBUG_MODE
342 : if (.false.) then
343 : call genPertPotDensGvecs( stars, cell, input, ngpqdp, ngpqdp2km, qpts%bk(1:3, iqpt), gpqdp )
344 : end if
345 : #endif
346 :
347 : ! Create final Eii2 matrix and temporary array for passing to the subroutine
348 0 : allocate( E2ndOrdII(3 * atoms%nat, 3 * atoms%nat) )
349 0 : allocate( E2ndOrdIIatFinQ(3 * atoms%nat, 3 * atoms%nat) )
350 0 : allocate( E2ndOrdIIatQ0(3 * atoms%nat, 3 * atoms%nat) )
351 0 : allocate( constTerm(3 * atoms%nat, 3 * atoms%nat) )
352 :
353 0 : E2ndOrdII = cmplx(0.0, 0.0)
354 0 : E2ndOrdIIatFinQ = cmplx(0.0, 0.0)
355 0 : E2ndOrdIIatQ0 = cmplx(0.0, 0.0)
356 :
357 : ! Call the routine for q = 0
358 0 : call CalcIIEnerg2MatElem(atoms, cell, [0.0,0.0,0.0], ngdp, gdp, E2ndOrdIIatQ0)
359 :
360 : ! Call the routine for finite q
361 0 : call CalcIIEnerg2MatElem(atoms, cell, -qpts%bk(1:3, iqpt), ngpqdp, gpqdp, E2ndOrdIIatFinQ)
362 :
363 0 : CALL getConstTerm(atoms, cell, constTerm)
364 :
365 : !write(4543,*) E2ndOrdIIatQ0
366 : !write(4544,*) E2ndOrdIIatFinQ
367 :
368 : #ifdef DEBUG_MODE
369 : if (.false.) then
370 : ! We get the same results for -q and q, probably because Eii2 is a real quantity
371 : call CalcIIEnerg2MatElem(atoms, cell, qpts%bk(1:3, iqpt), ngpqdp, gpqdp, E2ndOrdIIatFinQ)
372 : end if
373 : #endif
374 0 : iAatom = 0
375 0 : do iAtype = 1, atoms%ntype
376 0 : do iAeqat = 1, atoms%neq(iAtype)
377 0 : iAatom = iAatom + 1
378 0 : do iAdir = 1, 3
379 : iBatom = 0
380 0 : do iBtype = 1, atoms%ntype
381 0 : do iBeqat = 1, atoms%neq(iBtype)
382 0 : iBatom = iBatom + 1
383 0 : do iBdir = 1, 3
384 0 : IF (oldStuff) THEN
385 : E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) = &
386 : & E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) &
387 : & - E2ndOrdIIatQ0(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) &
388 : & + E2ndOrdIIatFinQ(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1))
389 : ELSE
390 : E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) = &
391 : & E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) + &
392 0 : & E2ndOrdIIatFinQ(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1))
393 0 : IF (iBatom/=iAatom) THEN
394 : ! E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) = &
395 : ! & E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) - &
396 : ! & constTerm(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1))
397 : ELSE
398 0 : DO iCtype = 1, atoms%ntype
399 : !DO iCdir = 1, 3
400 : E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) = &
401 : & E2ndOrdII(iBdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) - &
402 0 : & E2ndOrdIIatQ0(iBdir + 3 * (iCtype - 1), iAdir + 3 * (iAatom - 1))
403 : !END DO
404 : IF (iCtype==iBtype) CYCLE
405 0 : DO iCdir = 1, 3
406 : ! E2ndOrdII(iCdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) = &
407 : ! & E2ndOrdII(iCdir + 3 * (iBatom - 1), iAdir + 3 * (iAatom - 1)) + &
408 : ! & constTerm(iCdir + 3 * (iCtype - 1), iAdir + 3 * (iAatom - 1))
409 : END DO
410 : END DO
411 : END IF
412 : END IF
413 : end do ! iBdir
414 : end do ! iBeqat
415 : end do ! iBtype
416 : end do ! iAdir
417 : end do ! iAeqat
418 : end do ! iAtype
419 :
420 0 : end subroutine CalcIIEnerg2
421 :
422 : ! Creates interstitial and muffin-tin coefficients of the second-variation external potential required for the Dynamical Matrix using
423 : ! the Weinert method.
424 : ! Main subroutine for creating the second-order V_ext interstitial coefficients and muffin-tin coefficients using the Weinert method.
425 : ! Here the output is a matrix!!! So the factor of 2 is multiplied in the dynamical matrix not before multipliying the Qs.
426 0 : subroutine GenVext2(atoms, cell, ngdp, gdp, vExt2IR, vExt2MT, testMode)
427 :
428 : implicit none
429 :
430 : ! Type parameter
431 : type(t_atoms), intent(in) :: atoms
432 : type(t_cell), intent(in) :: cell
433 :
434 : ! Scalar parameter
435 : integer, intent(in) :: ngdp
436 : logical, intent(in) :: testMode
437 :
438 : ! Array parameter
439 : integer, intent(in) :: gdp(:, :)
440 : complex, allocatable, intent(out):: vExt2IR(:, :, :, :)
441 : complex, allocatable, intent(out):: vExt2MT(:, :, :, :)
442 :
443 : ! Scalar variables
444 : integer :: G0index
445 : integer :: iatom
446 : integer :: itype
447 : integer :: ieqat
448 : integer :: ii
449 : integer :: jj
450 : integer :: iG
451 :
452 : ! Array variables
453 0 : real, allocatable :: Gext(:, :)
454 0 : complex, allocatable :: psDens2ndOrd(:, :, :, :)
455 :
456 0 : allocate(vExt2IR(ngdp, 3, 3, atoms%nat))
457 0 : allocate(Gext(3, ngdp))
458 0 : vExt2IR = cmplx(0.0, 0.0)
459 :
460 0 : call GenPsDens2ndOrd(atoms, cell, ngdp, G0index, gdp, [0., 0., 0.], psDens2ndOrd, testMode)
461 :
462 0 : do iG = 1, ngdp
463 0 : Gext(:, iG) = matmul(gdp(:, iG), cell%bmat)
464 : end do
465 :
466 0 : iatom = 0
467 0 : do itype = 1, atoms%ntype
468 0 : do ieqat = 1, atoms%neq(itype)
469 0 : iatom = iatom + 1
470 0 : do jj = 1, 3
471 0 : do ii = 1, 3
472 0 : do iG = 1, ngdp
473 0 : if ( iG /= G0index ) then
474 0 : vExt2IR(iG, ii, jj, iatom) = fpi_const * psDens2ndOrd(iG, ii, iatom, jj ) / norm2(Gext(:, iG))**2
475 : end if
476 : end do
477 : end do
478 : end do
479 : end do
480 : end do ! iG
481 :
482 0 : call GenVext2MT(atoms, cell, ngdp, gdp, testMode, vExt2IR, vExt2MT)
483 :
484 0 : end subroutine GenVext2
485 :
486 : ! Generates second-order V_ext coefficients for the muffin-tin region.
487 0 : subroutine GenVext2MT(atoms, cell, ngdp, gdp, testMode, vExt2IR, vExt2MT)
488 :
489 : use m_gaunt, only : gaunt1
490 : use m_sphbes
491 : use m_jpGrVeff0, only : Phasy1nSym
492 :
493 : implicit none
494 :
495 : ! Type parameters
496 : type(t_atoms), intent(in) :: atoms
497 : type(t_cell), intent(in) :: cell
498 :
499 : ! Scalar parameters
500 : integer, intent(in) :: ngdp
501 :
502 : ! Array parameters
503 : integer, intent(in) :: gdp(:, :)
504 : logical, intent(in) :: testMode
505 : complex, intent(in) :: vExt2IR(:, :, :, :)
506 : complex, allocatable, intent(out) :: vExt2MT(:, :, :, :)
507 : ! Scalar variables
508 : integer :: itype
509 : integer :: ii
510 : integer :: jj
511 : integer :: oqn_l
512 : integer :: mqn_m
513 : integer :: t
514 : integer :: tPrime
515 : integer :: lm
516 : integer :: imesh
517 : integer :: iatom
518 : integer :: ieqat
519 : integer :: iG
520 : integer :: iDtype
521 : integer :: iDeqat
522 : integer :: iDatom
523 :
524 : !complex(kind=16) :: kindTest ?? This is not portable
525 :
526 : ! Array variables
527 : ! lmax is 2, so the upper limit of lm is (2 + 1)**3 = 9
528 : real :: Gext(3)
529 0 : real :: sbes(0:atoms%lmaxd)
530 0 : real, allocatable :: prfMesh(:, :)
531 0 : complex, allocatable :: pylm(:, :)
532 0 : complex, allocatable :: surfIntMat(:, :, :, :, :)
533 0 : complex, allocatable :: scalFac(:)
534 0 : real, allocatable :: recMesh(:, :)
535 : complex :: volIntMat(9, 3, 3)
536 : complex :: volIntMatTest(9, 3, 3)
537 :
538 : complex :: testMat(3, 3)
539 :
540 0 : allocate( recMesh(atoms%jmtd, atoms%ntype) )
541 0 : recMesh(:, :) = 0.
542 :
543 : ! Calculate the lm dependent matrix part of the 2nd order external potential
544 0 : volIntMat = cmplx(0.0, 0.0)
545 : volIntMatTest = cmplx(0.0, 0.0)
546 :
547 : ! l = 0, m = 0
548 : !if (testMode) then
549 : if (.false.) then
550 : volIntMat(1, 1, 1) = cmplx(sqrt(fpi_const), 0.)
551 : volIntMat(1, 2, 2) = cmplx(sqrt(fpi_const), 0.)
552 : volIntMat(1, 3, 3) = cmplx(sqrt(fpi_const), 0.)
553 : end if
554 :
555 : ! l = 1 has vanishing Gaunt coefficients.
556 :
557 : ! l = 2, m = -2
558 0 : volIntMat(5, 1, 1) = cmplx(sqrt(3 * tpi_const / 5), 0.)
559 0 : volIntMat(5, 2, 2) = cmplx(-sqrt(3 * tpi_const / 5), 0.)
560 0 : volIntMat(5, 1, 2) = cmplx(0., sqrt(3 * tpi_const / 5))
561 0 : volIntMat(5, 2, 1) = cmplx(0., sqrt(3 * tpi_const / 5))
562 :
563 : ! l = 2, m = -1
564 0 : volIntMat(6, 1, 3) = cmplx(sqrt(3 * tpi_const / 5), 0.)
565 0 : volIntMat(6, 3, 1) = cmplx(sqrt(3 * tpi_const / 5), 0.)
566 0 : volIntMat(6, 2, 3) = cmplx(0., sqrt(3 * tpi_const / 5))
567 0 : volIntMat(6, 3, 2) = cmplx(0., sqrt(3 * tpi_const / 5))
568 :
569 : ! l = 2, m = 0
570 0 : volIntMat(7, 1, 1) = cmplx(-sqrt(fpi_const / 5), 0.)
571 0 : volIntMat(7, 2, 2) = cmplx(-sqrt(fpi_const / 5), 0.)
572 0 : volIntMat(7, 3, 3) = cmplx(4 * sqrt(pi_const / 5), 0.)
573 :
574 : ! l = 2, m = 1
575 0 : volIntMat(8, 1, 3) = cmplx(-sqrt(3 * tpi_const / 5), 0.)
576 0 : volIntMat(8, 3, 1) = cmplx(-sqrt(3 * tpi_const / 5), 0.)
577 0 : volIntMat(8, 3, 2) = cmplx(0., sqrt(3 * tpi_const / 5))
578 0 : volIntMat(8, 2, 3) = cmplx(0., sqrt(3 * tpi_const / 5))
579 :
580 : ! l = 2, m = 2
581 0 : volIntMat(9, 1, 1) = cmplx(sqrt(3 * tpi_const / 5), 0.)
582 0 : volIntMat(9, 2, 2) = cmplx(-sqrt(3 * tpi_const / 5), 0.)
583 0 : volIntMat(9, 1, 2) = cmplx(0., -sqrt(3 * tpi_const / 5))
584 0 : volIntMat(9, 2, 1) = cmplx(0., -sqrt(3 * tpi_const / 5))
585 :
586 : if (.false.) then
587 : do jj = 1, 3
588 : do ii = 1, 3
589 : do oqn_l = 0, 2, 2
590 : do t = -1, 1
591 : do tPrime = -1, 1
592 : ! For l = 0, we only have m = 0
593 : if (abs(t + tPrime) > oqn_l) cycle
594 : lm = oqn_l * (oqn_l + 1) + 1 + t + tPrime
595 : ! + 2 because c_im array starts at 1 and not at - 1
596 : volIntMatTest(lm, ii, jj) = volIntMatTest(lm, ii, jj) + 3 * c_im(ii, t + 2) * c_im(jj, tPrime + 2) &
597 : & * gaunt1(oqn_l, 1, 1, t + tPrime, t, tPrime, 2)
598 : end do
599 : end do
600 : end do ! oqn_l
601 : end do ! ii
602 : end do ! jj
603 : ! if (.not.testMode) then
604 : ! volIntMat(1, :, :) = cmplx(0.0, 0.0)
605 : ! end if
606 :
607 : do lm = 1, 9
608 : do jj = 1, 3
609 : do ii = 1, 3
610 : write(1104, '(3(i8),2(f15.8))') jj, ii, lm, volIntMatTest(lm, ii, jj)
611 : write(1105, '(3(i8),2(es15.8))') jj, ii, lm, volIntMatTest(lm, ii, jj) - volIntMat(lm, ii, jj)
612 : end do ! lm
613 : end do ! ii
614 : write(1104, *)
615 : write(1105, *)
616 : end do ! jj
617 :
618 : do jj = 1, 3
619 : write(*, '(i8,1x,2(es15.8))') jj, volIntMat(1, jj, jj)
620 : end do ! lm
621 : write(*, *) 'shack'
622 :
623 : end if
624 :
625 :
626 : ! testMat = 0
627 : ! do lm = 1, 10
628 : ! testMat(:, :) = testMat(:, :) + volIntMat(lm, :, :)
629 : ! end do
630 : ! do ii = 1, 3
631 : ! do jj = 1, 3
632 : ! write(*, *) testMat(jj, ii)
633 : ! end do
634 : ! end do
635 : !
636 : ! NOstopNO
637 : ! do jj = 1, 9
638 : ! write(*, *) jj
639 : ! do ii = 1, 3
640 : ! write(*, '(3(2(f15.8),2x))') volIntMat(jj, ii, 1), volIntMat(jj, ii, 2), volIntMat(jj, ii, 3)
641 : ! end do
642 : ! write(*, *)
643 : ! end do
644 : ! ! do ii = 1, 3
645 : ! ! do jj = 1, 3
646 : ! ! write(*, *) jj, ii, c_im(jj, ii)
647 : ! ! end do
648 : ! ! end do
649 : ! ! NOstopNO
650 : !! write(*, *) 'gaunt'
651 : ! write(*, *) gaunt1(0, 1, 1, 0, 0, 0, 2)
652 : ! write(*, *) gaunt1(0, 1, 1, 0, -1, 1, 2)
653 : ! write(*, *) gaunt1(0, 1, 1, 0, 1, -1, 2)
654 : ! NOstopNO
655 : ! NOstopNO
656 : ! write(*, *) gaunt1(0, 1, 1, 0, -1, 1, 2)
657 : ! NOstopNO
658 : ! write(*, *) 1. / 4. / 3.14 ! NOstopNO
659 : ! do t = 1, 3
660 : ! do tPrime = 1, 3
661 : ! write(*, *) tPrime, t, c_im(tPrime, t)
662 : ! end do
663 : ! end do
664 : ! NOstopNO
665 :
666 :
667 : ! Calculate the mesh dependent prefactor
668 0 : allocate( prfMesh(atoms%jmtd, atoms%ntype) )
669 0 : prfMesh = cmplx(0.0, 0.0)
670 0 : do itype = 1, atoms%ntype
671 : ! Numerical optimization. At the MT boundary, we set the factor zero because it is analytically zero to prevent error propagation.
672 0 : do imesh = 1, atoms%jri(itype) - 1
673 0 : prfMesh(imesh, itype) = atoms%zatom(itype) * (1 - (atoms%rmsh(imesh, itype) / atoms%rmt(itype))**5)
674 0 : recMesh(imesh, itype) = atoms%rmsh(imesh, itype)**(-3)
675 : end do
676 : end do
677 :
678 : ! Calculate the surface integral part of the Dirichelet boundary problem
679 0 : allocate( surfIntMat((atoms%lmaxd + 1)**2, 3, 3, atoms%nat, atoms%nat) )
680 0 : allocate( pylm(( atoms%lmaxd + 1 )**2, atoms%nat) )
681 0 : allocate( scalFac(( atoms%lmaxd + 1 )**2) )
682 0 : surfIntMat(:, :, :, :, :) = cmplx(0.0, 0.0)
683 0 : pylm = cmplx(0.0, 0.0)
684 0 : scalFac = cmplx(0.0, 0.0)
685 : iatom = 0
686 0 : do itype = 1, atoms%ntype
687 0 : do ieqat = 1, atoms%neq(itype)
688 0 : iatom = iatom + 1
689 0 : do iG = 1, ngdp
690 0 : if ( all( gdp(:, iG) == 0 ) ) then
691 : !todo we should determine the G=0 vector once or put it to the beginning at all
692 : cycle
693 : end if
694 0 : Gext(1:3) = matmul( gdp(1:3, iG), cell%bmat(1:3, 1:3) )
695 0 : pylm(:, :) = cmplx(0.0, 0.0)
696 0 : scalFac(:) = cmplx(0.0, 0.0)
697 0 : sbes(:) = cmplx(0., 0.)
698 0 : call phasy1nSym( atoms, cell, gdp(:, iG), [0., 0., 0.], pylm )
699 0 : call Sphbes( atoms%lmax(itype), norm2( Gext ) * atoms%rmt(itype), sbes )
700 : ! If we precalcalculate the scalar factor we have 7 multiplications less
701 0 : do oqn_l = 0, atoms%lmax(itype)
702 0 : do mqn_m = -oqn_l, oqn_l
703 0 : lm = oqn_l * (oqn_l + 1) + mqn_m + 1
704 0 : scalFac(lm) = sbes(oqn_l) * pylm(lm, iatom)
705 : end do
706 : end do
707 : iDatom = 0
708 0 : do iDtype = 1, atoms%ntype
709 0 : do iDeqat = 1, atoms%neq(iDtype)
710 0 : iDatom = iDatom + 1
711 0 : do jj = 1, 3
712 0 : do ii = 1, 3
713 0 : do oqn_l = 0, atoms%lmax(itype)
714 0 : do mqn_m = -oqn_l, oqn_l
715 0 : lm = oqn_l * (oqn_l + 1) + mqn_m + 1
716 : surfIntMat(lm, ii, jj, iDatom, iatom) = surfIntMat(lm, ii, jj, iDatom, iatom) + vExt2IR(iG, ii, jj, iDatom) &
717 0 : & * scalFac(lm)
718 : end do
719 : end do
720 : end do
721 : end do
722 : end do
723 : end do
724 : end do
725 : end do
726 : end do
727 0 : deallocate(pylm, scalFac)
728 :
729 : ! Add the volume and surface integral contribution of the Dirichelet boundary problem for determing the MT coefficients of the 2nd
730 : ! order external potential
731 : ! We have chosen this loop structure sothat the if-statement is as outermost as possible
732 0 : allocate(vExt2MT(atoms%jmtd, (atoms%lmaxd + 1)**2, 3 * atoms%nat, 3 * atoms%nat ))
733 0 : vExt2MT = cmplx(0.0, 0.0)
734 : iDatom = 0
735 0 : do iDtype = 1, atoms%ntype
736 0 : do iDeqat = 1, atoms%neq(iDtype)
737 0 : iDatom = iDatom + 1
738 0 : do jj = 1, 3
739 : iatom = 0
740 0 : do itype = 1, atoms%ntype
741 0 : do ieqat = 1, atoms%neq(itype)
742 0 : iatom = iatom + 1
743 0 : if ( iatom == iDatom ) then
744 0 : do ii = 1, 3
745 : ! if trace is subtracted there is no l = 0 component
746 : !do oqn_l = 0, 2, 2
747 0 : do oqn_l = 2, 2
748 : !do oqn_l = 0, 2
749 0 : do mqn_m = -oqn_l, oqn_l
750 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
751 0 : do imesh = 1, atoms%jri(itype)
752 : vExt2MT(imesh, lm, ii + (iatom - 1) * 3, jj + (iDatom - 1) * 3) = &
753 : & vExt2MT(imesh, lm, ii + (iatom - 1) * 3, jj + (iDatom - 1) * 3) &
754 : ! We calculate prfMesh * volIntMat first due to numerical stability
755 0 : & - (prfMesh(imesh, itype) * volIntMat(lm, ii, jj)) * recMesh(imesh, itype)
756 : end do ! imesh
757 : end do ! mqn_m
758 : end do ! oqn_l
759 : end do ! ii
760 : end if ! within displaced muffin-tin
761 :
762 : ! Add surface integral contribution.
763 0 : do ii = 1, 3
764 0 : do oqn_l = 0, atoms%lmax(itype)
765 0 : do mqn_m = -oqn_l, oqn_l
766 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
767 0 : do imesh = 1, atoms%jri(itype)
768 : vExt2MT(imesh, lm, ii + (iatom - 1) * 3, jj + (iDatom - 1) * 3) = &
769 : & vExt2MT(imesh, lm, ii + (iatom - 1) * 3, jj + (iDatom - 1) * 3) &
770 0 : & + (atoms%rmsh(imesh, itype) / atoms%rmt(itype))**oqn_l * surfIntMat(lm, ii, jj, iDatom, iatom)
771 : end do ! imesh
772 : end do ! mqn_m
773 : end do ! oqn_l
774 : end do ! ii
775 :
776 : end do ! ieqat
777 : end do ! itype
778 : end do ! ii
779 : end do ! iDeqat
780 : end do ! iDtype
781 :
782 0 : end subroutine GenVext2MT
783 :
784 : !--------------------------------------------------------------------------------------------------------------------------------------
785 : !> @author
786 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
787 : !>
788 : !> @brief
789 : !> Auxiliary method which is similiar to FLEUR's phasy routine but does not implement symmetry.
790 : !>
791 : !>
792 : !> @details
793 : !> This method calculates for a given G-vector and q-vector: 4 π i^l exp((Gvec + qpoint) * taual) Y*_lm(Gvec + qpoint)
794 : !> Opposite to the original phasy1 routine, here, stars and any type of symmetry is not considered here.
795 : !>
796 : !> @param[in] atomsT : Contains atoms-related quantities; definition of its members in types.F90 file.
797 : !> @param[in] cellT : Contains unit-cell related quantities; definition of its members in types.F90 file.
798 : !> @param[in] Gvec : Current G-vector
799 : !> @param[in] qpoint : Current q-point
800 : !> @param[out] pylm : Result of the routine, is atom resolved due to lack of symmetry in first variation of the potentials.
801 : !--------------------------------------------------------------------------------------------------------------------------------------
802 : ! Deprecated
803 0 : subroutine phasy1lp2nSym(atomsT, cellT, Gvec, qptn, pylm)
804 :
805 : use m_ylm
806 : use m_types
807 :
808 : implicit none
809 :
810 : ! Scalar Type Arguments
811 : type(t_atoms), intent(in) :: atomsT
812 : type(t_cell), intent(in) :: cellT
813 :
814 : ! Array Arguments
815 : integer, intent(in) :: Gvec(:)
816 : real, intent(in) :: qptn(:)
817 : complex, intent(out) :: pylm((atomsT%lmaxd + 3)**2, atomsT%nat)
818 :
819 : !---------------------------------------------------------------------------------------------------------------------------------
820 : ! Local Scalar Variables
821 : ! iatom : runs over all atoms
822 : ! itype : runs over all types
823 : ! ineq : runs over all equivalent atoms of one atom type
824 : ! lm : encodes oqn_l and mqn_m
825 : ! sf : stores exponential function
826 : ! csf : stores exponential function times 4 π i^l
827 : ! x : stores argument of exponential function
828 : ! mqn_m : magnetic quantum number m
829 : ! oqn_l : orbital quantum number l
830 : ! ll1 : auxiliary variable to calculate lm
831 : !---------------------------------------------------------------------------------------------------------------------------------
832 : integer :: iatom
833 : integer :: itype
834 : integer :: ineq
835 : integer :: lm
836 : complex :: sf
837 : complex :: csf
838 : real :: x
839 : integer :: mqn_m
840 : integer :: oqn_l
841 : integer :: ll1
842 :
843 : !---------------------------------------------------------------------------------------------------------------------------------
844 : ! Local Array Variables
845 : ! fpiul: stores 4 π i^l
846 : ! Gqext: stores G + q in external coordinates
847 : ! ylm : stores Y_lm
848 : !---------------------------------------------------------------------------------------------------------------------------------
849 0 : complex :: fpiul(0:atomsT%lmaxd + 2)
850 : real :: Gqext(3)
851 0 : complex :: ylm((atomsT%lmaxd + 3)**2)
852 :
853 : ! calculates 4 π i^l resolved for every l, not divided by nop because no loop over symmetry operations
854 0 : do oqn_l = 0, atomsT%lmaxd + 2
855 0 : fpiul(oqn_l) = fpi_const * ImagUnit**oqn_l
856 : enddo
857 :
858 :
859 : ! calculates Y*_lm(\vec{G} + \vec{q}) for every l and m. The argument Gqext must be in external coordinates.
860 0 : Gqext = matmul(Gvec + qptn, cellT%bmat)
861 : !call Ylmnorm_init( atomsT%lmaxd + 2 )
862 0 : call ylm4(atomsT%lmaxd + 2, Gqext, ylm)
863 : !call Ylmnorm_init( atomsT%lmaxd)
864 0 : ylm = conjg(ylm)
865 :
866 :
867 : ! calculates first exp(i (G + q) tau) and multiplies recent factors before storing the final result to pylm
868 0 : iatom = 1
869 0 : pylm = cmplx(0.,0.)
870 0 : do itype = 1, atomsT%ntype
871 0 : do ineq = 1, atomsT%neq(itype)
872 0 : x = tpi_const * dot_product(Gvec + qptn, atomsT%taual(:, iatom))
873 0 : sf = exp(ImagUnit * x)
874 0 : do oqn_l = 0, atomsT%lmax(itype) + 2
875 0 : ll1 = oqn_l * (oqn_l + 1) + 1
876 0 : csf = fpiul(oqn_l) * sf
877 0 : do mqn_m = -oqn_l, oqn_l
878 0 : lm = ll1 + mqn_m
879 0 : pylm(lm, iatom) = csf * ylm(lm)
880 : enddo ! mqn_m
881 : enddo ! oqn_l
882 0 : iatom = iatom + 1
883 : enddo ! ineq
884 : enddo ! itype
885 :
886 0 : end subroutine phasy1lp2nSym
887 :
888 0 : function outerProduct(a, b)
889 :
890 : implicit none
891 :
892 : real, intent(in) :: a(:)
893 : real, intent(in) :: b(:)
894 :
895 : real :: outerProduct(size(a), size(b))
896 :
897 0 : outerProduct(:, :) = spread(a, dim=2, ncopies=size(b)) * spread(b, dim=1, ncopies=size(a))
898 :
899 : end function outerProduct
900 :
901 0 : function outerProductME(a, b, i, j)
902 :
903 : use m_juDFT_stop, only : juDFT_error
904 :
905 : implicit none
906 :
907 : complex, intent(in) :: a(:)
908 : complex, intent(in) :: b(:)
909 : integer, intent(in) :: i
910 : integer, intent(in) :: j
911 :
912 : complex :: outerProductME
913 :
914 0 : if (i > size(a) .or. j > size(b)) then
915 : call juDFT_error( 'Wished matrix element is out of scope of outer product matrix.', calledby='outerProductME', &
916 0 : & hint='Choose smaller indices.')
917 : end if
918 :
919 0 : outerProductME = a(i) * b(j)
920 :
921 0 : end function outerProductME
922 :
923 0 : subroutine genPertPotDensGvecs( stars, cell, input, ngpqdp, ngpqdp2km, qpoint, gpqdp )
924 :
925 : use m_types
926 :
927 : implicit none
928 :
929 : ! Type parameters
930 : type(t_stars), intent(in) :: stars
931 : type(t_cell), intent(in) :: cell
932 : type(t_input), intent(in) :: input
933 :
934 : ! Scalar parameters
935 : integer, intent(out) :: ngpqdp
936 : integer, intent(out) :: ngpqdp2km
937 :
938 : ! Array parameters
939 : real, intent(in) :: qpoint(:)
940 : integer, allocatable, intent(out) :: gpqdp(:, :)
941 :
942 : ! Scalar variables
943 : integer :: ngrest
944 : integer :: iGx
945 : integer :: iGy
946 : integer :: iGz
947 : integer :: iG
948 :
949 : ! Array variables
950 0 : integer, allocatable :: gpqdptemp2kmax(:, :)
951 0 : integer, allocatable :: gpqdptemprest(:, :)
952 : integer :: Gint(3)
953 : real :: Gpqext(3)
954 :
955 : allocate( gpqdptemp2kmax(3, (2 * stars%mx1 + 1) * (2 * stars%mx2 + 1) * (2 * stars%mx3 + 1)), &
956 0 : & gpqdptemprest(3, (2 * stars%mx1 + 1) * (2 * stars%mx2 + 1) * (2 * stars%mx3 + 1)) )
957 :
958 0 : ngpqdp = 0
959 0 : ngpqdp2km = 0
960 0 : ngrest = 0
961 0 : gpqdptemp2kmax(:, :) = 0
962 0 : gpqdptemprest(:, :) = 0
963 : ! From all possible G-vectors in a box, only these are accepted which are element of a sphere with radius gmax which is shifted.
964 : ! We need a little bit more than k*d because they are thought for a Gmax ball that is not shifted by a q, i.e. |G+q|<Gmax
965 0 : do iGx = -(stars%mx1 + 3), (stars%mx1 + 3)
966 0 : do iGy = -(stars%mx2 + 3), (stars%mx2 + 3)
967 0 : do iGz = -(stars%mx3 + 3), (stars%mx3 + 3)
968 0 : Gint = [iGx, iGy, iGz]
969 0 : Gpqext = matmul(real(Gint(1:3) + qpoint(1:3)),cell%bmat) !transform from internal to external coordinates
970 0 : if (norm2(Gpqext) <= input%gmax) then
971 0 : ngpqdp = ngpqdp + 1
972 : ! Sort G-vectors
973 0 : if ( norm2(Gpqext) <= 2 * input%rkmax ) then
974 0 : ngpqdp2km = ngpqdp2km + 1
975 0 : gpqdptemp2kmax(1:3, ngpqdp2km) = Gint(1:3)
976 : else
977 0 : ngrest = ngrest + 1
978 0 : gpqdptemprest(1:3, ngrest) = Gint(1:3)
979 : end if
980 : end if
981 : end do !iGz
982 : end do !iGy
983 : end do !iGx
984 0 : allocate(gpqdp(3, ngpqdp))
985 : ! Mapping array from G-vector to G-vector index
986 0 : gpqdp(:, :) = 0
987 0 : gpqdp(1:3, 1:ngpqdp2km) = gpqdptemp2kmax(1:3, 1:ngpqdp2km)
988 0 : gpqdp(1:3, ngpqdp2km + 1 : ngpqdp) = gpqdptemprest(1:3, 1:ngrest)
989 :
990 0 : end subroutine genPertPotDensGvecs
991 :
992 0 : SUBROUTINE dfpt_e2_madelung(atoms,jspins,grgrRho,grgrVC,e2_vm)
993 : USE m_intgr
994 : TYPE(t_atoms), INTENT(IN) :: atoms
995 : INTEGER, INTENT(IN) :: jspins
996 : REAL, INTENT(IN) :: grgrRho(:,:,:), grgrVC(:,:)
997 : REAL, INTENT(OUT) :: e2_vm(:)
998 :
999 0 : REAL :: mt(atoms%jmtd,atoms%nat), dpj(atoms%jmtd), rhs_arr(atoms%jmtd), totz_arr(atoms%jmtd)
1000 0 : REAL :: vmd(atoms%nat), zintn_r(atoms%nat)
1001 : INTEGER :: n, j
1002 :
1003 0 : e2_vm = 0.0
1004 0 : mt=0.0
1005 0 : DO n = 1,atoms%nat
1006 0 : DO j = 1,atoms%jri(n)
1007 0 : mt(j,n) = grgrRho(j,n,1) + grgrRho(j,n,jspins)
1008 : ENDDO
1009 : ENDDO
1010 0 : IF (jspins.EQ.1) mt=mt/2
1011 :
1012 0 : rhs_arr = 0.0
1013 0 : DO n = 1,atoms%nat
1014 0 : DO j = 1,atoms%jri(n)
1015 0 : dpj(j) = mt(j,n)/atoms%rmsh(j,n)
1016 : END DO
1017 : !CALL intgr3(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),rhs)
1018 0 : CALL sfint(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),rhs_arr)
1019 :
1020 0 : zintn_r(n) = atoms%zatom(n)*sfp_const*rhs_arr(atoms%jri(n))
1021 0 : e2_vm(n) = e2_vm(n) - zintn_r(n)
1022 :
1023 : !CALL intgr3(mt(1,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz)
1024 0 : CALL sfint(mt(:,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz_arr)
1025 0 : vmd(n) = atoms%rmt(n)*grgrVC(atoms%jri(n),n)/sfp_const - totz_arr(atoms%jri(n))*sfp_const
1026 0 : vmd(n) = -atoms%zatom(n)*vmd(n)/ atoms%rmt(n)
1027 0 : e2_vm(n) = e2_vm(n) + vmd(n)
1028 : END DO
1029 :
1030 0 : END SUBROUTINE dfpt_e2_madelung
1031 :
1032 : end module m_dfpt_eii2
|
