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

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

Generated by: LCOV version 1.13