auxiliary energies

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)
 

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

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

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, & 
+subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & 
-                                              Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD)
+                                              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'

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

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(:),                           & 
+                                          Eaux(:),ExDD(:),EcDD(:),ExcDD(:),E(:),                   & 
-                                          Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:))
+                                          Om(:),Omx(:),Omc(:),Omxc(:),Omaux,OmxDD(:),OmcDD(:),OmxcDD(:))
 
 end subroutine restricted_individual_energy