xcxal.f 5.15 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.
!--------------------------------------------------------------------------------

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
      MODULE m_xcxal
      use m_juDFT
!-----------------------------------------------------------------------
!     Called in case of icorr=0 : nonspin-polarized exchange-correlation 
!       potential with the X-alpha method
!
!     We caclulate  V_x = alpha V_{x,Slater}  where alpha = 2/3 
!     i.e. the output corresponds to the Gaspar result for the
!     local approximation to the Hartree-Fock method.
!
!     krla=1: Relativistic correction of exchange energy and potential 
!             related to Dirac kinetic energy, according to:
!             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
!
!     be careful: calculation in rydberg!
!
!     vxcxal   calculates the XC-potential and
!     excxal   calculates the XC-energy
!
!     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
!-----------------------------------------------------------------------

      USE m_constants, ONLY : pi_const
      USE m_relcor
      IMPLICIT NONE

      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)

      CONTAINS
!************************************************************************
      SUBROUTINE vxcxal(
     >                  krla,jspins,
     >                  mgrid,ngrid,rh,
     <                  vx,vxc)
!************************************************************************
!
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: jspins
      INTEGER, INTENT (IN) :: krla        !  switch for rel. corr.
      INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
!
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: vx (mgrid,jspins)
      REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
!
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
57 58 59 60 61

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


62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
!
!-----s Intrinsic Functions
      INTRINSIC alog,max

      thfpi  = three / ( four * pi_const )
!
!-----> evaluate relativistic corrections for exchange
!
      ALLOCATE ( psi(ngrid) )
      CALL relcor(
     >            mgrid,ngrid,jspins,krla,.true.,rh,
     <            psi)

      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
         rs = (thfpi/rho)**thrd
                                       ! exc:  exchange energy = -0.9163306/rs
         vxc(ir,1) = - psi(ir)*cvx/rs  ! paramagnetic exchange potential = 4/3 exc
         
         vx (ir,1) = - psi(ir)*cvx/rs
      
      ENDDO
      IF ( jspins .EQ. 2 ) THEN        ! spinpolarized calculation
         DO ir = 1,ngrid
            vxc(ir,jspins) = vxc(ir,1)
            
            vx(ir,jspins)  = vx(ir,1)
         ENDDO
      ENDIF

      DEALLOCATE (psi)
      RETURN

      END SUBROUTINE vxcxal
C***********************************************************************
      SUBROUTINE excxal(
     >                  iofile,krla,jspins,
     >                  mgrid,ngrid,rh,
     <                  exc)
C***********************************************************************
!
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: krla, jspins
      INTEGER, INTENT (IN) :: iofile      !  file number for read and write
      INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
!
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
!
!     .. Local Scalars ..
      REAL cex
118 119 120
      REAL ::  rs, rho, thfpi
      INTEGER :: ir

121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
!
!     .. Local Arrays ..
      REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
!
!-----> Intrinsic Functions
      INTRINSIC alog,max
!
      cex = three * cvx / four
      thfpi  = three / ( four * pi_const )

      ALLOCATE ( phi(ngrid) )
      CALL relcor(
     >            mgrid,ngrid,jspins,krla,.false.,rh,
     <            phi)

      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
        rs = (thfpi/rho)**thrd
                                       ! exc:  exchange energy = -0.9163306/rs
        exc(ir)     = - phi(ir)*cex / rs
      ENDDO
      IF ( jspins.EQ.2 ) WRITE (iofile,*)
     +                'WARNING: X-alpha XC & spin-polarized calculation'
      IF ( jspins.GT.2 ) THEN
          WRITE (iofile,'('' or error in jspins, jspins ='',i2)') jspins
           CALL juDFT_error("excxal",calledby="xcxal")
      END IF

      DEALLOCATE (phi)
      RETURN

      END SUBROUTINE excxal
      END MODULE m_xcxal