LCOV - code coverage report
Current view: top level - global - radsrdn.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 135 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 3 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_radsrdn
       8             :       use m_juDFT
       9             :       CONTAINS
      10           0 :       SUBROUTINE radsrdn(
      11           0 :      >                  e,l,vr,r0,h,jri,c,
      12           0 :      <                  udn,dudn,ddnn,nodedn,fn,fni,
      13           0 :      >                  f0,du0,n)
      14             : C*********************************************************************
      15             : C     Calculates the nth energy derivative of the scalar relativistic
      16             : C     wavefuction for energy e and angular momentum l.
      17             : C
      18             : C     The large and small components
      19             : C       of the nth derivative of f0 are returned in fn,
      20             : C       of the inhomogeneous part of the nth derivative of the Dirac
      21             : C       equation ( approx. n * (n-1)st derivative of f0 ) are returned
      22             : C       in fni.
      23             : C     For non-relativistic case, choose large c ( >> 137.0359895).
      24             : C
      25             : C     f0 (= primitive, 0th derivative) and
      26             : C     the radial derivative du0 of f0,
      27             : C     are required.
      28             : C
      29             : C     The nth derivative of the Dirac equation is solved in radsrdn1.
      30             : C
      31             : C                  C. Friedrich   Apr. 2005
      32             : C*********************************************************************
      33             :       USE m_intgr, ONLY : intgr0
      34             :       IMPLICIT NONE
      35             : C     ..
      36             : C     .. Scalar Arguments ..
      37             :       INTEGER, INTENT (IN) :: l,n
      38             :       INTEGER, INTENT (IN) :: jri
      39             :       INTEGER, INTENT (OUT):: nodedn
      40             :       REAL,    INTENT (IN) :: c
      41             :       REAL,    INTENT (IN) :: du0,e,h,r0
      42             :       REAL,    INTENT (OUT):: ddnn,dudn,udn
      43             : C     ..
      44             : C     .. Array Arguments ..
      45             :       REAL,    INTENT (IN) :: f0(:,:),vr(:)
      46             :       REAL,    INTENT (OUT):: fn(:,:),fni(:,:)
      47             : C     ..
      48             : C     .. Local Scalars ..
      49             :       INTEGER i,j,k
      50             :       REAL scprod,ovlap(9,9),fac
      51             : C     ..
      52             : C     .. Local Arrays ..
      53           0 :       REAL f(size(vr),2,0:n),help(size(vr))
      54             : C
      55             : C
      56           0 :       IF(n<1)  CALL juDFT_error("n<1",calledby ="radsrdn")
      57           0 :       f(:,:,0) = f0(:,:)
      58             : C
      59           0 :       DO i=1,n
      60             : C       Calculate
      61             : C           scprod = < f0(:,:) | f(:,:,i) >
      62             : C                  = - 1/2  SUM  ( i over j )  < f(:,:,j) | f(:,:,i-j) >   
      63             : C                          j=1,i ----fac-----  -----------ovlap---------
      64             : C           (This follows from d^n/de^n <f0|f0> = 0.)
      65           0 :         scprod=0
      66           0 :         DO j=1,int(i/2)
      67           0 :           help(:) = f(:,1,j)*f(:,1,i-j) + f(:,2,j)*f(:,2,i-j)
      68           0 :           CALL intgr0(help,r0,h,jri,ovlap(j,i-j))
      69           0 :           ovlap(i-j,j)=ovlap(j,i-j)
      70           0 :           fac=1        !
      71           0 :           DO k=i-j+1,i !
      72           0 :             fac=fac*k  ! Calculate
      73             :           ENDDO        ! fac = ( i over j ) = i! / j! / (i-j)!
      74           0 :           DO k=2,j     !
      75           0 :             fac=fac/k  !
      76             :           ENDDO        !
      77           0 :           IF(j.eq.i-j) fac=fac/2
      78           0 :           scprod=scprod-fac*ovlap(j,i-j)
      79             :         ENDDO
      80             : C
      81             : C        scprod=0 ! Instead of the nth derivative a linear combination
      82             :                   ! of derivatives of order <=n is calculated.
      83             :                   ! It can be shown, that the udot-coefficient in the LO
      84             :                   ! construction in setabc1lo is always zero with this choice.
      85             :                   ! This reduces the error arising from the inconsistency
      86             :                   ! between an "exact" skalar-relativistic udot and a
      87             :                   ! non-relativistic Hamiltonian. Uncomment this, if you
      88             :                   ! have set RELATIVISTIC_CORRECTIONS in radsrd. But note,
      89             :                   ! that then relativistic corrections are not fully treated
      90             :                   ! in the higher derivatives calculated here.
      91           0 :         CALL getfni(e,l,vr,r0,h,jri,c,i,f,fni)
      92             :         CALL radsrdn1(
      93             :      >             e,l,vr,r0,h,jri,c,
      94             :      <             udn,dudn,ddnn,nodedn,f(:,1,i),f(:,2,i),
      95           0 :      >             fni(:,1),fni(:,2),f0(:,1),f0(:,2),du0,scprod)
      96             :       ENDDO
      97           0 :       fn(:,:) = f(:,:,n)
      98           0 :       help(:) = fni(:,1)
      99           0 :       fni(:,1)=-fni(:,2)
     100           0 :       fni(:,2)= help(:)*c
     101           0 :       END SUBROUTINE radsrdn
     102             : 
     103             : C---------------------------------------------------------------------
     104             : 
     105           0 :       SUBROUTINE radsrdn1(
     106           0 :      >                  e,l,vr,r0,h,jri,c,
     107           0 :      <                  ud,dud,ddn,nodes,pe,qe,
     108           0 :      >                  pi,qi,phom,qhom,dus,scprod)
     109             : C*********************************************************************
     110             : C     Solves the "inhomogeneous" Dirac equation for energy e and angular
     111             : C     momentum l with pi (qi) on one side instead of zero.
     112             : C
     113             : C     The large (small) component is returned in pe (qe).
     114             : C     For non-relativistic case, choose large c ( >> 137.0359895).
     115             : C
     116             : C     The functions
     117             : C       pi,qi,
     118             : C       phom,qhom (= homogeneous solution, 0th derivative,
     119             : C                    as calculated in radsra),
     120             : C     as well as
     121             : C       the radial derivative dus of phom,
     122             : C       and scprod = < p(q)e | p(q)hom >
     123             : C     are required.
     124             : C
     125             : C     Modified from radsrd.    C. Friedrich   Apr. 2005
     126             : C*********************************************************************
     127             : C
     128             :       USE m_intgr, ONLY : intgr0
     129             :       IMPLICIT NONE
     130             : C     ..
     131             : C     .. Scalar Arguments ..
     132             :       INTEGER, INTENT (IN) :: l
     133             :       INTEGER, INTENT (IN) :: jri
     134             :       INTEGER, INTENT (OUT):: nodes
     135             :       REAL,    INTENT (IN) :: c
     136             :       REAL,    INTENT (IN) :: dus,e,h,r0,scprod
     137             :       REAL,    INTENT (OUT):: ddn,dud,ud
     138             : C     ..
     139             : C     .. Array Arguments ..
     140             :       REAL,    INTENT (IN) :: pi(:),qi(:),vr(:)
     141             :       REAL,    INTENT (OUT):: pe(:),qe(:)
     142             :       REAL,    INTENT (IN) :: phom(:),qhom(:)
     143             : C     ..
     144             : C     .. Local Scalars ..
     145             :       REAL dr,drh,eps,erp,erq,fl1,p0,p1,p1p,q0,q1,q1p,r,rh,rh2,rm,rve,
     146             :      +     sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,t,t1,t2,yn,zn,cin,cin2
     147             :       INTEGER i,it
     148             : C     ..
     149             : C     .. Local Arrays ..
     150           0 :       REAL pp(size(pi)),qp(size(pi))
     151             : C     ..
     152             : C     .. Intrinsic Functions ..
     153             :       INTRINSIC abs,exp,sign,sqrt
     154             : C     ..
     155             : C     .. Data statements ..
     156             :       DATA eps/1.e-06/
     157             : C     ..
     158           0 :       cin = 1.0/c
     159           0 :       cin2 = cin*cin
     160             : c
     161           0 :       IF (jri>size(pi))  CALL juDFT_error("dimension too small",  
     162           0 :      +     calledby ="radsrdn")
     163             : c--->    set up initial conditions
     164           0 :       fl1 = l* (l+1)
     165           0 :       pe(1) = 0
     166           0 :       qe(1) = 0
     167           0 :       pp(1) = r0*pi(1)
     168           0 :       qp(1) = r0*qi(1)
     169             : c--->    use 4th order runge-kutta to get first few mesh points
     170           0 :       dr = exp(h)
     171           0 :       drh = sqrt(dr)
     172           0 :       r = r0
     173           0 :       DO i = 1,5
     174           0 :          rh2 = drh*r
     175           0 :          rh = dr*r
     176           0 :          sk1 = h*pp(i)
     177           0 :          sl1 = h*qp(i)
     178           0 :          rve = 0.5* (vr(i)+vr(i+1)) - rh2*e
     179           0 :          rm = 2.*rh2 - cin2*rve
     180           0 :          yn = pe(i) + 0.5*sk1
     181           0 :          zn = qe(i) + 0.5*sl1
     182           0 :          t1 = 0.5*rh2*(pi(i)+pi(i+1))
     183           0 :          t2 = 0.5*rh2*(qi(i)+qi(i+1))
     184           0 :          sk2 = h* (rm*zn + yn              + t1)
     185           0 :          sl2 = h* (-zn   + (fl1/rm+rve)*yn + t2)
     186           0 :          yn = pe(i) + 0.5*sk2
     187           0 :          zn = qe(i) + 0.5*sl2
     188           0 :          sk3 = h* (rm*zn + yn              + t1)
     189           0 :          sl3 = h* (-zn   + (fl1/rm+rve)*yn + t2)
     190           0 :          rve = vr(i+1) - rh*e
     191           0 :          rm = 2.*rh - cin2*rve
     192           0 :          yn = pe(i) + sk3
     193           0 :          zn = qe(i) + sl3
     194           0 :          t1 = rh*pi(i+1)
     195           0 :          t2 = rh*qi(i+1)
     196           0 :          sk4 = h* (rm*zn + yn              + t1 )
     197           0 :          sl4 = h* (-zn   + (fl1/rm+rve)*yn + t2 )
     198           0 :          pe(i+1) = pe(i) + (sk1+2.*sk2+2.*sk3+sk4)/6.
     199           0 :          qe(i+1) = qe(i) + (sl1+2.*sl2+2.*sl3+sl4)/6.
     200           0 :          pp(i+1) = rm*qe(i+1) + pe(i+1)              + t1
     201           0 :          qp(i+1) = -qe(i+1)   + (fl1/rm+rve)*pe(i+1) + t2
     202           0 :          r = rh
     203             :       ENDDO
     204           0 :       nodes = 0
     205             : c--->    adams-bashforth-moulton predictor-corrector
     206           0 :       predictor: DO i = 6,jri - 1
     207           0 :          r = r*dr
     208             : c--->    predictor
     209             :          p0 = pe(i) + h* (4277.*pp(i)-7923.*pp(i-1)+9982.*pp(i-2)-
     210           0 :      +        7298.*pp(i-3)+2877.*pp(i-4)-475.*pp(i-5))/1440.
     211             :          q0 = qe(i) + h* (4277.*qp(i)-7923.*qp(i-1)+9982.*qp(i-2)-
     212           0 :      +        7298.*qp(i-3)+2877.*qp(i-4)-475.*qp(i-5))/1440.
     213             : c--->    evaluate derivatives at next point
     214           0 :          rve = vr(i+1) - r*e
     215           0 :          rm = 2.*r - cin2*rve
     216           0 :          t1 = r*pi(i+1)
     217           0 :          t2 = r*qi(i+1)
     218           0 :          p1p = rm*q0 + p0              + t1
     219           0 :          q1p = -q0   + (fl1/rm+rve)*p0 + t2
     220             : c--->    corrector
     221           0 :          corrector: DO it = 1,5
     222             :             p1 = pe(i) + h* (475.*p1p+1427.*pp(i)-798.*pp(i-1)+
     223           0 :      +           482.*pp(i-2)-173.*pp(i-3)+27.*pp(i-4))/1440.
     224             :             q1 = qe(i) + h* (475.*q1p+1427.*qp(i)-798.*qp(i-1)+
     225           0 :      +           482.*qp(i-2)-173.*qp(i-3)+27.*qp(i-4))/1440.
     226             : c--->    final evaluation
     227           0 :             p1p = rm*q1 + p1              + t1
     228           0 :             q1p = -q1   + (fl1/rm+rve)*p1 + t2
     229             : c--->    test quality of corrector and iterate if necessary
     230           0 :             erp = abs(p1-p0)/ (abs(p1)+abs(h*p1p))
     231           0 :             erq = abs(q1-q0)/ (abs(q1)+abs(h*p1p))
     232           0 :             IF (erp.LT.eps .AND. erq.LT.eps) EXIT corrector
     233           0 :             p0 = p1
     234           0 :             q0 = q1
     235             :          ENDDO corrector
     236           0 :          IF (it > 5)  CALL juDFT_error("Not converged.",calledby
     237           0 :      +        ="radsrdn")
     238             : c--->    store values
     239           0 :          pe(i+1) = p1
     240           0 :          qe(i+1) = q1
     241           0 :          pp(i+1) = p1p
     242           0 :          qp(i+1) = q1p
     243           0 :          nodes = nodes + 0.501*abs(sign(1.0,pe(i+1))-sign(1.0,pe(i)))
     244             :       ENDDO predictor
     245             : c--->    ensure < p(q)hom | p(q)e > = scprod
     246           0 :       DO i = 1,jri
     247           0 :          qe(i) = cin*qe(i)
     248             :       ENDDO
     249           0 :       DO i = 1,jri
     250           0 :          qp(i) = phom(i)*pe(i) + qhom(i)*qe(i)
     251             :       ENDDO
     252           0 :       CALL intgr0(qp,r0,h,jri,t)
     253           0 :       dud = (pp(jri)-pe(jri))/ (r*r)
     254           0 :       DO i = 1,jri
     255           0 :          pe(i) = pe(i) + (scprod-t)*phom(i)
     256           0 :          qe(i) = qe(i) + (scprod-t)*qhom(i)
     257             :       ENDDO
     258           0 :       ud = pe(jri)/r
     259           0 :       dud = dud + (scprod-t)*dus
     260           0 :       DO i = 1,jri
     261           0 :          qp(i) = pe(i)*pe(i) + qe(i)*qe(i)
     262             :       ENDDO
     263           0 :       CALL intgr0(qp,r0,h,jri,ddn)
     264           0 :       RETURN
     265             :       END SUBROUTINE radsrdn1
     266             :       
     267             : C---------------------------------------------------------------------
     268             : 
     269           0 :       SUBROUTINE getfni(e,l,vr,r0,h,jri,c,n,f,fni)
     270             : C*********************************************************************
     271             : C     Calculates the inhomogeneous part of the nth derivative of the
     272             : C     Dirac equation with energy e and angular momentum l.
     273             : C
     274             : C     It is
     275             : C
     276             : C     fni(:,1) = n*f(:,2,n-1)/c (Note that f(:,2)=Q/c)
     277             : C
     278             : C                                     l(l+1)
     279             : C     fni(:,2) =-n*f(:,1,n-1)* ( 1 + -------- )
     280             : C                                    (2Mrc)^2
     281             : C
     282             : C                  n      i           n! l(l+1)
     283             : C               + SUM (-1)   -----------------------------
     284             : C                 i=2        (n-i)! (2*M)^(n+1) r^2 c^(2n)
     285             : C
     286             : C     with M = 1 + (E-V)/(2*c^2).
     287             : C
     288             : C                  C. Friedrich   Apr. 2005
     289             : C*********************************************************************
     290             : C
     291             : 
     292             :       IMPLICIT NONE
     293             : C     ..
     294             : C     .. Scalar Arguments ..
     295             :       INTEGER, INTENT (IN) :: n,jri,l
     296             :       REAL,    INTENT (IN) :: c,e,h,r0
     297             : C     ..
     298             : C     .. Array Arguments ..
     299             :       REAL,    INTENT (IN) :: f(:,:,0:),vr(:)
     300             :       REAL,    INTENT (OUT):: fni(:,:)
     301             : C     ..
     302             : C     .. Local Scalars ..
     303             :       INTEGER i
     304             :       REAL r
     305             : C     ..
     306             : C     .. Local Arrays ..
     307           0 :       REAL hlp(jri),fac(0:n),rmsh(jri)
     308             : C
     309           0 :       IF(n.lt.1)  CALL juDFT_error("getfni: n<1!",calledby="radsrdn")
     310           0 :       r=exp(h)
     311           0 :       rmsh(1)=r0
     312           0 :       DO i=2,jri
     313           0 :         rmsh(i)=rmsh(i-1)*r
     314             :       ENDDO
     315           0 :       fac(0)=1
     316           0 :       DO i=1,n
     317           0 :         fac(i)=fac(i-1)*i
     318             :       ENDDO
     319           0 :       hlp(:)      = 2*rmsh(:) + (e*rmsh(:)-vr(:jri))/c**2
     320           0 :       fni(:jri,1) = n*f(:jri,2,n-1) /c
     321             :       fni(:jri,2) =-n*f(:jri,1,n-1)
     322           0 :      &              * ( 1 + l*(l+1)* (hlp(:)*c)**(-2) )
     323           0 :       hlp(:)      = hlp(:)/rmsh(:)
     324           0 :       DO i=2,n
     325           0 :         r = fac(n)/fac(n-i) * (-1)**i * l*(l+1) / c**(2*i)
     326             :         fni(:jri,2) = fni(:jri,2) +
     327           0 :      &                f(:jri,1,n-i) * r / hlp(:)**(i+1) / rmsh(:)**2
     328             :       ENDDO
     329           0 :       END SUBROUTINE getfni
     330             :       
     331             :       END MODULE m_radsrdn
     332             :       

Generated by: LCOV version 1.13