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

          Line data    Source code
       1             :       MODULE m_inwint
       2             :       CONTAINS
       3   109997280 :       SUBROUTINE inwint(
       4      178248 :      >                  e,fl,ki,fkap,cs,cis,s,z,h,dd,rn,rnot,msh,vr,
       5      178248 :      X                  a,b,
       6             :      <                  kj)
       7             : c
       8             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
       9             : c-----x  perform the  inward integration.   this routine is a     x----
      10             : c-----x  derivative of start2/diff2.   much unnecessary indexing  x----
      11             : c-----x  has been removed among other things.    dale koelling    x----
      12             : c
      13             : c     adams' procedure. more stable for GGA. Feb.20,1998. T.A.
      14             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
      15             :       use m_juDFT
      16             :       IMPLICIT NONE
      17             : 
      18             : c     .. Arguments ..
      19             :       INTEGER,INTENT (IN)  :: msh,ki
      20             :       INTEGER,INTENT (OUT) :: kj
      21             :       REAL,   INTENT (IN)  :: e,fl,fkap,cs,cis,s,z,h,dd,rn,rnot
      22             :       REAL,   INTENT (IN)  :: vr(msh)
      23             : c     ..
      24             :       REAL,   INTENT (INOUT) :: a(msh),b(msh)
      25             : c     ..
      26             : c     .. Local Scalars ..
      27             :       REAL atk,btk,d,da1,da2,da3,da4,db1,db2,db3,db4,eabsv,eps,f1,f2,f3,
      28             :      +     f4,f5,fg1,fg2,fg3,g,g1,g2,g3,g4,g5,h3,q11,q22,r,rp12,rp21
      29             :       INTEGER I,k,l,n
      30             :       REAL :: da(5),db(5)
      31             : c     ..
      32             : c     .. intrinsic functions ..
      33             :       INTRINSIC abs,mod,sqrt
      34             : c     ..
      35             : c     .. equivalences ..
      36             : cta+
      37             :       EQUIVALENCE (da1,da(1)), (da2,da(2)), (da3,da(3)), (da4,da(4)),
      38             :      &  (da5,da(5))
      39             :       EQUIVALENCE (db1,db(1)), (db2,db(2)), (db3,db(3)), (db4,db(4)),
      40             :      &  (db5,db(5))
      41             : 
      42             :       REAL c1,c5,c9,c19, t9,t37,t55,t59, da5,db5
      43             : cta-
      44             : c     ..
      45             : c     .. data statements ..
      46             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
      47             : c-----x this version modified to accept positive energies.        x----
      48             : c-----x                              d.d.koelling   7/7/77      x----x
      49             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
      50             : c---* eps    value used to determine the practical infinity
      51             :       DATA f1/1.045833333e0/,f2/2.691666666e0/,f3/1.1e0/
      52             :       DATA f4/.4416666666e0/,f5/.07916666666e0/
      53             :       DATA g1/.9666666666e0/,g2/4.133333333e0/,g3/.8e0/
      54             :       DATA g4/.1333333333e0/,g5/.03333333333e0/
      55             :       DATA fg1/.9333333333e0/,fg2/4.266666666e0/,fg3/1.6e0/
      56             :       DATA eps/75.0e0/
      57             : c     ..
      58      178248 :       h3 = -h/3.0e0
      59      178248 :       n = msh
      60             : cta+
      61      178248 :       t55=55.e0/8.e0
      62      178248 :       t59=59.e0/8.e0
      63      178248 :       t37=37.e0/8.e0
      64      178248 :       t9=9.e0/8.e0
      65             : 
      66      178248 :       c19=19.e0/8.e0
      67      178248 :       c9=9.e0/8.e0
      68      178248 :       c5=5.e0/8.e0
      69      178248 :       c1=1.e0/8.e0
      70             : cta-
      71             : 
      72             : c     e: energy
      73             : c     ki: the grid point corresponding to classical terning point
      74             : cc        at which potential exceeded the energy.
      75             : cc        ki=249 for Fe and l=0.
      76             : c     kj: undefined at the beginning. finally the farthest grid 
      77             : cc        point which can be considered as infinity for the 
      78             : cc        level to be obtained.
      79             : 
      80      178248 :       d = 1.0e0/dd
      81      178248 :       g = sqrt(fkap**2- (cis*z)**2)
      82      178248 :       q11 = s*fkap - g
      83      178248 :       q22 = -g - s*fkap
      84      178248 :       kj = ki + 10
      85    11410155 :       r = rnot* (dd**kj)
      86             : 
      87    11410155 :    10 kj = kj + 1
      88    11410155 :       IF (kj>msh)  CALL juDFT_error("inward integration failed",calledby
      89             :      +     ="inwint",hint
      90           0 :      +     ="This often indicates that the potential is unphysical")
      91             : 
      92             : c     eps is a big energy. thus if the next statement is fulfilled,
      93             : cc    the grid is considered to be infinitely far.
      94             : 
      95    11410155 :       IF (r* (vr(kj)-e*r).gt.eps) GOTO 20
      96    11231907 :       r = r*dd
      97    11231907 :       IF (kj.lt.n) go to 10
      98           0 :       b(kj) = mod(fl,2.0e0)
      99           0 :       a(kj) = 1.0e0 - b(kj)
     100           0 :       go to 30
     101      178248 :    20 a(kj) = 1.0e0
     102      178248 :       eabsv = abs(e)
     103      178248 :       b(kj) = s*sqrt(eabsv/ (2*cs*cs+e))
     104             :    30 continue
     105             : c  30 do 40 l = 1,4
     106             : 
     107             : c     kj=338 for Fe l=0 and e=-333.
     108             : c     a(kj)=1.0
     109             : 
     110     1604232 :       do 40 l = 1,4
     111      712992 :          k = kj - l
     112      712992 :          a(k) = a(kj)
     113      712992 :          b(k) = b(kj)
     114      178248 :    40 continue
     115             : 
     116             : c     a(k)=1.0 for k=kj-4,..,kj.
     117             : 
     118     2317224 :       do 70 i = 1,6
     119     1069488 :          k = kj + 1
     120     1069488 :          r = rn*d** (n-k)
     121     6416928 :          do 60 l = 1,5
     122     5347440 :             k = k - 1
     123     5347440 :             r = r*d
     124     5347440 :             rp21 = cis* (e*r-vr(k))
     125     5347440 :             rp12 = -rp21 - 2.0e0*cs*r
     126     5347440 :             da(l) = h3* (q11*a(k)+rp12*b(k))
     127     5347440 :             db(l) = h3* (q22*b(k)+rp21*a(k))
     128     1069488 :    60    continue
     129     1069488 :          l = kj - 1
     130     5347440 :          a(l) = a(kj) + f1*da1 + f2*da2 + f4*da4 - (f3*da3+f5*da(5))
     131     5347440 :          b(l) = b(kj) + f1*db1 + f2*db2 + f4*db4 - (f3*db3+f5*db(5))
     132     1069488 :          l = l - 1
     133     1069488 :          a(l) = a(kj) + g1*da1 + g2*da2 + g3*da3 + g4*da4 - g5*da(5)
     134     1069488 :          b(l) = b(kj) + g1*db1 + g2*db2 + g3*db3 + g4*db4 - g5*db(5)
     135     1069488 :          l = l - 1
     136             :          a(l) = a(kj) + 1.0125e0*da1 + 3.825e0*da2 + 2.7e0*da3 +
     137     1069488 :      +          1.575e0*da4 - .1125e0*da(5)
     138             :          b(l) = b(kj) + 1.0125e0*db1 + 3.825e0*db2 + 2.7e0*db3 +
     139     1069488 :      +          1.575e0*db4 - .1125e0*db(5)
     140     1069488 :          l = l - 1
     141             :          a(l) = a(kj) + fg1*da1 + fg2*da2 + fg3*da3 + fg2*da4 +
     142     1069488 :      +          fg1*da(5)
     143             :          b(l) = b(kj) + fg1*db1 + fg2*db2 + fg3*db3 + fg2*db4 +
     144     1069488 :      +          fg1*db(5)
     145      178248 :    70 continue
     146             : 
     147             : c     a(kj)=1.0, a(kj-1),..,a(kj-4) have been obtained.
     148             : 
     149             : c     adams' modification here on.
     150             : 
     151      178248 :       k = kj - 3
     152    12657891 :       r = rn*d** (n-k)
     153    12657891 :    80 k = k - 1
     154    12657891 :       r = r*d
     155    12657891 :       rp21 = cis* (e*r-vr(k))
     156    12657891 :       rp12 = -rp21 - 2.0e0*cs*r
     157             : 
     158             : cta+
     159             : c     adams' predictor.
     160    63289455 :       atk=a(k+1) + t55*da4-t59*da3+t37*da2-t9*da1
     161    63289455 :       btk=b(k+1) + t55*db4-t59*db3+t37*db2-t9*db1
     162             : 
     163    12657891 :       da1 = da2
     164    12657891 :       da2 = da3
     165    12657891 :       da3 = da4
     166             : 
     167    12657891 :       db1 = db2
     168    12657891 :       db2 = db3
     169    12657891 :       db3 = db4
     170             : 
     171             : cta-
     172    12657891 :       da4=h3*(q11*atk+rp12*btk)
     173    12657891 :       db4=h3*(q22*btk+rp21*atk)
     174             : 
     175             : cta+
     176             : c     adams interpolation formula 25.5.5 (Handbook of math.functions).
     177    12657891 :       a(k)=a(k+1) + c9*da4+c19*da3-c5*da2+c1*da1
     178    12657891 :       b(k)=b(k+1) + c9*db4+c19*db3-c5*db2+c1*db1
     179             : cta-
     180             : 
     181    12657891 :       da4 = h3* (q11*a(k)+rp12*b(k))
     182    12657891 :       db4 = h3* (q22*b(k)+rp21*a(k))
     183             : 
     184    12657891 :       if(k.gt.ki) go to 80
     185             : 
     186      178248 :       return
     187             :       END SUBROUTINE
     188             :       end
     189             : 
     190             : 

Generated by: LCOV version 1.13