xcwgn.f90 4.94 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
MODULE m_xcwgn
8
!-----------------------------------------------------------------------
Matthias Redies's avatar
Matthias Redies committed
9
!     Called in case of icorr=1 : non-spinpolarized exchange-correlation
10
!                                     from wigners interpolation formula.
Matthias Redies's avatar
Matthias Redies committed
11 12

!     krla=1: Relativistic correction of exchange energy and potential
13 14
!             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
15

16
!     be careful: calculation in rydberg!
Matthias Redies's avatar
Matthias Redies committed
17

18 19
!     vxc   calculates the XC-potential and
!     exc   calculates the XC-energy
Matthias Redies's avatar
Matthias Redies committed
20

21 22 23
!     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
!-----------------------------------------------------------------------

Matthias Redies's avatar
Matthias Redies committed
24 25 26
   USE m_constants, ONLY : pi_const
   USE m_relcor
   IMPLICIT NONE
27

Matthias Redies's avatar
Matthias Redies committed
28 29 30
   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 :: cex = 0.91633058742  ! 3/2 * ( 3/(2*pi) )^(2/3)
31

Matthias Redies's avatar
Matthias Redies committed
32
CONTAINS
33
!************************************************************************
Matthias Redies's avatar
Matthias Redies committed
34 35 36 37
   SUBROUTINE vxcwgn( &
      krla,jspins, &
      mgrid,ngrid,rh, &
      vx,vxc)
38
!************************************************************************
Matthias Redies's avatar
Matthias Redies committed
39

40 41 42 43
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: jspins
      INTEGER, INTENT (IN) :: krla        !  run mode parameters
      INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
Matthias Redies's avatar
Matthias Redies committed
44

45 46 47 48
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: vx(mgrid,jspins)      ! x  potential
      REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
Matthias Redies's avatar
Matthias Redies committed
49

50
!     .. Local Scalars ..
Matthias Redies's avatar
Matthias Redies committed
51
      REAL :: fothrd,vxp,vcp
52 53
      REAL ::  rs, rho, thfpi, exp, ecp
      INTEGER :: ir
Matthias Redies's avatar
Matthias Redies committed
54

55 56
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
Matthias Redies's avatar
Matthias Redies committed
57

58 59 60 61 62
!-----s Intrinsic Functions
      INTRINSIC alog,max

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

64
!-----> evaluate relativistic corrections for exchange
Matthias Redies's avatar
Matthias Redies committed
65

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

      DO ir = 1,ngrid
Matthias Redies's avatar
Matthias Redies committed
72 73
         IF (jspins == 1) THEN
            rho = max(d_15,rh(ir,1))
74
         ELSE
Matthias Redies's avatar
Matthias Redies committed
75
            rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
76 77 78 79 80 81 82 83 84
         ENDIF
         rs  = (thfpi/rho)**thrd
         exp = - cex / rs              ! exchange energy = -0.9163306/rs
         ecp = - 0.88 / (rs + 7.8)     ! correlation energy = -0.88/(rs + 7.8)

         vxp = fothrd * exp            ! exchange potential = 4/3 exp
         vcp = fothrd * ecp + 2.288 / ( rs + 7.8 ) ** 2

         vxc(ir,1) = vcp + vxp * psi(ir)
Matthias Redies's avatar
Matthias Redies committed
85

86 87
         vx(ir,1)  = vxp * psi(ir)
      ENDDO
Matthias Redies's avatar
Matthias Redies committed
88
      IF ( jspins == 2 ) THEN        ! spinpolarized calculation
89 90
         DO ir = 1,ngrid
            vxc(ir,jspins) = vxc(ir,1)
Matthias Redies's avatar
Matthias Redies committed
91

92 93 94 95 96 97 98
            vx(ir,jspins)  = vx(ir,1)
         ENDDO
      ENDIF

      DEALLOCATE (psi)
      RETURN

Matthias Redies's avatar
Matthias Redies committed
99 100 101 102 103 104 105 106
   END SUBROUTINE vxcwgn
!***********************************************************************
   SUBROUTINE excwgn( &
      iofile,krla,jspins, &
      mgrid,ngrid,rh, &
      exc)
!***********************************************************************

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

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

117 118
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
119 120 121

      REAL ::  rs, rho, thfpi, exp, ecp
      INTEGER :: ir
Matthias Redies's avatar
Matthias Redies committed
122

123 124 125
      thfpi  = three / ( four * pi_const )

      ALLOCATE ( phi(ngrid) )
Matthias Redies's avatar
Matthias Redies committed
126 127 128
      CALL relcor( &
         mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
         phi)
129 130

      DO ir = 1,ngrid
Matthias Redies's avatar
Matthias Redies committed
131 132
         IF (jspins == 1) THEN
            rho = max(d_15,rh(ir,1))
133
         ELSE
Matthias Redies's avatar
Matthias Redies committed
134
            rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
135 136 137 138 139 140 141
         ENDIF
         rs  = (thfpi/rho)**thrd
         exp = - cex / rs              ! exchange energy = -0.9163306/rs
         ecp = - 0.88 / (rs + 7.8)     ! correlation energy = -0.88/(rs + 7.8)

         exc(ir) = ecp + exp * phi(ir)
      ENDDO
Matthias Redies's avatar
Matthias Redies committed
142 143 144
      IF (jspins == 2) THEN
         WRITE (iofile,'('' WARNING: Wigner correlation !'', &
                ''only applicable for non-spinpolarized calculations'')')
145 146 147 148 149
      ENDIF

      DEALLOCATE (phi)
      RETURN

Matthias Redies's avatar
Matthias Redies committed
150 151
   END SUBROUTINE excwgn
END MODULE m_xcwgn