diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90
index 49ac42d..e6f0d1a 100644
--- a/src/GW/GGW_ppBSE.f90
+++ b/src/GW/GGW_ppBSE.f90
@@ -115,6 +115,11 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
       call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                           Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
+  !----------------!
+  ! Upfolded ppBSE !
+  !----------------!
+
+  call GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
 
   deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
 
diff --git a/src/GW/GGW_ppBSE_dynamic_kernel_D.f90 b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90
index 73dbe43..ba1a670 100644
--- a/src/GW/GGW_ppBSE_dynamic_kernel_D.f90
+++ b/src/GW/GGW_ppBSE_dynamic_kernel_D.f90
@@ -53,25 +53,25 @@ subroutine GGW_ppBSE_dynamic_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,eGW,Om,
   
           do m=1,nS
 
-            dem = OmBSE - eGW(k) + Om(m) - eGW(i)
+            dem = OmBSE - eGW(k) + Om(m) - eGW(j)
             num = rho(i,k,m)*rho(j,l,m)
 
-            KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
-            ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+            KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+            ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
 
-            dem = OmBSE - eGW(k) + Om(m) - eGW(j)
+            dem = OmBSE - eGW(k) + Om(m) - eGW(i)
             num = rho(j,k,m)*rho(i,l,m)
 
             KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
             ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
 
-            dem = OmBSE - eGW(l) + Om(m) - eGW(j)
+            dem = OmBSE - eGW(l) + Om(m) - eGW(i)
             num = rho(i,k,m)*rho(j,l,m)
 
-            KD_dyn(ij,kl) = KD_dyn(ij,kl) - num*dem/(dem**2 + eta**2)
-            ZD_dyn(ij,kl) = ZD_dyn(ij,kl) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+            KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
+            ZD_dyn(ij,kl) = ZD_dyn(ij,kl) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
 
-            dem = OmBSE - eGW(l) + Om(m) - eGW(i)
+            dem = OmBSE - eGW(l) + Om(m) - eGW(j)
             num = rho(j,k,m)*rho(i,l,m)
 
             KD_dyn(ij,kl) = KD_dyn(ij,kl) + num*dem/(dem**2 + eta**2)
diff --git a/src/GW/GGW_ppBSE_static_kernel_D.f90 b/src/GW/GGW_ppBSE_static_kernel_D.f90
index 37c0672..57e5f10 100644
--- a/src/GW/GGW_ppBSE_static_kernel_D.f90
+++ b/src/GW/GGW_ppBSE_static_kernel_D.f90
@@ -55,7 +55,7 @@ subroutine GGW_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,r
                       + rho(i,l,m)*rho(j,k,m)*Om(m)/eps
           end do
 
-          KD(ij,kl) = 2d0*lambda*chi
+          KD(ij,kl) = lambda*chi
 
         end do
       end do
diff --git a/src/GW/RGW_ppBSE_upfolded.f90 b/src/GW/RGW_ppBSE_upfolded.f90
index 3555425..d5e8deb 100644
--- a/src/GW/RGW_ppBSE_upfolded.f90
+++ b/src/GW/RGW_ppBSE_upfolded.f90
@@ -76,7 +76,7 @@ subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
   n3h1p = n2h*n1h1p
   n3p1h = n2p*n1h1p
 
-  nH = n1h1p + 4*n3h1p 
+  nH = n2h + 4*n3h1p 
 
 ! Memory allocation
 
@@ -220,8 +220,6 @@ subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
 ! Diagonalize supermatrix !
 !-------------------------!
 
-  call matout(nH,nH,H)
-
   call diagonalize_general_matrix(nH,H,OmBSE,X)
 
   do s=1,nH