Line data    Source code

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