Line data Source code
1 : MODULE m_radsrd
2 : use m_juDFT
3 :
4 : c*********************************************************************
5 : c calculates the energy derivative of the scalar relativistic
6 : c wavefuction for energy e and angular momentum l.
7 : c the large (small) component is returned in pe (qe). for non-
8 : c relativistic case, set cin=0.
9 : c the functions p,q and the radial derivative dus are required.
10 : c based on code by m. weinert
11 : c*********************************************************************
12 :
13 : c*********************************************************************
14 : c
15 : c Uncomment the following, if you want to calculate an "exact" udot
16 : c including relativistic corrections in the inhomogeneous part of the
17 : c Dirac-DGL. Note that these corrections are inconsistent with the
18 : c construction of the Hamiltonian in hssphn, where the inhomogeneous
19 : c part is assumed to be non-relativistic.
20 : c
21 : c C. Friedrich June 2009
22 : c*********************************************************************
23 : # define RELATIVISTIC_CORRECTIONS
24 :
25 : # ifdef RELATIVISTIC_CORRECTIONS
26 : # define CORR *(1.+cin2*fl1/rm**2)
27 : # else
28 : # define CORR
29 : # endif
30 :
31 : CONTAINS
32 160958 : SUBROUTINE radsrd(
33 160958 : > e,l,vr,r0,h,jri,c,
34 80479 : < ud,dud,ddn,nodes,pe,qe,
35 160958 : > p,q,dus)
36 :
37 : USE m_intgr, ONLY : intgr0
38 : IMPLICIT NONE
39 : C ..
40 : C .. Scalar Arguments ..
41 : INTEGER, INTENT (IN) :: l
42 : INTEGER, INTENT (IN) :: jri
43 : INTEGER, INTENT (OUT):: nodes
44 : REAL, INTENT (IN) :: c
45 : REAL, INTENT (IN) :: dus,e,h,r0
46 : REAL, INTENT (OUT):: ddn,dud,ud
47 : C ..
48 : C .. Array Arguments ..
49 : REAL, INTENT (IN) :: p(:),q(:),vr(:)
50 : REAL, INTENT (OUT):: pe(:),qe(:)
51 : C ..
52 : C .. Local Scalars ..
53 : REAL dr,drh,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,
54 : + sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,t1,t2,yn,zn,cin,cin2
55 : INTEGER i,it
56 : C ..
57 : C .. Local Arrays ..
58 80479 : REAL pp(size(p)),qp(size(p))
59 : C ..
60 : C .. Data statements ..
61 : REAL,PARAMETER :: eps=1.e-06
62 : C ..
63 80479 : cin = 1.0/c
64 80479 : cin2 = cin*cin
65 : c
66 80479 : IF(jri>size(p)) CALL juDFT_error("dimension too small",
67 0 : + calledby="radsrd")
68 : c---> set up initial conditions
69 80479 : fl1 = l* (l+1)
70 80479 : rm = 2.*r0 + cin2* (r0*e-vr(1))
71 80479 : pe(1) = 0
72 80479 : qe(1) = 0
73 80479 : pp(1) = r0*cin*q(1)
74 80479 : qp(1) = - r0*p(1) CORR
75 : c---> use 4th order runge-kutta to get first few mesh points
76 80479 : dr = exp(h)
77 80479 : drh = sqrt(dr)
78 80479 : r = r0
79 482874 : DO i = 1,5
80 402395 : rh2 = drh*r
81 402395 : rh = dr*r
82 402395 : sk1 = h*pp(i)
83 402395 : sl1 = h*qp(i)
84 402395 : rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
85 402395 : rm = 2.*rh2 - cin2*rve
86 402395 : yn = pe(i) + 0.5*sk1
87 402395 : zn = qe(i) + 0.5*sl1
88 402395 : t1 = 0.5*rh2*cin* (q(i)+q(i+1))
89 402395 : t2 = 0.5*rh2* (p(i)+p(i+1)) CORR
90 402395 : sk2 = h* (rm*zn+yn+t1)
91 402395 : sl2 = h* ((fl1/rm+rve)*yn-zn-t2)
92 402395 : yn = pe(i) + 0.5*sk2
93 402395 : zn = qe(i) + 0.5*sl2
94 402395 : sk3 = h* (rm*zn+yn+t1)
95 402395 : sl3 = h* ((fl1/rm+rve)*yn-zn-t2)
96 402395 : rve = vr(i+1) - rh*e
97 402395 : rm = 2.*rh - cin2*rve
98 402395 : yn = pe(i) + sk3
99 402395 : zn = qe(i) + sl3
100 402395 : sk4 = h* (rm*zn+yn+rh*cin*q(i+1))
101 402395 : sl4 = h* ((fl1/rm+rve)*yn-zn-rh*p(i+1) CORR )
102 402395 : pe(i+1) = pe(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
103 402395 : qe(i+1) = qe(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
104 402395 : pp(i+1) = rm*qe(i+1) + pe(i+1) + rh*cin*q(i+1)
105 : qp(i+1) = (fl1/rm+rve)*pe(i+1) - qe(i+1) -
106 402395 : + rh*p(i+1) CORR
107 482874 : r = rh
108 : ENDDO
109 80479 : nodes = 0
110 :
111 : c---> adams-bashforth-moulton predictor-corrector
112 58290974 : predictor: DO i = 6,jri - 1
113 58210495 : r = r*dr
114 : c---> predictor
115 : p0 = pe(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
116 58210495 : + 7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
117 : q0 = qe(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
118 58210495 : + 7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
119 : c---> evaluate derivatives at next point
120 58210495 : rve = vr(i+1) - r*e
121 58210495 : rm = 2.*r - cin2*rve
122 58210495 : t1 = cin*r*q(i+1)
123 58210495 : t2 = r*p(i+1) CORR
124 58210495 : p1p = rm*q0 + p0 + t1
125 58210495 : q1p = (fl1/rm+rve)*p0 - q0 - t2
126 : c---> corrector
127 63732141 : corrector: DO it = 1,5
128 : p1 = pe(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
129 63732141 : + 482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
130 : q1 = qe(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
131 63732141 : + 482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
132 : c---> final evaluation
133 63732141 : p1p = rm*q1 + p1 + t1
134 63732141 : q1p = (fl1/rm+rve)*p1 - q1 - t2
135 : c---> test quality of corrector and iterate if necessary
136 63732141 : erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
137 63732141 : erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
138 63732141 : IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
139 5521646 : p0 = p1
140 63732141 : q0 = q1
141 : ENDDO corrector
142 : c---> store values
143 58210495 : pe(i+1) = p1
144 58210495 : qe(i+1) = q1
145 58210495 : pp(i+1) = p1p
146 58210495 : qp(i+1) = q1p
147 58290974 : nodes = nodes + 0.501*abs(sign(1.0,pe(i+1))-sign(1.0,pe(i)))
148 : ENDDO predictor
149 : c---> ensure orthogonality
150 58773848 : DO i = 1,jri
151 58773848 : qe(i) = cin*qe(i)
152 : ENDDO
153 58773848 : DO i = 1,jri
154 58773848 : qp(i) = p(i)*pe(i) + q(i)*qe(i)
155 : ENDDO
156 80479 : CALL intgr0(qp,r0,h,jri,t)
157 80479 : dud = (pp(jri)-pe(jri))/ (r*r)
158 58773848 : DO i = 1,jri
159 58693369 : pe(i) = pe(i) - t*p(i)
160 58773848 : qe(i) = qe(i) - t*q(i)
161 : ENDDO
162 80479 : ud = pe(jri)/r
163 80479 : dud = dud - t*dus
164 58773848 : DO i = 1,jri
165 58773848 : qp(i) = pe(i)*pe(i) + qe(i)*qe(i)
166 : ENDDO
167 80479 : CALL intgr0(qp,r0,h,jri,ddn)
168 :
169 80479 : END SUBROUTINE radsrd
170 : END MODULE m_radsrd
|