vacuz.f 4.12 KB
 Markus Betzinger committed Apr 26, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 `````` MODULE m_vacuz CONTAINS SUBROUTINE vacuz( > e,vz,vz0,nmz,dz, < uz,duz,u) use m_juDFT c********************************************************************* c integrates the vacuum wavefunction fo energy e<0 inward from the c last mesh point using the schrodinger equation and assuming that c this point is past the classical turning point. ! (vz(i),i = 1,nmz) contains the potential on a linear mesh. c vz0 is the vacuum zero. c dz is the spacing of the mesh points. c uz,duz are the value and outward normal derivative (i.e. into c the bulk) evaluated at the first mesh point. c (u(i),i=1,nmz) contains the normalized wavefunction c based on code by m. weinert c********************************************************************* USE m_intgr, ONLY : intgz0 IMPLICIT NONE C .. C .. Scalar Arguments .. INTEGER, INTENT (IN) :: nmz REAL, INTENT (IN) :: dz,e,vz0 REAL, INTENT (OUT):: duz,uz C .. C .. Array Arguments .. REAL, INTENT (IN) :: vz(nmz) REAL, INTENT (OUT):: u(nmz) C .. C .. Local Scalars .. REAL a,eps,eru,erw,h,sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,tnorm,u1,u10, + u1p,vme,w1,w10,w1p,yn,zn,znm INTEGER i,it,n,n1 LOGICAL tail C .. C .. Local Arrays .. REAL w(nmz),wp(nmz) C .. C .. Intrinsic Functions .. INTRINSIC abs,exp,sqrt C .. C .. Data statements .. DATA eps/1.e-06/ C .. IF (e.GE.vz0) THEN WRITE (6,*) 'e>vz0; e=',e,'vz0=',vz0 CALL juDFT_error("e>vz0",calledby ="vacuz",hint + ="Vacuum energy parameter too high") ENDIF c--- > set up initial conditions znm = (nmz-1)*dz a = sqrt(2.* (vz0-e)) u(nmz) = exp(-a*znm) w(nmz) = a*u(nmz) wp(nmz) = a*a*u(nmz) c--- > use 4nd order runge-kutta to first few mesh points h = dz DO i = 1,4 n = nmz + 1 - i yn = u(n) zn = w(n) sk1 = h*zn sl1 = h*wp(n) sk2 = h* (zn+0.5*sl1) sl2 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk1)) sk3 = h* (zn+0.5*sl2) sl3 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk2)) sk4 = h* (zn+sl3) sl4 = h* (2.* (vz(n-1)-e)* (yn+sk3)) u(n-1) = yn + (sk1+2.*sk2+2.*sk3+sk4)/6. w(n-1) = zn + (sl1+2.*sl2+2.*sl3+sl4)/6. wp(n-1) = 2.* (vz(n-1)-e)*u(n-1) ENDDO c--- > use adams-bashforth-moulton predictor-corrector for rest DO i = 5,nmz - 1 c--- > adams-bashforth predictor (order h**5) n = nmz + 1 - i n1 = n - 1 u10 = u(n) + h* (55.*w(n)-59.*w(n+1)+37.*w(n+2)-9.*w(n+3))/24. w10 = w(n) + h* (55.*wp(n)-59.*wp(n+1)+37.*wp(n+2)-9.*wp(n+3))/ + 24. c--- > evaluate derivative at next point( n-1 corresponds to m+1) vme = 2.* (vz(n-1)-e) u1p = w10 w1p = vme*u10 c--- > adams-moulton correctors (order h**6) DO it = 1,10 u1 = u(n) + h* (251.*u1p+646.*w(n)-264.*w(n+1)+106.*w(n+2)- + 19.*w(n+3))/720. w1 = w(n) + h* (251.*w1p+646.*wp(n)-264.*wp(n+1)+ + 106.*wp(n+2)-19.*wp(n+3))/720. c--- > final evaluation u1p = w1 w1p = vme*u1 c--- > test quality of corrector and iterate if necessary eru = abs(u1-u10)/ (abs(u1)+abs(h*u1p)) erw = abs(w1-w10)/ (abs(w1)+abs(h*w1p)) IF (eru.LT.eps .AND. erw.LT.eps) EXIT u10 = u1 w10 = w1 ENDDO IF (it>10) THEN WRITE (6,FMT=8000) n1,eru,erw 8000 FORMAT (' ***vacuz - step may be too big - mesh point',i5, + ', eru,erw=',1p,2e16.7) ENDIF c--- > store values u(n1) = u1 w(n1) = w1 wp(n1) = w1p ENDDO c--- > normalize function tail = .true. DO i = 1,nmz wp(i) = u(nmz+1-i)*u(nmz+1-i) ENDDO CALL intgz0(wp,dz,nmz,tnorm,tail) tnorm = 1./sqrt(tnorm) DO i = 1,nmz u(i) = tnorm*u(i) ENDDO uz = u(1) duz = tnorm*w(1) RETURN END SUBROUTINE END ``````