LCOV - code coverage report
Current view: top level - global - radsrd.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 86 87 98.9 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_radsrd
       2             :       use m_juDFT
       3             : 
       4             : c*********************************************************************
       5             : c     calculates the energy derivative of the scalar relativistic
       6             : c     wavefuction for energy e and angular momentum l.
       7             : c     the large (small) component is returned in pe (qe). for non-
       8             : c     relativistic case, set cin=0.
       9             : c     the functions p,q and the radial derivative dus are required.
      10             : c                 based on code by m. weinert  
      11             : c*********************************************************************
      12             : 
      13             : c*********************************************************************
      14             : c
      15             : c Uncomment the following, if you want to calculate an "exact" udot
      16             : c including relativistic corrections in the inhomogeneous part of the
      17             : c Dirac-DGL. Note that these corrections are inconsistent with the
      18             : c construction of the Hamiltonian in hssphn, where the inhomogeneous
      19             : c part is assumed to be non-relativistic.
      20             : c
      21             : c C. Friedrich June 2009
      22             : c*********************************************************************
      23             : # define RELATIVISTIC_CORRECTIONS
      24             : 
      25             : # ifdef RELATIVISTIC_CORRECTIONS
      26             : #   define CORR *(1.+cin2*fl1/rm**2)
      27             : # else
      28             : #   define CORR
      29             : # endif
      30             : 
      31             :       CONTAINS
      32       29930 :       SUBROUTINE radsrd(
      33       29930 :      >                  e,l,vr,r0,h,jri,c,
      34       29930 :      <                  ud,dud,ddn,nodes,pe,qe,
      35       59860 :      >                  p,q,dus)
      36             : 
      37             :       USE m_intgr, ONLY : intgr0
      38             :       IMPLICIT NONE
      39             : C     ..
      40             : C     .. Scalar Arguments ..
      41             :       INTEGER, INTENT (IN) :: l
      42             :       INTEGER, INTENT (IN) :: jri
      43             :       INTEGER, INTENT (OUT):: nodes
      44             :       REAL,    INTENT (IN) :: c
      45             :       REAL,    INTENT (IN) :: dus,e,h,r0
      46             :       REAL,    INTENT (OUT):: ddn,dud,ud
      47             : C     ..
      48             : C     .. Array Arguments ..
      49             :       REAL,    INTENT (IN) :: p(:),q(:),vr(:)
      50             :       REAL,    INTENT (OUT):: pe(:),qe(:)
      51             : C     ..
      52             : C     .. Local Scalars ..
      53             :       REAL dr,drh,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,
      54             :      +     sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,t1,t2,yn,zn,cin,cin2
      55             :       INTEGER i,it
      56             : C     ..
      57             : C     .. Local Arrays ..
      58       59860 :       REAL pp(size(p)),qp(size(p))
      59             : C     ..
      60             : C     .. Data statements ..
      61             :       REAL,PARAMETER :: eps=1.e-06
      62             : C     ..
      63       29930 :       cin = 1.0/c
      64       29930 :       cin2 = cin*cin
      65             : c
      66       29930 :       IF(jri>size(p))  CALL juDFT_error("dimension too small",
      67           0 :      +     calledby="radsrd")
      68             : c--->    set up initial conditions
      69       29930 :       fl1 = l* (l+1)
      70       29930 :       rm = 2.*r0 + cin2* (r0*e-vr(1))
      71       29930 :       pe(1) = 0
      72       29930 :       qe(1) = 0
      73       29930 :       pp(1) =   r0*cin*q(1)
      74       29930 :       qp(1) = - r0*p(1) CORR
      75             : c--->    use 4th order runge-kutta to get first few mesh points
      76       29930 :       dr = exp(h)
      77       29930 :       drh = sqrt(dr)
      78       29930 :       r = r0
      79      179580 :       DO i = 1,5
      80      149650 :          rh2 = drh*r
      81      149650 :          rh = dr*r
      82      149650 :          sk1 = h*pp(i)
      83      149650 :          sl1 = h*qp(i)
      84      149650 :          rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
      85      149650 :          rm = 2.*rh2 - cin2*rve
      86      149650 :          yn = pe(i) + 0.5*sk1
      87      149650 :          zn = qe(i) + 0.5*sl1
      88      149650 :          t1 = 0.5*rh2*cin* (q(i)+q(i+1))
      89      149650 :          t2 = 0.5*rh2* (p(i)+p(i+1)) CORR
      90      149650 :          sk2 = h* (rm*zn+yn+t1)
      91      149650 :          sl2 = h* ((fl1/rm+rve)*yn-zn-t2)
      92      149650 :          yn = pe(i) + 0.5*sk2
      93      149650 :          zn = qe(i) + 0.5*sl2
      94      149650 :          sk3 = h* (rm*zn+yn+t1)
      95      149650 :          sl3 = h* ((fl1/rm+rve)*yn-zn-t2)
      96      149650 :          rve = vr(i+1) - rh*e
      97      149650 :          rm = 2.*rh - cin2*rve
      98      149650 :          yn = pe(i) + sk3
      99      149650 :          zn = qe(i) + sl3
     100      149650 :          sk4 = h* (rm*zn+yn+rh*cin*q(i+1))
     101      149650 :          sl4 = h* ((fl1/rm+rve)*yn-zn-rh*p(i+1) CORR )
     102      149650 :          pe(i+1) = pe(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
     103      149650 :          qe(i+1) = qe(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
     104      149650 :          pp(i+1) = rm*qe(i+1) + pe(i+1) + rh*cin*q(i+1)
     105             :          qp(i+1) = (fl1/rm+rve)*pe(i+1) - qe(i+1) -
     106      149650 :      +             rh*p(i+1) CORR
     107      179580 :          r = rh
     108             :       ENDDO
     109       29930 :       nodes = 0
     110             : 
     111             : c--->    adams-bashforth-moulton predictor-corrector
     112    17859588 :       predictor: DO i = 6,jri - 1
     113    17829658 :          r = r*dr
     114             : c--->    predictor
     115             :          p0 = pe(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
     116    17829658 :      +        7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
     117             :          q0 = qe(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
     118    17829658 :      +        7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
     119             : c--->    evaluate derivatives at next point
     120    17829658 :          rve = vr(i+1) - r*e
     121    17829658 :          rm = 2.*r - cin2*rve
     122    17829658 :          t1 = cin*r*q(i+1)
     123    17829658 :          t2 = r*p(i+1) CORR
     124    17829658 :          p1p = rm*q0 + p0 + t1
     125    17829658 :          q1p = (fl1/rm+rve)*p0 - q0 - t2
     126             : c--->    corrector
     127    21621899 :          corrector: DO it = 1,5
     128             :             p1 = pe(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
     129    21621899 :      +           482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
     130             :             q1 = qe(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
     131    21621899 :      +           482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
     132             : c--->    final evaluation
     133    21621899 :             p1p = rm*q1 + p1 + t1
     134    21621899 :             q1p = (fl1/rm+rve)*p1 - q1 - t2
     135             : c--->    test quality of corrector and iterate if necessary
     136    21621899 :             erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
     137    21621899 :             erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
     138    21621899 :             IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
     139     3792241 :             p0 = p1
     140     3792241 :             q0 = q1
     141             :          ENDDO corrector
     142             : c--->    store values
     143    17829658 :          pe(i+1) = p1
     144    17829658 :          qe(i+1) = q1
     145    17829658 :          pp(i+1) = p1p
     146    17829658 :          qp(i+1) = q1p
     147    17859588 :          nodes = nodes + 0.501*abs(sign(1.0,pe(i+1))-sign(1.0,pe(i)))
     148             :       ENDDO predictor
     149             : c--->    ensure orthogonality
     150    36048406 :       DO i = 1,jri
     151    18039168 :          qe(i) = cin*qe(i)
     152             :       ENDDO
     153    36048406 :       DO i = 1,jri
     154    18039168 :          qp(i) = p(i)*pe(i) + q(i)*qe(i)
     155             :       ENDDO
     156       29930 :       CALL intgr0(qp,r0,h,jri,t)
     157       29930 :       dud = (pp(jri)-pe(jri))/ (r*r)
     158    18039168 :       DO i = 1,jri
     159    18009238 :          pe(i) = pe(i) - t*p(i)
     160    18039168 :          qe(i) = qe(i) - t*q(i)
     161             :       ENDDO
     162       29930 :       ud = pe(jri)/r
     163       29930 :       dud = dud - t*dus
     164    18039168 :       DO i = 1,jri
     165    18039168 :          qp(i) = pe(i)*pe(i) + qe(i)*qe(i)
     166             :       ENDDO
     167       29930 :       CALL intgr0(qp,r0,h,jri,ddn)
     168             : 
     169       29930 :       END SUBROUTINE radsrd
     170             :       END MODULE m_radsrd

Generated by: LCOV version 1.13