fix problem with individual energies

input/dft

@@ -1,25 +1,25 @@
 # Restricted or unrestricted KS calculation
   eDFT-UKS
 # exchange rung:
-# Hartree      = 0
+# Hartree      = 0: H
-# LDA          = 1: RS51,RMFL20
+# LDA          = 1: S51,CC
-# GGA          = 2: RB88
+# GGA          = 2: B88
-# Hybrid       = 4
+# Hybrid       = 4:
-# Hartree-Fock = 666
+# Hartree-Fock = 666: HF
-  666 HF
+  1 S51
 # correlation rung: 
-# Hartree      = 0
+# Hartree      = 0: H
-# LDA          = 1: RVWN5,RMFL20
+# LDA          = 1: VWN5,MFL20
 # GGA          = 2: 
 # Hybrid       = 4: 
-# Hartree-Fock = 666
+# Hartree-Fock = 666: HF
   0 H 
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
   3
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.0 0.0
+  0.000000 0.0
 # Parameters for CC weight-dependent exchange functional
 0.000000 0.0000000 0.000000
 0.000000 0.0000000 0.0000000

input/int

@@ -1,13 +1,13 @@
 # Debuggin mode?
   F
 # Chemist notation for two-electron integral?
-  F
+  T
 # Exposant of the Slater geminal
-  0.5
+  1.0
 # One-electron integrals: Ov Kin Nuc
                           T  T   T
 # Two-electron integrals: ERI F12 Yuk Erf
-                          T   F   F   T
+                          T   F   F   F
 # Three-electron integrals: Type1 Type2 Type3
                             F     F     F
 # Four-electron integrals: Type1 Type2 Type3

input/options

@@ -9,7 +9,7 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.00367493 3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta        COHSEX SOSEX TDA_W G0W GW0 
-         256    0.00001   T    5     T  0.00367493  F      F     T     F   F  
+         256    0.00001   T    5     T  0.00367493  F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn

src/QuAcK/G0W0.f90

@@ -104,6 +104,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
   call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & 
                        rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
+  if(print_W) call print_excitation('RPA@HF      ',ispin,nS,Omega(:,ispin))
+
 ! Compute correlation part of the self-energy 
 
   call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
@@ -143,8 +145,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
 
 ! Dump results
 
-  if(print_W) call print_excitation('RPA@G0W0    ',ispin,nS,Omega(:,ispin))
-
   call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
 
 ! Compute the RPA correlation energy

src/eDFT/UCC_lda_exchange_individual_energy.f90

@@ -86,11 +86,18 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
     r  = max(0d0,rhow(iG))
     rI = max(0d0,rho(iG))
 
-    if(r > threshold .and. rI > threshold) then
+    if(r > threshold) then
 
       e_p  =         Cx*r**(1d0/3d0)
       dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
-      Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
+
+      Ex = Ex - weight(iG)*dedr*r*r
+
+      if(rI > threshold) then
+
+      Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
+
+      endif
 
     endif
 

src/eDFT/UVWN5_lda_correlation_individual_energy.f90

@@ -61,7 +61,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
    
 !   spin-up contribution
    
-    if(ra > threshold .and. raI > threshold) then
+    if(ra > threshold .or. raI > threshold) then
    
       r = ra + rb
       rI = raI + rbI
@@ -131,7 +131,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
 
 !   spin-down contribution
    
-    if(rb > threshold .and. rbI > threshold) then
+    if(rb > threshold .or. rbI > threshold) then
    
       r = ra + rb
       rI = raI + rbI

src/eDFT/print_UKS.f90

@@ -17,18 +17,31 @@ subroutine print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
   double precision,intent(in)        :: Ec(nsp)
   double precision,intent(in)        :: Ew
 
+  integer                            :: ispin
   integer                            :: HOMO(nspin)
   integer                            :: LUMO(nspin)
   double precision                   :: Gap(nspin)
 
 ! HOMO and LUMO
 
-  HOMO(:) = nO(:)
 
-  LUMO(:) = HOMO(:) + 1
+  do ispin=1,nspin
 
-  Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1)
+    if(nO(ispin) > 0) then 
-  Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2)
+
+      HOMO(ispin) = nO(ispin)
+      LUMO(ispin) = HOMO(ispin) + 1
+      Gap(ispin)  = eps(LUMO(ispin),ispin) - eps(HOMO(ispin),ispin)
+
+    else
+
+      HOMO(ispin) = 0
+      LUMO(ispin) = 0
+      Gap(ispin)  = 0d0
+
+    end if
+
+  end do
 
 ! Dump results
 

src/eDFT/print_unrestricted_individual_energy.f90

@@ -137,7 +137,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
 !------------------------------------------------------------------------
 
   write(*,'(A60)') '-------------------------------------------------'
-  write(*,'(A60)') ' IP and EA FROM AUXILIARY ENERGIES '
+  write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES '
   write(*,'(A60)') '-------------------------------------------------'
 
     write(*,'(A43,F16.10,A4)') ' Ionization Potential  1 -> 2:',Omaux(2)+OmxcDD(2),' au'