Line data Source code
1 : MODULE m_differ
2 : use m_juDFT
3 : c
4 : cdiffer subroutine for the differential equations
5 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
6 : c-----x this version modified to accept energies up to emax=2 x----
7 : c-----x d.d.koelling 7/7/77 x----x
8 : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
9 : CONTAINS
10 976914 : SUBROUTINE differ(
11 46673 : > fn,fl,fj,c,z,h,rnot,rn,d,msh,vr,
12 : X e,
13 46673 : < a,b,ierr)
14 : USE m_constants
15 : USE m_inwint
16 : USE m_outint
17 :
18 : IMPLICIT NONE
19 :
20 : C .. Scalar Arguments ..
21 : INTEGER,INTENT(IN) :: msh
22 : INTEGER,INTENT(OUT) :: ierr
23 : REAL, INTENT(IN) :: fj,fl,fn,c,z,h,rnot,rn,d
24 : REAL, INTENT(IN) :: vr(msh)
25 : REAL, INTENT(OUT) :: a(msh),b(msh)
26 : REAL, INTENT(INOUT) :: e
27 :
28 : C .. Local Scalars ..
29 : REAL cis,cs,de,del,dg,eabsv,emax,emin,fkap,g,h3,
30 : + qcoef,qqqq,r,ra,rb,rg,rj,s,w,wmin
31 : INTEGER k,ki,kj,n,nodes,nqnt,ntimes
32 : LOGICAL dbl
33 : CHARACTER(LEN=150) hintString
34 :
35 : C .. Local Arrays ..
36 : REAL a0(5),b0(5)
37 :
38 41729390 : a = 0.0
39 41729390 : b = 0.0
40 46673 : ierr = 0
41 46673 : nqnt = fn - fl - 0.99e0
42 46673 : n = msh
43 46673 : del = 1.e-7
44 46673 : emax = 2
45 46673 : emin = -z*z/fn**2 - 10.0
46 46673 : IF ((e.GE.emax) .OR. (e.LE.emin)) e = 0.5e0* (emax+emin)
47 46673 : s = 2.0e0* (fj-fl)
48 46673 : cs = c*s
49 46673 : cis = 1.0e0/cs
50 46673 : fkap = fj + 0.5e0
51 46673 : g = sqrt(fkap**2- (cis*z)**2)
52 46673 : h3 = h/3.0e0
53 46673 : ntimes = 0
54 520234 : 10 ntimes = ntimes + 1
55 520234 : IF (ntimes.GT.200) GO TO 140
56 : CALL outint(
57 : > msh,e,fkap,cs,cis,s,vr,z,rn,rnot,h,d,
58 : < a0,b0,a,b,
59 520234 : < ki,nodes)
60 520234 : IF (nqnt-nodes.EQ.0) GOTO 30
61 110227 : IF (nqnt-nodes.GT.0) GOTO 130
62 : c**** too many nodes
63 69437 : 20 IF (e.LT.emax) emax = e
64 69437 : e = 0.5e0* (e+emin)
65 69437 : IF ((emax-emin).GT.del) GO TO 10
66 0 : WRITE (oUnit,FMT=8010) nodes,fn,fl,fj,emin,e,emax
67 0 : WRITE (oUnit,FMT=8030) vr
68 0 : WRITE (oUnit,FMT=8030) a
69 0 : WRITE (hintString,'(a,i0,a,i0,a,i0,a)') "The n=",NINT(fn)," l=",
70 0 : + NINT(fl), " state of an atom with Z=", NINT(z),
71 0 : + " seems to be below the reasonable energy range."
72 : CALL juDFT_error("differ 1: problems with solving dirac equation"
73 410007 : + ,calledby ="differ", hint=TRIM(hintString))
74 : c**** correct number of nodes
75 410007 : 30 ra = a(ki)
76 410007 : rb = b(ki)
77 : CALL inwint(
78 : > e,fl,ki,fkap,cs,cis,s,z,h,d,rn,rnot,msh,vr,
79 : X a,b,
80 410007 : < kj)
81 410007 : ra = ra/a(ki)
82 410007 : rb = rb/b(ki)
83 35874563 : DO 40 k = ki,kj
84 35464556 : a(k) = a(k)*ra
85 35464556 : b(k) = ra*b(k)
86 410007 : 40 CONTINUE
87 410007 : dg = exp(h*g)
88 410007 : rg = rnot**g
89 318292276 : DO 50 k = 1,kj
90 317882269 : a(k) = a(k)*rg
91 317882269 : b(k) = rg*b(k)
92 317882269 : rg = rg*dg
93 410007 : 50 CONTINUE
94 410007 : eabsv = abs(e)
95 410007 : qcoef = sqrt(eabsv+eabsv)
96 410007 : rj = rnot*d** (kj-1)
97 410007 : r = rj
98 410007 : rg = a(kj)
99 410007 : dg = b(kj)
100 410007 : qqqq = min(abs(rg),abs(dg))
101 410007 : IF (qqqq.LT.1.0e-25) GO TO 70
102 410007 : wmin = (1.e-35)/qqqq
103 410007 : k = kj
104 36612120 : 60 k = k + 1
105 36612120 : IF (k.GT.n) GO TO 90
106 36402975 : r = r*d
107 36402975 : w = exp(qcoef* (rj-r))
108 36402975 : IF (w.LT.wmin) GO TO 80
109 36202113 : a(k) = w*rg
110 36202113 : b(k) = w*dg
111 36202113 : GO TO 60
112 200862 : 70 k = kj + 1
113 12328574 : 80 IF (k.GT.n) GO TO 90
114 12127712 : a(k) = 0.0e0
115 12127712 : b(k) = 0.0e0
116 12127712 : k = k + 1
117 12537719 : GO TO 80
118 410007 : 90 r = rn
119 410007 : w = r* (a(n)**2+b(n)**2)
120 410007 : k = n
121 410007 : r = r + r
122 410007 : rj = 1.0e0/d
123 410007 : dbl = .false.
124 365392080 : 100 k = k - 1
125 365392080 : r = r*rj
126 365392080 : dbl = .NOT. dbl
127 365392080 : rg = r* (a(k)**2+b(k)**2)
128 365392080 : w = w + rg
129 365392080 : IF (dbl) w = w + rg
130 365392080 : IF (k.GT.2) GO TO 100
131 410007 : w = h3* (w+rnot* (a(1)**2+b(1)**2))
132 410007 : de = cs*a(ki)*b(ki)* (ra-rb)/ (ra*w)
133 410007 : IF (de.GT.0.0e0) emin = e
134 410007 : IF (de.LT.0.0e0) emax = e
135 410007 : e = e + de
136 410007 : IF (abs(de).LT.del) GO TO 110
137 363334 : IF ((e.GE.emax) .OR. (e.LE.emin)) e = 0.5e0* (emax+emin)
138 46673 : GO TO 10
139 46673 : 110 w = 1.0e0/sqrt(w)
140 : !
141 : ! for a consistent definition of the small component in
142 : ! the Dirac and scalar-relativistic approximation (SRA)
143 : ! we have to multiply the small component with -1;
144 : ! then the small component of the Dirac equation
145 : ! also fulfills the SRA for s-states
146 : ! M.B. (July, 2012)
147 41729390 : DO 120 k = 1,n
148 41682717 : a(k) = a(k)*w
149 41682717 : b(k) = -b(k)*w
150 46673 : 120 CONTINUE
151 46673 : ierr = 0
152 46673 : RETURN
153 : c**** too few nodes
154 40790 : 130 IF (e.GT.emin) emin = e
155 40790 : e = 0.5e0* (e+emax)
156 40790 : IF ((emax-emin).GT.del) GO TO 10
157 0 : WRITE (oUnit,FMT=8020) nodes,fn,fl,fj,emin,e,emax
158 0 : WRITE (oUnit,FMT=8030) vr
159 0 : WRITE (oUnit,FMT=8030) a
160 0 : WRITE (hintString,'(a,i0,a,i0,a,i0,a)') "The n=",NINT(fn)," l=",
161 0 : + NINT(fl), " state of an atom with Z=", NINT(z),
162 0 : + " seems to be above the reasonable energy range."
163 : CALL juDFT_error("differ 2: problems with solving dirac equation"
164 0 : + ,calledby ="differ", hint=TRIM(hintString))
165 0 : 140 WRITE (oUnit,FMT=8000)
166 0 : WRITE (oUnit,FMT=8030) fn,fl,fj,emin,emax
167 0 : WRITE (oUnit,FMT=8030) e,de
168 0 : WRITE (oUnit,FMT=8030) ra,rb,w,a(ki),b(ki)
169 0 : WRITE (oUnit,FMT=8000)
170 0 : WRITE (oUnit,FMT=8030) vr
171 0 : WRITE (oUnit,FMT=8030) a
172 0 : ierr = 1
173 : 8000 FORMAT (/,/,/,/,10x,' too many tries required')
174 : 8010 FORMAT (/,/,/,/,10x,' too many nodes.',i5,3f4.1,3e15.7)
175 : 8020 FORMAT (/,/,/,/,10x,' too few nodes. ',i5,3f4.1,3e15.7)
176 : 8030 FORMAT (10x,5e14.4)
177 : END SUBROUTINE differ
178 : END MODULE m_differ
|
