LCOV - code coverage report
Current view: top level - global - radsra.f (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 79 80 98.8 %
Date: 2024-03-29 04:21:46 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      752592 :       SUBROUTINE radsra(
      11      376296 :      >                  e,l,vr,r0,h,jri,c,
      12      376296 :      <                  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      376296 :       REAL pp(jri),qp(jri)
      36             : C     ..
      37             : C     .. Data statements ..
      38             :       REAL,PARAMETER  :: eps=1.e-06
      39             : C     ..
      40      376296 :       cin = 1.0/c
      41      376296 :       cin2 = cin*cin
      42             : 
      43      376296 :       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      376296 :       fl1 = l* (l+1)
      49      376296 :       s = sqrt(fl1+1-cin2*vr(1)*vr(1)) - 1.
      50      376296 :       rm = 2.*r0 + cin2* (r0*e-vr(1))
      51      376296 :       p(1) = r0** (l+1)
      52      376296 :       q(1) = s*p(1)/rm
      53      376296 :       pp(1) = rm*q(1) + p(1)
      54      376296 :       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      376296 :       dr = exp(h)
      57      376296 :       drh = sqrt(dr)
      58      376296 :       r = r0
      59     2257776 :       DO i = 1,5
      60     1881480 :          rh2 = drh*r
      61     1881480 :          rh = dr*r
      62     1881480 :          sk1 = h*pp(i)
      63     1881480 :          sl1 = h*qp(i)
      64     1881480 :          rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
      65     1881480 :          rm = 2.*rh2 - cin2*rve
      66     1881480 :          yn = p(i) + 0.5*sk1
      67     1881480 :          zn = q(i) + 0.5*sl1
      68     1881480 :          sk2 = h* (rm*zn+yn)
      69     1881480 :          sl2 = h* ((fl1/rm+rve)*yn-zn)
      70     1881480 :          yn = p(i) + 0.5*sk2
      71     1881480 :          zn = q(i) + 0.5*sl2
      72     1881480 :          sk3 = h* (rm*zn+yn)
      73     1881480 :          sl3 = h* ((fl1/rm+rve)*yn-zn)
      74     1881480 :          rve = vr(i+1) - rh*e
      75     1881480 :          rm = 2.*rh - cin2*rve
      76     1881480 :          yn = p(i) + sk3
      77     1881480 :          zn = q(i) + sl3
      78     1881480 :          sk4 = h* (rm*zn+yn)
      79     1881480 :          sl4 = h* ((fl1/rm+rve)*yn-zn)
      80     1881480 :          p(i+1) = p(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
      81     1881480 :          q(i+1) = q(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
      82     1881480 :          pp(i+1) = rm*q(i+1) + p(i+1)
      83     1881480 :          qp(i+1) = (fl1/rm+rve)*p(i+1) - q(i+1)
      84     2257776 :          r = rh
      85             :       ENDDO
      86      376296 :       nodes = 0
      87             : c--->    adams-bashforth-moulton predictor-corrector
      88   278551298 :       predictor: DO i = 6,jri - 1
      89   278175002 :          r = r*dr
      90             : c--->    predictor
      91             :          p0 = p(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
      92   278175002 :      +        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   278175002 :      +        7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
      95             : c--->    evaluate derivatives at next point
      96   278175002 :          rve = vr(i+1) - r*e
      97   278175002 :          rm = 2.*r - cin2*rve
      98   278175002 :          p1p = rm*q0 + p0
      99   278175002 :          q1p = (fl1/rm+rve)*p0 - q0
     100             : c--->    corrector
     101   321265190 :          corrector: DO it = 1,5
     102             :             p1 = p(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
     103   318327024 :      +           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   318327024 :      +           482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
     106             : c--->       final evaluation
     107   318327024 :             p1p = rm*q1 + p1
     108   318327024 :             q1p = (fl1/rm+rve)*p1 - q1
     109             : c--->       test quality of corrector and iterate if necessary
     110   318327024 :             erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
     111   318327024 :             erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
     112   318327024 :             IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
     113    43090188 :             p0 = p1
     114   321265190 :             q0 = q1
     115             :          ENDDO corrector
     116             : c--->    store values
     117   278175002 :          p(i+1) = p1
     118   278175002 :          q(i+1) = q1
     119   278175002 :          pp(i+1) = p1p
     120   278175002 :          qp(i+1) = q1p
     121   278551298 :          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   280809074 :       DO i = 1,jri
     126   280809074 :          q(i) = cin*q(i)
     127             :       ENDDO
     128   280809074 :       DO i = 1,jri
     129   280809074 :          qp(i) = p(i)*p(i) + q(i)*q(i)
     130             :       ENDDO
     131      376296 :       CALL intgr0(qp,r0,h,jri,t)
     132      376296 :       t = 1.0/sqrt(t)
     133      376296 :       us = t*p(jri)/r
     134      376296 :       dus = t* (pp(jri)-p(jri))/ (r*r)
     135   280809074 :       DO i = 1,jri
     136   280432778 :          p(i) = t*p(i)
     137   280809074 :          q(i) = t*q(i)
     138             :       ENDDO
     139             : 
     140      376296 :       END SUBROUTINE radsra
     141             :       END MODULE m_radsra

Generated by: LCOV version 1.14