LCOV - code coverage report
Current view: top level - global - vacuz.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 58 61 95.1 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_vacuz
       2             :       CONTAINS
       3        3694 :       SUBROUTINE vacuz(
       4        3694 :      >     e,vz,vz0,nmz,dz,
       5        3694 :      <     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_intgr, ONLY : intgz0
      21             :       IMPLICIT NONE
      22             : C     ..
      23             : C     .. Scalar Arguments ..
      24             :       INTEGER, INTENT (IN) :: nmz
      25             :       REAL,    INTENT (IN) :: dz,e,vz0
      26             :       REAL,    INTENT (OUT):: duz,uz
      27             : C     ..
      28             : C     .. Array Arguments ..
      29             :       REAL, INTENT (IN) :: vz(nmz)
      30             :       REAL, INTENT (OUT)::  u(nmz)
      31             : C     ..
      32             : C     .. Local Scalars ..
      33             :       REAL a,eps,eru,erw,h,sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,tnorm,u1,u10,
      34             :      +     u1p,vme,w1,w10,w1p,yn,zn,znm
      35             :       INTEGER i,it,n,n1
      36             :       LOGICAL tail
      37             : C     ..
      38             : C     .. Local Arrays ..
      39        7388 :       REAL w(nmz),wp(nmz)
      40             : C     ..
      41             : C     .. Intrinsic Functions ..
      42             :       INTRINSIC abs,exp,sqrt
      43             : C     ..
      44             : C     .. Data statements ..
      45             :       DATA eps/1.e-06/
      46             : C     ..
      47        3694 :       IF (e.GE.vz0) THEN
      48           0 :          WRITE (6,*) 'e>vz0; e=',e,'vz0=',vz0
      49             :          CALL juDFT_error("e>vz0",calledby ="vacuz",hint
      50           0 :      +        ="Vacuum energy parameter too high")
      51             :       ENDIF
      52             : c---  >    set up initial conditions
      53        3694 :       znm = (nmz-1)*dz
      54        3694 :       a = sqrt(2.* (vz0-e))
      55        3694 :       u(nmz) = exp(-a*znm)
      56        3694 :       w(nmz) = a*u(nmz)
      57        3694 :       wp(nmz) = a*a*u(nmz)
      58             : c---  >    use 4nd order runge-kutta to first few mesh points
      59        3694 :       h = dz
      60       18470 :       DO i = 1,4
      61       14776 :          n = nmz + 1 - i
      62       14776 :          yn = u(n)
      63       14776 :          zn = w(n)
      64       14776 :          sk1 = h*zn
      65       14776 :          sl1 = h*wp(n)
      66       14776 :          sk2 = h* (zn+0.5*sl1)
      67       14776 :          sl2 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk1))
      68       14776 :          sk3 = h* (zn+0.5*sl2)
      69       14776 :          sl3 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk2))
      70       14776 :          sk4 = h* (zn+sl3)
      71       14776 :          sl4 = h* (2.* (vz(n-1)-e)* (yn+sk3))
      72       14776 :          u(n-1) = yn + (sk1+2.*sk2+2.*sk3+sk4)/6.
      73       14776 :          w(n-1) = zn + (sl1+2.*sl2+2.*sl3+sl4)/6.
      74       18470 :          wp(n-1) = 2.* (vz(n-1)-e)*u(n-1)
      75             :       ENDDO
      76             : c---  >    use adams-bashforth-moulton predictor-corrector for rest
      77     1813754 :       DO  i = 5,nmz - 1
      78             : c---  >    adams-bashforth predictor (order h**5)
      79      905030 :          n = nmz + 1 - i
      80      905030 :          n1 = n - 1
      81      905030 :          u10 = u(n) + h* (55.*w(n)-59.*w(n+1)+37.*w(n+2)-9.*w(n+3))/24.
      82             :          w10 = w(n) + h* (55.*wp(n)-59.*wp(n+1)+37.*wp(n+2)-9.*wp(n+3))/
      83      905030 :      +        24.
      84             : c---  >    evaluate derivative at next point( n-1 corresponds to m+1)
      85      905030 :          vme = 2.* (vz(n-1)-e)
      86      905030 :          u1p = w10
      87      905030 :          w1p = vme*u10
      88             : c---  >    adams-moulton correctors (order h**6)
      89     3348834 :          DO it = 1,10
      90             :             u1 = u(n) + h* (251.*u1p+646.*w(n)-264.*w(n+1)+106.*w(n+2)-
      91     3348834 :      +           19.*w(n+3))/720.
      92             :             w1 = w(n) + h* (251.*w1p+646.*wp(n)-264.*wp(n+1)+
      93     3348834 :      +           106.*wp(n+2)-19.*wp(n+3))/720.
      94             : c---  >    final evaluation
      95     3348834 :             u1p = w1
      96     3348834 :             w1p = vme*u1
      97             : c---  >    test quality of corrector and iterate if necessary
      98     3348834 :             eru = abs(u1-u10)/ (abs(u1)+abs(h*u1p))
      99     3348834 :             erw = abs(w1-w10)/ (abs(w1)+abs(h*w1p))
     100     3348834 :             IF (eru.LT.eps .AND. erw.LT.eps) EXIT
     101     2443804 :             u10 = u1
     102     2443804 :             w10 = w1
     103             :          ENDDO
     104      905030 :          IF (it>10) THEN
     105           0 :             WRITE (6,FMT=8000) n1,eru,erw
     106             :  8000       FORMAT (' ***vacuz - step may be too big - mesh point',i5,
     107             :      +           ', eru,erw=',1p,2e16.7)
     108             :          ENDIF
     109             : c---  >    store values
     110      905030 :          u(n1) = u1
     111      905030 :          w(n1) = w1
     112      908724 :          wp(n1) = w1p
     113             :       ENDDO
     114             : c---  >    normalize function
     115        3694 :       tail = .true.
     116      927194 :       DO i = 1,nmz
     117      927194 :          wp(i) = u(nmz+1-i)*u(nmz+1-i)
     118             :       ENDDO
     119        3694 :       CALL intgz0(wp,dz,nmz,tnorm,tail)
     120        3694 :       tnorm = 1./sqrt(tnorm)
     121      927194 :       DO  i = 1,nmz
     122      927194 :          u(i) = tnorm*u(i)
     123             :       ENDDO
     124        3694 :       uz = u(1)
     125        3694 :       duz = tnorm*w(1)
     126        3694 :       RETURN
     127             :       END SUBROUTINE
     128             :       END 

Generated by: LCOV version 1.13