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_radsrdn
8 : use m_juDFT
9 : CONTAINS
10 38 : SUBROUTINE radsrdn(
11 38 : > e,l,vr,r0,h,jri,c,
12 38 : < udn,dudn,ddnn,nodedn,fn,fni,
13 38 : > f0,du0,n)
14 : C*********************************************************************
15 : C Calculates the nth energy derivative of the scalar relativistic
16 : C wavefuction for energy e and angular momentum l.
17 : C
18 : C The large and small components
19 : C of the nth derivative of f0 are returned in fn,
20 : C of the inhomogeneous part of the nth derivative of the Dirac
21 : C equation ( approx. n * (n-1)st derivative of f0 ) are returned
22 : C in fni.
23 : C For non-relativistic case, choose large c ( >> 137.0359895).
24 : C
25 : C f0 (= primitive, 0th derivative) and
26 : C the radial derivative du0 of f0,
27 : C are required.
28 : C
29 : C The nth derivative of the Dirac equation is solved in radsrdn1.
30 : C
31 : C C. Friedrich Apr. 2005
32 : C*********************************************************************
33 : USE m_intgr, ONLY : intgr0
34 : IMPLICIT NONE
35 : C ..
36 : C .. Scalar Arguments ..
37 : INTEGER, INTENT (IN) :: l,n
38 : INTEGER, INTENT (IN) :: jri
39 : INTEGER, INTENT (OUT):: nodedn
40 : REAL, INTENT (IN) :: c
41 : REAL, INTENT (IN) :: du0,e,h,r0
42 : REAL, INTENT (OUT):: ddnn,dudn,udn
43 : C ..
44 : C .. Array Arguments ..
45 : REAL, INTENT (IN) :: f0(:,:),vr(:)
46 : REAL, INTENT (OUT):: fn(:,:),fni(:,:)
47 : C ..
48 : C .. Local Scalars ..
49 : INTEGER i,j,k
50 : REAL scprod,ovlap(9,9),fac
51 : C ..
52 : C .. Local Arrays ..
53 38 : REAL f(size(vr),2,0:n),help(size(vr))
54 : C
55 : C
56 38 : IF(n<1) CALL juDFT_error("n<1",calledby ="radsrdn")
57 57646 : f(:,:,0) = f0(:,:)
58 : C
59 114 : DO i=1,n
60 : C Calculate
61 : C scprod = < f0(:,:) | f(:,:,i) >
62 : C = - 1/2 SUM ( i over j ) < f(:,:,j) | f(:,:,i-j) >
63 : C j=1,i ----fac----- -----------ovlap---------
64 : C (This follows from d^n/de^n <f0|f0> = 0.)
65 76 : scprod=0
66 114 : DO j=1,int(i/2)
67 28804 : help(:) = f(:,1,j)*f(:,1,i-j) + f(:,2,j)*f(:,2,i-j)
68 38 : CALL intgr0(help,r0,h,jri,ovlap(j,i-j))
69 38 : ovlap(i-j,j)=ovlap(j,i-j)
70 38 : fac=1 !
71 76 : DO k=i-j+1,i !
72 76 : fac=fac*k ! Calculate
73 : ENDDO ! fac = ( i over j ) = i! / j! / (i-j)!
74 38 : DO k=2,j !
75 38 : fac=fac/k !
76 : ENDDO !
77 38 : IF(j.eq.i-j) fac=fac/2
78 114 : scprod=scprod-fac*ovlap(j,i-j)
79 : ENDDO
80 : C
81 : C scprod=0 ! Instead of the nth derivative a linear combination
82 : ! of derivatives of order <=n is calculated.
83 : ! It can be shown, that the udot-coefficient in the LO
84 : ! construction in setabc1lo is always zero with this choice.
85 : ! This reduces the error arising from the inconsistency
86 : ! between an "exact" skalar-relativistic udot and a
87 : ! non-relativistic Hamiltonian. Uncomment this, if you
88 : ! have set RELATIVISTIC_CORRECTIONS in radsrd. But note,
89 : ! that then relativistic corrections are not fully treated
90 : ! in the higher derivatives calculated here.
91 76 : CALL getfni(e,l,vr,r0,h,jri,c,i,f,fni)
92 : CALL radsrdn1(
93 : > e,l,vr,r0,h,jri,c,
94 : < udn,dudn,ddnn,nodedn,f(:,1,i),f(:,2,i),
95 114 : > fni(:,1),fni(:,2),f0(:,1),f0(:,2),du0,scprod)
96 : ENDDO
97 57646 : fn(:,:) = f(:,:,n)
98 28804 : help(:) = fni(:,1)
99 28804 : fni(:,1)=-fni(:,2)
100 28804 : fni(:,2)= help(:)*c
101 38 : END SUBROUTINE radsrdn
102 :
103 : C---------------------------------------------------------------------
104 :
105 152 : SUBROUTINE radsrdn1(
106 152 : > e,l,vr,r0,h,jri,c,
107 76 : < ud,dud,ddn,nodes,pe,qe,
108 152 : > pi,qi,phom,qhom,dus,scprod)
109 : C*********************************************************************
110 : C Solves the "inhomogeneous" Dirac equation for energy e and angular
111 : C momentum l with pi (qi) on one side instead of zero.
112 : C
113 : C The large (small) component is returned in pe (qe).
114 : C For non-relativistic case, choose large c ( >> 137.0359895).
115 : C
116 : C The functions
117 : C pi,qi,
118 : C phom,qhom (= homogeneous solution, 0th derivative,
119 : C as calculated in radsra),
120 : C as well as
121 : C the radial derivative dus of phom,
122 : C and scprod = < p(q)e | p(q)hom >
123 : C are required.
124 : C
125 : C Modified from radsrd. C. Friedrich Apr. 2005
126 : C*********************************************************************
127 : C
128 : USE m_intgr, ONLY : intgr0
129 : IMPLICIT NONE
130 : C ..
131 : C .. Scalar Arguments ..
132 : INTEGER, INTENT (IN) :: l
133 : INTEGER, INTENT (IN) :: jri
134 : INTEGER, INTENT (OUT):: nodes
135 : REAL, INTENT (IN) :: c
136 : REAL, INTENT (IN) :: dus,e,h,r0,scprod
137 : REAL, INTENT (OUT):: ddn,dud,ud
138 : C ..
139 : C .. Array Arguments ..
140 : REAL, INTENT (IN) :: pi(:),qi(:),vr(:)
141 : REAL, INTENT (OUT):: pe(:),qe(:)
142 : REAL, INTENT (IN) :: phom(:),qhom(:)
143 : C ..
144 : C .. Local Scalars ..
145 : REAL dr,drh,eps,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,
146 : + sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,t1,t2,yn,zn,cin,cin2
147 : INTEGER i,it
148 : C ..
149 : C .. Local Arrays ..
150 76 : REAL pp(size(pi)),qp(size(pi))
151 : C ..
152 : C .. Intrinsic Functions ..
153 : INTRINSIC abs,exp,sign,sqrt
154 : C ..
155 : C .. Data statements ..
156 : DATA eps/1.e-06/
157 : C ..
158 76 : cin = 1.0/c
159 76 : cin2 = cin*cin
160 : c
161 76 : IF (jri>size(pi)) CALL juDFT_error("dimension too small",
162 0 : + calledby ="radsrdn")
163 : c---> set up initial conditions
164 76 : fl1 = l* (l+1)
165 76 : pe(1) = 0
166 76 : qe(1) = 0
167 76 : pp(1) = r0*pi(1)
168 76 : qp(1) = r0*qi(1)
169 : c---> use 4th order runge-kutta to get first few mesh points
170 76 : dr = exp(h)
171 76 : drh = sqrt(dr)
172 76 : r = r0
173 456 : DO i = 1,5
174 380 : rh2 = drh*r
175 380 : rh = dr*r
176 380 : sk1 = h*pp(i)
177 380 : sl1 = h*qp(i)
178 380 : rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
179 380 : rm = 2.*rh2 - cin2*rve
180 380 : yn = pe(i) + 0.5*sk1
181 380 : zn = qe(i) + 0.5*sl1
182 380 : t1 = 0.5*rh2*(pi(i)+pi(i+1))
183 380 : t2 = 0.5*rh2*(qi(i)+qi(i+1))
184 380 : sk2 = h* (rm*zn + yn + t1)
185 380 : sl2 = h* (-zn + (fl1/rm+rve)*yn + t2)
186 380 : yn = pe(i) + 0.5*sk2
187 380 : zn = qe(i) + 0.5*sl2
188 380 : sk3 = h* (rm*zn + yn + t1)
189 380 : sl3 = h* (-zn + (fl1/rm+rve)*yn + t2)
190 380 : rve = vr(i+1) - rh*e
191 380 : rm = 2.*rh - cin2*rve
192 380 : yn = pe(i) + sk3
193 380 : zn = qe(i) + sl3
194 380 : t1 = rh*pi(i+1)
195 380 : t2 = rh*qi(i+1)
196 380 : sk4 = h* (rm*zn + yn + t1 )
197 380 : sl4 = h* (-zn + (fl1/rm+rve)*yn + t2 )
198 380 : pe(i+1) = pe(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
199 380 : qe(i+1) = qe(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
200 380 : pp(i+1) = rm*qe(i+1) + pe(i+1) + t1
201 380 : qp(i+1) = -qe(i+1) + (fl1/rm+rve)*pe(i+1) + t2
202 456 : r = rh
203 : ENDDO
204 76 : nodes = 0
205 : c---> adams-bashforth-moulton predictor-corrector
206 57152 : predictor: DO i = 6,jri - 1
207 57076 : r = r*dr
208 : c---> predictor
209 : p0 = pe(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
210 57076 : + 7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
211 : q0 = qe(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
212 57076 : + 7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
213 : c---> evaluate derivatives at next point
214 57076 : rve = vr(i+1) - r*e
215 57076 : rm = 2.*r - cin2*rve
216 57076 : t1 = r*pi(i+1)
217 57076 : t2 = r*qi(i+1)
218 57076 : p1p = rm*q0 + p0 + t1
219 57076 : q1p = -q0 + (fl1/rm+rve)*p0 + t2
220 : c---> corrector
221 61940 : corrector: DO it = 1,5
222 : p1 = pe(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
223 61940 : + 482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
224 : q1 = qe(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
225 61940 : + 482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
226 : c---> final evaluation
227 61940 : p1p = rm*q1 + p1 + t1
228 61940 : q1p = -q1 + (fl1/rm+rve)*p1 + t2
229 : c---> test quality of corrector and iterate if necessary
230 61940 : erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
231 61940 : erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
232 61940 : IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
233 4864 : p0 = p1
234 61940 : q0 = q1
235 : ENDDO corrector
236 57076 : IF (it > 5) CALL juDFT_error("Not converged.",calledby
237 0 : + ="radsrdn")
238 : c---> store values
239 57076 : pe(i+1) = p1
240 57076 : qe(i+1) = q1
241 57076 : pp(i+1) = p1p
242 57076 : qp(i+1) = q1p
243 57152 : nodes = nodes + 0.501*abs(sign(1.0,pe(i+1))-sign(1.0,pe(i)))
244 : ENDDO predictor
245 : c---> ensure < p(q)hom | p(q)e > = scprod
246 57608 : DO i = 1,jri
247 57608 : qe(i) = cin*qe(i)
248 : ENDDO
249 57608 : DO i = 1,jri
250 57608 : qp(i) = phom(i)*pe(i) + qhom(i)*qe(i)
251 : ENDDO
252 76 : CALL intgr0(qp,r0,h,jri,t)
253 76 : dud = (pp(jri)-pe(jri))/ (r*r)
254 57608 : DO i = 1,jri
255 57532 : pe(i) = pe(i) + (scprod-t)*phom(i)
256 57608 : qe(i) = qe(i) + (scprod-t)*qhom(i)
257 : ENDDO
258 76 : ud = pe(jri)/r
259 76 : dud = dud + (scprod-t)*dus
260 57608 : DO i = 1,jri
261 57608 : qp(i) = pe(i)*pe(i) + qe(i)*qe(i)
262 : ENDDO
263 76 : CALL intgr0(qp,r0,h,jri,ddn)
264 76 : RETURN
265 : END SUBROUTINE radsrdn1
266 :
267 : C---------------------------------------------------------------------
268 :
269 76 : SUBROUTINE getfni(e,l,vr,r0,h,jri,c,n,f,fni)
270 : C*********************************************************************
271 : C Calculates the inhomogeneous part of the nth derivative of the
272 : C Dirac equation with energy e and angular momentum l.
273 : C
274 : C It is
275 : C
276 : C fni(:,1) = n*f(:,2,n-1)/c (Note that f(:,2)=Q/c)
277 : C
278 : C l(l+1)
279 : C fni(:,2) =-n*f(:,1,n-1)* ( 1 + -------- )
280 : C (2Mrc)^2
281 : C
282 : C n i n! l(l+1)
283 : C + SUM (-1) -----------------------------
284 : C i=2 (n-i)! (2*M)^(n+1) r^2 c^(2n)
285 : C
286 : C with M = 1 + (E-V)/(2*c^2).
287 : C
288 : C C. Friedrich Apr. 2005
289 : C*********************************************************************
290 : C
291 :
292 : IMPLICIT NONE
293 : C ..
294 : C .. Scalar Arguments ..
295 : INTEGER, INTENT (IN) :: n,jri,l
296 : REAL, INTENT (IN) :: c,e,h,r0
297 : C ..
298 : C .. Array Arguments ..
299 : REAL, INTENT (IN) :: f(:,:,0:),vr(:)
300 : REAL, INTENT (OUT):: fni(:,:)
301 : C ..
302 : C .. Local Scalars ..
303 : INTEGER i
304 : REAL r
305 : C ..
306 : C .. Local Arrays ..
307 76 : REAL hlp(jri),fac(0:n),rmsh(jri)
308 : C
309 76 : IF(n.lt.1) CALL juDFT_error("getfni: n<1!",calledby="radsrdn")
310 76 : r=exp(h)
311 76 : rmsh(1)=r0
312 57532 : DO i=2,jri
313 57532 : rmsh(i)=rmsh(i-1)*r
314 : ENDDO
315 76 : fac(0)=1
316 190 : DO i=1,n
317 190 : fac(i)=fac(i-1)*i
318 : ENDDO
319 57608 : hlp(:) = 2*rmsh(:) + (e*rmsh(:)-vr(:jri))/c**2
320 57608 : fni(:jri,1) = n*f(:jri,2,n-1) /c
321 : fni(:jri,2) =-n*f(:jri,1,n-1)
322 57608 : & * ( 1 + l*(l+1)* (hlp(:)*c)**(-2) )
323 57608 : hlp(:) = hlp(:)/rmsh(:)
324 114 : DO i=2,n
325 38 : r = fac(n)/fac(n-i) * (-1)**i * l*(l+1) / c**(2*i)
326 : fni(:jri,2) = fni(:jri,2) +
327 28880 : & f(:jri,1,n-i) * r / hlp(:)**(i+1) / rmsh(:)**2
328 : ENDDO
329 76 : END SUBROUTINE getfni
330 :
331 : END MODULE m_radsrdn
332 :
|
