diff --git a/include/parameters.h b/include/parameters.h index e5bd88d..279fe1e 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -18,6 +18,5 @@ double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) -! double precision,parameter :: Cx1 = - 0.904d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0) double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) diff --git a/input/dft b/input/dft index 35b19de..4a629a5 100644 --- a/input/dft +++ b/input/dft @@ -6,19 +6,19 @@ # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 1 RMFL20 + 666 HF # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 1 RMFL20 + 666 HF # quadrature grid SG-n 1 # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.0000 0.00000 + 0.0000 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index 0c70040..ca5dc13 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -337,7 +337,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS call 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(:,:), & + eps(:),Pw(:,:),rhow(:),drhow(:,:),J(:,:),Fx(:,:),FxHF(:,:), & Fc(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:)) end subroutine GOK_RKS diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 index 847aa17..9f2657a 100644 --- a/src/eDFT/print_restricted_individual_energy.f90 +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -1,5 +1,5 @@ -subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & - Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD) +subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & + Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD) ! Print individual energies for eDFT calculation @@ -16,8 +16,10 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec, 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) :: Eaux(nEns) double precision,intent(in) :: ExDD(nEns), EcDD(nEns), ExcDD(nEns) double precision,intent(in) :: Omx(nEns), Omc(nEns), Omxc(nEns) + double precision,intent(in) :: Omaux(nEns) double precision,intent(in) :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) double precision,intent(in) :: E(nEns) double precision,intent(in) :: Om(nEns) @@ -103,6 +105,19 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec, write(*,'(A60)') '-------------------------------------------------' write(*,*) +!------------------------------------------------------------------------ +! Auxiliary energies +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' AUXILIARY KS ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') 'Auxiliary KS energy state ',iEns,': ',Eaux(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + !------------------------------------------------------------------------ ! Compute derivative discontinuities !------------------------------------------------------------------------ @@ -124,7 +139,32 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec, !------------------------------------------------------------------------ write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL AND EXCITATION ENERGIES' + write(*,'(A60)') ' EXCITATION ENERGIES FROM AUXILIARY ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Auxiliary excitation energy 1 ->',iEns,': ',Omaux(iEns)+OmxcDD(iEns),' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Auxiliary excitation energy 1 ->',iEns,': ',(Om(iEns)+OmxcDD(iEns))*HaToeV,' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' EXCITATION ENERGIES FROM INDIVIDUAL ENERGIES ' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' diff --git a/src/eDFT/restricted_auxiliary_energy.f90 b/src/eDFT/restricted_auxiliary_energy.f90 new file mode 100644 index 0000000..7e3f251 --- /dev/null +++ b/src/eDFT/restricted_auxiliary_energy.f90 @@ -0,0 +1,42 @@ +subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) + +! Compute the auxiliary KS energies + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nEns + integer,intent(in) :: nO + double precision,intent(in) :: eps(nBas) + +! Local variables + + integer :: iEns + +! Output variables + + double precision,intent(out) :: Eaux(nEns) + +! Ground state density matrix + + iEns = 1 + + Eaux(iEns) = 2d0*sum(eps(1:nO)) + +! Doubly-excited state density matrix + + iEns = 2 + + if(nO > 1) then + Eaux(iEns) = 2d0*sum(eps(1:nO-1)) + else + Eaux(iEns) = 0d0 + end if + + Eaux(iEns) = Eaux(iEns) + 2d0*eps(nO+1) + +end subroutine restricted_auxiliary_energy diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 0f97e30..0191299 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -1,5 +1,5 @@ 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,Ew,EwGIC,E,Om) + nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,EwGIC,E,Om) ! Compute individual energies as well as excitation energies @@ -18,12 +18,14 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri double precision,intent(in) :: AO(nBas,nGrid) double precision,intent(in) :: dAO(ncart,nBas,nGrid) - integer,intent(in) :: nO,nV + integer,intent(in) :: nO + integer,intent(in) :: 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) :: eps(nBas) double precision,intent(in) :: Pw(nBas,nBas) double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) @@ -45,8 +47,10 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri double precision :: EV(nEns) double precision :: EJ(nEns) double precision :: Ex(nEns), Ec(nEns), Exc(nEns) + double precision :: Eaux(nEns) double precision :: ExDD(nEns), EcDD(nEns), ExcDD(nEns) double precision :: Omx(nEns), Omc(nEns), Omxc(nEns) + double precision :: Omaux(nEns) double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) double precision,external :: trace_matrix @@ -106,6 +110,15 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri end do + + + +!------------------------------------------------------------------------ +! Compute auxiliary energies +!------------------------------------------------------------------------ + + call restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) + !------------------------------------------------------------------------ ! Compute derivative discontinuities !------------------------------------------------------------------------ @@ -147,6 +160,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri Omc(iEns) = Ec(iEns) - Ec(1) Omxc(iEns) = Exc(iEns) - Exc(1) + Omaux(iEns) = Eaux(iENs) - Eaux(1) + OmxDD(iEns) = ExDD(iEns) - ExDD(1) OmcDD(iEns) = EcDD(iEns) - EcDD(1) OmxcDD(iEns) = ExcDD(iEns) - ExcDD(1) @@ -158,7 +173,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ call print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & - ExDD(:),EcDD(:),ExcDD(:),E(:), & - Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:)) + Eaux(:),ExDD(:),EcDD(:),ExcDD(:),E(:), & + Om(:),Omx(:),Omc(:),Omxc(:),Omaux,OmxDD(:),OmcDD(:),OmxcDD(:)) end subroutine restricted_individual_energy