radsra.f 4.45 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
      MODULE m_radsra
      use m_juDFT
c*********************************************************************
c     calculates the scalar relativistic wavefuction for energy e and
c     angular momentum l for the potential vr by integrating outward.
c     the large (small) component is returned in p (q). for non-
c     relativistic case, set cin=0.        based on code by m. weinert 
c*********************************************************************
      CONTAINS
      SUBROUTINE radsra(
     >                  e,l,vr,r0,h,jri,c,
     <                  us,dus,nodes,p,q)

      USE m_intgr, ONLY : intgr0
      IMPLICIT NONE
C     ..
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: l
      INTEGER, INTENT (IN) :: jri
      INTEGER, INTENT (OUT):: nodes
      REAL,    INTENT (IN) :: e,h,r0
      REAL,    INTENT (IN) :: c
      REAL,    INTENT (OUT):: dus,us
C     ..
C     .. Array Arguments ..
      REAL,    INTENT (IN) :: vr(:)
      REAL,    INTENT (OUT):: p(:),q(:)
C     ..
C     .. Local Scalars ..
      REAL dr,drh,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,s,
     +     sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,yn,zn,cin,cin2
      INTEGER i,it
C     ..
C     .. Local Arrays ..
35
      REAL pp(jri),qp(jri)
36 37 38 39 40 41 42
C     ..
C     .. Data statements ..
      REAL,PARAMETER  :: eps=1.e-06
C     ..
      cin = 1.0/c
      cin2 = cin*cin

43 44 45
      IF (jri>SIZE(vr).OR.jri>SIZE(p))  CALL juDFT_error
     +    ("BUG: data dimension in radsra too small",calledby ="radsra")
   
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

c--->    set up initial conditions
      fl1 = l* (l+1)
      s = sqrt(fl1+1-cin2*vr(1)*vr(1)) - 1.
      rm = 2.*r0 + cin2* (r0*e-vr(1))
      p(1) = r0** (l+1)
      q(1) = s*p(1)/rm
      pp(1) = rm*q(1) + p(1)
      qp(1) = (fl1/rm+vr(1)-r0*e)*p(1) - q(1)
c--->    use 4th order runge-kutta to get first few mesh points
      dr = exp(h)
      drh = sqrt(dr)
      r = r0
      DO i = 1,5
         rh2 = drh*r
         rh = dr*r
         sk1 = h*pp(i)
         sl1 = h*qp(i)
         rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
         rm = 2.*rh2 - cin2*rve
         yn = p(i) + 0.5*sk1
         zn = q(i) + 0.5*sl1
         sk2 = h* (rm*zn+yn)
         sl2 = h* ((fl1/rm+rve)*yn-zn)
         yn = p(i) + 0.5*sk2
         zn = q(i) + 0.5*sl2
         sk3 = h* (rm*zn+yn)
         sl3 = h* ((fl1/rm+rve)*yn-zn)
         rve = vr(i+1) - rh*e
         rm = 2.*rh - cin2*rve
         yn = p(i) + sk3
         zn = q(i) + sl3
         sk4 = h* (rm*zn+yn)
         sl4 = h* ((fl1/rm+rve)*yn-zn)
         p(i+1) = p(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
         q(i+1) = q(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
         pp(i+1) = rm*q(i+1) + p(i+1)
         qp(i+1) = (fl1/rm+rve)*p(i+1) - q(i+1)
         r = rh
      ENDDO
      nodes = 0
c--->    adams-bashforth-moulton predictor-corrector
      predictor: DO i = 6,jri - 1
         r = r*dr
c--->    predictor
         p0 = p(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
     +        7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
         q0 = q(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
     +        7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
c--->    evaluate derivatives at next point
         rve = vr(i+1) - r*e
         rm = 2.*r - cin2*rve
         p1p = rm*q0 + p0
         q1p = (fl1/rm+rve)*p0 - q0
c--->    corrector
         corrector: DO it = 1,5
            p1 = p(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
     +           482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
            q1 = q(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
     +           482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
c--->       final evaluation
            p1p = rm*q1 + p1
            q1p = (fl1/rm+rve)*p1 - q1
c--->       test quality of corrector and iterate if necessary
            erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
            erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
            IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
            p0 = p1
            q0 = q1
         ENDDO corrector
c--->    store values
         p(i+1) = p1
         q(i+1) = q1
         pp(i+1) = p1p
         qp(i+1) = q1p
         nodes = nodes + 0.501*abs(sign(1.0,p(i+1))-sign(1.0,p(i)))
      ENDDO predictor

c--->    normalize function
      DO i = 1,jri
         q(i) = cin*q(i)
      ENDDO
      DO i = 1,jri
         qp(i) = p(i)*p(i) + q(i)*q(i)
      ENDDO
      CALL intgr0(qp,r0,h,jri,t)
      t = 1.0/sqrt(t)
      us = t*p(jri)/r
      dus = t* (pp(jri)-p(jri))/ (r*r)
      DO i = 1,jri
         p(i) = t*p(i)
         q(i) = t*q(i)
      ENDDO

      END SUBROUTINE radsra
      END MODULE m_radsra