xcxal.f90 4.87 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

Matthias Redies's avatar
Matthias Redies committed
7 8
MODULE m_xcxal
   use m_juDFT
9
!-----------------------------------------------------------------------
Matthias Redies's avatar
Matthias Redies committed
10
!     Called in case of icorr=0 : nonspin-polarized exchange-correlation
11
!       potential with the X-alpha method
Matthias Redies's avatar
Matthias Redies committed
12 13

!     We caclulate  V_x = alpha V_{x,Slater}  where alpha = 2/3
14 15
!     i.e. the output corresponds to the Gaspar result for the
!     local approximation to the Hartree-Fock method.
Matthias Redies's avatar
Matthias Redies committed
16 17

!     krla=1: Relativistic correction of exchange energy and potential
18 19
!             related to Dirac kinetic energy, according to:
!             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
Matthias Redies's avatar
Matthias Redies committed
20

21
!     be careful: calculation in rydberg!
Matthias Redies's avatar
Matthias Redies committed
22

23 24
!     vxcxal   calculates the XC-potential and
!     excxal   calculates the XC-energy
Matthias Redies's avatar
Matthias Redies committed
25

26 27 28
!     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
!-----------------------------------------------------------------------

Matthias Redies's avatar
Matthias Redies committed
29 30 31
   USE m_constants, ONLY : pi_const
   USE m_relcor
   IMPLICIT NONE
32

Matthias Redies's avatar
Matthias Redies committed
33 34 35
   REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
   REAL, PARAMETER, PRIVATE :: thrd = one/three , d_15 = 1.e-15
   REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422  ! 2 * ( 3/(2*pi) )^(2/3)
36

Matthias Redies's avatar
Matthias Redies committed
37
CONTAINS
38
!************************************************************************
Matthias Redies's avatar
Matthias Redies committed
39 40 41 42
   SUBROUTINE vxcxal( &
      krla,jspins, &
      mgrid,ngrid,rh, &
      vx,vxc)
43
!************************************************************************
Matthias Redies's avatar
Matthias Redies committed
44

45 46 47 48
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: jspins
      INTEGER, INTENT (IN) :: krla        !  switch for rel. corr.
      INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
Matthias Redies's avatar
Matthias Redies committed
49

50 51 52 53
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: vx (mgrid,jspins)
      REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
Matthias Redies's avatar
Matthias Redies committed
54

55 56
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
57 58 59 60

      REAL ::  rs, rho, thfpi
      INTEGER :: ir

61 62 63 64
!-----s Intrinsic Functions
      INTRINSIC alog,max

      thfpi  = three / ( four * pi_const )
Matthias Redies's avatar
Matthias Redies committed
65

66
!-----> evaluate relativistic corrections for exchange
Matthias Redies's avatar
Matthias Redies committed
67

68
      ALLOCATE ( psi(ngrid) )
Matthias Redies's avatar
Matthias Redies committed
69 70 71
      CALL relcor( &
         mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
         psi)
72 73

      DO ir = 1,ngrid
Matthias Redies's avatar
Matthias Redies committed
74 75
         IF (jspins == 1) THEN
            rho = max(d_15,rh(ir,1))
76
         ELSE
Matthias Redies's avatar
Matthias Redies committed
77
            rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
78 79
         ENDIF
         rs = (thfpi/rho)**thrd
Matthias Redies's avatar
Matthias Redies committed
80
         ! exc:  exchange energy = -0.9163306/rs
81
         vxc(ir,1) = - psi(ir)*cvx/rs  ! paramagnetic exchange potential = 4/3 exc
Matthias Redies's avatar
Matthias Redies committed
82

83
         vx (ir,1) = - psi(ir)*cvx/rs
Matthias Redies's avatar
Matthias Redies committed
84

85
      ENDDO
Matthias Redies's avatar
Matthias Redies committed
86
      IF ( jspins == 2 ) THEN        ! spinpolarized calculation
87 88
         DO ir = 1,ngrid
            vxc(ir,jspins) = vxc(ir,1)
Matthias Redies's avatar
Matthias Redies committed
89

90 91 92 93 94 95 96
            vx(ir,jspins)  = vx(ir,1)
         ENDDO
      ENDIF

      DEALLOCATE (psi)
      RETURN

Matthias Redies's avatar
Matthias Redies committed
97 98 99 100 101 102 103 104
   END SUBROUTINE vxcxal
!***********************************************************************
   SUBROUTINE excxal( &
      iofile,krla,jspins, &
      mgrid,ngrid,rh, &
      exc)
!***********************************************************************

105 106 107 108
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: krla, jspins
      INTEGER, INTENT (IN) :: iofile      !  file number for read and write
      INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
Matthias Redies's avatar
Matthias Redies committed
109

110 111 112
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
Matthias Redies's avatar
Matthias Redies committed
113

114
!     .. Local Scalars ..
Matthias Redies's avatar
Matthias Redies committed
115
      REAL :: cex
116 117 118
      REAL ::  rs, rho, thfpi
      INTEGER :: ir

119 120
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
Matthias Redies's avatar
Matthias Redies committed
121

122 123
!-----> Intrinsic Functions
      INTRINSIC alog,max
Matthias Redies's avatar
Matthias Redies committed
124

125 126 127 128
      cex = three * cvx / four
      thfpi  = three / ( four * pi_const )

      ALLOCATE ( phi(ngrid) )
Matthias Redies's avatar
Matthias Redies committed
129 130 131
      CALL relcor( &
         mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
         phi)
132 133

      DO ir = 1,ngrid
Matthias Redies's avatar
Matthias Redies committed
134 135 136 137 138 139 140 141
         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
         rs = (thfpi/rho)**thrd
         ! exc:  exchange energy = -0.9163306/rs
         exc(ir)     = - phi(ir)*cex / rs
142
      ENDDO
Matthias Redies's avatar
Matthias Redies committed
143 144 145 146 147
      IF ( jspins == 2 ) WRITE (iofile,*) &
         'WARNING: X-alpha XC & spin-polarized calculation'
      IF ( jspins > 2 ) THEN
         WRITE (iofile,'('' or error in jspins, jspins ='',i2)') jspins
         CALL juDFT_error("excxal",calledby="xcxal")
148 149 150 151 152
      END IF

      DEALLOCATE (phi)
      RETURN

Matthias Redies's avatar
Matthias Redies committed
153 154
   END SUBROUTINE excxal
END MODULE m_xcxal