From 4607acc649dfb8767001892c63b8d52ff65e6cc5 Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Wed, 15 Jan 2020 22:29:43 +0100
Subject: [PATCH] clean up

---
 examples/molecule.H2          |   2 +-
 examples/molecule.N2          |   2 +-
 extract.sh                    | 215 ++++++++++++++++++++--------------
 input/methods                 |   4 +-
 input/molecule                |   2 +-
 input/options                 |   6 +-
 scan_H2.sh                    |   6 +-
 scan_N2.sh                    |   6 +-
 src/QuAcK/CIS.f90             |  19 ++-
 src/QuAcK/G0W0.f90            |  47 ++++----
 src/QuAcK/MP2.f90             |  20 ++--
 src/QuAcK/RHF.f90             |   2 +-
 src/QuAcK/RPA.f90             |  16 +--
 src/QuAcK/evGW.f90            |  40 ++++---
 src/QuAcK/huckel_guess.f90    |  37 +-----
 src/QuAcK/linear_response.f90 |   8 +-
 src/QuAcK/mo_guess.f90        |   8 +-
 src/QuAcK/print_G0W0.f90      |  10 +-
 src/QuAcK/print_RHF.f90       |  36 +++---
 src/QuAcK/print_evGW.f90      |  10 +-
 src/QuAcK/print_qsGW.f90      |  13 +-
 src/QuAcK/qsGW.f90            |  40 ++++---
 22 files changed, 287 insertions(+), 262 deletions(-)

diff --git a/examples/molecule.H2 b/examples/molecule.H2
index 81c624a..579cbe6 100644
--- a/examples/molecule.H2
+++ b/examples/molecule.H2
@@ -2,4 +2,4 @@
  2   1     1   0     0
 # Znuc   x            y           z
   H     0.           0.         0.
-  H     0.           0.         1.4
+  H     0.           0.         3
diff --git a/examples/molecule.N2 b/examples/molecule.N2
index 6905a9c..58ce30b 100644
--- a/examples/molecule.N2
+++ b/examples/molecule.N2
@@ -2,4 +2,4 @@
  2   7     7   0     0
 # Znuc   x            y           z
   N     0.           0.         0.
-  N     0.           0.         2.
+  N     0.           0.         3
diff --git a/extract.sh b/extract.sh
index f72e21a..95e4400 100755
--- a/extract.sh
+++ b/extract.sh
@@ -10,123 +10,164 @@ INPUT=$1
 
   echo 
   echo '*** WFT information ***'
+  echo 
   grep "Hartree-Fock energy" $INPUT
-  EHF=`grep "Hartree-Fock energy" $INPUT | cut -f2 -d":"`
+  EHF=`grep "Hartree-Fock energy" $INPUT | cut -f2 -d"="`
   grep "MP2 correlation energy" $INPUT
+  EcMP2=`grep "MP2 correlation energy" $INPUT | cut -f2 -d"="`
   grep "Ec(MP2) =" $INPUT
+  grep "Ec(CCD) =" $INPUT
   grep "Ec(CCSD) =" $INPUT
   grep "Ec(CCSD(T)) =" $INPUT
 
-  echo 
-  echo '*** Gap information: HF, G0F2, GF2, G0W0 & evGW ***'
-  HF=`grep "HF HOMO-LUMO gap    (eV):" $INPUT | cut -f2 -d":"`
-  G0F2=`grep "GF2  HOMO-LUMO gap    (eV):" $INPUT | head -1 | cut -f2 -d":"`
-  GF2=`grep "GF2  HOMO-LUMO gap    (eV):" $INPUT | tail -1 | cut -f2 -d":"`
-  G0W0=`grep "G0W0 HOMO-LUMO gap    (eV):" $INPUT | cut -f2 -d":"`
-  evGW=`grep "evGW HOMO-LUMO gap    (eV):" $INPUT | tail -1 | cut -f2 -d":"`
+#  echo 
+#  echo '*** Gap information: HF, G0F2, GF2, G0W0 & evGW ***'
+#  HF=`grep "HF HOMO-LUMO gap    (eV):" $INPUT | cut -f2 -d":"`
+#  G0F2=`grep "GF2  HOMO-LUMO gap    (eV):" $INPUT | head -1 | cut -f2 -d":"`
+#  GF2=`grep "GF2  HOMO-LUMO gap    (eV):" $INPUT | tail -1 | cut -f2 -d":"`
+#  G0W0=`grep "G0W0 HOMO-LUMO gap    (eV):" $INPUT | cut -f2 -d":"`
+#  evGW=`grep "evGW HOMO-LUMO gap    (eV):" $INPUT | tail -1 | cut -f2 -d":"`
 
-  echo -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW
+#  echo -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW
 
   echo 
-  echo '*** RPA@TDHF information: RPA, RPA1 & RPA3  ***'
-  RPA_TDHF=`grep "RPA@TDHF correlation energy           =" $INPUT| cut -f2 -d"="`
-  RPA1_TDHF=`grep "RPA@TDHF correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
-  RPA3_TDHF=`grep "RPA@TDHF correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  echo '*** RPA information: Tr@RPA (singlet), Tr@RPA (triplet), AC@RPA (singlet), AC@RPA (triplet) ***'
+  echo 
+  Tr_RPA_1=`grep "Tr@RPA  correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_RPA_3=`grep "Tr@RPA  correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  AC_RPA_1=`grep "AC@RPA  correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  AC_RPA_3=`grep "AC@RPA  correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
 
-  echo -e "\t" $RPA_TDHF "\t" $RPA1_TDHF "\t" $RPA3_TDHF
+  echo -e "\t" $Tr_RPA_1 "\t" $Tr_RPA_3 "\t" $AC_RPA_1 "\t" $AC_RPA_3
 
   echo 
-  echo '*** Ec@G0W0 information: RPA, GM, BSE1 & BSE3 ***'
-  RPA_G0W0=`grep "RPA@G0W0 correlation energy =" $INPUT| cut -f2 -d"="`
-  GM_G0W0=`grep "GM@G0W0  correlation energy =" $INPUT| cut -f2 -d"="`
-  BSE1_G0W0=`grep "BSE@G0W0 correlation energy (singlet)" $INPUT| cut -f2 -d"="`
-  BSE3_G0W0=`grep "BSE@G0W0 correlation energy (triplet)" $INPUT| cut -f2 -d"="`
+  echo '*** RPAx information: Tr@RPAx (singlet), Tr@RPAx (triplet), AC@RPAx (singlet), AC@RPAx (triplet) ***'
+  echo 
+  Tr_RPAx_1=`grep "Tr@RPAx correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_RPAx_3=`grep "Tr@RPAx correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  AC_RPAx_1=`grep "AC@RPAx correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  AC_RPAx_3=`grep "AC@RPAx correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
 
-  echo -e "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0
+  echo -e "\t" $Tr_RPAx_1 "\t" $Tr_RPAx_3 "\t" $AC_RPAx_1 "\t" $AC_RPAx_3
 
   echo 
-  echo '*** Ec@evGW information: RPA, GM, BSE1 & BSE3 ***'
-  RPA_evGW=`grep "RPA@evGW correlation energy =" $INPUT | tail -1| cut -f2 -d"="`
-  GM_evGW=`grep "GM@evGW  correlation energy =" $INPUT | tail -1 | cut -f2 -d"="`
-  BSE1_evGW=`grep "BSE@evGW correlation energy (singlet)" $INPUT | cut -f2 -d"="`
-  BSE3_evGW=`grep "BSE@evGW correlation energy (triplet)" $INPUT | cut -f2 -d"="`
+  echo '*** G0W0 information: Tr@RPA (singlet), Tr@RPA (triplet), Tr@BSE (singlet), Tr@BSE (triplet), AC@BSE (singlet), AC@BSE (triplet) ***'
+  echo 
+  Tr_RPA_G0W0_1=`grep "Tr@RPA@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_RPA_G0W0_3=`grep "Tr@RPA@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  Tr_BSE_G0W0_1=`grep "Tr@BSE@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_BSE_G0W0_3=`grep "Tr@BSE@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  AC_BSE_G0W0_1=`grep "AC@BSE@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  AC_BSE_G0W0_3=`grep "AC@BSE@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
 
-  echo -e "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW
+  echo -e "\t" $Tr_RPA_G0W0_1 "\t" $Tr_RPA_G0W0_3 "\t" $Tr_BSE_G0W0_1 "\t" $Tr_BSE_G0W0_3 "\t" $AC_BSE_G0W0_1 "\t" $AC_BSE_G0W0_3
 
   echo 
-  echo '*** Ec@qsGW information: RPA, GM, BSE1 & BSE3 ***'
-  RPA_qsGW=`grep "RPA@qsGW correlation energy =" $INPUT | tail -1| cut -f2 -d"="`
-  GM_qsGW=`grep "GM@qsGW  correlation energy =" $INPUT | tail -1 | cut -f2 -d"="`
-  BSE1_qsGW=`grep "BSE@qsGW correlation energy (singlet)" $INPUT | cut -f2 -d"="`
-  BSE3_qsGW=`grep "BSE@qsGW correlation energy (triplet)" $INPUT | cut -f2 -d"="`
+  echo '*** evGW information: Tr@RPA (singlet), Tr@RPA (triplet), Tr@BSE (singlet), Tr@BSE (triplet), AC@BSE (singlet), AC@BSE (triplet) ***'
+  echo 
+  Tr_RPA_evGW_1=`grep "Tr@RPA@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_RPA_evGW_3=`grep "Tr@RPA@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  Tr_BSE_evGW_1=`grep "Tr@BSE@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  Tr_BSE_evGW_3=`grep "Tr@BSE@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
+  AC_BSE_evGW_1=`grep "AC@BSE@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
+  AC_BSE_evGW_3=`grep "AC@BSE@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
 
-  echo -e "\t" $RPA_qsGW "\t" $GM_qsGW "\t" $BSE1_qsGW "\t" $BSE3_qsGW
+  echo -e "\t" $Tr_RPA_evGW_1 "\t" $Tr_RPA_evGW_3 "\t" $Tr_BSE_evGW_1 "\t" $Tr_BSE_evGW_3 "\t" $AC_BSE_evGW_1 "\t" $AC_BSE_evGW_3
 
   echo 
-  echo '*** CIS and TDHF excitation energy (singlet & triplet) ***'
-  CIS1_1=`grep "|     1 |" $INPUT | head -1 | cut -f3 -d"|"`
-  CIS3_1=`grep "|     1 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
-  TDHF1_1=`grep "|     1 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
-  TDHF3_1=`grep "|     1 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
-  echo -e "\t" $CIS1_1 "\t" $CIS3_1 "\t" $TDHF1_1 "\t" $TDHF3_1
+  echo '*** CIS excitation energy (singlet & triplet) ***'
+  echo 
 
-  CIS1_2=`grep "|     2 |" $INPUT | head -1 | cut -f3 -d"|"`
-  CIS3_2=`grep "|     2 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
-  TDHF1_2=`grep "|     2 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
-  TDHF3_2=`grep "|     2 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
-  echo -e "\t" $CIS1_2 "\t" $CIS3_2 "\t" $TDHF1_2 "\t" $TDHF3_2
+  CIS_1_1=`grep "|     1 |" $INPUT | head -1 | cut -f3 -d"|"`
+  CIS_1_2=`grep "|     2 |" $INPUT | head -1 | cut -f3 -d"|"`
+  CIS_1_3=`grep "|     3 |" $INPUT | head -1 | cut -f3 -d"|"`
+  CIS_1_4=`grep "|     4 |" $INPUT | head -1 | cut -f3 -d"|"`
+  CIS_1_5=`grep "|     5 |" $INPUT | head -1 | cut -f3 -d"|"`
 
-  CIS1_3=`grep "|     3 |" $INPUT | head -1 | cut -f3 -d"|"`
-  CIS3_3=`grep "|     3 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
-  TDHF1_3=`grep "|     3 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
-  TDHF3_3=`grep "|     3 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
-  echo -e "\t" $CIS1_3 "\t" $CIS3_3 "\t" $TDHF1_3 "\t" $TDHF3_3
+  CIS_3_1=`grep "|     1 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
+  CIS_3_2=`grep "|     2 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
+  CIS_3_3=`grep "|     3 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
+  CIS_3_4=`grep "|     4 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
+  CIS_3_5=`grep "|     5 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
 
-  CIS1_4=`grep "|     4 |" $INPUT | head -1 | cut -f3 -d"|"`
-  CIS3_4=`grep "|     4 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
-  TDHF1_4=`grep "|     4 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
-  TDHF3_4=`grep "|     4 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
-  echo -e "\t" $CIS1_4 "\t" $CIS3_4 "\t" $TDHF1_4 "\t" $TDHF3_4
-
-  CIS1_5=`grep "|     5 |" $INPUT | head -1 | cut -f3 -d"|"`
-  CIS3_5=`grep "|     5 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
-  TDHF1_5=`grep "|     5 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
-  TDHF3_5=`grep "|     5 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
-  echo -e "\t" $CIS1_5 "\t" $CIS3_5 "\t" $TDHF1_5 "\t" $TDHF3_5
+  echo -e "\t" $CIS_1_1 "\t" $CIS_3_1 
+  echo -e "\t" $CIS_1_2 "\t" $CIS_3_2
+  echo -e "\t" $CIS_1_3 "\t" $CIS_3_3
+  echo -e "\t" $CIS_1_4 "\t" $CIS_3_4
+  echo -e "\t" $CIS_1_5 "\t" $CIS_3_5
 
   echo 
-  echo '*** BSE@G0W0 and BSE@evGW excitation energy (singlet & triplet) ***'
-  G0W01_1=`grep "|     1 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
-  G0W03_1=`grep "|     1 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
-  evGW1_1=`grep "|     1 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
-  evGW3_1=`grep "|     1 |" $INPUT | tail -1 | cut -f3 -d"|"`
-  echo -e "\t" $G0W01_1 "\t" $G0W03_1 "\t" $evGW1_1 "\t" $evGW3_1
+  echo '*** RPA excitation energy (singlet & triplet) ***'
+  echo 
 
-  G0W01_2=`grep "|     2 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
-  G0W03_2=`grep "|     2 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
-  evGW1_2=`grep "|     2 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
-  evGW3_2=`grep "|     2 |" $INPUT | tail -1 | cut -f3 -d"|"`
-  echo -e "\t" $G0W01_2 "\t" $G0W03_2 "\t" $evGW1_2 "\t" $evGW3_2
+  RPA_1_1=`grep "|     1 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
+  RPA_1_2=`grep "|     2 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
+  RPA_1_3=`grep "|     3 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
+  RPA_1_4=`grep "|     4 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
+  RPA_1_5=`grep "|     5 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
 
-  G0W01_3=`grep "|     3 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
-  G0W03_3=`grep "|     3 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
-  evGW1_3=`grep "|     3 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
-  evGW3_3=`grep "|     3 |" $INPUT | tail -1 | cut -f3 -d"|"`
-  echo -e "\t" $G0W01_3 "\t" $G0W03_3 "\t" $evGW1_3 "\t" $evGW3_3
+  RPA_3_1=`grep "|     1 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
+  RPA_3_2=`grep "|     2 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
+  RPA_3_3=`grep "|     3 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
+  RPA_3_4=`grep "|     4 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
+  RPA_3_5=`grep "|     5 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
 
-  G0W01_4=`grep "|     4 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
-  G0W03_4=`grep "|     4 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
-  evGW1_4=`grep "|     4 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
-  evGW3_4=`grep "|     4 |" $INPUT | tail -1 | cut -f3 -d"|"`
-  echo -e "\t" $G0W01_4 "\t" $G0W03_4 "\t" $evGW1_4 "\t" $evGW3_4
+  echo -e "\t" $RPA_1_1 "\t" $RPA_3_1 
+  echo -e "\t" $RPA_1_2 "\t" $RPA_3_2
+  echo -e "\t" $RPA_1_3 "\t" $RPA_3_3
+  echo -e "\t" $RPA_1_4 "\t" $RPA_3_4
+  echo -e "\t" $RPA_1_5 "\t" $RPA_3_5
 
-  G0W01_5=`grep "|     5 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
-  G0W03_5=`grep "|     5 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
-  evGW1_5=`grep "|     5 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
-  evGW3_5=`grep "|     5 |" $INPUT | tail -1 | cut -f3 -d"|"`
-  echo -e "\t" $G0W01_5 "\t" $G0W03_5 "\t" $evGW1_5 "\t" $evGW3_5
+  echo 
+  echo '*** RPAx excitation energy (singlet & triplet) ***'
+  echo 
 
+  RPAx_1_1=`grep "|     1 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
+  RPAx_1_2=`grep "|     2 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
+  RPAx_1_3=`grep "|     3 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
+  RPAx_1_4=`grep "|     4 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
+  RPAx_1_5=`grep "|     5 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
+
+  RPAx_3_1=`grep "|     1 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
+  RPAx_3_2=`grep "|     2 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
+  RPAx_3_3=`grep "|     3 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
+  RPAx_3_4=`grep "|     4 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
+  RPAx_3_5=`grep "|     5 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
+
+  echo -e "\t" $RPAx_1_1 "\t" $RPAx_3_1 
+  echo -e "\t" $RPAx_1_2 "\t" $RPAx_3_2
+  echo -e "\t" $RPAx_1_3 "\t" $RPAx_3_3
+  echo -e "\t" $RPAx_1_4 "\t" $RPAx_3_4
+  echo -e "\t" $RPAx_1_5 "\t" $RPAx_3_5
+
+  echo 
+  echo '*** BSE@G0W0 excitation energy (singlet & triplet) ***'
+  echo 
+
+  G0W0_1_1=`grep "|     1 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
+  G0W0_1_2=`grep "|     2 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
+  G0W0_1_3=`grep "|     3 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
+  G0W0_1_4=`grep "|     4 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
+  G0W0_1_5=`grep "|     5 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
+
+  G0W0_3_1=`grep "|     1 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
+  G0W0_3_2=`grep "|     2 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
+  G0W0_3_3=`grep "|     3 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
+  G0W0_3_4=`grep "|     4 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
+  G0W0_3_5=`grep "|     5 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
+
+  echo -e "\t" $G0W0_1_1 "\t" $G0W0_3_1 
+  echo -e "\t" $G0W0_1_2 "\t" $G0W0_3_2
+  echo -e "\t" $G0W0_1_3 "\t" $G0W0_3_3
+  echo -e "\t" $G0W0_1_4 "\t" $G0W0_3_4
+  echo -e "\t" $G0W0_1_5 "\t" $G0W0_3_5
+
+  echo 
   echo '*** MATHEMATICA OUTPUT ***'
-  echo -e "\t" $EHF "\t" $CIS1_1 "\t" $CIS1_2 "\t" $CIS1_3 "\t" $CIS1_4 "\t" $CIS1_5 "\t" $CIS3_1 "\t" $CIS3_2 "\t" $CIS3_3 "\t" $CIS3_4 "\t" $CIS3_5 "\t" $RPA_TDHF "\t" $RPA1_TDHF "\t" $RPA3_TDHF "\t" $TDHF1_1 "\t" $TDHF1_2 "\t" $TDHF1_3 "\t" $TDHF1_4 "\t" $TDHF1_5 "\t" $TDHF3_1 "\t" $TDHF3_2 "\t" $TDHF3_3 "\t" $TDHF3_4 "\t" $TDHF3_5 "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0 "\t" $G0W01_1 "\t" $G0W01_2 "\t" $G0W01_3 "\t" $G0W01_4 "\t" $G0W01_5 "\t" $G0W03_1 "\t" $G0W03_2 "\t" $G0W03_3 "\t" $G0W03_4 "\t" $G0W03_5 "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW "\t" $evGW1_1 "\t" $evGW1_2 "\t" $evGW1_3 "\t" $evGW1_4 "\t" $evGW1_5 "\t" $evGW3_1 "\t" $evGW3_2 "\t" $evGW3_3 "\t" $evGW3_4 "\t" $evGW3_5 
-  echo '*** DONE ***'
+  echo 
+  echo -e "\t" $EHF "\t" $EcMP2 "\t" $Tr_RPA_1 "\t" $Tr_RPA_3 "\t" $AC_RPA_1 "\t" $AC_RPA_3 "\t" $Tr_RPAx_1 "\t" $Tr_RPAx_3 "\t" $AC_RPAx_1 "\t" $AC_RPAx_3 "\t" $Tr_RPA_G0W0_1 "\t" $Tr_RPA_G0W0_3 "\t" $Tr_BSE_G0W0_1 "\t" $Tr_BSE_G0W0_3 "\t" $AC_BSE_G0W0_1 "\t" $AC_BSE_G0W0_3 "\t" $CIS_1_1 "\t" $CIS_1_2 "\t" $CIS_1_3 "\t" $CIS_1_4 "\t" $CIS_1_5 "\t" $CIS_3_1 "\t" $CIS_3_2 "\t" $CIS_3_3 "\t" $CIS_3_4 "\t" $CIS_3_5 "\t" $RPA_1_1 "\t" $RPA_1_2 "\t" $RPA_1_3 "\t" $RPA_1_4 "\t" $RPA_1_5 "\t" $RPA_3_1 "\t" $RPA_3_2 "\t" $RPA_3_3 "\t" $RPA_3_4 "\t" $RPA_3_5 "\t" $RPAx_1_1 "\t" $RPAx_1_2 "\t" $RPAx_1_3 "\t" $RPAx_1_4 "\t" $RPAx_1_5 "\t" $RPAx_3_1 "\t" $RPAx_3_2 "\t" $RPAx_3_3 "\t" $RPAx_3_4 "\t" $RPAx_3_5 "\t" $G0W0_1_1 "\t" $G0W0_1_2 "\t" $G0W0_1_3 "\t" $G0W0_1_4 "\t" $G0W0_1_5 "\t" $G0W0_3_1 "\t" $G0W0_3_2 "\t" $G0W0_3_3 "\t" $G0W0_3_4 "\t" $G0W0_3_5
+
+  echo 
+  echo '*** DONE ***'
+  echo 
 
diff --git a/input/methods b/input/methods
index a97c716..5b88684 100644
--- a/input/methods
+++ b/input/methods
@@ -1,11 +1,11 @@
 # RHF UHF MOM
   T   F   F
 # MP2 MP3 MP2-F12
-  F   F   F
+  T   F   F
 # CCD CCSD CCSD(T) ringCCD ladderCCD
   F   F    F       F       F
 # CIS RPA RPAx ppRPA ADC
-  F   T   F    F     F
+  T   T   T    F     F
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW 
diff --git a/input/molecule b/input/molecule
index 6905a9c..58ce30b 100644
--- a/input/molecule
+++ b/input/molecule
@@ -2,4 +2,4 @@
  2   7     7   0     0
 # Znuc   x            y           z
   N     0.           0.         0.
-  N     0.           0.         2.
+  N     0.           0.         3
diff --git a/input/options b/input/options
index bdf63b6..b8b7c10 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # RHF: maxSCF thresh   DIIS n_diis guess_type ortho_type
-       64    0.0000001  T     5      2          1
+       64    0.0000001  T     5      1          1
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -7,10 +7,10 @@
 # CIS/TDHF/BSE: singlet triplet
                 T       T
 # GF: maxSCF thresh  DIIS n_diis renormalization
-       1     0.00001 T    10      3 
+       64    0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta     
       256    0.00001    T  5      F      F     T   F   F   F   F  0.000  
 # ACFDT: AC XBS
-         T  F
+         T  T
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/scan_H2.sh b/scan_H2.sh
index cbc6e29..798b62c 100755
--- a/scan_H2.sh
+++ b/scan_H2.sh
@@ -1,10 +1,10 @@
 #! /bin/bash
 
 MOL="H2"
-BASIS="VDZ"
+BASIS="VTZ"
 R_START=0.5
 R_END=3.0
-DR=0.01
+DR=0.05
 
 for R in $(seq $R_START $DR $R_END)
 do
@@ -14,6 +14,6 @@ do
   echo "  H     0.           0.         0."  >> examples/molecule.$MOL
   echo "  H     0.           0.         $R"  >> examples/molecule.$MOL
   ./GoDuck $MOL $BASIS > ${MOL}_${BASIS}_${R}.out
-  echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -2 | head -1`
+  echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -4 | head -1`
 done
 
diff --git a/scan_N2.sh b/scan_N2.sh
index 86a445e..35f50e2 100755
--- a/scan_N2.sh
+++ b/scan_N2.sh
@@ -3,8 +3,8 @@
 MOL="N2"
 BASIS="VDZ"
 R_START=1.5
-R_END=3.5
-DR=0.1
+R_END=3.0
+DR=0.05
 
 for R in $(seq $R_START $DR $R_END)
 do
@@ -14,6 +14,6 @@ do
   echo "  N     0.           0.         0."  >> examples/molecule.$MOL
   echo "  N     0.           0.         $R"  >> examples/molecule.$MOL
   ./GoDuck $MOL $BASIS > ${MOL}_${BASIS}_${R}.out
-  echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -2 | head -1`
+  echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -4 | head -1`
 done
 
diff --git a/src/QuAcK/CIS.f90 b/src/QuAcK/CIS.f90
index 2026a36..f3bf07f 100644
--- a/src/QuAcK/CIS.f90
+++ b/src/QuAcK/CIS.f90
@@ -8,13 +8,14 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
 
 ! Input variables
 
-  logical,intent(in)            :: singlet_manifold,triplet_manifold
+  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: triplet_manifold
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas),eHF(nBas)
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
-  logical                       :: dRPA
   logical                       :: dump_matrix = .false.
   logical                       :: dump_trans = .false.
   integer                       :: ispin
@@ -33,10 +34,6 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
 
   lambda = 1d0
 
-! Switch on exchange for CIS
-
-  dRPA = .false.
-
 ! Memory allocation
 
   allocate(A(nS,nS),Omega(nS))
@@ -46,7 +43,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   if(singlet_manifold) then
 
     ispin = 1
-    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (singlet state)'
@@ -55,7 +52,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
     endif
 
     call diagonalize_matrix(nS,A,Omega)
-    call print_excitation('CIS  ',ispin,nS,Omega)
+    call print_excitation('CIS   ',ispin,nS,Omega)
  
     if(dump_trans) then
       print*,'Singlet CIS transition vectors'
@@ -68,7 +65,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   if(triplet_manifold) then
 
     ispin = 2
-    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
+    call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (triplet state)'
@@ -77,7 +74,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
     endif
  
     call diagonalize_matrix(nS,A,Omega)
-    call print_excitation('CIS  ',ispin,nS,Omega)
+    call print_excitation('CIS   ',ispin,nS,Omega)
 
     if(dump_trans) then
       print*,'Triplet CIS transition vectors'
diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90
index cf9a4e4..1e86fe7 100644
--- a/src/QuAcK/G0W0.f90
+++ b/src/QuAcK/G0W0.f90
@@ -97,16 +97,30 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
 ! Solve the quasi-particle equation
 
   eGW(:) = eHF(:) + Z(:)*SigC(:)
-
+ 
 ! Dump results
 
-  call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
+! call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
   call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
 
+! Compute the RPA correlation energy
+
+  call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & 
+                       rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy           =',EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
 ! Plot stuff
 
-  call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
- 
+! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
+
 ! Perform BSE calculation
 
   if(BSE) then
@@ -116,19 +130,10 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -153,10 +158,10 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
 
       write(*,*)
       write(*,*)'-------------------------------------------------------------------------------'
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy           =',EcAC(1) + EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy           =',EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)
 
diff --git a/src/QuAcK/MP2.f90 b/src/QuAcK/MP2.f90
index 3e38149..3f48a11 100644
--- a/src/QuAcK/MP2.f90
+++ b/src/QuAcK/MP2.f90
@@ -56,16 +56,16 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
   EcMP2(1) = EcMP2(2) + EcMP2(3)
 
   write(*,*)
-  write(*,'(A32)')           '-----------------------'
-  write(*,'(A32)')           ' MP2 calculation       '
-  write(*,'(A32)')           '-----------------------'
-  write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
-  write(*,'(A32,1X,F16.10)') ' Direct part           ',EcMP2(2)
-  write(*,'(A32,1X,F16.10)') ' Exchange part         ',EcMP2(3)
-  write(*,'(A32)')           '-----------------------'
-  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy',EHF + EcMP2(1)
-  write(*,'(A32,1X,F16.10)') ' MP2 total       energy',ENuc + EHF + EcMP2(1)
-  write(*,'(A32)')           '-----------------------'
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32)')           ' MP2 calculation          '
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2(1)
+  write(*,'(A32,1X,F16.10)') ' Direct part            = ',EcMP2(2)
+  write(*,'(A32,1X,F16.10)') ' Exchange part          = ',EcMP2(3)
+  write(*,'(A32)')           '--------------------------'
+  write(*,'(A32,1X,F16.10)') ' MP2 electronic  energy = ',EHF + EcMP2(1)
+  write(*,'(A32,1X,F16.10)') ' MP2 total       energy = ',ENuc + EHF + EcMP2(1)
+  write(*,'(A32)')           '--------------------------'
   write(*,*)
 
 end subroutine MP2
diff --git a/src/QuAcK/RHF.f90 b/src/QuAcK/RHF.f90
index ef86d79..ec9ad6a 100644
--- a/src/QuAcK/RHF.f90
+++ b/src/QuAcK/RHF.f90
@@ -61,7 +61,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH
 
 ! Guess coefficients and eigenvalues
 
-  call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
+  call mo_guess(nBas,nO,guess_type,S,Hc,ERI,X,cp,F,Fp,e,c,P)
 
 ! ON(:) = 0d0
 ! do i=1,nO
diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90
index c7f256d..da6afe4 100644
--- a/src/QuAcK/RPA.f90
+++ b/src/QuAcK/RPA.f90
@@ -76,10 +76,10 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A40,F15.6)') 'Tr@RPA  correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A40,F15.6)') 'Tr@RPA  correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A40,F15.6)') 'Tr@RPA  correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A40,F15.6)') 'Tr@RPA  total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA  correlation energy (singlet) =',EcRPA(1)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA  correlation energy (triplet) =',EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA  correlation energy           =',EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA  total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
@@ -98,10 +98,10 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A40,F15.6)') 'AC@RPA  correlation energy (singlet) =',EcAC(1)
-  write(*,'(2X,A40,F15.6)') 'AC@RPA  correlation energy (triplet) =',EcAC(2)
-  write(*,'(2X,A40,F15.6)') 'AC@RPA  correlation energy           =',EcAC(1) + EcAC(2)
-  write(*,'(2X,A40,F15.6)') 'AC@RPA  total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
+  write(*,'(2X,A50,F15.6)') 'AC@RPA  correlation energy (singlet) =',EcAC(1)
+  write(*,'(2X,A50,F15.6)') 'AC@RPA  correlation energy (triplet) =',EcAC(2)
+  write(*,'(2X,A50,F15.6)') 'AC@RPA  correlation energy           =',EcAC(1) + EcAC(2)
+  write(*,'(2X,A50,F15.6)') 'AC@RPA  total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90
index 0648467..126f712 100644
--- a/src/QuAcK/evGW.f90
+++ b/src/QuAcK/evGW.f90
@@ -161,7 +161,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
     ! Print results
 
-    call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
+!   call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
     call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
 
     ! Linear mixing or DIIS extrapolation
@@ -196,7 +196,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
 ! Plot stuff
 
-  call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
+! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
 
 ! Did it actually converge?
 
@@ -212,6 +212,17 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
   endif
 
+! Dump the RPA correlation energy
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy           =',EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
 ! Perform BSE calculation
 
   if(BSE) then
@@ -221,19 +232,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -258,10 +260,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
       write(*,*)
       write(*,*)'-------------------------------------------------------------------------------'
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy           =',EcAC(1) + EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy           =',EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)
 
diff --git a/src/QuAcK/huckel_guess.f90 b/src/QuAcK/huckel_guess.f90
index 0bd9508..ed18bd8 100644
--- a/src/QuAcK/huckel_guess.f90
+++ b/src/QuAcK/huckel_guess.f90
@@ -1,4 +1,4 @@
-subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
+subroutine huckel_guess(nBas,nO,S,Hc,ERI,X,cp,F,Fp,e,c)
 
 !  Hickel guess of the molecular orbitals for HF calculation
 
@@ -11,13 +11,11 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
   double precision,intent(in)   :: S(nBas,nBas)
   double precision,intent(in)   :: Hc(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(inout):: J(nBas,nBas)
-  double precision,intent(inout):: K(nBas,nBas)
   double precision,intent(in)   :: X(nBas,nBas)
   double precision,intent(inout):: cp(nBas,nBas)
+  double precision,intent(inout):: F(nBas,nBas)
   double precision,intent(inout):: Fp(nBas,nBas)
   double precision,intent(inout):: e(nBas)
-  double precision,intent(inout):: P(nBas,nBas)
 
 ! Local variables
 
@@ -30,43 +28,20 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
 
   a = 1.75d0
 
-  Fp(:,:) = Hc(:,:)
+  F(:,:) = Hc(:,:)
 
   do mu=1,nBas
     do nu=mu+1,nBas
 
-      Fp(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
-      Fp(nu,mu) = Fp(mu,nu)
+      F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
+      F(nu,mu) = F(mu,nu)
 
     enddo
   enddo
   
-  Fp(:,:) = matmul(transpose(X(:,:)),matmul(Fp(:,:),X(:,:)))
+  Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
   cp(:,:) = Fp(:,:)
   call diagonalize_matrix(nBas,cp,e)
   c(:,:) = matmul(X(:,:),cp(:,:))
-  P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
-
-! call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
-! call exchange_matrix_AO_basis(nBas,P,ERI,K)
-
-! do mu=1,nBas
-
-!   Fp(mu,mu) = Hc(mu,mu) + J(mu,mu) + 0.5d0*K(mu,mu)
-
-!   do nu=mu+1,nBas
-
-!     Fp(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
-!     Fp(nu,mu) = Fp(mu,nu)
-
-!   enddo
-
-! enddo
-
-!   Fp = matmul(transpose(X),matmul(Fp,X))
-
-!   cp(:,:) = Fp(:,:)
-!   call diagonalize_matrix(nBas,cp,e)
-!   c = matmul(X,cp)
 
 end subroutine
diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90
index 989f410..5c865e7 100644
--- a/src/QuAcK/linear_response.f90
+++ b/src/QuAcK/linear_response.f90
@@ -59,6 +59,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
   if(minval(Omega) < 0d0) &
     call print_warning('You may have instabilities in linear response!!')
 
+  do ia=1,nS
+    if(Omega(ia) < 0d0) Omega(ia) = 0d0
+  end do
+
   call ADAt(nS,AmB,sqrt(Omega),AmBSq)
 
   Z = matmul(AmBSq,matmul(ApB,AmBSq))
@@ -74,11 +78,11 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
 
   Omega = sqrt(Omega)
   XpY = matmul(transpose(Z),AmBSq)
-  call DA(nS,1d0/sqrt(abs(Omega)),XpY)
+  call DA(nS,1d0/sqrt(Omega),XpY)
 
   call ADAt(nS,AmB,1d0/sqrt(Omega),AmBSq)
   XmY = matmul(transpose(Z),AmBSq)
-  call DA(nS,sqrt(abs(Omega)),XmY)
+  call DA(nS,sqrt(Omega),XmY)
 
 ! Compute the RPA correlation energy
 
diff --git a/src/QuAcK/mo_guess.f90 b/src/QuAcK/mo_guess.f90
index f7b416b..8c0e143 100644
--- a/src/QuAcK/mo_guess.f90
+++ b/src/QuAcK/mo_guess.f90
@@ -1,4 +1,4 @@
-subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
+subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,X,cp,F,Fp,e,c,P)
 
 !  Guess of the molecular orbitals for HF calculation
 
@@ -12,10 +12,9 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
   double precision,intent(in)   :: S(nBas,nBas)
   double precision,intent(in)   :: Hc(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(inout):: J(nBas,nBas)
-  double precision,intent(inout):: K(nBas,nBas)
   double precision,intent(in)   :: X(nBas,nBas)
   double precision,intent(inout):: cp(nBas,nBas)
+  double precision,intent(inout):: F(nBas,nBas)
   double precision,intent(inout):: Fp(nBas,nBas)
   double precision,intent(inout):: e(nBas)
   double precision,intent(inout):: P(nBas,nBas)
@@ -37,7 +36,7 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
 
   elseif(guess_type == 2) then
 
-    call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
+    call huckel_guess(nBas,nO,S,Hc,ERI,X,cp,F,Fp,e,c)
 
   elseif(guess_type == 3) then
 
@@ -52,5 +51,4 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
 
   P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
 
-
 end subroutine
diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90
index 782643b..95d1051 100644
--- a/src/QuAcK/print_G0W0.f90
+++ b/src/QuAcK/print_G0W0.f90
@@ -40,11 +40,11 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
   write(*,'(2X,A30,F15.6)') 'G0W0 LUMO      energy (eV):',eGW(LUMO)*HaToeV
   write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap    (eV):',Gap*HaToeV
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy       =',ENuc + EHF + EcRPA
-  write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
-  write(*,'(2X,A30,F15.6)') 'GM@G0W0  total energy       =',ENuc + EHF + EcGM
-  write(*,'(2X,A30,F15.6)') 'GM@G0W0  correlation energy =',EcGM
-  write(*,*)'-------------------------------------------------------------------------------'
+! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy       =',ENuc + EHF + EcRPA
+! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA
+! write(*,'(2X,A30,F15.6)') 'GM@G0W0  total energy       =',ENuc + EHF + EcGM
+! write(*,'(2X,A30,F15.6)') 'GM@G0W0  correlation energy =',EcGM
+! write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
 end subroutine print_G0W0
diff --git a/src/QuAcK/print_RHF.f90 b/src/QuAcK/print_RHF.f90
index bd0aa51..0cea950 100644
--- a/src/QuAcK/print_RHF.f90
+++ b/src/QuAcK/print_RHF.f90
@@ -21,25 +21,25 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF)
 
 
   write(*,*)
-  write(*,'(A50)')           '---------------------------------------'
+  write(*,'(A50)')           '-----------------------------------------'
   write(*,'(A32)')           ' Summary              '
-  write(*,'(A50)')           '---------------------------------------'
-  write(*,'(A32,1X,F16.10)') ' One-electron energy: ',ET + EV
-  write(*,'(A32,1X,F16.10)') ' Kinetic      energy: ',ET
-  write(*,'(A32,1X,F16.10)') ' Potential    energy: ',EV
-  write(*,'(A50)')           '---------------------------------------'
-  write(*,'(A32,1X,F16.10)') ' Two-electron energy: ',EJ + EK
-  write(*,'(A32,1X,F16.10)') ' Coulomb      energy: ',EJ
-  write(*,'(A32,1X,F16.10)') ' Exchange     energy: ',EK
-  write(*,'(A50)')           '---------------------------------------'
-  write(*,'(A32,1X,F16.10)') ' Electronic   energy: ',ERHF
-  write(*,'(A32,1X,F16.10)') ' Nuclear   repulsion: ',ENuc
-  write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy: ',ERHF + ENuc
-  write(*,'(A50)')           '---------------------------------------'
-  write(*,'(A36,F13.6)')  ' HF HOMO      energy (eV):',eHF(HOMO)*HaToeV
-  write(*,'(A36,F13.6)')  ' HF LUMO      energy (eV):',eHF(LUMO)*HaToeV
-  write(*,'(A36,F13.6)')  ' HF HOMO-LUMO gap    (eV):',Gap*HaToeV
-  write(*,'(A50)')           '---------------------------------------'
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A32,1X,F16.10)') ' One-electron energy = ',ET + EV
+  write(*,'(A32,1X,F16.10)') ' Kinetic      energy = ',ET
+  write(*,'(A32,1X,F16.10)') ' Potential    energy = ',EV
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A32,1X,F16.10)') ' Two-electron energy = ',EJ + EK
+  write(*,'(A32,1X,F16.10)') ' Coulomb      energy = ',EJ
+  write(*,'(A32,1X,F16.10)') ' Exchange     energy = ',EK
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A32,1X,F16.10)') ' Electronic   energy = ',ERHF
+  write(*,'(A32,1X,F16.10)') ' Nuclear   repulsion = ',ENuc
+  write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy = ',ERHF + ENuc
+  write(*,'(A50)')           '-----------------------------------------'
+  write(*,'(A36,F13.6)')  ' HF HOMO      energy (eV) = ',eHF(HOMO)*HaToeV
+  write(*,'(A36,F13.6)')  ' HF LUMO      energy (eV) = ',eHF(LUMO)*HaToeV
+  write(*,'(A36,F13.6)')  ' HF HOMO-LUMO gap    (eV) = ',Gap*HaToeV
+  write(*,'(A50)')           '-----------------------------------------'
   write(*,*)
 
 ! Print results
diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90
index fe61210..7596655 100644
--- a/src/QuAcK/print_evGW.f90
+++ b/src/QuAcK/print_evGW.f90
@@ -47,11 +47,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
   write(*,'(2X,A30,F15.6)') 'evGW LUMO      energy (eV):',eGW(LUMO)*HaToeV
   write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap    (eV):',Gap*HaToeV
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy       =',ENuc + EHF + EcRPA
-  write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
-  write(*,'(2X,A30,F15.6)') 'GM@evGW  total energy       =',ENuc + EHF + EcGM
-  write(*,'(2X,A30,F15.6)') 'GM@evGW  correlation energy =',EcGM
-  write(*,*)'-------------------------------------------------------------------------------'
+! write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy       =',ENuc + EHF + EcRPA
+! write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA
+! write(*,'(2X,A30,F15.6)') 'GM@evGW  total energy       =',ENuc + EHF + EcGM
+! write(*,'(2X,A30,F15.6)') 'GM@evGW  correlation energy =',EcGM
+! write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 
 end subroutine print_evGW
diff --git a/src/QuAcK/print_qsGW.f90 b/src/QuAcK/print_qsGW.f90
index 22a709a..b61bdd3 100644
--- a/src/QuAcK/print_qsGW.f90
+++ b/src/QuAcK/print_qsGW.f90
@@ -1,4 +1,4 @@
-subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigmaC,Z,EcRPA,EcGM,EqsGW)
+subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW)
 
 
 ! Print one-electron energies and other stuff for qsGW
@@ -13,7 +13,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
   double precision,intent(in)        :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) 
   double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
   double precision,intent(in)        :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
-  double precision,intent(in)        :: Z(nBas),SigmaC(nBas,nBas)
+  double precision,intent(in)        :: Z(nBas),SigC(nBas,nBas)
 
 ! Local variables
 
@@ -35,8 +35,9 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
   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
+! Ec = 0.5d0*trace_matrix(nBas,matmul(P,SigC))
   Ec = 0d0
+  EqsGW = ET + EV + EJ + Ex + Ec
 
 ! Dump results
 
@@ -69,8 +70,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
   write(*,'(2X,A30,F15.6)') 'qsGW GM total    energy   =',EqsGW + ENuc + EcGM
   write(*,'(2X,A30,F15.6)') 'qsGW exchange    energy   =',Ex
   write(*,'(2X,A30,F15.6)') 'qsGW correlation energy   =',Ec
-  write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA
-  write(*,'(2X,A30,F15.6)') 'GM@qsGW  correlation energy =',EcGM
+! write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA
+! write(*,'(2X,A30,F15.6)') 'GM@qsGW  correlation energy =',EcGM
   write(*,*)'-------------------------------------------'
   write(*,*)
 
@@ -94,7 +95,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
     write(*,'(A32,1X,F16.10)') ' Electronic   energy  ',EqsGW
     write(*,'(A32,1X,F16.10)') ' Nuclear   repulsion  ',ENuc
     write(*,'(A32,1X,F16.10)') ' qsGW         energy  ',ENuc + EqsGW
-    write(*,'(A32,1X,F16.10)') ' RPA corr.    energy  ',EcRPA
+!   write(*,'(A32,1X,F16.10)') ' RPA corr.    energy  ',EcRPA
     write(*,'(A50)')           '---------------------------------------'
     write(*,*)
  
diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90
index 74e5e93..d2a10b0 100644
--- a/src/QuAcK/qsGW.f90
+++ b/src/QuAcK/qsGW.f90
@@ -204,7 +204,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
     ! Print results
 
-    call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
+!   call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
     call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW)
 
     ! Increment
@@ -218,7 +218,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
 ! Compute second-order correction of the Hermitization error
 
-  call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
+! call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
 
 ! Compute the overlap between HF and GW orbitals
 
@@ -242,6 +242,17 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
   endif
 
+! Dump RPA correlation energy
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy           =',EcRPA(1) + EcRPA(2)
+  write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW total energy                 =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
 ! Perform BSE calculation
 
   if(BSE) then
@@ -251,19 +262,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy           =',EcRPA(1) + EcRPA(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW total energy                 =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy           =',EcBSE(1) + EcBSE(2)
-    write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW total energy                 =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW total energy                 =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
@@ -288,10 +290,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
 
       write(*,*)
       write(*,*)'-------------------------------------------------------------------------------'
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy           =',EcAC(1) + EcAC(2)
-      write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW total energy                 =',ENuc + EqsGW + EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy           =',EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW total energy                 =',ENuc + EqsGW + EcAC(1) + EcAC(2)
       write(*,*)'-------------------------------------------------------------------------------'
       write(*,*)