Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 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_tlo
8 : USE m_juDFT
9 : !***********************************************************************
10 : ! sets up the extra t-matrix elements due to the local orbitals.
11 : ! only non=zero elements are calculated
12 : !
13 : ! p.kurz jul. 1996
14 : !***********************************************************************
15 : CONTAINS
16 2332 : SUBROUTINE tlo(atoms,sym,sphhar,iSpinPr,iSpin,jsp,ntyp,enpara,lh0,input,vr,&
17 1166 : na,flo,f,g,usdus, tlmplm, one, l_dfpt, l_V1)
18 : ! Contruct the additional local Hamiltonian matrices
19 : ! t_{L'L}^{\mu} = \sum_{lh} \int dV u_{l',order'}^{\mu}(r)Y_{l'}^{m'*}(\Omega)
20 : ! * V_{lh}(r)Y_{lh}(\Omega)
21 : ! * u_{l,order}^{\mu}(r)Y_{l}^{m}(\Omega)
22 : ! * i^{l-l'}
23 : ! + \int dV u_{l',order'}^{\mu}(r) H_{sph}
24 : ! * u_{l,order}^{\mu}(r)Y_{l}^{m}(\Omega)
25 : ! of a real valued potential V(\bm{r}). The superindex L is defined as
26 : ! L := (l,m,order)
27 : ! with order = 0 refering to radial functions u and order = 1 denoting
28 : ! their energy derivatives (udot). This construction is not k-dependent
29 : ! and therefore executed only once each scf iteration.
30 :
31 : ! Abbreviations:
32 : ! tuulo: t-matrix element of an LO and the APW radial fuction
33 : ! tdulo: t-matrix element of an LO and the energy derivative of
34 : ! the APW radial fuction
35 : ! tulou: t-matrix element of the APW radial fuction and an LO
36 : ! tulod: t-matrix element of the APW radial fuction derivative and an LO
37 : ! tuloulo: t-matrix element of two LOs
38 :
39 : USE m_intgr, ONLY : intgr3
40 : USE m_gaunt, ONLY: gaunt1
41 : USE m_types
42 : USE m_constants
43 :
44 : IMPLICIT NONE
45 :
46 : TYPE(t_input), INTENT(IN) :: input
47 : TYPE(t_sphhar), INTENT(IN) :: sphhar
48 : TYPE(t_atoms), INTENT(IN) :: atoms
49 : TYPE(t_sym), INTENT(IN) :: sym
50 : TYPE(t_usdus), INTENT(IN) :: usdus
51 : TYPE(t_tlmplm), INTENT(INOUT) :: tlmplm
52 : TYPE(t_enpara), INTENT(IN) :: enpara
53 :
54 : ! Scalar Arguments
55 : INTEGER, INTENT (IN) :: iSpinPr,iSpin,jsp,ntyp ,lh0,na
56 :
57 : ! Array Arguments
58 : REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd)
59 : REAL, INTENT (IN) :: f(:,:,0:,:),g(:,:,0:,:) !(atoms%jmtd,2,0:atoms%lmaxd,spins)
60 : REAL, INTENT (IN) :: flo(:,:,:,:)!(atoms%jmtd,2,atoms%nlod,spins)
61 :
62 : COMPLEX, INTENT(IN) :: one
63 :
64 : LOGICAL, INTENT(IN) :: l_dfpt, l_V1
65 :
66 : ! Local Scalars
67 : COMPLEX :: cil
68 : INTEGER :: i,l,lh,lm ,lmin,lmp,lo,lop,loplo,lp,lpmax,lpmax0,lpmin,lpmin0,lpp ,mem,mp,mpp,m,lmx,mlo,mlolo,s
69 : INTEGER :: loplo_new, mlolo_new
70 :
71 : ! Local Arrays
72 1166 : REAL :: x(atoms%jmtd),ulovulo(atoms%nlod*(atoms%nlod+1)/2,lh0:sphhar%nlhd)
73 1166 : REAL :: uvulo(atoms%nlod,0:atoms%lmaxd,lh0:sphhar%nlhd),dvulo(atoms%nlod,0:atoms%lmaxd,lh0:sphhar%nlhd)
74 :
75 3112 : DO lo = 1,atoms%nlo(ntyp)
76 1946 : l = atoms%llo(lo,ntyp)
77 21446 : DO lp = 0,atoms%lmax(ntyp)
78 18334 : lmin = ABS(lp-l)
79 18334 : lmx = lp + l
80 717354 : DO lh = lh0,sphhar%nlh(sym%ntypsy(na))
81 697074 : lpp = sphhar%llh(lh,sym%ntypsy(na))
82 715408 : IF ((MOD(l+lp+lpp,2).EQ.1).OR.(lpp.LT.lmin).OR.(lpp.GT.lmx)) THEN
83 590444 : uvulo(lo,lp,lh) = 0.0
84 590444 : dvulo(lo,lp,lh) = 0.0
85 : ELSE
86 82764164 : DO i = 1,atoms%jri(ntyp)
87 82764164 : x(i) = (f(i,1,lp,iSpinPr)*flo(i,1,lo,iSpin)+ f(i,2,lp,iSpinPr)*flo(i,2,lo,iSpin))*vr(i,lh)
88 : END DO
89 106630 : CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),uvulo(lo,lp,lh))
90 :
91 82764164 : DO i = 1,atoms%jri(ntyp)
92 82764164 : x(i) = (g(i,1,lp,iSpinPr)*flo(i,1,lo,iSpin)+ g(i,2,lp,iSpinPr)*flo(i,2,lo,iSpin))*vr(i,lh)
93 : END DO
94 106630 : CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),dvulo(lo,lp,lh))
95 : END IF
96 : END DO
97 : END DO
98 : END DO
99 1166 : loplo = 0
100 3112 : DO lop = 1,atoms%nlo(ntyp)
101 1946 : lp = atoms%llo(lop,ntyp)
102 5854 : DO lo = 1,lop
103 2742 : l = atoms%llo(lo,ntyp)
104 2742 : loplo = loplo + 1
105 2742 : IF (loplo>SIZE(ulovulo,1)) CALL juDFT_error("loplo too large!!!" ,calledby ="tlo")
106 111674 : DO lh = lh0,sphhar%nlh(sym%ntypsy(na))
107 106986 : lpp = sphhar%llh(lh,sym%ntypsy(na))
108 106986 : lmin = ABS(lp - l)
109 106986 : lmx = lp + l
110 109728 : IF ((MOD(l+lp+lpp,2).EQ.1).OR.(lpp.LT.lmin).OR.(lpp.GT.lmx)) THEN
111 102596 : ulovulo(loplo,lh) = 0.0
112 : ELSE
113 3372236 : DO i = 1,atoms%jri(ntyp)
114 3372236 : x(i) = (flo(i,1,lop,iSpinPr)*flo(i,1,lo,iSpin)+flo(i,2,lop,iSpinPr)*flo(i,2,lo,iSpin))*vr(i,lh)
115 : END DO
116 4390 : CALL intgr3(x,atoms%rmsh(:,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ulovulo(loplo,lh))
117 : END IF
118 : END DO
119 : END DO
120 : END DO
121 :
122 : ! Generate the t-matrices. For optimal performance consider only those
123 : ! combinations of l,l',l'',m,m',m'' that satisfy the three conditions for
124 : ! non-zero Gaunt coefficients, i.e.
125 : ! |l - l''| <= l' <= l + l'' (triangular condition)
126 : ! m' = m + m''
127 : ! l + l' + l'' even
128 :
129 : ! Loop over the local orbitals
130 1688 : mlo=SUM(atoms%nlo(:ntyp-1))
131 1166 : s=tlmplm%h_loc2_nonsph(ntyp)
132 3112 : DO lo = 1,atoms%nlo(ntyp)
133 1946 : l = atoms%llo(lo,ntyp)
134 7742 : DO m = -l,l
135 4630 : lm=l*(l+1)+m
136 : ! Loop over the lattice harmonics
137 165602 : DO lh = lh0,sphhar%nlh(sym%ntypsy(na))
138 159026 : lpp = sphhar%llh(lh,sym%ntypsy(na))
139 159026 : lpmin0 = ABS(l-lpp)
140 159026 : lpmax0 = l + lpp
141 : ! Check that lpmax is smaller than the max l of the
142 : ! wavefunction expansion at this atom
143 159026 : lpmax = MIN(lpmax0,atoms%lnonsph(ntyp))
144 : ! Make sure that l + l'' + lpmax is even
145 159026 : lpmax = lpmax - MOD(l+lpp+lpmax,2)
146 470914 : DO mem = 1,sphhar%nmem(lh,sym%ntypsy(na))
147 307258 : mpp = sphhar%mlh(mem,lh,sym%ntypsy(na))
148 307258 : mp = m + mpp
149 307258 : lpmin = MAX(lpmin0,ABS(mp))
150 : !- Make sure that l + l'' + lpmin is even
151 307258 : lpmin = lpmin + MOD(ABS(lpmax-lpmin),2)
152 : ! Loop over l'
153 466284 : DO lp = lpmin,lpmax,2
154 292990 : lmp = lp* (lp+1) + mp
155 : cil = ImagUnit**(l-lp) * sphhar%clnu(mem,lh,sym%ntypsy(na)) &
156 292990 : & * gaunt1(lp,lpp,l,mp,mpp,m,atoms%lmaxd)
157 : tlmplm%tuulo(lmp,m,lo+mlo,iSpinPr,iSpin) = &
158 292990 : & tlmplm%tuulo(lmp,m,lo+mlo,iSpinPr,iSpin) + one * cil * uvulo(lo,lp,lh)
159 : tlmplm%tdulo(lmp,m,lo+mlo,iSpinPr,iSpin) = &
160 292990 : & tlmplm%tdulo(lmp,m,lo+mlo,iSpinPr,iSpin) + one * cil * dvulo(lo,lp,lh)
161 292990 : tlmplm%h_LO(lmp,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lmp,m,lo+mlo,iSpinPr,iSpin) + one * cil * uvulo(lo,lp,lh)
162 292990 : tlmplm%h_LO(lmp+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lmp+s,m,lo+mlo,iSpinPr,iSpin) + one * cil * dvulo(lo,lp,lh)
163 : tlmplm%tulou(lmp,m,lo+mlo,iSpinPr,iSpin) = &
164 292990 : & tlmplm%tulou(lmp,m,lo+mlo,iSpinPr,iSpin) + one * CONJG(cil*uvulo(lo,lp,lh))
165 : tlmplm%tulod(lmp,m,lo+mlo,iSpinPr,iSpin) = &
166 292990 : & tlmplm%tulod(lmp,m,lo+mlo,iSpinPr,iSpin) + one * CONJG(cil*dvulo(lo,lp,lh))
167 292990 : tlmplm%h_LO2(lmp,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lmp,m,lo+mlo,iSpinPr,iSpin) + one * CONJG(cil*uvulo(lo,lp,lh))
168 389408 : tlmplm%h_LO2(lmp+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lmp+s,m,lo+mlo,iSpinPr,iSpin) + one * CONJG(cil*dvulo(lo,lp,lh))
169 : END DO
170 : END DO
171 : END DO
172 : END DO
173 : END DO
174 :
175 : ! Generate the t-matrix including two local orbitals for LO <= LO'
176 : ! Loop over LO'
177 1688 : mlolo = DOT_PRODUCT(atoms%nlo(:ntyp-1),atoms%nlo(:ntyp-1)+1)/2
178 : mlolo_new = DOT_PRODUCT(atoms%nlo(:ntyp-1),atoms%nlo(:ntyp-1))
179 3112 : DO lop = 1,atoms%nlo(ntyp)
180 1946 : lp = atoms%llo(lop,ntyp)
181 7742 : DO mp = -lp,lp
182 : ! Loop over the lattice harmonics
183 165602 : DO lh = lh0,sphhar%nlh(sym%ntypsy(na))
184 159026 : lpp = sphhar%llh(lh,sym%ntypsy(na))
185 470914 : DO mem = 1,sphhar%nmem(lh,sym%ntypsy(na))
186 307258 : mpp = sphhar%mlh(mem,lh,sym%ntypsy(na))
187 307258 : m = mp - mpp
188 : ! Loop over LO
189 962342 : DO lo = 1,lop
190 496058 : l = atoms%llo(lo,ntyp)
191 496058 : loplo = ((lop-1)*lop)/2 + lo
192 496058 : loplo_new = (lop-1) * atoms%nlo(ntyp) + lo
193 803316 : IF ((ABS(l-lpp).LE.lp).AND.(lp.LE.(l+lpp)).AND.(MOD(l+lp+lpp,2).EQ.0).AND.(ABS(m).LE.l)) THEN
194 : cil = ImagUnit**(l-lp) * sphhar%clnu(mem,lh,sym%ntypsy(na)) &
195 11970 : & * gaunt1(lp,lpp,l,mp,mpp,m,atoms%lmaxd)
196 : tlmplm%tuloulo(mp,m,loplo+mlolo,iSpinPr,iSpin) = &
197 11970 : & tlmplm%tuloulo(mp,m,loplo+mlolo,iSpinPr,iSpin) + one * cil * ulovulo(loplo,lh)
198 : ! tlmplm%tuloulo_new(mp,m,mlolo_new+loplo_new,iSpinPr,iSpin) = &
199 : ! & tlmplm%tuloulo_new(mp,m,mlolo_new+loplo_new,iSpinPr,iSpin) + one * cil * ulovulo(loplo,lh)
200 : tlmplm%tuloulo_newer(mp,m,lop,lo,ntyp,iSpinPr,iSpin) = &
201 11970 : & tlmplm%tuloulo_newer(mp,m,lop,lo,ntyp,iSpinPr,iSpin) + one * cil * ulovulo(loplo,lh)
202 11970 : IF (lop.NE.lo) THEN
203 : !loplo = ((lo-1)*lo)/2 + lop
204 : !loplo_new = (lo-1) * atoms%nlo(ntyp) + lop
205 : ! tlmplm%tuloulo_new(m,mp,mlolo_new+loplo_new,iSpinPr,iSpin) = &
206 : ! & tlmplm%tuloulo_new(m,mp,mlolo_new+loplo_new,iSpinPr,iSpin) + one * CONJG(cil * ulovulo(loplo,lh))
207 : tlmplm%tuloulo_newer(m,mp,lo,lop,ntyp,iSpinPr,iSpin) = &
208 1688 : & tlmplm%tuloulo_newer(m,mp,lo,lop,ntyp,iSpinPr,iSpin) + one * CONJG(cil * ulovulo(loplo,lh))
209 : END IF
210 : END IF
211 : END DO
212 : END DO
213 : END DO
214 : END DO
215 : END DO
216 :
217 : ! TODO: Do *not* symmetrize the local hamiltonian in dfpt.
218 : ! Add the diagonal terms from the spherical Hamiltonian. These terms have
219 : ! to be made Hermitian. If second variation is switched on, the t-matrices
220 : ! contain only the contributions from the non-spherical Hamiltonian.
221 1166 : IF (.NOT.input%secvar.AND.iSpinPr==iSpin.AND..NOT.l_V1) THEN
222 2752 : DO lo = 1,atoms%nlo(ntyp)
223 1706 : l = atoms%llo(lo,ntyp)
224 6902 : DO m = -l,l
225 4150 : lm = l* (l+1) + m
226 : !IF (.NOT.l_dfpt) THEN
227 1706 : IF (.TRUE.) THEN
228 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) &
229 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr) &
230 4150 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) )
231 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin) &
232 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr) &
233 4150 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) )
234 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) &
235 : + 0.5 * usdus%dulon(lo,ntyp,iSpinPr) &
236 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) ) &
237 4150 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr)
238 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin) &
239 : + 0.5 * usdus%dulon(lo,ntyp,iSpinPr) &
240 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) ) &
241 4150 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr)
242 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) &
243 : & + 0.5 * usdus%uulon(lo,ntyp,iSpinPr) &
244 4150 : & * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) )
245 : tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin) &
246 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr) &
247 4150 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) )
248 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) &
249 : & + 0.5 * usdus%dulon(lo,ntyp,iSpinPr) &
250 : & * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) ) &
251 4150 : & + 0.5 * usdus%uulon(lo,ntyp,iSpinPr)
252 : tlmplm%h_LO2(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lm+s,m,lo+mlo,iSpinPr,iSpin) &
253 : + 0.5 * usdus%dulon(lo,ntyp,iSpinPr) &
254 : * ( enpara%el0(l,ntyp,iSpinPr)+enpara%ello0(lo,ntyp,iSpinPr) ) &
255 4150 : + 0.5 * usdus%uulon(lo,ntyp,iSpinPr)
256 4150 : IF (atoms%ulo_der(lo,ntyp).GE.1) THEN
257 20 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr)
258 20 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr)
259 20 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr)
260 20 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr)
261 :
262 20 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr)
263 20 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr)
264 20 : tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr)
265 20 : tlmplm%h_LO2(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lm+s,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr)
266 : END IF
267 : !+apw_lo
268 4150 : IF (atoms%l_dulo(lo,ntyp)) THEN
269 0 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
270 0 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
271 0 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)= 0.0
272 0 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = 0.0
273 0 : tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO2(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
274 0 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
275 0 : tlmplm%h_LO2(lm+s,m,lo+mlo,iSpinPr,iSpin)= 0.0
276 0 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = 0.0
277 : END IF
278 : !+apw_lo
279 : ELSE
280 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) &
281 : + usdus%uulon(lo,ntyp,iSpinPr) &
282 : * enpara%ello0(lo,ntyp,iSpinPr)
283 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)&
284 : + usdus%uulon(lo,ntyp,iSpinPr) &
285 : * enpara%ello0(lo,ntyp,iSpinPr)
286 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) &
287 : + usdus%dulon(lo,ntyp,iSpinPr) &
288 : * enpara%ello0(lo,ntyp,iSpinPr) &
289 : + 0.0 * usdus%uulon(lo,ntyp,iSpinPr)
290 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)&
291 : + usdus%dulon(lo,ntyp,iSpinPr) &
292 : * enpara%ello0(lo,ntyp,iSpinPr) &
293 : + 0.0 * usdus%uulon(lo,ntyp,iSpinPr)
294 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) &
295 : & + usdus%uulon(lo,ntyp,iSpinPr) &
296 : & * enpara%el0(l,ntyp,iSpinPr)
297 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) &
298 : & + usdus%dulon(lo,ntyp,iSpinPr) &
299 : & * enpara%el0(l,ntyp,iSpinPr) &
300 : & + 1.0 * usdus%uulon(lo,ntyp,iSpinPr)
301 : ! TODO: Implement boundary term.
302 : IF (atoms%ulo_der(lo,ntyp).GE.1) THEN
303 : CALL juDFT_error("ulo_der>0 for DFPT" ,calledby ="tlo")
304 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr) !TODO: 1.0 or 0.0?
305 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr) !TODO: 1.0 or 0.0?
306 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr)
307 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr)
308 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%uuilon(lo,ntyp,iSpinPr) !TODO: 1.0 or 0.0?
309 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5 * usdus%duilon(lo,ntyp,iSpinPr) !TODO: 1.0 or 0.0?
310 : END IF
311 : IF (atoms%l_dulo(lo,ntyp)) THEN
312 : CALL juDFT_error("l_dulo for DFPT" ,calledby ="tlo")
313 : tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tuulo(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
314 : tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)=tlmplm%h_LO(lm,m,lo+mlo,iSpinPr,iSpin)+0.5
315 : tlmplm%h_LO(lm+s,m,lo+mlo,iSpinPr,iSpin)=0.0
316 : tlmplm%tdulo(lm,m,lo+mlo,iSpinPr,iSpin) = 0.0
317 : tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) = tlmplm%tulou(lm,m,lo+mlo,iSpinPr,iSpin) + 0.5
318 : tlmplm%tulod(lm,m,lo+mlo,iSpinPr,iSpin) = 0.0
319 : END IF
320 : END IF
321 : END DO
322 : END DO
323 2752 : DO lop = 1,atoms%nlo(ntyp)
324 1706 : lp = atoms%llo(lop,ntyp)
325 4462 : DO lo = atoms%lo1l(lp,ntyp),lop
326 1710 : loplo = ((lop-1)*lop)/2 + lo
327 1710 : loplo_new = (lop-1) * atoms%nlo(ntyp) + lo
328 7578 : DO m = -lp,lp
329 : !IF (.NOT.l_dfpt) THEN
330 1710 : IF (.TRUE.) THEN
331 : tlmplm%tuloulo(m,m,loplo+mlolo,iSpinPr,iSpin) = tlmplm%tuloulo(m,m,loplo+mlolo,iSpinPr,iSpin) &
332 : & + 0.5 * (enpara%ello0(lop,ntyp,iSpinPr) &
333 : & + enpara%ello0(lo,ntyp,iSpinPr)) &
334 : & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
335 : & + 0.5 * (usdus%ulouilopn(lop,lo,ntyp,iSpinPr) &
336 4162 : & + usdus%ulouilopn(lo,lop,ntyp,iSpinPr))
337 : !tlmplm%tuloulo_new(m,m,mlolo_new+loplo_new,iSpinPr,iSpin) = tlmplm%tuloulo_new(m,m,mlolo_new+loplo_new,iSpinPr,iSpin) &
338 : ! & + 0.5 * (enpara%ello0(lop,ntyp,iSpinPr) &
339 : ! & + enpara%ello0(lo,ntyp,iSpinPr)) &
340 : ! & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
341 : ! & + 0.5 * (usdus%ulouilopn(lop,lo,ntyp,iSpinPr) &
342 : ! & + usdus%ulouilopn(lo,lop,ntyp,iSpinPr))
343 : tlmplm%tuloulo_newer(m,m,lop,lo,ntyp,iSpinPr,iSpin) = tlmplm%tuloulo_newer(m,m,lop,lo,ntyp,iSpinPr,iSpin) &
344 : & + 0.5 * (enpara%ello0(lop,ntyp,iSpinPr) &
345 : & + enpara%ello0(lo,ntyp,iSpinPr)) &
346 : & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
347 : & + 0.5 * (usdus%ulouilopn(lop,lo,ntyp,iSpinPr) &
348 4162 : & + usdus%ulouilopn(lo,lop,ntyp,iSpinPr))
349 4162 : IF (.NOT.lop.EQ.lo) THEN
350 : !loplo_new = (lo-1) * atoms%nlo(ntyp) + lop
351 : !tlmplm%tuloulo_new(m,m,mlolo_new+loplo_new,iSpinPr,iSpin) = tlmplm%tuloulo_new(m,m,mlolo_new+loplo_new,iSpinPr,iSpin) &
352 : ! & + 0.5 * (enpara%ello0(lo,ntyp,iSpinPr) &
353 : ! & + enpara%ello0(lop,ntyp,iSpinPr)) &
354 : ! & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
355 : ! & + 0.5 * (usdus%ulouilopn(lo,lop,ntyp,iSpinPr) &
356 : ! & + usdus%ulouilopn(lop,lo,ntyp,iSpinPr))
357 : tlmplm%tuloulo_newer(m,m,lo,lop,ntyp,iSpinPr,iSpin) = tlmplm%tuloulo_newer(m,m,lo,lop,ntyp,iSpinPr,iSpin) &
358 : & + 0.5 * (enpara%ello0(lo,ntyp,iSpinPr) &
359 : & + enpara%ello0(lop,ntyp,iSpinPr)) &
360 : & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
361 : & + 0.5 * (usdus%ulouilopn(lo,lop,ntyp,iSpinPr) &
362 12 : & + usdus%ulouilopn(lop,lo,ntyp,iSpinPr))
363 : END IF
364 : ELSE
365 : !tlmplm%tuloulo(m,m,loplo+mlolo,iSpinPr,iSpin) = tlmplm%tuloulo(m,m,loplo+mlolo,iSpinPr,iSpin) &
366 : ! & + enpara%ello0(lo,ntyp,iSpinPr) &
367 : ! & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
368 : ! & + usdus%ulouilopn(lop,lo,ntyp,iSpinPr)
369 : tlmplm%tuloulo_newer(m,m,lop,lo,ntyp,iSpinPr,iSpin) = tlmplm%tuloulo_newer(m,m,lop,lo,ntyp,iSpinPr,iSpin) &
370 : & + enpara%ello0(lo,ntyp,iSpinPr) &
371 : & * usdus%uloulopn(lop,lo,ntyp,iSpinPr) &
372 : & + usdus%ulouilopn(lop,lo,ntyp,iSpinPr)
373 : IF (.NOT.lop.EQ.lo) THEN
374 : tlmplm%tuloulo_newer(m,m,lo,lop,ntyp,iSpinPr,iSpin) = tlmplm%tuloulo_newer(m,m,lo,lop,ntyp,iSpinPr,iSpin) &
375 : & + enpara%ello0(lop,ntyp,iSpinPr) &
376 : & * usdus%uloulopn(lo,lop,ntyp,iSpinPr) &
377 : & + usdus%ulouilopn(lo,lop,ntyp,iSpinPr)
378 : END IF
379 : END IF
380 : END DO
381 : END DO
382 : END DO
383 : END IF
384 :
385 1166 : END SUBROUTINE tlo
386 : END MODULE m_tlo
|
