relcor.f90 2.62 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
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
MODULE m_relcor
!************************************************************************

! calculate the relativistic  corrections for exchange as formulated by
!              A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
! either phi (for xc-energy) or psi (for xc-potential) are calculated
!             (if l_psi=.true. we call from vxc.. and psi is evaluated)

!************************************************************************
CONTAINS
   SUBROUTINE relcor( &
      mgrid,ngrid,jspins,krla,l_psi,rh, &
      phsi)
      IMPLICIT NONE

!     .. Scalar Arguments ..
      INTEGER, INTENT (IN)  :: mgrid,krla,ngrid,jspins
      LOGICAL, INTENT (IN)  :: l_psi

!     .. Array Arguments ..
      REAL,    INTENT (IN)  :: rh(mgrid,jspins)
      REAL,    INTENT (OUT) :: phsi(ngrid)

!     .. Local Parameters ..
      REAL, PARAMETER :: betac = 2.2576918e-2  ! alpha * (3 * pi)^(1/3)
      REAL, PARAMETER :: d_15 = 1.e-15 , d_3 = 1.e-3
      REAL, PARAMETER :: one = 1.0 , three = 3.0 , half = 0.5
      REAL, PARAMETER :: thrhalf = three * half , thrd = one/three
      REAL, PARAMETER :: bs1 = 0.75 , bs2 = 0.45 ,  bf2 = 0.4
      REAL, PARAMETER :: bf1 = 2*thrd

!     .. Locals ..
      INTEGER :: ir
      REAL :: beta              ! Fermi velocity devided by speed of light
      REAL :: rho,eta,xi,betasq

      INTRINSIC max,sqrt

      IF (krla == 1) THEN      !  evaluate relativistic corrections for exchange

         DO ir = 1,ngrid
            IF (jspins == 1) THEN
               rho = max(d_15,rh(ir,1))
            ELSE
               rho = max(d_15,rh(ir,1))+max(d_15,rh(ir,jspins))
            ENDIF
            beta = betac * rho**thrd
            betasq = beta*beta
            eta = sqrt(one+betasq)
            xi = beta + eta

            !----->    If beta.LT.10**(-3) use taylor series of psi,phi with respect to
            !          beta, because of accuracy considerations. Taylor series
            !          implemented is exact up to  beta**5 (see notes S.B.)

            IF (l_psi) THEN
               IF (beta < d_3) THEN
                  phsi(ir) = one - betasq + bs1*beta*betasq - bs2*betasq**2
               ELSE
                  phsi(ir) = half* (-one+three*alog(xi)/ (beta*eta))
               END IF
            ELSE
               IF (beta < d_3) THEN
                  phsi(ir) = one - bf1*betasq + bf2*betasq*betasq
               ELSE
                  phsi(ir) = one - thrhalf*((beta*eta-alog(xi))/betasq)**2
               END IF
            ENDIF
         ENDDO

      ELSE
         DO ir = 1,ngrid
            phsi(ir) = one
         ENDDO
      ENDIF

      RETURN
   END SUBROUTINE relcor
END MODULE m_relcor