From 52fe50f992b281e7736c4facaffceb7eee98e824 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 17 Mar 2020 11:29:29 +0100 Subject: [PATCH] DD eLDA --- ..._lda_exchange_derivative_discontinuity.f90 | 52 +++++++++++++++++ .../RMFL20_lda_exchange_individual_energy.f90 | 57 +++++++++++++++++++ .../RS51_lda_exchange_individual_energy.f90 | 48 ++++++++++++++++ .../restricted_lda_correlation_energy.f90 | 46 +++++++++++++++ .../restricted_lda_correlation_potential.f90 | 48 ++++++++++++++++ 5 files changed, 251 insertions(+) create mode 100644 src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/RMFL20_lda_exchange_individual_energy.f90 create mode 100644 src/eDFT/RS51_lda_exchange_individual_energy.f90 create mode 100644 src/eDFT/restricted_lda_correlation_energy.f90 create mode 100644 src/eDFT/restricted_lda_correlation_potential.f90 diff --git a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..ab6b337 --- /dev/null +++ b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 @@ -0,0 +1,52 @@ +subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD) + +! Compute the restricted version of the eLDA exchange 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 :: iG + double precision :: Cx0 + double precision :: Cx1 + double precision :: rw + double precision :: dExdw + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Weight-dependent Cx coefficient for RMFL20 exchange functional + + Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) + Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) + +! Compute correlation energy for ground, singly-excited and doubly-excited states + + dExdw = 0d0 + + do iG=1,nGrid + + rw = max(0d0,rhow(iG)) + + if(rw > threshold) then + + dExdw = dExdw + weight(iG)*(Cx1 - Cx0)*rw**(4d0/3d0) + + end if + + end do + + ExDD(1) = - wEns(2) *dExdw + ExDD(2) = (1d0 - wEns(2))*dExdw + +end subroutine RMFL20_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..54f663b --- /dev/null +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -0,0 +1,57 @@ +subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex) + +! Compute the restricted version of the Marut-Fromager-Loos 2020 weight-dependent exchange functional + + 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 + + integer :: iG + double precision :: Cx0 + double precision :: Cx1 + double precision :: CxLDA + double precision :: Cxw + double precision :: r,rI + double precision :: e,dedr + +! Output variables + + double precision,intent(out) :: Ex + +! Weight-dependent Cx coefficient for RMFL20 exchange functional + + Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) + Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) + CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + + Cxw = CxLDA + wEns(2)*(Cx1 - Cx0) + +! Compute LDA exchange matrix in the AO basis + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .and. rI > threshold) then + + e = Cxw*r**(1d0/3d0) + dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + + endif + + enddo + +end subroutine RMFL20_lda_exchange_individual_energy diff --git a/src/eDFT/RS51_lda_exchange_individual_energy.f90 b/src/eDFT/RS51_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..4e43078 --- /dev/null +++ b/src/eDFT/RS51_lda_exchange_individual_energy.f90 @@ -0,0 +1,48 @@ +subroutine RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) + +! Compute the restricted version of Slater's LDA exchange individual energy + + 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 :: Cx + double precision :: r,rI + double precision :: e,dedr + +! Output variables + + double precision,intent(out) :: Ex + +! Cx coefficient for Slater LDA exchange + + Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + +! Compute LDA exchange matrix in the AO basis + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .and. rI > threshold) then + + e = Cx*r**(1d0/3d0) + dedr = 1d0/3d0*Cx*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + + endif + + enddo + +end subroutine RS51_lda_exchange_individual_energy diff --git a/src/eDFT/restricted_lda_correlation_energy.f90 b/src/eDFT/restricted_lda_correlation_energy.f90 new file mode 100644 index 0000000..0c7aa9e --- /dev/null +++ b/src/eDFT/restricted_lda_correlation_energy.f90 @@ -0,0 +1,46 @@ +subroutine restricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) + +! Select LDA correlation functional + + 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) :: Ec + +! Select correlation functional + + select case (DFA) + +! Hartree-Fock + + case ('HF') + + Ec = 0d0 + + case ('RVWN5') + + call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),Ec) + + case ('RMFL20') + + call RMFL20_lda_correlation_energy(nEns,wEns(:),nGrid,weight(:),rho(:),Ec) + + case default + + call print_warning('!!! LDA correlation functional not available !!!') + stop + + end select + +end subroutine restricted_lda_correlation_energy diff --git a/src/eDFT/restricted_lda_correlation_potential.f90 b/src/eDFT/restricted_lda_correlation_potential.f90 new file mode 100644 index 0000000..3c3f171 --- /dev/null +++ b/src/eDFT/restricted_lda_correlation_potential.f90 @@ -0,0 +1,48 @@ +subroutine restricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) + +! Select LDA correlation potential + + 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) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision,intent(out) :: Fc(nBas,nBas) + +! Select correlation functional + + select case (DFA) + +! Hartree-Fock + + case ('HF') + + Fc(:,:) = 0d0 + + case ('RVWN5') + + call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) + + case ('RMFL20') + + call RMFL20_lda_correlation_potential(nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) + + case default + + call print_warning('!!! LDA correlation functional not available !!!') + stop + + end select + +end subroutine restricted_lda_correlation_potential