differ.f90 4.7 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 35 36 37 38 39 40 41 42 43 44 45 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 142 143 144 145 146 147 148 149 150 151
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_differ
  USE m_juDFT
  !c
  !cdiffer   subroutine for the differential equations
  !c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
  !c-----x this version modified to accept energies up to emax=2     x----
  !c-----x                              d.d.koelling   7/7/77      x----x
  !c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
CONTAINS
  SUBROUTINE differ(fn,fl,fj,c,z,h,rnot,rn,d,msh,vr,e,a,b,ierr)
    USE m_inwint
    USE m_outint
    IMPLICIT NONE

    !     .. Scalar Arguments ..
    INTEGER,INTENT(IN)  :: msh
    INTEGER,INTENT(OUT) :: ierr
    REAL,   INTENT(IN)  :: fj,fl,fn,c,z,h,rnot,rn,d
    REAL,   INTENT(IN)  :: vr(msh)
    REAL,   INTENT(OUT) :: a(msh),b(msh)
    REAL,   INTENT(INOUT) :: e
    !     .. Local Scalars ..
    REAL cis,cs,de,del,dg,eabsv,emax,emin,fkap,g,h3,&
         qcoef,qqqq,r,ra,rb,rg,rj,s,w,wmin
    INTEGER k,ki,kj,nodes,nqnt,ntimes
    LOGICAL dbl
    !     .. Local Arrays ..
    REAL a0(5),b0(5)
    !     .. Intrinsic Functions ..
    INTRINSIC abs,exp,min,sqrt

    a = 0.0
    b = 0.0
    nqnt = fn - fl - 0.99e0
    del = 5.e-7
    emax = 2
    emin = -z*z/fn**2 - 10.0
    IF ((e.GE.emax) .OR. (e.LE.emin)) e = 0.5e0* (emax+emin)
    s      = 2.0e0* (fj-fl)
    cs     = c*s
    cis    = 1.0e0/cs
    fkap   = fj + 0.5e0
    g      = SQRT(fkap**2- (cis*z)**2)
    h3     = h/3.0e0
    ntimes_loop:DO ntimes=1,200
       CALL outint( msh,e,fkap,cs,cis,s,vr,z,rn,rnot,h,d,a0,b0,a,b,ki,nodes)
       IF (nqnt-nodes.GT.0) THEN
          IF (e.GT.emin) emin = e
          e = 0.5e0* (e+emax)
          IF ((emax-emin).GT.del) CYCLE ntimes_loop
          WRITE (6,FMT=8020) nodes,fn,fl,fj,emin,e,emax
          WRITE (6,FMT=8030) vr
          WRITE (6,FMT=8030) a
          CALL juDFT_error("differ 2: problems with solving dirac equation"&
               ,calledby ="differ")
       ENDIF
       IF (nqnt-nodes<0) THEN
          IF (e.LT.emax) emax = e
          e = 0.5e0* (e+emin)
          IF ((emax-emin).GT.del) CYCLE ntimes_loop
          WRITE (6,FMT=8010) nodes,fn,fl,fj,emin,e,emax
          WRITE (6,FMT=8030) vr
          WRITE (6,FMT=8030) a
          CALL juDFT_error("differ 1: problems with solving dirac equation"&
               ,calledby ="differ")
       ENDIF
       !**** correct number of nodes
       ra = a(ki)
       rb = b(ki)
       CALL inwint(e,fl,ki,fkap,cs,cis,s,z,h,d,rn,rnot,msh,vr,a,b,kj)
       ra = ra/a(ki)
       rb = rb/b(ki)
       a(ki:kj) = a(ki:kj)*ra
       b(ki:kj) = ra*b(ki:kj)
       dg = EXP(h*g)
       rg = rnot**g
       DO  k = 1,kj
          a(k) = a(k)*rg
          b(k) = rg*b(k)
          rg   = rg*dg
       ENDDO
       eabsv = ABS(e)
       qcoef = SQRT(eabsv+eabsv)
       rj = rnot*d** (kj-1)
       r = rj
       rg = a(kj)
       dg = b(kj)
       qqqq = MIN(ABS(rg),ABS(dg))
       IF (qqqq.LT.1.0e-25) THEN
          IF (kj+1.LE.SIZE(a)) THEN
             a(kj+1:)=0.0
             b(kj+1:)=0.0
          ENDIF
       ELSE
          wmin = (1.e-35)/qqqq
          DO k=kj+1,SIZE(a)
             r = r*d
             w = EXP(qcoef* (rj-r))
             IF (w.LT.wmin) THEN
                a(k)=0.0
                b(k)=0.0
                EXIT
             ENDIF
             a(k) = w*rg
             b(k) = w*dg
          ENDDO
       ENDIF
       r = rn
       w = r* (a(msh)**2+b(msh)**2)
       r = r + r
       rj = 1.0e0/d
       dbl = .FALSE.
       DO k=msh-1,2,-1
          r = r*rj
          dbl = .NOT. dbl
          rg = r* (a(k)**2+b(k)**2)
          w = w + rg
          IF (dbl) w = w + rg
       ENDDO
       w = h3* (w+rnot* (a(1)**2+b(1)**2))
       de = cs*a(ki)*b(ki)* (ra-rb)/ (ra*w)
       IF (de.GT.0.0e0) emin = e
       IF (de.LT.0.0e0) emax = e
       IF (ABS(de).LT.del) THEN
          e = e+de
          a =  a/SQRT(w)
          b = -b/SQRT(w)
          ierr = 0
          RETURN
       END IF
       e = 0.5e0* (emax+emin)
    END DO ntimes_loop
    WRITE (6,FMT=8000)
    WRITE (6,FMT=8030) fn,fl,fj,emin,emax
    WRITE (6,FMT=8030) e,de
    WRITE (6,FMT=8030) ra,rb,w,a(ki),b(ki)
    WRITE (6,FMT=8000)
    WRITE (6,FMT=8030) vr
    WRITE (6,FMT=8030) a
    ierr = 1
8000 FORMAT (/,/,/,/,10x,' too many tries required')
8010 FORMAT (/,/,/,/,10x,' too many nodes.',i5,3f4.1,3e15.7)
8020 FORMAT (/,/,/,/,10x,' too few nodes. ',i5,3f4.1,3e15.7)
8030 FORMAT (10x,5e14.4)
  END SUBROUTINE differ
END MODULE m_differ