Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! Copyright (c) 2022 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !-------------------------------------------------------------------------------
6 :
7 : MODULE m_dfpt_rhonmt
8 : !! Module adapted from rhonmt21, which in turn was adapted from rhonmt, and
9 : !! from rhonmtlo for the lo part.
10 : !!
11 : !! It contains subroutines for the non-LO and LO non-spherical contributions
12 : !! to the density coefficients used to construct \(\rho_{L}^{\alpha}(r)\).
13 : USE m_gaunt,ONLY:gaunt1
14 : USE m_types_setup
15 : USE m_types_cdnval
16 : USE m_constants
17 :
18 : IMPLICIT NONE
19 :
20 : CONTAINS
21 :
22 8382 : SUBROUTINE dfpt_rhonmt(atoms,sphhar,we,we1,ne,ilSpinPr,ilSpin,qpoint,l_dfpt,l_less_effort,sym,eigVecCoeffs,eigVecCoeffs1,denCoeffs,eigVecCoeffs1m)
23 : !! Subroutine to construct all non-spherical MT density coefficients (for a
24 : !! density perturbation) without LOs in one routine. The spin input dictates,
25 : !! which element is gonna be built.
26 : !! The coefficients are of the form:
27 : !! $$\begin{aligned}
28 : !! d_{l',l,L,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha} &= \sum_{\nu\boldsymbol{k}}\sum_{m\mu(L)}\\
29 : !! &* c_{L,\mu}^{*}G_{l,l''(L),l'}^{m,m''(\mu),m-m''(\mu)}A_{l',m-m''(\mu),\lambda'}^{\sigma_{\alpha}',\nu\boldsymbol{k}*}\\
30 : !! &* (2\tilde{f}_{\nu\boldsymbol{k}}A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{q},j,\beta~(1)}+\tilde{f}_{\nu\boldsymbol{k}\boldsymbol{q}}^{j,\beta~(1)}A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}})
31 : !! \end{aligned}$$
32 : !! The k-point loop is performed outside this routine. In contrast to older
33 : !! routines, the arrays uunmt etc. and uunmt21 etc. are merged into one
34 : !! spin and u-order dependent array.
35 : !!
36 : !! \(G_{l,l'',l'}^{m,m'',m'}\): Gaunt coefficient
37 : !!
38 : !! \(\sigma_{\alpha}(')\): local spin indices \(\rightarrow\) ilSpinPr, ilSpin
39 : !!
40 : !! \(\lambda(')\): order of the radial basis function (0: u, 1: d)
41 : !!
42 : !! \(L\): Lattice harmonic index
43 : !!
44 : !! \(\mu(L)\): Lattice harmonic member index
45 : !!
46 : !! \(c_{L,\mu}\): Lattice harmonic coefficients
47 : !!
48 : !! \(\nu\boldsymbol{k}\): State index (k-point and number of state)
49 : !!
50 : !! \(\boldsymbol{q},j,\beta\): Perturbation index; q-point, direction and atom
51 : !!
52 : !! \(\tilde{f}_{\nu\boldsymbol{k}}\): (Smeared) occupation number [perturbed for \(\tilde{f}^{(1)}\)]
53 : !!
54 : !! \(A\): Summed matching coefficients and eigenvectors [perturbed for \(A^{(1)}\)]
55 :
56 : TYPE(t_sym), INTENT(IN) :: sym
57 : TYPE(t_sphhar), INTENT(IN) :: sphhar
58 : TYPE(t_atoms), INTENT(IN) :: atoms
59 :
60 : INTEGER, INTENT(IN) :: ne
61 : INTEGER, INTENT(IN) :: ilSpinPr !! \(\sigma_{\alpha}^{'}\)
62 : INTEGER, INTENT(IN) :: ilSpin !! \(\sigma_{\alpha}\)
63 : REAL, INTENT(IN) :: we(ne) !! \(\tilde{f}_{\nu\boldsymbol{k}}\)
64 : REAL, INTENT(IN) :: we1(ne) !! \(\tilde{f}_{\nu\boldsymbol{k}\boldsymbol{q}}^{j,\beta~(1)}\)
65 : REAL, INTENT(IN) :: qpoint(3) !! \(\boldsymbol{q}\)
66 : LOGICAL, INTENT(IN) :: l_dfpt, l_less_effort
67 :
68 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}}\)
69 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs1 !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{q},j,\beta~(1)}\)
70 :
71 : TYPE(t_denCoeffs), INTENT(INOUT) :: denCoeffs !! \(d_{l',l,L,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}\)
72 :
73 : TYPE(t_eigVecCoeffs), OPTIONAL, INTENT(IN) :: eigVecCoeffs1m !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{-q},j,\beta~(1)}\)
74 :
75 : COMPLEX :: coef, cil, cmv
76 8382 : COMPLEX :: temp(ne)
77 :
78 : INTEGER :: jmem,l,lh,llp,llpmax,lm,lmp,lp,lv,m, mp,mv,na,natom,nn,ns,nt,lphi,lplow
79 :
80 : LOGICAL :: l_minusq
81 :
82 8382 : l_minusq = PRESENT(eigVecCoeffs1m)
83 :
84 19044 : DO ns=1,sym%nsymt
85 : !$OMP parallel do default(none) &
86 : !$OMP private(lh,lp,l,lv,mp,m,mv,lm,lmp,llp,llpmax,lphi,lplow) &
87 : !$OMP private(cil,jmem,cmv,coef,temp,na,nt,nn,natom) &
88 19044 : !$OMP shared(sym,we,we1,ne,ns,atoms,sphhar,eigVecCoeffs,eigVecCoeffs1,denCoeffs,ilSpinPr,ilSpin,l_dfpt,l_less_effort,qpoint,l_minusq,eigVecCoeffs1m)
89 : ! !$OMP collapse(2)
90 : DO lh = 1, sphhar%nlh(ns)
91 : lv = sphhar%llh(lh,ns)
92 : DO jmem = 1, sphhar%nmem(lh,ns)
93 : mv = sphhar%mlh(jmem,lh,ns)
94 : cmv = conjg(sphhar%clnu(jmem,lh,ns))
95 : DO l = 0, atoms%lmaxd
96 : m_loop: DO m = -l,l
97 : lm= l*(l+1) + m
98 : mp = m - mv
99 :
100 : !maximum value of lp
101 : lphi = l + lv
102 : !---> check that lphi is smaller than the max l of the
103 : !---> wavefunction expansion
104 : lphi = MIN(lphi,atoms%lmaxd)
105 : !---> make sure that l + l'' + lphi is even
106 : lphi = lphi - MOD(l+lv+lphi,2)
107 :
108 : IF(l_less_effort) lphi = l - mod(lv,2)
109 :
110 : lplow = abs(l-lv)
111 : lplow = MAX(lplow,ABS(mp))
112 : !---> make sure that l + l'' + lplow is even
113 : lplow = lplow + MOD(ABS(lphi-lplow),2)
114 :
115 : IF (lplow.GT.lphi) CYCLE m_loop
116 :
117 : DO lp = lplow, lphi,2
118 : cil = ImagUnit**(l-lp)
119 : lmp = lp*(lp+1) + mp
120 : IF (lmp>lm.AND.l_less_effort) CYCLE m_loop
121 :
122 : coef= cmv * cil * gaunt1(l,lv,lp,m,mv,mp,atoms%lmaxd)
123 : !IF (ABS(coef) .LT. 1e-12 ) CYCLE
124 : natom= 0
125 : typeloop: DO nn=1,atoms%ntype
126 : IF (l_less_effort) THEN
127 : llp = (l* (l+1))/2 + lp
128 : llpmax = (atoms%lmax(nn)* (atoms%lmax(nn)+3))/2
129 : ELSE
130 : llp= lp*(atoms%lmax(nn)+1)+l
131 : llpmax = ((atoms%lmax(nn)+1)**2)-1
132 : END IF
133 :
134 : IF(llp.GT.llpmax) THEN
135 : natom = natom + atoms%neq(nn)
136 : CYCLE typeloop
137 : END IF
138 :
139 : nt= natom
140 : DO na= 1,atoms%neq(nn)
141 : nt= nt+1
142 : IF (sym%ntypsy(nt)==ns) THEN
143 : ! uu/du
144 : temp(:) = coef * we(:) * eigVecCoeffs1%abcof(:,lm,0,nt,ilSpin) ! If not DFPT, this is the base case for rhonmt(21)
145 : IF (lmp/=lm.AND.l_less_effort) temp(:) = temp(:) * 2.0
146 : IF (l_dfpt) THEN
147 : IF (.NOT.l_minusq) temp(:) = temp(:) * 2.0
148 : IF (norm2(qpoint)<=1e-8) THEN
149 : temp(:) = temp(:) + coef * we1(:) * eigVecCoeffs%abcof(:,lm,0,nt,ilSpin)
150 : END IF
151 : END IF
152 : denCoeffs%nmt_coeff(llp,lh,nn,0,0,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,0,0,ilSpinPr,ilSpin) &
153 : & + dot_product(eigVecCoeffs%abcof(:ne,lmp,0,nt,ilSpinPr),temp(:ne))
154 : denCoeffs%nmt_coeff(llp,lh,nn,1,0,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,1,0,ilSpinPr,ilSpin) &
155 : & + dot_product(eigVecCoeffs%abcof(:ne,lmp,1,nt,ilSpinPr),temp(:ne))
156 :
157 : IF (l_minusq) THEN
158 : denCoeffs%nmt_coeff(llp,lh,nn,0,0,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,0,0,ilSpinPr,ilSpin) &
159 : & + dot_product(eigVecCoeffs1m%abcof(:ne,lmp,0,nt,ilSpinPr), &
160 : & eigVecCoeffs%abcof(:,lm,0,nt,ilSpin))
161 : denCoeffs%nmt_coeff(llp,lh,nn,1,0,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,1,0,ilSpinPr,ilSpin) &
162 : & + dot_product(eigVecCoeffs1m%abcof(:ne,lmp,1,nt,ilSpinPr), &
163 : & eigVecCoeffs%abcof(:,lm,0,nt,ilSpin))
164 : END IF
165 :
166 : ! dd/ud
167 : temp(:) = coef * we(:) * eigVecCoeffs1%abcof(:,lm,1,nt,ilSpin)
168 : IF (lmp/=lm.AND.l_less_effort) temp(:) = temp(:) * 2.0
169 : IF (l_dfpt) THEN
170 : IF (.NOT.l_minusq) temp(:) = temp(:) * 2.0
171 : IF (norm2(qpoint)<=1e-8) THEN
172 : temp(:) = temp(:) + coef * we1(:) * eigVecCoeffs%abcof(:,lm,1,nt,ilSpin)
173 : END IF
174 : END IF
175 : denCoeffs%nmt_coeff(llp,lh,nn,1,1,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,1,1,ilSpinPr,ilSpin) &
176 : & + dot_product(eigVecCoeffs%abcof(:ne,lmp,1,nt,ilSpinPr),temp(:ne))
177 : denCoeffs%nmt_coeff(llp,lh,nn,0,1,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,0,1,ilSpinPr,ilSpin) &
178 : & + dot_product(eigVecCoeffs%abcof(:ne,lmp,0,nt,ilSpinPr),temp(:ne))
179 : IF (l_minusq) THEN
180 : denCoeffs%nmt_coeff(llp,lh,nn,1,1,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,1,1,ilSpinPr,ilSpin) &
181 : & + dot_product(eigVecCoeffs1m%abcof(:ne,lmp,1,nt,ilSpinPr), &
182 : & eigVecCoeffs%abcof(:,lm,1,nt,ilSpin))
183 : denCoeffs%nmt_coeff(llp,lh,nn,0,1,ilSpinPr,ilSpin) = denCoeffs%nmt_coeff(llp,lh,nn,0,1,ilSpinPr,ilSpin) &
184 : & + dot_product(eigVecCoeffs1m%abcof(:ne,lmp,0,nt,ilSpinPr), &
185 : & eigVecCoeffs%abcof(:,lm,1,nt,ilSpin))
186 : END IF
187 : ENDIF ! (sym%ntypsy(nt)==ns)
188 : ENDDO ! na
189 : natom = natom + atoms%neq(nn)
190 : END DO typeloop! nn
191 : END DO
192 : END DO m_loop ! m
193 : END DO ! jmem
194 : END DO ! l
195 : END DO ! lh
196 : !$OMP end parallel do
197 : END DO ! ns
198 8382 : END SUBROUTINE dfpt_rhonmt
199 :
200 8382 : SUBROUTINE dfpt_rhonmtlo(atoms,sphhar,sym,ne,we,we1,eigVecCoeffs,eigVecCoeffs1,denCoeffs,ilSpinPr,ilSpin,l_dfpt,qpoint,eigVecCoeffs1m)
201 : !! This is a complementary routine to the one above for \(\lambda(')\ge 2\),
202 : !! i.e. mixed or pure LO contributions.
203 : USE m_gaunt,ONLY:gaunt1
204 : USE m_types
205 : use m_constants
206 :
207 : IMPLICIT NONE
208 :
209 : TYPE(t_sphhar), INTENT(IN) :: sphhar
210 : TYPE(t_atoms), INTENT(IN) :: atoms
211 : TYPE(t_sym), INTENT(IN) :: sym
212 :
213 : INTEGER, INTENT (IN) :: ne, ilSpinPr, ilSpin
214 :
215 : REAL, INTENT (IN) :: we(:),we1(:)!(nobd)
216 :
217 : REAL, INTENT(IN) :: qpoint(3)
218 : LOGICAL, INTENT(IN) :: l_dfpt
219 :
220 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
221 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs1
222 :
223 : TYPE(t_denCoeffs), INTENT(INOUT) :: denCoeffs
224 :
225 : TYPE(t_eigVecCoeffs), OPTIONAL, INTENT(IN) :: eigVecCoeffs1m
226 :
227 : COMPLEX :: cmv,fact,cf1, cf2
228 : INTEGER :: i,jmem,l,lh,lmp,lo,lop,lp,lpmax,lpmax0,lpmin,lpmin0,m,lpp ,mp,mpp,na,neqat0,nn,ntyp
229 :
230 : LOGICAL :: l_minusq
231 :
232 8382 : l_minusq = PRESENT(eigVecCoeffs1m)
233 :
234 8382 : neqat0 = 0
235 26088 : DO ntyp = 1,atoms%ntype
236 450220 : DO lh = 1,sphhar%nlh(sym%ntypsy(neqat0+1))
237 432514 : lpp = sphhar%llh(lh,sym%ntypsy(neqat0+1))
238 1255476 : DO jmem = 1,sphhar%nmem(lh,sym%ntypsy(neqat0+1))
239 805256 : mpp = sphhar%mlh(jmem,lh,sym%ntypsy(neqat0+1))
240 805256 : cmv = CONJG(sphhar%clnu(jmem,lh,sym%ntypsy(neqat0+1)))
241 2568556 : DO lo = 1,atoms%nlo(ntyp)
242 1330786 : l = atoms%llo(lo,ntyp)
243 1330786 : lpmin0 = ABS(l-lpp)
244 1330786 : lpmax0 = l + lpp
245 :
246 1330786 : lpmax = MIN(lpmax0,atoms%lmax(ntyp))
247 1330786 : lpmax = lpmax - MOD(l+lpp+lpmax,2)
248 4912140 : DO m = -l,l
249 2776098 : mp = m - mpp
250 2776098 : lpmin = MAX(lpmin0,ABS(mp))
251 2776098 : lpmin = lpmin + MOD(l+lpp+lpmin,2)
252 6765508 : DO lp = lpmin,lpmax,2
253 3989410 : lmp = lp* (lp+1) + mp
254 3989410 : fact = cmv* (ImagUnit** (l-lp))*gaunt1(l,lp,lpp,m,mp,mpp,atoms%lmaxd)
255 3989410 : na = neqat0
256 10793628 : DO nn = 1,atoms%neq(ntyp)
257 4028120 : na = na + 1
258 72360090 : DO i = 1,ne
259 64342560 : cf1 = fact * we(i) * eigVecCoeffs1%ccof(m,i,lo,na,ilSpin)! If not DFPT, this is the base case for rhonmtlo
260 64342560 : IF (l_dfpt) THEN
261 0 : cf1 = cf1 * 2.0
262 0 : IF (norm2(qpoint)<=1e-8) THEN
263 0 : cf1 = cf1 + fact * we1(i) * eigVecCoeffs%ccof(m,i,lo,na,ilSpin)
264 : END IF
265 : END IF
266 : denCoeffs%nmt_ulo_coeff(lp,lo,lh,ntyp,0,ilSpinPr,ilSpin) = denCoeffs%nmt_ulo_coeff(lp,lo,lh,ntyp,0,ilSpinPr,ilSpin) &
267 64342560 : & + CONJG(eigVecCoeffs%abcof(i,lmp,0,na,ilSpinPr)) * cf1
268 : denCoeffs%nmt_ulo_coeff(lp,lo,lh,ntyp,1,ilSpinPr,ilSpin) = denCoeffs%nmt_ulo_coeff(lp,lo,lh,ntyp,1,ilSpinPr,ilSpin) &
269 68370680 : & + CONJG(eigVecCoeffs%abcof(i,lmp,1,na,ilSpinPr)) * cf1
270 : END DO
271 : END DO
272 : END DO
273 :
274 2776098 : mp = m + mpp
275 2776098 : lpmin = MAX(lpmin0,ABS(mp))
276 2776098 : lpmin = lpmin + MOD(l+lpp+lpmin,2)
277 6765508 : DO lp = lpmin,lpmax,2
278 3989410 : lmp = lp* (lp+1) + mp
279 3989410 : fact = cmv* (ImagUnit** (lp-l))*gaunt1(lp,l,lpp,mp,m,mpp,atoms%lmaxd)
280 3989410 : na = neqat0
281 10793628 : DO nn = 1,atoms%neq(ntyp)
282 4028120 : na = na + 1
283 72360090 : DO i = 1,ne
284 64342560 : cf1 = fact * we(i) * eigVecCoeffs1%abcof(i,lmp,0,na,ilSpin)! If not DFPT, this is the base case for rhonmtlo
285 64342560 : cf2 = fact * we(i) * eigVecCoeffs1%abcof(i,lmp,1,na,ilSpin)
286 64342560 : IF (l_dfpt) THEN
287 0 : cf1 = cf1 * 2.0
288 0 : cf2 = cf2 * 2.0
289 0 : IF (norm2(qpoint)<=1e-8) THEN
290 0 : cf1 = cf1 + fact * we1(i) * eigVecCoeffs%abcof(i,lmp,0,na,ilSpin)
291 0 : cf2 = cf2 + fact * we1(i) * eigVecCoeffs%abcof(i,lmp,1,na,ilSpin)
292 : END IF
293 : END IF
294 : denCoeffs%nmt_lou_coeff(lp,lo,lh,ntyp,0,ilSpinPr,ilSpin) = denCoeffs%nmt_lou_coeff(lp,lo,lh,ntyp,0,ilSpinPr,ilSpin) &
295 64342560 : & + CONJG(eigVecCoeffs%ccof(m,i,lo,na,ilSpinPr)) * cf1
296 : denCoeffs%nmt_lou_coeff(lp,lo,lh,ntyp,1,ilSpinPr,ilSpin) = denCoeffs%nmt_lou_coeff(lp,lo,lh,ntyp,1,ilSpinPr,ilSpin) &
297 68370680 : & + CONJG(eigVecCoeffs%ccof(m,i,lo,na,ilSpinPr)) * cf2
298 : END DO
299 : END DO
300 : END DO
301 :
302 9359062 : DO lop = 1,atoms%nlo(ntyp)
303 5252178 : lp = atoms%llo(lop,ntyp)
304 5252178 : mp = m - mpp
305 8028276 : IF ((ABS(l-lpp).LE.lp) .AND.(lp.LE. (l+lpp)) .AND.(MOD(l+lp+lpp,2).EQ.0) .AND.(ABS(mp).LE.lp)) THEN
306 117196 : fact = cmv* (ImagUnit** (l-lp))*gaunt1(l,lp,lpp,m,mp,mpp,atoms%lmaxd)
307 117196 : na = neqat0
308 235446 : DO nn = 1,atoms%neq(ntyp)
309 118250 : na = na + 1
310 2096747 : DO i = 1,ne
311 1861301 : cf1 = fact * we(i) * eigVecCoeffs1%ccof(m,i,lo,na,ilSpin)! If not DFPT, this is the base case for rhonmtlo
312 1861301 : IF (l_dfpt) THEN
313 0 : cf1 = cf1 * 2.0
314 0 : IF (norm2(qpoint)<=1e-8) THEN
315 0 : cf1 = cf1 + fact * we1(i) * eigVecCoeffs%ccof(m,i,lo,na,ilSpin)
316 : END IF
317 : END IF
318 : denCoeffs%nmt_lolo_coeff(lop,lo,lh,ntyp,ilSpinPr,ilSpin) = denCoeffs%nmt_lolo_coeff(lop,lo,lh,ntyp,ilSpinPr,ilSpin) &
319 1979551 : & + CONJG(eigVecCoeffs%ccof(mp,i,lop,na,ilSpinPr)) * cf1
320 : END DO
321 : END DO
322 : END IF
323 : END DO
324 : END DO
325 : END DO
326 : END DO
327 : END DO
328 26088 : neqat0 = neqat0 + atoms%neq(ntyp)
329 : END DO
330 8382 : END SUBROUTINE dfpt_rhonmtlo
331 : END MODULE m_dfpt_rhonmt
|
