Line data Source code
1 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
2 : ! This file is part of FLEUR and available as free software under the conditions
3 : ! of the MIT license as expressed in the LICENSE file in more detail.
4 : !--------------------------------------------------------------------------------
5 :
6 : MODULE m_cdnmt
7 : !! Archived comment:
8 : !!
9 : !! This subroutine calculates the spherical and non-spherical charge-
10 : !! density and the orbital moment inside the muffin-tin spheres.
11 : !! Philipp Kurz 2000-02-03
12 : CONTAINS
13 :
14 509 : SUBROUTINE cdnmt(jspd,input,atoms,sym,sphhar,noco,jsp_start,jsp_end,enpara,banddos,&
15 509 : vr,denCoeffs,usdus,orb,denCoeffsOffdiag,rho,hub1inp,moments,jDOS,hub1data,rhoIm)
16 : !! Current situation:
17 : !!
18 : !! This routine calculates density contributions
19 : !! $$\rho_{L}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}(r)=
20 : !! \sum_{l',l,\lambda',\lambda,s}d_{l',l,L,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}
21 : !! u_{l',\lambda',s}^{\sigma_{\alpha}',\alpha}(r)u_{l,\lambda,s}^{\sigma_{\alpha},\alpha}(r)$$
22 : !! \(s\) is the index for the big/small components yielded by the
23 : !! scalar-relativistic Schrödinger equation.
24 :
25 : USE m_types
26 : USE m_constants
27 : USE m_cdnmtlo
28 : USE m_radfun
29 : USE m_orbmom2
30 : USE m_xmlOutput
31 : USE m_types_orbcomp
32 : USE m_types_jDOS
33 : USE m_types_mcd
34 :
35 : IMPLICIT NONE
36 :
37 : TYPE(t_input), INTENT(IN) :: input
38 : TYPE(t_usdus), INTENT(INOUT) :: usdus !in fact only the lo part is intent(in)
39 : TYPE(t_noco), INTENT(IN) :: noco
40 : TYPE(t_sphhar), INTENT(IN) :: sphhar
41 : TYPE(t_atoms), INTENT(IN) :: atoms
42 : TYPE(t_sym), INTENT(IN) :: sym
43 : TYPE(t_enpara), INTENT(IN) :: enpara
44 : TYPE(t_banddos), INTENT(IN) :: banddos
45 : TYPE(t_hub1inp), INTENT(IN) :: hub1inp
46 : TYPE (t_orb), INTENT(IN) :: orb
47 : TYPE (t_denCoeffs), INTENT(IN) :: denCoeffs
48 : TYPE (t_denCoeffsOffdiag), INTENT(IN) :: denCoeffsOffdiag
49 :
50 : TYPE(t_jDOS), OPTIONAL, INTENT(IN) :: jDOS
51 : TYPE(t_moments), OPTIONAL, INTENT(INOUT) :: moments
52 : TYPE(t_hub1data), OPTIONAL, INTENT(INOUT) :: hub1data
53 :
54 : INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
55 :
56 : REAL, INTENT (IN) :: vr(atoms%jmtd,atoms%ntype,jspd)
57 : REAL, INTENT (INOUT) :: rho(:,0:,:,:)!(toms%jmtd,0:sphhar%nlhd,atoms%ntype,jspd)
58 : REAL, OPTIONAL, INTENT(INOUT) :: rhoIm(:,0:,:,:)
59 :
60 : INTEGER, PARAMETER :: lcf=3
61 :
62 : INTEGER :: itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu,llpb,natom,jj,n_dos
63 : INTEGER :: ilo,ilop,i,i_hia,i_exc
64 : REAL :: wronk,qmtt
65 : COMPLEX :: cs, rho21
66 : LOGICAL :: l_hia,l_performSpinavg
67 :
68 509 : REAL :: qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd),vrTmp(atoms%jmtd)
69 509 : REAL,ALLOCATABLE :: vr0(:,:)
70 : CHARACTER(LEN=20) :: attributes(6)
71 :
72 509 : REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:)
73 :
74 509 : IF (noco%l_mperp) THEN
75 29 : IF (denCoeffsOffdiag%l_fmpl) THEN
76 4297649 : rho(:,:,:,3:4) = CMPLX(0.0,0.0)
77 : END IF
78 : END IF
79 :
80 509 : l_performSpinavg = .FALSE.
81 509 : IF(PRESENT(hub1data)) l_performSpinavg = hub1data%l_performSpinavg
82 :
83 17404 : qmtl = 0
84 : !$OMP PARALLEL DEFAULT(none) &
85 : !$OMP SHARED(usdus,rho,moments,qmtl,hub1inp,hub1data,rhoIm) &
86 : !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar,l_performSpinavg)&
87 : !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd,input,sym)&
88 : !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,qmtllo,qmtt,nd,lh,lp,llp,llpb,cs)&
89 509 : !$OMP PRIVATE(l_hia,vrTmp,vr0)
90 : IF (noco%l_mperp) THEN
91 : ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
92 : ELSE
93 : ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
94 : ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
95 : END IF
96 :
97 : !$OMP DO
98 : DO itype = 1,atoms%ntype
99 : if (atoms%l_nonpolbas(itype)) THEN
100 : if (.not.allocated(vr0)) allocate(vr0(atoms%jmtd,jspd))
101 : vr0(:,1)=(vr(:,itype,1)+vr(:,itype,2))/2
102 : vr0(:,2)=vr0(:,1)
103 : else
104 : vr0=vr(:,itype,:)
105 : ENDIF
106 :
107 : na = atoms%firstAtom(itype)
108 :
109 : DO ispin = jsp_start,jsp_end
110 : !Spherical component
111 : CALL timestart("cdnmt spherical diagonal")
112 : DO l = 0,atoms%lmax(itype)
113 : !Check if the orbital is treated with Hubbard 1
114 : l_hia=.FALSE.
115 : DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
116 : IF(atoms%lda_u(i)%atomType==itype.AND.atoms%lda_u(i)%l==l) THEN
117 : l_hia=.TRUE.
118 : END IF
119 : END DO
120 :
121 : !In the case of a spin-polarized calculation with Hubbard 1 we want to treat
122 : !the correlated orbitals with a non-spin-polarized basis
123 : IF(l_hia.AND.jspd==2 .AND. l_performSpinavg) THEN
124 : vrTmp = (vr0(:,1) + vr0(:,2))/2.0
125 : ELSE
126 : vrTmp = vr0(:,ispin)
127 : END IF
128 :
129 : CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vrTmp,atoms,&
130 : f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
131 : llp = (l*(l+1))/2 + l
132 :
133 : DO j = 1,atoms%jri(itype)
134 : cs = denCoeffs%mt_coeff(l,itype,0,0,ispin,ispin)*(f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin)) &
135 : + denCoeffs%mt_coeff(l,itype,0,1,ispin,ispin)*(f(j,1,l,ispin)*g(j,1,l,ispin)+f(j,2,l,ispin)*g(j,2,l,ispin)) &
136 : + denCoeffs%mt_coeff(l,itype,1,0,ispin,ispin)*(g(j,1,l,ispin)*f(j,1,l,ispin)+g(j,2,l,ispin)*f(j,2,l,ispin)) &
137 : + denCoeffs%mt_coeff(l,itype,1,1,ispin,ispin)*(g(j,1,l,ispin)*g(j,1,l,ispin)+g(j,2,l,ispin)*g(j,2,l,ispin))
138 : rho(j,0,itype,ispin) = rho(j,0,itype,ispin) + REAL(cs)/(atoms%neq(itype)*sfp_const)
139 : IF (l<=input%lResMax.AND.PRESENT(moments)) THEN !DFT case
140 : moments%rhoLRes(j,0,llp,itype,ispin) = moments%rhoLRes(j,0,llp,itype,ispin) + REAL(cs)/(atoms%neq(itype)*sfp_const)
141 : ELSE IF (.NOT.PRESENT(moments)) THEN
142 : rhoIm(j,0,itype,ispin) = rhoIm(j,0,itype,ispin) + AIMAG(cs)/(atoms%neq(itype)*sfp_const)
143 : END IF
144 : IF(PRESENT(hub1data).AND.l<=lmaxU_const) THEN
145 : hub1data%cdn_atomic(j,l,itype,ispin) = hub1data%cdn_atomic(j,l,itype,ispin) &
146 : + REAL(denCoeffs%mt_coeff(l,itype,0,0,ispin,ispin)) &
147 : * (f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin)) &
148 : * 1.0/(atoms%neq(itype)*sfp_const)
149 : END IF
150 : END DO
151 : END DO
152 : CALL timestop("cdnmt spherical diagonal")
153 :
154 : !Add the contribution of LO-LO and LAPW-LO cross-terms to rho and
155 : !qmtl. The latter are stored in qmtllo.
156 : DO l = 0,atoms%lmaxd
157 : qmtllo(l) = 0.0
158 : END DO
159 : IF (PRESENT(moments)) THEN !DFT case
160 : CALL timestart("cdnmt LO diagonal")
161 : CALL cdnmtlo(itype,ispin,ispin,input,atoms,sphhar,sym,usdus,noco,&
162 : enpara%ello0(:,itype,:),vr0(:,:),denCoeffs,&
163 : f(:,:,0:,ispin),g(:,:,0:,ispin),&
164 : rho(:,0:,itype,ispin),qmtllo,moments=moments)
165 : CALL timestop("cdnmt LO diagonal")
166 :
167 : !l-decomposed density for each atom type
168 : qmtt = 0.0
169 : DO l = 0,atoms%lmax(itype)
170 : qmtl(l,ispin,itype) = REAL(denCoeffs%mt_coeff(l,itype,0,0,ispin,ispin)+denCoeffs%mt_coeff(l,itype,1,1,ispin,ispin) &
171 : * usdus%ddn(l,itype,ispin))/atoms%neq(itype) + qmtllo(l)
172 : qmtt = qmtt + qmtl(l,ispin,itype)
173 : END DO
174 : moments%chmom(itype,ispin) = qmtt
175 : ELSE !DFPT case
176 : CALL timestart("cdnmt LO diagonal")
177 : CALL cdnmtlo(itype,ispin,ispin,input,atoms,sphhar,sym,usdus,noco,&
178 : enpara%ello0(:,itype,:),vr0(:,:),denCoeffs,&
179 : f(:,:,0:,ispin),g(:,:,0:,ispin),&
180 : rho(:,0:,itype,ispin),qmtllo,&
181 : rhoIm=rhoIm(:,0:,itype,ispin), f2=f(:,:,0:,ispin), g2=g(:,:,0:,ispin))
182 : CALL timestop("cdnmt LO diagonal")
183 : END IF
184 :
185 : !Get the magnetic moment for the shells where we defined additional exchange splittings for DFT+Hubbard 1
186 : IF(PRESENT(hub1data)) THEN
187 : DO i_hia = 1, atoms%n_hia
188 : IF(atoms%lda_u(atoms%n_u+i_hia)%atomType/=itype) CYCLE
189 : DO i_exc = 1, hub1inp%n_exc(i_hia)
190 : hub1data%mag_mom(i_hia,i_exc) = hub1data%mag_mom(i_hia,i_exc) + (-1)**(ispin-1) &
191 : * qmtl(hub1inp%exc_l(i_hia,i_exc),ispin,itype)
192 : END DO
193 : END DO
194 : END IF
195 :
196 : !Spherical angular component for the SOC contribution
197 : IF (noco%l_soc.AND.PRESENT(moments)) THEN
198 : CALL orbmom2(atoms,itype,ispin,usdus%ddn(0:,itype,ispin),&
199 : orb,usdus%uulon(:,itype,ispin),usdus%dulon(:,itype,ispin),&
200 : usdus%uloulopn(:,:,itype,ispin),moments%clmom(:,itype,ispin))!keep
201 : END IF
202 :
203 : !Non-spherical components
204 : CALL timestart("cdnmt non-spherical diagonal")
205 : nd = sym%ntypsy(na)
206 : DO lh = 1,sphhar%nlh(nd)
207 : DO l = 0,atoms%lmax(itype)
208 : DO lp = 0,MERGE(l,atoms%lmax(itype),PRESENT(moments))
209 : llp = (l* (l+1))/2 + lp
210 : IF (.NOT.PRESENT(moments)) llp = lp*(atoms%lmax(itype)+1)+l
211 : IF(atoms%l_outputCFpot(itype).AND.atoms%l_outputCFremove4f(itype)&
212 : .AND.(l==lcf.AND.lp==lcf)) CYCLE !Exclude non-spherical contributions for CF
213 :
214 : DO j = 1,atoms%jri(itype)
215 : cs = 0.0
216 : cs = denCoeffs%nmt_coeff(llp,lh,itype,0,0,ispin,ispin)*(f(j,1,lp,ispin)*f(j,1,l,ispin)+ f(j,2,lp,ispin)*f(j,2,l,ispin)) &
217 : + denCoeffs%nmt_coeff(llp,lh,itype,0,1,ispin,ispin)*(f(j,1,lp,ispin)*g(j,1,l,ispin)+ f(j,2,lp,ispin)*g(j,2,l,ispin)) &
218 : + denCoeffs%nmt_coeff(llp,lh,itype,1,0,ispin,ispin)*(g(j,1,lp,ispin)*f(j,1,l,ispin)+ g(j,2,lp,ispin)*f(j,2,l,ispin)) &
219 : + denCoeffs%nmt_coeff(llp,lh,itype,1,1,ispin,ispin)*(g(j,1,lp,ispin)*g(j,1,l,ispin)+ g(j,2,lp,ispin)*g(j,2,l,ispin))
220 : rho(j,lh,itype,ispin) = rho(j,lh,itype,ispin)+ REAL(cs)/atoms%neq(itype)
221 : IF ((l<=input%lResMax).AND.(lp<=input%lResMax).AND.PRESENT(moments)) THEN !DFT case
222 : moments%rhoLRes(j,lh,llp,itype,ispin) = moments%rhoLRes(j,lh,llp,itype,ispin) + REAL(cs)/atoms%neq(itype)
223 : ELSE IF (.NOT.PRESENT(moments)) THEN
224 : rhoIm(j,lh,itype,ispin) = rhoIm(j,lh,itype,ispin)+ AIMAG(cs)/atoms%neq(itype)
225 : END IF
226 : END DO
227 : END DO
228 : END DO
229 : END DO
230 : CALL timestop("cdnmt non-spherical diagonal")
231 : END DO ! end of spin loop (ispin = jsp_start,jsp_end)
232 :
233 : IF (noco%l_mperp) THEN
234 : IF (PRESENT(moments)) THEN
235 : !Calculate the off-diagonal integrated density
236 : DO l = 0, atoms%lmax(itype)
237 : moments%qa21(itype) = moments%qa21(itype) + CONJG( &
238 : denCoeffsOffdiag%uu21(l,itype) * denCoeffsOffdiag%uu21n(l,itype) + &
239 : denCoeffsOffdiag%ud21(l,itype) * denCoeffsOffdiag%ud21n(l,itype) + &
240 : denCoeffsOffdiag%du21(l,itype) * denCoeffsOffdiag%du21n(l,itype) + &
241 : denCoeffsOffdiag%dd21(l,itype) * denCoeffsOffdiag%dd21n(l,itype) )/atoms%neq(itype)
242 : END DO
243 : DO ilo = 1, atoms%nlo(itype)
244 : moments%qa21(itype) = moments%qa21(itype) + CONJG( &
245 : denCoeffsOffdiag%ulou21(ilo,itype) * denCoeffsOffdiag%ulou21n(ilo,itype) + &
246 : denCoeffsOffdiag%ulod21(ilo,itype) * denCoeffsOffdiag%ulod21n(ilo,itype) + &
247 : denCoeffsOffdiag%uulo21(ilo,itype) * denCoeffsOffdiag%uulo21n(ilo,itype) + &
248 : denCoeffsOffdiag%dulo21(ilo,itype) * denCoeffsOffdiag%dulo21n(ilo,itype) )/atoms%neq(itype)
249 : DO ilop = 1, atoms%nlo(itype)
250 : moments%qa21(itype) = moments%qa21(itype) + CONJG( &
251 : denCoeffsOffdiag%uloulop21(ilo,ilop,itype) * &
252 : denCoeffsOffdiag%uloulop21n(ilo,ilop,itype) )/atoms%neq(itype)
253 : END DO
254 : END DO
255 : END IF
256 :
257 : !The following part can be used to calculate the full spherical
258 : !and non-spherical parts of the off-diagonal magnetization
259 : !density.
260 : IF (denCoeffsOffdiag%l_fmpl) THEN
261 : !Spherical components for the off-diagonal density
262 : CALL timestart("cdnmt spherical off-diagonal")
263 : DO l = 0,atoms%lmax(itype)
264 : llp = (l* (l+1))/2 + l
265 : DO j = 1, atoms%jri(itype)
266 : cs = denCoeffs%mt_coeff(l,itype,0,0,2,1)*(f(j,1,l,2)*f(j,1,l,1)+f(j,2,l,2)*f(j,2,l,1)) &
267 : + denCoeffs%mt_coeff(l,itype,0,1,2,1)*(f(j,1,l,2)*g(j,1,l,1)+f(j,2,l,2)*g(j,2,l,1)) &
268 : + denCoeffs%mt_coeff(l,itype,1,0,2,1)*(g(j,1,l,2)*f(j,1,l,1)+g(j,2,l,2)*f(j,2,l,1)) &
269 : + denCoeffs%mt_coeff(l,itype,1,1,2,1)*(g(j,1,l,2)*g(j,1,l,1)+g(j,2,l,2)*g(j,2,l,1))
270 : rho21 = cs/(atoms%neq(itype)*sfp_const)
271 : rho(j,0,itype,3) = rho(j,0,itype,3) + REAL(rho21)
272 : IF (PRESENT(moments)) THEN
273 : rho(j,0,itype,4) = rho(j,0,itype,4) + AIMAG(rho21)
274 : ELSE
275 : rhoIm(j,0,itype,3) = rhoIm(j,0,itype,3) + AIMAG(rho21)
276 : cs = denCoeffs%mt_coeff(l,itype,0,0,1,2)*(f(j,1,l,1)*f(j,1,l,2)+f(j,2,l,1)*f(j,2,l,2)) &
277 : + denCoeffs%mt_coeff(l,itype,0,1,1,2)*(f(j,1,l,1)*g(j,1,l,2)+f(j,2,l,1)*g(j,2,l,2)) &
278 : + denCoeffs%mt_coeff(l,itype,1,0,1,2)*(g(j,1,l,1)*f(j,1,l,2)+g(j,2,l,1)*f(j,2,l,2)) &
279 : + denCoeffs%mt_coeff(l,itype,1,1,1,2)*(g(j,1,l,1)*g(j,1,l,2)+g(j,2,l,1)*g(j,2,l,2))
280 : rho21 = cs/(atoms%neq(itype)*sfp_const)
281 : rho(j,0,itype,4) = rho(j,0,itype,4) + REAL(rho21)
282 : rhoIm(j,0,itype,4) = rhoIm(j,0,itype,4) + AIMAG(rho21)
283 : END IF
284 : IF (l<=input%lResMax.AND.PRESENT(moments)) THEN
285 : moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3)+ REAL(cs/(atoms%neq(itype)*sfp_const))
286 : moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4)+ AIMAG(cs/(atoms%neq(itype)*sfp_const))
287 : END IF
288 : END DO
289 : END DO
290 : CALL timestop("cdnmt spherical off-diagonal")
291 :
292 : !New feature: LOs for the offdiagonal density.
293 : !Add the contribution of LO-LO and LAPW-LO cross-terms to rho for
294 : !the offdiagonal magnetism.
295 : IF (PRESENT(moments)) THEN !DFT case
296 : CALL timestart("cdnmt LO off-diagonal")
297 : CALL cdnmtlo(itype,2,1,input,atoms,sphhar,sym,usdus,noco,&
298 : enpara%ello0(:,itype,:),vr0(:,:),denCoeffs,&
299 : f(:,:,0:,1),g(:,:,0:,1),&
300 : rho(:,0:,itype,3),qmtllo,moments=moments,&
301 : rhoIm=rho(:,0:,itype,4), f2=f(:,:,0:,2), g2=g(:,:,0:,2))
302 : !Note: qmtllo is irrelevant here
303 : CALL timestop("cdnmt LO off-diagonal")
304 : ELSE
305 : CALL timestart("cdnmt LO off-diagonal")
306 : CALL cdnmtlo(itype,2,1,input,atoms,sphhar,sym,usdus,noco,&
307 : enpara%ello0(:,itype,:),vr0(:,:),denCoeffs,&
308 : f(:,:,0:,1),g(:,:,0:,1),&
309 : rho(:,0:,itype,3),qmtllo,&
310 : rhoIm=rhoIm(:,0:,itype,3), f2=f(:,:,0:,2), g2=g(:,:,0:,2))
311 : CALL cdnmtlo(itype,1,2,input,atoms,sphhar,sym,usdus,noco,&
312 : enpara%ello0(:,itype,:),vr0(:,:),denCoeffs,&
313 : f(:,:,0:,2),g(:,:,0:,2),&
314 : rho(:,0:,itype,4),qmtllo,&
315 : rhoIm=rhoIm(:,0:,itype,4), f2=f(:,:,0:,1), g2=g(:,:,0:,1))
316 : !Note: qmtllo is irrelevant here
317 : CALL timestop("cdnmt LO off-diagonal")
318 : END IF
319 :
320 : !Non-spherical components for the off-diagonal density
321 : CALL timestart("cdnmt non-spherical off-diagonal")
322 : nd = sym%ntypsy(na)
323 : DO lh = 1,sphhar%nlh(nd)
324 : DO l = 0,atoms%lmax(itype)
325 : DO lp = 0,atoms%lmax(itype)
326 : llp = lp*(atoms%lmax(itype)+1)+l
327 : llpb = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
328 : DO j = 1,atoms%jri(itype)
329 : cs = denCoeffs%nmt_coeff(llp,lh,itype,0,0,2,1)*(f(j,1,lp,2)*f(j,1,l,1)+f(j,2,lp,2)*f(j,2,l,1)) &
330 : + denCoeffs%nmt_coeff(llp,lh,itype,0,1,2,1)*(f(j,1,lp,2)*g(j,1,l,1)+f(j,2,lp,2)*g(j,2,l,1)) &
331 : + denCoeffs%nmt_coeff(llp,lh,itype,1,0,2,1)*(g(j,1,lp,2)*f(j,1,l,1)+g(j,2,lp,2)*f(j,2,l,1)) &
332 : + denCoeffs%nmt_coeff(llp,lh,itype,1,1,2,1)*(g(j,1,lp,2)*g(j,1,l,1)+g(j,2,lp,2)*g(j,2,l,1))
333 : rho21 = cs/atoms%neq(itype)
334 : rho(j,lh,itype,3) = rho(j,lh,itype,3) + REAL(rho21)
335 : rho(j,lh,itype,4) = rho(j,lh,itype,4) + AIMAG(rho21)
336 : IF ((l<=input%lResMax).AND.(lp<=input%lResMax).AND.PRESENT(moments)) THEN
337 : moments%rhoLRes(j,lh,llpb,itype,3)= moments%rhoLRes(j,lh,llpb,itype,3) + REAL(cs/atoms%neq(itype))
338 : moments%rhoLRes(j,lh,llpb,itype,4)= moments%rhoLRes(j,lh,llpb,itype,4) + AIMAG(cs/atoms%neq(itype))
339 : END IF
340 : END DO
341 : END DO
342 : END DO
343 : END DO
344 : CALL timestop("cdnmt non-spherical off-diagonal")
345 : END IF ! denCoeffsOffdiag%l_fmpl
346 : END IF ! noco%l_mperp
347 : END DO ! end of loop over atom types
348 : !$OMP END DO
349 : DEALLOCATE (f,g)
350 : !$OMP END PARALLEL
351 :
352 509 : IF (.NOT.PRESENT(moments)) RETURN
353 : !if (size(rho,4)>3) rho(:,:,:,4)=-rho(:,:,:,4)
354 509 : WRITE (oUnit,FMT=8000)
355 : 8000 FORMAT (/,5x,'l-like charge',/,t6,'atom',t15,'s',t24,'p',&
356 : t33,'d',t42,'f',t51,'total')
357 :
358 1406 : DO itype = 1,atoms%ntype
359 2341 : DO ispin = jsp_start,jsp_end
360 935 : WRITE (oUnit,FMT=8100) itype, (qmtl(l,ispin,itype),l=0,3),moments%chmom(itype,ispin)
361 : 8100 FORMAT (' -->',i3,2x,4f9.5,2x,f9.5)
362 6545 : attributes = ''
363 935 : WRITE(attributes(1),'(i0)') itype
364 935 : WRITE(attributes(2),'(f12.7)') moments%chmom(itype,ispin)
365 935 : WRITE(attributes(3),'(f12.7)') qmtl(0,ispin,itype)
366 935 : WRITE(attributes(4),'(f12.7)') qmtl(1,ispin,itype)
367 935 : WRITE(attributes(5),'(f12.7)') qmtl(2,ispin,itype)
368 935 : WRITE(attributes(6),'(f12.7)') qmtl(3,ispin,itype)
369 : CALL writeXMLElementForm('mtCharge',(/'atomType','total ','s ','p ','d ','f '/),attributes(:6),&
370 7442 : reshape((/8,5,1,1,1,1,6,12,12,12,12,12/),(/6,2/)))
371 : END DO
372 : END DO
373 :
374 509 : IF(banddos%l_jDOS) THEN
375 1 : IF(PRESENT(jDOS)) THEN
376 1 : WRITE(oUnit,8200)
377 : 8200 FORMAT(/,5x,'j-decomposed charge',/,t6,'atom',t15,'s',t24,'p1/2',t33,'p3/2',&
378 : t42,'d3/2',t51,'d5/2',t60,'f5/2',t69,'f7/2')
379 2 : DO itype = 1, atoms%ntype
380 1 : natom = atoms%firstAtom(itype)
381 1 : IF (.NOT.banddos%dos_atom(natom)) CYCLE
382 : !Find index for dos.
383 1 : DO n_dos = 1, size(banddos%dos_atomlist)
384 1 : IF (banddos%dos_atomlist(n_dos)==natom) EXIT
385 : END DO
386 :
387 10 : WRITE(oUnit,8300) itype, jDOS%occ(0,1,n_dos), ((jDOS%occ(l,jj,n_dos),jj = 1, 2),l = 1, 3)
388 : 8300 FORMAT(' -->',i3,2x,f9.5,2x,6f9.5,/)
389 :
390 3 : CALL openXMLElementPoly('mtJcharge',['atomType'],[itype])
391 :
392 7 : attributes = ''
393 1 : WRITE(attributes(1),'(f12.7)') jDOS%occ(1,1,n_dos)
394 1 : WRITE(attributes(2),'(f12.7)') jDOS%occ(2,1,n_dos)
395 1 : WRITE(attributes(3),'(f12.7)') jDOS%occ(3,1,n_dos)
396 4 : CALL writeXMLElementForm('lowJ',['p','d','f'],attributes(:3),reshape([1,1,1,12,12,12],[3,2]))
397 :
398 7 : attributes = ''
399 1 : WRITE(attributes(1),'(f12.7)') jDOS%occ(1,2,n_dos)
400 1 : WRITE(attributes(2),'(f12.7)') jDOS%occ(2,2,n_dos)
401 1 : WRITE(attributes(3),'(f12.7)') jDOS%occ(3,2,n_dos)
402 4 : CALL writeXMLElementForm('highJ',['p','d','f'],attributes(:3),reshape([1,1,1,12,12,12],[3,2]))
403 :
404 2 : CALL closeXMLElement('mtJcharge')
405 :
406 : END DO
407 : END IF
408 : END IF
409 509 : END SUBROUTINE cdnmt
410 : END MODULE m_cdnmt
|
