From 310f4ab132c55b5be343a73c198e13039d17faa1 Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Thu, 6 Jan 2022 10:23:48 +0100
Subject: [PATCH] problem print HOMO/LUMO

---
 input/methods          |  8 ++---
 src/eDFT/eDFT_UKS.f90  |  2 +-
 src/eDFT/print_UKS.f90 | 75 +++++++++++++++++++++++++-----------------
 3 files changed, 50 insertions(+), 35 deletions(-)

diff --git a/input/methods b/input/methods
index dcb210a..8124aba 100644
--- a/input/methods
+++ b/input/methods
@@ -1,5 +1,5 @@
 # RHF UHF KS MOM 
-   T   F   F  F   
+   F   F   T  F   
 # MP2* MP3 MP2-F12
   F   F   F
 # CCD pCCD DCD CCSD CCSD(T) 
@@ -11,11 +11,11 @@
 # RPA* RPAx* crRPA ppRPA 
   F    F     F      F
 # G0F2* evGF2* qsGF2* G0F3 evGF3
-  T     F       F      F    F
+  F     F       F      F    F
 # G0W0* evGW* qsGW* ufG0W0 ufGW
-  T    F    F       F      F
+  F    F    F       F      F
 # G0T0 evGT qsGT 
-  T    F    F
+  F    F    F
 # MCMP2
   F
 # * unrestricted version available
diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90
index 154c3c0..f4a20e8 100644
--- a/src/eDFT/eDFT_UKS.f90
+++ b/src/eDFT/eDFT_UKS.f90
@@ -367,7 +367,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
 ! Compute final KS energy
 
   call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_UKS(nBas,nEns,nO,S,wEns,eKS,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
+  call print_UKS(nBas,nEns,occnum,S,wEns,eKS,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
 
 ! Compute Vxc for post-HF calculations
 
diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90
index 2615291..76b68ff 100644
--- a/src/eDFT/print_UKS.f90
+++ b/src/eDFT/print_UKS.f90
@@ -1,4 +1,4 @@
-subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
+subroutine print_UKS(nBas,nEns,occnum,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
 
 ! Print one- and two-electron energies and other stuff for KS calculation
 
@@ -7,39 +7,54 @@ subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
 
 ! Input variables
 
-  integer,intent(in)                 :: nBas
-  integer,intent(in)                 :: nEns
-  integer,intent(in)                 :: nO(nspin)
-  double precision,intent(in)        :: Ov(nBas,nBas)
-  double precision,intent(in)        :: wEns(nEns)
-  double precision,intent(in)        :: eps(nBas,nspin)
-  double precision,intent(in)        :: c(nBas,nBas,nspin)
-  double precision,intent(in)        :: ENuc
-  double precision,intent(in)        :: ET(nspin)
-  double precision,intent(in)        :: EV(nspin)
-  double precision,intent(in)        :: EH(nsp)
-  double precision,intent(in)        :: Ex(nspin)
-  double precision,intent(in)        :: Ec(nsp)
-  double precision,intent(in)        :: Ew
-  double precision,intent(in)        :: dipole(ncart)
+  integer,intent(in)           :: nBas
+  integer,intent(in)           :: nEns
+  double precision,intent(in)  :: occnum(nBas,nspin,nEns)
+  double precision,intent(in)  :: Ov(nBas,nBas)
+  double precision,intent(in)  :: wEns(nEns)
+  double precision,intent(in)  :: eps(nBas,nspin)
+  double precision,intent(in)  :: c(nBas,nBas,nspin)
+  double precision,intent(in)  :: ENuc
+  double precision,intent(in)  :: ET(nspin)
+  double precision,intent(in)  :: EV(nspin)
+  double precision,intent(in)  :: EH(nsp)
+  double precision,intent(in)  :: Ex(nspin)
+  double precision,intent(in)  :: Ec(nsp)
+  double precision,intent(in)  :: Ew
+  double precision,intent(in)  :: dipole(ncart)
 
 ! Local variables
 
-  integer                            :: ixyz
-  integer                            :: ispin
-  integer                            :: iEns
-  integer                            :: iBas
-  integer                            :: HOMO(nspin)
-  integer                            :: LUMO(nspin)
-  double precision                   :: Gap(nspin)
-  double precision                   :: S_exact,S2_exact
-  double precision                   :: S,S2
+  integer                      :: ixyz
+  integer                      :: ispin
+  integer                      :: iEns
+  integer                      :: iBas
+  integer                      :: HOMO(nspin)
+  integer                      :: LUMO(nspin)
+  double precision             :: Gap(nspin)
+  double precision             :: S_exact,S2_exact
+  double precision             :: S,S2
+
+  double precision             :: nO(nspin)
+
+! Compute the number of spin-up and spin-down electrons
+
+  nO(:) = 0d0
+  do ispin=1,nspin
+    do iEns=1,nEns
+      do iBas=1,nBas
+        nO(ispin) = nO(ispin) + wEns(iEns)*occnum(iBas,ispin,iEns) 
+      end do
+    end do
+  end do
+
+  print*,'nO = ',nO
 
 ! HOMO and LUMO
 
   do ispin=1,nspin
 
-      HOMO(ispin) = nO(ispin)
+      HOMO(ispin) = ceiling(nO(ispin))
       LUMO(ispin) = HOMO(ispin) + 1
       Gap(ispin)  = eps(LUMO(ispin),ispin) - eps(HOMO(ispin),ispin)
 
@@ -47,11 +62,11 @@ subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
 
 ! Spin comtamination
 
-  S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
-  S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
+! S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
+! S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
 
-  S_exact = 0.5d0*dble(nO(1) - nO(2))
-  S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
+! S_exact = 0.5d0*dble(nO(1) - nO(2))
+! S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
 
 ! Dump results