IffGit has a new shared runner for building Docker images in GitLab CI. Visit https://iffgit.fz-juelich.de/examples/ci-docker-in-docker for more details.

relcor.f 2.62 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
      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.EQ.1) THEN      !  evaluate relativistic corrections for exchange

         DO ir = 1,ngrid
           IF (jspins.EQ.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.LT.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.LT.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