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

          Line data    Source code
       1             :       MODULE m_radsra
       2             :       use m_juDFT
       3             : c*********************************************************************
       4             : c     calculates the scalar relativistic wavefuction for energy e and
       5             : c     angular momentum l for the potential vr by integrating outward.
       6             : c     the large (small) component is returned in p (q). for non-
       7             : c     relativistic case, set cin=0.        based on code by m. weinert 
       8             : c*********************************************************************
       9             :       CONTAINS
      10     2135490 :       SUBROUTINE radsra(
      11     2135490 :      >                  e,l,vr,r0,h,jri,c,
      12     2135490 :      <                  us,dus,nodes,p,q)
      13             : 
      14             :       USE m_intgr, ONLY : intgr0
      15             :       IMPLICIT NONE
      16             : C     ..
      17             : C     .. Scalar Arguments ..
      18             :       INTEGER, INTENT (IN) :: l
      19             :       INTEGER, INTENT (IN) :: jri
      20             :       INTEGER, INTENT (OUT):: nodes
      21             :       REAL,    INTENT (IN) :: e,h,r0
      22             :       REAL,    INTENT (IN) :: c
      23             :       REAL,    INTENT (OUT):: dus,us
      24             : C     ..
      25             : C     .. Array Arguments ..
      26             :       REAL,    INTENT (IN) :: vr(:)
      27             :       REAL,    INTENT (OUT):: p(:),q(:)
      28             : C     ..
      29             : C     .. Local Scalars ..
      30             :       REAL dr,drh,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,s,
      31             :      +     sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,yn,zn,cin,cin2
      32             :       INTEGER i,it
      33             : C     ..
      34             : C     .. Local Arrays ..
      35     4270980 :       REAL pp(jri),qp(jri)
      36             : C     ..
      37             : C     .. Data statements ..
      38             :       REAL,PARAMETER  :: eps=1.e-06
      39             : C     ..
      40     2135490 :       cin = 1.0/c
      41     2135490 :       cin2 = cin*cin
      42             : 
      43     2135490 :       IF (jri>SIZE(vr).OR.jri>SIZE(p))  CALL juDFT_error
      44           0 :      +    ("BUG: data dimension in radsra too small",calledby ="radsra")
      45             :    
      46             : 
      47             : c--->    set up initial conditions
      48     2135490 :       fl1 = l* (l+1)
      49     2135490 :       s = sqrt(fl1+1-cin2*vr(1)*vr(1)) - 1.
      50     2135490 :       rm = 2.*r0 + cin2* (r0*e-vr(1))
      51     2135490 :       p(1) = r0** (l+1)
      52     2135490 :       q(1) = s*p(1)/rm
      53     2135490 :       pp(1) = rm*q(1) + p(1)
      54     2135490 :       qp(1) = (fl1/rm+vr(1)-r0*e)*p(1) - q(1)
      55             : c--->    use 4th order runge-kutta to get first few mesh points
      56     2135490 :       dr = exp(h)
      57     2135490 :       drh = sqrt(dr)
      58     2135490 :       r = r0
      59    12812940 :       DO i = 1,5
      60    10677450 :          rh2 = drh*r
      61    10677450 :          rh = dr*r
      62    10677450 :          sk1 = h*pp(i)
      63    10677450 :          sl1 = h*qp(i)
      64    10677450 :          rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
      65    10677450 :          rm = 2.*rh2 - cin2*rve
      66    10677450 :          yn = p(i) + 0.5*sk1
      67    10677450 :          zn = q(i) + 0.5*sl1
      68    10677450 :          sk2 = h* (rm*zn+yn)
      69    10677450 :          sl2 = h* ((fl1/rm+rve)*yn-zn)
      70    10677450 :          yn = p(i) + 0.5*sk2
      71    10677450 :          zn = q(i) + 0.5*sl2
      72    10677450 :          sk3 = h* (rm*zn+yn)
      73    10677450 :          sl3 = h* ((fl1/rm+rve)*yn-zn)
      74    10677450 :          rve = vr(i+1) - rh*e
      75    10677450 :          rm = 2.*rh - cin2*rve
      76    10677450 :          yn = p(i) + sk3
      77    10677450 :          zn = q(i) + sl3
      78    10677450 :          sk4 = h* (rm*zn+yn)
      79    10677450 :          sl4 = h* ((fl1/rm+rve)*yn-zn)
      80    10677450 :          p(i+1) = p(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
      81    10677450 :          q(i+1) = q(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
      82    10677450 :          pp(i+1) = rm*q(i+1) + p(i+1)
      83    10677450 :          qp(i+1) = (fl1/rm+rve)*p(i+1) - q(i+1)
      84    12812940 :          r = rh
      85             :       ENDDO
      86     2135490 :       nodes = 0
      87             : c--->    adams-bashforth-moulton predictor-corrector
      88  1318933968 :       predictor: DO i = 6,jri - 1
      89  1316798478 :          r = r*dr
      90             : c--->    predictor
      91             :          p0 = p(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
      92  1316798478 :      +        7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
      93             :          q0 = q(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
      94  1316798478 :      +        7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
      95             : c--->    evaluate derivatives at next point
      96  1316798478 :          rve = vr(i+1) - r*e
      97  1316798478 :          rm = 2.*r - cin2*rve
      98  1316798478 :          p1p = rm*q0 + p0
      99  1316798478 :          q1p = (fl1/rm+rve)*p0 - q0
     100             : c--->    corrector
     101  1392058640 :          corrector: DO it = 1,5
     102             :             p1 = p(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
     103  1392058640 :      +           482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
     104             :             q1 = q(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
     105  1392058640 :      +           482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
     106             : c--->       final evaluation
     107  1392058640 :             p1p = rm*q1 + p1
     108  1392058640 :             q1p = (fl1/rm+rve)*p1 - q1
     109             : c--->       test quality of corrector and iterate if necessary
     110  1392058640 :             erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
     111  1392058640 :             erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
     112  1392058640 :             IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
     113    75260162 :             p0 = p1
     114    75260162 :             q0 = q1
     115             :          ENDDO corrector
     116             : c--->    store values
     117  1316798478 :          p(i+1) = p1
     118  1316798478 :          q(i+1) = q1
     119  1316798478 :          pp(i+1) = p1p
     120  1316798478 :          qp(i+1) = q1p
     121  1318933968 :          nodes = nodes + 0.501*abs(sign(1.0,p(i+1))-sign(1.0,p(i)))
     122             :       ENDDO predictor
     123             : 
     124             : c--->    normalize function
     125  2661358326 :       DO i = 1,jri
     126  1331746908 :          q(i) = cin*q(i)
     127             :       ENDDO
     128  2661358326 :       DO i = 1,jri
     129  1331746908 :          qp(i) = p(i)*p(i) + q(i)*q(i)
     130             :       ENDDO
     131     2135490 :       CALL intgr0(qp,r0,h,jri,t)
     132     2135490 :       t = 1.0/sqrt(t)
     133     2135490 :       us = t*p(jri)/r
     134     2135490 :       dus = t* (pp(jri)-p(jri))/ (r*r)
     135  1331746908 :       DO i = 1,jri
     136  1329611418 :          p(i) = t*p(i)
     137  1331746908 :          q(i) = t*q(i)
     138             :       ENDDO
     139             : 
     140     2135490 :       END SUBROUTINE radsra
     141             :       END MODULE m_radsra

Generated by: LCOV version 1.13