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

          Line data    Source code
       1             :       MODULE m_differ
       2             :       use m_juDFT
       3             : c
       4             : cdiffer   subroutine for the differential equations
       5             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
       6             : c-----x this version modified to accept energies up to emax=2     x----
       7             : c-----x                              d.d.koelling   7/7/77      x----x
       8             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
       9             :       CONTAINS
      10       20833 :       SUBROUTINE differ(
      11       20833 :      >                  fn,fl,fj,c,z,h,rnot,rn,d,msh,vr,
      12             :      X                  e,
      13       20833 :      <                  a,b,ierr)
      14             :       USE m_inwint
      15             :       USE m_outint
      16             : 
      17             :       IMPLICIT NONE
      18             : C
      19             : C     .. Scalar Arguments ..
      20             :       INTEGER,INTENT(IN)  :: msh
      21             :       INTEGER,INTENT(OUT) :: ierr
      22             :       REAL,   INTENT(IN)  :: fj,fl,fn,c,z,h,rnot,rn,d
      23             :       REAL,   INTENT(IN)  :: vr(msh)
      24             :       REAL,   INTENT(OUT) :: a(msh),b(msh)
      25             :       REAL,   INTENT(INOUT) :: e
      26             : C     ..
      27             : C     .. Local Scalars ..
      28             :       REAL cis,cs,de,del,dg,eabsv,emax,emin,fkap,g,h3,
      29             :      +     qcoef,qqqq,r,ra,rb,rg,rj,s,w,wmin
      30             :       INTEGER k,ki,kj,n,nodes,nqnt,ntimes
      31             :       LOGICAL dbl
      32             : C     ..
      33             : C     .. Local Arrays ..
      34             :       REAL a0(5),b0(5)
      35             : C     ..
      36             : C     ..
      37             : C     .. Intrinsic Functions ..
      38             :       INTRINSIC abs,exp,min,sqrt
      39             : C     ..
      40    16209587 :       a = 0.0
      41    16209587 :       b = 0.0
      42       20833 :       ierr = 0
      43       20833 :       nqnt = fn - fl - 0.99e0
      44       20833 :       n = msh
      45       20833 :       del = 5.e-7
      46       20833 :       emax = 2
      47       20833 :       emin = -z*z/fn**2 - 10.0
      48       20833 :       IF ((e.GE.emax) .OR. (e.LE.emin)) e = 0.5e0* (emax+emin)
      49       20833 :       s      = 2.0e0* (fj-fl)
      50       20833 :       cs     = c*s
      51       20833 :       cis    = 1.0e0/cs
      52       20833 :       fkap   = fj + 0.5e0
      53       20833 :       g      = sqrt(fkap**2- (cis*z)**2)
      54       20833 :       h3     = h/3.0e0
      55       37415 :       ntimes = 0
      56      225766 :    10 ntimes = ntimes + 1
      57      225766 :       IF (ntimes.GT.200) GO TO 140
      58             :       CALL outint(
      59             :      >            msh,e,fkap,cs,cis,s,vr,z,rn,rnot,h,d,
      60             :      <            a0,b0,a,b,
      61      225766 :      <            ki,nodes)
      62      225766 :       IF (nqnt-nodes.EQ.0) GOTO 30
      63       47518 :       IF (nqnt-nodes.GT.0) GOTO 130
      64             : c**** too many nodes
      65       30936 :    20 IF (e.LT.emax) emax = e
      66       30936 :       e = 0.5e0* (e+emin)
      67       30936 :       IF ((emax-emin).GT.del) GO TO 10
      68           0 :       WRITE (6,FMT=8010) nodes,fn,fl,fj,emin,e,emax
      69           0 :       WRITE (6,FMT=8030) vr
      70           0 :       WRITE (6,FMT=8030) a
      71             :       CALL juDFT_error("differ 1: problems with solving dirac equation"
      72      178248 :      +     ,calledby ="differ")
      73             : c**** correct number of nodes
      74      178248 :    30 ra = a(ki)
      75      178248 :       rb = b(ki)
      76             :       CALL inwint(
      77             :      >            e,fl,ki,fkap,cs,cis,s,z,h,d,rn,rnot,msh,vr,
      78             :      X            a,b,
      79      178248 :      <            kj)
      80      178248 :       ra = ra/a(ki)
      81      178248 :       rb = rb/b(ki)
      82    13549131 :       DO 40 k = ki,kj
      83    13370883 :          a(k) = a(k)*ra
      84    13370883 :          b(k) = ra*b(k)
      85      178248 :    40 CONTINUE
      86      178248 :       dg = exp(h*g)
      87      178248 :       rg = rnot**g
      88   121474631 :       DO 50 k = 1,kj
      89   121296383 :          a(k) = a(k)*rg
      90   121296383 :          b(k) = rg*b(k)
      91   121296383 :          rg   = rg*dg
      92      178248 :    50 CONTINUE
      93      178248 :       eabsv = abs(e)
      94      178248 :       qcoef = sqrt(eabsv+eabsv)
      95      178248 :       rj = rnot*d** (kj-1)
      96      178248 :       r = rj
      97      178248 :       rg = a(kj)
      98      178248 :       dg = b(kj)
      99      178248 :       qqqq = min(abs(rg),abs(dg))
     100      178248 :       IF (qqqq.LT.1.0e-25) GO TO 70
     101      178248 :       wmin = (1.e-35)/qqqq
     102      178248 :       k = kj
     103    13437027 :    60 k = k + 1
     104    13437027 :       IF (k.GT.n) GO TO 90
     105    13336117 :       r = r*d
     106    13336117 :       w = exp(qcoef* (rj-r))
     107    13336117 :       IF (w.LT.wmin) GO TO 80
     108    13258779 :       a(k) = w*rg
     109    13258779 :       b(k) = w*dg
     110    13258779 :       GO TO 60
     111       77338 :    70 k = kj + 1
     112     3980817 :    80 IF (k.GT.n) GO TO 90
     113     3903479 :       a(k) = 0.0e0
     114     3903479 :       b(k) = 0.0e0
     115     3903479 :       k = k + 1
     116     3980817 :       GO TO 80
     117      178248 :    90 r = rn
     118      178248 :       w = r* (a(n)**2+b(n)**2)
     119      178248 :       k = n
     120      178248 :       r = r + r
     121      178248 :       rj = 1.0e0/d
     122   138102145 :       dbl = .false.
     123   138102145 :   100 k = k - 1
     124   138102145 :       r = r*rj
     125   138102145 :       dbl = .NOT. dbl
     126   138102145 :       rg = r* (a(k)**2+b(k)**2)
     127   138102145 :       w = w + rg
     128   138102145 :       IF (dbl) w = w + rg
     129   138102145 :       IF (k.GT.2) GO TO 100
     130      178248 :       w = h3* (w+rnot* (a(1)**2+b(1)**2))
     131      178248 :       de = cs*a(ki)*b(ki)* (ra-rb)/ (ra*w)
     132      178248 :       IF (de.GT.0.0e0) emin = e
     133      178248 :       IF (de.LT.0.0e0) emax = e
     134      178248 :       e = e + de
     135      178248 :       IF (abs(de).LT.del) GO TO 110
     136      157415 :       IF ((e.GE.emax) .OR. (e.LE.emin)) e = 0.5e0* (emax+emin)
     137       20833 :       GO TO 10
     138       20833 :   110 w = 1.0e0/sqrt(w)
     139             :       !
     140             :       ! for a consistent definition of the small component in
     141             :       ! the Dirac and scalar-relativistic approximation (SRA)
     142             :       ! we have to multiply the small component with -1;
     143             :       ! then the small component of the Dirac equation
     144             :       ! also fulfills the SRA for s-states
     145             :       !                                               M.B. (July, 2012)
     146    16209587 :       DO 120 k = 1,n
     147    16188754 :          a(k) =  a(k)*w
     148    16188754 :          b(k) = -b(k)*w
     149       20833 :   120 CONTINUE
     150       20833 :       ierr = 0
     151       20833 :       RETURN
     152             : c**** too few nodes
     153       16582 :   130 IF (e.GT.emin) emin = e
     154       16582 :       e = 0.5e0* (e+emax)
     155       16582 :       IF ((emax-emin).GT.del) GO TO 10
     156           0 :       WRITE (6,FMT=8020) nodes,fn,fl,fj,emin,e,emax
     157           0 :       WRITE (6,FMT=8030) vr
     158           0 :       WRITE (6,FMT=8030) a
     159             :       CALL juDFT_error("differ 2: problems with solving dirac equation"
     160           0 :      +     ,calledby ="differ")
     161           0 :   140 WRITE (6,FMT=8000)
     162           0 :       WRITE (6,FMT=8030) fn,fl,fj,emin,emax
     163           0 :       WRITE (6,FMT=8030) e,de
     164           0 :       WRITE (6,FMT=8030) ra,rb,w,a(ki),b(ki)
     165           0 :       WRITE (6,FMT=8000)
     166           0 :       WRITE (6,FMT=8030) vr
     167           0 :       WRITE (6,FMT=8030) a
     168           0 :       ierr = 1
     169             :  8000 FORMAT (/,/,/,/,10x,' too many tries required')
     170             :  8010 FORMAT (/,/,/,/,10x,' too many nodes.',i5,3f4.1,3e15.7)
     171             :  8020 FORMAT (/,/,/,/,10x,' too few nodes. ',i5,3f4.1,3e15.7)
     172             :  8030 FORMAT (10x,5e14.4)
     173             :       END SUBROUTINE differ
     174             :       END MODULE m_differ

Generated by: LCOV version 1.13