From 642a21756493b71e9d76196a848c7f295ef4a784 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 15 Mar 2020 22:37:40 +0100 Subject: [PATCH] restricted formalism --- input/dft | 2 +- ...20_lda_correlation_Levy_Zahariev_shift.f90 | 78 +++++++++++++++++++ ...a_correlation_derivative_discontinuity.f90 | 59 ++++++++++++++ ...FL20_lda_correlation_individual_energy.f90 | 73 +++++++++++++++++ ...N5_lda_correlation_Levy_Zahariev_shift.f90 | 74 ++++++++++++++++++ ...VWN5_lda_correlation_individual_energy.f90 | 75 ++++++++++++++++++ ...ricted_correlation_Levy_Zahariev_shift.f90 | 69 ++++++++++++++++ ...d_correlation_derivative_discontinuity.f90 | 65 ++++++++++++++++ ...stricted_correlation_individual_energy.f90 | 69 ++++++++++++++++ ...d_elda_correlation_Levy_Zahariev_shift.f90 | 48 ++++++++++++ .../restricted_elda_correlation_energy.f90 | 44 +++++++++++ ...ted_elda_correlation_individual_energy.f90 | 52 +++++++++++++ src/eDFT/restricted_individual_energy.f90 | 8 +- ...ed_lda_correlation_Levy_Zahariev_shift.f90 | 45 +++++++++++ ...a_correlation_derivative_discontinuity.f90 | 54 +++++++++++++ ...cted_lda_correlation_individual_energy.f90 | 45 +++++++++++ src/eDFT/select_rung.f90 | 2 +- 17 files changed, 857 insertions(+), 5 deletions(-) create mode 100644 src/eDFT/RMFL20_lda_correlation_Levy_Zahariev_shift.f90 create mode 100644 src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/RMFL20_lda_correlation_individual_energy.f90 create mode 100644 src/eDFT/RVWN5_lda_correlation_Levy_Zahariev_shift.f90 create mode 100644 src/eDFT/RVWN5_lda_correlation_individual_energy.f90 create mode 100644 src/eDFT/restricted_correlation_Levy_Zahariev_shift.f90 create mode 100644 src/eDFT/restricted_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/restricted_correlation_individual_energy.f90 create mode 100644 src/eDFT/restricted_elda_correlation_Levy_Zahariev_shift.f90 create mode 100644 src/eDFT/restricted_elda_correlation_energy.f90 create mode 100644 src/eDFT/restricted_elda_correlation_individual_energy.f90 create mode 100644 src/eDFT/restricted_lda_correlation_Levy_Zahariev_shift.f90 create mode 100644 src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/restricted_lda_correlation_individual_energy.f90 diff --git a/input/dft b/input/dft index f6badeb..1f60135 100644 --- a/input/dft +++ b/input/dft @@ -9,6 +9,6 @@ # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.50000 0.00000 + 0.00000 0.00000 # eKS: maxSCF thresh DIIS n_diis guess_type ortho_type 64 0.0000001 T 15 1 1 diff --git a/src/eDFT/RMFL20_lda_correlation_Levy_Zahariev_shift.f90 b/src/eDFT/RMFL20_lda_correlation_Levy_Zahariev_shift.f90 new file mode 100644 index 0000000..b3dbadf --- /dev/null +++ b/src/eDFT/RMFL20_lda_correlation_Levy_Zahariev_shift.f90 @@ -0,0 +1,78 @@ +subroutine RMFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight,rho,EcLZ) + +! Compute the restricted Marut-Fromager-Loos LDA correlation contribution to Levy-Zahariev shift + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + logical :: LDA_centered = .true. + integer :: iEns + double precision :: EcLZLDA + double precision,allocatable :: aMFL(:,:) + double precision,allocatable :: EcLZeLDA(:) + +! Output variables + + double precision,intent(out) :: EcLZ + + +! Allocation + + allocate(aMFL(3,nEns),EcLZeLDA(nEns)) + +! Parameters for weight-dependent LDA correlation functional + + aMFL(1,1) = -0.0238184d0 + aMFL(2,1) = +0.00540994d0 + aMFL(3,1) = +0.0830766d0 + + aMFL(1,2) = -0.0144633d0 + aMFL(2,2) = -0.0506019d0 + aMFL(3,2) = +0.0331417d0 + +! Compute correlation energy for ground, singly-excited and doubly-excited states + + do iEns=1,nEns + + call restricted_elda_correlation_Levy_Zahariev_shift(nEns,aMFL(:,iEns),nGrid,weight(:),rho(:),EcLZeLDA(iEns)) + + end do + +! LDA-centered functional + + EcLZLDA = 0d0 + + if(LDA_centered) then + + call RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:),EcLZLDA) + + do iEns=1,nEns + + EcLZeLDA(iEns) = EcLZeLDA(iEns) + EcLZLDA - EcLZeLDA(1) + + end do + + end if + +! Weight-denpendent functional for ensembles + + EcLZ = 0d0 + + do iEns=1,nEns + + EcLZ = EcLZ + wEns(iEns)*EcLZeLDA(iEns) + + enddo + +end subroutine RMFL20_lda_correlation_Levy_Zahariev_shift diff --git a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..5b2b0e6 --- /dev/null +++ b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 @@ -0,0 +1,59 @@ +subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,Ec) + +! Compute the restricted version of the eLDA correlation part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + +! Local variables + + integer :: iEns,jEns + double precision,allocatable :: aMFL(:,:) + double precision :: dEc(nEns) + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: Ec(nEns) + +! Allocation + + allocate(aMFL(3,nEns)) + +! Parameters for weight-dependent LDA correlation functional + + aMFL(1,1) = -0.0238184d0 + aMFL(2,1) = +0.00540994d0 + aMFL(3,1) = +0.0830766d0 + + aMFL(1,2) = -0.0144633d0 + aMFL(2,2) = -0.0506019d0 + aMFL(3,2) = +0.0331417d0 + +! Compute correlation energy for ground, singly-excited and doubly-excited states + + do iEns=1,nEns + + call restricted_elda_correlation_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:),dEc(iEns)) + + end do + + Ec(:) = 0d0 + + do iEns=1,nEns + do jEns=1,nEns + + Ec(iEns) = Ec(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dEc(jEns) - dEc(1)) + + end do + end do + +end subroutine RMFL20_lda_correlation_derivative_discontinuity diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 new file mode 100644 index 0000000..ef06070 --- /dev/null +++ b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 @@ -0,0 +1,73 @@ +subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ec) + +! Compute eLDA correlation energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + logical :: LDA_centered = .true. + integer :: iEns + double precision :: EcLDA + double precision,allocatable :: aMFL(:,:) + double precision,allocatable :: EceLDA(:) + +! Output variables + + double precision :: Ec + +! Allocation + + allocate(aMFL(3,nEns),EceLDA(nEns)) + +! Parameters for weight-dependent LDA correlation functional + + aMFL(1,1) = -0.0238184d0 + aMFL(2,1) = +0.00540994d0 + aMFL(3,1) = +0.0830766d0 + + aMFL(1,2) = -0.0144633d0 + aMFL(2,2) = -0.0506019d0 + aMFL(3,2) = +0.0331417d0 + +! Compute correlation energy for ground, singly-excited and doubly-excited states + + do iEns=1,nEns + + call restricted_elda_correlation_individual_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:),rho(:),EceLDA(iEns)) + + end do + +! LDA-centered functional + + EcLDA = 0d0 + + if(LDA_centered) then + + call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) + + do iEns=1,nEns + EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1) + end do + + end if + +! Weight-denpendent functional for ensembles + + Ec = 0d0 + + do iEns=1,nEns + Ec = Ec + wEns(iEns)*EceLDA(iEns) + enddo + +end subroutine RMFL20_lda_correlation_individual_energy diff --git a/src/eDFT/RVWN5_lda_correlation_Levy_Zahariev_shift.f90 b/src/eDFT/RVWN5_lda_correlation_Levy_Zahariev_shift.f90 new file mode 100644 index 0000000..da04ad4 --- /dev/null +++ b/src/eDFT/RVWN5_lda_correlation_Levy_Zahariev_shift.f90 @@ -0,0 +1,74 @@ +subroutine RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ) + +! Compute the restricted VNW5 LDA correlation contribution to Levy-Zahariev shift + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r,rs,x + double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p + double precision :: dxdrs,dxdx_p,decdx_p + double precision :: drsdra,decdra_p + + double precision :: ec_p + +! Output variables + + double precision,intent(out) :: EcLZ + +! Parameters of the functional + + a_p = +0.0621814D0/2D0 + x0_p = -0.10498d0 + b_p = +3.72744d0 + c_p = +12.9352d0 + +! Initialization + + EcLZ = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + rs = (4d0*pi*r/3d0)**(-1d0/3d0) + x = sqrt(rs) + + x_p = x*x + b_p*x + c_p + + xx0_p = x0_p*x0_p + b_p*x0_p + c_p + + q_p = sqrt(4d0*c_p - b_p*b_p) + + ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & + - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) + + drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) + dxdrs = 0.5d0/sqrt(rs) + + dxdx_p = 2d0*x + b_p + + decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & + - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) + + decdra_p = drsdra*dxdrs*decdx_p + + EcLZ = EcLZ + weight(iG)*decdra_p*r*r + + end if + + end do + +end subroutine RVWN5_lda_correlation_Levy_Zahariev_shift diff --git a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 new file mode 100644 index 0000000..df3ef33 --- /dev/null +++ b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 @@ -0,0 +1,75 @@ +subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) + +! Compute the restricted VWN5 LDA correlation individual energies + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r,rI,rs,x + double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p + double precision :: dxdrs,dxdx_p,decdx_p + double precision :: drsdra,decdra_p + double precision :: ec_p + +! Output variables + + double precision :: Ec + +! Parameters of the functional + + a_p = +0.0621814D0/2D0 + x0_p = -0.10498d0 + b_p = +3.72744d0 + c_p = +12.9352d0 + +! Initialization + + Ec = 0d0 + + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + +! spin-up contribution + + if(r > threshold .and. rI > threshold) then + + rs = (4d0*pi*r/3d0)**(-1d0/3d0) + x = sqrt(rs) + + x_p = x*x + b_p*x + c_p + xx0_p = x0_p*x0_p + b_p*x0_p + c_p + q_p = sqrt(4d0*c_p - b_p*b_p) + + ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & + - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) + + drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) + dxdrs = 0.5d0/sqrt(rs) + + dxdx_p = 2d0*x + b_p + + decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & + - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) + + decdra_p = drsdra*dxdrs*decdx_p + + Ec = Ec + weight(iG)*(ec_p + decdra_p*r)*rI + + end if + + end do + +end subroutine RVWN5_lda_correlation_individual_energy diff --git a/src/eDFT/restricted_correlation_Levy_Zahariev_shift.f90 b/src/eDFT/restricted_correlation_Levy_Zahariev_shift.f90 new file mode 100644 index 0000000..695ab2d --- /dev/null +++ b/src/eDFT/restricted_correlation_Levy_Zahariev_shift.f90 @@ -0,0 +1,69 @@ +subroutine restricted_correlation_Levy_Zahariev_shift(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,EcLZ) + +! Compute the correlation part of the Levy-Zahariev shift + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision :: EcLZLDA,EcLZGGA + double precision :: aC + +! Output variables + + double precision,intent(out) :: EcLZ + + select case (rung) + +! Hartree calculation + + case(0) + + EcLZ = 0d0 + +! LDA functionals + + case(1) + + call restricted_lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),EcLZ) + +! GGA functionals + + case(2) + + call print_warning('!!! Individual energies NYI for GGAs !!!') + stop + +! Hybrid functionals + + case(4) + + call print_warning('!!! Individual energies NYI for hybrids !!!') + stop + + aC = 0.81d0 + + EcLZ = EcLZLDA + aC*(EcLZGGA - EcLZLDA) + +! Hartree-Fock calculation + + case(666) + + EcLZ = 0d0 + + end select + +end subroutine restricted_correlation_Levy_Zahariev_shift diff --git a/src/eDFT/restricted_correlation_derivative_discontinuity.f90 b/src/eDFT/restricted_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..bb10b3c --- /dev/null +++ b/src/eDFT/restricted_correlation_derivative_discontinuity.f90 @@ -0,0 +1,65 @@ +subroutine restricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) + +! Compute the correlation part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + +! Local variables + + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec(nEns) + + select case (rung) + +! Hartree calculation + + case(0) + + Ec(:) = 0d0 + +! LDA functionals + + case(1) + + call restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),Ec(:)) + +! GGA functionals + + case(2) + + call print_warning('!!! derivative discontinuity NYI for GGAs !!!') + stop + +! Hybrid functionals + + case(4) + + call print_warning('!!! derivative discontinuity NYI for hybrids !!!') + stop + + aC = 0.81d0 + +! Hartree-Fock calculation + + case(666) + + Ec(:) = 0d0 + + end select + +end subroutine restricted_correlation_derivative_discontinuity diff --git a/src/eDFT/restricted_correlation_individual_energy.f90 b/src/eDFT/restricted_correlation_individual_energy.f90 new file mode 100644 index 0000000..e6cfdb3 --- /dev/null +++ b/src/eDFT/restricted_correlation_individual_energy.f90 @@ -0,0 +1,69 @@ +subroutine restricted_correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) + +! Compute the correlation energy of individual states + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision :: EcLDA + double precision :: EcGGA + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec + + select case (rung) + +! Hartree calculation + + case(0) + + Ec = 0d0 + +! LDA functionals + + case(1) + + call restricted_lda_correlation_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ec) + +! GGA functionals + + case(2) + + call print_warning('!!! Individual energies NYI for GGAs !!!') + stop + +! Hybrid functionals + + case(4) + + call print_warning('!!! Individual energies NYI for hybrids !!!') + stop + + aC = 0.81d0 + +! Hartree-Fock calculation + + case(666) + + Ec = 0d0 + + end select + +end subroutine restricted_correlation_individual_energy diff --git a/src/eDFT/restricted_elda_correlation_Levy_Zahariev_shift.f90 b/src/eDFT/restricted_elda_correlation_Levy_Zahariev_shift.f90 new file mode 100644 index 0000000..c572a0c --- /dev/null +++ b/src/eDFT/restricted_elda_correlation_Levy_Zahariev_shift.f90 @@ -0,0 +1,48 @@ +subroutine restricted_elda_correlation_Levy_Zahariev_shift(nEns,aLF,nGrid,weight,rho,EcLZ) + +! Compute the restricted Levy-Zahariev LDA correlation shift of 2-glomium for various states + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: aLF(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r,ec_p + double precision :: dFcdr + +! Output variables + + double precision,intent(out) :: EcLZ + +! Compute Levy-Zahariev eLDA correlation shift + + EcLZ = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) + dFcdr = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0) + dFcdr = dFcdr/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) + dFcdr = ec_p*dFcdr/6d0 + dFcdr = ec_p + dFcdr + + EcLZ = EcLZ - weight(iG)*r*r*dFcdr + + end if + + end do + +end subroutine restricted_elda_correlation_Levy_Zahariev_shift diff --git a/src/eDFT/restricted_elda_correlation_energy.f90 b/src/eDFT/restricted_elda_correlation_energy.f90 new file mode 100644 index 0000000..8986566 --- /dev/null +++ b/src/eDFT/restricted_elda_correlation_energy.f90 @@ -0,0 +1,44 @@ +subroutine restricted_elda_correlation_energy(nEns,aLF,nGrid,weight,rho,Ec) + +! Compute the restricted LDA correlation energy of 2-glomium for various states + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: aLF(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r,ec_p + +! Output variables + + double precision,intent(out) :: Ec + + +! Compute eLDA correlation energy + + Ec = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) + + Ec = Ec + weight(iG)*ec_p*r + + end if + + end do + +end subroutine restricted_elda_correlation_energy diff --git a/src/eDFT/restricted_elda_correlation_individual_energy.f90 b/src/eDFT/restricted_elda_correlation_individual_energy.f90 new file mode 100644 index 0000000..dd61945 --- /dev/null +++ b/src/eDFT/restricted_elda_correlation_individual_energy.f90 @@ -0,0 +1,52 @@ +subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec) + +! Compute LDA correlation individual energy of 2-glomium for various states + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: aLF(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r + double precision :: rI + double precision :: ec_p,dFcdr + +! Output variables + + double precision,intent(out) :: Ec + +! Compute eLDA correlation potential + + Ec = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .and. rI > threshold) then + + ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) + + dFcdr = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0) + dFcdr = dFcdr/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) + dFcdr = ec_p*dFcdr/(6d0*r) + dFcdr = ec_p + dFcdr*r + + Ec = Ec + weight(iG)*rI*dFcdr + + end if + + end do + +end subroutine restricted_elda_correlation_individual_energy diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index f351af0..2a57ef1 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -93,13 +93,15 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri end do + Ex(:) = 0.5d0*Ex(:) + !------------------------------------------------------------------------ ! Correlation energy !------------------------------------------------------------------------ do iEns=1,nEns - call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & + call restricted_correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & rho(:,iEns),drho(:,:,iEns),Ec(iEns)) end do @@ -108,14 +110,14 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Compute Levy-Zahariev shift !------------------------------------------------------------------------ - call correlation_Levy_Zahariev_shift(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:), & + call restricted_correlation_Levy_Zahariev_shift(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:), & ExLZ,EcLZ,ExcLZ) !------------------------------------------------------------------------ ! Compute derivative discontinuities !------------------------------------------------------------------------ - call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & + call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & ExDD(:),EcDD(:),ExcDD(:)) !------------------------------------------------------------------------ diff --git a/src/eDFT/restricted_lda_correlation_Levy_Zahariev_shift.f90 b/src/eDFT/restricted_lda_correlation_Levy_Zahariev_shift.f90 new file mode 100644 index 0000000..d666ec0 --- /dev/null +++ b/src/eDFT/restricted_lda_correlation_Levy_Zahariev_shift.f90 @@ -0,0 +1,45 @@ +subroutine restricted_lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns,nGrid,weight,rho,EcLZ) + +! Compute the lda correlation part of the Levy-Zahariev shift + + implicit none + + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision,intent(out) :: EcLZ + +! Select correlation functional + + select case (DFA) + +! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 + + case ('RVWN5') + + call RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:),EcLZ) + +! Marut-Fromager-Loos weight-dependent correlation functional + + case ('RMFL20') + + call RMFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight(:),rho(:),EcLZ) + + case default + + call print_warning('!!! LDA correlation functional not available !!!') + stop + + end select + +end subroutine restricted_lda_correlation_Levy_Zahariev_shift diff --git a/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..d6531f3 --- /dev/null +++ b/src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90 @@ -0,0 +1,54 @@ +subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + +! Compute the correlation LDA part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + +! Local variables + + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec(nEns) + +! Select correlation functional + + select case (DFA) + +! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678 + + case ('W38') + + Ec(:) = 0d0 + +! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 + + case ('VWN5') + + Ec(:) = 0d0 + +! Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation) + + case ('LF19') + + call RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),Ec(:)) + + case default + + call print_warning('!!! LDA correlation functional not available !!!') + stop + + end select + +end subroutine restricted_lda_correlation_derivative_discontinuity diff --git a/src/eDFT/restricted_lda_correlation_individual_energy.f90 b/src/eDFT/restricted_lda_correlation_individual_energy.f90 new file mode 100644 index 0000000..a72cad8 --- /dev/null +++ b/src/eDFT/restricted_lda_correlation_individual_energy.f90 @@ -0,0 +1,45 @@ +subroutine restricted_lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ec) + +! Compute LDA correlation energy for individual states + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision :: Ec + +! Select correlation functional + + select case (DFA) + +! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 + + case ('RVWN5') + + call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ec) + +! marut-Fromager-Loos weight-dependent correlation functional + + case ('RMFL20') + + call RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ec) + + case default + + call print_warning('!!! LDA correlation functional not available !!!') + stop + + end select + +end subroutine restricted_lda_correlation_individual_energy diff --git a/src/eDFT/select_rung.f90 b/src/eDFT/select_rung.f90 index 8ae7343..2160590 100644 --- a/src/eDFT/select_rung.f90 +++ b/src/eDFT/select_rung.f90 @@ -34,7 +34,7 @@ subroutine select_rung(rung,DFA) ! Hartree-Fock calculation case(666) - write(*,*) "* rung 666: Hartree-Fock calculation *" + write(*,*) "* rung 666: Hartree-Fock calculation *" ! Default case default