Line data Source code
1 : !----------------------------------------------------------------------------------------------------------------------------------------
2 : ! Forschungszentrum Jülich, juPhon Plugin for the FLEUR program
3 : !----------------------------------------------------------------------------------------------------------------------------------------
4 : !
5 : ! MODULE: Gradient of unperturbed effective potential
6 : !
7 : !> @author
8 : !> Christian-Roman Gerhorst
9 : !>
10 : !> @note
11 : !> Many routines are very similiar to the routines of FLEUR.
12 : !>
13 : !> @brief
14 : !> This module contains routines related to the calculation of the gradient of the unperturbed effective potential.
15 : !>
16 : !> @todo
17 : !> Complete documentation
18 : !> Put in dimensions of array?
19 : !>
20 : !> @note
21 : !> Additional information and formulas pointing out the routines of this module can be found within this
22 : !> <a href='jpGrVeff0.pdf'>document</a>.
23 : !----------------------------------------------------------------------------------------------------------------------------------------
24 : module m_jpGrVeff0
25 :
26 : USE m_constants
27 : implicit none
28 :
29 : contains
30 :
31 : !--------------------------------------------------------------------------------------------------------------------------------------
32 : !> @author
33 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
34 : !>
35 : !> @brief
36 : !> Main routine to calculate the gradient of the unperturbed effective potential.
37 : !>
38 : !> @details
39 : !> The first variation of the Coulomb potential is calculated using the Weinert method without using symmetry by not exploiting
40 : !> stars and lattice harmonics according to M.Weinert, J.Math.Phys. 22(11) (1981) p.2434 and the PhD thesis of Aaron Klüppelberg.
41 : !> For the xc-potential, the functional derivative of the kernel with respect to the density is connected with the gradient of the
42 : !> density. In the interstitial region, this is done by a convolution, whereas in the muffin-tin a projection onto a spherical
43 : !> harmonic is evaluted by using Gauss quadrature.
44 : !>
45 : !> @note
46 : !> At the moment, only the x-alpha kernel functional derivative is implemented, but its calculation is in seperate routines, that
47 : !> can be replaced by the libxc library later.
48 : !>
49 : !> @param[in] atoms : Atoms type, see types.f90
50 : !> @param[in] cell : Unit cell type, see types.f90
51 : !> @param[in] dimensions : Dimension type, see types.f90
52 : !> @param[in] stars : Stars type, see types.f90
53 : !> @param[in] ngdp : Number of G-vectors for potentials and densities
54 : !> @param[in] harSw : Switch to enable Hartree contribution
55 : !> @param[in] extSw : Switch to enable external contribution
56 : !> @param[in] xcSw : Switch to enable xc contribution
57 : !> @param[in] testGoldstein : Switch needed by a test using the Goldstone condition
58 : !> @param[in] grRhoTermSw : Switch to activate volume term contribution of muffin-tin boundary problem
59 : !> @param[in] gdp : G-vectors of potentials and densities
60 : !> @param[in] rho0IRpw : Plane-wave coefficients of the unperturbed and converged interstitial density parsed from Fleur
61 : !> @param[in] rho0MTsh : Spherical harmonic coefficients of unperturbed and converged muffin-tin density parsed from Fleur
62 : !> @param[in] grRho0IR : Plane-wave coefficients of the gradient of the interstitial unperturbed density
63 : !> @param[in] grRho0MT : Spherical harmonic coefficients of the gradient of the muffin-tin unperturbed density
64 : !> @param[in] grVxcIRKern : Plane-wave interstitial coefficients of the functional derivative of the xc-kernel with respect to
65 : !> the density
66 : !> @param[in] gausWts : Gauss weights for Gauss quadrature created in m_jppotdenshelper::calcmtdvxckern
67 : !> @param[in] ylm : Set of spherical harmonics whose arguments are the unit vectors of the Gauss mesh points up to lmax
68 : !> created in m_jppotdenshelper::calcmtdvxckern
69 : !> @param[in] dKernMTGpts : Spherical harmonic muffin-tin coefficients of the functional derivative of the xc-kernel with
70 : !> respect to the density created in m_jppotdenshelper::calcmtdvxckern
71 : !> @param[out] grVeff0IR : Plane-wave insterstitial coefficients of the gradient of the unperturbed effective potential
72 : !> @param[out] grVeff0MT : Spherical harmonic muffin-tin coefficients of the gradient of the unperturbed effective potential
73 : !--------------------------------------------------------------------------------------------------------------------------------------
74 0 : subroutine GenGrVeff0( atoms, cell, stars, ngdp, harSw, extSw, xcSw, gdp, rho0IRpw, rho0MTsh, grRho0IR, grRho0MT, &
75 0 : & gausWts, ylm, dKernMTGPts, grVxcIRKern, testGoldstein, grRhoTermSw, grVeff0IR, grVeff0MT )
76 :
77 : USE m_constants
78 : use m_types, only : t_atoms, t_cell, t_stars
79 :
80 : implicit none
81 :
82 : ! Type parameters
83 : type(t_atoms), intent(in) :: atoms
84 : type(t_cell), intent(in) :: cell
85 : type(t_stars), intent(in) :: stars
86 :
87 : ! Scalar parameter
88 : integer, intent(in) :: ngdp
89 : logical, intent(in) :: harSw
90 : logical, intent(in) :: extSw
91 : logical, intent(in) :: xcSw
92 : logical, intent(in) :: testGoldstein
93 : logical, intent(in) :: grRhoTermSw
94 :
95 : ! Array parameters
96 : integer, intent(in) :: gdp(:, :)
97 : complex, intent(in) :: rho0IRpw(:, :)
98 : complex, intent(in) :: rho0MTsh(:, :, :, :)
99 : complex, intent(in) :: grRho0IR(:, :)
100 : complex, intent(in) :: grRho0MT(:, :, :, :)
101 : complex, intent(in) :: grVxcIRKern(:)
102 : real, intent(in) :: gausWts(:)
103 : complex, intent(in) :: ylm(:, :)
104 : real, intent(in) :: dKernMTGPts(:, :, :)
105 : complex, allocatable, intent(out) :: grVeff0IR(:, :)
106 : complex, allocatable, intent(out) :: grVeff0MT(:, :, :, :)
107 :
108 : ! Local scalar variables
109 : real :: normGext
110 : integer :: iG
111 : integer :: idir
112 : integer :: iatom
113 : integer :: itype
114 : integer :: ieqat
115 : integer :: lm
116 : integer :: oqn_l
117 : integer :: mqn_m
118 : integer :: nfftx
119 : integer :: nffty
120 : integer :: nfftz
121 : integer :: nfftxy
122 : integer :: GxFFT
123 : integer :: GyFFT
124 : integer :: GzFFT
125 : integer :: imesh
126 : logical :: didvext
127 :
128 : ! Local array variables
129 0 : complex, allocatable :: qlmGrVc0(:, :, :)
130 0 : complex, allocatable :: qlmGrVh0Vol(:, :, :)
131 0 : complex, allocatable :: psqGrVc0(:, :)
132 0 : complex, allocatable :: grVxc0IR(:, :)
133 0 : complex, allocatable :: grVxc0MT(:, :, :, :)
134 0 : integer, allocatable :: pdG2FouM(:)
135 0 : complex, allocatable :: qlmGrVh0Surf(:, :)
136 : real :: Gext(3)
137 :
138 :
139 0 : if (harSw .or. extSw) then
140 0 : allocate( qlmGrVc0(( atoms%lmaxd + 1 )**2, atoms%nat, 3) )
141 0 : qlmGrVc0(:, :, :) = cmplx(0., 0.)
142 : end if ! harSw .or. extSw
143 :
144 0 : if (harSw) then
145 : ! Multipole moments to calculate the gradient of unperturbed Hartree potential using formulas 7.58, 7.59, 7.60, 4.17 and 3.30
146 : ! from PhDthesAK.
147 : ! NOTE: in Equation 4.17 within PhDthesAK, there are typos and errors, one should consult PhDthesCRG
148 0 : allocate( qlmGrVh0Vol(( atoms%lmaxd + 1 )**2, atoms%nat, 3) )
149 0 : qlmGrVh0Vol(:, :, :) = cmplx(0., 0.)
150 0 : call CalcQlmGrVh0Vol( atoms, cell, ngdp, gdp, grRho0IR, grRho0MT, qlmGrVh0Vol )
151 :
152 0 : do idir = 1, 3
153 0 : iatom = 0
154 0 : do itype = 1, atoms%ntype
155 0 : do ieqat = 1, atoms%neq(itype)
156 0 : iatom = iatom + 1
157 0 : do oqn_l = 0, atoms%lmax(itype)
158 0 : do mqn_m = -oqn_l, oqn_l
159 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
160 0 : qlmGrVc0(lm, iatom, idir) = qlmGrVc0(lm, iatom, idir) + qlmGrVh0Vol(lm, iatom, idir)
161 : end do ! mqn_m
162 : end do ! oqn_l
163 : end do ! ieqat
164 : end do !itype
165 : end do ! idir
166 :
167 : ! Additional multipole moment dependent on discontinuity of density at the MT boundary that is shared for gradient of
168 : ! unperturbed Hartree potential and the linear variation of the Hartree potential using formulas 7.54, 3.30 and 7.56 from
169 : ! PhDthesAK
170 0 : iatom = 0
171 0 : do itype = 1, atoms%ntype
172 0 : do ieqat = 1, atoms%neq(itype)
173 0 : iatom = iatom + 1
174 :
175 0 : allocate( qlmGrVh0Surf( ( atoms%lmax(itype) + 1)**2, 3 ) )
176 0 : qlmGrVh0Surf(:, :) = cmplx(0., 0.)
177 0 : call CalcQlmHarSurf( atoms, cell, itype, iatom, ngdp, gdp, rho0IRpw, rho0MTsh, qlmGrVh0Surf )
178 :
179 0 : do idir = 1, 3
180 0 : do oqn_l = 0, atoms%lmax(itype)
181 0 : do mqn_m = -oqn_l, oqn_l
182 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
183 0 : qlmgrvc0(lm, iatom, idir) = qlmgrvc0(lm, iatom, idir) - qlmGrVh0Surf(lm, idir)
184 : end do ! mqn_m
185 : end do ! oqn_l
186 : end do ! idir
187 0 : deallocate(qlmGrVh0Surf)
188 :
189 : end do ! ieqat
190 : end do ! itype
191 :
192 : end if ! harSw
193 :
194 : ! N_check: qlmgrvc0(lm, iatom, idir)-=\sum{\alpha} qlmGrVh0Surf(lm, idir)[\alpha]
195 : !
196 : ! Consistent with Aaron (7.62) and (2.41) in CRG 16.12.2020-19:00
197 :
198 : ! the atom loop is over all atoms alpha which are displaced this has to be considered later when adding to other quantities
199 0 : if ( extSw ) then
200 0 : do idir = 1, 3
201 0 : iatom = 0
202 0 : do itype = 1, atoms%ntype
203 0 : do ieqat = 1, atoms%neq(itype)
204 0 : iatom = iatom + 1
205 0 : do lm = 2, 4
206 : ! This is + instead of - because we have a - from 7.45 which can be drawn until here.
207 0 : qlmGrVc0(lm, iatom, idir) = qlmGrVc0(lm, iatom, idir) + atoms%zatom(itype) * 3 / 4 / pi_const * c_im(idir, lm - 1)
208 : end do ! lm
209 : end do ! ieqat
210 : end do ! itype
211 : end do ! idir
212 : end if ! extSw
213 :
214 : ! N_check: qlmGrVc0(1m, iatom, idir)+=\frac{3Z_{\alpha}}{4\pi}c_im(idir,m)
215 : !
216 : ! Consistent with Aaron (7.38) and (2.21) in CRG 16.12.2020-19:00 [- omitted; also in (2.27a)]
217 :
218 : ! Here, first the plane-wave contribution to the pseudocharge of the gradient of the unperturbed Hartree potential is
219 : ! evaluated. Then the pseudo charges are calculated for the Hartree case, as well as for the gradient of the unperturbed
220 : ! external potential. Finally, the pseudocharges (in the case of the external potential already summed over all atoms)
221 : ! are used to get the interstitial contribution of the potentials discussed here. According to 7.45b the external potential part
222 : ! needs a minus.
223 :
224 : ! Calculate the pseudo-charge for the Coulomb potential in general or the Hartree potential or the external potential
225 0 : if (harSw .or. extSw) then
226 :
227 0 : allocate(psqGrVc0(ngdp, 3))
228 0 : psqGrVc0(:, :) = cmplx(0., 0.)
229 :
230 0 : call psqpwVeclp1(atoms, cell, ngdp, grRho0IR, harSw, extSw, gdp, qlmGrVc0, psqGrVc0)
231 :
232 0 : deallocate(qlmGrVc0)
233 :
234 : end if ! harSw .or. extSw
235 :
236 0 : allocate(grVeff0IR(ngdp, 3))
237 0 : grVeff0IR(:, :) = cmplx(0., 0.)
238 0 : if ( harSW .or. extSw ) then
239 0 : do idir = 1, 3
240 0 : do iG = 1, ngdp
241 0 : Gext = matmul(cell%bmat, gdp(:, iG))
242 0 : normGext = norm2(Gext)
243 0 : if (normGext == 0) then
244 : cycle
245 : end if
246 0 : grVeff0IR(iG, idir) = fpi_const * psqGrVc0(iG, idir) / normGext**2
247 : end do ! iG
248 : end do ! idir
249 0 : deallocate(psqGrVc0)
250 : end if ! harSW .or. extSw
251 :
252 0 : allocate(grVeff0MT(atoms%jmtd, (atoms%lmaxd + 1)**2, 3, atoms%nat)) ! todo change order of direction and nat
253 0 : grVeff0MT(:, :, :, :) = cmplx(0., 0.)
254 0 : if ( harSW .or. extSw ) then
255 :
256 0 : iatom = 0
257 0 : do itype = 1, atoms%ntype
258 0 : do ieqat = 1, atoms%neq(itype)
259 0 : iatom = iatom + 1
260 : ! Note: We only have to subtract the external volume integral contribution in the displaced atom. However, the gradient of
261 : ! the external potential or the effective potential is only used in the muffin-tin that is displaced, except for the
262 : ! the case when it is plotted, but this does not justify another atom dimension. For plotting we have to create a
263 : ! new array in the plotting routine where grRhoTermSw is enabled in the displaced atom and disabled anywhere else.
264 : !
265 : ! Tests have shown that for kind=8 accuracy, it is not important whether one calculates the external and the
266 : ! Hartree potential seperately or not. By default and to have consistency with Fleur the interstitial boundary value
267 : ! is the Coulomb potential, i.e. the Hartree and the external potential.
268 : call vmtsCoul(atoms, cell, ngdp, iatom, itype, harSw, extSW, grVeff0IR, grRho0MT, gdp, &
269 0 : & grVeff0MT(:, :, :, iatom), testGoldstein, grRhoTermSw )
270 : end do ! ieqat
271 : end do ! itype
272 :
273 : end if ! harSW .or. extSw
274 :
275 0 : if ( xcSw ) then
276 : ! Prepare mapping array between a G-vector and its respective point on the FFT mesh, according to how the FFT routine requires it.
277 0 : allocate(pdG2FouM(ngdp))
278 0 : pdG2FouM(:) = 0
279 0 : nfftx = 3 * stars%mx1
280 0 : nffty = 3 * stars%mx2
281 0 : nfftz = 3 * stars%mx3
282 0 : nfftxy= 9 * stars%mx1 * stars%mx2
283 :
284 0 : do iG = 1, ngdp
285 0 : GxFFT = gdp(1, iG)
286 0 : GyFFT = gdp(2, iG)
287 0 : GzFFT = gdp(3, iG)
288 0 : if (GxFFT < 0) GxFFT = GxFFT + nfftx
289 0 : if (GyFFT < 0) GyFFT = GyFFT + nffty
290 0 : if (GzFFT < 0) GzFFT = GzFFT + nfftz
291 0 : pdG2FouM(iG) = GxFFT + GyFFT * nfftx + GzFFT * nfftxy
292 : end do
293 0 : allocate(grVxc0IR(ngdp, 3))
294 0 : grVxc0IR(:, :) = cmplx(0., 0.)
295 0 : do idir = 1, 3
296 0 : call convolGrRhoKern(stars, ngdp, ngdp, idir, grRho0IR(:, idir), grVxcIRKern, pdG2FouM, pdG2FouM, grVxc0IR(:, idir), 1)
297 : end do
298 0 : deallocate(pdG2FouM)
299 :
300 : ! Add x-alpha xc-potential contribution for IR
301 0 : do idir = 1, 3
302 0 : do iG = 1, ngdp
303 0 : grVeff0IR(iG, idir) = grVeff0IR(iG, idir) + grVxc0IR(iG, idir)
304 : end do ! iG
305 : end do ! idir
306 0 : deallocate(grVxc0IR)
307 :
308 : ! During the normal Sternheimer iterations the density variation is passed to the routines without the gradient of the density
309 : ! coming from the basis set variation. Therefore, we do not calculate terms accounting for them an switch off terms canceling
310 : ! anyway.
311 0 : if ( grRhoTermSw ) then
312 : ! Add x-alpha xc-potential contribution for MT
313 0 : allocate( grVxc0MT(atoms%jmtd, ( atoms%lmaxd + 1)**2, atoms%nat, 3) )
314 0 : grVxc0MT(:, :, :, :) = cmplx(0., 0.)
315 :
316 : ! Note: we leave the lower block of dead code, because in near future, there will be an optimization, after which it will
317 : ! become required to divide out the r^2 again.
318 : !do idir = 1, 3
319 : ! iatom = 0
320 : ! do itype = 1, atoms%ntype
321 : ! do ieqat = 1, atoms%neq(itype)
322 : ! iatom = iatom + 1
323 : ! do oqn_l = 0, atoms%lmax(itype) + 1
324 : ! do mqn_m = -oqn_l, oqn_l
325 : ! lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
326 : ! do imesh = 1, atoms%jri(itype)
327 : ! grRho0MT(imesh, lm, iatom, idir) = grRho0MT(imesh, lm, iatom, idir) / atoms%rmsh(imesh, itype) &
328 : ! & / atoms%rmsh(imesh, itype)
329 : ! end do
330 : ! end do
331 : ! end do
332 : ! end do
333 : ! end do
334 : !end do
335 :
336 0 : call convolMTgrVeff0dKern(atoms, grRho0MT, dKernMTGPts, gausWts, ylm, grVxc0MT)
337 :
338 0 : iatom = 0
339 0 : do itype = 1, atoms%ntype
340 0 : do ieqat = 1, atoms%neq(itype)
341 0 : iatom = iatom + 1
342 0 : do idir = 1, 3
343 : ! This is not really valid with the l + 1 and l+ 2 because the grVxc0MT is not calculated correctly with the gauß mesh
344 0 : do lm = 1, (atoms%lmax(itype) + 1)**2
345 0 : do imesh = 1, atoms%jri(itype)
346 0 : grVeff0MT(imesh, lm, idir, iatom) = grVeff0MT(imesh, lm, idir, iatom) + grVxc0MT(imesh, lm, iatom, idir)
347 : end do
348 : end do
349 : end do
350 : end do
351 : end do
352 : end if
353 : end if
354 :
355 0 : end subroutine genGrVeff0
356 :
357 : !--------------------------------------------------------------------------------------------------------------------------------------
358 : !> @author
359 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
360 : !>
361 : !> @brief
362 : !> Calculates the multipole moments of the gradient of the Hartree potential incorporating a gradient of the density as real
363 : !> charge.
364 : !>
365 : !> @param[in] atoms : Atoms type, see types.f90
366 : !> @param[in] cell : Unit cell type, see types.f90
367 : !> @param[in] ngdp : Number of G-vectors for potentials and densities
368 : !> @param[in] gdp : G-vectors of potentials and densities
369 : !> @param[in] grRho0IR : Plane-wave coefficients of the gradient of the interstitial unperturbed density
370 : !> @param[in] grRho0MT : Spherical harmonic coefficients of the gradient of the muffin-tin unperturbed density
371 : !> @param[out] qlmGrVc0 : Multipole coefficients of the gradient of the Hartree potential incorporating a volume integral over the
372 : !> gradient of the unperturbed density
373 : !--------------------------------------------------------------------------------------------------------------------------------------
374 0 : subroutine CalcQlmGrVh0Vol( atoms, cell, ngdp, gdp, grRho0IR, grRho0MT, qlmGrVc0 )
375 :
376 : use m_types, only : t_atoms, t_cell
377 : use m_intgr, only : intgr3LinIntp
378 : use m_sphbes, only : Sphbes
379 :
380 : implicit none
381 :
382 : ! Type Parameters
383 : type(t_atoms), intent(in) :: atoms
384 : type(t_cell), intent(in) :: cell
385 :
386 : ! Scalar Parameters
387 : integer, intent(in) :: ngdp
388 :
389 : ! Array Parameters
390 : integer, intent(in) :: gdp(:, :)
391 : complex, intent(in) :: grRho0IR(:, :)
392 : complex, intent(in) :: grRho0MT(:, :, :, :)
393 : complex, intent(out) :: qlmGrVc0(:, :, :)
394 :
395 : ! Local Variables:
396 : !
397 : ! iatom : runs over all atoms
398 : ! lm : encodes oqn_l and mqn_m
399 : ! itype : runs over all atom types
400 : ! idir : runs over 3 directions the atom can be displaced to
401 : ! tempGaunt1 : auxillary variable to store a Gaunt coefficient
402 : ! tempGaunt2 : auxillary variable to store a Gaunt coefficient
403 : ! imesh : runs over mesh points of current grid
404 : ! rl2 : stores R^(l + 2)
405 : ! fint : stores integral
406 : ! ll1 : auxillary variable to calculate lm
407 : ! mqn_m : magnetic quantum number m
408 : ! sk3r : stores argument of spherical Bessel function
409 : ! mqn_mpp : magnetic quantum number m", also used for indexing 3 directions the atom can be displaced to
410 : ! oqn_l : orbital quantum number l
411 : ! iGvec : indexes current G-vector
412 : ! gradrho0PWi : stores the the second part of equation 7.58 using equation 7.60 from PhD thesis Aaron Klüppelberg
413 : ! pylm : contains 4π i i^l G / |G| exp(i G τ) Y*_lm(G / |G|)
414 : ! sbes : stores the Bessel function
415 : ! cil : stores everything except for pylm of the result of this routine
416 : ! f : stores the integrand of the first term in (7.58, PhD thesis Aaron Klüppelberg)
417 :
418 : ! Local Scalar Variables
419 : integer :: iatom
420 : integer :: lm
421 : integer :: itype
422 : integer :: idir
423 : integer :: imesh
424 : real :: rl2
425 : real :: intgrResR
426 : real :: intgrResI
427 : real :: ll1
428 : integer :: mqn_m
429 : integer :: ieqat
430 : real :: sk3r
431 : integer :: oqn_l
432 : integer :: iG
433 : real :: normGext
434 : complex :: cil
435 : #ifdef DEBUG_MODE
436 : !todo review these variables
437 : real :: analyticalInt
438 : logical :: expensiveDebug
439 : #endif
440 :
441 : ! Local Array Variables
442 0 : complex, allocatable :: grRhoIRlm(:, :, :)
443 0 : complex, allocatable :: pylm(:, :)
444 0 : real, allocatable :: sbes(:)
445 0 : real, allocatable :: intgrR(:)
446 0 : real, allocatable :: intgrI(:)
447 : real :: Gext(3)
448 :
449 : #ifdef DEBUG_MODE
450 : !todo review these variables
451 : complex, allocatable :: gradrho0PWiTest(:, :, :)
452 : real, allocatable :: integrandTest(:)
453 : real, allocatable :: sbesIntegr(:, :, :)
454 : complex, allocatable :: pylmOld(:, :, :)
455 : #endif
456 :
457 : ! Initialize dummy assumed shape array
458 0 : allocate( intgrR(atoms%jmtd), intgrI(atoms%jmtd) )
459 0 : intgrR(:) = 0.
460 0 : intgrI(:) = 0.
461 0 : qlmGrVc0(:, :, :) = cmplx(0., 0.)
462 :
463 : ! Calculate the first term in 7.58 using 7.59 (in 7.59 there are mistakes, please refer to dissCRG).
464 0 : do idir = 1, 3
465 0 : iatom = 0
466 0 : do itype = 1, atoms%ntype
467 0 : do ieqat = 1, atoms%neq(itype)
468 0 : iatom = iatom + 1
469 0 : do oqn_l = 0, atoms%lmax(itype)
470 0 : do mqn_m = -oqn_l, oqn_l
471 0 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
472 0 : intgrR = 0
473 0 : intgrI = 0
474 0 : do imesh = 1, atoms%jri(itype)
475 0 : intgrR(imesh) = atoms%rmsh(imesh, itype)**(oqn_l + 2) * real( grRho0MT(imesh, lm, iatom, idir) )
476 0 : intgrI(imesh) = atoms%rmsh(imesh, itype)**(oqn_l + 2) * aimag( grRho0MT(imesh, lm, iatom, idir) )
477 : end do ! imesh
478 0 : call Intgr3LinIntp( intgrR(1:atoms%jri(itype)), atoms%rmsh(1:atoms%jri(itype), itype), atoms%dx(itype), atoms%jri(itype), intgrResR, 1 )
479 0 : call Intgr3LinIntp( intgrI(1:atoms%jri(itype)), atoms%rmsh(1:atoms%jri(itype), itype), atoms%dx(itype), atoms%jri(itype), intgrResI, 1 )
480 0 : qlmGrVc0(lm, iatom, idir) = cmplx( intgrResR, intgrResI )
481 : end do ! mqn_m
482 : end do ! oqn_l
483 : end do ! ieqat
484 : end do ! itype
485 : end do ! idir
486 :
487 : ! N_check: qlmGrVc0(lm,iatom,idir)=\int_{0}^{R_{MT^{\alpha}}}dr r^{l+2} grRho0MT(r, lm, iatom, idir)
488 : ! Consistent with Aaron (7.58) part 1 and (2.35) in CRG 16.12.2020-19:00
489 :
490 : ! This block calculates the second part of equation 7.58, where the integrand is rewritten using using equation 7.60 in PhDthesAK
491 : ! and then can be evaluated in similiar manner as in 7.57 PhDthesAK analytically.
492 : ! 4 π i^l i \sum_G G / |G| ρ^(0)(G) Y*_lm(G / |G|) exp(i G τ_β) R_β^(l + 2) j_(l + 1)(|G|R_β) .
493 :
494 0 : allocate( grRhoIRlm(( atoms%lmaxd + 1 )**2, atoms%nat, 3) )
495 0 : allocate( pylm(( atoms%lmaxd + 1 )**2, atoms%nat ) )
496 0 : allocate( sbes(0:atoms%lmaxd + 1) )
497 0 : grRhoIRlm(:, :, :) = cmplx(0., 0.)
498 0 : pylm(:, :) = cmplx(0., 0.)
499 0 : sbes(:) = 0.
500 :
501 0 : do idir = 1, 3
502 0 : do iG = 1, ngdp
503 : !call Phasy1nSymVeclp1( atoms, cell, gdp(1:3, iG), pylm )
504 : ! calculates 4 π i^l exp(i Gvec * taual) Y*_lm(Gvec / |Gvec|)
505 0 : pylm(:, :) = cmplx(0., 0.)
506 0 : call Phasy1nSym( atoms, cell, gdp(1:3, iG), [0., 0., 0.], pylm)
507 0 : Gext(1:3) = matmul( cell%bmat(1:3, 1:3), gdp(1:3, iG) )
508 0 : normGext = norm2( Gext(1:3) )
509 0 : if (normGext < 1e-12) cycle
510 0 : iatom = 0
511 0 : do itype = 1, atoms%ntype
512 0 : do ieqat = 1, atoms%neq(itype)
513 0 : iatom = iatom + 1
514 0 : sk3r = normGext * atoms%rmt(itype)
515 0 : sbes(:) = 0.
516 : ! calculates spherical bessel function with argument sk3r
517 0 : call Sphbes( atoms%lmax(itype) + 1, sk3r, sbes )
518 0 : rl2 = atoms%rmt(itype)**2
519 0 : do oqn_l = 0, atoms%lmax(itype)
520 0 : cil = sbes(oqn_l + 1) * rl2 * grRho0IR(iG, idir) / normGext
521 0 : ll1 = oqn_l * (oqn_l + 1) + 1
522 0 : do mqn_m = -oqn_l, oqn_l
523 0 : lm = ll1 + mqn_m
524 0 : grRhoIRlm(lm, iatom, idir) = grRhoIRlm(lm, iatom, idir) + cil * pylm(lm, iatom)
525 : end do ! m
526 0 : rl2 = rl2 * atoms%rmt(itype)
527 : end do ! l
528 : end do !ieqat
529 : end do ! itype
530 : end do ! iG
531 : end do ! idir
532 :
533 : ! N_check: grRhoIRlm(lm, iatom, idir)=\sum_{\bm{G}\neq\bm{0}} \frac{j_{l+1}(GR_{\alpha})R_{\alpha}^{l+2}}{G}grRho0IR (G, idir)*
534 : ! 4\pi i^{l} e^{i2\pi \bm{G}_{int}\cdot\bm{tau}_{\alpha,int}}Y_{lm}^{*}(\har\bm{G})
535 : ! Consistent with Aaron (7.58) part 2 and (2.35) in CRG 16.12.2020-19:00
536 :
537 : ! Unit the MT and the IR contribution of 7.58 in PhDthesAK
538 0 : do idir = 1, 3
539 0 : iatom = 0
540 0 : do itype = 1, atoms%ntype
541 0 : do ieqat = 1, atoms%neq(itype)
542 0 : iatom = iatom + 1
543 0 : do oqn_l = 0, atoms%lmax(itype)
544 0 : do mqn_m = -oqn_l, oqn_l
545 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
546 0 : qlmGrVc0(lm, iatom, idir) = qlmGrVc0(lm, iatom, idir) - grRhoIRlm(lm, iatom, idir)
547 : end do ! mqn_m
548 : end do ! oqn_l
549 : end do ! ieqat
550 : end do ! itype
551 : end do ! idir
552 :
553 : ! N_check:
554 : ! Consistent with (2.35) in CRG 16.12.2020-19:00
555 :
556 : if (.false.) then
557 : ! Debug output for maintenance
558 : do idir = 1, 3
559 : do oqn_l = 0, atoms%lmax(1)
560 : do mqn_m = -oqn_l, oqn_l
561 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
562 : write(2022, '(2(i8),2(f15.8))') lm, idir, qlmGrVc0(lm, 1, idir)
563 : end do ! mqn_m
564 : end do ! oqn_l
565 : end do ! idir
566 : end if
567 :
568 : !#ifdef DEBUG_MODE
569 : !
570 : ! ! This tests whether the numerical and the analytical integral of r^(l+2) j_l(GR) is the same
571 : ! !TODO not sure if this is working anymore, is this even important to test anymore, will not be formatted now?
572 : ! !TODO format this and review
573 : ! expensiveDebug = .false.
574 : ! if (expensiveDebug) then
575 : ! allocate( gradRho0PWiTest((atoms%lmaxd + 2)**2, atoms%nat, 3) )
576 : ! allocate( integrandTest(atoms%jmtd) )
577 : ! allocate( sbesIntegr(0:atoms%lmaxd + 1, ngdp, atoms%ntype) )
578 : ! gradRho0PWiTest = cmplx(0., 0.)
579 : ! integrandTest = 0.
580 : ! sbesIntegr(:, :, :) = 0.
581 : !
582 : ! write (*, '(a)') 'l, numerical integral of spherical lth Bessel function, analytical version and differences'
583 : ! do itype = 1, atoms%ntype
584 : ! write (*, '(a, i1)') 'atom type ', itype
585 : ! do iG = 1, ngdp
586 : ! Gext = matmul(cell%bmat, gdp(:, iG))
587 : ! normGext = norm2(Gext)
588 : ! if (normGext == 0) cycle
589 : ! do oqn_l = 0, atoms%lmax(itype) + 1
590 : ! integrandTest = 0
591 : ! do imesh = 1, atoms%jri(itype)
592 : ! sk3r = normGext * atoms%rmsh(imesh, itype)
593 : ! call sphbes(atoms%lmax(itype) + 2, sk3r, sbes)
594 : ! integrandTest(imesh) = sbes(oqn_l) * atoms%rmsh(imesh, itype)**(oqn_l + 2)
595 : ! end do ! imesh
596 : ! !call intgr3(integrandTest, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), intgrResR)
597 : ! call intgr3LinIntp(integrandTest, atoms%rmsh(:, itype), atoms%dx(itype), atoms%jri(itype), intgrResR, 1)
598 : ! sbesIntegr(oqn_l, iG, itype) = intgrResR
599 : ! analyticalInt = sbes(oqn_l + 1) * atoms%rmt(itype)**(oqn_l + 3) / atoms%rmt(itype)**(1) / normGext
600 : ! if (iG == 1) then
601 : ! if (abs(analyticalInt - intgrResR) > 1e-9) then
602 : ! write (*, '(i2,2x,3(es15.8,2x))') oqn_l, intgrResR, analyticalInt, abs(analyticalInt - intgrResR)
603 : ! end if
604 : ! end if
605 : ! end do
606 : ! end do
607 : ! end do
608 : !
609 : !
610 : ! !todo optimize loop structure
611 : ! allocate(pylmOld((atoms%lmaxd + 1)**2, atoms%nat, 3))
612 : ! do idir = 1, 3
613 : ! do iG = 1, ngdp
614 : ! iatom = 0
615 : ! call Phasy1nSymUVeclp1(atoms, cell, gdp(:, iG), pylmOld) ! todo is this method really necessary and not yet there?, review it! Can be replaced by normal phasy routine, by dividing through norm of G-vector
616 : ! do itype = 1, atoms%ntype
617 : ! do ieqat = 1, atoms%neq(itype)
618 : ! iatom = iatom + 1
619 : ! do oqn_l = 0, atoms%lmax(itype) + 1
620 : ! cil = sbesIntegr(oqn_l, iG, itype) * rho0IRpw(iG)
621 : ! ll1 = oqn_l * (oqn_l + 1) + 1
622 : ! do mqn_m = -oqn_l, oqn_l
623 : ! lm = ll1 + mqn_m
624 : ! gradRho0PWiTest(lm, iatom, idir) = gradRho0PWiTest(lm, iatom, idir) + cil * pylmOld(lm, iatom, idir)
625 : ! end do
626 : ! end do
627 : ! end do
628 : ! end do
629 : ! end do
630 : ! end do
631 : !
632 : ! deallocate(pylmOld)
633 : !
634 : ! write (*, '(a)') 'Analytical version and numerical end result and its difference is displayed if difference not numerically zero'
635 : ! do idir = 1, 3
636 : ! iatom = 0
637 : ! do itype = 1, atoms%ntype
638 : ! do ieqat = 1, atoms%neq(itype)
639 : ! iatom = iatom + 1
640 : ! do oqn_l = 0, atoms%lmax(itype) + 1
641 : ! do mqn_m = -oqn_l, oqn_l
642 : ! lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
643 : ! if ( abs(real(grRhoIRlm(lm, iatom, idir) - gradRho0PWiTest(lm, iatom, idir)) > 1e-9) .or. &
644 : ! &abs( aimag(grRhoIRlm(lm, iatom, idir) - gradRho0PWiTest(lm, iatom, idir)) > 1e-9) ) then
645 : ! ! write (*, *) 'idir= ', idir, 'iatom= ', iatom, 'oqn_l= ', oqn_l, 'mqn_m= ', mqn_m
646 : ! write (*, *) idir, oqn_l, mqn_m
647 : ! write (*, '(3(es15.8,2x))') real(grRhoIRlm(lm, iatom, idir)), real(gradRho0PWiTest(lm, iatom, idir)), &
648 : ! &abs(real(grRhoIRlm(lm, iatom, idir) - gradRho0PWiTest(lm, iatom, idir)))
649 : ! write (*, '(3(es15.8,2x))') aimag(grRhoIRlm(lm, iatom, idir)), aimag(gradRho0PWiTest(lm, iatom, idir)), &
650 : ! &aimag(grRhoIRlm(lm, iatom, idir) - gradRho0PWiTest(lm, iatom, idir))
651 : ! write (*, *)
652 : ! end if
653 : ! end do
654 : ! end do
655 : ! end do
656 : ! end do
657 : ! end do
658 : ! end if
659 : !
660 : !
661 : ! do idir = 1, 1
662 : ! iatom = 0
663 : ! do itype = 1, atoms%ntype
664 : ! do ieqat = 1, atoms%neq(itype)
665 : ! iatom = iatom + 1
666 : ! do oqn_l = 0, atoms%lmax(itype)
667 : ! do mqn_m = -oqn_l, oqn_l
668 : ! lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
669 : ! if ( ((abs(grRhoIRlm(lm, iatom, idir)) < 1e-9) .and. (abs(qlmGrVc0(lm, 1, idir)) > 1e-9))&
670 : ! .or. ((abs(grRhoIRlm(lm, iatom, idir)) > 1e-9) .and. (abs(qlmGrVc0(lm, 1, idir)) < 1e-9))) then
671 : ! ! lm channels where there is a MT part but no PW part
672 : ! write (*, '(a)') 'No muffin-tin and plane-wave part for lm channels:'
673 : ! write(*, '(i2, 2x, i3, 2x, 2(2(es12.5),2x))') oqn_l, mqn_m, grRhoIRlm(lm, iatom, idir), qlmGrVc0(lm, iatom, idir)
674 : ! end if
675 : ! end do
676 : ! end do
677 : ! end do
678 : ! end do
679 : ! end do
680 : !
681 : !#endif
682 :
683 0 : end subroutine CalcQlmGrVh0Vol
684 :
685 : !--------------------------------------------------------------------------------------------------------------------------------------
686 : !> @author
687 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
688 : !>
689 : !> @brief
690 : !> Calculates the pseudocharges of the gradient of the Hartree, external or Coulomb unperturbed potential.
691 : !>
692 : !> @details
693 : !> This routine generates the fourier coefficients of the pseudo charge density. It
694 : !> is very similar to psqpw.F90 but lacks any consideration of symmetry. It is based
695 : !> on equation 28, 30 in M.Weinert J.Math.Phys. 22(11) (1981) p.2434, where a factor
696 : !> 1 / R^l has been forgotten which can be seen starting from the unsimplified equation
697 : !> above equation 28 in the Weinert paper. As this routine has been modified to deliver
698 : !> the multipole moments to calculate interstitial contributions of potentials for the
699 : !> determination of the first variation of the Coulomb potential due to a phonon
700 : !> pertubation the result psq of the method is 3-dimensional to mirror the 3 directions
701 : !> the atoms can be displaced due to a phonon mode and the phonon q-vector is considered.
702 : !> The respective formula can be found in equation 7.63 in the PhD thesis of Aaron Klüppelberg,
703 : !> where equations 28 and 30 have been united.
704 : !>
705 : !> @param[in] atoms : Atoms type, see types.f90
706 : !> @param[in] cell : Unit cell type, see types.f90
707 : !> @param[in] dimensions : Dimension type, see types.f90
708 : !> @param[in] ngdp : Number of G-vectors for potentials and densities
709 : !> @param[in] harSw : Switch to enable Hartree contribution
710 : !> @param[in] extSw : Switch to enable external contribution
711 : !> @param[in] grRho0IR : Plane-wave coefficients of the gradient of the interstitial unperturbed density
712 : !> @param[in] gdp : G-vectors of potentials and densities
713 : !> @param[in] qlmGrVc0 : Multipole coefficients of the gradient of the Hartree, external or Coulomb unperturbed potential
714 : !> @param[out] psqGrVc0 : Plane-wave coefficients of the pseudocharge to construct the gradient of the Hartree, external or
715 : !> Coulomb potential depending on what multipole moments are passed in.
716 : !--------------------------------------------------------------------------------------------------------------------------------------
717 0 : subroutine PsqpwVeclp1(atoms, cell, ngdp, grRho0IR, harSw, extSw, gdp, qlmGrVc0, psqGrVc0)
718 :
719 : use m_types, only : t_atoms, t_cell
720 : use m_sphbes
721 :
722 : implicit none
723 :
724 : ! Variables:
725 : !
726 : ! atoms : atoms type defined in m_types
727 : ! cell : unit cell type defined in m_types
728 : ! dimensions : dimension type defined in m_types
729 : ! ngdp : number of G-vectors used for the potential and the density
730 : ! grRho0IR : planewave density in the interstitial region
731 : ! gdp : array of G-vectors used for the potential and density
732 : ! qlmGrVc0 : multipole moments to construct the pseudo charge
733 : ! psq : resulting 3-D pseudocharge of this routine
734 : ! itype : runs over all types
735 : ! rmtl : contains R^l for the prefactor
736 : ! oqn_l : orbital quantum number l
737 : ! p : auxiliary variable to calculate prefactor
738 : ! nc : loop variable to calculate prefactor
739 : ! fpo : auxiliary variable which contains 1 / Ω
740 : ! iGvec : runs over G-vectors which are used for the potential and the density
741 : ! idir : runs over directions the atoms can be displaced to
742 : ! iatom : runs over all atoms
743 : ! ncvn : contains n + l for a certain atom type, while n is taken from a table in the Weinert paper
744 : ! ieqat : runs over all equal atoms
745 : ! n1 : auxiliary variable which contains n - l + 1
746 : ! ll1 : auxiliary variable to calculate lm
747 : ! mqn_m : magnetic quantum number
748 : ! lm : variable encoding oqn_l and mqn_m
749 : ! pn : variable which stores prefactor with double faculties times 1 / R^l
750 : ! pylm : contains 4 π i^l exp((Gvec) * taual) Y*_lm(Gvec)
751 : ! Gext : G-vector in external coordinates
752 : ! sbes : contains spherical bessel function
753 : ! sl : auxiliary variable to calculate psq
754 : ! sm : auxiliary variable to calculate psq
755 : ! sa : auxiliary variable to calculate psq
756 :
757 : ! .. Scalar Arguments ..
758 : type(t_atoms), intent(in) :: atoms
759 : type(t_cell), intent(in) :: cell
760 :
761 : integer, intent(in) :: ngdp
762 : logical, intent(in) :: harSw
763 : logical, intent(in) :: extSw
764 :
765 : ! .. Array Arguments ..
766 : complex, intent(in) :: grRho0IR(:, :)
767 : integer, intent(in) :: gdp(:, :)
768 : complex, intent(in) :: qlmGrVc0(:, :, :)
769 : complex, intent(out) :: psqGrVc0(:, :)
770 :
771 : ! .. Local Scalars ..
772 : integer :: itype
773 : real :: rmtl
774 : integer :: oqn_l
775 : real :: p
776 : integer :: nc
777 : real :: fpo
778 : integer :: iG
779 : integer :: idir
780 : integer :: iatom
781 : integer :: ncvn
782 : integer :: ieqat
783 : integer :: n1
784 : integer :: ll1
785 : integer :: mqn_m
786 : integer :: lm
787 : complex :: sl ! sum over l
788 : complex :: sm ! sum over m
789 : complex :: sa ! sum over atoms
790 :
791 : ! .. Local Arrays ..
792 0 : real, allocatable :: pn(:, :)
793 0 : complex, allocatable :: pylm(:, :)
794 0 : real, allocatable :: sbes(:)
795 : real :: Gext(3)
796 :
797 : ! Init assumed shape array dummy
798 0 : psqGrVc0(:, :) = cmplx(0., 0.)
799 :
800 : ! This block calculates pn(l,itype) = (2l + 2nc(itype) + 3)!! / ((2l + 1)!! R^l) from equation 28 in the Weinert paper cited above or the
801 : ! second term of equation 7.63 (e.g.) in the PhD thesis of Aaron Klüppelberg. ncv(itype) is according to Weinert's paper still
802 : ! constant and defined as n + l and the formula to be calculated is simplified in the following
803 : ! way: Assume l + n = l', so the factorials cancel each other until l' = l - 1 since 3 is 1 + 2 and gives the next odd number.
804 : ! Since 2l + 2nc(itype) + 3 is odd we then multiply over all odd numbers from (2 l + 3) until (2ncv(itype) +3) which can be done
805 : ! by incrementing l until ncv(itype). That's why the loop starts at nc=l.
806 : ! The value of ncv per type is taken from Table 1 in the Weinert paper and given by the Fleur input .
807 0 : allocate( pn(0:atoms%lmaxd, atoms%ntype) )
808 0 : pn(:, :) = 0
809 :
810 0 : do itype = 1, atoms%ntype
811 0 : rmtl = 1.0
812 0 : do oqn_l = 0, atoms%lmax(itype)
813 0 : if (oqn_l < atoms%ncv(itype)) then ! N>=1
814 : p = 1.
815 0 : do nc = oqn_l, atoms%ncv(itype)
816 0 : p = p * (2 * nc + 3)
817 : end do ! nc
818 0 : pn(oqn_l, itype) = p / rmtl
819 : end if ! oqn_l < atoms%ncv(itype)
820 0 : rmtl = rmtl * atoms%rmt(itype)
821 : end do ! oqn_l
822 : end do ! itype
823 :
824 : ! N_check: valid.
825 :
826 : ! This block calculates the G /= 0 term in equation 28 of the Weinert paper cited above:
827 : ! \tilde \rho_s (K) = 4 π / Ω \sum_{lmi} (-i)^l (2l + 2nc(n) + 3)!! / ((2l + 1)!! R^l) j_{l + n + 1}(|G + q|R_i) / (|G + q| R_i)^{n-l+1}
828 : ! qlmGrVc0 \exp{-i (G + q) τ} Y_{lm}(G + q), which is the second term in equation 7.63. It is similiar to equation 28 in the Weinert
829 : ! paper cited above.
830 : ! The G = 0 term is not needed for calculating the interstitial contributions of the potentials and should be zero anyway.
831 0 : allocate( pylm((atoms%lmaxd + 1)**2, atoms%nat) )
832 0 : allocate( sbes(0:MAXVAL(atoms%ncv) + 1))
833 0 : pylm(:, :) = cmplx(0., 0.)
834 0 : sbes(:) = 0.
835 :
836 0 : fpo = 1. / cell%omtil
837 :
838 0 : do iG = 1, ngdp
839 0 : Gext(1:3) = matmul( cell%bmat(1:3, 1:3), real(gdp(1:3, iG)) )
840 0 : if ( norm2( Gext(1:3) ) < 1e-12 ) cycle
841 0 : pylm(:, :) = cmplx(0., 0.)
842 : ! Return 4 pi i^l exp(iG tau) Y_lm^{*}(G)
843 0 : call Phasy1nSym( atoms, cell, gdp(1:3, iG), [0., 0., 0.], pylm )
844 0 : do idir = 1, 3
845 : ! variable accumulating sum over l resetted
846 : sa = cmplx(0., 0.)
847 : iatom = 0
848 0 : do itype = 1, atoms%ntype
849 0 : ncvn = atoms%ncv(itype)
850 : ! We need the spherical Bessel functions up to ncvn + 1 + 1 because, ncvn only assumes lmax not lmax + 1
851 0 : sbes(:) = 0.
852 0 : call sphbes(ncvn + 1, norm2(Gext(1:3)) * atoms%rmt(itype), sbes)
853 0 : do ieqat = 1, atoms%neq(itype)
854 0 : iatom = iatom + 1
855 : ! variable accumulating sum over l resetted
856 0 : sl = cmplx(0., 0.)
857 0 : do oqn_l = 0, atoms%lmax(itype)
858 0 : if ( oqn_l >= ncvn ) cycle ! in earlier versions of FLEUR pn was set zero in this case, so zero was added to sl
859 0 : n1 = ncvn - oqn_l + 1
860 0 : ll1 = oqn_l * (oqn_l + 1) + 1
861 : ! variable accumulating sum over m resetted
862 0 : sm = cmplx(0., 0.)
863 0 : do mqn_m = -oqn_l, oqn_l
864 0 : lm = ll1 + mqn_m
865 0 : sm = sm + qlmGrVc0(lm, iatom, idir) * conjg( pylm(lm, iatom) )
866 : end do ! mqn_m
867 0 : sl = sl + sbes(ncvn + 1) * sm * (pn(oqn_l, itype) / ( ( norm2( Gext(1:3) ) * atoms%rmt(itype) )**n1 ))
868 : end do ! oqn_l
869 0 : sa = sa + sl
870 : end do ! ieqat
871 : end do ! itype
872 0 : if (harSw) then
873 0 : psqGrVc0(iG, idir) = grRho0IR(iG, idir) + fpo * sa
874 0 : else if (.not.harSw .and. extSw) then
875 0 : psqGrVc0(iG, idir) = fpo * sa
876 : end if ! special treatment if only external potential required
877 : end do ! idir
878 : end do ! iG
879 :
880 : ! N_check: psqGrVc0(iG, idir)=[grRho0IR(iG, idir) +] frac{1}{\Omega}\sum_{\alpha l m} qlmGrVc0(lm, iatom, idir)*
881 : ! 4\pi (-i)^{l} e^{-i2\pi \bm{G}_{int}\cdot\bm{tau}_{\alpha,int}}Y_{lm}(\hat{\bm{G})*
882 : ! \frac{j_{l+N+1}(|GR_{\alpha})}{(GR_{\alpha})^{N+1}}*
883 : ! \frac{(2l + 2N + 3)!!}{(2l + 1)!! R_{\alpha}^{l})}
884 : !
885 : ! Consistent with Aaron (7.64) and (2.40) in CRG 16.12.2020-19:00
886 :
887 0 : end subroutine psqpwVeclp1
888 :
889 : !--------------------------------------------------------------------------------------------------------------------------------------
890 : !> @author
891 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
892 : !>
893 : !> @brief
894 : !> Solves a boundary problem according to Weinert to gain the muffin-tin contribution of the Hartree, external or
895 : !> Coulomb potential.
896 : !>
897 : !> @param[in] atoms : Atoms type, see types.f90
898 : !> @param[in] cell : Unit cell type, see types.f90
899 : !> @param[in] ngdp : Number of G-vectors for potentials and densities
900 : !> @param[in] iatom : Index of atom in whose muffin-tin the calculation is to be performed
901 : !> @param[in] itype : Index of atom type to which atom belongs in whose muffin-tin the calculation is to be performed
902 : !> @param[in] harSw : Switch to enable Hartree contribution
903 : !> @param[in] extSw : Switch to enable external contribution
904 : !> @param[in] testGoldstein : Switch needed by a test using the Goldstone condition
905 : !> @param[in] grRhoTermSw : Switch to activate volume term contribution of muffin-tin boundary problem
906 : !> @param[in] grVc0IR : Plane-wave insterstitial coefficients of the gradient of the unperturbed effective potential
907 : !> @param[in] grRho0MT : Spherical harmonic coefficients of the gradient of the muffin-tin unperturbed density
908 : !> @param[in] gdp : G-vectors of potentials and densities
909 : !> @param[out] grVc0MT : Spherical harmonic muffin-tin coefficients of the gradient of the unperturbed effective potential
910 : !--------------------------------------------------------------------------------------------------------------------------------------
911 0 : subroutine vmtsCoul( atoms, cell, ngdp, iatom, itype, harSw, extSw, grVc0IR, grRho0MT, gdp, grVc0MT, testGoldstein, grRhoTermSw)
912 :
913 : use m_types_atoms
914 : use m_types_cell
915 : use m_intgr, only : intgr2LinIntp
916 : use m_sphbes
917 :
918 : implicit none
919 :
920 : ! Type Parameter
921 : type(t_atoms), intent(in) :: atoms
922 : type(t_cell), intent(in) :: cell
923 :
924 : ! Scalar Parameter
925 : integer, intent(in) :: ngdp
926 : integer, intent(in) :: iatom
927 : integer, intent(in) :: itype
928 : logical, intent(in) :: harSw
929 : logical, intent(in) :: extSw
930 : logical, intent(in) :: testGoldstein
931 : logical, intent(in) :: grRhoTermSw
932 :
933 : ! Array Parameter
934 : complex, intent(in) :: grVc0IR(:, :)
935 : complex, intent(in) :: grRho0MT(:,:,:,:)
936 : integer, intent(in) :: gdp(:, :)
937 : complex, intent(out) :: grVc0MT(:, :, :)
938 :
939 : !----------------------------------------------------------------------------------------------------------------------------------
940 : ! Local Scalar Variables
941 : ! oqn_l : orbital quantum number l
942 : ! mqn_m : magnetic quantum number m
943 : ! l21 : contains 2 l + 1
944 : ! fpl21 : contains prefactor of volume integral
945 : ! rmtl : auxiliary variable to calculate the volume integral
946 : ! rmt2l : auxiliary variable to calculate the volume integral
947 : ! rrlr : auxiliary variable to calculate the volume integral
948 : ! ror : auxiliary variable to calculate the volume integral
949 : ! lm : encodes orbital quantum number l and magnetic quantum number m
950 : ! iGvec : runs over all G-vectors used here
951 : ! idir : runs over all directions the atoms can be displaced to
952 : ! imesh : runs over mesh points of local sphere mesh
953 : !----------------------------------------------------------------------------------------------------------------------------------
954 : integer :: oqn_l
955 : integer :: mqn_m
956 : integer :: l21
957 : real :: fpl21
958 : real :: rmtl
959 : real :: rmt2l
960 : integer :: lm
961 : integer :: iG
962 : integer :: idir
963 : integer :: imesh
964 : real :: rrlr
965 : real :: ror
966 : real :: prefactor
967 :
968 : !----------------------------------------------------------------------------------------------------------------------------------
969 : ! Local Array Variables
970 : ! Gqext : cartesian representation of G + q
971 : ! sbf : contains spherical Bessel functions
972 : ! pylm : contains 4 π i^l exp(i (Gvec) * taual) Y*_lm(Gvec )
973 : ! vtl : contains the surface integral
974 : ! rrl : auxiliary variable to calculate the volume integral
975 : ! rrl1 : auxiliary variable to calculate the volume integral
976 : ! f1r : contains the integrand of the first type of integral occuring for the volume integral part
977 : ! f2r : contains the integrand of the second type of integral occuring for the volume integral part
978 : ! x1r : contains the result of the the first type of integral occuring for the volume integral part
979 : ! x2r : contains the result of the the first type of integral occuring for the volume integral part
980 : !todo comment imag
981 : !----------------------------------------------------------------------------------------------------------------------------------
982 : complex, allocatable :: pylm(:, :)
983 0 : real, allocatable :: sbf(:)
984 0 : complex, allocatable :: vtl(:, :)
985 0 : real, allocatable :: rrl(:)
986 0 : real, allocatable :: rrl1(:)
987 0 : complex, allocatable :: f1r(:)
988 0 : complex, allocatable :: f2r(:)
989 0 : real, allocatable :: f1rReal(:)
990 0 : real, allocatable :: f2rReal(:)
991 0 : real, allocatable :: f1rImag(:)
992 0 : real, allocatable :: f2rImag(:)
993 0 : real, allocatable :: x1rReal(:)
994 0 : real, allocatable :: x2rReal(:)
995 0 : real, allocatable :: x1rImag(:)
996 0 : real, allocatable :: x2rImag(:)
997 : real :: Gext(3)
998 :
999 : ! For a given iatom and itype the second term in equation (7.66 / 7.67, PhD thesis Aaron Klüppelberg) is evaluated, beginning from the sum
1000 0 : allocate( vtl((atoms%lmaxd + 1)**2, 3) )
1001 0 : allocate( sbf(0:atoms%lmaxd) )
1002 0 : allocate( pylm(( atoms%lmaxd + 1 )**2, atoms%nat) )
1003 0 : vtl(:, :) = cmplx(0., 0.)
1004 0 : sbf(:) = 0.
1005 0 : pylm(:, :) = cmplx(0., 0.)
1006 :
1007 : ! Init assumed shape dummy array
1008 0 : grVc0MT(:, :, :) = cmplx(0., 0.)
1009 :
1010 0 : do iG = 1, ngdp
1011 0 : Gext(1:3) = matmul( cell%bmat(1:3, 1:3), gdp(1:3, iG) )
1012 0 : if ( norm2( Gext(1:3) ) < 1e-12 ) cycle
1013 : ! Return 4 pi i^l exp(iG tau) Y_lm^{*}(G)
1014 0 : pylm(:, :) = cmplx(0., 0.)
1015 0 : call Phasy1nSym( atoms, cell, gdp(1:3, iG), [0., 0., 0.], pylm )
1016 0 : sbf(:) = 0.
1017 0 : call Sphbes( atoms%lmax(itype), norm2( Gext(1:3) ) * atoms%rmt(itype), sbf )
1018 0 : do idir = 1, 3
1019 0 : do oqn_l = 0, atoms%lmax(itype)
1020 0 : do mqn_m = -oqn_l, oqn_l
1021 0 : lm = oqn_l * (oqn_l + 1) + mqn_m + 1
1022 0 : vtl(lm, idir) = vtl(lm, idir) + grVc0IR(iG, idir) * sbf(oqn_l) * pylm(lm, iatom)
1023 : end do ! mqn_m
1024 : end do ! oqn_l
1025 : end do ! idir
1026 : end do ! iG
1027 0 : deallocate(pylm)
1028 0 : deallocate(sbf)
1029 :
1030 : ! This is not similar to the code beneath because we return.
1031 0 : if ( .not.grRhoTermSw ) then
1032 0 : do idir = 1, 3
1033 0 : do oqn_l = 0, atoms%lmax(itype)
1034 0 : rmtl = 1. / atoms%rmt(itype)**oqn_l ! 1 / R^l
1035 0 : do mqn_m = -oqn_l, oqn_l
1036 0 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
1037 0 : do imesh = 1, atoms%jri(itype)
1038 0 : ror = atoms%rmsh(imesh, itype)**oqn_l * rmtl ! r^l / R^l
1039 0 : grVc0MT(imesh, lm, idir) = ror * vtl(lm, idir)
1040 : end do
1041 : end do
1042 : end do
1043 : end do
1044 :
1045 : ! N_check: vtl(lm, idir)=\sum_{\bm{G}\new\bm{0}} grVc0IR(G, idir)*
1046 : ! 4\pi i^{l} e^{i2\pi (\bm{G}_{int}+\bm{q}_{int})\cdot\bm{tau}_{\alpha,int}}Y_{lm}^{*}(\hat{\bm{G}})*
1047 : ! j_{l}(G|R_{\alpha})+(\frac{r_{\alpha}}{\bm{R}_{\alpha}})^{l}
1048 : !
1049 : ! Consistent with Aaron (7.67) line 2 and (2.44) line 2/3 without lm-sum in CRG 16.12.2020-19:00
1050 :
1051 : ! Maintenance
1052 : if (.false.) then
1053 : do idir = 1, 3
1054 : do oqn_l = 0, atoms%lmax(itype)
1055 : rmtl = 1. / atoms%rmt(itype)**oqn_l ! 1 / R^l
1056 : do mqn_m = -oqn_l, oqn_l
1057 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
1058 : do imesh = 1, atoms%jri(itype)
1059 : write(2130, '(3(i8),2(f15.8))') imesh, lm, idir, grVc0MT(imesh, lm, idir)
1060 : end do
1061 : end do
1062 : end do
1063 : end do
1064 : end if
1065 : return
1066 : end if
1067 :
1068 : ! For a given iatom and itype the first term in equation (7.66 / 7.67, PhD thesis Aaron Klüppelberg) is calculated first and then
1069 : ! the second surface integral contribution inclucing the correct prefactor is added so that the complete Dirichlet boundary value
1070 : ! problem is solved.
1071 : ! Concerning the first term, basically, the integral is splitted into four integrals with respect to the definitions of r_> and r_<.
1072 : ! This gives two types of integrands x1r = r^(l + 2) * ρ(r) dr and x2r = r^(1 - l) * ρ(r) dr. Considering the bounds of the
1073 : ! integrals they can be rearranged to:
1074 : ! 4π / (2 l + 1) * (1 / r^(l + 1) \int_0^r s^(l + 2) ρ ds - r^l / R^(2 * l + 1) \int_0^R ρ s^(l + 2)
1075 : ! + r^l (\int_0^R 1 / s^(l + 1) ρ ds - \int_0^r 1 / s^(l + 1) ρ))
1076 :
1077 : ! Prefactor is for external contribution
1078 0 : prefactor = atoms%zatom(itype) * 4 * pi_const / 3
1079 0 : grVc0MT = cmplx(0., 0.)
1080 0 : if (harSw) then
1081 : allocate( rrl1(atoms%jmtd), rrl(atoms%jmtd), x1rReal(atoms%jmtd), x2rReal(atoms%jmtd), x1rImag(atoms%jmtd), &
1082 : & x2rImag(atoms%jmtd), f1rReal(atoms%jmtd), f1rImag(atoms%jmtd), f1r(atoms%jmtd), f2rReal(atoms%jmtd), &
1083 0 : & f2rImag(atoms%jmtd), f2r(atoms%jmtd) )
1084 : end if
1085 0 : do idir = 1, 3
1086 : ! grRhoTermSw = true implicitely
1087 0 : if (harSw) then
1088 0 : rrl1 = 0.
1089 0 : rrl = 0.
1090 0 : x1rReal = 0.
1091 0 : x2rReal = 0.
1092 0 : x1rImag = 0.
1093 0 : x2rImag = 0.
1094 0 : f1rReal = 0.
1095 0 : f1rImag = 0.
1096 0 : f1r = cmplx(0., 0.)
1097 0 : f2rReal = 0.
1098 0 : f2rImag = 0.
1099 0 : f2r = cmplx(0., 0.)
1100 0 : do oqn_l = 0, atoms%lmax(itype)
1101 0 : l21 = 2 * oqn_l + 1
1102 0 : fpl21 = fpi_const / l21 ! prefactor 1st term
1103 0 : rmtl = 1. / atoms%rmt(itype)**oqn_l ! 1 / R^l
1104 0 : rmt2l = 1. / atoms%rmt(itype)**l21 ! 1 / R^(2 l + 1)
1105 0 : do imesh = 1, atoms%jri(itype)
1106 0 : rrl(imesh) = atoms%rmsh(imesh, itype)**oqn_l ! r^l
1107 0 : rrl1(imesh) = 1. / (rrl(imesh) * atoms%rmsh(imesh, itype)) ! 1 / r^(l + 1)
1108 : end do ! imesh
1109 0 : do mqn_m = -oqn_l, oqn_l
1110 0 : lm = oqn_l * (oqn_l + 1) + mqn_m + 1
1111 0 : do imesh = 1, atoms%jri(itype)
1112 0 : x1rReal(imesh) = rrl(imesh) * real( grRho0MT(imesh, lm, iatom, idir) ) * atoms%rmsh(imesh, itype)**2
1113 0 : x2rReal(imesh) = rrl1(imesh) * real( grRho0MT(imesh, lm, iatom, idir) ) * atoms%rmsh(imesh, itype)**2
1114 0 : x1rImag(imesh) = rrl(imesh) * aimag( grRho0MT(imesh, lm, iatom, idir) )* atoms%rmsh(imesh, itype)**2
1115 0 : x2rImag(imesh) = rrl1(imesh) * aimag( grRho0MT(imesh, lm, iatom, idir) )* atoms%rmsh(imesh, itype)**2
1116 : end do ! mqn_m
1117 0 : call intgr2LinIntp(x1rReal, atoms%rmsh(1, itype), atoms%dx(itype), atoms%jri(itype), f1rReal)
1118 0 : call intgr2LinIntp(x2rReal, atoms%rmsh(1, itype), atoms%dx(itype), atoms%jri(itype), f2rReal)
1119 0 : call intgr2LinIntp(x1rImag, atoms%rmsh(1, itype), atoms%dx(itype), atoms%jri(itype), f1rImag)
1120 0 : call intgr2LinIntp(x2rImag, atoms%rmsh(1, itype), atoms%dx(itype), atoms%jri(itype), f2rImag)
1121 0 : x1rReal(:) = 0.
1122 0 : x2rReal(:) = 0.
1123 0 : x1rImag(:) = 0.
1124 0 : x2rImag(:) = 0.
1125 0 : f1r(1:atoms%jri(itype)) = cmplx( f1rReal(1:atoms%jri(itype)), f1rImag(1:atoms%jri(itype)) )
1126 0 : f2r(1:atoms%jri(itype)) = cmplx( f2rReal(1:atoms%jri(itype)), f2rImag(1:atoms%jri(itype)) )
1127 0 : f1rReal(:) = 0.
1128 0 : f2rReal(:) = 0.
1129 0 : f1rImag(:) = 0.
1130 0 : f2rImag(:) = 0.
1131 0 : do imesh = 1, atoms%jri(itype)
1132 0 : rrlr = rrl(imesh) * rmt2l ! r^l / R^(2 l + 1)
1133 0 : ror = rrl(imesh) * rmtl ! r^l / R^l
1134 : grVc0MT(imesh, lm, idir) = fpl21 * (rrl1(imesh) * f1r(imesh) - rrlr * f1r(atoms%jri(itype)) &
1135 0 : & + rrl(imesh) * (f2r(atoms%jri(itype)) - f2r(imesh))) + ror * vtl(lm, idir)
1136 : end do ! imesh
1137 0 : f1r = cmplx(0., 0.)
1138 0 : f2r = cmplx(0., 0.)
1139 : end do ! mqn_m
1140 0 : rrl = 0.
1141 0 : rrl1 = 0.
1142 : end do ! oqn_l
1143 : ! .and. grRhoTermSw implicitely
1144 0 : else if (.not.harSw .and. extSw) then
1145 : ! This part can be outsourced to a seperate routine
1146 0 : do oqn_l = 0, atoms%lmax(itype)
1147 0 : rmtl = 1. / atoms%rmt(itype)**oqn_l ! 1 / R^l
1148 0 : do mqn_m = -oqn_l, oqn_l
1149 0 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
1150 0 : do imesh = 1, atoms%jri(itype)
1151 0 : ror = atoms%rmsh(imesh, itype)**oqn_l * rmtl ! r^l / R^l
1152 0 : grVc0MT(imesh, lm, idir) = ror * vtl(lm, idir)
1153 : end do ! imesh
1154 : end do ! mqn_m
1155 : end do ! oqn_l
1156 :
1157 : ! N_check: grVc0MT(r, lm, idir)=\frac{4\pi}{2l+1}\int_{0}^{R_{\alpha}}ds_{\alpha} s_{\alpha}^{2}\frac{r_{<}^{l}}{r_{>}^{l+1}}
1158 : ! grRho0MT(r, lm, iatom, idir)*[1-(\frac{r_{>}}{\bm{R}_{\alpha}})^{2l+1}]+boundary term
1159 : !
1160 : ! Consistent with Aaron (7.67) line 1 and (2.44) line 1 in CRG 16.12.2020-19:00
1161 :
1162 : ! Maintenance
1163 : if (.false.) then
1164 : do oqn_l = 0, atoms%lmax(itype)
1165 : rmtl = 1. / atoms%rmt(itype)**oqn_l ! 1 / R^l
1166 : do mqn_m = -oqn_l, oqn_l
1167 : lm = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
1168 : do imesh = 1, atoms%jri(itype)
1169 : write(2020, '(3(i8),2(f15.8))') imesh, lm, idir, grVc0MT(imesh, lm, idir)
1170 : end do
1171 : end do
1172 : end do
1173 : end if ! maintenance
1174 :
1175 : end if
1176 : ! .and. grRhoTermSw implicitely
1177 : ! This has to stand here because it is valid for both ifs
1178 0 : if ( extSw ) then
1179 0 : do lm = 2, 4
1180 0 : do imesh = 1, atoms%jri(itype)
1181 : grVc0MT(imesh, lm, idir) = grVc0MT(imesh, lm, idir) + prefactor / atoms%rmsh(imesh, itype)**2 &
1182 0 : & * ( 1 - (atoms%rmsh(imesh, itype) / atoms%rmt(itype))**3) * 3 / 4 / pi_const * c_im(idir, lm - 1)
1183 : end do
1184 : end do
1185 : end if
1186 : end do
1187 :
1188 : ! N_check: grVc0MT(r, lm, idir)+=\frac{Z_{\alpha}4\pi}{3}r_{\alpha}^{-2}*
1189 : ! [1-(\frac{r_{\alpha}}{\bm{R}_{\alpha}})^{3}]\frac{3}{4\pi}c_im(idir, m)\delta_{l1}+boundary term
1190 : !
1191 : ! Consistent with Aaron (7.48) and (2.32)/(2.33) without lm-sum in CRG 16.12.2020-19:00
1192 :
1193 : ! Maintenance
1194 : if (.false.) then
1195 : do idir = 1, 3
1196 : do oqn_l = 0, atoms%lmax(itype)
1197 : do mqn_m = -oqn_l, oqn_l
1198 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
1199 : do imesh = 1, atoms%jri(itype)
1200 : write(2027, '(3(i8),2(f15.8))') imesh, lm, idir, grVc0MT(imesh, lm, idir)
1201 : end do ! imesh
1202 : end do ! mqn_m
1203 : end do ! oqn_l
1204 : end do ! idir
1205 : end if
1206 :
1207 0 : end subroutine vmtsCoul
1208 :
1209 : !--------------------------------------------------------------------------------------------------------------------------------------
1210 : !> @author
1211 : !> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
1212 : !>
1213 : !> @brief
1214 : !> Alternative way to calculate the external potential in the interstitial according to Equation 7.43 in PhdAK
1215 : !>
1216 : !> @note
1217 : !> This routine is out of order at the moment and has to be reviewed before being reactivated. But it has run once in the past.
1218 : !--------------------------------------------------------------------------------------------------------------------------------------
1219 0 : subroutine psqpwVecExt(atoms, cell, gdp, ngdp, psq)
1220 : ! *************************************************************************************
1221 : ! This routine calculates the Fourier coefficients of the pseudo-charge to gain the
1222 : ! gradient of the unperturbed external potential in the interstitial region. For reasons
1223 : ! of convenience, the sum over all atoms is already evaluated.
1224 : ! *************************************************************************************
1225 :
1226 : use m_sphbes
1227 : use m_types
1228 :
1229 : implicit none
1230 :
1231 : ! Variables:
1232 : !
1233 : ! atoms : atoms type defined in m_types
1234 : ! cell : unit cell type defined in m_types
1235 : ! dimens : dimension type defined in m_types
1236 : ! ngdp : number of G-vectors used for the potential and the density
1237 : ! gdp : array of G-vectors used for the potential and density
1238 : ! psq : resulting 3-D pseudocharge of this routine
1239 : ! itype : runs over all types
1240 : ! oqn_l : orbital quantum number l
1241 : ! df : loop variable to calculate prefactor
1242 : ! pfac : auxiliary variable which contains 1 / Ω
1243 : ! iGvec : runs over G-vectors which are used for the potential and the density
1244 : ! idirec : runs over directions the atoms can be displaced to
1245 : ! iatom : runs over all atoms
1246 : ! ncvn : contains n + l for a certain atom type, while n is taken from a table in the Weinert paper
1247 : ! ineq : runs over all equal atoms
1248 : ! pn : variable which stores prefactor with double faculties times 1 / R^l
1249 : ! Gext : G-vector in external coordinates
1250 : ! sbes : contains spherical bessel function
1251 :
1252 : ! .. Scalar Arguments ..
1253 : type(t_atoms), intent(in) :: atoms
1254 : type(t_cell), intent(in) :: cell
1255 : integer, intent(in) :: ngdp
1256 :
1257 : ! .. Array Arguments ..
1258 : integer, intent(in) :: gdp(:, :)
1259 : complex, intent(out) :: psq(:, :)
1260 :
1261 : ! .. Local Scalars ..
1262 : integer :: itype
1263 : integer :: oqn_l
1264 : integer :: df
1265 : complex :: pfac
1266 : integer :: iGvec
1267 : integer :: idirec
1268 : integer :: iatom
1269 : integer :: ncvn
1270 : integer :: ineq
1271 : integer :: ufacultb
1272 : integer :: nc
1273 :
1274 : ! .. Local Arrays ..
1275 0 : real :: pn(atoms%ntype)
1276 : real :: Gext(3)
1277 : real :: nGext
1278 0 : real :: sbes(0:atoms%lmaxd+MAXVAL(atoms%ncv)+1)
1279 :
1280 : ! This block calculates pn(itype) = (2 n + 5)!! in equation 7.43a PhD thesis of Aaron Klüppelberg.
1281 :
1282 : ! pn = 0
1283 : ! do itype = 1, atoms%ntype
1284 : ! oqn_l = 1
1285 : ! if (oqn_l < atoms%ncv(itype)) then
1286 : ! pn(itype) = 1.
1287 : ! do nc = oqn_l, atoms%ncv(itype)
1288 : ! pn(itype) = pn(itype) * (2 * nc + 3)
1289 : ! enddo
1290 : ! end if
1291 : ! enddo
1292 : ! write (*, *) pn * 3
1293 :
1294 0 : do itype = 1, atoms%ntype
1295 0 : oqn_l = 1
1296 0 : if (oqn_l >= atoms%ncv(itype)) then
1297 0 : pn(itype) = 0.0
1298 : else
1299 0 : pn(itype) = 1.
1300 0 : ufacultb = 2 * atoms%ncv(itype) + 3 ! + 2l is hidden in atoms%ncv, as this is n + l
1301 0 : do df = ufacultb, 1, -2
1302 0 : pn(itype) = pn(itype) * df
1303 : enddo
1304 : endif
1305 : enddo
1306 :
1307 :
1308 : ! This block calculates equation 7.43a in the PhD thesis of Aaron Klüppelberg:
1309 : ! There, the N is the bare n without the l in ncvn, l is already 1 and added to the 1 which was already there in l + n + 1
1310 : ! Thus, since ncvn is constant for all l per type we only add 1 instead of two to ncvn.
1311 0 : psq = cmplx(0., 0.)
1312 0 : do iGvec = 1, ngdp
1313 0 : Gext = matmul(cell%bmat, gdp(:, iGvec))
1314 0 : nGext = norm2(Gext)
1315 0 : if ( nGext == 0 ) cycle
1316 0 : do idirec = 1, 3
1317 : iatom = 1
1318 0 : do itype = 1, atoms%ntype
1319 0 : pfac = ImagUnit * cmplx(atoms%zatom(itype), 0) / cmplx(cell%omtil, 0)
1320 0 : ncvn = atoms%ncv(itype)
1321 0 : call sphbes(ncvn + 1, norm2(Gext) * atoms%rmt(itype), sbes) ! again one +1 is hidden in ncvn
1322 0 : do ineq = 1, atoms%neq(itype)
1323 : psq(iGvec, idirec) = psq(iGvec, idirec) + pfac * pn(itype)* sbes(ncvn + 1) / (nGext * atoms%rmt(itype))**(ncvn + 1)&
1324 0 : &* exp(-ImagUnit * tpi_const * dot_product(gdp(:, iGvec), atoms%taual(:, iatom))) * Gext(idirec)
1325 0 : iatom = iatom + 1
1326 : enddo
1327 : enddo
1328 : enddo
1329 : enddo
1330 :
1331 0 : end subroutine psqpwVecExt
1332 :
1333 0 : subroutine CalcQlmHarSurf( atoms, cell, iDtype, iDatom, ngdp, gdp, rho0IRpw, rho0MTsh, qlmHartSurf )
1334 :
1335 : use m_types, only : t_atoms, t_cell
1336 :
1337 : implicit none
1338 :
1339 :
1340 : ! Type parameters
1341 : type(t_atoms), intent(in) :: atoms
1342 : type(t_cell), intent(in) :: cell
1343 :
1344 : ! Scalar parameter
1345 : integer, intent(in) :: iDtype
1346 : integer, intent(in) :: iDatom
1347 : integer, intent(in) :: ngdp
1348 :
1349 : ! Array parameters
1350 : integer, intent(in) :: gdp(:, :)
1351 : complex, intent(in) :: rho0IRpw(:, :)
1352 : complex, intent(in) :: rho0MTsh(:, :, :, :)
1353 : complex, intent(out) :: qlmHartSurf(:, :)
1354 :
1355 : ! Scalar variables
1356 : integer :: oqn_l
1357 : integer :: mqn_m
1358 : integer :: lm
1359 : integer :: idir
1360 :
1361 : ! Array variables
1362 0 : complex, allocatable :: qlmHartSurfIR(:, :)
1363 0 : complex, allocatable :: qlmHartSurfMT(:, :)
1364 :
1365 0 : qlmHartSurf(:, :) = cmplx(0., 0.)
1366 :
1367 0 : allocate( qlmHartSurfIR(3, (atoms%lmax(iDtype) + 1)**2) )
1368 0 : allocate( qlmHartSurfMT(3, (atoms%lmax(iDtype) + 1)**2) )
1369 0 : qlmHartSurfIR(:, :) = cmplx(0., 0.)
1370 0 : qlmHartSurfMT(:, :) = cmplx(0., 0.)
1371 :
1372 0 : call CalcQlmHarSurfIR( atoms, cell, ngdp, iDtype, iDatom, gdp, rho0IRpw(:, 1), qlmHartSurfIR )
1373 0 : call CalcQlmHarSurMT(atoms, iDtype, iDatom, rho0MTsh, qlmHartSurfMT)
1374 :
1375 0 : do oqn_l = 0, atoms%lmax(iDtype)
1376 0 : do mqn_m = -oqn_l, oqn_l
1377 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
1378 0 : do idir = 1, 3
1379 0 : qlmHartSurf(lm, idir) = qlmHartSurfMT(idir, lm)
1380 0 : qlmHartSurf(lm, idir) = qlmHartSurf(lm, idir) - qlmHartSurfIR(idir, lm)
1381 0 : if (.false.) then
1382 : write(2401, '(3i8,2f15.8)') oqn_l, mqn_m, idir, qlmHartSurfIR(idir, lm)
1383 : write(2402, '(3i8,2f15.8)') oqn_l, mqn_m, idir, qlmHartSurfMT(idir, lm)
1384 : write(2403, '(3i8,2f15.8)') oqn_l, mqn_m, idir, qlmHartSurf(lm, idir)
1385 : end if
1386 : end do ! idir
1387 : end do ! mqn_m
1388 : end do ! oqn_l
1389 :
1390 0 : end subroutine CalcQlmHarSurf
1391 :
1392 0 : subroutine CalcQlmHarSurfIR( atoms, cell, ngdp, iDtype, iDatom, gdp, rho0IRpw, qlmHartSurfIR )
1393 :
1394 : use m_types, only : t_atoms, t_cell
1395 : use m_ylm, only : ylm4
1396 : use m_sphbes, only : sphbes
1397 : use m_gaunt, only : Gaunt1
1398 :
1399 : implicit none
1400 :
1401 : ! Type parameter
1402 : type(t_atoms), intent(in) :: atoms
1403 : type(t_cell), intent(in) :: cell
1404 :
1405 : ! Scalar parameter
1406 : integer, intent(in) :: ngdp
1407 : integer, intent(in) :: iDtype
1408 : integer, intent(in) :: iDatom
1409 :
1410 : ! Array parameter
1411 : integer, intent(in) :: gdp(:, :)
1412 : complex, intent(in) :: rho0IRpw(:)
1413 : complex, intent(out) :: qlmHartSurfIR(:, :)
1414 :
1415 : ! Scalar variables
1416 : integer :: idir
1417 : integer :: iG
1418 : integer :: oqn_l
1419 : integer :: mqn_m
1420 : integer :: oqn_l1p
1421 : integer :: mqn_m1p
1422 : integer :: mqn_m2p
1423 : integer :: lm
1424 : integer :: lm1p
1425 : complex :: phaseFac
1426 : complex :: pref
1427 : complex :: temp1
1428 : complex :: temp2
1429 : complex :: temp3
1430 : real :: gauntFactor
1431 :
1432 : ! Array variables
1433 0 : complex, allocatable :: ylm(:)
1434 0 : real, allocatable :: sbes(:)
1435 : real :: gExt(3)
1436 :
1437 : ! Init assumed-shape dummy array
1438 0 : qlmHartSurfIR(:, :) = cmplx(0., 0.)
1439 :
1440 0 : allocate( ylm( (atoms%lmax(iDtype) + 1)**2 ) )
1441 0 : allocate( sbes(0:atoms%lmax(iDtype)) )
1442 0 : ylm(:) = cmplx(0., 0.)
1443 0 : sbes(:) = 0.
1444 :
1445 0 : pref = fpi_const * atoms%rmt(iDtype) * atoms%rmt(iDtype)
1446 0 : do iG = 1, ngdp
1447 0 : gExt(1:3) = matmul( cell%bmat(1:3, 1:3), real(gdp(1:3, iG)) )
1448 :
1449 0 : ylm(:) = cmplx(0., 0.)
1450 0 : call ylm4( atoms%lmax(iDtype), gExt(1:3), ylm )
1451 :
1452 0 : sbes(:) = 0
1453 0 : call sphbes(atoms%lmax(iDtype), norm2(gExt(1:3)) * atoms%rmt(iDtype), sbes)
1454 :
1455 0 : phaseFac = exp( ImagUnit * tpi_const * dot_product(gdp(1:3, iG), atoms%taual(1:3, iDatom)))
1456 :
1457 0 : do oqn_l = 0, atoms%lmax(iDtype)
1458 0 : temp1 = pref * phaseFac * atoms%rmt(iDtype)**oqn_l * rho0IRpw(iG)
1459 0 : do mqn_m = -oqn_l, oqn_l
1460 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
1461 0 : do oqn_l1p = 0, atoms%lmax(iDtype)
1462 0 : temp2 = temp1 * sbes(oqn_l1p) * ImagUnit**oqn_l1p
1463 0 : do mqn_m1p = -oqn_l1p, oqn_l1p
1464 0 : lm1p = oqn_l1p * (oqn_l1p + 1) + 1 + mqn_m1p
1465 0 : temp3 = temp2 * conjg(ylm(lm1p))
1466 0 : do mqn_m2p = -1, 1
1467 0 : gauntFactor = Gaunt1( oqn_l, oqn_l1p, 1, mqn_m, mqn_m1p, mqn_m2p, atoms%lmax(iDtype))
1468 0 : do idir = 1, 3
1469 0 : qlmHartSurfIR(idir, lm) = qlmHartSurfIR(idir, lm) + c_im(idir, mqn_m2p + 2) * gauntFactor * temp3
1470 : end do ! idir
1471 : end do ! mqn_mpp
1472 : end do ! mqn_mp
1473 : end do ! oqn_lp
1474 : end do ! mqn_m
1475 : end do ! oqn_l
1476 : end do ! iG
1477 :
1478 0 : end subroutine CalcQlmHarSurfIR
1479 :
1480 0 : subroutine CalcQlmHarSurMT(atoms, iDtype, iDatom, rho0MTsh, qlmHartSurfMT)
1481 :
1482 : use m_types, only : t_atoms
1483 : use m_gaunt, only : Gaunt1
1484 :
1485 : implicit none
1486 :
1487 : ! Type parameters
1488 : type(t_atoms), intent(in) :: atoms
1489 :
1490 : ! Scalar parameter
1491 : integer, intent(in) :: iDtype
1492 : integer, intent(in) :: iDatom
1493 :
1494 : ! Array parameters
1495 : complex, intent(in) :: rho0MTsh(:, :, :, :)
1496 : complex, intent(out) :: qlmHartSurfMT(:, :)
1497 :
1498 : ! Scalar variables
1499 : integer :: oqn_l
1500 : integer :: oqn_l1p
1501 : integer :: mqn_m1p
1502 : integer :: mqn_m2p
1503 : integer :: mqn_m
1504 : integer :: lm
1505 : integer :: lm1p
1506 : integer :: idir
1507 : real :: gauntFactor
1508 :
1509 0 : qlmHartSurfMT(:, :) = cmplx(0., 0.)
1510 :
1511 0 : do oqn_l = 0, atoms%lmax(iDtype)
1512 0 : do mqn_m = -oqn_l, oqn_l
1513 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
1514 0 : do oqn_l1p = 0, atoms%lmax(iDtype)
1515 0 : do mqn_m1p = -oqn_l1p, oqn_l1p
1516 0 : lm1p = oqn_l1p * (oqn_l1p + 1) + 1 + mqn_m1p
1517 0 : do mqn_m2p = -1, 1
1518 0 : gauntFactor = Gaunt1( oqn_l, oqn_l1p, 1, mqn_m, mqn_m1p, mqn_m2p, atoms%lmax(iDtype))
1519 0 : do idir = 1, 3
1520 : qlmHartSurfMT(idir, lm) = qlmHartSurfMT(idir, lm) + c_im(idir, mqn_m2p + 2) * atoms%rmt(iDtype)**(oqn_l + 2) &
1521 0 : & * rho0MTsh(atoms%jri(iDtype), lm1p, iDatom, 1) * gauntFactor
1522 : end do ! idir
1523 : end do ! mqn_m2p
1524 : end do ! mqn_m1p
1525 : end do ! oqn_l1p
1526 : end do ! mqn_m
1527 : end do ! oqn_l
1528 :
1529 0 : end subroutine CalcQlmHarSurMT
1530 :
1531 0 : subroutine phasy1nSym(atoms, cell, Gvec, qptn, pylm)
1532 :
1533 : use m_ylm
1534 : use m_types_atoms
1535 : use m_types_cell
1536 :
1537 : implicit none
1538 :
1539 : ! Scalar Type Arguments
1540 : type(t_atoms), intent(in) :: atoms
1541 : type(t_cell), intent(in) :: cell
1542 :
1543 : ! Array Arguments
1544 : integer, intent(in) :: Gvec(:)
1545 : real, intent(in) :: qptn(:)
1546 : complex, intent(out) :: pylm(:, :)
1547 :
1548 : !-------------------------------------------------------------------------------------------------------------------------------
1549 : ! Local Scalar Variables
1550 : ! iatom : runs over all atoms
1551 : ! itype : runs over all types
1552 : ! ineq : runs over all equivalent atoms of one atom type
1553 : ! lm : encodes oqn_l and mqn_m
1554 : ! sf : stores exponential function
1555 : ! csf : stores exponential function times 4 π i^l
1556 : ! x : stores argument of exponential function
1557 : ! mqn_m : magnetic quantum number m
1558 : ! oqn_l : orbital quantum number l
1559 : ! ll1 : auxiliary variable to calculate lm
1560 : !-------------------------------------------------------------------------------------------------------------------------------
1561 : integer :: iatom
1562 : integer :: itype
1563 : integer :: ineq
1564 : integer :: lm
1565 : complex :: sf
1566 : complex :: csf
1567 : real :: x
1568 : integer :: mqn_m
1569 : integer :: oqn_l
1570 : integer :: ll1
1571 :
1572 : !-------------------------------------------------------------------------------------------------------------------------------
1573 : ! Local Array Variables
1574 : ! fpiul: stores 4 π i^l
1575 : ! Gqext: stores G + q in external coordinates
1576 : ! ylm : stores Y_lm
1577 : !-------------------------------------------------------------------------------------------------------------------------------
1578 0 : complex, allocatable :: fpiul(:)
1579 0 : complex, allocatable :: ylm(:)
1580 : real :: Gqext(3)
1581 :
1582 0 : allocate(fpiul(0:atoms%lmaxd))
1583 0 : fpiul(:) = cmplx(0., 0.)
1584 : ! calculates 4 π i^l resolved for every l, not divided by nop because no loop over symmetry operations
1585 0 : do oqn_l = 0, atoms%lmaxd
1586 0 : fpiul(oqn_l) = fpi_const * ImagUnit**oqn_l
1587 : enddo
1588 :
1589 :
1590 : ! calculates Y*_lm(\vec{G} + \vec{q}) for every l and m. The argument Gqext must be in external coordinates.
1591 0 : allocate(ylm((atoms%lmaxd + 1)**2))
1592 0 : ylm(:) = cmplx(0., 0.)
1593 0 : Gqext(1:3) = matmul(cell%bmat(1:3, 1:3), real(Gvec(1:3) + qptn(1:3)))
1594 0 : call ylm4(atoms%lmaxd, Gqext(1:3), ylm)
1595 0 : ylm = conjg(ylm)
1596 :
1597 :
1598 : ! calculates first exp(i (G + q) tau) and multiplies recent factors before storing the final result to pylm
1599 0 : iatom = 1
1600 0 : pylm(:, :) = cmplx(0.,0.)
1601 0 : do itype = 1, atoms%ntype
1602 0 : do ineq = 1, atoms%neq(itype)
1603 0 : x = tpi_const * dot_product(real(Gvec(1:3) + qptn(1:3)), atoms%taual(1:3, iatom))
1604 0 : sf = exp(ImagUnit * x)
1605 0 : do oqn_l = 0, atoms%lmax(itype)
1606 0 : ll1 = oqn_l * (oqn_l + 1) + 1
1607 0 : csf = fpiul(oqn_l) * sf
1608 0 : do mqn_m = -oqn_l, oqn_l
1609 0 : lm = ll1 + mqn_m
1610 0 : pylm(lm, iatom) = csf * ylm(lm)
1611 : enddo ! mqn_m
1612 : enddo ! oqn_l
1613 0 : iatom = iatom + 1
1614 : enddo ! ineq
1615 : enddo ! itype
1616 :
1617 0 : end subroutine phasy1nSym
1618 :
1619 0 : subroutine convolMTgrVeff0dKern(atoms, grRho0MT, dKernMTGPts, gWghts, ylm, grVxc0MT)
1620 :
1621 : use m_types_atoms
1622 :
1623 : implicit none
1624 :
1625 : ! Type parameter
1626 : type(t_atoms), intent(in) :: atoms
1627 :
1628 : ! Array parameter
1629 : complex, intent(in) :: grRho0MT(:, :, :, :)
1630 : complex, intent(in) :: ylm(:, :)
1631 : real, intent(in) :: gWghts(:) ! gaussian weights belonging to gausPts
1632 : real, intent(in) :: dKernMTGPts(:, :, :)
1633 : complex, intent(out) :: grVxc0MT(:, :, :, :)
1634 :
1635 : ! Local scalar variables
1636 : integer :: idir
1637 : integer :: iatom
1638 : integer :: itype
1639 : integer :: ieqat
1640 : integer :: igmesh ! Loop variable over sampling points of spherical Gauss mesh
1641 : integer :: irmesh ! Loop variable over radial MT mesh
1642 : integer :: oqn_l
1643 : integer :: lm_lonly !reduce multiplication when calculating lm
1644 : integer :: mqn_m
1645 : integer :: lm
1646 : complex :: grVxcMTKernAdd
1647 :
1648 : ! Local allocatable variables
1649 0 : complex, allocatable :: grRhoMTGpts(:, :) !grRhoMT on Gauss mesh
1650 0 : real, allocatable :: grVxcMTKernGPts(:, :)
1651 :
1652 0 : grVxc0MT(:, :, :, :) = cmplx(0., 0.)
1653 :
1654 0 : allocate( grRhoMTGpts(atoms%nsp(), atoms%jmtd), grVxcMTKernGPts(atoms%nsp(), atoms%jmtd) )
1655 0 : grRhoMTGpts(:, :) = cmplx(0., 0.)
1656 0 : grVxcMTKernGPts(:, :) = 0.
1657 :
1658 0 : do idir = 1, 3
1659 0 : iatom = 0
1660 0 : do itype = 1, atoms%ntype
1661 0 : do ieqat = 1, atoms%neq(itype)
1662 0 : iatom = iatom + 1
1663 0 : grRhoMTGpts(:, :) = cmplx(0., 0.)
1664 0 : grVxcMTKernGPts(:, :) = 0.
1665 : ! Density's gradient has l + 1 entries
1666 0 : do oqn_l = 0, atoms%lmax(itype)
1667 0 : lm_lonly = oqn_l * ( oqn_l + 1 ) + 1
1668 0 : do mqn_m = -oqn_l, oqn_l
1669 0 : lm = lm_lOnly + mqn_m
1670 : ! Evaluate grRho on spherical Gauss mesh in order to apply Gauss quadrature.
1671 0 : do irmesh = 1, atoms%jri(itype)
1672 0 : do igmesh = 1, atoms%nsp()
1673 0 : grRhoMTGpts(igmesh, irmesh) = grRhoMTGpts(igmesh, irmesh) + grRho0MT(irmesh, lm, iatom, idir) * ylm(igmesh, lm)
1674 : end do ! igmesh
1675 : end do ! irmesh
1676 : end do ! mqn_m
1677 : end do ! oqn_l
1678 : #if DEBUG_MODE
1679 : if ( any( aimag( grRhoMTGpts ) > 1e-7 ) ) then
1680 : write(*, *) 'Warning rhoMTGpts has imaginary components.'
1681 : end if
1682 : #endif
1683 : ! On the spherical Gauss mesh the integral reduces to a weighted (gWghts) sum (over all sampling points on the Gauss mesh)
1684 : ! of the MT exchange-correlation kernel and either the density's gradient or the first variation of the gradient.
1685 0 : do irmesh = 1, atoms%jri(itype)
1686 0 : do igmesh = 1, atoms%nsp()
1687 : grVxcMTKernGPts(igmesh, irmesh) = grVxcMTKernGpts(igmesh, irmesh) + real(grRhoMTGpts(igmesh, irmesh)) &
1688 0 : & * dKernMTGPts(igmesh, irmesh, iatom) * gWghts(igmesh)
1689 : end do
1690 : end do
1691 0 : do oqn_l = 0, atoms%lmax(itype)
1692 0 : lm_lonly = oqn_l * ( oqn_l + 1 ) + 1
1693 0 : do mqn_m = -oqn_l, oqn_l
1694 0 : lm = lm_lOnly + mqn_m
1695 0 : do irmesh = 1, atoms%jri(itype)
1696 : ! Back-transformation of the MT coefficients. Now they are expansion coefficients of the MT grid.
1697 0 : grVxcMTKernAdd = dot_product( grVxcMTKernGPts(:atoms%nsp(), irmesh), conjg(ylm(: atoms%nsp(), lm)) )
1698 : ! Add this contribution to MT exchange-correlation contribution to the potential
1699 0 : grVxc0MT(irmesh, lm, iatom, idir) = grVxc0MT(irmesh, lm, iatom, idir) + grVxcMTKernAdd
1700 : end do ! irmesh
1701 : end do ! mqn_m
1702 : end do ! oqn_l
1703 : end do ! ieqat
1704 : end do ! itype
1705 : end do ! itype
1706 :
1707 0 : end subroutine convolMTgrVeff0dKern
1708 :
1709 0 : subroutine convolGrRhoKern(stars, ngdp, ngpqdp, iDir, f1IR, f2IR, pdG2FouM, pdG2FouMv, f3IR, iqpt)
1710 :
1711 : use m_cfft
1712 : use m_types
1713 : use m_npy
1714 :
1715 : implicit none
1716 :
1717 : ! Type parameter
1718 : ! -------------
1719 : type(t_stars), intent(in) :: stars
1720 :
1721 : ! Scalar parameter
1722 : integer, intent(in) :: ngdp
1723 : integer, intent(in) :: ngpqdp
1724 : integer, intent(in) :: iDir
1725 : integer, intent(in) :: iqpt
1726 :
1727 : ! Array parameter
1728 : !----------------
1729 : ! 2 quantities to convolute
1730 : complex, intent(in) :: f1IR(:)
1731 : complex, intent(in) :: f2IR(:)
1732 : ! maps from a G to respective mesh entry of FFT mesh
1733 : integer, intent(in) :: pdG2FouM(:)
1734 : integer, intent(in) :: pdG2FouMv(:)
1735 : ! Convoluted quantity
1736 : complex, intent(out) :: f3IR(:)
1737 :
1738 : ! Maximal length of FFT mesh
1739 : integer :: ifftd
1740 : ! Loop index
1741 : integer :: iG
1742 : ! Scaling factor for canceling FFT artefacts.
1743 : real :: scaling
1744 : integer :: ii
1745 :
1746 : ! Real and imaginary parts of quantities to be convoluted (1, 2) and convoluted quantity 3
1747 0 : real, allocatable :: rf1(:)
1748 0 : real, allocatable :: if1(:)
1749 0 : real, allocatable :: rf2(:)
1750 0 : real, allocatable :: if2(:)
1751 0 : real, allocatable :: rf3(:)
1752 0 : real, allocatable :: if3(:)
1753 :
1754 :
1755 : ! Length of FFT mesh. The cartesian dimensions are stored sequentially. The G's expand maximally from -k_i to k_i (from this
1756 : ! cube a ball of radius gmax is cut off) so 3 * k_i should be in principle large enough for all given G's components. And we
1757 : ! have some free space to avoid aliasing. However, this factor of 3 is only working well for symmetrized quantities which is
1758 : ! why we have to expand the mx1, mx2, mx3 if we expand our quantities in plane waves.
1759 0 : ifftd = 27 * stars%mx1 * stars%mx2 * stars%mx3
1760 :
1761 : ! All quantities have the size of the FFT mesh.
1762 0 : allocate( rf1(0:ifftd - 1), if1(0:ifftd - 1), rf2(0:ifftd - 1), if2(0:ifftd - 1), rf3(0:ifftd - 1), if3(0:ifftd - 1) )
1763 :
1764 0 : rf1(:) = 0
1765 0 : rf2(:) = 0
1766 0 : rf3(:) = 0
1767 0 : if1(:) = 0
1768 0 : if2(:) = 0
1769 0 : if3(:) = 0
1770 :
1771 : ! Extract the real and imaginary part of the expansion coefficients of the quantities to be convoluted for every G-vector and
1772 : ! map them with pdG2FouM to its respective mesh point.
1773 0 : do iG = 1, ngpqdp
1774 0 : rf1(pdG2FouMv(iG)) = real(f1IR(iG))
1775 0 : if1(pdG2FouMv(iG)) = aimag(f1IR(iG))
1776 : end do
1777 :
1778 0 : do iG = 1, ngdp ! is it kimax - 1 or kimax?
1779 0 : rf2(pdG2FouM(iG)) = real(f2IR(iG))
1780 0 : if2(pdG2FouM(iG)) = aimag(f2IR(iG))
1781 : end do
1782 :
1783 : ! Complex FFT of 1st quantity into real space, this is done as it is done in FLEUR fft3d
1784 0 : call Cfft(rf1, if1, ifftd, 3 * stars%mx1, 3 * stars%mx1, 1)
1785 0 : call Cfft(rf1, if1, ifftd, 3 * stars%mx2, 9 * stars%mx1 * stars%mx2, 1)
1786 0 : call Cfft(rf1, if1, ifftd, 3 * stars%mx3, ifftd, 1)
1787 :
1788 : ! Complex FFT of 2nd quantity into real space
1789 0 : call Cfft(rf2, if2, ifftd, 3 * stars%mx1, 3 * stars%mx1, 1)
1790 0 : call Cfft(rf2, if2, ifftd, 3 * stars%mx2, 9 * stars%mx1 * stars%mx2, 1)
1791 0 : call Cfft(rf2, if2, ifftd, 3 * stars%mx3, ifftd, 1)
1792 :
1793 : ! Exploiting the convolution theorem
1794 0 : do ii = 0, ifftd - 1
1795 0 : rf3(ii) = real( cmplx(rf1(ii),if1(ii)) * cmplx(rf2(ii), if2(ii)))
1796 0 : if3(ii) = aimag( cmplx(rf1(ii),if1(ii)) * cmplx(rf2(ii), if2(ii)))
1797 : end do
1798 :
1799 : # ifdef DEBUG_MODE
1800 : ! It is sufficient to test it only for iqpt == 1 because the gradient of rho is calculated outside the q-loop and the density
1801 : ! variation has only to be real for q = 0
1802 : if(iqpt == 1) then
1803 : if ( any(abs(if3) > 1e-7 )) then
1804 : write(*, *) 'Convolution FFT has imaginary components'
1805 : !NOstopNO'Convolution FFT has imaginary components'
1806 : end if
1807 : if3 = 0
1808 : end if
1809 : #endif
1810 :
1811 : ! Complex FFT of convoluted quantity into reciprocal space
1812 0 : call Cfft(rf3, if3, ifftd, 3 * stars%mx1, 3 * stars%mx1, -1)
1813 0 : call Cfft(rf3, if3, ifftd, 3 * stars%mx2, 9 * stars%mx1 * stars%mx2, -1)
1814 0 : call Cfft(rf3, if3, ifftd, 3 * stars%mx3, ifftd, -1)
1815 :
1816 : ! We have to care for the artefacts of this FFT
1817 0 : scaling = 1. / ifftd
1818 0 : f3IR = cmplx( 0.0, 0.0 )
1819 : ! Map convoluted quantity from FFT mesh to plane-wave expansion coefficient representation.
1820 0 : do iG = 1, ngpqdp !kimax is max G-vector
1821 0 : f3IR(iG) = f3IR(iG) + cmplx(rf3(pdG2FouMv(iG)), if3(pdG2FouMv(iG))) * scaling
1822 : end do
1823 :
1824 0 : end subroutine convolGrRhoKern
1825 :
1826 0 : end module m_jpGrVeff0
|
