From ed2992bc8515bdb4b2281f6b1da2b6c16e26fd8a Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Tue, 23 Mar 2021 22:35:46 +0100
Subject: [PATCH] clean up qsGW print

---
 input/methods           |  4 ++--
 input/options           |  2 +-
 mol/h2.xyz              |  2 +-
 src/MBPT/print_qsGW.f90 | 15 +++++++--------
 src/MBPT/qsGW.f90       | 30 +++++++++++++++++++++++++++++-
 5 files changed, 40 insertions(+), 13 deletions(-)

diff --git a/input/methods b/input/methods
index 098169e..24288c1 100644
--- a/input/methods
+++ b/input/methods
@@ -11,9 +11,9 @@
 # RPA* RPAx* ppRPA 
   F   F    F    
 # G0F2* evGF2* qsGF2* G0F3 evGF3
-  F     F     T      F    F
+  F     F     F      F    F
 # G0W0* evGW* qsGW*
-  F    F    F
+  F    F    T
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2
diff --git a/input/options b/input/options
index d2044fb..9ce0271 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess stability
-       128    0.0000001  T    5      1          1         F         F
+       128    0.0000001  T    5      1          1         T         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
diff --git a/mol/h2.xyz b/mol/h2.xyz
index 7ab70eb..85810e4 100644
--- a/mol/h2.xyz
+++ b/mol/h2.xyz
@@ -1,4 +1,4 @@
 2
 
 H 0.0 0.0 0.0
-H 0.0 0.0 0.741
+H 0.0 0.0 2.0
diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90
index ffcd712..c3995b0 100644
--- a/src/MBPT/print_qsGW.f90
+++ b/src/MBPT/print_qsGW.f90
@@ -1,4 +1,5 @@
-subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcGM,EcRPA,EqsGW,dipole)
+subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, & 
+                      ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
 
 ! Print one-electron energies and other stuff for qsGW
 
@@ -11,6 +12,10 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
   integer,intent(in)                 :: nO
   integer,intent(in)                 :: nSCF
   double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ET
+  double precision,intent(in)        :: EV
+  double precision,intent(in)        :: EJ
+  double precision,intent(in)        :: Ex
   double precision,intent(in)        :: EcGM
   double precision,intent(in)        :: EcRPA
   double precision,intent(in)        :: Conv
@@ -27,7 +32,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
 ! Local variables
 
   integer                            :: x,ixyz,HOMO,LUMO
-  double precision                   :: Gap,ET,EV,EJ,Ex
+  double precision                   :: Gap
   double precision,external          :: trace_matrix
 
 ! Output variables
@@ -42,12 +47,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z
 
 ! Compute energies
 
-  ET = trace_matrix(nBas,matmul(P,T))
-  EV = trace_matrix(nBas,matmul(P,V))
-  EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
-  Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
-  EqsGW = ET + EV + EJ + Ex 
-
 ! Dump results
 
   write(*,*)'-------------------------------------------------------------------------------'
diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90
index 7802f04..e143ebe 100644
--- a/src/MBPT/qsGW.f90
+++ b/src/MBPT/qsGW.f90
@@ -55,6 +55,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   integer                       :: nBasSq
   integer                       :: ispin
   integer                       :: n_diis
+  double precision              :: ET
+  double precision              :: EV
+  double precision              :: EJ
+  double precision              :: Ex
   double precision              :: EqsGW
   double precision              :: EcRPA
   double precision              :: EcBSE(nspin)
@@ -234,10 +238,34 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
     P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
 
+    !------------------------------------------------------------------------
+    !   Compute total energy
+    !------------------------------------------------------------------------
+
+    ! Kinetic energy
+
+    ET = trace_matrix(nBas,matmul(P,T))
+
+    ! Potential energy
+
+    EV = trace_matrix(nBas,matmul(P,V))
+
+    ! Coulomb energy
+
+    EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
+
+    ! Exchange energy
+
+    Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
+
+    ! Total energy
+
+    EqsGW = ET + EV + EJ + Ex 
+
     ! Print results
 
     call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcGM,EcRPA,EqsGW,dipole)
+    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
 
   enddo
 !------------------------------------------------------------------------