Line data Source code
1 : MODULE m_inwint
2 : CONTAINS
3 36694577 : SUBROUTINE inwint(
4 410007 : > e,fl,ki,fkap,cs,cis,s,z,h,dd,rn,rnot,msh,vr,
5 410007 : X a,b,
6 : < kj)
7 : c
8 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
9 : c-----x perform the inward integration. this routine is a x----
10 : c-----x derivative of start2/diff2. much unnecessary indexing x----
11 : c-----x has been removed among other things. dale koelling x----
12 : c
13 : c adams' procedure. more stable for GGA. Feb.20,1998. T.A.
14 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
15 : use m_juDFT
16 : IMPLICIT NONE
17 :
18 : c .. Arguments ..
19 : INTEGER,INTENT (IN) :: msh,ki
20 : INTEGER,INTENT (OUT) :: kj
21 : REAL, INTENT (IN) :: e,fl,fkap,cs,cis,s,z,h,dd,rn,rnot
22 : REAL, INTENT (IN) :: vr(msh)
23 : c ..
24 : REAL, INTENT (INOUT) :: a(msh),b(msh)
25 : c ..
26 : c .. Local Scalars ..
27 : REAL atk,btk,d,da1,da2,da3,da4,db1,db2,db3,db4,eabsv,eps,f1,f2,f3,
28 : + f4,f5,fg1,fg2,fg3,g,g1,g2,g3,g4,g5,h3,q11,q22,r,rp12,rp21
29 : INTEGER I,k,l,n
30 : REAL :: da(5),db(5)
31 : c ..
32 : c .. intrinsic functions ..
33 : INTRINSIC abs,mod,sqrt
34 : c ..
35 : c .. equivalences ..
36 : cta+
37 : EQUIVALENCE (da1,da(1)), (da2,da(2)), (da3,da(3)), (da4,da(4)),
38 : & (da5,da(5))
39 : EQUIVALENCE (db1,db(1)), (db2,db(2)), (db3,db(3)), (db4,db(4)),
40 : & (db5,db(5))
41 :
42 : REAL c1,c5,c9,c19, t9,t37,t55,t59, da5,db5
43 : cta-
44 : c ..
45 : c .. data statements ..
46 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
47 : c-----x this version modified to accept positive energies. x----
48 : c-----x d.d.koelling 7/7/77 x----x
49 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
50 : c---* eps value used to determine the practical infinity
51 : DATA f1/1.045833333e0/,f2/2.691666666e0/,f3/1.1e0/
52 : DATA f4/.4416666666e0/,f5/.07916666666e0/
53 : DATA g1/.9666666666e0/,g2/4.133333333e0/,g3/.8e0/
54 : DATA g4/.1333333333e0/,g5/.03333333333e0/
55 : DATA fg1/.9333333333e0/,fg2/4.266666666e0/,fg3/1.6e0/
56 : DATA eps/75.0e0/
57 : c ..
58 410007 : h3 = -h/3.0e0
59 410007 : n = msh
60 : cta+
61 410007 : t55=55.e0/8.e0
62 410007 : t59=59.e0/8.e0
63 410007 : t37=37.e0/8.e0
64 410007 : t9=9.e0/8.e0
65 :
66 410007 : c19=19.e0/8.e0
67 410007 : c9=9.e0/8.e0
68 410007 : c5=5.e0/8.e0
69 410007 : c1=1.e0/8.e0
70 : cta-
71 :
72 : c e: energy
73 : c ki: the grid point corresponding to classical terning point
74 : cc at which potential exceeded the energy.
75 : cc ki=249 for Fe and l=0.
76 : c kj: undefined at the beginning. finally the farthest grid
77 : cc point which can be considered as infinity for the
78 : cc level to be obtained.
79 :
80 410007 : d = 1.0e0/dd
81 410007 : g = sqrt(fkap**2- (cis*z)**2)
82 410007 : q11 = s*fkap - g
83 410007 : q22 = -g - s*fkap
84 410007 : kj = ki + 10
85 30954479 : r = rnot* (dd**kj)
86 :
87 30954479 : 10 kj = kj + 1
88 30954479 : IF (kj>msh) CALL juDFT_error("inward integration failed",calledby
89 : + ="inwint",hint
90 0 : + ="This often indicates that the potential is unphysical")
91 :
92 : c eps is a big energy. thus if the next statement is fulfilled,
93 : cc the grid is considered to be infinitely far.
94 :
95 30954479 : IF (r* (vr(kj)-e*r).gt.eps) GOTO 20
96 30544472 : r = r*dd
97 30544472 : IF (kj.lt.n) go to 10
98 0 : b(kj) = mod(fl,2.0e0)
99 0 : a(kj) = 1.0e0 - b(kj)
100 410007 : go to 30
101 410007 : 20 a(kj) = 1.0e0
102 410007 : eabsv = abs(e)
103 410007 : b(kj) = s*sqrt(eabsv/ (2*cs*cs+e))
104 : 30 continue
105 : c 30 do 40 l = 1,4
106 :
107 : c kj=338 for Fe l=0 and e=-333.
108 : c a(kj)=1.0
109 :
110 2050035 : do 40 l = 1,4
111 1640028 : k = kj - l
112 1640028 : a(k) = a(kj)
113 1640028 : b(k) = b(kj)
114 410007 : 40 continue
115 :
116 : c a(k)=1.0 for k=kj-4,..,kj.
117 :
118 2870049 : do 70 i = 1,6
119 2460042 : k = kj + 1
120 2460042 : r = rn*d** (n-k)
121 14760252 : do 60 l = 1,5
122 12300210 : k = k - 1
123 12300210 : r = r*d
124 12300210 : rp21 = cis* (e*r-vr(k))
125 12300210 : rp12 = -rp21 - 2.0e0*cs*r
126 12300210 : da(l) = h3* (q11*a(k)+rp12*b(k))
127 12300210 : db(l) = h3* (q22*b(k)+rp21*a(k))
128 2460042 : 60 continue
129 2460042 : l = kj - 1
130 2460042 : a(l) = a(kj) + f1*da1 + f2*da2 + f4*da4 - (f3*da3+f5*da(5))
131 2460042 : b(l) = b(kj) + f1*db1 + f2*db2 + f4*db4 - (f3*db3+f5*db(5))
132 2460042 : l = l - 1
133 2460042 : a(l) = a(kj) + g1*da1 + g2*da2 + g3*da3 + g4*da4 - g5*da(5)
134 2460042 : b(l) = b(kj) + g1*db1 + g2*db2 + g3*db3 + g4*db4 - g5*db(5)
135 2460042 : l = l - 1
136 : a(l) = a(kj) + 1.0125e0*da1 + 3.825e0*da2 + 2.7e0*da3 +
137 2460042 : + 1.575e0*da4 - .1125e0*da(5)
138 : b(l) = b(kj) + 1.0125e0*db1 + 3.825e0*db2 + 2.7e0*db3 +
139 2460042 : + 1.575e0*db4 - .1125e0*db(5)
140 2460042 : l = l - 1
141 : a(l) = a(kj) + fg1*da1 + fg2*da2 + fg3*da3 + fg2*da4 +
142 2460042 : + fg1*da(5)
143 : b(l) = b(kj) + fg1*db1 + fg2*db2 + fg3*db3 + fg2*db4 +
144 2460042 : + fg1*db(5)
145 410007 : 70 continue
146 :
147 : c a(kj)=1.0, a(kj-1),..,a(kj-4) have been obtained.
148 :
149 : c adams' modification here on.
150 :
151 410007 : k = kj - 3
152 33824528 : r = rn*d** (n-k)
153 33824528 : 80 k = k - 1
154 33824528 : r = r*d
155 33824528 : rp21 = cis* (e*r-vr(k))
156 33824528 : rp12 = -rp21 - 2.0e0*cs*r
157 :
158 : cta+
159 : c adams' predictor.
160 33824528 : atk=a(k+1) + t55*da4-t59*da3+t37*da2-t9*da1
161 33824528 : btk=b(k+1) + t55*db4-t59*db3+t37*db2-t9*db1
162 :
163 33824528 : da1 = da2
164 33824528 : da2 = da3
165 33824528 : da3 = da4
166 :
167 33824528 : db1 = db2
168 33824528 : db2 = db3
169 33824528 : db3 = db4
170 :
171 : cta-
172 33824528 : da4=h3*(q11*atk+rp12*btk)
173 33824528 : db4=h3*(q22*btk+rp21*atk)
174 :
175 : cta+
176 : c adams interpolation formula 25.5.5 (Handbook of math.functions).
177 33824528 : a(k)=a(k+1) + c9*da4+c19*da3-c5*da2+c1*da1
178 33824528 : b(k)=b(k+1) + c9*db4+c19*db3-c5*db2+c1*db1
179 : cta-
180 :
181 33824528 : da4 = h3* (q11*a(k)+rp12*b(k))
182 33824528 : db4 = h3* (q22*b(k)+rp21*a(k))
183 :
184 33824528 : if(k.gt.ki) go to 80
185 :
186 410007 : return
187 : END SUBROUTINE
188 : end
189 :
190 :
|
