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 : MODULE m_sorad
7 : USE m_juDFT
8 : !*********************************************************************
9 : ! generates radial spin-orbit matrix elements
10 : !*********************************************************************
11 : CONTAINS
12 204 : SUBROUTINE sorad(atoms,input,ntyp,vr,enpara,spav,rsoc,usdus,hub1data)
13 :
14 : USE m_constants, ONLY : c_light
15 : USE m_intgr, ONLY : intgr0
16 : USE m_sointg
17 : USE m_radsra
18 : USE m_radsrd
19 : USE m_radsrdn
20 : USE m_types
21 : IMPLICIT NONE
22 : TYPE(t_enpara),INTENT(IN) :: enpara
23 : TYPE(t_input),INTENT(IN) :: input
24 : TYPE(t_atoms),INTENT(IN) :: atoms
25 : TYPE(t_usdus),INTENT(INOUT) :: usdus
26 : TYPE(t_rsoc),INTENT(INOUT) :: rsoc
27 : TYPE(t_hub1data),OPTIONAL,INTENT(IN) :: hub1data
28 : ! ..
29 : ! .. Scalar Arguments ..
30 : INTEGER, INTENT (IN) :: ntyp
31 : LOGICAL, INTENT (IN) :: spav ! if T, spin-averaged pot is used
32 : ! ..
33 : ! .. Array Arguments ..
34 : REAL, INTENT (IN) :: vr(:,:)!(atoms%jmtd,input%jspins),
35 : ! ..
36 : ! .. Local Scalars ..
37 : REAL ddn1,e ,ulops,dulops,duds1
38 : INTEGER i,j,ir,jspin,l,noded,nodeu,ilo,ilop
39 : LOGICAL l_hia
40 : ! ..
41 : ! .. Local Arrays ..
42 204 : REAL, ALLOCATABLE :: p(:,:),pd(:,:),q(:,:),qd(:,:),plo(:,:)
43 204 : REAL, ALLOCATABLE :: plop(:,:),glo(:,:),fint(:),pqlo(:,:)
44 204 : REAL, ALLOCATABLE :: filo(:,:)
45 204 : REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:),vrTmp(:)
46 : ! ..
47 :
48 204 : IF (atoms%jri(ntyp)>atoms%jmtd) CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
49 : ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
50 2856 : & qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2),vrTmp(atoms%jmtd) )
51 :
52 1322712 : p = 0.0 ; pd = 0.0 ; q = 0.0 ; plo = 0.0 ; fint = 0.0
53 1469476 : qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0; vrTmp = 0.0
54 :
55 :
56 2144 : DO l = 0,atoms%lmax(ntyp)
57 1940 : l_hia=.FALSE.
58 1940 : DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
59 1940 : IF(atoms%lda_u(i)%atomType.EQ.ntyp.AND.atoms%lda_u(i)%l.EQ.l) THEN
60 0 : l_hia=.TRUE.
61 : ENDIF
62 : ENDDO
63 5374 : DO jspin = 1,input%jspins
64 : !TODO: here genMTBasis should be used
65 2432390 : vrTmp = vr(:,jspin)
66 221738 : if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
67 3434 : IF(l_hia.AND.input%jspins.EQ.2) THEN
68 0 : IF(PRESENT(hub1data)) THEN
69 0 : IF(hub1data%l_performSpinavg) vrTmp = (vr(:,1)+vr(:,2))/2.0
70 : ENDIF
71 : ENDIF
72 : !
73 : !---> calculate normalized function at e: p and q
74 : !
75 3434 : e = enpara%el0(l,ntyp,jspin)
76 : CALL radsra(&
77 : e,l,vrTmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
78 : usdus%us(l,ntyp,jspin),usdus%dus(l,ntyp,jspin),&
79 3434 : nodeu,p(:,jspin),q(:,jspin))
80 : !
81 : !---> calculate orthogonal energy derivative at e : pd and qd
82 : !
83 : CALL radsrd(&
84 : e,l,vrTmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
85 : usdus%uds(l,ntyp,jspin),usdus%duds(l,ntyp,jspin),&
86 : usdus%ddn(l,ntyp,jspin),noded,pd(:,jspin),qd(:,jspin),&
87 5374 : p(:,jspin),q(:,jspin),usdus%dus(l,ntyp,jspin))
88 :
89 : END DO ! end of spin loop
90 : !
91 : !---> in case of jspins=1
92 : !
93 :
94 1940 : IF (input%jspins.EQ.1) THEN
95 382036 : DO i = 1,atoms%jri(ntyp)
96 381590 : p(i,2) = p(i,1)
97 382036 : pd(i,2) = pd(i,1)
98 : ENDDO
99 : ENDIF
100 : !
101 : !---> common spin-orbit integrant V (average spin directions)
102 : ! SO
103 1416360 : v0(:) = 0.0
104 1940 : IF (input%jspins.EQ.1) THEN
105 382036 : v0(1:atoms%jri(ntyp)) = vr(1:atoms%jri(ntyp),1)
106 446 : e = enpara%el0(l,ntyp,1)
107 : ELSE
108 1012596 : DO i = 1,atoms%jri(ntyp)
109 1012596 : v0(i) = (vr(i,1)+vr(i,input%jspins))/2.
110 : END DO
111 1494 : e = (enpara%el0(l,ntyp,1)+enpara%el0(l,ntyp,input%jspins))/2.
112 : END IF
113 :
114 2682112 : CALL sointg(ntyp,e,vr,v0,atoms,input,vso)
115 1940 : IF (spav) THEN
116 109152 : DO i= 1,atoms%jmtd
117 109008 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
118 109152 : vso(i,2)= vso(i,1)
119 : ENDDO
120 : ENDIF
121 :
122 : ! s s' .s s'
123 : !--> radial integrals <u |V |u > = rsopp, <u |V |u > = rsopdp etc.
124 : ! SO SO
125 :
126 1940 : IF (l.GT.0) THEN ! there is no spin-orbit for s-states
127 5208 : DO i = 1, 2
128 12152 : DO j = 1, 2
129 5009248 : rsoc%rsopp(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
130 5009248 : rsoc%rsopdp(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
131 5009248 : rsoc%rsoppd(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
132 5012720 : rsoc%rsopdpd(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),(vso(:atoms%jri(ntyp),i)+vso(:atoms%jri(ntyp),j))*0.5,atoms%dx(ntyp),atoms%rmsh(1,ntyp))
133 : ENDDO
134 : ENDDO
135 : ENDIF ! l>0
136 : !
137 : !---> Check for local orbitals with same l
138 : !
139 4044 : DO ilo = 1, atoms%nlo(ntyp)
140 3840 : IF (atoms%llo(ilo,ntyp).EQ.l) THEN
141 :
142 550 : DO jspin = 1,input%jspins
143 281150 : vrTmp = vr(:,jspin)
144 48866 : if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
145 354 : e = enpara%ello0(ilo,ntyp,jspin)
146 : CALL radsra(&
147 : e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
148 : usdus%ulos(ilo,ntyp,jspin),usdus%dulos(ilo,ntyp,jspin),&
149 354 : nodeu,plo(:,jspin),qlo(:,jspin))
150 :
151 : !+apw+lo
152 550 : IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate energy derivative (of order atoms%ulo_der) at e
153 0 : ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
154 0 : glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
155 0 : pqlo(1:atoms%jri(ntyp),1)=plo(1:atoms%jri(ntyp),jspin)
156 0 : pqlo(1:atoms%jri(ntyp),2)=qlo(1:atoms%jri(ntyp),jspin)
157 0 : i = atoms%ulo_der(ilo,ntyp)
158 0 : IF(atoms%l_dulo(ilo,ntyp)) i=1
159 : CALL radsrdn(&
160 : e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
161 : usdus%ulos(ilo,ntyp,jspin),duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
162 0 : pqlo,usdus%dulos(ilo,ntyp,jspin),i)
163 0 : ddn1 = SQRT(ddn1)
164 0 : IF(atoms%l_dulo(ilo,ntyp)) ddn1=1.0
165 0 : plo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)/ddn1
166 0 : qlo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),2)/ddn1
167 0 : usdus%dulos(ilo,ntyp,jspin) = duds1/ddn1
168 0 : DEALLOCATE (glo,pqlo,filo)
169 : ENDIF
170 : !-apw+lo
171 : ENDDO
172 :
173 196 : IF (input%jspins.EQ.1) THEN
174 33860 : plo(1:atoms%jri(ntyp),2) = plo(1:atoms%jri(ntyp),1)
175 38 : e = (enpara%ello0(ilo,ntyp,1) + enpara%el0(l,ntyp,1) )/2
176 : ELSE
177 : e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilo,ntyp,input%jspins) +&
178 158 : enpara%el0(l,ntyp,1) + enpara%el0(l,ntyp,input%jspins) )/4
179 : END IF
180 335488 : CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
181 196 : IF (spav) THEN
182 24256 : DO i= 1,atoms%jmtd
183 24224 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
184 24256 : vso(i,2)= vso(i,1)
185 : ENDDO
186 : ENDIF
187 :
188 588 : DO i = 1, 2
189 1372 : DO j = 1, 2
190 784 : rsoc%rsoplop (ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),p (:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
191 784 : rsoc%rsoplopd(ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
192 784 : rsoc%rsopplo (ntyp,ilo,i,j) = radso(p (:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
193 1176 : rsoc%rsopdplo(ntyp,ilo,i,j) = radso(pd(:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
194 : ENDDO
195 : ENDDO
196 :
197 550 : DO i = 1,input%jspins
198 280796 : fint(:) = plo(:,i) * p(:,i) + qlo(:,i) * q(:,i)
199 354 : CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,i))
200 280796 : fint(:) = plo(:,i) * pd(:,i) + qlo(:,i) * qd(:,i)
201 550 : CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,i))
202 : ENDDO
203 :
204 548 : DO ilop = 1, atoms%nlo(ntyp)
205 548 : IF (atoms%llo(ilop,ntyp).EQ.l) THEN
206 :
207 550 : DO jspin = 1,input%jspins
208 281150 : vrTmp = vr(:,jspin)
209 48866 : if (atoms%l_nonpolbas(ntyp)) vrTmp = (vr(:,1)+vr(:,2))/2.0
210 354 : e = enpara%ello0(ilop,ntyp,jspin)
211 : CALL radsra(&
212 : e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
213 354 : ulops,dulops,nodeu,plop(:,jspin),q(:,1))
214 : !+apw+lo
215 550 : IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate orthogonal energy derivative at e
216 0 : ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
217 0 : glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
218 0 : pqlo(1:atoms%jri(ntyp),1)=plop(1:atoms%jri(ntyp),jspin)
219 0 : pqlo(1:atoms%jri(ntyp),2)=q(1:atoms%jri(ntyp),1)
220 0 : i = atoms%ulo_der(ilo,ntyp)
221 0 : IF(atoms%l_dulo(ilo,ntyp)) i=1
222 : CALL radsrdn(&
223 : e,l,vrtmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
224 : ulops,duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
225 0 : pqlo,dulops,i)
226 0 : plop(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)
227 0 : DEALLOCATE (glo,pqlo,filo)
228 : ENDIF
229 : !-apw+lo
230 : ENDDO
231 :
232 196 : IF (input%jspins.EQ.1) THEN
233 33860 : plop(1:atoms%jri(ntyp),2) = plop(1:atoms%jri(ntyp),1)
234 38 : e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilop,ntyp,1) )/2
235 : ELSE
236 : e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilo,ntyp,input%jspins) + &
237 158 : enpara%ello0(ilop,ntyp,1) + enpara%ello0(ilop,ntyp,input%jspins) )/4
238 : END IF
239 335488 : CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
240 196 : IF (spav) THEN
241 24256 : DO i= 1,atoms%jmtd
242 24224 : vso(i,1)= (vso(i,1)+vso(i,2))/2.
243 24256 : vso(i,2)= vso(i,1)
244 : ENDDO
245 : ENDIF
246 :
247 588 : DO i = 1, 2
248 1372 : DO j = 1, 2
249 : rsoc%rsoploplop(ntyp,ilo,ilop,i,j) =&
250 1176 : radso(plo(:atoms%jri(ntyp),i),plop(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
251 : ENDDO
252 : ENDDO
253 :
254 : ENDIF
255 : ENDDO
256 :
257 : ENDIF
258 : ENDDO ! end of lo-loop
259 :
260 : ENDDO ! end of l-loop
261 :
262 204 : DEALLOCATE ( p,pd,q,qd,plo,plop,qlo,fint,v0,vso )
263 : ! rsoplop (:,:,:,:) = 0.0
264 : ! rsoplopd(:,:,:,:) = 0.0
265 : ! rsopplo (:,:,:,:) = 0.0
266 : ! rsopdplo(:,:,:,:) = 0.0
267 : ! rsoploplop(:,:,:,:,:) = 0.0
268 :
269 204 : END SUBROUTINE sorad
270 : !--------------------------------------------------------------------
271 31696 : REAL FUNCTION radso(a,b,vso,dx,r0)
272 : !
273 : ! compute radial spin-orbit integrals
274 : !
275 : USE m_intgr, ONLY : intgr0
276 : IMPLICIT NONE
277 : !
278 : ! .. Scalar Arguments ..
279 : REAL, INTENT (IN) :: r0,dx
280 : ! ..
281 : ! .. Array Arguments ..
282 : REAL, INTENT (IN) :: a(:),b(:),vso(:)
283 : ! ..
284 : ! .. Local Arrays ..
285 31696 : REAL q(size(a))
286 : ! ..
287 23149056 : q = a*b*vso
288 31696 : CALL intgr0(q,r0,dx,size(a),radso)
289 :
290 31696 : RETURN
291 : END FUNCTION radso
292 : !--------------------------------------------------------------------
293 : END MODULE m_sorad
|