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         1159 :    SUBROUTINE relcor( &
      12         1159 :       mgrid,ngrid,jspins,krla,l_psi,rh, &
      13         1159 :       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         1159 :       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     44967364 :          DO ir = 1,ngrid
      73     44967364 :             phsi(ir) = one
      74              :          ENDDO
      75              :       ENDIF
      76              : 
      77         1159 :       RETURN
      78              :    END SUBROUTINE relcor
      79              : END MODULE m_relcor