Commit 8cc55f94 authored by Matthias Redies's avatar Matthias Redies

xcxal moved

parent 32c41c5d
...@@ -23,11 +23,11 @@ xc-pot/xcbh.f90 ...@@ -23,11 +23,11 @@ xc-pot/xcbh.f90
xc-pot/xcpz.f90 xc-pot/xcpz.f90
xc-pot/xcvwn.f90 xc-pot/xcvwn.f90
xc-pot/xcwgn.f90 xc-pot/xcwgn.f90
xc-pot/xcxal.f90
) )
set(fleur_F77 ${fleur_F77} set(fleur_F77 ${fleur_F77}
xc-pot/xch91.F xc-pot/xch91.F
xc-pot/xcxal.f
) )
set(inpgen_F90 ${inpgen_F90} set(inpgen_F90 ${inpgen_F90}
...@@ -50,11 +50,11 @@ xc-pot/xcbh.f90 ...@@ -50,11 +50,11 @@ xc-pot/xcbh.f90
xc-pot/xcpz.f90 xc-pot/xcpz.f90
xc-pot/xcvwn.f90 xc-pot/xcvwn.f90
xc-pot/xcwgn.f90 xc-pot/xcwgn.f90
xc-pot/xcxal.f90
) )
set(inpgen_F77 ${inpgen_F77} set(inpgen_F77 ${inpgen_F77}
xc-pot/xch91.F xc-pot/xch91.F
xc-pot/xcxal.f
) )
...@@ -4,91 +4,89 @@ ...@@ -4,91 +4,89 @@
! of the MIT license as expressed in the LICENSE file in more detail. ! of the MIT license as expressed in the LICENSE file in more detail.
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
MODULE m_xcxal MODULE m_xcxal
use m_juDFT use m_juDFT
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
! Called in case of icorr=0 : nonspin-polarized exchange-correlation ! Called in case of icorr=0 : nonspin-polarized exchange-correlation
! potential with the X-alpha method ! potential with the X-alpha method
!
! We caclulate V_x = alpha V_{x,Slater} where alpha = 2/3 ! We caclulate V_x = alpha V_{x,Slater} where alpha = 2/3
! i.e. the output corresponds to the Gaspar result for the ! i.e. the output corresponds to the Gaspar result for the
! local approximation to the Hartree-Fock method. ! local approximation to the Hartree-Fock method.
!
! krla=1: Relativistic correction of exchange energy and potential ! krla=1: Relativistic correction of exchange energy and potential
! related to Dirac kinetic energy, according to: ! related to Dirac kinetic energy, according to:
! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979) ! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
!
! be careful: calculation in rydberg! ! be careful: calculation in rydberg!
!
! vxcxal calculates the XC-potential and ! vxcxal calculates the XC-potential and
! excxal calculates the XC-energy ! excxal calculates the XC-energy
!
! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96 ! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
USE m_constants, ONLY : pi_const USE m_constants, ONLY : pi_const
USE m_relcor USE m_relcor
IMPLICIT NONE IMPLICIT NONE
REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0 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 :: thrd = one/three , d_15 = 1.e-15
REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422 ! 2 * ( 3/(2*pi) )^(2/3) REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422 ! 2 * ( 3/(2*pi) )^(2/3)
CONTAINS CONTAINS
!************************************************************************ !************************************************************************
SUBROUTINE vxcxal( SUBROUTINE vxcxal( &
> krla,jspins, krla,jspins, &
> mgrid,ngrid,rh, mgrid,ngrid,rh, &
< vx,vxc) vx,vxc)
!************************************************************************ !************************************************************************
!
! .. Scalar Arguments .. ! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jspins INTEGER, INTENT (IN) :: jspins
INTEGER, INTENT (IN) :: krla ! switch for rel. corr. INTEGER, INTENT (IN) :: krla ! switch for rel. corr.
INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
!
! .. Array Arguments .. ! .. Array Arguments ..
REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
REAL, INTENT (OUT) :: vx (mgrid,jspins) REAL, INTENT (OUT) :: vx (mgrid,jspins)
REAL, INTENT (OUT) :: vxc(mgrid,jspins) ! xc potential REAL, INTENT (OUT) :: vxc(mgrid,jspins) ! xc potential
!
! .. Local Arrays .. ! .. Local Arrays ..
REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr. REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr.
REAL :: rs, rho, thfpi REAL :: rs, rho, thfpi
INTEGER :: ir INTEGER :: ir
!
!-----s Intrinsic Functions !-----s Intrinsic Functions
INTRINSIC alog,max INTRINSIC alog,max
thfpi = three / ( four * pi_const ) thfpi = three / ( four * pi_const )
!
!-----> evaluate relativistic corrections for exchange !-----> evaluate relativistic corrections for exchange
!
ALLOCATE ( psi(ngrid) ) ALLOCATE ( psi(ngrid) )
CALL relcor( CALL relcor( &
> mgrid,ngrid,jspins,krla,.true.,rh, mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
< psi) psi)
DO ir = 1,ngrid DO ir = 1,ngrid
IF (jspins.EQ.1) THEN IF (jspins == 1) THEN
rho = max(d_15,rh(ir,1)) rho = max(d_15,rh(ir,1))
ELSE ELSE
rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins)) rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
ENDIF ENDIF
rs = (thfpi/rho)**thrd rs = (thfpi/rho)**thrd
! exc: exchange energy = -0.9163306/rs ! exc: exchange energy = -0.9163306/rs
vxc(ir,1) = - psi(ir)*cvx/rs ! paramagnetic exchange potential = 4/3 exc vxc(ir,1) = - psi(ir)*cvx/rs ! paramagnetic exchange potential = 4/3 exc
vx (ir,1) = - psi(ir)*cvx/rs vx (ir,1) = - psi(ir)*cvx/rs
ENDDO ENDDO
IF ( jspins .EQ. 2 ) THEN ! spinpolarized calculation IF ( jspins == 2 ) THEN ! spinpolarized calculation
DO ir = 1,ngrid DO ir = 1,ngrid
vxc(ir,jspins) = vxc(ir,1) vxc(ir,jspins) = vxc(ir,1)
vx(ir,jspins) = vx(ir,1) vx(ir,jspins) = vx(ir,1)
ENDDO ENDDO
ENDIF ENDIF
...@@ -96,62 +94,61 @@ ...@@ -96,62 +94,61 @@
DEALLOCATE (psi) DEALLOCATE (psi)
RETURN RETURN
END SUBROUTINE vxcxal END SUBROUTINE vxcxal
C*********************************************************************** !***********************************************************************
SUBROUTINE excxal( SUBROUTINE excxal( &
> iofile,krla,jspins, iofile,krla,jspins, &
> mgrid,ngrid,rh, mgrid,ngrid,rh, &
< exc) exc)
C*********************************************************************** !***********************************************************************
!
! .. Scalar Arguments .. ! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: krla, jspins INTEGER, INTENT (IN) :: krla, jspins
INTEGER, INTENT (IN) :: iofile ! file number for read and write INTEGER, INTENT (IN) :: iofile ! file number for read and write
INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
!
! .. Array Arguments .. ! .. Array Arguments ..
REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
REAL, INTENT (OUT) :: exc(mgrid) ! xc energy REAL, INTENT (OUT) :: exc(mgrid) ! xc energy
!
! .. Local Scalars .. ! .. Local Scalars ..
REAL cex REAL :: cex
REAL :: rs, rho, thfpi REAL :: rs, rho, thfpi
INTEGER :: ir INTEGER :: ir
!
! .. Local Arrays .. ! .. Local Arrays ..
REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct. REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct.
!
!-----> Intrinsic Functions !-----> Intrinsic Functions
INTRINSIC alog,max INTRINSIC alog,max
!
cex = three * cvx / four cex = three * cvx / four
thfpi = three / ( four * pi_const ) thfpi = three / ( four * pi_const )
ALLOCATE ( phi(ngrid) ) ALLOCATE ( phi(ngrid) )
CALL relcor( CALL relcor( &
> mgrid,ngrid,jspins,krla,.false.,rh, mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
< phi) phi)
DO ir = 1,ngrid DO ir = 1,ngrid
IF (jspins.EQ.1) THEN IF (jspins == 1) THEN
rho = max(d_15,rh(ir,1)) rho = max(d_15,rh(ir,1))
ELSE ELSE
rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins)) rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
ENDIF ENDIF
rs = (thfpi/rho)**thrd rs = (thfpi/rho)**thrd
! exc: exchange energy = -0.9163306/rs ! exc: exchange energy = -0.9163306/rs
exc(ir) = - phi(ir)*cex / rs exc(ir) = - phi(ir)*cex / rs
ENDDO ENDDO
IF ( jspins.EQ.2 ) WRITE (iofile,*) IF ( jspins == 2 ) WRITE (iofile,*) &
+ 'WARNING: X-alpha XC & spin-polarized calculation' 'WARNING: X-alpha XC & spin-polarized calculation'
IF ( jspins.GT.2 ) THEN IF ( jspins > 2 ) THEN
WRITE (iofile,'('' or error in jspins, jspins ='',i2)') jspins WRITE (iofile,'('' or error in jspins, jspins ='',i2)') jspins
CALL juDFT_error("excxal",calledby="xcxal") CALL juDFT_error("excxal",calledby="xcxal")
END IF END IF
DEALLOCATE (phi) DEALLOCATE (phi)
RETURN RETURN
END SUBROUTINE excxal END SUBROUTINE excxal
END MODULE m_xcxal END MODULE m_xcxal
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment