vacudz.f 4.16 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 `````` MODULE m_vacudz CONTAINS SUBROUTINE vacudz( > e,vz,vz0,nmz,dz, < udz,dudz,ddnv,ue, > duz,u) c********************************************************************* c integrates the vacuum energy derivative wavefunction ue for the c energy e>0 using the schrodinger equation. the wavefunction u c is necessary (obtained from a previous call to vacuz). the c conventions are the same as in vacuz. ue will be orthogonal to c u. udz and dudz are the value and normal derivative at the c vacuum boundary. c based on code by m. weinert c********************************************************************* use m_juDFT USE m_intgr, ONLY : intgz0 IMPLICIT NONE C .. C .. Scalar Arguments .. INTEGER, INTENT (IN) :: nmz REAL, INTENT (IN) :: duz,dz,e,vz0 REAL, INTENT (OUT):: ddnv,dudz,udz C .. C .. Array Arguments .. REAL, INTENT (IN) :: u(nmz),vz(nmz) REAL, INTENT (OUT):: ue(nmz) C .. C .. Local Scalars .. REAL a,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 we(nmz),wp(nmz) C .. C .. Data statements .. REAL,PARAMETER:: eps=1.e-06 C .. IF (e.GE.vz0) CALL + juDFT_error("e >= vz0",calledby ="vacudz",hint = + 'Vacuum energy-parameter too high') c--- > set up initial conditions znm = (nmz-1)*dz a = sqrt(2.* (vz0-e)) ue(nmz) = (a*znm-0.5)*u(nmz)/ (a*a) we(nmz) = (a*znm-1.5)*u(nmz)/a wp(nmz) = 2.* (vz(nmz)-e)*ue(nmz) - 2.*u(nmz) c--- > use 4nd order runge-kutta to first few mesh points h = dz DO i = 1,4 n = nmz + 1 - i yn = ue(n) zn = we(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)- (u(n)+u(n-1))) sk3 = h* (zn+0.5*sl2) sl3 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk2)- (u(n)+u(n-1))) sk4 = h* (zn+sl3) sl4 = h* (2.* (vz(n-1)-e)* (yn+sk3)-2.*u(n-1)) ue(n-1) = yn + (sk1+2.*sk2+2.*sk3+sk4)/6. we(n-1) = zn + (sl1+2.*sl2+2.*sl3+sl4)/6. wp(n-1) = 2.* (vz(n-1)-e)*ue(n-1) - 2.*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 = ue(n) + h* (55.*we(n)-59.*we(n+1)+37.*we(n+2)- + 9.*we(n+3))/24. w10 = we(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 - 2.*u(n-1) c--- > adams-moulton correctors: (order h**6) DO it = 1,10 u1 = ue(n) + h* (251.*u1p+646.*we(n)-264.*we(n+1)+ + 106.*we(n+2)-19.*we(n+3))/720. w1 = we(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 - 2.*u(n-1) 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 (' ***vacudz - step may be too big - mesh point',i5, + ', eru,erw=',1p,2e16.7) ENDIF c--- > store values ue(n1) = u1 we(n1) = w1 wp(n1) = w1p ENDDO c--- > ensure orthogonality tail = .true. DO i = 1,nmz wp(i) = ue(nmz+1-i)*u(nmz+1-i) ENDDO CALL intgz0(wp,dz,nmz,tnorm,tail) ue(:) = ue(:) - tnorm*u(:) udz = ue(1) dudz = we(1) - tnorm*duz c--- > determine normalization DO i = 1,nmz wp(i) = ue(nmz+1-i)*ue(nmz+1-i) enddo CALL intgz0(wp,dz,nmz,ddnv,tail) RETURN END subroutine END``````