diff --git a/basis_cap/install.py b/basis_cap/install.py
new file mode 100755
index 0000000..a6db931
--- /dev/null
+++ b/basis_cap/install.py
@@ -0,0 +1,107 @@
+#!/usr/bin/env python3
+
+"""Installs the default basis set files from the BSE command-line tool in the required format for """
+
+import basis_set_exchange as bse
+conversion = {
+    "3-21g":  "3-21G",
+    "4-31g":  "4-31G",
+    "6-311g":  "6-311G",
+    "6-311++g_2d_2p":  "6-311++G(2d,2p)",
+    "6-311g_2df_2pd":  "6-311G(2df,2pd)",
+    "6-311++g_3df_3pd":  "6-311++G(3df,3pd)",
+    "6-311+g_star":  "6-311+G*",
+    "6-311G_star":  "6-311G*",
+    "6-311++g_star_star":  "6-311++G**",
+    "6-311G_star_star":  "6-311G**",
+    "6-31g":  "6-31G",
+    "6-31+g":  "6-31+G",
+    "6-31++g":  "6-31++G",
+    "6-31g_3df_3pd":  "6-31G(3df,3pd)",
+    "6-31++g_star":  "6-31++G*",
+    "6-31+g_star":  "6-31+G*",
+    "6-31g_star":  "6-31G*",
+    "6-31++g_star_star":  "6-31++G**",
+    "6-31g_star_star":  "6-31G**",
+    "ano2_ames":  "NASA Ames ANO2",
+    "ano_ames":  "NASA Ames ANO",
+    "ano-rcc":  "ANO-RCC",
+    "apr-cc-pv_q+d_z":  "apr-cc-pV(Q+d)Z",
+    "aug-cc-pcv5z":  "aug-cc-pCV5Z",
+    "aug-cc-pcvdz":  "aug-cc-pCVDZ",
+    "aug-cc-pcvqz":  "aug-cc-pCVQZ",
+    "aug-cc-pcvtz":  "aug-cc-pCVTZ",
+    "aug-cc-pv5z":  "aug-cc-pV5Z",
+    "aug-cc-pv6z":  "aug-cc-pV6Z",
+    "aug-cc-pvdz":  "aug-cc-pVDZ",
+    "aug-cc-pvqz":  "aug-cc-pVQZ",
+    "aug-cc-pvtz":  "aug-cc-pVTZ",
+    "aug-cc-pwcv5z":  "aug-cc-pwCV5Z",
+    "aug-cc-pwcvdz":  "aug-cc-pwCVDZ",
+    "aug-cc-pwcvqz":  "aug-cc-pwCVQZ",
+    "aug-cc-pwcvtz":  "aug-cc-pwCVTZ",
+    "cc-pcv5z":  "cc-pCV5Z",
+    "cc-pcvdz":  "cc-pCVDZ",
+    "cc-pcvqz":  "cc-pCVQZ",
+    "cc-pcvtz":  "cc-pCVTZ",
+    "cc-pv5z":  "cc-pV5Z",
+    "cc-pv6z":  "cc-pV6Z",
+    "cc-pv8z":  "cc-pV8Z",
+    "cc-pv9z":  "cc-pV9Z",
+    "cc-pvdz":  "cc-pVDZ",
+    "cc-pvqz":  "cc-pVQZ",
+    "cc-pvtz":  "cc-pVTZ",
+    "cc-pwcv5z":  "cc-pwCV5Z",
+    "cc-pwcvdz":  "cc-pwCVDZ",
+    "cc-pwcvqz":  "cc-pwCVQZ",
+    "cc-pwcvtz":  "cc-pwCVTZ",
+    "d-aug-cc-pv5z":  "d-aug-cc-pV5Z",
+    "d-aug-cc-pv6z":  "d-aug-cc-pV6Z",
+    "d-aug-cc-pvdz":  "d-aug-cc-pVDZ",
+    "d-aug-cc-pvqz":  "d-aug-cc-pVQZ",
+    "d-aug-cc-pvtz":  "d-aug-cc-pVTZ",
+    "def2-qzvp":  "def2-QZVP",
+    "def2-qzvpd":  "def2-QZVPD",
+    "def2-qzvpp":  "def2-QZVPP",
+    "def2-qzvppd":  "def2-QZVPPD",
+    "def2-svp":  "def2-SVP",
+    "def2-sv_p":  "def2-SV(P)",
+    "def2-svpd":  "def2-SVPD",
+    "def2-tzvp":  "def2-TZVP",
+    "def2-tzvpd":  "def2-TZVPD",
+    "def2-tzvpp":  "def2-TZVPP",
+    "def2-tzvppd":  "def2-TZVPPD",
+    "jul-cc-pv_d+d_z":  "jul-cc-pV(D+d)Z",
+    "jul-cc-pv_q+d_z":  "jul-cc-pV(Q+d)Z",
+    "jul-cc-pv_t+d_z":  "jul-cc-pV(T+d)Z",
+    "jun-cc-pv_d+d_z":  "jun-cc-pV(D+d)Z",
+    "jun-cc-pv_q+d_z":  "jun-cc-pV(Q+d)Z",
+    "jun-cc-pv_t+d_z":  "jun-cc-pV(T+d)Z",
+    "lanl08":  "LANL08",
+    "lanl08+":  "LANL08+",
+    "lanl08d":  "LANL08(d)",
+    "lanl08f":  "LANL08(f)",
+    "m6-31g":  "m6-31G",
+    "may-cc-pv_q+d_z":  "may-cc-pV(Q+d)Z",
+    "may-cc-pv_t+d_z":  "may-cc-pV(T+d)Z",
+    "midi_bang":  "MIDI!",
+    "midi_huzinaga":  "MIDI",
+    "mini":  "MINI",
+    "pv6z":  "pV6Z",
+    "pv7z":  "pV7Z",
+    "pvdz_ahlrichs":  "Ahlrichs pVDZ",
+    "sto-2g":  "STO-2G",
+    "sto-3g":  "STO-3G",
+    "sto-4g":  "STO-4G",
+    "sto-5g":  "STO-5G",
+    "sto-6g":  "STO-6G",
+    "sto-3g_star":  "STO-3G*",
+}
+
+
+for filename in conversion.keys():
+    print(filename)
+    data = bse.get_basis(
+        conversion[filename], fmt='psi4', uncontract_general=True, uncontract_spdf=True)
+    with open(filename, 'w') as f:
+        f.write(data)
diff --git a/cap_onsets/CO b/cap_onsets/CO
new file mode 100644
index 0000000..efd25f6
--- /dev/null
+++ b/cap_onsets/CO
@@ -0,0 +1 @@
+2.76 2.76 4.97
diff --git a/cap_onsets/N2 b/cap_onsets/N2
new file mode 100644
index 0000000..2caa419
--- /dev/null
+++ b/cap_onsets/N2
@@ -0,0 +1 @@
+2.76 2.76 4.88
diff --git a/mol/COX.xyz b/mol/COX.xyz
new file mode 100644
index 0000000..0571bda
--- /dev/null
+++ b/mol/COX.xyz
@@ -0,0 +1,5 @@
+3
+
+C  0.0000  0.0000  0.0000
+O  0.0000  0.0000  1.1335 
+X  0.0000  0.0000  0.56675
diff --git a/mol/N2X.xyz b/mol/N2X.xyz
new file mode 100644
index 0000000..5a6ac60
--- /dev/null
+++ b/mol/N2X.xyz
@@ -0,0 +1,5 @@
+3
+
+N  0.0000  0.0000  -0.5475132
+N  0.0000  0.0000  0.5475132
+X  0.0000  0.0000  0.000000
diff --git a/src/GW/cRGW_Im_SigC.f90 b/src/GW/cRGW_Im_SigC.f90
new file mode 100644
index 0000000..d038c23
--- /dev/null
+++ b/src/GW/cRGW_Im_SigC.f90
@@ -0,0 +1,56 @@
+double precision function cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: Re_w
+  double precision,intent(in)   :: Im_w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Re_e(nBas)
+  double precision,intent(in)   :: Im_e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps,eta_tilde
+
+! Initialize 
+
+  cRGW_Im_SigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = Re_w - Re_e(i) + Om(m)
+      eta_tilde = eta - Im_w + Im_e(i)
+      num = 2d0*rho(p,i,m)**2
+      cRGW_Im_SigC = cRGW_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2)
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = Re_w - Re_e(a) - Om(m)
+      eta_tilde = eta + Im_w - Im_e(a)
+      num = 2d0*rho(p,a,m)**2
+      cRGW_Im_SigC =cRGW_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2)
+    end do
+  end do
+
+end function 
diff --git a/src/GW/cRGW_Im_dSigC.f90 b/src/GW/cRGW_Im_dSigC.f90
new file mode 100644
index 0000000..13c8fce
--- /dev/null
+++ b/src/GW/cRGW_Im_dSigC.f90
@@ -0,0 +1,56 @@
+double precision function cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
+
+! Compute the derivative of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: Re_w
+  double precision,intent(in)   :: Im_w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Re_e(nBas)
+  double precision,intent(in)   :: Im_e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps,eta_tilde
+
+! Initialize 
+
+  cRGW_Im_dSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = Re_w - Re_e(i) + Om(m)
+      eta_tilde = eta - Im_w + Im_e(i)
+      num = 2d0*rho(p,i,m)**2
+      cRGW_Im_dSigC = cRGW_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = Re_w - Re_e(a) - Om(m)
+      eta_tilde = eta + Im_w - Im_e(a)
+      num = 2d0*rho(p,a,m)**2
+      cRGW_Im_dSigC = cRGW_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2
+    end do
+  end do
+
+end function 
diff --git a/src/GW/cRGW_QP_graph.f90 b/src/GW/cRGW_QP_graph.f90
new file mode 100644
index 0000000..bb8ac6b
--- /dev/null
+++ b/src/GW/cRGW_QP_graph.f90
@@ -0,0 +1,107 @@
+subroutine cRGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin, &
+        Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_Z)
+
+! Compute the graphical solution of the QP equation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+
+  logical,intent(in)            :: doSRG
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: flow
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: e_cap(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+  double precision,intent(in)   :: Re_eGWlin(nBas)
+  double precision,intent(in)   :: Im_eGWlin(nBas)
+  double precision,external     :: cRGW_Re_SigC,cRGW_Re_dSigC
+  double precision,external     :: cRGW_Im_SigC,cRGW_Im_dSigC
+  double precision,intent(in)   :: Re_eOld(nBas)
+  double precision,intent(in)   :: Im_eOld(nBas)
+
+! Local variables
+
+  integer                       :: p
+  integer                       :: nIt
+  integer,parameter             :: maxIt = 64
+  double precision,parameter    :: thresh = 1d-6
+  double precision              :: Re_SigC,Re_dSigC
+  double precision              :: Im_SigC,Im_dSigC
+  double precision              :: Re_f,Im_f,Re_df,Im_df
+  double precision              :: Re_w
+  double precision              :: Im_w
+
+! Output variables
+
+  double precision,intent(out)  :: Re_eGW(nBas),Im_eGW(nBas)
+  double precision,intent(out)  :: Re_Z(nBas),Im_Z(nBas)
+
+! Run Newton's algorithm to find the root
+
+  write(*,*)'-----------------------------------------------------'
+  write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Re(e_GWlin) (eV)','Re(e_GW (eV))','Re(Z)'
+  write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Im(e_GWlin) (eV)','Im(e_GW (eV))','Im(Z)'
+  write(*,*)'-----------------------------------------------------'
+
+  do p=nC+1,nBas-nR
+
+    Re_w = Re_eGWlin(p)
+    Im_w = Im_eGWlin(p)
+    nIt = 0
+    Re_f = 1d0
+    Im_f = 1d0
+    
+    do while (sqrt(Re_f**2+Im_f**2) > thresh .and. nIt < maxIt)
+    
+      nIt = nIt + 1
+
+
+      Re_SigC  = cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
+      Im_SigC  = cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
+      Re_dSigC = cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
+      Im_dSigC = cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho)
+
+
+      Re_f  = Re_w - eHF(p) - Re_SigC
+      Im_f  = Im_w - e_cap(p) - Im_SigC
+      Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
+      Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2)
+      Re_w = Re_w - Re_df*Re_f + Im_df*Im_f
+      Im_w = Im_w - Re_f*Im_df - Re_df*Im_f
+    
+    end do
+ 
+    if(nIt == maxIt) then 
+
+      Re_eGW(p) = Re_eGWlin(p)
+      write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p),'Cvg Failed!'
+
+    else
+
+      Re_eGW(p) = Re_w
+      Im_eGW(p) = Im_w
+      Re_Z(p)   = Re_df
+      Im_Z(p)   = Im_df
+
+      write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p)
+      write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGWlin(p)*HaToeV,Im_eGW(p)*HaToeV,Im_Z(p)
+
+    end if
+          
+  end do
+
+  write(*,*)'-----------------------------------------------------'
+  write(*,*)
+
+end subroutine 
diff --git a/src/GW/cRGW_Re_SigC.f90 b/src/GW/cRGW_Re_SigC.f90
new file mode 100644
index 0000000..ca85e39
--- /dev/null
+++ b/src/GW/cRGW_Re_SigC.f90
@@ -0,0 +1,57 @@
+double precision function cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: Re_w
+  double precision,intent(in)   :: Im_w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Re_e(nBas)
+  double precision,intent(in)   :: Im_e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps
+  double precision              :: eta_tilde
+
+! Initialize 
+
+  cRGW_Re_SigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = Re_w - Re_e(i) + Om(m)
+      eta_tilde = eta - Im_w + Im_e(i)
+      num = 2d0*rho(p,i,m)**2
+      cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = Re_w - Re_e(a) - Om(m)
+      eta_tilde = eta + Im_w - Im_e(a)
+      num = 2d0*rho(p,a,m)**2
+      cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2)
+    end do
+  end do
+
+end function 
diff --git a/src/GW/cRGW_Re_dSigC.f90 b/src/GW/cRGW_Re_dSigC.f90
new file mode 100644
index 0000000..33ba960
--- /dev/null
+++ b/src/GW/cRGW_Re_dSigC.f90
@@ -0,0 +1,56 @@
+double precision function cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho)
+
+! Compute the derivative of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: Re_w
+  double precision,intent(in)   :: Im_w
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Re_e(nBas)
+  double precision,intent(in)   :: Im_e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+
+! Local variables
+
+  integer                       :: i,a,m
+  double precision              :: num,eps,eta_tilde
+
+! Initialize 
+
+  cRGW_Re_dSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = Re_w - Re_e(i) + Om(m)
+      eta_tilde = eta - Im_w + Im_e(i)
+      num = 2d0*rho(p,i,m)**2
+      cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = Re_w - Re_e(a) - Om(m)
+      eta_tilde = eta + Im_w - Im_e(a)
+      num = 2d0*rho(p,a,m)**2
+      cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
+    end do
+  end do
+
+end function 
diff --git a/src/GW/cRGW_self_energy_diag.f90 b/src/GW/cRGW_self_energy_diag.f90
new file mode 100644
index 0000000..eb0717c
--- /dev/null
+++ b/src/GW/cRGW_self_energy_diag.f90
@@ -0,0 +1,101 @@
+subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z,CAP)
+
+! Compute diagonal of the correlation part of the self-energy and the renormalization factor
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: CAP(nOrb,nOrb)
+
+! Local variables
+
+  integer                       :: i,a,p,m
+  double precision              :: num,eps
+  double precision              :: eta_tilde
+  double precision              :: Re_DS(nBas)
+  double precision              :: Im_DS(nBas)
+  
+! Output variables
+
+  double precision,intent(out)  :: Re_Sig(nBas)
+  double precision,intent(out)  :: Im_Sig(nBas)
+  double precision,intent(out)  :: Re_Z(nBas)
+  double precision,intent(out)  :: Im_Z(nBas)
+  double precision,intent(out)  :: EcGM
+
+! Initialize 
+  Re_Sig(:) = 0d0
+  Im_Sig(:) = 0d0
+  Re_DS(:)   = 0d0
+  Im_DS(:)   = 0d0
+
+!----------------!
+! GW self-energy !
+!----------------!
+ 
+! Occupied part of the correlation self-energy
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do m=1,nS
+        eps = e(p) - e(i) + Om(m)
+        eta_tilde = eta  - CAP(p,p) + CAP(i,i) 
+        num = 2d0*rho(p,i,m)**2
+        Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2)
+        Im_Sig(p) = Im_Sig(p) + num*eta_tilde/(eps**2 + eta_tilde**2)
+        Re_DS(p)   = Re_DS(p)   - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
+        Im_DS(p)   = Im_DS(p)    - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
+
+      end do
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      do m=1,nS
+
+        eps = e(p) - e(a) - Om(m)
+        eta_tilde = eta  + CAP(p,p) - CAP(a,a)
+        num = 2d0*rho(p,a,m)**2
+        Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2)
+        Im_Sig(p) = Im_Sig(p) - num*eta_tilde/(eps**2 + eta_tilde**2)
+        Re_DS(p)   = Re_DS(p)   - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2
+        Im_DS(p)   = Im_DS(p)    + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2
+      end do
+    end do
+  end do
+
+! Galitskii-Migdal correlation energy
+! MAYBE MODIFY THIS  TERM
+  EcGM = 0d0
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      do m=1,nS
+
+        eps = e(a) - e(i) + Om(m)
+        num = 4d0*rho(a,i,m)**2
+        EcGM = EcGM - num*eps/(eps**2 + eta**2)
+
+      end do
+    end do
+  end do
+
+! Compute renormalization factor from derivative 
+  Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
+  Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
+
+end subroutine 
diff --git a/src/GW/print_cRG0W0.f90 b/src/GW/print_cRG0W0.f90
new file mode 100644
index 0000000..0121e08
--- /dev/null
+++ b/src/GW/print_cRG0W0.f90
@@ -0,0 +1,64 @@
+subroutine print_cRG0W0(nBas,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP)
+
+! Print one-electron energies and other stuff for G0W0
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas,nO
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ERHF
+  double precision,intent(in)        :: EcRPA
+  double precision,intent(in)        :: EcGM
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: Re_SigC(nBas)
+  double precision,intent(in)        :: Im_SigC(nBas)
+  double precision,intent(in)        :: Re_Z(nBas)
+  double precision,intent(in)        :: Im_Z(nBas)
+  double precision,intent(in)        :: Re_eGW(nBas)
+  double precision,intent(in)        :: Im_eGW(nBas)
+  double precision,intent(in)        :: CAP(nBas,nBas)
+
+  integer                            :: p
+  double precision                   :: eHOMO,eLUMO,Gap
+
+! HOMO and LUMO
+
+  eHOMO = maxval(Re_eGW(1:nO))
+  eLUMO = minval(Re_eGW(nO+1:nBas))
+  Gap = eLUMO-eHOMO
+
+! Dump results
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)' G0W0@RHF calculation '
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(1X,A11,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Re(Sig_GW) (eV)','|','Re(Z)','|','Re(e_GW) (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+  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,F15.6,1X,A1,1X)') &
+    '|',p,'|',eHF(p)*HaToeV,'|',Re_SigC(p)*HaToeV,'|',Re_Z(p),'|',Re_eGW(p)*HaToeV,'|'
+  end do
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(1X,A11,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','CAP (eV)','|','Im(Sig_GW) (eV)','|','Im(Z)','|','Im(e_GW) (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+  do p=1,nBas
+    write(*,'(1X,A11,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',p,'|',CAP(p,p)*HaToeV,'|',Im_SigC(p)*HaToeV,'|',Im_Z(p),'|',Im_eGW(p)*HaToeV,'|'
+  end do
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A60,F15.6,A3)') 'G0W0@RHF HOMO      energy = ',eHOMO*HaToeV,' eV'
+  write(*,'(2X,A60,F15.6,A3)') 'G0W0@RHF LUMO      energy = ',eLUMO*HaToeV,' eV'
+  write(*,'(2X,A60,F15.6,A3)') 'G0W0@RHF HOMO-LUMO gap    = ',Gap*HaToeV,' eV'
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A60,F15.6,A3)') 'phRPA@G0W0@RHF total       energy = ',ENuc + ERHF + EcRPA,' au'
+  write(*,'(2X,A60,F15.6,A3)') 'phRPA@G0W0@RHF correlation energy = ',EcRPA,' au'
+  write(*,'(2X,A60,F15.6,A3)') '   GM@G0W0@RHF total       energy = ',ENuc + ERHF + EcGM,' au'
+  write(*,'(2X,A60,F15.6,A3)') '   GM@G0W0@RHF correlation energy = ',EcGM,' au'
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+end subroutine 
diff --git a/src/cRGW_Im_dSigC.f90 b/src/cRGW_Im_dSigC.f90
new file mode 100644
index 0000000..e69de29