Line data Source code
1 : MODULE m_radsra
2 : use m_juDFT
3 : c*********************************************************************
4 : c calculates the scalar relativistic wavefuction for energy e and
5 : c angular momentum l for the potential vr by integrating outward.
6 : c the large (small) component is returned in p (q). for non-
7 : c relativistic case, set cin=0. based on code by m. weinert
8 : c*********************************************************************
9 : CONTAINS
10 752592 : SUBROUTINE radsra(
11 376296 : > e,l,vr,r0,h,jri,c,
12 376296 : < us,dus,nodes,p,q)
13 :
14 : USE m_intgr, ONLY : intgr0
15 : IMPLICIT NONE
16 : C ..
17 : C .. Scalar Arguments ..
18 : INTEGER, INTENT (IN) :: l
19 : INTEGER, INTENT (IN) :: jri
20 : INTEGER, INTENT (OUT):: nodes
21 : REAL, INTENT (IN) :: e,h,r0
22 : REAL, INTENT (IN) :: c
23 : REAL, INTENT (OUT):: dus,us
24 : C ..
25 : C .. Array Arguments ..
26 : REAL, INTENT (IN) :: vr(:)
27 : REAL, INTENT (OUT):: p(:),q(:)
28 : C ..
29 : C .. Local Scalars ..
30 : REAL dr,drh,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,s,
31 : + sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,yn,zn,cin,cin2
32 : INTEGER i,it
33 : C ..
34 : C .. Local Arrays ..
35 376296 : REAL pp(jri),qp(jri)
36 : C ..
37 : C .. Data statements ..
38 : REAL,PARAMETER :: eps=1.e-06
39 : C ..
40 376296 : cin = 1.0/c
41 376296 : cin2 = cin*cin
42 :
43 376296 : IF (jri>SIZE(vr).OR.jri>SIZE(p)) CALL juDFT_error
44 0 : + ("BUG: data dimension in radsra too small",calledby ="radsra")
45 :
46 :
47 : c---> set up initial conditions
48 376296 : fl1 = l* (l+1)
49 376296 : s = sqrt(fl1+1-cin2*vr(1)*vr(1)) - 1.
50 376296 : rm = 2.*r0 + cin2* (r0*e-vr(1))
51 376296 : p(1) = r0** (l+1)
52 376296 : q(1) = s*p(1)/rm
53 376296 : pp(1) = rm*q(1) + p(1)
54 376296 : qp(1) = (fl1/rm+vr(1)-r0*e)*p(1) - q(1)
55 : c---> use 4th order runge-kutta to get first few mesh points
56 376296 : dr = exp(h)
57 376296 : drh = sqrt(dr)
58 376296 : r = r0
59 2257776 : DO i = 1,5
60 1881480 : rh2 = drh*r
61 1881480 : rh = dr*r
62 1881480 : sk1 = h*pp(i)
63 1881480 : sl1 = h*qp(i)
64 1881480 : rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
65 1881480 : rm = 2.*rh2 - cin2*rve
66 1881480 : yn = p(i) + 0.5*sk1
67 1881480 : zn = q(i) + 0.5*sl1
68 1881480 : sk2 = h* (rm*zn+yn)
69 1881480 : sl2 = h* ((fl1/rm+rve)*yn-zn)
70 1881480 : yn = p(i) + 0.5*sk2
71 1881480 : zn = q(i) + 0.5*sl2
72 1881480 : sk3 = h* (rm*zn+yn)
73 1881480 : sl3 = h* ((fl1/rm+rve)*yn-zn)
74 1881480 : rve = vr(i+1) - rh*e
75 1881480 : rm = 2.*rh - cin2*rve
76 1881480 : yn = p(i) + sk3
77 1881480 : zn = q(i) + sl3
78 1881480 : sk4 = h* (rm*zn+yn)
79 1881480 : sl4 = h* ((fl1/rm+rve)*yn-zn)
80 1881480 : p(i+1) = p(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
81 1881480 : q(i+1) = q(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
82 1881480 : pp(i+1) = rm*q(i+1) + p(i+1)
83 1881480 : qp(i+1) = (fl1/rm+rve)*p(i+1) - q(i+1)
84 2257776 : r = rh
85 : ENDDO
86 376296 : nodes = 0
87 : c---> adams-bashforth-moulton predictor-corrector
88 278551298 : predictor: DO i = 6,jri - 1
89 278175002 : r = r*dr
90 : c---> predictor
91 : p0 = p(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
92 278175002 : + 7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
93 : q0 = q(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
94 278175002 : + 7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
95 : c---> evaluate derivatives at next point
96 278175002 : rve = vr(i+1) - r*e
97 278175002 : rm = 2.*r - cin2*rve
98 278175002 : p1p = rm*q0 + p0
99 278175002 : q1p = (fl1/rm+rve)*p0 - q0
100 : c---> corrector
101 321265190 : corrector: DO it = 1,5
102 : p1 = p(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
103 318327024 : + 482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
104 : q1 = q(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
105 318327024 : + 482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
106 : c---> final evaluation
107 318327024 : p1p = rm*q1 + p1
108 318327024 : q1p = (fl1/rm+rve)*p1 - q1
109 : c---> test quality of corrector and iterate if necessary
110 318327024 : erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
111 318327024 : erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
112 318327024 : IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
113 43090188 : p0 = p1
114 321265190 : q0 = q1
115 : ENDDO corrector
116 : c---> store values
117 278175002 : p(i+1) = p1
118 278175002 : q(i+1) = q1
119 278175002 : pp(i+1) = p1p
120 278175002 : qp(i+1) = q1p
121 278551298 : nodes = nodes + 0.501*abs(sign(1.0,p(i+1))-sign(1.0,p(i)))
122 : ENDDO predictor
123 :
124 : c---> normalize function
125 280809074 : DO i = 1,jri
126 280809074 : q(i) = cin*q(i)
127 : ENDDO
128 280809074 : DO i = 1,jri
129 280809074 : qp(i) = p(i)*p(i) + q(i)*q(i)
130 : ENDDO
131 376296 : CALL intgr0(qp,r0,h,jri,t)
132 376296 : t = 1.0/sqrt(t)
133 376296 : us = t*p(jri)/r
134 376296 : dus = t* (pp(jri)-p(jri))/ (r*r)
135 280809074 : DO i = 1,jri
136 280432778 : p(i) = t*p(i)
137 280809074 : q(i) = t*q(i)
138 : ENDDO
139 :
140 376296 : END SUBROUTINE radsra
141 : END MODULE m_radsra
|