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_wann_rad_twf
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_rad_twf(
11 0 : > nwfs,jmtd,natd,ind,rwf,zona,regio,
12 0 : > us,dus,uds,duds,ff,gg,lmaxd,ikpt,
13 0 : > ntypd,ntp,jri,rmsh,dx,
14 0 : > nlod,flo,llo,nlo,
15 0 : < rads)
16 : c*********************************************************************
17 : c..this subroutine calculates on the radial grid the radial function
18 : c..corresponding to the twf, its mt sphere parameters, diffusivity
19 : c..and region
20 : c
21 : c..Y. Mokrousov
22 : c
23 : c..tidied-up && extension of rwf-parameter
24 : c..Frank Freimuth
25 : c*********************************************************************
26 :
27 : use m_intgr, only : intgr3
28 :
29 : implicit none
30 :
31 : integer, intent (in) :: nwfs,natd,jmtd,ntypd,lmaxd,ikpt,nlod
32 : integer, intent (in) :: rwf(nwfs),ind(nwfs),ntp(natd),jri(ntypd)
33 : real, intent (in) :: zona(nwfs),regio(nwfs),rmsh(jmtd,ntypd)
34 : real, intent (in) :: us(0:lmaxd,ntypd),dus(0:lmaxd,ntypd)
35 : real, intent (in) :: uds(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
36 : real, intent (in) :: ff(ntypd,jmtd,2,0:lmaxd)
37 : real, intent (in) :: gg(ntypd,jmtd,2,0:lmaxd),dx(ntypd)
38 : real, intent (in) :: flo(ntypd,jmtd,2,nlod)
39 : integer, intent (in) :: llo(nlod,ntypd)
40 : integer, intent (in) :: nlo(ntypd)
41 : real, intent (out) :: rads(nwfs,0:3,jmtd,2)
42 :
43 : integer :: nwf,nat,ntyp,j,n,l,lo
44 : real :: rr,alpha,const,fact,wronk,radi
45 0 : real :: acft,bcft,aa,bb,a1,b1,rho,radf(jmtd)
46 :
47 0 : do nwf = 1,nwfs
48 :
49 0 : if (abs(regio(nwf)-1.00000).ge.0.00001) then
50 : CALL juDFT_error("no projections outside muffin tins",calledby
51 0 : + ="wann_rad_twf")
52 : endif
53 :
54 0 : nat = ind(nwf)
55 0 : ntyp = ntp(nat)
56 0 : alpha = zona(nwf)
57 :
58 0 : rads (nwf,:,:,:) = 0.
59 :
60 0 : if (rwf(nwf).eq.-5) then
61 : c*********************************************************************
62 : c..For rwf.eq.-5 the following way of implementing is realized.
63 : c..Suppose, the twf has its center on some atom. The
64 : c..self-consistent potential V_l(r) and the energy parameters
65 : c..E_l are known. A linear combination of the u_l and \dot{u}_l
66 : c..is constructed, with the coefficients A and B such that
67 : c..on the mt boundary A*u_l + B*\dot{u}_l is smoothly continuous
68 : c..with the 1s function 2*((zona)**(1.5))*exp(-r*zona)
69 : c*********************************************************************
70 0 : aa = 2.*(zona(nwf)**(1.5))
71 0 : bb = zona(nwf)
72 :
73 0 : a1 = aa*exp(-bb*rmsh(jri(ntyp),ntyp))
74 0 : b1 = -aa*bb*exp(-bb*rmsh(jri(ntyp),ntyp))
75 :
76 0 : do l = 0,3
77 :
78 0 : wronk = us(l,ntyp)*duds(l,ntyp)-uds(l,ntyp)*dus(l,ntyp)
79 :
80 0 : acft = (a1*duds(l,ntyp) - b1*uds(l,ntyp))/wronk
81 0 : bcft = (b1*us(l,ntyp) - a1*dus(l,ntyp))/wronk
82 :
83 0 : do j = 1,jri(ntyp)
84 : rads(nwf,l,j,:) = acft*ff(ntyp,j,:,l) +
85 0 : + bcft*gg(ntyp,j,:,l)
86 : enddo
87 :
88 : enddo
89 0 : elseif(rwf(nwf).eq.0) then
90 : c***********************************************************
91 : c..For rwf.eq.0 the zona parameter mixes the radial function
92 : c..and its energy derivative linearly.
93 : c***********************************************************
94 0 : do l = 0,3
95 0 : do j = 1,jri(ntyp)
96 : rads(nwf,l,j,:) = ff(ntyp,j,:,l)*(1.-abs(zona(nwf))) +
97 0 : + zona(nwf)*gg(ntyp,j,:,l)
98 : enddo
99 : enddo
100 : c**************************************
101 : c..project onto local orbitals
102 : c**************************************
103 0 : elseif(rwf(nwf).eq.-6)then
104 0 : IF(nlod<1) CALL juDFT_error("nlod<1",calledby="wann_rad_twf"
105 0 : + )
106 0 : do l=0,3
107 0 : do j=1,jri(ntyp)
108 0 : rads(nwf,l,j,:)=flo(ntyp,j,:,1)
109 : enddo!j
110 : enddo!l
111 0 : elseif(rwf(nwf).eq.-7)then
112 0 : IF(nlod<2) CALL juDFT_error("nlod<2",calledby="wann_rad_twf"
113 0 : + )
114 0 : do l=0,3
115 0 : do j=1,jri(ntyp)
116 0 : rads(nwf,l,j,:)=flo(ntyp,j,:,2)
117 : enddo!j
118 : enddo!l
119 0 : elseif(rwf(nwf).eq.-8)then
120 0 : IF(nlod<3) CALL juDFT_error("nlod<3",calledby ="wann_rad_twf"
121 0 : + )
122 0 : do l=0,3
123 0 : do j=1,jri(ntyp)
124 0 : rads(nwf,l,j,:)=flo(ntyp,j,:,3)
125 : enddo!j
126 : enddo!l
127 0 : elseif(rwf(nwf).eq.-9)then
128 0 : rads(nwf,:,:,:)=0
129 0 : do lo=1,nlo(ntyp)
130 0 : l=llo(lo,ntyp)
131 0 : do j=1,jri(ntyp)
132 0 : rads(nwf,l,j,:)=flo(ntyp,j,:,lo)
133 : enddo!j
134 : enddo !lo
135 :
136 0 : elseif((-5.lt.rwf(nwf)).and.(rwf(nwf).lt.0)) then
137 : c************************************************************
138 : c..use the radial with l=abs(rwf)-1 for all components of the
139 : c..trial wave function
140 : c************************************************************
141 0 : do l = 0,3
142 0 : do j = 1,jri(ntyp)
143 : rads(nwf,l,j,:) =
144 : = ff(ntyp,j,:,abs(rwf(nwf))-1)*(1.-abs(zona(nwf))) +
145 0 : + zona(nwf)*gg(ntyp,j,:,abs(rwf(nwf))-1)
146 : enddo
147 : enddo
148 :
149 : c**************************************************************
150 : c..Hydrogen-like test orbitals.
151 : c..here we have tabulated the radial functions from orbitron.
152 : c..for 1 <= rwf <= 6 different components of the test orbital
153 : c..are constructed from different radials.
154 : c For Hybrides (e.g. sp3) it might be better to use the same
155 : c radials for all components. See below (7 <= rwf <= 12)
156 : c**************************************************************
157 0 : elseif (rwf(nwf).eq.1) then
158 : c..n=1
159 0 : do j = 1,jri(ntyp)
160 0 : rho = 2.*alpha*rmsh(j,ntyp)
161 : rads(nwf,0,j,1) = 2.*((alpha)**(1.5))*
162 0 : * exp(-rho/2.)
163 : enddo
164 :
165 0 : elseif (rwf(nwf).eq.2) then
166 : c..n=2
167 0 : do j = 1,jri(ntyp)
168 0 : rho = alpha*(rmsh(j,ntyp))
169 : rads(nwf,0,j,1) = (1./(2.*sqrt(2.)))*
170 : * (2.-rho)*
171 0 : * ((alpha)**(1.5))*exp(-rho/2.)
172 : rads(nwf,1,j,1) = (1./(2.*sqrt(6.)))*
173 : * rho*
174 0 : * ((alpha)**(1.5))*exp(-rho/2.)
175 : enddo
176 :
177 0 : elseif (rwf(nwf).eq.3) then
178 : c..n=3
179 0 : do j = 1,jri(ntyp)
180 0 : rho = 2*alpha*rmsh(j,ntyp)/3.
181 : rads(nwf,0,j,1) = (1./(9.*sqrt(3.)))*
182 : * (6. - 6.*rho + rho**2)*
183 0 : * ((alpha)**(1.5))*exp(-rho/2.)
184 : rads(nwf,1,j,1) = (1./(9.*sqrt(6.)))*
185 : * (4. - rho)*rho*
186 0 : * ((alpha)**(1.5))*exp(-rho/2.)
187 : rads(nwf,2,j,1) = (1./(9.*sqrt(30.)))*
188 : * rho*rho*
189 0 : * ((alpha)**(1.5))*exp(-rho/2.)
190 : enddo
191 :
192 0 : elseif (rwf(nwf).eq.4) then
193 : c..n=4
194 0 : do j = 1,jri(ntyp)
195 0 : rho = 2*alpha*rmsh(j,ntyp)/4.
196 : rads(nwf,0,j,1) = (1./96.)*
197 : * (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
198 0 : * ((alpha)**(1.5))*exp(-rho/2.)
199 : rads(nwf,1,j,1) = (1./(32.*sqrt(15.)))*
200 : * rho*(20. - 10.*rho + rho**2)*
201 0 : * ((alpha)**(1.5))*exp(-rho/2.)
202 : rads(nwf,2,j,1) = (1./(96.*sqrt(5.)))*
203 : * rho*rho*(6. - rho)*
204 0 : * ((alpha)**(1.5))*exp(-rho/2.)
205 : rads(nwf,3,j,1) = (1./(96.*sqrt(35.)))*
206 : * rho*rho*rho*
207 0 : * ((alpha)**(1.5))*exp(-rho/2.)
208 : enddo
209 :
210 0 : elseif (rwf(nwf).eq.5) then
211 0 : do j = 1,jri(ntyp)
212 0 : rho = 2*alpha*rmsh(j,ntyp)/5.
213 : rads(nwf,0,j,1) = (1./(300.*sqrt(5.)))*
214 : * (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
215 : + +(rho**4))*
216 0 : * ((alpha)**(1.5))*exp(-rho/2.)
217 : rads(nwf,1,j,1) = (1./(150.*sqrt(30.)))*
218 : * rho*(120. - 90.*rho + 18.*(rho**2) - (rho**3))*
219 0 : * ((alpha)**(1.5))*exp(-rho/2.)
220 : rads(nwf,2,j,1) = (1./(150.*sqrt(70.)))*
221 : * rho*rho*(42. - 14.*rho + (rho*rho))*
222 0 : * ((alpha)**(1.5))*exp(-rho/2.)
223 : rads(nwf,3,j,1) = (1./(300.*sqrt(70.)))*
224 : * rho*rho*rho*(8. - rho)*
225 0 : * ((alpha)**(1.5))*exp(-rho/2.)
226 : enddo
227 0 : elseif (rwf(nwf).eq.6) then
228 :
229 0 : do j = 1,jri(ntyp)
230 0 : rho = 2*alpha*rmsh(j,ntyp)/6.
231 : rads(nwf,0,j,1) = (1./2160.*(sqrt(6.)))*
232 : * (720. - 1800.*rho + 1200.*rho*rho-
233 : - 300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
234 0 : * ((alpha)**(1.5))*exp(-rho/2.)
235 : rads(nwf,1,j,1) = (1./432.*(sqrt(210.)))*
236 : * rho*(840. - 840.*rho + 252.*rho*rho-
237 : - 28.*rho*rho*rho + (rho**4))*
238 0 : * ((alpha)**(1.5))*exp(-rho/2.)
239 : enddo
240 :
241 0 : elseif (rwf(nwf).eq.7) then
242 : c..n=1
243 0 : do j = 1,jri(ntyp)
244 0 : rho = 2.*alpha*rmsh(j,ntyp)
245 0 : do l=0,3
246 : rads(nwf,l,j,1) = 2.*((alpha)**(1.5))*
247 0 : * exp(-rho/2.)
248 : enddo
249 : enddo
250 :
251 0 : elseif (rwf(nwf).eq.8) then
252 : c..n=2
253 0 : do j = 1,jri(ntyp)
254 0 : rho = alpha*(rmsh(j,ntyp))
255 0 : do l=0,3
256 : rads(nwf,l,j,1) = (1./(2.*sqrt(2.)))*
257 : * (2.-rho)*
258 0 : * ((alpha)**(1.5))*exp(-rho/2.)
259 : enddo
260 : enddo
261 :
262 0 : elseif (rwf(nwf).eq.9) then
263 : c..n=3
264 0 : do j = 1,jri(ntyp)
265 0 : rho = 2*alpha*rmsh(j,ntyp)/3.
266 0 : do l=0,3
267 : rads(nwf,l,j,1) = (1./(9.*sqrt(3.)))*
268 : * (6. - 6.*rho + rho**2)*
269 0 : * ((alpha)**(1.5))*exp(-rho/2.)
270 : enddo
271 : enddo
272 :
273 0 : elseif (rwf(nwf).eq.10) then
274 : c..n=4
275 0 : do j = 1,jri(ntyp)
276 0 : rho = 2*alpha*rmsh(j,ntyp)/4.
277 0 : do l=0,3
278 : rads(nwf,l,j,1) = (1./96.)*
279 : * (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
280 0 : * ((alpha)**(1.5))*exp(-rho/2.)
281 : enddo
282 : enddo
283 0 : elseif (rwf(nwf).eq.11) then
284 :
285 0 : do j = 1,jri(ntyp)
286 0 : rho = 2*alpha*rmsh(j,ntyp)/5.
287 0 : do l=0,3
288 : rads(nwf,l,j,1) = (1./(300.*sqrt(5.)))*
289 : * (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
290 : + +(rho**4))*
291 0 : * ((alpha)**(1.5))*exp(-rho/2.)
292 : enddo
293 : enddo
294 :
295 0 : elseif (rwf(nwf).eq.12) then
296 :
297 0 : do j = 1,jri(ntyp)
298 0 : rho = 2*alpha*rmsh(j,ntyp)/6.
299 0 : do l=0,3
300 : rads(nwf,0,j,1) = 0*(1./2160.*(sqrt(6.)))*
301 : * (720. - 1800.*rho + 1200.*rho*rho-
302 : - 300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
303 0 : * ((alpha)**(1.5))*exp(-rho/2.)
304 : enddo
305 : enddo
306 :
307 : else
308 :
309 : CALL juDFT_error("radial function is not tabulated" ,calledby
310 0 : + ="wann_rad_twf")
311 :
312 :
313 : endif
314 :
315 0 : if (ikpt.eq.1) then
316 0 : do l = 0,3
317 0 : do j = 1,jri(ntyp)
318 : radf(j) = rmsh(j,ntyp)*rmsh(j,ntyp)*rads(nwf,l,j,1)*
319 0 : * rads(nwf,l,j,1)
320 : c radf(j) = rads(nwf,l,j,1)*rads(nwf,l,j,1)
321 : enddo
322 0 : call intgr3(radf,rmsh(1,ntyp),dx(ntyp),jri(ntyp),radi)
323 0 : write (6,*)
324 0 : write (6,*) 'Wannier Function N:',nwf
325 0 : write (6,*) 'angular momentum',l
326 0 : write (6,*) 'radial function at the MT boundary:',
327 0 : & rads(nwf,l,jri(ntyp),1)
328 0 : write (6,*) 'norma =',radi
329 : enddo
330 : endif
331 :
332 : enddo ! by twfs
333 :
334 0 : return
335 :
336 : end subroutine wann_rad_twf
337 : end module m_wann_rad_twf
|
