diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 new file mode 100644 index 0000000..eb1657f --- /dev/null +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -0,0 +1,136 @@ +subroutine print_restricted_individual_energy(nEns,ET,EV,EJ,Ex,Ec,Exc,ExLZ,EcLZ,ExcLZ,ExDD,EcDD,ExcDD,E,Om) + +! Print individual energies for eDFT calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: ET(nEns) + double precision,intent(in) :: EV(nEns) + double precision,intent(in) :: EJ(nEns) + double precision,intent(in) :: Ex(nEns),Ec(nEns),Exc(nEns) + double precision,intent(in) :: ExLZ,EcLZ,ExcLZ + double precision,intent(in) :: ExDD(nEns),EcDD(nEns),ExcDD(nEns) + double precision,intent(in) :: E(nEns) + double precision,intent(in) :: Om(nEns) + +! Local variables + + integer :: iEns + +!------------------------------------------------------------------------ +! Kinetic energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual Kinetic energies' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',ET(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Potential energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual Potential energies' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',EV(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Hartree energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual Hartree energies' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',EJ(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Exchange energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual exchange energies' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',Ex(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Correlation energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual correlation energies' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',Ec(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Compute Levy-Zahariev shift +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,2X,2X,F16.10,A3)') ' x Levy-Zahariev shifts: ',ExLZ, ' au' + write(*,'(A40,2X,2X,F16.10,A3)') ' c Levy-Zahariev shifts: ',EcLZ, ' au' + write(*,'(A40,2X,2X,F16.10,A3)') ' xc Levy-Zahariev shifts: ',ExcLZ,' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Compute derivative discontinuities +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Derivative discontinuities (DD) ' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative ',iEns,': ',ExDD(iEns), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative ',iEns,': ',EcDD(iEns), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative ',iEns,': ',ExcDD(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Total and Excitation energies +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A50)') ' Individual and excitation energies ' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +end subroutine print_restricted_individual_energy diff --git a/src/eDFT/restricted_correlation_energy.f90 b/src/eDFT/restricted_correlation_energy.f90 new file mode 100644 index 0000000..9af7789 --- /dev/null +++ b/src/eDFT/restricted_correlation_energy.f90 @@ -0,0 +1,68 @@ +subroutine restricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + +! Compute the correlation energy + + 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 :: 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 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) + +! 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) + +! Hartree-Fock calculation + + case(666) + + Ec = 0d0 + + end select + +end subroutine restricted_correlation_energy diff --git a/src/eDFT/restricted_correlation_potential.f90 b/src/eDFT/restricted_correlation_potential.f90 new file mode 100644 index 0000000..ccac4fc --- /dev/null +++ b/src/eDFT/restricted_correlation_potential.f90 @@ -0,0 +1,75 @@ +subroutine restricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + +! Compute the correlation potential + + 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) + 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) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision,allocatable :: FcLDA(:,:) + double precision,allocatable :: FcGGA(:,:) + double precision :: aC + +! Output variables + + double precision,intent(out) :: Fc(nBas,nBas) + +! Memory allocation + + select case (rung) + +! Hartree calculation + + case(0) + + Fc(:,:) = 0d0 + +! LDA functionals + + case(1) + + call 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(:,:)) + +! 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(:,:)) + +! Hartree-Fock calculation + + case(666) + + Fc(:,:) = 0d0 + + end select + +end subroutine restricted_correlation_potential diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 new file mode 100644 index 0000000..f351af0 --- /dev/null +++ b/src/eDFT/restricted_individual_energy.f90 @@ -0,0 +1,146 @@ +subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO, & + nO,nV,T,V,ERI,ENuc,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,E,Om) + +! Compute individual energies as well as excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: x_rung,c_rung + character(len=12),intent(in) :: x_DFA,c_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) + + integer,intent(in) :: nO,nV + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ENuc + + double precision,intent(in) :: Pw(nBas,nBas) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + + double precision,intent(in) :: P(nBas,nBas,nEns) + double precision,intent(in) :: rho(nGrid,nEns) + double precision,intent(in) :: drho(ncart,nGrid,nEns) + + double precision,intent(in) :: J(nBas,nBas) + double precision,intent(in) :: Fx(nBas,nBas) + double precision,intent(in) :: FxHF(nBas,nBas) + double precision,intent(in) :: Fc(nBas,nBas) + +! Local variables + + double precision :: ET(nEns) + double precision :: EV(nEns) + double precision :: EJ(nEns) + double precision :: Ex(nEns),Ec(nEns),Exc(nEns) + double precision :: ExLZ,EcLZ,ExcLZ + double precision :: ExDD(nEns),EcDD(nEns),ExcDD(nEns) + + double precision,external :: trace_matrix + + integer :: iEns + +! Output variables + + double precision,intent(out) :: E(nEns) + double precision,intent(out) :: Om(nEns) + +!------------------------------------------------------------------------ +! Kinetic energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + ET(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),T(:,:))) + end do + +!------------------------------------------------------------------------ +! Potential energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + EV(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),V(:,:))) + end do + +!------------------------------------------------------------------------ +! Hartree energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + call hartree_coulomb(nBas,P(:,:,iEns),ERI,J(:,:)) + EJ(iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,iEns),J(:,:))) + end do + +!------------------------------------------------------------------------ +! Exchange energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + + call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),ERI(:,:,:,:), & + AO(:,:),dAO(:,:,:),rho(:,iEns),drho(:,:,iEns),Fx(:,:),FxHF(:,:)) + call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),FxHF(:,:), & + rho(:,iEns),drho(:,:,iEns),Ex(iEns)) + + end do + +!------------------------------------------------------------------------ +! Correlation energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + + call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & + rho(:,iEns),drho(:,:,iEns),Ec(iEns)) + + end do + +!------------------------------------------------------------------------ +! Compute Levy-Zahariev shift +!------------------------------------------------------------------------ + + call 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(:,:), & + ExDD(:),EcDD(:),ExcDD(:)) + +!------------------------------------------------------------------------ +! Total energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + Exc(iEns) = Ex(iEns) + Ec(iEns) + E(iEns) = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) & + + Ex(iEns) + Ec(iEns) + ExcLZ + ExcDD(iEns) + end do + +!------------------------------------------------------------------------ +! Excitation energies +!------------------------------------------------------------------------ + + do iEns=1,nEns + Om(iEns) = E(iEns) - E(1) + end do + +!------------------------------------------------------------------------ +! Dump results +!------------------------------------------------------------------------ + + call print_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:),ExLZ,EcLZ,ExcLZ, & + ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:)) + +end subroutine restricted_individual_energy