Line data Source code
1 : MODULE m_vacudz
2 : CONTAINS
3 57052 : SUBROUTINE vacudz(
4 28526 : > e,vz,vz0,nmz,dz,
5 28526 : < udz,dudz,ddnv,ue,
6 28526 : > duz,u)
7 : c*********************************************************************
8 : c integrates the vacuum energy derivative wavefunction ue for the
9 : c energy e>0 using the schrodinger equation. the wavefunction u
10 : c is necessary (obtained from a previous call to vacuz). the
11 : c conventions are the same as in vacuz. ue will be orthogonal to
12 : c u. udz and dudz are the value and normal derivative at the
13 : c vacuum boundary.
14 : c based on code by m. weinert
15 : c*********************************************************************
16 : use m_juDFT
17 : USE m_constants
18 : USE m_intgr, ONLY : intgz0
19 : IMPLICIT NONE
20 : C ..
21 : C .. Scalar Arguments ..
22 : INTEGER, INTENT (IN) :: nmz
23 : REAL, INTENT (IN) :: duz,dz,e,vz0
24 : REAL, INTENT (OUT):: ddnv,dudz,udz
25 : C ..
26 : C .. Array Arguments ..
27 : REAL, INTENT (IN) :: u(nmz),vz(nmz)
28 : REAL, INTENT (OUT):: ue(nmz)
29 : C ..
30 : C .. Local Scalars ..
31 : REAL a,eru,erw,h,sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,tnorm,u1,u10,
32 : + u1p,vme,w1,w10,w1p,yn,zn,znm
33 : INTEGER i,it,n,n1
34 : LOGICAL tail
35 : C ..
36 : C .. Local Arrays ..
37 28526 : REAL we(nmz),wp(nmz)
38 : C ..
39 : C .. Data statements ..
40 : REAL,PARAMETER:: eps=1.e-06
41 : C ..
42 28526 : IF (e.GE.vz0) CALL
43 : + juDFT_error("e >= vz0",calledby ="vacudz",hint =
44 0 : + 'Vacuum energy-parameter too high')
45 : c--- > set up initial conditions
46 28526 : znm = (nmz-1)*dz
47 28526 : a = sqrt(2.* (vz0-e))
48 28526 : ue(nmz) = (a*znm-0.5)*u(nmz)/ (a*a)
49 28526 : we(nmz) = (a*znm-1.5)*u(nmz)/a
50 28526 : wp(nmz) = 2.* (vz(nmz)-e)*ue(nmz) - 2.*u(nmz)
51 : c--- > use 4nd order runge-kutta to first few mesh points
52 28526 : h = dz
53 142630 : DO i = 1,4
54 114104 : n = nmz + 1 - i
55 114104 : yn = ue(n)
56 114104 : zn = we(n)
57 114104 : sk1 = h*zn
58 114104 : sl1 = h*wp(n)
59 114104 : sk2 = h* (zn+0.5*sl1)
60 114104 : sl2 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk1)- (u(n)+u(n-1)))
61 114104 : sk3 = h* (zn+0.5*sl2)
62 114104 : sl3 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk2)- (u(n)+u(n-1)))
63 114104 : sk4 = h* (zn+sl3)
64 114104 : sl4 = h* (2.* (vz(n-1)-e)* (yn+sk3)-2.*u(n-1))
65 114104 : ue(n-1) = yn + (sk1+2.*sk2+2.*sk3+sk4)/6.
66 114104 : we(n-1) = zn + (sl1+2.*sl2+2.*sl3+sl4)/6.
67 142630 : wp(n-1) = 2.* (vz(n-1)-e)*ue(n-1) - 2.*u(n-1)
68 : ENDDO
69 : c--- > use adams-bashforth-moulton predictor-corrector for rest
70 7017396 : DO i = 5,nmz - 1
71 : c--- > adams-bashforth predictor (order h**5)
72 6988870 : n = nmz + 1 - i
73 6988870 : n1 = n - 1
74 : u10 = ue(n) + h* (55.*we(n)-59.*we(n+1)+37.*we(n+2)-
75 6988870 : + 9.*we(n+3))/24.
76 : w10 = we(n) + h* (55.*wp(n)-59.*wp(n+1)+37.*wp(n+2)-
77 6988870 : + 9.*wp(n+3))/24.
78 : c--- > evaluate derivative at next point( n-1 corresponds to m+1)
79 6988870 : vme = 2.* (vz(n-1)-e)
80 6988870 : u1p = w10
81 6988870 : w1p = vme*u10 - 2.*u(n-1)
82 : c--- > adams-moulton correctors: (order h**6)
83 26147998 : DO it = 1,10
84 : u1 = ue(n) + h* (251.*u1p+646.*we(n)-264.*we(n+1)+
85 26147998 : + 106.*we(n+2)-19.*we(n+3))/720.
86 : w1 = we(n) + h* (251.*w1p+646.*wp(n)-264.*wp(n+1)+
87 26147998 : + 106.*wp(n+2)-19.*wp(n+3))/720.
88 : c--- > final evaluation
89 26147998 : u1p = w1
90 26147998 : w1p = vme*u1 - 2.*u(n-1)
91 : c--- > test quality of corrector and iterate if necessary
92 26147998 : eru = abs(u1-u10)/ (abs(u1)+abs(h*u1p))
93 26147998 : erw = abs(w1-w10)/ (abs(w1)+abs(h*w1p))
94 26147998 : IF (eru.LT.eps .AND. erw.LT.eps) EXIT
95 19159128 : u10 = u1
96 26147998 : w10 = w1
97 : ENDDO
98 6988870 : IF (it>10) THEN
99 0 : WRITE (oUnit,FMT=8000) n1,eru,erw
100 : 8000 FORMAT (' ***vacudz - step may be too big - mesh point',i5,
101 : + ', eru,erw=',1p,2e16.7)
102 : ENDIF
103 : c--- > store values
104 6988870 : ue(n1) = u1
105 6988870 : we(n1) = w1
106 7017396 : wp(n1) = w1p
107 : ENDDO
108 : c--- > ensure orthogonality
109 28526 : tail = .true.
110 7160026 : DO i = 1,nmz
111 7160026 : wp(i) = ue(nmz+1-i)*u(nmz+1-i)
112 : ENDDO
113 28526 : CALL intgz0(wp,dz,nmz,tnorm,tail)
114 7160026 : ue(:) = ue(:) - tnorm*u(:)
115 28526 : udz = ue(1)
116 28526 : dudz = we(1) - tnorm*duz
117 : c--- > determine normalization
118 7160026 : DO i = 1,nmz
119 7160026 : wp(i) = ue(nmz+1-i)*ue(nmz+1-i)
120 : enddo
121 28526 : CALL intgz0(wp,dz,nmz,ddnv,tail)
122 28526 : RETURN
123 : END subroutine
124 : END
|
