Line data    Source code

       1             : MODULE m_relcor
       2             : !************************************************************************
       3             : 
       4             : ! calculate the relativistic  corrections for exchange as formulated by
       5             : !              A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
       6             : ! either phi (for xc-energy) or psi (for xc-potential) are calculated
       7             : !             (if l_psi=.true. we call from vxc.. and psi is evaluated)
       8             : 
       9             : !************************************************************************
      10             : CONTAINS
      11         246 :    SUBROUTINE relcor( &
      12         246 :       mgrid,ngrid,jspins,krla,l_psi,rh, &
      13         246 :       phsi)
      14             :       IMPLICIT NONE
      15             : 
      16             : !     .. Scalar Arguments ..
      17             :       INTEGER, INTENT (IN)  :: mgrid,krla,ngrid,jspins
      18             :       LOGICAL, INTENT (IN)  :: l_psi
      19             : 
      20             : !     .. Array Arguments ..
      21             :       REAL,    INTENT (IN)  :: rh(mgrid,jspins)
      22             :       REAL,    INTENT (OUT) :: phsi(ngrid)
      23             : 
      24             : !     .. Local Parameters ..
      25             :       REAL, PARAMETER :: betac = 2.2576918e-2  ! alpha * (3 * pi)^(1/3)
      26             :       REAL, PARAMETER :: d_15 = 1.e-15 , d_3 = 1.e-3
      27             :       REAL, PARAMETER :: one = 1.0 , three = 3.0 , half = 0.5
      28             :       REAL, PARAMETER :: thrhalf = three * half , thrd = one/three
      29             :       REAL, PARAMETER :: bs1 = 0.75 , bs2 = 0.45 ,  bf2 = 0.4
      30             :       REAL, PARAMETER :: bf1 = 2*thrd
      31             : 
      32             : !     .. Locals ..
      33             :       INTEGER :: ir
      34             :       REAL :: beta              ! Fermi velocity devided by speed of light
      35             :       REAL :: rho,eta,xi,betasq
      36             : 
      37             :       INTRINSIC max,sqrt
      38             : 
      39         246 :       IF (krla == 1) THEN      !  evaluate relativistic corrections for exchange
      40             : 
      41           0 :          DO ir = 1,ngrid
      42           0 :             IF (jspins == 1) THEN
      43           0 :                rho = max(d_15,rh(ir,1))
      44             :             ELSE
      45           0 :                rho = max(d_15,rh(ir,1))+max(d_15,rh(ir,jspins))
      46             :             ENDIF
      47           0 :             beta = betac * rho**thrd
      48           0 :             betasq = beta*beta
      49           0 :             eta = sqrt(one+betasq)
      50           0 :             xi = beta + eta
      51             : 
      52             :             !----->    If beta.LT.10**(-3) use taylor series of psi,phi with respect to
      53             :             !          beta, because of accuracy considerations. Taylor series
      54             :             !          implemented is exact up to  beta**5 (see notes S.B.)
      55             : 
      56           0 :             IF (l_psi) THEN
      57           0 :                IF (beta < d_3) THEN
      58           0 :                   phsi(ir) = one - betasq + bs1*beta*betasq - bs2*betasq**2
      59             :                ELSE
      60           0 :                   phsi(ir) = half* (-one+three*alog(xi)/ (beta*eta))
      61             :                END IF
      62             :             ELSE
      63           0 :                IF (beta < d_3) THEN
      64           0 :                   phsi(ir) = one - bf1*betasq + bf2*betasq*betasq
      65             :                ELSE
      66           0 :                   phsi(ir) = one - thrhalf*((beta*eta-alog(xi))/betasq)**2
      67             :                END IF
      68             :             ENDIF
      69             :          ENDDO
      70             : 
      71             :       ELSE
      72    23057698 :          DO ir = 1,ngrid
      73         246 :             phsi(ir) = one
      74             :          ENDDO
      75             :       ENDIF
      76             : 
      77         246 :       RETURN
      78             :    END SUBROUTINE relcor
      79             : END MODULE m_relcor