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