From de4927aad4bab4162189874843157b210872b887 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 14 Feb 2021 22:52:17 +0100 Subject: [PATCH] hybrids --- input/dft | 6 +- ...d_correlation_derivative_discontinuity.f90 | 6 +- ...cted_exchange_derivative_discontinuity.f90 | 3 +- ...a_correlation_derivative_discontinuity.f90 | 44 +++++++++++ ...d_correlation_derivative_discontinuity.f90 | 48 ++++++++++++ ...unrestricted_hybrid_correlation_energy.f90 | 56 +++++++++++++ ...estricted_hybrid_correlation_potential.f90 | 65 ++++++++++++++++ ...brid_exchange_derivative_discontinuity.f90 | 54 +++++++++++++ .../unrestricted_hybrid_exchange_energy.f90 | 69 ++++++++++++++++ ...unrestricted_hybrid_exchange_potential.f90 | 78 +++++++++++++++++++ ...a_correlation_derivative_discontinuity.f90 | 34 ++++++++ .../unrestricted_mgga_correlation_energy.f90 | 36 +++++++++ 12 files changed, 492 insertions(+), 7 deletions(-) create mode 100644 src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/unrestricted_hybrid_correlation_energy.f90 create mode 100644 src/eDFT/unrestricted_hybrid_correlation_potential.f90 create mode 100644 src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/unrestricted_hybrid_exchange_energy.f90 create mode 100644 src/eDFT/unrestricted_hybrid_exchange_potential.f90 create mode 100644 src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 create mode 100644 src/eDFT/unrestricted_mgga_correlation_energy.f90 diff --git a/input/dft b/input/dft index bfcaa1b..5c2cad4 100644 --- a/input/dft +++ b/input/dft @@ -6,14 +6,14 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4: HF,B3,PBE - 4 B3 + 4 HF # correlation rung: # Hartree = 0: H # LDA = 1: VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: -# Hybrid = 4: HF,B88,PBE - 4 LYP +# Hybrid = 4: HF,LYP,PBE + 4 HF # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 index 7c7c33c..98e8841 100644 --- a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 @@ -42,19 +42,19 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns, case(2) - call print_warning('!!! derivative discontinuity NYI for GGAs !!!') + call unrestricted_gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! MGGA functionals case(3) - call print_warning('!!! derivative discontinuity NYI for MGGAs !!!') + call unrestricted_mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Hybrid functionals case(4) - call print_warning('!!! derivative discontinuity NYI for hybrids !!!') + call unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) end select diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 index 382e0aa..0a5308a 100644 --- a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 @@ -59,7 +59,8 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC case(4) - call print_warning('!!! exchange part of derivative discontinuity NYI for hybrids !!!') + call unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),& + rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) end select diff --git a/src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..6c83d30 --- /dev/null +++ b/src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 @@ -0,0 +1,44 @@ +subroutine unrestricted_gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + +! Compute the correlation GGA 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,nspin) + +! Local variables + + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec(nsp,nEns) + +! Select correlation functional + + select case (DFA) + + case ('LYP') + + Ec(:,:) = 0d0 + + case ('PBE') + + Ec(:,:) = 0d0 + + case default + + call print_warning('!!! GGA correlation functional not available !!!') + stop + + end select + +end subroutine unrestricted_gga_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..ed74ddd --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 @@ -0,0 +1,48 @@ +subroutine unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + +! Compute the correlation hybrid 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,nspin) + +! Local variables + + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec(nsp,nEns) + +! Select correlation functional + + select case (DFA) + + case ('HF') + + Ec(:,:) = 0d0 + + case ('LYP') + + Ec(:,:) = 0d0 + + case ('PBE') + + Ec(:,:) = 0d0 + + case default + + call print_warning('!!! Hybrid correlation functional not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_hybrid_correlation_energy.f90 b/src/eDFT/unrestricted_hybrid_correlation_energy.f90 new file mode 100644 index 0000000..0e030d0 --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_correlation_energy.f90 @@ -0,0 +1,56 @@ +subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + +! Compute the unrestricted version of the correlation energy for hybrid functionals + + 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,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) + +! Local variables + + double precision :: EcLDA(nsp) + double precision :: EcGGA(nsp) + double precision :: aC + +! Output variables + + double precision,intent(out) :: Ec(nsp) + + select case (DFA) + + case('HF') + + Ec(:) = 0d0 + + case('LYP') + + aC = 0.81d0 + + call unrestricted_lda_correlation_energy('VWN5 ',nEns,wEns,nGrid,weight,rho,EcLDA) + call unrestricted_gga_correlation_energy('LYP ',nEns,wEns,nGrid,weight,rho,drho,EcGGA) + + Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) + + case('PBE') + + call unrestricted_gga_correlation_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,EcGGA) + + Ec(:) = EcGGA(:) + + case default + + call print_warning('!!! Hybrid correlation energy not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_correlation_energy diff --git a/src/eDFT/unrestricted_hybrid_correlation_potential.f90 b/src/eDFT/unrestricted_hybrid_correlation_potential.f90 new file mode 100644 index 0000000..41afc8f --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_correlation_potential.f90 @@ -0,0 +1,65 @@ +subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + +! Compute the correlation potential for hybrid functionals + + 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) :: dAO(ncart,nBas,nGrid) + double precision,intent(in) :: rho(nGrid,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) + +! Local variables + + double precision,allocatable :: FcLDA(:,:,:) + double precision,allocatable :: FcGGA(:,:,:) + double precision :: aC + +! Output variables + + double precision,intent(out) :: Fc(nBas,nBas,nspin) + +! Memory allocation + + select case (DFA) + + case('HF') + + Fc(:,:,:) = 0d0 + + case('LYP') + + allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin)) + + aC = 0.81d0 + + call unrestricted_lda_correlation_potential('VWN5 ',nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA) + call unrestricted_gga_correlation_potential('LYP ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) + + Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:)) + + case('PBE') + + allocate(FcGGA(nBas,nBas,nspin)) + + call unrestricted_gga_correlation_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) + + Fc(:,:,:) = FcGGA(:,:,:) + + case default + + call print_warning('!!! Hybrid correlation potential not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_correlation_potential diff --git a/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..7e5bde8 --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 @@ -0,0 +1,54 @@ +subroutine unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,& + Cx_choice,doNcentered,kappa,ExDD) + +! Compute the exchange part of the derivative discontinuity for hybrid functionals + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + integer,intent(in) :: Cx_choice + logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) + +! Local variables + + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Select exchange functional + + select case (DFA) + + case ('HF') + + ExDD(:) = 0d0 + + case ('B3') + + ExDD(:) = 0d0 + + case ('PBE') + + ExDD(:) = 0d0 + + case default + + call print_warning('!!! Hybrid exchange derivative discontinuity not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_hybrid_exchange_energy.f90 b/src/eDFT/unrestricted_hybrid_exchange_energy.f90 new file mode 100644 index 0000000..d0365b4 --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_exchange_energy.f90 @@ -0,0 +1,69 @@ +subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF, & + rho,drho,Ex,Cx_choice) + +! Compute the exchange energy for hybrid functionals + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: FxHF(nBas,nBas) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + integer,intent(in) :: Cx_choice + +! Local variables + + double precision :: ExLDA,ExGGA,ExHF + double precision :: a0,aX + +! Output variables + + double precision,intent(out) :: Ex + + select case (DFA) + + case ('HF') + + call unrestricted_fock_exchange_energy(nBas,P,FxHF,Ex) + + case ('B3') + + a0 = 0.20d0 + aX = 0.72d0 + + call unrestricted_lda_exchange_energy('S51 ',LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,& + rho,ExLDA,Cx_choice) + call unrestricted_gga_exchange_energy('B88 ',nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + + Ex = ExLDA & + + a0*(ExHF - ExLDA) & + + aX*(ExGGA - ExLDA) + + case ('PBE') + + call unrestricted_gga_exchange_energy('PBE ',nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + + Ex = 0.25d0*ExHF + 0.75d0*ExGGA + + case default + + call print_warning('!!! Hybrid exchange energy not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_exchange_energy diff --git a/src/eDFT/unrestricted_hybrid_exchange_potential.f90 b/src/eDFT/unrestricted_hybrid_exchange_potential.f90 new file mode 100644 index 0000000..f7c5d37 --- /dev/null +++ b/src/eDFT/unrestricted_hybrid_exchange_potential.f90 @@ -0,0 +1,78 @@ +subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, & + ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice) + +! Compute the exchange potential for hybrid functionals + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(ncart,nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + integer,intent(in) :: Cx_choice + +! Local variables + + double precision,allocatable :: FxLDA(:,:),FxGGA(:,:) + double precision :: a0,aX + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas) + +! Memory allocation + + select case (DFA) + + case('HF') + + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + Fx(:,:) = FxHF(:,:) + + case('B3') + + allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas)) + + a0 = 0.20d0 + aX = 0.72d0 + + call unrestricted_lda_exchange_potential('S51 ',LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & + nBas,AO,rho,FxLDA,Cx_choice) + call unrestricted_gga_exchange_potential('B88 ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + + Fx(:,:) = FxLDA(:,:) & + + a0*(FxHF(:,:) - FxLDA(:,:)) & + + aX*(FxGGA(:,:) - FxLDA(:,:)) + + case('PBE') + + allocate(FxGGA(nBas,nBas)) + + call unrestricted_gga_exchange_potential('PBE ',nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + + Fx(:,:) = 0.25d0*FxHF(:,:) + 0.75d0*FxGGA(:,:) + + case default + + call print_warning('!!! Hybrid exchange potential not available !!!') + stop + + end select + +end subroutine unrestricted_hybrid_exchange_potential diff --git a/src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 new file mode 100644 index 0000000..4a870fe --- /dev/null +++ b/src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 @@ -0,0 +1,34 @@ +subroutine unrestricted_mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + +! Compute the correlation MGGA 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,nspin) + +! Local variables + +! Output variables + + double precision,intent(out) :: Ec(nsp,nEns) + +! Select correlation functional + + select case (DFA) + + case default + + call print_warning('!!! MGGA correlation functional not available !!!') + stop + + end select + +end subroutine unrestricted_mgga_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_mgga_correlation_energy.f90 b/src/eDFT/unrestricted_mgga_correlation_energy.f90 new file mode 100644 index 0000000..ae797ac --- /dev/null +++ b/src/eDFT/unrestricted_mgga_correlation_energy.f90 @@ -0,0 +1,36 @@ +subroutine unrestricted_mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + +! Compute unrestricted MGGA correlation energy + + 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,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) + +! Local variables + + integer :: iG + double precision :: ra,rb,ga,gb + +! Output variables + + double precision :: Ec(nsp) + + select case (DFA) + + case default + + call print_warning('!!! MGGA correlation energy not available !!!') + stop + + end select + +end subroutine unrestricted_mgga_correlation_energy