diff --git a/input/dft b/input/dft index ce4e31b..9537434 100644 --- a/input/dft +++ b/input/dft @@ -3,7 +3,7 @@ # exchange rung: Hartree = 0, LDA = 1 (RS51,S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 1 RS51 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 - 0 HF + 1 RVWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/src/eDFT/exchange_individual_energy.f90 b/src/eDFT/exchange_individual_energy.f90 index f305d03..f67a0e0 100644 --- a/src/eDFT/exchange_individual_energy.f90 +++ b/src/eDFT/exchange_individual_energy.f90 @@ -1,6 +1,7 @@ -subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex) +subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, & + ERI,P,FxHF,rhow,drhow,rho,drho,Ex) -! Compute the exchange energy of individual states +! Compute the exchange individual energy implicit none include 'parameters.h' @@ -13,6 +14,10 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: FxHF(nBas,nBas) double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: rho(nGrid) @@ -21,7 +26,6 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow ! Local variables double precision :: ExLDA - double precision :: ExGGA double precision :: ExHF ! Output variables @@ -40,7 +44,9 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow case(1) - call lda_exchange_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ex) + call lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,ExLDA) + + Ex = ExLDA ! GGA functionals @@ -53,14 +59,16 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow case(4) - call print_warning('!!! Individual energies NYI for hybrids !!!') + call print_warning('!!! Individual energies NYI for Hybrids !!!') stop ! Hartree-Fock calculation case(666) -! call fock_exchange_individual_energy(nEns,wEns(:),nBas,P,FxHF,ExHF) + call fock_exchange_potential(nBas,P,ERI,FxHF) + call fock_exchange_energy(nBas,P,FxHF,ExHF) + Ex = ExHF end select diff --git a/src/eDFT/lda_correlation_energy.f90 b/src/eDFT/lda_correlation_energy.f90 index b82100f..31008fc 100644 --- a/src/eDFT/lda_correlation_energy.f90 +++ b/src/eDFT/lda_correlation_energy.f90 @@ -44,12 +44,6 @@ subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) call C16_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:)) -! Restricted version of the Marut-Fromager-Loos weight-dependent correlation functional - - case ('RMFL20') - - call RMFL20_lda_correlation_energy(nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:)) - case default call print_warning('!!! LDA correlation functional not available !!!') diff --git a/src/eDFT/lda_correlation_potential.f90 b/src/eDFT/lda_correlation_potential.f90 index 6172ce5..b2d8d2a 100644 --- a/src/eDFT/lda_correlation_potential.f90 +++ b/src/eDFT/lda_correlation_potential.f90 @@ -48,12 +48,6 @@ include 'parameters.h' call C16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) -! Restricted version of the Marut-Fromager-Loos weight-dependent correlation functional - - 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 !!!') diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index 8296c7e..c607139 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -29,6 +29,10 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow ExDD(:) = 0d0 + case ('RS51') + + ExDD(:) = 0d0 + case ('RMFL20') ! call MFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 1be7d48..6fbb190 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -23,13 +23,13 @@ subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex select case (DFA) - case ('S51') + case ('RS51') -! call S51_lda_exchange_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ex) + call RS51_lda_exchange_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ex) - case ('MFL20') + case ('RMFL20') -! call MFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) +! call RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) case default diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index 91c5810..19a0a52 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -39,8 +39,8 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th method = 'GOK-RKS' x_rung = 1 c_rung = 1 - x_DFA = 'S51' - c_DFA = 'W38' + x_DFA = 'RS51' + c_DFA = 'RVWN5' SGn = 0 wEns(:) = 0d0 diff --git a/src/eDFT/restricted_correlation_energy.f90 b/src/eDFT/restricted_correlation_energy.f90 index 9af7789..ec8491b 100644 --- a/src/eDFT/restricted_correlation_energy.f90 +++ b/src/eDFT/restricted_correlation_energy.f90 @@ -38,24 +38,20 @@ subroutine restricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drh case(1) - call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),Ec) + call restricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),Ec) ! GGA functionals case(2) - call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),drho(:,:),Ec) + call print_warning('!!! restricted correlation energies NYI for GGAs !!!') + stop -! Hybrid functionals case(4) - aC = 0.81d0 - - call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),EcLDA) - call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),drho(:,:),EcGGA) - - Ec = EcLDA + aC*(EcGGA - EcLDA) + call print_warning('!!! restricted correlation energies NYI for Hybrids !!!') + stop ! Hartree-Fock calculation diff --git a/src/eDFT/restricted_correlation_potential.f90 b/src/eDFT/restricted_correlation_potential.f90 index ccac4fc..a731203 100644 --- a/src/eDFT/restricted_correlation_potential.f90 +++ b/src/eDFT/restricted_correlation_potential.f90 @@ -22,8 +22,6 @@ subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas ! Local variables double precision,allocatable :: FcLDA(:,:) - double precision,allocatable :: FcGGA(:,:) - double precision :: aC ! Output variables @@ -43,26 +41,21 @@ subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas case(1) - call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) + call restricted_lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:)) ! GGA functionals case(2) - call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:),drho(:,:),Fc(:,:)) + call print_warning('!!! restricted correlation potentials NYI for GGAs !!!') + stop ! Hybrid functionals case(4) - allocate(FcLDA(nBas,nBas),FcGGA(nBas,nBas)) - - aC = 0.81d0 - - call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),FcLDA(:,:)) - call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:),drho(:,:),FcGGA(:,:)) - - Fc(:,:) = FcLDA(:,:) + aC*(FcGGA(:,:) - FcLDA(:,:)) + call print_warning('!!! restricted correlation potentials NYI for Hybrids !!!') + stop ! Hartree-Fock calculation diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 7e58b82..d67a2bd 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -75,7 +75,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ do iEns=1,nEns - call hartree_coulomb(nBas,P(:,:,iEns),ERI,J(:,:)) + call hartree_coulomb(nBas,P(:,:,iEns),ERI(:,:,:,:),J(:,:)) EJ(iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,iEns),J(:,:))) end do @@ -85,8 +85,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri do iEns=1,nEns - call exchange_energy_individual_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),FxHF(:,:), & - rho(:,iEns),drho(:,:,iEns),Ex(iEns)) + call exchange_individual_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & + P(:,:,iEns),FxHF(:,:),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) end do @@ -97,7 +97,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri do iEns=1,nEns call restricted_correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & - rho(:,iEns),drho(:,:,iEns),Ec(iEns)) + rho(:,iEns),drho(:,:,iEns),Ec(iEns)) end do @@ -133,7 +133,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Dump results !------------------------------------------------------------------------ - call print_restricted_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:),ExLZ,EcLZ,ExcLZ, & + call print_restricted_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:)) end subroutine restricted_individual_energy