Line data Source code
1 : module m_mpmom
2 : ! ***********************************************************
3 : ! calculation of the multipole moments of the original charge
4 : ! density minus the interstitial charge density
5 : ! for each atom type
6 : !
7 : ! For yukawa_residual = .true. the subroutines calculate the
8 : ! multipole moments for the Yukawa potential instead of the
9 : ! Coulomb potential. This is used in the preconditioning of
10 : ! the SCF iteration for metallic systems.
11 : !
12 : ! qlmo(m,l,n) : mult.mom. of the muffin-tin charge density
13 : ! qlmp(m,l,n) : mult.mom. of the plane-wave charge density
14 : ! qlm (m,l,n) : (output) difference of the former quantities
15 : !
16 : ! references:
17 : ! for both the Coulomb and the Yukawa pseudo charge density:
18 : ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
19 : ! or see the original paper for the normal Coulomb case only:
20 : ! M. Weinert: J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15)
21 : ! ***********************************************************
22 :
23 : contains
24 :
25 1392 : subroutine mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell, qpw, rho, potdenType, qlm,l_coreCharge,&
26 1392 : & rhoimag, stars2, iDtype, iDir, rho0, qpw0, iDir2, mat2ord )
27 :
28 : use m_types
29 : USE m_constants
30 : implicit none
31 :
32 : type(t_input), intent(in) :: input
33 : type(t_mpi), intent(in) :: fmpi
34 :
35 : type(t_sym), intent(in) :: sym
36 : type(t_juphon), intent(in) :: juphon
37 : type(t_stars), intent(in) :: stars
38 : type(t_cell), intent(in) :: cell
39 : type(t_sphhar), intent(in) :: sphhar
40 : type(t_atoms), intent(in) :: atoms
41 : real, intent(in) :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
42 : complex, intent(in) :: qpw(:) !(stars%ng3)
43 : integer, intent(in) :: potdenType
44 : complex, intent(out) :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
45 : LOGICAL, OPTIONAL, INTENT(IN) :: l_coreCharge
46 :
47 : REAL, OPTIONAL, INTENT(IN) :: rhoimag(:,0:,:), rho0(:,0:,:)
48 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir ! DFPT: Type and direction of displaced atom
49 : COMPLEX, OPTIONAL, INTENT(IN) :: qpw0(:)
50 : TYPE(t_stars), OPTIONAL, INTENT(IN) :: stars2
51 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
52 : COMPLEX, OPTIONAL, INTENT(IN) :: mat2ord(5,3,3)
53 :
54 : integer :: j, jm, lh, mb, mem, mems, n, nd, l, nat, m
55 696 : complex :: qlmo(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
56 696 : complex :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
57 696 : complex :: qlmp_SF(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
58 :
59 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
60 :
61 696 : l_dfptvgen = PRESENT(stars2)
62 :
63 : ! multipole moments of original charge density
64 696 : if ( fmpi%irank == 0 ) then
65 114949 : qlmo = 0.0
66 : ! call mt_moments( input, atoms, sphhar, rho(:,:,:), potdenType, qlmo )
67 348 : IF (.NOT.l_dfptvgen) THEN
68 348 : call mt_moments( input, atoms, sym, juphon, sphhar, rho(:,:,:), potdenType,qlmo,l_coreCharge)
69 0 : ELSE IF (.NOT.PRESENT(iDir2)) THEN
70 : ! qlmo for the real part of rho1:
71 0 : call mt_moments( input, atoms, sym, juphon, sphhar, rho(:,:,:), potdenType,qlmo,l_coreCharge=.FALSE.)
72 : ! qlmo for the imaginary part of rho1 and the perturbation of vExt in the displaced atom:
73 0 : call mt_moments( input, atoms, sym, juphon, sphhar, rhoimag(:,:,:), potdenType,qlmo,l_coreCharge=.TRUE.,l_rhoimag=.TRUE.,iDtype=iDtype,iDir=iDir)
74 0 : IF (juphon%l_phonon) CALL dfpt_mt_moments_SF(atoms, sym, sphhar, iDtype, iDir, rho0(:,:,:), qlmo)
75 : ELSE
76 0 : call mt_moments( input, atoms, sym, juphon, sphhar, rho(:,:,:), potdenType,qlmo,l_coreCharge=.TRUE.,l_rhoimag=.FALSE.,iDtype=iDtype,iDir=iDir,iDir2=iDir2,mat2ord=mat2ord)
77 : END IF
78 : end if
79 :
80 : ! multipole moments of the interstitial charge density in the spheres
81 696 : call pw_moments( input, fmpi, stars, atoms, cell, sym, qpw(:), potdenType, qlmp , l_dfptvgen)
82 :
83 696 : IF (l_dfptvgen.AND..NOT.PRESENT(iDir2).AND.juphon%l_phonon) THEN
84 0 : CALL dfpt_pw_moments_SF( fmpi, stars2, atoms, cell, sym, iDtype, iDir, qpw0(:), qlmp_SF )
85 0 : qlmp = qlmp + qlmp_SF
86 : END IF
87 :
88 696 : if ( fmpi%irank == 0 ) then
89 : ! see (A14)
90 114949 : qlm = qlmo - qlmp
91 : ! output section
92 975 : do n = 1, atoms%ntype
93 627 : nat = atoms%firstAtom(n)
94 627 : write(oUnit, fmt=8000 ) n
95 6770 : do l = 0, atoms%lmax(n)
96 61025 : do m = -l, l
97 60398 : if ( qlmo(m,l,n)/=CMPLX(0.0) .or. qlmp(m,l,n)/=CMPLX(0.0) ) then
98 50668 : write(oUnit, fmt=8010 ) l, m, qlmo(m,l,n), qlmp(m,l,n)
99 : end if
100 : end do
101 : end do
102 : end do
103 :
104 : 8000 FORMAT (/,10x,'multipole moments for atom type=',i5,/,/,t3,'l',t7,&
105 : & 'm',t27,'original',t57,'plane wave')
106 : 8010 FORMAT (1x,i2,2x,i3,2x,2 (5x,2e15.5))
107 : !
108 : end if ! fmpi%irank == 0
109 :
110 696 : end subroutine mpmom
111 :
112 :
113 : ! subroutine mt_moments( input, atoms, sphhar, rho, potdenType, qlmo )
114 410 : subroutine mt_moments( input, atoms, sym, juphon, sphhar, rho, potdenType,qlmo,l_coreCharge,l_rhoimag,iDtype,iDir,iDir2,mat2ord)
115 : ! multipole moments of original charge density
116 : ! see (A15) (Coulomb case) or (A17) (Yukawa case)
117 :
118 : use m_intgr, only: intgr3
119 : use m_constants, only: sfp_const, POTDEN_TYPE_POTYUK, POTDEN_TYPE_CRYSTALFIELD
120 : use m_types
121 : use m_DoubleFactorial
122 : use m_SphBessel
123 : use m_juDFT
124 : implicit none
125 :
126 : type(t_input), intent(in) :: input
127 : type(t_sphhar), intent(in) :: sphhar
128 : type(t_atoms), intent(in) :: atoms
129 : type(t_sym), intent(in) :: sym
130 : type(t_juphon), intent(in) :: juphon
131 : real, intent(in) :: rho(: ,0:, :)
132 : integer, intent(in) :: potdenType
133 : complex, intent(inout) :: qlmo(-atoms%lmaxd:,0:,:)
134 : LOGICAL, OPTIONAL, INTENT(IN) :: l_coreCharge,l_rhoimag
135 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir ! DFPT: Type and direction of displaced atom
136 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
137 : COMPLEX, OPTIONAL, INTENT(IN) :: mat2ord(5,3,3)
138 :
139 : integer :: n, ns, jm, nl, l, j, mb, m, nat, i, imax, lmax
140 : real :: fint
141 1107 : real :: f( maxval( atoms%jri ) )
142 410 : real, allocatable, dimension(:,:) :: il, kl
143 : LOGICAL :: l_subtractCoreCharge
144 :
145 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
146 :
147 410 : l_dfptvgen = PRESENT(iDtype)
148 :
149 410 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
150 18 : allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
151 : end if
152 :
153 410 : l_subtractCoreCharge = .TRUE.
154 410 : if ( potdenType == POTDEN_TYPE_POTYUK ) l_subtractCoreCharge = .FALSE.
155 410 : IF(PRESENT(l_coreCharge)) l_subtractCoreCharge = l_coreCharge
156 :
157 410 : nat = 1
158 1107 : do n = 1, atoms%ntype
159 697 : ns = sym%ntypsy(nat)
160 697 : jm = atoms%jri(n)
161 697 : imax = atoms%jri(n)
162 697 : lmax = sphhar%llh(sphhar%nlh(ns), ns)
163 697 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
164 : !do concurrent (i = 1:imax)
165 11370 : do i = 1,imax
166 11370 : call ModSphBessel( il(:, i), kl(:, i), input%preconditioning_param * atoms%rmsh(i, n), lmax )
167 : end do
168 : end if
169 22198 : do nl = 0, sphhar%nlh(ns)
170 21501 : l = sphhar%llh(nl,ns)
171 21501 : if(jm < 2) call juDFT_error("This would be uninit in intgr3.")
172 15313898 : do j = 1, jm
173 15313898 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
174 15135698 : f(j) = atoms%rmsh(j,n) ** l * rho(j,nl,n)
175 : else
176 156699 : f(j) = il(l, j) * rho(j,nl,n)
177 : end if
178 : end do
179 21501 : call intgr3( f, atoms%rmsh(:,n), atoms%dx(n), jm, fint )
180 21501 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
181 207 : fint = fint * DoubleFactorial( l ) / input%preconditioning_param ** l
182 : end if
183 64942 : do mb = 1, sphhar%nmem(nl,ns)
184 42744 : m = sphhar%mlh(mb,nl,ns)
185 64245 : IF (.NOT.PRESENT(l_rhoimag)) THEN
186 42744 : qlmo(m,l,n) = qlmo(m,l,n) + sphhar%clnu(mb,nl,ns) * fint
187 : ELSE
188 0 : qlmo(m,l,n) = qlmo(m,l,n) + ImagUnit*sphhar%clnu(mb,nl,ns) * fint
189 : END IF
190 : end do
191 : end do
192 : ! if ( potdenType /= POTDEN_TYPE_POTYUK ) then
193 697 : if (l_subtractCoreCharge) then
194 612 : IF (.NOT.l_dfptvgen) THEN
195 612 : qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n) / sfp_const
196 0 : ELSE IF (.NOT.PRESENT(iDir2).AND.juphon%l_phonon) THEN
197 0 : IF ((n.EQ.iDtype)) THEN
198 0 : qlmo(-1:1,1,n) = qlmo(-1:1,1,n) - 3.0 / fpi_const * atoms%zatom(n) * c_im(iDir, :)
199 0 : ELSE IF ((0.EQ.iDtype)) THEN
200 0 : qlmo(-1:1,1,n) = qlmo(-1:1,1,n) + 3.0 / fpi_const * atoms%zatom(n) * c_im(iDir, :)
201 : END IF
202 0 : ELSE IF (juphon%l_phonon) THEN
203 : !IF ((n.EQ.iDtype).OR.(0.EQ.iDtype)) qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n) * (-0.2660214309643778) ! TODO: What the hell is this value???
204 0 : IF ((n.EQ.iDtype).OR.(0.EQ.iDtype)) qlmo(-2:2,2,n) = qlmo(-2:2,2,n) - 5.0 / fpi_const * atoms%zatom(n) * mat2ord(:,iDir2,iDir)
205 : END IF
206 : end if
207 1107 : nat = nat + atoms%neq(n)
208 : end do
209 :
210 410 : end subroutine mt_moments
211 :
212 :
213 : ! subroutine pw_moments( input, fmpi, stars, atoms, cell, sym, qpw, potdenType, qlmp_out )
214 696 : subroutine pw_moments( input, fmpi, stars, atoms, cell, sym, qpw_in, potdenType, qlmp_out, l_dfptvgen )
215 : ! multipole moments of the interstitial charge in the spheres
216 :
217 : use m_mpi_bc_tool
218 : use m_mpi_reduce_tool
219 : use m_phasy1
220 : use m_sphbes
221 :
222 : use m_constants, only: sfp_const, POTDEN_TYPE_POTYUK
223 : use m_types
224 : use m_DoubleFactorial
225 : use m_SphBessel
226 : implicit none
227 :
228 : type(t_input), intent(in) :: input
229 : type(t_mpi), intent(in) :: fmpi
230 :
231 : type(t_sym), intent(in) :: sym
232 : type(t_stars), intent(in) :: stars
233 : type(t_cell), intent(in) :: cell
234 : type(t_atoms), intent(in) :: atoms
235 : complex, intent(in) :: qpw_in(:)
236 : integer, intent(in) :: potdenType
237 : complex, intent(out) :: qlmp_out(-atoms%lmaxd:,0:,:)
238 : LOGICAL, INTENT(IN) :: l_dfptvgen
239 :
240 : integer :: n, k, l, ll1, lm, ierr, m
241 : integer :: maxBunchSize, numBunches, iBunch, firstStar, iTempArray
242 : complex :: sk3i, cil, nqpw, temp
243 1950 : complex :: pylm(( maxval( atoms%lmax ) + 1 ) ** 2, atoms%ntype)
244 : real :: sk3r, rl2
245 1950 : real :: aj(0:maxval( atoms%lmax ) + 1 )
246 696 : complex, ALLOCATABLE :: qpw(:)
247 1950 : real :: il(0:maxval( atoms%lmax ) + 1 )
248 1950 : real :: kl(0:maxval( atoms%lmax ) + 1 )
249 696 : complex :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
250 696 : complex, ALLOCATABLE :: qlmpStars(:,:,:,:)
251 : TYPE(t_parallelLoop) :: mpiLoop, ompLoop
252 :
253 2088 : ALLOCATE(qpw(stars%ng3))
254 1563080 : qpw = qpw_in(:stars%ng3)
255 229898 : qlmp = 0.0
256 :
257 696 : call mpi_bc(qpw,0,fmpi%mpi_comm)
258 :
259 2784 : firstStar = MERGE(2,1,norm2(stars%center)<=1e-8)
260 696 : maxBunchSize = 2*getNumberOfThreads() ! This bunch size is kind of a magic number detrmined from some
261 : ! naive performance tests for a 64 atom unit cell
262 696 : CALL calcNumberComputationBunches(firstStar, stars%ng3, maxBunchSize, numBunches)
263 :
264 4176 : ALLOCATE(qlmpStars(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype,maxBunchSize))
265 920288 : qlmpStars = CMPLX(0.0,0.0)
266 :
267 696 : CALL mpiLoop%init(fmpi%irank,fmpi%isize,0,numBunches-1)
268 195887 : DO iBunch = mpiLoop%bunchMinIndex, mpiLoop%bunchMaxIndex
269 195191 : CALL ompLoop%init(iBunch,numBunches,firstStar,stars%ng3)
270 : !$OMP parallel do default( NONE ) &
271 : !$OMP SHARED(input,atoms,stars,sym,cell,ompLoop,qpw,qlmpStars,potdenType) &
272 195887 : !$OMP private(iTempArray,pylm,nqpw,n,sk3r,aj,rl2,il,kl,sk3i,l,cil,ll1,m,lm)
273 : do k = ompLoop%bunchMinIndex, ompLoop%bunchMaxIndex
274 : iTempArray = k - ompLoop%bunchMinIndex + 1
275 : call phasy1( atoms, stars, sym, cell, k, pylm )
276 : nqpw = qpw(k) * stars%nstr(k)
277 : do n = 1, atoms%ntype
278 : sk3r = stars%sk3(k) * atoms%rmt(n)
279 : call sphbes( atoms%lmax(n) + 1, sk3r, aj )
280 : rl2 = atoms%rmt(n) ** 2
281 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
282 : call ModSphBessel( il(0:atoms%lmax(n)+1), kl(0:atoms%lmax(n)+1), input%preconditioning_param * atoms%rmt(n), atoms%lmax(n) + 1 )
283 : sk3i = nqpw / ( stars%sk3(k) ** 2 + input%preconditioning_param ** 2 ) * rl2
284 : else
285 : sk3i = nqpw / stars%sk3(k)
286 : end if
287 : do l = 0, atoms%lmax(n)
288 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
289 : cil = ( stars%sk3(k) * il(l) * aj(l+1) + input%preconditioning_param * il(l+1) * aj(l) ) * ( DoubleFactorial( l ) / input%preconditioning_param ** l ) * sk3i
290 : else
291 : cil = aj(l+1) * sk3i * rl2
292 : rl2 = rl2 * atoms%rmt(n)
293 : end if
294 : ll1 = l * ( l + 1 ) + 1
295 : do m = -l, l
296 : lm = ll1 + m
297 : qlmpStars(m,l,n,iTempArray) = qlmpStars(m,l,n,iTempArray) + cil * pylm(lm,n)
298 : end do
299 : end do ! l = 0, atoms%lmax(n)
300 : end do ! n = 1, atoms%ntype
301 : end do
302 : !$OMP END PARALLEL DO
303 : END DO
304 : !$OMP parallel do default( NONE ) &
305 : !$OMP SHARED(atoms,maxBunchSize,qlmpStars,qlmp) &
306 696 : !$OMP private(l,m,temp,iTempArray)
307 : DO n = 1, atoms%ntype
308 : DO l = 0, atoms%lmax(n)
309 : DO m = -l, l
310 : temp = CMPLX(0.0,0.0)
311 : DO iTempArray = 1, maxBunchSize
312 : temp = temp + qlmpStars(m,l,n,iTempArray)
313 : END DO
314 : qlmp(m,l,n) = temp
315 : END DO
316 : END DO
317 : END DO
318 : !$OMP END PARALLEL DO
319 :
320 696 : if ( fmpi%irank == 0 .AND. .NOT. l_dfptvgen) then
321 : ! q=0 term: see (A19) (Coulomb case) or (A20) (Yukawa case)
322 975 : do n = 1, atoms%ntype
323 975 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
324 612 : qlmp(0,0,n) = qlmp(0,0,n) + qpw(1) * stars%nstr(1) * atoms%volmts(n) / sfp_const
325 : else
326 15 : call ModSphBessel( il(0:1), kl(0:1), input%preconditioning_param * atoms%rmt(n), 1 )
327 15 : qlmp(0,0,n) = qlmp(0,0,n) + qpw(1) * stars%nstr(1) * sfp_const * atoms%rmt(n) ** 2 * il(1) / input%preconditioning_param
328 : end if
329 : end do
330 : end if
331 :
332 696 : CALL mpi_sum_reduce(qlmp, qlmp_out, fmpi%mpi_comm)
333 :
334 1392 : end subroutine pw_moments
335 :
336 0 : SUBROUTINE dfpt_mt_moments_SF(atoms, sym, sphhar, iDtype, iDir, rho0, qlmo)
337 : USE m_types
338 : USE m_gaunt, only : Gaunt1
339 : USE m_constants
340 :
341 : IMPLICIT NONE
342 :
343 : TYPE(t_atoms), INTENT(IN) :: atoms
344 : TYPE(t_sym), INTENT(IN) :: sym
345 : TYPE(t_sphhar), INTENT(IN) :: sphhar
346 : INTEGER, INTENT(IN) :: iDtype, iDir
347 : REAL, INTENT(IN) :: rho0(:, 0:, :)
348 : COMPLEX, INTENT(INOUT) :: qlmo(-atoms%lmaxd:,0:, :)
349 :
350 : INTEGER :: mb, n, nat, nl, ns, jm, l, lp, m, mp, mVec, pref
351 : REAL :: fint, gauntFactor
352 :
353 0 : pref = -1
354 :
355 0 : IF (iDtype.NE.0) THEN
356 0 : nat = iDtype
357 0 : pref = 1
358 : END IF
359 :
360 0 : DO n = MERGE(1,iDtype,iDtype.EQ.0), MERGE(atoms%ntype,iDtype,iDtype.EQ.0)
361 0 : nat = atoms%firstAtom(n)
362 0 : ns = sym%ntypsy(nat)
363 0 : jm = atoms%jri(n)
364 0 : DO nl = 0, sphhar%nlh(ns)
365 0 : lp = sphhar%llh(nl,ns)
366 0 : DO l = MERGE(1, lp - 1, lp.EQ.0), MERGE(1, lp + 1, lp.EQ.0), 2 ! Gaunt selection
367 0 : IF (l.GT.atoms%lmax(n)) CYCLE
368 0 : fint = atoms%rmt(n)**l * rho0(jm,nl,n)
369 0 : DO mb = 1, sphhar%nmem(nl,ns)
370 0 : mp = sphhar%mlh(mb,nl,ns)
371 0 : DO mVec = -1, 1
372 0 : m = mVec + mp ! Gaunt selection
373 0 : IF (ABS(m).GT.l) CYCLE
374 0 : gauntFactor = Gaunt1(l, 1, lp, m, mVec, mp, atoms%lmax(n))
375 : qlmo(m, l, n) = qlmo(m, l, n) + c_im(iDir, mVec + 2) * gauntFactor * &
376 0 : & sphhar%clnu(mb,nl,ns) * fint * pref
377 : END DO ! mVec
378 : END DO ! mb
379 : END DO ! l
380 : END DO ! nl
381 : END DO ! n
382 :
383 0 : END SUBROUTINE dfpt_mt_moments_SF
384 :
385 0 : SUBROUTINE dfpt_pw_moments_SF( fmpi, stars, atoms, cell, sym, iDtype, iDir, qpw_in, qlmp_SF )
386 :
387 : use m_mpi_bc_tool
388 : use m_mpi_reduce_tool
389 : use m_phasy1
390 : use m_sphbes
391 : use m_constants
392 : use m_types
393 : USE m_gaunt, only : Gaunt1
394 :
395 : implicit none
396 :
397 : type(t_mpi), intent(in) :: fmpi
398 : type(t_sym), intent(in) :: sym
399 : type(t_stars), intent(in) :: stars
400 : type(t_cell), intent(in) :: cell
401 : type(t_atoms), intent(in) :: atoms
402 : INTEGER, INTENT(IN) :: iDtype, iDir
403 : complex, intent(in) :: qpw_in(:)
404 : complex, intent(out) :: qlmp_SF(-atoms%lmaxd:,0:,:)
405 :
406 : integer :: n, k, l, ll1p, lmp, ierr, m, lp, mp, mVec, pref
407 : complex :: cil, nqpw
408 0 : complex :: pylm(( maxval( atoms%lmax ) + 1 ) ** 2, atoms%ntype)
409 : real :: sk3r, rl2
410 0 : real :: aj(0:maxval( atoms%lmax ) + 1 )
411 0 : complex, ALLOCATABLE :: qpw(:)
412 0 : complex :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
413 :
414 : ! TYPE(t_atoms), INTENT(IN) :: atoms
415 : ! TYPE(t_cell), INTENT(IN) :: cell
416 : ! INTEGER, INTENT(IN) :: ngdp
417 :
418 : ! INTEGER, INTENT(IN) :: gdp(:, :)
419 : ! COMPLEX, INTENT(IN) :: rho0IRpw(:)
420 : ! COMPLEX, INTENT(INOUT) :: qlmp(:, :,:)
421 :
422 : ! INTEGER :: iG, l, m, lp, mp, m2p, lm, lmp, iDir
423 : ! COMPLEX :: pref, phaseFac, temp1, temp2, temp3
424 : REAL :: gauntFactor
425 :
426 0 : ALLOCATE (qpw(stars%ng3))
427 0 : qpw = qpw_in(:stars%ng3)
428 0 : qlmp = 0.0
429 :
430 0 : pref = -1
431 0 : IF (iDtype.NE.0) pref = 1
432 :
433 0 : if ( fmpi%irank == 0 ) then
434 0 : do n = MERGE(1,iDtype,iDtype.EQ.0), MERGE(atoms%ntype,iDtype,iDtype.EQ.0)
435 0 : DO mVec = -1, 1
436 0 : qlmp(mVec,1,n) = pref * c_im(iDir, mVec + 2) * qpw(1) * stars%nstr(1) * atoms%rmt(n)**3
437 : END DO
438 : end do
439 : end if
440 :
441 0 : call mpi_bc(qpw,0,fmpi%mpi_comm)
442 :
443 0 : do k = fmpi%irank+2, stars%ng3, fmpi%isize
444 0 : call phasy1( atoms, stars, sym, cell, k, pylm )
445 :
446 0 : nqpw = qpw(k) * stars%nstr(k)
447 :
448 0 : do n = MERGE(1,iDtype,iDtype.EQ.0), MERGE(atoms%ntype,iDtype,iDtype.EQ.0)
449 :
450 0 : sk3r = stars%sk3(k) * atoms%rmt(n)
451 0 : call sphbes( atoms%lmax(n), sk3r, aj )
452 0 : rl2 = atoms%rmt(n) ** 2
453 :
454 0 : DO lp = 0, atoms%lmax(n)
455 0 : cil = aj(lp) * nqpw * rl2
456 0 : ll1p = lp * ( lp + 1 ) + 1
457 0 : DO l = MERGE(1, lp - 1, lp.EQ.0), MERGE(1, lp + 1, lp.EQ.0), 2 ! Gaunt selection
458 0 : IF (l.GT.atoms%lmax(n)) CYCLE
459 0 : DO mp = -lp, lp
460 0 : lmp = ll1p + mp
461 0 : DO mVec = -1, 1
462 0 : m = mVec + mp ! Gaunt selection
463 0 : IF (ABS(m).GT.l) CYCLE
464 0 : gauntFactor = Gaunt1(l, 1, lp, m, mVec, mp, atoms%lmax(n))
465 : qlmp(m,l,n) = qlmp(m,l,n) + c_im(iDir, mVec + 2) * gauntFactor * &
466 0 : & cil * atoms%rmt(n)**l * pylm(lmp,n) * pref
467 : END DO ! mVec
468 : END DO ! mp
469 : END DO ! l
470 : END DO ! lp
471 : END DO ! n = 1, atoms%ntype
472 : END DO ! k = 2, stars%ng3
473 :
474 0 : CALL mpi_sum_reduce(qlmp, qlmp_SF, fmpi%mpi_comm)
475 :
476 : ! pref = fpi_const * atoms%rmt(iDtype) * atoms%rmt(iDtype)
477 : ! DO iG = 1, ngdp
478 : ! gext = matmul(cell%bmat, real(gdp(:, iG)))
479 : ! gnorm = norm2(gExt)
480 :
481 : ! call ylm4( atoms%lmax(iDtype), gExt(1:3), ylm )
482 : ! call sphbes(atoms%lmax(iDtype), gnorm * atoms%rmt(iDtype), sbes)
483 :
484 : ! phaseFac = exp( ImagUnit * tpi_const * dot_product(gdp(:, iG), atoms%taual(:, iDatom)))
485 :
486 : ! DO l = 0, atoms%lmax(iDtype)
487 : ! temp1 = pref * phaseFac * atoms%rmt(iDtype)**l * rho0IRpw(iG)
488 : ! DO m = -l, l
489 : ! lm = l * (l + 1) + m + 1
490 : ! DO lp = 0, atoms%lmax(iDtype)
491 : ! temp2 = temp1 * sbes(lp) * ImagUnit**lp
492 : ! DO mp = -lp, lp
493 : ! lmp = lp * (lp + 1) + mp + 1
494 : ! temp3 = temp2 * conjg(ylm(lmp))
495 : ! DO m2p = -1, 1
496 : ! gauntFactor = Gaunt1( l, lp, 1, m, mp, m2p, atoms%lmax(iDtype))
497 : ! DO iDir = 1, 3
498 : ! qlmp(lm, iDatom, iDir) = qlmp(lm, iDatom, iDir) + c_im(iDir, m2p + 2) * gauntFactor * temp3
499 : ! END DO ! iDir
500 : ! END DO ! m2p
501 : ! END DO ! mp
502 : ! END DO ! lp
503 : ! END DO ! m
504 : ! END DO ! l
505 : ! END DO ! iG
506 :
507 0 : END SUBROUTINE dfpt_pw_moments_SF
508 :
509 : end module m_mpmom
|
