From 941ab79dc10812eb76e0fc88a36abbfb7737869f Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Sat, 21 Mar 2020 16:31:39 +0100
Subject: [PATCH] numgrid

---
 examples/molecule.CO                         |   2 +-
 examples/molecule.F2                         |   2 +-
 examples/molecule.H2                         |   2 +-
 examples/molecule.LiF                        |   2 +-
 input/basis                                  | 121 ++++++++++----
 input/methods                                |   4 +-
 input/molecule                               |   5 +-
 input/molecule.xyz                           |   5 +-
 input/options                                |   4 +-
 input/weight                                 | 121 ++++++++++----
 scan_H2.sh                                   |   2 +-
 scan_LiF.sh                                  |   8 +-
 scan_LiH.sh                                  |   2 +-
 src/QuAcK/G0F2.f90                           |  31 ++--
 src/QuAcK/G0T0.f90                           |  10 +-
 src/QuAcK/QuAcK.f90                          |   5 +-
 src/QuAcK/evGF2.f90                          |  14 +-
 src/QuAcK/excitation_density_Tmatrix.f90     |  14 +-
 src/QuAcK/linear_response_C_pp.f90           |   1 +
 src/QuAcK/linear_response_D_pp.f90           |   1 +
 src/QuAcK/linear_response_pp.f90             |   7 +-
 src/QuAcK/ppRPA.f90                          |  11 +-
 src/QuAcK/print_G0F2.f90                     |  21 +--
 src/QuAcK/print_G0T0.f90                     |  11 +-
 src/QuAcK/print_evGF2.f90                    |  43 +++--
 src/QuAcK/print_evGF3.f90                    |   6 +-
 src/QuAcK/renormalization_factor_Tmatrix.f90 |   8 +-
 src/QuAcK/self_energy_Tmatrix_diag.f90       |   8 +-
 src/QuAcK/soG0T0.f90                         |  11 +-
 src/QuAcK/sort_ppRPA.f90                     |  31 ++--
 src/eDFT/Makefile                            |  21 +--
 src/eDFT/numgrid.f90                         | 113 ++++++++++++++
 src/eDFT/test_explicit.f90                   | 156 +++++++++++++++++++
 33 files changed, 623 insertions(+), 180 deletions(-)
 create mode 100644 src/eDFT/numgrid.f90
 create mode 100644 src/eDFT/test_explicit.f90

diff --git a/examples/molecule.CO b/examples/molecule.CO
index 4254094..cbb5734 100644
--- a/examples/molecule.CO
+++ b/examples/molecule.CO
@@ -2,4 +2,4 @@
  2   7     7   0     0
 # Znuc   x            y           z
   C     0.           0.         0.
-  O     0.           0.         2.134
+  O     0.           0.         2.8
diff --git a/examples/molecule.F2 b/examples/molecule.F2
index 9c99be9..8f011d0 100644
--- a/examples/molecule.F2
+++ b/examples/molecule.F2
@@ -2,4 +2,4 @@
  2   9     9   0     0
 # Znuc   x            y           z
   F     0.           0.         0.
-  F     0.           0.         2
+  F     0.           0.         2.640
diff --git a/examples/molecule.H2 b/examples/molecule.H2
index 93d91e8..8b4dc85 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.402
+  H     0.           0.         2.3
diff --git a/examples/molecule.LiF b/examples/molecule.LiF
index 648b5da..381ba52 100644
--- a/examples/molecule.LiF
+++ b/examples/molecule.LiF
@@ -2,4 +2,4 @@
  2   6     6   0     0
 # Znuc   x            y           z
   Li    0.           0.         0.
-  F     0.           0.         2.974
+  F     0.           0.         2.965
diff --git a/input/basis b/input/basis
index 6f3d2a9..76e69dd 100644
--- a/input/basis
+++ b/input/basis
@@ -1,39 +1,98 @@
-1  10
-S   8
-  1  24350.0000000              0.0005020        
-  2   3650.0000000              0.0038810        
-  3    829.6000000              0.0199970        
-  4    234.0000000              0.0784180        
-  5     75.6100000              0.2296760        
-  6     26.7300000              0.4327220        
-  7      9.9270000              0.3506420        
-  8      1.1020000             -0.0076450        
-S   8
-  1  24350.0000000             -0.0001180        
-  2   3650.0000000             -0.0009150        
-  3    829.6000000             -0.0047370        
-  4    234.0000000             -0.0192330        
-  5     75.6100000             -0.0603690        
-  6     26.7300000             -0.1425080        
-  7      9.9270000             -0.1777100        
-  8      1.1020000              0.6058360        
+1  15
+S   9
+  1   6601.0000000              0.0001170        
+  2    989.7000000              0.0009110        
+  3    225.7000000              0.0047280        
+  4     64.2900000              0.0191970        
+  5     21.1800000              0.0630470        
+  6      7.7240000              0.1632080        
+  7      3.0030000              0.3148270        
+  8      1.2120000              0.3939360        
+  9      0.4930000              0.1969180        
+S   9
+  1   6601.0000000             -0.0000180        
+  2    989.7000000             -0.0001420        
+  3    225.7000000             -0.0007410        
+  4     64.2900000             -0.0030200        
+  5     21.1800000             -0.0101230        
+  6      7.7240000             -0.0270940        
+  7      3.0030000             -0.0573590        
+  8      1.2120000             -0.0938950        
+  9      0.4930000             -0.1210910        
 S   1
-  1      2.8360000              1.0000000        
+  1      0.0951500              1.0000000        
 S   1
-  1      0.3782000              1.0000000        
+  1      0.0479100              1.0000000        
+S   1
+  1      0.0222000              1.0000000        
 P   3
-  1     54.7000000              0.0171510        
-  2     12.4300000              0.1076560        
-  3      3.6790000              0.3216810        
+  1      6.2500000              0.0033880        
+  2      1.3700000              0.0193160        
+  3      0.3672000              0.0791040        
 P   1
-  1      1.1430000              1.0000000        
+  1      0.1192000              1.0000000        
 P   1
-  1      0.3300000              1.0000000        
+  1      0.0447400              1.0000000        
+P   1
+  1      0.0179500              1.0000000        
 D   1
-  1      4.0140000              1.0000000        
+  1      0.3440000              1.0000000        
 D   1
-  1      1.0960000              1.0000000        
+  1      0.1530000              1.0000000        
+D   1
+  1      0.0680000              1.0000000        
 F   1
-  1      2.5440000              1.0000000
-
-
+  1      0.2460000              1.0000000        
+F   1
+  1      0.1292000              1.0000000        
+G   1
+  1      0.2380000              1.0000000
+2  15
+S   9
+  1  74530.0000000              0.0000950        
+  2  11170.0000000              0.0007380        
+  3   2543.0000000              0.0038580        
+  4    721.0000000              0.0159260        
+  5    235.9000000              0.0542890        
+  6     85.6000000              0.1495130        
+  7     33.5500000              0.3082520        
+  8     13.9300000              0.3948530        
+  9      5.9150000              0.2110310        
+S   9
+  1  74530.0000000             -0.0000220        
+  2  11170.0000000             -0.0001720        
+  3   2543.0000000             -0.0008910        
+  4    721.0000000             -0.0037480        
+  5    235.9000000             -0.0128620        
+  6     85.6000000             -0.0380610        
+  7     33.5500000             -0.0862390        
+  8     13.9300000             -0.1558650        
+  9      5.9150000             -0.1109140        
+S   1
+  1      1.8430000              1.0000000        
+S   1
+  1      0.7124000              1.0000000        
+S   1
+  1      0.2637000              1.0000000        
+P   3
+  1     80.3900000              0.0063470        
+  2     18.6300000              0.0442040        
+  3      5.6940000              0.1685140        
+P   1
+  1      1.9530000              1.0000000        
+P   1
+  1      0.6702000              1.0000000        
+P   1
+  1      0.2166000              1.0000000        
+D   1
+  1      5.0140000              1.0000000        
+D   1
+  1      1.7250000              1.0000000        
+D   1
+  1      0.5860000              1.0000000        
+F   1
+  1      3.5620000              1.0000000        
+F   1
+  1      1.1480000              1.0000000        
+G   1
+  1      2.3760000              1.0000000
diff --git a/input/methods b/input/methods
index 2b2d735..5ee8f11 100644
--- a/input/methods
+++ b/input/methods
@@ -7,10 +7,10 @@
 # CIS RPA RPAx ppRPA ADC
   F   F   F    F     F
 # G0F2 evGF2 G0F3 evGF3
-  T    F     F    F
+  F    F     F    F
 # G0W0 evGW qsGW 
   F    F    F
 # G0T0 evGT qsGT 
-  F    F    F
+  T    F    F
 # MCMP2
   F
diff --git a/input/molecule b/input/molecule
index edeba31..381ba52 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   5     5   0     0
+ 2   6     6   0     0
 # Znuc   x            y           z
-  Ne    0.0          0.0         0.0
+  Li    0.           0.         0.
+  F     0.           0.         2.965
diff --git a/input/molecule.xyz b/input/molecule.xyz
index 1c70680..87babf0 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,3 +1,4 @@
-  1
+  2
 
- Ne     0.0000000000    0.0000000000    0.0000000000
+ Li     0.0000000000    0.0000000000    0.0000000000
+ F      0.0000000000    0.0000000000    1.5690105433
diff --git a/input/options b/input/options
index b94b9a1..2206914 100644
--- a/input/options
+++ b/input/options
@@ -7,10 +7,10 @@
 # CIS/TDHF/BSE: singlet triplet
                 T       T
 # GF: maxSCF thresh  DIIS n_diis lin renorm
-      256    0.00001 T    5      F   3 
+      256    0.00001 T    5      T   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   T    0.000   
 # ACFDT: AC Kx  XBS
-         T  F   F
+         T  F   T
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/input/weight b/input/weight
index 6f3d2a9..76e69dd 100644
--- a/input/weight
+++ b/input/weight
@@ -1,39 +1,98 @@
-1  10
-S   8
-  1  24350.0000000              0.0005020        
-  2   3650.0000000              0.0038810        
-  3    829.6000000              0.0199970        
-  4    234.0000000              0.0784180        
-  5     75.6100000              0.2296760        
-  6     26.7300000              0.4327220        
-  7      9.9270000              0.3506420        
-  8      1.1020000             -0.0076450        
-S   8
-  1  24350.0000000             -0.0001180        
-  2   3650.0000000             -0.0009150        
-  3    829.6000000             -0.0047370        
-  4    234.0000000             -0.0192330        
-  5     75.6100000             -0.0603690        
-  6     26.7300000             -0.1425080        
-  7      9.9270000             -0.1777100        
-  8      1.1020000              0.6058360        
+1  15
+S   9
+  1   6601.0000000              0.0001170        
+  2    989.7000000              0.0009110        
+  3    225.7000000              0.0047280        
+  4     64.2900000              0.0191970        
+  5     21.1800000              0.0630470        
+  6      7.7240000              0.1632080        
+  7      3.0030000              0.3148270        
+  8      1.2120000              0.3939360        
+  9      0.4930000              0.1969180        
+S   9
+  1   6601.0000000             -0.0000180        
+  2    989.7000000             -0.0001420        
+  3    225.7000000             -0.0007410        
+  4     64.2900000             -0.0030200        
+  5     21.1800000             -0.0101230        
+  6      7.7240000             -0.0270940        
+  7      3.0030000             -0.0573590        
+  8      1.2120000             -0.0938950        
+  9      0.4930000             -0.1210910        
 S   1
-  1      2.8360000              1.0000000        
+  1      0.0951500              1.0000000        
 S   1
-  1      0.3782000              1.0000000        
+  1      0.0479100              1.0000000        
+S   1
+  1      0.0222000              1.0000000        
 P   3
-  1     54.7000000              0.0171510        
-  2     12.4300000              0.1076560        
-  3      3.6790000              0.3216810        
+  1      6.2500000              0.0033880        
+  2      1.3700000              0.0193160        
+  3      0.3672000              0.0791040        
 P   1
-  1      1.1430000              1.0000000        
+  1      0.1192000              1.0000000        
 P   1
-  1      0.3300000              1.0000000        
+  1      0.0447400              1.0000000        
+P   1
+  1      0.0179500              1.0000000        
 D   1
-  1      4.0140000              1.0000000        
+  1      0.3440000              1.0000000        
 D   1
-  1      1.0960000              1.0000000        
+  1      0.1530000              1.0000000        
+D   1
+  1      0.0680000              1.0000000        
 F   1
-  1      2.5440000              1.0000000
-
-
+  1      0.2460000              1.0000000        
+F   1
+  1      0.1292000              1.0000000        
+G   1
+  1      0.2380000              1.0000000
+2  15
+S   9
+  1  74530.0000000              0.0000950        
+  2  11170.0000000              0.0007380        
+  3   2543.0000000              0.0038580        
+  4    721.0000000              0.0159260        
+  5    235.9000000              0.0542890        
+  6     85.6000000              0.1495130        
+  7     33.5500000              0.3082520        
+  8     13.9300000              0.3948530        
+  9      5.9150000              0.2110310        
+S   9
+  1  74530.0000000             -0.0000220        
+  2  11170.0000000             -0.0001720        
+  3   2543.0000000             -0.0008910        
+  4    721.0000000             -0.0037480        
+  5    235.9000000             -0.0128620        
+  6     85.6000000             -0.0380610        
+  7     33.5500000             -0.0862390        
+  8     13.9300000             -0.1558650        
+  9      5.9150000             -0.1109140        
+S   1
+  1      1.8430000              1.0000000        
+S   1
+  1      0.7124000              1.0000000        
+S   1
+  1      0.2637000              1.0000000        
+P   3
+  1     80.3900000              0.0063470        
+  2     18.6300000              0.0442040        
+  3      5.6940000              0.1685140        
+P   1
+  1      1.9530000              1.0000000        
+P   1
+  1      0.6702000              1.0000000        
+P   1
+  1      0.2166000              1.0000000        
+D   1
+  1      5.0140000              1.0000000        
+D   1
+  1      1.7250000              1.0000000        
+D   1
+  1      0.5860000              1.0000000        
+F   1
+  1      3.5620000              1.0000000        
+F   1
+  1      1.1480000              1.0000000        
+G   1
+  1      2.3760000              1.0000000
diff --git a/scan_H2.sh b/scan_H2.sh
index 2b7a66a..0d4103a 100755
--- a/scan_H2.sh
+++ b/scan_H2.sh
@@ -1,7 +1,7 @@
 #! /bin/bash
 
 MOL="H2"
-BASIS="cc-pvqz"
+BASIS="cc-pvdz"
 R_START=1.0
 R_END=2.4
 DR=0.1
diff --git a/scan_LiF.sh b/scan_LiF.sh
index c3f11d7..f98b3d4 100755
--- a/scan_LiF.sh
+++ b/scan_LiF.sh
@@ -1,10 +1,10 @@
 #! /bin/bash
 
 MOL="LiF"
-BASIS="cc-pvqz"
-R_START=2.965
-R_END=2.965
-DR=0.001
+BASIS="cc-pvdz"
+R_START=2.5
+R_END=3.6
+DR=0.1
 
 for R in $(seq $R_START $DR $R_END)
 do
diff --git a/scan_LiH.sh b/scan_LiH.sh
index 38759d5..c4229cb 100755
--- a/scan_LiH.sh
+++ b/scan_LiH.sh
@@ -1,7 +1,7 @@
 #! /bin/bash
 
 MOL="LiH"
-BASIS="cc-pvqz"
+BASIS="cc-pvdz"
 R_START=2.5
 R_END=3.6
 DR=0.1
diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90
index bc7511c..4f03cb6 100644
--- a/src/QuAcK/G0F2.f90
+++ b/src/QuAcK/G0F2.f90
@@ -21,7 +21,7 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
   double precision              :: eps
   double precision              :: VV
   double precision,allocatable  :: eGF2(:)
-  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Sig(:)
   double precision,allocatable  :: Z(:)
 
   integer                       :: i,j,a,b,p
@@ -36,13 +36,19 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
 
 ! Memory allocation
 
-  allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas))
+  allocate(Sig(nBas),Z(nBas),eGF2(nBas))
 
+  if(linearize) then 
+  
+     write(*,*) '*** Quasiparticle equation will be linearized ***'
+     write(*,*)
 
-  ! Frequency-dependent second-order contribution
+  end  if
 
-  Bpp(:,:) = 0d0
-  Z(:)     = 0d0
+! Frequency-dependent second-order contribution
+
+  Sig(:) = 0d0
+  Z(:)   = 0d0
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
@@ -51,8 +57,8 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
 
           eps = e0(p) + e0(a) - e0(i) - e0(j)
           VV  = (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)
-          Bpp(p,1) = Bpp(p,1) + VV/eps
-          Z(p)     = Z(p)     + VV/eps**2
+          Sig(p) = Sig(p) + VV/eps
+          Z(p)   = Z(p)   + VV/eps**2
 
         end do
       end do
@@ -66,8 +72,8 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
 
           eps = e0(p) + e0(i) - e0(a) - e0(b)
           VV  = (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)
-          Bpp(p,2) = Bpp(p,2) + VV/eps
-          Z(p)     = Z(p)     + VV/eps**2
+          Sig(p) = Sig(p) + VV/eps
+          Z(p)   = Z(p)   + VV/eps**2
 
         end do
       end do
@@ -78,15 +84,16 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
 
   if(linearize) then
 
-    eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
+    eGF2(:) = e0(:) + Z(:)*Sig(:)
 
   else
 
-    eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
+    eGF2(:) = e0(:) + Sig(:)
 
   end if
+
   ! Print results
 
-  call print_G0F2(nBas,nO,e0,eGF2,Z)
+  call print_G0F2(nBas,nO,e0,Sig,eGF2,Z)
 
 end subroutine G0F2
diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90
index 7a541d2..b6f5caa 100644
--- a/src/QuAcK/G0T0.f90
+++ b/src/QuAcK/G0T0.f90
@@ -71,7 +71,7 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, & 
+  call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, & 
                           nOOs,nVVs,eHF(:),ERI(:,:,:,:),  & 
                           Omega1s(:),X1s(:,:),Y1s(:,:),   & 
                           Omega2s(:),X2s(:,:),Y2s(:,:),   & 
@@ -96,7 +96,7 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, & 
+  call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, & 
                           nOOt,nVVt,eHF(:),ERI(:,:,:,:),  & 
                           Omega1t(:),X1t(:,:),Y1t(:,:),   & 
                           Omega2t(:),X2t(:,:),Y2t(:,:),   & 
@@ -113,6 +113,9 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
                                   X1t(:,:),Y1t(:,:),rho1t(:,:,:),                & 
                                   X2t(:,:),Y2t(:,:),rho2t(:,:,:))
 
+  rho2s(:,:,:) = 0d0
+  rho2t(:,:,:) = 0d0
+
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
@@ -133,7 +136,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 ! Solve the quasi-particle equation
 !----------------------------------------------
 
-  eG0T0(:) = eHF(:) + Z(:)*SigT(:)
+  eG0T0(:) = eHF(:) + SigT(:)
+! eG0T0(:) = eHF(:) + Z(:)*SigT(:)
 
 !----------------------------------------------
 ! Dump results
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 909e20d..a10d588 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -481,7 +481,7 @@ program QuAcK
   if(doppRPA) then
 
     call cpu_time(start_ppRPA)
-    call ppRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
+    call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
     call cpu_time(end_ppRPA)
 
     t_ppRPA = end_ppRPA - start_ppRPA
@@ -636,7 +636,8 @@ program QuAcK
     
     call cpu_time(start_G0T0)
     call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
-!   call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
+    call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
+
     call cpu_time(end_G0T0)
   
     t_G0T0 = end_G0T0 - start_G0T0
diff --git a/src/QuAcK/evGF2.f90 b/src/QuAcK/evGF2.f90
index e0abc77..27539f0 100644
--- a/src/QuAcK/evGF2.f90
+++ b/src/QuAcK/evGF2.f90
@@ -28,7 +28,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
   double precision              :: rcond
   double precision,allocatable  :: eGF2(:)
   double precision,allocatable  :: eOld(:)
-  double precision,allocatable  :: Bpp(:,:)
+  double precision,allocatable  :: Sig(:)
   double precision,allocatable  :: Z(:)
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: e_diis(:,:)
@@ -45,7 +45,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
 ! Memory allocation
 
-  allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
+  allocate(Sig(nBas),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
 
 ! Initialization
 
@@ -65,7 +65,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
     ! Frequency-dependent second-order contribution
 
-    Bpp(:,:) = 0d0
+    Sig(:) = 0d0
 
     do p=nC+1,nBas-nR
       do i=nC+1,nO
@@ -74,7 +74,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
             eps = eGF2(p) + e0(a) - e0(i) - e0(j)
 
-            Bpp(p,1) = Bpp(p,1) &
+            Sig(p) = Sig(p) &
                      + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
 
           end do
@@ -89,7 +89,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
             eps = eGF2(p) + e0(i) - e0(a) - e0(b)
 
-            Bpp(p,2) = Bpp(p,2) &
+            Sig(p) = Sig(p) &
                      + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
 
           end do
@@ -133,11 +133,11 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
 
     if(linearize) then
 
-      eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
+      eGF2(:) = e0(:) + Z(:)*Sig(:)
 
     else
 
-      eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
+      eGF2(:) = e0(:) + Sig(:)
 
     end if
 
diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90
index b9dfdb6..f541da9 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -49,6 +49,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 !                          + ERI(p,i,c,d)*X1(cd,ab) 
                            + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
                            /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
+
             end do
           end do
 
@@ -60,6 +61,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 !                          + ERI(p,i,k,l)*Y1(kl,ab)
                            + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
                            /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
+
             end do
           end do
 
@@ -86,7 +88,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k,nO
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-!                          + ERI(p,nO+a,k,l)*Y2(kl,ij)  
+!                          + ERI(p,nO+a,k,l)*Y2(kl,ij) 
                            + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
                            /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
 
@@ -116,6 +118,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
               cd = cd + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
                            + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
+              print*,rho1(p,i,ab),ERI(p,i,c,d),X1(cd,ab)
+
             end do
           end do
 
@@ -124,7 +128,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k+1,nO
               kl = kl + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
-                           + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
+                           + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) 
+              print*,rho1(p,i,ab),ERI(p,i,k,l),Y1(kl,ab)
             end do
           end do
 
@@ -139,7 +144,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do d=c+1,nBas-nR
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
+                           + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) 
+
             end do
           end do
 
@@ -148,7 +154,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k+1,nO
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
+                           + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) 
             end do
           end do
 
diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90
index b82f1a3..ae7d399 100644
--- a/src/QuAcK/linear_response_C_pp.f90
+++ b/src/QuAcK/linear_response_C_pp.f90
@@ -25,6 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 ! Define the chemical potential
  
   eF = e(nO) + e(nO+1)
+  eF = 0d0
 
 ! Build C matrix for the singlet manifold
 
diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90
index 1b38222..002a6bc 100644
--- a/src/QuAcK/linear_response_D_pp.f90
+++ b/src/QuAcK/linear_response_D_pp.f90
@@ -25,6 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 ! Define the chemical potential
 
   eF = e(nO) + e(nO+1)
+  eF = 0d0
  
 ! Build the D matrix for the singlet manifold
 
diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90
index 4a42f94..f178b7f 100644
--- a/src/QuAcK/linear_response_pp.f90
+++ b/src/QuAcK/linear_response_pp.f90
@@ -1,4 +1,4 @@
-subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & 
+subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & 
                               e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
 
 ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
@@ -8,6 +8,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
 
 ! Input variables
 
+  logical,intent(in)            :: ortho_eigvec
   logical,intent(in)            :: BSE
   integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR
   integer,intent(in)            :: nOO
@@ -127,11 +128,13 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
 
 ! Split the various quantities in p-p and h-h parts
 
-  call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
+  call sort_ppRPA(ortho_eigvec,nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
 
 ! Compute the RPA correlation energy
 
 ! EcppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
+! print*,+sum(Omega1(:)),- trace_matrix(nVV,C(:,:))
+! print*,-sum(Omega2(:)),- trace_matrix(nOO,D(:,:))
   EcppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
 ! EcppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
 
diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90
index 0d45d81..642fc23 100644
--- a/src/QuAcK/ppRPA.f90
+++ b/src/QuAcK/ppRPA.f90
@@ -1,4 +1,4 @@
-subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
+subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
 
 ! Perform pp-RPA calculation
 
@@ -7,7 +7,6 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
 
 ! Input variables
 
-  logical,intent(in)            :: doACFDT
   logical,intent(in)            :: singlet_manifold
   logical,intent(in)            :: triplet_manifold
   integer,intent(in)            :: nBas
@@ -62,7 +61,7 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
     allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
              Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
 
-    call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & 
+    call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & 
                             Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),  & 
                             Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),  & 
                             Ec_ppRPA(ispin))
@@ -91,9 +90,9 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
              Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
 
 
-    call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
-                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),  & 
-                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),  & 
+    call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
+                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),          & 
+                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),          & 
                             Ec_ppRPA(ispin))
 
     call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
diff --git a/src/QuAcK/print_G0F2.f90 b/src/QuAcK/print_G0F2.f90
index 1e80ac1..e560db8 100644
--- a/src/QuAcK/print_G0F2.f90
+++ b/src/QuAcK/print_G0F2.f90
@@ -1,4 +1,4 @@
-subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
+subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
 
 ! Print one-electron energies and other stuff for G0F2
 
@@ -8,6 +8,7 @@ subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
   integer,intent(in)                 :: nBas
   integer,intent(in)                 :: nO
   double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: Sig(nBas)
   double precision,intent(in)        :: eGF2(nBas)
   double precision,intent(in)        :: Z(nBas)
 
@@ -24,23 +25,23 @@ subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
 
 ! Dump results
 
-  write(*,*)'--------------------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,*)' Frequency-dependent G0F2 calculation'
-  write(*,*)'--------------------------------------------------------'
-  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
-            '|','#','|','e_HF (eV)','|','e_G0F2 (eV)','|','Z','|'
-  write(*,*)'--------------------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Sigma (eV)','|','e_G0F2 (eV)','|','Z','|'
+  write(*,*)'--------------------------------------------------------------------------'
 
   do p=1,nBas
-    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
-    '|',p,'|',eHF(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
+    '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
   enddo
 
-  write(*,*)'--------------------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,'(2X,A27,F15.6)') 'G0F2  HOMO      energy (eV):',eGF2(HOMO)*HaToeV
   write(*,'(2X,A27,F15.6)') 'G0F2  LUMO      energy (eV):',eGF2(LUMO)*HaToeV
   write(*,'(2X,A27,F15.6)') 'G0F2  HOMO-LUMO gap    (eV):',Gap*HaToeV
-  write(*,*)'--------------------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,*)
 
 end subroutine print_G0F2
diff --git a/src/QuAcK/print_G0T0.f90 b/src/QuAcK/print_G0T0.f90
index 6b8bcee..165fa9b 100644
--- a/src/QuAcK/print_G0T0.f90
+++ b/src/QuAcK/print_G0T0.f90
@@ -21,7 +21,16 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA)
 
   HOMO = nO
   LUMO = HOMO + 1
-  Gap = eGW(LUMO)-eGW(HOMO)
+  if(nBas >= LUMO) then
+
+    Gap = eGW(LUMO) - eGW(HOMO)
+
+  else
+
+    Gap = 0d0
+
+  end if
+
 
 ! Dump results
 
diff --git a/src/QuAcK/print_evGF2.f90 b/src/QuAcK/print_evGF2.f90
index 050dc7c..fdab254 100644
--- a/src/QuAcK/print_evGF2.f90
+++ b/src/QuAcK/print_evGF2.f90
@@ -1,44 +1,53 @@
-subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
+subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
 
-! Print one-electron energies and other stuff for GF2
+! Print one-electron energies and other stuff for G0F2
 
   implicit none
   include 'parameters.h'
 
-  integer,intent(in)                 :: nBas,nO,nSCF
-  double precision,intent(in)        :: Conv,eHF(nBas),eGF2(nBas)
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO
+  integer,intent(in)                 :: nSCF
+  double precision,intent(in)        :: Conv
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: Sig(nBas)
+  double precision,intent(in)        :: eGF2(nBas)
+  double precision,intent(in)        :: Z(nBas)
 
-  integer                            :: x,HOMO,LUMO
+  integer                            :: p
+  integer                            :: HOMO
+  integer                            :: LUMO
   double precision                   :: Gap
 
 ! HOMO and LUMO
 
   HOMO = nO
   LUMO = HOMO + 1
-  Gap = eGF2(LUMO)-eGF2(HOMO)
+  Gap = eGF2(LUMO) - eGF2(HOMO)
 
 ! Dump results
 
-  write(*,*)'-------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,*)' Frequency-dependent evGF2 calculation'
-  write(*,*)'-------------------------------------------'
-  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
-            '|','#','|','e_HF (eV)','|','e_evGF2 (eV)','|'
-  write(*,*)'-------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Sigma (eV)','|','e_evGF2 (eV)','|','Z','|'
+  write(*,*)'--------------------------------------------------------------------------'
 
-  do x=1,nBas
-    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-    '|',x,'|',eHF(x)*HaToeV,'|',eGF2(x)*HaToeV,'|'
+  do p=1,nBas
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
+    '|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
   enddo
 
-  write(*,*)'-------------------------------------------'
+  write(*,*)'-------------------------------------------------------------'
   write(*,'(2X,A10,I3)')   'Iteration ',nSCF
   write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
-  write(*,*)'-------------------------------------------'
+
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,'(2X,A27,F15.6)') 'evGF2  HOMO      energy (eV):',eGF2(HOMO)*HaToeV
   write(*,'(2X,A27,F15.6)') 'evGF2  LUMO      energy (eV):',eGF2(LUMO)*HaToeV
   write(*,'(2X,A27,F15.6)') 'evGF2  HOMO-LUMO gap    (eV):',Gap*HaToeV
-  write(*,*)'-------------------------------------------'
+  write(*,*)'--------------------------------------------------------------------------'
   write(*,*)
 
 end subroutine print_evGF2
diff --git a/src/QuAcK/print_evGF3.f90 b/src/QuAcK/print_evGF3.f90
index 3e7f740..6747058 100644
--- a/src/QuAcK/print_evGF3.f90
+++ b/src/QuAcK/print_evGF3.f90
@@ -8,7 +8,7 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
   integer,intent(in)                 :: nBas,nO,nSCF
   double precision,intent(in)        :: Conv,eHF(nBas),eGF3(nBas),Z(nBas)
 
-  integer                            :: x,HOMO,LUMO
+  integer                            :: p,HOMO,LUMO
   double precision                   :: Gap
 
 ! HOMO and LUMO
@@ -26,9 +26,9 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
             '|','#','|','e_HF (eV)','|','Z','|','e_evGF3 (eV)','|'
   write(*,*)'-------------------------------------------------------------'
 
-  do x=1,nBas
+  do p=1,nBas
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-    '|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|'
+    '|',p,'|',eHF(p)*HaToeV,'|',Z(p),'|',eGF3(p)*HaToeV,'|'
   enddo
 
   write(*,*)'-------------------------------------------------------------'
diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90
index 1eb76fe..ce4be02 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix.f90
@@ -42,7 +42,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do i=nC+1,nO
       do cd=1,nVVs
         eps = e(p) + e(i) - Omega1s(cd)
-        Z(p) = Z(p) + 1d0*(rho1s(p,i,cd)/eps)**2
+        Z(p) = Z(p) + (rho1s(p,i,cd)/eps)**2
       enddo
     enddo
   enddo
@@ -53,7 +53,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do a=1,nV-nR
       do kl=1,nOOs
         eps = e(p) + e(nO+a) - Omega2s(kl)
-        Z(p) = Z(p) + 1d0*(rho2s(p,a,kl)/eps)**2
+        Z(p) = Z(p) + (rho2s(p,a,kl)/eps)**2
       enddo
     enddo
   enddo
@@ -68,7 +68,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do i=nC+1,nO
       do cd=1,nVVt
         eps = e(p) + e(i) - Omega1t(cd)
-        Z(p) = Z(p) + 1d0*(rho1t(p,i,cd)/eps)**2
+        Z(p) = Z(p) + (rho1t(p,i,cd)/eps)**2
       enddo
     enddo
   enddo
@@ -79,7 +79,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
     do a=1,nV-nR
       do kl=1,nOOt
         eps = e(p) + e(nO+a) - Omega2t(kl)
-        Z(p) = Z(p) + 1d0*(rho2t(p,a,kl)/eps)**2
+        Z(p) = Z(p) + (rho2t(p,a,kl)/eps)**2
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90
index 72d7963..5e509ba 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag.f90
@@ -46,7 +46,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do i=nC+1,nO
       do cd=1,nVVs
         eps = e(p) + e(i) - Omega1s(cd)
-        SigT(p) = SigT(p) + 1d0*rho1s(p,i,cd)**2/eps
+        SigT(p) = SigT(p) + rho1s(p,i,cd)**2/eps
       enddo
     enddo
   enddo
@@ -57,7 +57,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do a=1,nV-nR
       do kl=1,nOOs
         eps = e(p) + e(nO+a) - Omega2s(kl)
-        SigT(p) = SigT(p) + 1d0*rho2s(p,a,kl)**2/eps
+        SigT(p) = SigT(p) + rho2s(p,a,kl)**2/eps
       enddo
     enddo
   enddo
@@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do i=nC+1,nO
       do cd=1,nVVt
         eps = e(p) + e(i) - Omega1t(cd)
-        SigT(p) = SigT(p) + 1d0*rho1t(p,i,cd)**2/eps
+        SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps
       enddo
     enddo
   enddo
@@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do a=1,nV-nR
       do kl=1,nOOt
         eps = e(p) + e(nO+a) - Omega2t(kl)
-        SigT(p) = SigT(p) + 1d0*rho2t(p,a,kl)**2/eps
+        SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90
index 03b1ba9..bf57a4f 100644
--- a/src/QuAcK/soG0T0.f90
+++ b/src/QuAcK/soG0T0.f90
@@ -79,8 +79,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, & 
-                          seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2,         & 
+  call linear_response_pp(ispin,.true.,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, & 
+                          seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2,                & 
                           EcRPA)
 
   call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
@@ -92,6 +92,10 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
                                   X1(:,:),Y1(:,:),rho1(:,:,:), &
                                   X2(:,:),Y2(:,:),rho2(:,:,:))
 
+  
+
+  rho2(:,:,:) = 0d0
+
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
@@ -112,7 +116,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 ! Solve the quasi-particle equation
 !----------------------------------------------
 
-  eG0T0(:) = seHF(:) + Z(:)*SigT(:)
+  eG0T0(:) = seHF(:) + SigT(:)
+! eG0T0(:) = seHF(:) + Z(:)*SigT(:)
 
 !----------------------------------------------
 ! Dump results
diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90
index 90ba9c0..c3e7162 100644
--- a/src/QuAcK/sort_ppRPA.f90
+++ b/src/QuAcK/sort_ppRPA.f90
@@ -1,4 +1,4 @@
-subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
+subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
 ! Compute the metric matrix for pp-RPA
 
@@ -7,6 +7,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
 ! Input variables
 
+  logical,intent(in)            :: ortho_eigvec
   integer,intent(in)            :: nOO
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: Omega(nOO+nVV)
@@ -91,24 +92,28 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
   if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
 
-  write(*,*) 'pp-RPA positive excitation energies'
-  call matout(nVV,1,Omega1(:))
-  write(*,*)
+! write(*,*) 'pp-RPA positive excitation energies'
+! call matout(nVV,1,Omega1(:))
+! write(*,*)
 
-  write(*,*) 'pp-RPA negative excitation energies'
-  call matout(nOO,1,Omega2(:))
-  write(*,*)
+! write(*,*) 'pp-RPA negative excitation energies'
+! call matout(nOO,1,Omega2(:))
+! write(*,*)
 
 ! Orthogonalize eigenvectors
 
-  S1 = + matmul(transpose(Z1),matmul(M,Z1))
-  S2 = - matmul(transpose(Z2),matmul(M,Z2))
+  if(ortho_eigvec) then
 
-  call orthogonalization_matrix(1,nVV,S1,O1)
-  call orthogonalization_matrix(1,nOO,S2,O2)
+    S1 = + matmul(transpose(Z1),matmul(M,Z1))
+    S2 = - matmul(transpose(Z2),matmul(M,Z2))
 
-  Z1 = matmul(Z1,O1)
-  Z2 = matmul(Z2,O2)
+    if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
+    if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
+
+    Z1 = matmul(Z1,O1)
+    Z2 = matmul(Z2,O2)
+
+  end if
 
   X1(1:nVV,1:nVV) = Z1(    1:    nVV,1:nVV)
   Y1(1:nOO,1:nVV) = Z1(nVV+1:nOO+nVV,1:nVV)
diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile
index 448a2e1..159acd5 100644
--- a/src/eDFT/Makefile
+++ b/src/eDFT/Makefile
@@ -14,8 +14,8 @@ ifeq ($(PROFIL),1)
   FC += -p -fno-inline
 endif
 
-LIBS = ../../lib/*.a
-LIBS += -lblas -llapack
+LIBS = -lstdc++ ../../lib/*.a
+LIBS += -lblas -llapack 
 
 SRCF90 = $(wildcard *.f90)
 
@@ -23,15 +23,18 @@ SRC = $(wildcard *.F)
 
 OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.F,$(ODIR)/%.o,$(SRC))
 
-$(ODIR)/%.o: %.f90 
-	$(FC) -c -o $@ $< $(FFLAGS) 
-
-$(ODIR)/%.o: %.F
-	$(FC) -c -o $@ $< $(FFLAGS) 
-
 $(BDIR)/eDFT: $(OBJ)
 	$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
 
+numgrid.mod: numgrid.f90 
+	$(FC) -c -o $@ $< $(FFLAGS) 
+
+$(ODIR)/%.o: %.f90 numgrid.mod
+	$(FC) -c -o $@ $< $(FFLAGS) 
+
+$(ODIR)/%.o: %.F numgrid.mod
+	$(FC) -c -o $@ $< $(FFLAGS) 
+
 debug: 
 	DEBUG=1 make $(BDIR)/eDFT
 
@@ -39,4 +42,4 @@ profil:
 	PROFIL=1 make $(BDIR)/eDFT
 
 clean:
-	rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug
+	rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug numgrid.mod
diff --git a/src/eDFT/numgrid.f90 b/src/eDFT/numgrid.f90
new file mode 100644
index 0000000..2717e6b
--- /dev/null
+++ b/src/eDFT/numgrid.f90
@@ -0,0 +1,113 @@
+module numgrid
+
+   use, intrinsic :: iso_c_binding, only: c_ptr, c_double, c_int
+
+   implicit none
+
+   public numgrid_new_atom_grid
+   public numgrid_free_atom_grid
+   public numgrid_get_num_grid_points
+   public numgrid_get_num_radial_grid_points
+   public numgrid_get_grid
+   public numgrid_get_radial_grid
+   public numgrid_get_angular_grid
+
+   private
+
+   interface numgrid_new_atom_grid
+      function numgrid_new_atom_grid(radial_precision,       &
+                                     min_num_angular_points, &
+                                     max_num_angular_points, &
+                                     proton_charge,          &
+                                     alpha_max,              &
+                                     max_l_quantum_number,   &
+                                     alpha_min) result(context) bind (c)
+         import :: c_ptr, c_double, c_int
+         type(c_ptr) :: context
+         real(c_double), intent(in), value :: radial_precision
+         integer(c_int), intent(in), value :: min_num_angular_points
+         integer(c_int), intent(in), value :: max_num_angular_points
+         integer(c_int), intent(in), value :: proton_charge
+         real(c_double), intent(in), value :: alpha_max
+         integer(c_int), intent(in), value :: max_l_quantum_number
+         real(c_double), intent(in)        :: alpha_min(*)
+      end function
+   end interface
+
+   interface numgrid_free_atom_grid
+      subroutine numgrid_free_atom_grid(context) bind (c)
+         import :: c_ptr
+         type(c_ptr), value :: context
+      end subroutine
+   end interface
+
+   interface numgrid_get_num_grid_points
+      function numgrid_get_num_grid_points(context) result(num_grid_points) bind (c)
+         import :: c_ptr, c_int
+         type(c_ptr), value :: context
+         integer(c_int)     :: num_grid_points
+      end function
+   end interface
+
+   interface numgrid_get_num_radial_grid_points
+      function numgrid_get_num_radial_grid_points(context) result(num_radial_grid_points) bind (c)
+         import :: c_ptr, c_int
+         type(c_ptr), value :: context
+         integer(c_int)     :: num_radial_grid_points
+      end function
+   end interface
+
+   interface numgrid_get_grid
+      subroutine numgrid_get_grid(context,          &
+                                  num_centers,      &
+                                  center_index,     &
+                                  x_coordinates_bohr, &
+                                  y_coordinates_bohr, &
+                                  z_coordinates_bohr, &
+                                  proton_charges,   &
+                                  grid_x_bohr,        &
+                                  grid_y_bohr,        &
+                                  grid_z_bohr,        &
+                                  grid_w) bind (c)
+         import :: c_ptr, c_double, c_int
+         type(c_ptr), value                :: context
+         integer(c_int), intent(in), value :: num_centers
+         integer(c_int), intent(in), value :: center_index
+         real(c_double), intent(in)        :: x_coordinates_bohr(*)
+         real(c_double), intent(in)        :: y_coordinates_bohr(*)
+         real(c_double), intent(in)        :: z_coordinates_bohr(*)
+         integer(c_int), intent(in)        :: proton_charges(*)
+         real(c_double), intent(inout)     :: grid_x_bohr(*)
+         real(c_double), intent(inout)     :: grid_y_bohr(*)
+         real(c_double), intent(inout)     :: grid_z_bohr(*)
+         real(c_double), intent(inout)     :: grid_w(*)
+      end subroutine
+   end interface
+
+   interface numgrid_get_radial_grid
+      subroutine numgrid_get_radial_grid(context,  &
+                                         radial_grid_r_bohr,&
+                                         radial_grid_w) bind (c)
+         import :: c_ptr, c_double, c_int
+         type(c_ptr), value            :: context
+         real(c_double), intent(inout) :: radial_grid_r_bohr(*)
+         real(c_double), intent(inout) :: radial_grid_w(*)
+      end subroutine
+   end interface
+
+   interface numgrid_get_angular_grid
+      subroutine numgrid_get_angular_grid(num_angular_grid_points, &
+                                          angular_grid_x_bohr,       &
+                                          angular_grid_y_bohr,       &
+                                          angular_grid_z_bohr,       &
+                                          angular_grid_w) bind (c)
+         import :: c_double, c_int
+         integer(c_int), intent(in), value :: num_angular_grid_points
+         real(c_double), intent(inout)     :: angular_grid_x_bohr(*)
+         real(c_double), intent(inout)     :: angular_grid_y_bohr(*)
+         real(c_double), intent(inout)     :: angular_grid_z_bohr(*)
+         real(c_double), intent(inout)     :: angular_grid_w(*)
+      end subroutine
+   end interface
+
+end module
diff --git a/src/eDFT/test_explicit.f90 b/src/eDFT/test_explicit.f90
new file mode 100644
index 0000000..2214109
--- /dev/null
+++ b/src/eDFT/test_explicit.f90
@@ -0,0 +1,156 @@
+subroutine test_numgrid
+
+   ! in this test we compute the grid for O
+   ! in the presence of two H in the H2O molecule
+
+   use numgrid
+   use, intrinsic :: iso_c_binding, only: c_ptr
+
+   implicit none
+
+   real(8)              :: radial_precision
+   integer              :: min_num_angular_points
+   integer              :: max_num_angular_points
+   integer              :: num_centers
+   integer              :: proton_charges(3)
+   real(8)              :: x_coordinates_bohr(3)
+   real(8)              :: y_coordinates_bohr(3)
+   real(8)              :: z_coordinates_bohr(3)
+   real(8)              :: alpha_max
+   integer              :: max_l_quantum_number
+   real(8)              :: alpha_min(3)
+   integer              :: num_points
+   integer              :: num_radial_points
+   integer              :: num_angular_points
+   real(8), allocatable :: grid_x_bohr(:)
+   real(8), allocatable :: grid_y_bohr(:)
+   real(8), allocatable :: grid_z_bohr(:)
+   real(8), allocatable :: grid_w(:)
+   real(8), allocatable :: angular_grid_x_bohr(:)
+   real(8), allocatable :: angular_grid_y_bohr(:)
+   real(8), allocatable :: angular_grid_z_bohr(:)
+   real(8), allocatable :: angular_grid_w(:)
+   real(8), allocatable :: radial_grid_r_bohr(:)
+   real(8), allocatable :: radial_grid_w(:)
+   integer              :: center_index
+   integer, parameter   :: io_unit = 13
+   real(8)              :: ref(4), own(4)
+   integer              :: ipoint
+   integer              :: j
+   real(8)              :: error
+   type(c_ptr)          :: context
+
+   radial_precision = 1.0d-12
+   min_num_angular_points = 86
+   max_num_angular_points = 302
+
+   num_centers = 3
+
+   proton_charges(1) = 8
+   proton_charges(2) = 1
+   proton_charges(3) = 1
+
+   x_coordinates_bohr(1) = 0.0d0
+   x_coordinates_bohr(2) = 1.43d0
+   x_coordinates_bohr(3) = -1.43d0
+
+   y_coordinates_bohr(1) = 0.0d0
+   y_coordinates_bohr(2) = 0.0d0
+   y_coordinates_bohr(3) = 0.0d0
+
+   z_coordinates_bohr(1) = 0.0d0
+   z_coordinates_bohr(2) = 1.1d0
+   z_coordinates_bohr(3) = 1.1d0
+
+   ! oxygen
+   alpha_max = 11720.0d0
+   max_l_quantum_number = 2
+   alpha_min(1) = 0.3023d0
+   alpha_min(2) = 0.2753d0
+   alpha_min(3) = 1.185d0
+
+   open(unit=io_unit, file='../test/reference_grid/cc-pVDZ.txt', access='sequential', action='read')
+
+   context = numgrid_new_atom_grid(radial_precision,       &
+                                   min_num_angular_points, &
+                                   max_num_angular_points, &
+                                   proton_charges(1),      &
+                                   alpha_max,              &
+                                   max_l_quantum_number,   &
+                                   alpha_min)
+
+   num_points = numgrid_get_num_grid_points(context)
+   if (num_points /= 16364) stop 1
+
+   allocate(grid_x_bohr(num_points))
+   allocate(grid_y_bohr(num_points))
+   allocate(grid_z_bohr(num_points))
+   allocate(grid_w(num_points))
+
+   center_index = 0
+
+   call numgrid_get_grid(context,          &
+                         num_centers,      &
+                         center_index,     &
+                         x_coordinates_bohr, &
+                         y_coordinates_bohr, &
+                         z_coordinates_bohr, &
+                         proton_charges,   &
+                         grid_x_bohr,        &
+                         grid_y_bohr,        &
+                         grid_z_bohr,        &
+                         grid_w)
+
+   do ipoint = 1, num_points
+      read(io_unit, *) ref(1), ref(2), ref(3), ref(4)
+      own(1) = grid_x_bohr(ipoint)
+      own(2) = grid_y_bohr(ipoint)
+      own(3) = grid_z_bohr(ipoint)
+      own(4) = grid_w(ipoint)
+      do j = 1, 4
+         error = own(j) - ref(j)
+         if (dabs(ref(j)) > 1.0e-15) error = error/ref(j)
+         if (dabs(error) > 1.0e-5) stop 1
+      end do
+   end do
+   close(io_unit)
+
+   deallocate(grid_x_bohr)
+   deallocate(grid_y_bohr)
+   deallocate(grid_z_bohr)
+   deallocate(grid_w)
+
+   num_radial_points = numgrid_get_num_radial_grid_points(context)
+   if (num_radial_points /= 106) stop 1
+
+   allocate(radial_grid_r_bohr(num_radial_points))
+   allocate(radial_grid_w(num_radial_points))
+
+   call numgrid_get_radial_grid(context,          &
+                                radial_grid_r_bohr, &
+                                radial_grid_w)
+
+   deallocate(radial_grid_r_bohr)
+   deallocate(radial_grid_w)
+
+   num_angular_points = 14
+
+   allocate(angular_grid_x_bohr(num_angular_points))
+   allocate(angular_grid_y_bohr(num_angular_points))
+   allocate(angular_grid_z_bohr(num_angular_points))
+   allocate(angular_grid_w(num_angular_points))
+
+   call numgrid_get_angular_grid(num_angular_points, &
+                                 angular_grid_x_bohr,  &
+                                 angular_grid_y_bohr,  &
+                                 angular_grid_z_bohr,  &
+                                 angular_grid_w)
+
+   deallocate(angular_grid_x_bohr)
+   deallocate(angular_grid_y_bohr)
+   deallocate(angular_grid_z_bohr)
+   deallocate(angular_grid_w)
+
+   call numgrid_free_atom_grid(context)
+
+end subroutine