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_cdnmtlo
8 : !! Archived comment:
9 : !!
10 : !! Add the local orbital contributions to the charge density. The
11 : !! corresponding summation of the pure apw contribuions is done in
12 : !! cdnval.
13 : !! Philipp Kurz 99/04
14 : CONTAINS
15 972 : SUBROUTINE cdnmtlo(itype,ilSpinPr,ilSpin,input,atoms,sphhar,sym,usdus,noco, &
16 972 : ello,vr,denCoeffs,f,g,rho,qmtllo,moments,rhoIm,f2,g2)
17 : !! Current situation:
18 : !!
19 : !! Renamed the routine from rhosphnlo to cdnmtlo and adjusted it to handle
20 : !! variables wrapped up in types while expanding the functionality to the
21 : !! previously missing case of offdiagonal density matrix elements. Hence,
22 : !! this routine can now calculate density contributions
23 : !! $$\rho_{L}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}(r)=
24 : !! \sum_{l',l,\lambda',\lambda,s}d_{l',l,L,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}
25 : !! u_{l',\lambda',s}^{\sigma_{\alpha}',\alpha}(r)u_{l,\lambda,s}^{\sigma_{\alpha},\alpha}(r)$$
26 : !! where one of the \(\lambda\) describes a local orbital. \(s\) is the
27 : !! index for the big/small components yielded by the scalar-relativistic
28 : !! Schrödinger equation.
29 : USE m_constants, ONLY : c_light,sfp_const
30 : USE m_types
31 : USE m_radsra
32 : USE m_radsrdn
33 :
34 : IMPLICIT NONE
35 :
36 : TYPE(t_input), INTENT(IN) :: input
37 : TYPE(t_sphhar), INTENT(IN) :: sphhar
38 : TYPE(t_atoms), INTENT(IN) :: atoms
39 : TYPE(t_sym), INTENT(IN) :: sym
40 : TYPE(t_usdus), INTENT(IN) :: usdus
41 : TYPE(t_noco), INTENT(IN) :: noco
42 : TYPE (t_denCoeffs), INTENT(IN) :: denCoeffs
43 :
44 : INTEGER, INTENT (IN) :: itype, ilSpinPr, ilSpin
45 : REAL, INTENT (IN) :: vr(:,:)
46 : REAL, INTENT (IN) :: ello(:,:)
47 : REAL, INTENT (IN) :: f(:,:,0:), g(:,:,0:)
48 : REAL, INTENT (INOUT) :: qmtllo(0:)
49 : REAL, INTENT (INOUT) :: rho(:,0:)
50 :
51 : REAL, OPTIONAL, INTENT(INOUT) :: rhoIm(:,0:)
52 : REAL, OPTIONAL, INTENT(IN) :: f2(:,:,0:), g2(:,:,0:)
53 :
54 : TYPE(t_moments), OPTIONAL, INTENT(INOUT) :: moments
55 :
56 : INTEGER, PARAMETER :: lcf=3
57 :
58 : INTEGER :: j,l,lh,lo,lop,lp,nodedum,llp,iSpin,jsp_start,jsp_end
59 : REAL :: dsdum,usdum,c_1,c_2,dus,ddn,c
60 : COMPLEX :: ctemp
61 :
62 972 : REAL :: filo(atoms%jmtd,2)
63 :
64 972 : REAL, ALLOCATABLE :: flo(:,:,:,:),glo(:,:)
65 972 : REAL, ALLOCATABLE :: fPr(:,:,:),gPr(:,:,:)
66 :
67 972 : c = c_light(1.0)
68 972 : c_1 = 1.0 / atoms%neq(itype)
69 972 : c_2 = 1.0 /(atoms%neq(itype)*sfp_const)
70 :
71 3888 : ALLOCATE(fPr(atoms%jmtd,2,0:atoms%lmaxd))
72 2916 : ALLOCATE(gPr(atoms%jmtd,2,0:atoms%lmaxd))
73 :
74 972 : IF (PRESENT(f2)) THEN
75 507329 : fPr = f2
76 507329 : gPr = g2
77 : ELSE
78 12901413 : fPr = f
79 12901413 : gPr = g
80 : END IF
81 :
82 972 : IF (ilSpinPr==ilSpin) THEN
83 1784 : DO lo = 1,atoms%nlo(itype)
84 849 : l = atoms%llo(lo,itype)
85 : ctemp = (denCoeffs%mt_ulo_coeff(lo,itype,0,ilSpinPr,ilSpin)+denCoeffs%mt_lou_coeff(lo,itype,0,ilSpinPr,ilSpin)) &
86 : & * usdus%uulon(lo,itype,ilSpin) &
87 : & + (denCoeffs%mt_ulo_coeff(lo,itype,1,ilSpinPr,ilSpin)+denCoeffs%mt_lou_coeff(lo,itype,1,ilSpinPr,ilSpin)) &
88 849 : & * usdus%dulon(lo,itype,ilSpin)
89 849 : qmtllo(l) = qmtllo(l) + REAL(ctemp) * c_1
90 3305 : DO lop = 1,atoms%nlo(itype)
91 2370 : IF (atoms%llo(lop,itype).EQ.l) THEN
92 853 : ctemp = denCoeffs%mt_lolo_coeff(lop,lo,itype,ilSpinPr,ilSpin) * usdus%uloulopn(lop,lo,itype,ilSpin)
93 853 : qmtllo(l) = qmtllo(l) + REAL(ctemp) * c_1
94 : END IF
95 : END DO
96 : END DO
97 : END IF
98 :
99 972 : jsp_start = MERGE(ilSpin,1,ilSpinPr==ilSpin)
100 972 : jsp_end = MERGE(ilSpin,2,ilSpinPr==ilSpin)
101 :
102 972 : IF (noco%l_mperp) THEN
103 555 : ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,input%jspins) )
104 : ELSE
105 4305 : ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,jsp_start:jsp_end) )
106 : END IF
107 :
108 2916 : ALLOCATE ( glo(atoms%jmtd,2) )
109 :
110 : ! Calculate the local orbital radial functions
111 1981 : DO iSpin = jsp_start,jsp_end
112 2970 : DO lo = 1,atoms%nlo(itype)
113 989 : l = atoms%llo(lo,itype)
114 : CALL radsra(ello(lo,iSpin),l,vr(:,iSpin),atoms%rmsh(1,itype),atoms%dx(itype), &
115 989 : atoms%jri(itype),c,usdum,dus,nodedum,flo(:,1,lo,iSpin),flo(:,2,lo,iSpin))
116 :
117 1998 : IF (atoms%l_dulo(lo,itype).OR.atoms%ulo_der(lo,itype)>=1) THEN
118 : ! Calculate orthogonal energy derivative at E
119 2 : j = atoms%ulo_der(lo,itype)
120 2 : IF(atoms%l_dulo(lo,itype)) j = 1
121 : CALL radsrdn(ello(lo,iSpin),l,vr(:,iSpin),atoms%rmsh(1,itype),atoms%dx(itype), &
122 2 : atoms%jri(itype),c,usdum,dsdum,ddn,nodedum,glo,filo,flo(:,:,lo,iSpin),dus,j)
123 : ! filo is a dummy array
124 1516 : DO j=1,atoms%jri(itype)
125 1514 : flo(j,1,lo,iSpin) = glo(j,1)
126 1516 : flo(j,2,lo,iSpin) = glo(j,2)
127 : END DO
128 2 : ddn = sqrt(ddn)
129 2 : IF(atoms%l_dulo(lo,itype)) ddn=1.0
130 3034 : flo(:,:,lo,iSpin) = flo(:,:,lo,iSpin)/ddn ! Normalize ulo (flo) if APW+lo is not used
131 : END IF
132 : END DO
133 : END DO
134 :
135 : ! Add the contribution of LO-LAPW and LO-LO cross-terms to the spherical
136 : ! charge density inside the Muffin Tins.
137 1891 : DO lo = 1,atoms%nlo(itype)
138 919 : l = atoms%llo(lo,itype)
139 919 : llp = (l* (l+1))/2 + l
140 708866 : DO j = 1,atoms%jri(itype)
141 707947 : IF (.NOT.PRESENT(rhoIm)) THEN
142 : ! Base case for diagonal density components.
143 : ctemp = (denCoeffs%mt_ulo_coeff(lo,itype,0,ilSpinPr,ilSpin)+denCoeffs%mt_lou_coeff(lo,itype,0,ilSpinPr,ilSpin)) &
144 : * (fPr(j,1,l)*flo(j,1,lo,ilSpin)+fPr(j,2,l)*flo(j,2,lo,ilSpin)) &
145 : + (denCoeffs%mt_ulo_coeff(lo,itype,1,ilSpinPr,ilSpin)+denCoeffs%mt_lou_coeff(lo,itype,1,ilSpinPr,ilSpin)) &
146 655057 : * (gPr(j,1,l)*flo(j,1,lo,ilSpin)+gPr(j,2,l)*flo(j,2,lo,ilSpin))
147 655057 : rho(j,0) = rho(j,0) + c_2 * REAL(ctemp)
148 : ELSE
149 : ! If the local spins are not the same, the radial functions differ
150 : ! and we need to treat u-ulo and ulo-u terms separately and save
151 : ! the imaginary part as well.
152 : ctemp = denCoeffs%mt_ulo_coeff(lo,itype,0,ilSpinPr,ilSpin) &
153 : * (fPr(j,1,l)*flo(j,1,lo,ilSpin)+fPr(j,2,l)*flo(j,2,lo,ilSpin)) &
154 : + denCoeffs%mt_ulo_coeff(lo,itype,1,ilSpinPr,ilSpin) &
155 : * (gPr(j,1,l)*flo(j,1,lo,ilSpin)+gPr(j,2,l)*flo(j,2,lo,ilSpin)) &
156 : + denCoeffs%mt_lou_coeff(lo,itype,0,ilSpinPr,ilSpin) &
157 : * (f(j,1,l)*flo(j,1,lo,ilSpinPr)+f(j,2,l)*flo(j,2,lo,ilSpinPr)) &
158 : + denCoeffs%mt_lou_coeff(lo,itype,1,ilSpinPr,ilSpin) &
159 52890 : * (g(j,1,l)*flo(j,1,lo,ilSpinPr)+g(j,2,l)*flo(j,2,lo,ilSpinPr))
160 52890 : rho(j,0) = rho(j,0) + c_2 * REAL(ctemp)
161 52890 : rhoIm(j,0) = rhoIm(j,0) + c_2 * AIMAG(ctemp)
162 : END IF
163 708866 : IF (l.LE.input%lResMax.AND.PRESENT(moments)) THEN
164 0 : IF (ilSpinPr==ilSpin) THEN
165 0 : moments%rhoLRes(j,0,llp,itype,ilSpin) = moments%rhoLRes(j,0,llp,itype,ilSpin) + c_2 * REAL(ctemp)
166 : ELSE
167 0 : moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3) + c_2 * REAL(ctemp)
168 0 : moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4) + c_2 * AIMAG(ctemp)
169 : END IF
170 : END IF
171 : END DO
172 3552 : DO lop = 1,atoms%nlo(itype)
173 2580 : IF (atoms%llo(lop,itype)==l) THEN
174 711898 : DO j = 1,atoms%jri(itype)
175 : ctemp = c_2 * denCoeffs%mt_lolo_coeff(lop,lo,itype,ilSpinPr,ilSpin) &
176 710975 : * (flo(j,1,lop,ilSpinPr)*flo(j,1,lo,ilSpin)+flo(j,2,lop,ilSpinPr)*flo(j,2,lo,ilSpin))
177 710975 : rho(j,0) = rho(j,0) + REAL(ctemp)
178 710975 : IF (PRESENT(rhoIm)) rhoIm(j,0) = rhoIm(j,0) + AIMAG(ctemp)
179 711898 : IF (l<=input%lResMax.AND.PRESENT(moments)) THEN
180 0 : IF (ilSpinPr==ilSpin) THEN
181 0 : moments%rhoLRes(j,0,llp,itype,ilSpin) = moments%rhoLRes(j,0,llp,itype,ilSpin) + REAL(ctemp)
182 : ELSE
183 0 : moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3) + REAL(ctemp)
184 0 : moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4) + AIMAG(ctemp)
185 : END IF
186 : END IF
187 : END DO
188 : END IF
189 : END DO
190 : END DO
191 :
192 : ! Add the contribution of LO-LAPW and LO-LO cross-terms to the non-spherical
193 : ! charge density inside the Muffin Tins.
194 28208 : DO lh = 1,sphhar%nlh(sym%ntypsy(atoms%firstAtom(itype)))
195 283926 : DO lp = 0,atoms%lmax(itype)
196 588033 : DO lo = 1,atoms%nlo(itype)
197 304107 : l = atoms%llo(lo,itype)
198 :
199 : ! Exclude non-spherical contributions for CF in the diagonal case.
200 : IF (atoms%l_outputCFpot(itype).AND.atoms%l_outputCFremove4f(itype) &
201 304107 : .AND.(l==lcf.AND.lp==lcf).AND.ilSpinPr==ilSpin) CYCLE
202 :
203 304107 : llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
204 237823352 : DO j = 1,atoms%jri(itype)
205 237262555 : IF (.NOT.PRESENT(rhoIm)) THEN
206 : ctemp = c_1 &
207 : * ((denCoeffs%nmt_ulo_coeff(lp,lo,lh,itype,0,ilSpinPr,ilSpin)+denCoeffs%nmt_lou_coeff(lp,lo,lh,itype,0,ilSpinPr,ilSpin)) &
208 : * (fPr(j,1,lp)*flo(j,1,lo,ilSpin)+fPr(j,2,lp)*flo(j,2,lo,ilSpin)) &
209 : + (denCoeffs%nmt_ulo_coeff(lp,lo,lh,itype,1,ilSpinPr,ilSpin)+denCoeffs%nmt_lou_coeff(lp,lo,lh,itype,1,ilSpinPr,ilSpin)) &
210 201125449 : * (gPr(j,1,lp)*flo(j,1,lo,ilSpin)+gPr(j,2,lp)*flo(j,2,lo,ilSpin)))
211 201125449 : rho(j,lh) = rho(j,lh) + REAL(ctemp)
212 : ELSE
213 : ctemp = c_1 &
214 : * (denCoeffs%nmt_ulo_coeff(lp,lo,lh,itype,0,ilSpinPr,ilSpin) &
215 : * (fPr(j,1,lp)*flo(j,1,lo,ilSpin)+fPr(j,2,lp)*flo(j,2,lo,ilSpin)) &
216 : + denCoeffs%nmt_ulo_coeff(lp,lo,lh,itype,1,ilSpinPr,ilSpin) &
217 : * (gPr(j,1,lp)*flo(j,1,lo,ilSpin)+gPr(j,2,lp)*flo(j,2,lo,ilSpin)) &
218 : + denCoeffs%nmt_lou_coeff(lp,lo,lh,itype,0,ilSpinPr,ilSpin) &
219 : * (f(j,1,lp)*flo(j,1,lo,ilSpinPr)+f(j,2,lp)*flo(j,2,lo,ilSpinPr)) &
220 : + denCoeffs%nmt_lou_coeff(lp,lo,lh,itype,1,ilSpinPr,ilSpin) &
221 36137106 : * (g(j,1,lp)*flo(j,1,lo,ilSpinPr)+g(j,2,lp)*flo(j,2,lo,ilSpinPr)))
222 36137106 : rho(j,lh) = rho(j,lh) + REAL(ctemp)
223 36137106 : rhoIm(j,lh) = rhoIm(j,lh) + AIMAG(ctemp)
224 : END IF
225 237566662 : IF ((l.LE.input%lResMax).AND.(lp.LE.input%lResMax).AND.PRESENT(moments)) THEN
226 0 : IF (ilSpinPr==ilSpin) THEN
227 0 : moments%rhoLRes(j,lh,llp,itype,ilSpin) = moments%rhoLRes(j,lh,llp,itype,ilSpin) + REAL(ctemp)
228 : ELSE
229 0 : moments%rhoLRes(j,lh,llp,itype,3) = moments%rhoLRes(j,lh,llp,itype,3) + REAL(ctemp)
230 0 : moments%rhoLRes(j,lh,llp,itype,4) = moments%rhoLRes(j,lh,llp,itype,4) + AIMAG(ctemp)
231 : END IF
232 : END IF
233 : END DO
234 : END DO
235 : END DO
236 60453 : DO lo = 1,atoms%nlo(itype)
237 32245 : l = atoms%llo(lo,itype)
238 119740 : DO lop = 1,atoms%nlo(itype)
239 60259 : lp = atoms%llo(lop,itype)
240 60259 : llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
241 46593671 : DO j = 1,atoms%jri(itype)
242 : ctemp = c_1 * denCoeffs%nmt_lolo_coeff(lop,lo,lh,itype,ilSpinPr,ilSpin) &
243 46501167 : * (flo(j,1,lop,ilSpinPr)*flo(j,1,lo,ilSpin)+flo(j,2,lop,ilSpinPr)*flo(j,2,lo,ilSpin))
244 46501167 : rho(j,lh) = rho(j,lh) + REAL(ctemp)
245 46501167 : IF (PRESENT(rhoIm)) rhoIm(j,lh) = rhoIm(j,lh) + AIMAG(ctemp)
246 46561426 : IF ((l.LE.input%lResMax).AND.(lp.LE.input%lResMax).AND.PRESENT(moments)) THEN
247 0 : IF (ilSpinPr==ilSpin) THEN
248 0 : moments%rhoLRes(j,lh,llp,itype,ilSpin) = moments%rhoLRes(j,lh,llp,itype,ilSpin) + REAL(ctemp)
249 : ELSE
250 0 : moments%rhoLRes(j,lh,llp,itype,3) = moments%rhoLRes(j,lh,llp,itype,3) + REAL(ctemp)
251 0 : moments%rhoLRes(j,lh,llp,itype,4) = moments%rhoLRes(j,lh,llp,itype,4) + AIMAG(ctemp)
252 : END IF
253 : END IF
254 : END DO
255 : END DO
256 : END DO
257 : END DO
258 :
259 972 : DEALLOCATE (flo,glo)
260 :
261 972 : END SUBROUTINE cdnmtlo
262 : END MODULE m_cdnmtlo
|
