diff --git a/src/GW/GW_ImSigC.f90 b/src/GW/GW_ImSigC.f90
new file mode 100644
index 0000000..414acb6
--- /dev/null
+++ b/src/GW/GW_ImSigC.f90
@@ -0,0 +1,52 @@
+double precision function GW_ImSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,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)   :: 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)   :: 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
+
+! Initialize 
+
+  GW_ImSigC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do m=1,nS
+      eps = w - e(i) + Om(m)
+      num = 2d0*rho(p,i,m)**2
+      GW_ImSigC = GW_ImSigC + num*eta/(eps**2 + eta**2)
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do m=1,nS
+      eps = w - e(a) - Om(m)
+      num = 2d0*rho(p,a,m)**2
+      GW_ImSigC = GW_ImSigC - num*eta/(eps**2 + eta**2)
+    end do
+  end do
+
+end function 
diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90
index 3f72071..2ec503c 100644
--- a/src/GW/GW_phBSE.f90
+++ b/src/GW/GW_phBSE.f90
@@ -120,6 +120,12 @@ subroutine GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,
     call print_excitation_energies('phBSE@GW@RHF','singlet',nS,OmBSE)
     call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
+    !--------------------!
+    ! Cumulant expansion !
+    !--------------------!
+
+    call RGWC(.false.,nBas,nC,nO,nR,nS,OmBSE,rho_RPA,eGW)
+
     !----------------------------------------------------!
     ! Compute the dynamical screening at the phBSE level !
     !----------------------------------------------------!