diff --git a/input/methods.default b/input/methods.default
index 685092a..addaf23 100644
--- a/input/methods.default
+++ b/input/methods.default
@@ -12,8 +12,8 @@
   F     F      F     F 
 # G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
   F    F     F     F    F    F
-# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
-  F    F    F    F        F      F
+# G0W0 evGW qsGW ufG0W0 ufGW
+  F    F    F    F      F
 # G0T0pp evGTpp qsGTpp ufG0T0pp
   F      F      F      F
 # G0T0eh evGTeh qsGTeh
diff --git a/input/options.default b/input/options.default
index 8a86013..0c6c2db 100644
--- a/input/options.default
+++ b/input/options.default
@@ -4,7 +4,7 @@
       F
 # CC: maxSCF thresh   DIIS
       64     0.00001  5
-# spin: TDA singlet triplet
+# LR: TDA singlet triplet
         F   T       T
 # GF: maxSCF thresh  DIIS lin eta renorm reg
       256    0.00001 5    F   0.0 0      F
diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90
index 20451d5..038df13 100644
--- a/src/GW/GG0W0.f90
+++ b/src/GW/GG0W0.f90
@@ -1,5 +1,5 @@
 subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+                 linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 ! Perform G0W0 calculation
 
   implicit none
@@ -22,7 +22,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
   logical,intent(in)            :: dTDA
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -77,11 +77,13 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Memory allocation
@@ -110,9 +112,11 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Compute GW self-energy !
 !------------------------!
 
-  if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
-
-  call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+  if(doSRG) then 
+    call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+  else
+    call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+  end if
 
 !-----------------------------------!
 ! Solve the quasi-particle equation !
diff --git a/src/GW/GGW.f90 b/src/GW/GGW.f90
index 52e8ab7..adacd18 100644
--- a/src/GW/GGW.f90
+++ b/src/GW/GGW.f90
@@ -1,5 +1,5 @@
-subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,  & 
-               TDA_W,TDA,dBSE,dTDA,linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc, & 
+subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & 
+               TDA_W,TDA,dBSE,dTDA,linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,     & 
                ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,eHF)
 
 ! GW module
@@ -30,7 +30,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -71,7 +71,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
     
     call wall_time(start_GW)
     call GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
-              linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+              linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
     call wall_time(end_GW)
   
     t_GW = end_GW - start_GW
@@ -88,7 +88,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
 
     call wall_time(start_GW)
     call evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
-               linearize,eta,regularize,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+               linearize,eta,doSRG,nBas2,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
     call wall_time(end_GW)
 
     t_GW = end_GW - start_GW
@@ -105,7 +105,7 @@ subroutine GGW(dotest,doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT,exchan
 
     call wall_time(start_GW)
     call qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-               eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI,                  & 
+               eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,S,X,T,V,Hc,ERI_AO,ERI,                       & 
                dipole_int_AO,dipole_int,PHF,cHF,eHF)
     call wall_time(end_GW)
 
diff --git a/src/GW/GGW_SRG_self_energy.f90 b/src/GW/GGW_SRG_self_energy.f90
new file mode 100644
index 0000000..a0c7c70
--- /dev/null
+++ b/src/GW/GGW_SRG_self_energy.f90
@@ -0,0 +1,144 @@
+subroutine GGW_SRG_self_energy(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+
+! Compute correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  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(nOrb)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  integer                       :: p,q
+  integer                       :: m
+  double precision              :: Dpim,Dqim,Dpam,Dqam,Diam
+  double precision              :: s
+  
+! Output variables
+
+  double precision,intent(out)  :: EcGM
+  double precision,intent(out)  :: SigC(nOrb,nOrb)
+  double precision,intent(out)  :: Z(nOrb)
+
+!--------------------!
+! SRG flow parameter !
+!--------------------!
+
+  s = 500d0
+
+!--------------------!
+! SRG-GW self-energy !
+!--------------------!
+
+  SigC(:,:) = 0d0
+
+  ! Occupied part of the correlation self-energy
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) &
+  !$OMP PRIVATE(m,i,q,p,Dpim,Dqim) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO 
+  do q=nC+1,nOrb-nR
+    do p=nC+1,nOrb-nR
+      do m=1,nS
+        do i=nC+1,nO
+
+          Dpim = e(p) - e(i) + Om(m)
+          Dqim = e(q) - e(i) + Om(m)
+          SigC(p,q) = SigC(p,q)                  &
+                    + rho(p,i,m)*rho(q,i,m)*     &
+                      (1d0-exp(-s*Dpim*Dpim)*exp(-s*Dqim*Dqim))*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
+
+        end do
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Virtual part of the correlation self-energy
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) &
+  !$OMP PRIVATE(m,a,q,p,Dpam,Dqam) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do q=nC+1,nOrb-nR
+    do p=nC+1,nOrb-nR
+      do m=1,nS
+        do a=nO+1,nOrb-nR
+ 
+           Dpam = e(p) - e(a) - Om(m)
+           Dqam = e(q) - e(a) - Om(m)
+           SigC(p,q) = SigC(p,q)                  & 
+                     + rho(p,a,m)*rho(q,a,m)*     &
+                       (1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam))*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
+ 
+        end do
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+ 
+!------------------------!
+! Renormalization factor !
+!------------------------!
+
+  Z(:)  = 0d0
+
+  ! Occupied part of the renormlization factor
+
+  do p=nC+1,nOrb-nR
+    do i=nC+1,nO
+      do m=1,nS
+
+        Dpim = e(p) - e(i) + Om(m)
+        Z(p) = Z(p) - rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2
+
+      end do
+    end do
+  end do
+
+  ! Virtual part of the renormlization factor
+
+  do p=nC+1,nOrb-nR
+     do a=nO+1,nOrb-nR
+        do m=1,nS
+           Dpam = e(p) - e(a) - Om(m)
+           Z(p) = Z(p)  - rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2
+        end do
+     end do
+  end do
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+!-------------------------------------!
+! Galitskii-Migdal correlation energy !
+!-------------------------------------!
+
+  EcGM = 0d0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      do m=1,nS
+
+        Diam = e(a) - e(i) + Om(m)
+        EcGM = EcGM - rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*s*Diam*Diam))/Diam 
+
+      end do
+    end do
+  end do
+
+end subroutine 
diff --git a/src/GW/GGW_SRG_self_energy_diag.f90 b/src/GW/GGW_SRG_self_energy_diag.f90
new file mode 100644
index 0000000..c63b2b3
--- /dev/null
+++ b/src/GW/GGW_SRG_self_energy_diag.f90
@@ -0,0 +1,134 @@
+subroutine GGW_SRG_self_energy_diag(nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z)
+
+! Compute correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  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(nOrb)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  integer                       :: p
+  integer                       :: m
+  double precision              :: Dpim,Dpam,Diam
+  double precision              :: s
+  
+! Output variables
+
+  double precision,intent(out)  :: EcGM
+  double precision,intent(out)  :: SigC(nOrb)
+  double precision,intent(out)  :: Z(nOrb)
+
+!--------------------!
+! SRG flow parameter !
+!--------------------!
+
+  s = 500d0
+
+!--------------------!
+! SRG-GW self-energy !
+!--------------------!
+
+  SigC(:) = 0d0
+
+  ! Occupied part of the correlation self-energy
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(SigC,rho,s,nS,nC,nO,nOrb,nR,e,Om) &
+  !$OMP PRIVATE(m,i,p,Dpim) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO 
+  do p=nC+1,nOrb-nR
+    do m=1,nS
+      do i=nC+1,nO
+
+        Dpim = e(p) - e(i) + Om(m)
+        SigC(p) = SigC(p) + rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim
+
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  ! Virtual part of the correlation self-energy
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nOrb,e,Om) &
+  !$OMP PRIVATE(m,a,p,Dpam) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do p=nC+1,nOrb-nR
+    do m=1,nS
+      do a=nO+1,nOrb-nR
+ 
+         Dpam = e(p) - e(a) - Om(m)
+         SigC(p) = SigC(p) + rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam
+ 
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+ 
+!------------------------!
+! Renormalization factor !
+!------------------------!
+
+  Z(:)  = 0d0
+
+  ! Occupied part of the renormlization factor
+
+  do p=nC+1,nOrb-nR
+    do i=nC+1,nO
+      do m=1,nS
+
+        Dpim = e(p) - e(i) + Om(m)
+        Z(p) = Z(p) - rho(p,i,m)**2*(1d0-exp(-2d0*s*Dpim*Dpim))/Dpim**2
+
+      end do
+    end do
+  end do
+
+  ! Virtual part of the renormlization factor
+
+  do p=nC+1,nOrb-nR
+     do a=nO+1,nOrb-nR
+        do m=1,nS
+           Dpam = e(p) - e(a) - Om(m)
+           Z(p) = Z(p)  - rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2
+        end do
+     end do
+  end do
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+!-------------------------------------!
+! Galitskii-Migdal correlation energy !
+!-------------------------------------!
+
+  EcGM = 0d0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      do m=1,nS
+
+        Diam = e(a) - e(i) + Om(m)
+        EcGM = EcGM - rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*s*Diam*Diam))/Diam 
+
+      end do
+    end do
+  end do
+
+end subroutine 
diff --git a/src/GW/GGW_phBSE.f90 b/src/GW/GGW_phBSE.f90
index 5e15087..bee353e 100644
--- a/src/GW/GGW_phBSE.f90
+++ b/src/GW/GGW_phBSE.f90
@@ -73,9 +73,21 @@ subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,di
   call GGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
   call GGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
 
-  EcBSE = 0d0
 
-  ! Compute BSE excitation energies
+!-----!
+! TDA !
+!-----!
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
+    write(*,*)
+  end if
+
+!---------------------------------!
+! Compute BSE excitation energies !
+!---------------------------------!
+
+  EcBSE = 0d0
 
                call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
   if(.not.TDA) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
diff --git a/src/GW/GGW_self_energy.f90 b/src/GW/GGW_self_energy.f90
index 3d2bc73..135068d 100644
--- a/src/GW/GGW_self_energy.f90
+++ b/src/GW/GGW_self_energy.f90
@@ -1,4 +1,4 @@
-subroutine GGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
+subroutine GGW_self_energy(eta,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 
 ! Compute correlation part of the self-energy and the renormalization factor 
 
@@ -8,15 +8,15 @@ subroutine GGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 ! 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)   :: e(nOrb)
   double precision,intent(in)   :: Om(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
 
 ! Local variables
 
@@ -27,71 +27,118 @@ subroutine GGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
 ! Output variables
 
   double precision,intent(out)  :: EcGM
-  double precision,intent(out)  :: Sig(nBas,nBas)
-  double precision,intent(out)  :: Z(nBas)
-
-! Initialize 
-
-  Sig(:,:) = 0d0
-  Z(:)     = 0d0
+  double precision,intent(out)  :: Sig(nOrb,nOrb)
+  double precision,intent(out)  :: Z(nOrb)
 
 !----------------!
 ! GW self-energy !
 !----------------!
 
+  Sig(:,:) = 0d0
+
 ! Occupied part of the correlation self-energy
 
-!$OMP PARALLEL &
-!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
-!$OMP PRIVATE(m,i,q,p,eps,num) &
-!$OMP DEFAULT(NONE)
-!$OMP DO
-  do q=nC+1,nBas-nR
-     do p=nC+1,nBas-nR
-        do m=1,nS
-           do i=nC+1,nO
+  !$OMP PARALLEL &
+  !$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
+  !$OMP PRIVATE(m,i,q,p,eps,num) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO
+  do q=nC+1,nOrb-nR
+    do p=nC+1,nOrb-nR
+      do m=1,nS
+        do i=nC+1,nO
  
-              eps = e(p) - e(i) + Om(m)
-              num = rho(p,i,m)*rho(q,i,m)
-              Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
-              if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+          eps = e(p) - e(i) + Om(m)
+          num = rho(p,i,m)*rho(q,i,m)
+          Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
  
-           end do
         end do
-     end do
+      end do
+    end do
   end do
   !$OMP END DO
   !$OMP END PARALLEL
 
 ! Virtual part of the correlation self-energy
 
+  !$OMP PARALLEL &
+  !$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
+  !$OMP PRIVATE(m,a,q,p,eps,num) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO  
+  do q=nC+1,nOrb-nR
+    do p=nC+1,nOrb-nR
+      do m=1,nS
+        do a=nO+1,nOrb-nR
+ 
+          eps = e(p) - e(a) - Om(m)
+          num = rho(p,a,m)*rho(q,a,m)
+          Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
+ 
+        end do
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+!------------------------!
+! Renormalization factor !
+!------------------------!
+
+  Z(:)     = 0d0
+
+! Occupied part of the renormalization factor
+
 !$OMP PARALLEL &
-!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
-!$OMP PRIVATE(m,a,q,p,eps,num) &
+!$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
+!$OMP PRIVATE(m,i,q,p,eps,num) &
 !$OMP DEFAULT(NONE)
-!$OMP DO  
-  do q=nC+1,nBas-nR
-     do p=nC+1,nBas-nR
-        do m=1,nS
-           do a=nO+1,nBas-nR
+!$OMP DO
+  do p=nC+1,nOrb-nR
+     do m=1,nS
+        do i=nC+1,nO
  
-              eps = e(p) - e(a) - Om(m)
-              num = rho(p,a,m)*rho(q,a,m)
-              Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2)
-              if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+           eps = e(p) - e(i) + Om(m)
+           num = rho(p,i,m)*rho(q,i,m)
+           Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
  
-           end do
         end do
      end do
   end do
-!$OMP END DO
-!$OMP END PARALLEL
+  !$OMP END DO
+  !$OMP END PARALLEL
 
-! Galitskii-Migdal correlation energy
+! Virtual part of the renormalization factor
+
+  !$OMP PARALLEL &
+  !$OMP SHARED(Sig,Z,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
+  !$OMP PRIVATE(m,a,q,p,eps,num) &
+  !$OMP DEFAULT(NONE)
+  !$OMP DO  
+  do p=nC+1,nOrb-nR
+    do m=1,nS
+      do a=nO+1,nOrb-nR
+ 
+        eps = e(p) - e(a) - Om(m)
+        num = rho(p,a,m)*rho(q,a,m)
+        Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+ 
+      end do
+    end do
+  end do
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  Z(:) = 1d0/(1d0 - Z(:))
+
+!-------------------------------------!
+! Galitskii-Migdal correlation energy !
+!-------------------------------------!
 
   EcGM = 0d0
   do m=1,nS
-    do a=nO+1,nBas-nR
+    do a=nO+1,nOrb-nR
       do i=nC+1,nO
 
         eps = e(a) - e(i) + Om(m)
@@ -102,8 +149,4 @@ subroutine GGW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
     end do
   end do
 
-! Compute renormalization factor from derivative 
-
-  Z(:) = 1d0/(1d0 - Z(:))
-
 end subroutine 
diff --git a/src/GW/GW_regularization.f90 b/src/GW/GW_regularization.f90
deleted file mode 100644
index 33dead2..0000000
--- a/src/GW/GW_regularization.f90
+++ /dev/null
@@ -1,53 +0,0 @@
-subroutine GW_regularization(nBas,nC,nO,nV,nR,nS,e,Om,rho)
-
-! Regularize GW excitation densities via SRG
-
-  implicit none
-
-! 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
-  integer,intent(in)            :: e(nBas)
-  integer,intent(in)            :: Om(nS)
-
-! Local variables
-
-  integer                       :: p,i,a,m
-  double precision              :: s
-  double precision              :: kappa
-  double precision              :: Dpim,Dpam
-
-! Output variables
-
-  double precision,intent(inout):: rho(nBas,nBas,nS)
-
-! SRG flow parameter 
-
-  s = 500d0
-
-! Regularize excitation densities
-
-  do p=nC+1,nBas-nR
-    do m=1,nS
-
-      do i=nC+1,nO
-        Dpim = e(p) - e(i) + Om(m)
-        kappa = 1d0 - exp(-Dpim*Dpim*s)
-        rho(p,i,m) = kappa*rho(p,i,m)
-      end do
-
-      do a=nO+1,nBas-nR
-        Dpam = e(p) - e(a) - Om(m)
-        kappa = 1d0 - exp(-Dpam*Dpam*s)
-        rho(p,a,m) = kappa*rho(p,a,m)
-      end do
-
-    end do
-  end do
-
-end subroutine 
diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90
index 897969b..ae1a340 100644
--- a/src/GW/RGW.f90
+++ b/src/GW/RGW.f90
@@ -1,6 +1,6 @@
-subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
-               exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,   &
-               linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,           &
+subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT,         &
+               exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
+               linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,              &
                S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Restricted GW module
@@ -17,7 +17,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
   logical,intent(in)            :: doqsGW
   logical,intent(in)            :: doufG0W0
   logical,intent(in)            :: doufGW
-  logical,intent(in)            :: doSRGqsGW
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -36,7 +35,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -77,7 +76,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
     
     call wall_time(start_GW)
     call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
-               linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
+               linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
     call wall_time(end_GW)
   
     t_GW = end_GW - start_GW
@@ -94,7 +93,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
     call wall_time(start_GW)
     call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
-               singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
+               singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
     call wall_time(end_GW)
 
     t_GW = end_GW - start_GW
@@ -111,7 +110,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
     call wall_time(start_GW)
     call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, &
-               TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,    &
+               TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,nNuc,ZNuc,rNuc,         &
                ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,                  &
                dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
     call wall_time(end_GW)
@@ -122,26 +121,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
   end if
 
-!------------------------------------------------------------------------
-! Perform SRG-qsGW calculation
-!------------------------------------------------------------------------
-
-  if(doSRGqsGW) then 
-
-    call wall_time(start_GW)
-    call SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
-                   dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,    &
-                   nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,                &
-                   ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,   &
-                   PHF,cHF,eHF)
-    call wall_time(end_GW)
-
-    t_GW = end_GW - start_GW
-    write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
-    write(*,*)
-
-  end if
-
 !------------------------------------------------------------------------
 ! Perform ufG0W0 calculatiom
 !------------------------------------------------------------------------
diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90
deleted file mode 100644
index 04404cc..0000000
--- a/src/GW/SRG_qsUGW.f90
+++ /dev/null
@@ -1,419 +0,0 @@
-subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
-                     eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
-                     dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
-
-! Perform a quasiparticle self-consistent GW calculation
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  logical,intent(in)            :: dotest
-
-  integer,intent(in)            :: maxSCF
-  integer,intent(in)            :: max_diis
-  double precision,intent(in)   :: thresh
-  logical,intent(in)            :: doACFDT
-  logical,intent(in)            :: exchange_kernel
-  logical,intent(in)            :: doXBS
-  logical,intent(in)            :: BSE
-  logical,intent(in)            :: TDA_W
-  logical,intent(in)            :: TDA
-  logical,intent(in)            :: dBSE
-  logical,intent(in)            :: dTDA
-  logical,intent(in)            :: spin_conserved
-  logical,intent(in)            :: spin_flip
-  double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
-
-  integer,intent(in)            :: nNuc
-  double precision,intent(in)   :: ZNuc(nNuc)
-  double precision,intent(in)   :: rNuc(nNuc,ncart)
-  double precision,intent(in)   :: ENuc
-
-  integer,intent(in)            :: nBas
-  integer,intent(in)            :: nC(nspin)
-  integer,intent(in)            :: nO(nspin)
-  integer,intent(in)            :: nV(nspin)
-  integer,intent(in)            :: nR(nspin)
-  integer,intent(in)            :: nS(nspin)
-
-  double precision,intent(in)   :: EUHF
-  double precision,intent(in)   :: eHF(nBas,nspin)
-  double precision,intent(in)   :: cHF(nBas,nBas,nspin)
-  double precision,intent(in)   :: PHF(nBas,nBas,nspin)
-  double precision,intent(in)   :: S(nBas,nBas)
-  double precision,intent(in)   :: T(nBas,nBas)
-  double precision,intent(in)   :: V(nBas,nBas)
-  double precision,intent(in)   :: Hc(nBas,nBas)
-  double precision,intent(in)   :: X(nBas,nBas)
-  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
-  double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas)
-  double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas)
-  double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: dipole_int_AO(nBas,nBas,ncart)
-  double precision,intent(inout):: dipole_int_aa(nBas,nBas,ncart)
-  double precision,intent(inout):: dipole_int_bb(nBas,nBas,ncart)
-
-! Local variables
-
-  logical                       :: dRPA
-  integer                       :: nSCF
-  integer                       :: nBasSq
-  integer                       :: ispin
-  integer                       :: ixyz
-  integer                       :: is
-  integer                       :: n_diis
-  integer                       :: nSa,nSb,nSt
-  double precision              :: dipole(ncart)
-
-  double precision              :: ET(nspin)
-  double precision              :: EV(nspin)
-  double precision              :: EJ(nsp)
-  double precision              :: EK(nspin)
-  double precision              :: EcRPA(nspin)
-  double precision              :: EcGM(nspin)
-  double precision              :: EqsGW
-  double precision              :: EcBSE(nspin)
-  double precision              :: Conv
-  double precision              :: rcond(nspin)
-  double precision,external     :: trace_matrix
-  double precision,allocatable  :: err_diis(:,:,:)
-  double precision,allocatable  :: F_diis(:,:,:)
-
-  double precision,allocatable  :: Aph(:,:)
-  double precision,allocatable  :: Bph(:,:)
-  double precision,allocatable  :: Om(:)
-  double precision,allocatable  :: XpY(:,:)
-  double precision,allocatable  :: XmY(:,:)
-  double precision,allocatable  :: rho(:,:,:,:)
-  double precision,allocatable  :: c(:,:,:)
-  double precision,allocatable  :: cp(:,:,:)
-  double precision,allocatable  :: eGW(:,:)
-  double precision,allocatable  :: P(:,:,:)
-  double precision,allocatable  :: F(:,:,:)
-  double precision,allocatable  :: Fp(:,:,:)
-  double precision,allocatable  :: J(:,:,:)
-  double precision,allocatable  :: K(:,:,:)
-  double precision,allocatable  :: SigC(:,:,:)
-  double precision,allocatable  :: SigCp(:,:,:)
-  double precision,allocatable  :: Z(:,:)
-  double precision,allocatable  :: err(:,:,:)
-
-! Hello world
-
-  write(*,*)
-  write(*,*)'*************************************'
-  write(*,*)'* Unrestricted SRG-qsGW Calculation *'
-  write(*,*)'*************************************'
-  write(*,*)
-
-! Warning 
-
-  write(*,*) '!! ERIs in MO basis will be overwritten in qsUGW !!'
-  write(*,*)
-
-! Stuff 
-
-  nBasSq = nBas*nBas
-  dRPA = .true.
-
-! TDA for W
-
-  if(TDA_W) then 
-    write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
-    write(*,*)
-  end if
-
-! TDA 
-
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
-    write(*,*)
-  end if
-
-! Memory allocation
-
-  nSa = nS(1)
-  nSb = nS(2)
-  nSt = nSa + nSb
-
-  allocate(Aph(nSt,nSt),Bph(nSt,nSt),eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin), &
-           F(nBas,nBas,nspin),Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),                        & 
-           SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),Z(nBas,nspin),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt),        &  
-           rho(nBas,nBas,nSt,nspin),err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
-
-! Initialization
-  
-  nSCF            = -1
-  n_diis          = 0
-  ispin           = 1
-  Conv            = 1d0
-  P(:,:,:)        = PHF(:,:,:)
-  eGW(:,:)        = eHF(:,:)
-  c(:,:,:)        = cHF(:,:,:)
-  F_diis(:,:,:)   = 0d0
-  err_diis(:,:,:) = 0d0
-  rcond(:)        = 0d0
-
-!------------------------------------------------------------------------
-! Main loop
-!------------------------------------------------------------------------
-
-  do while(Conv > thresh .and. nSCF < maxSCF)
-
-    ! Increment
-
-    nSCF = nSCF + 1
-
-    ! Buid Hartree matrix
-
-    do is=1,nspin
-      call Hartree_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is))
-    end do
-
-    ! Compute exchange part of the self-energy 
-
-    do is=1,nspin
-      call exchange_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),K(:,:,is))
-    end do
-
-    !--------------------------------------------------
-    ! AO to MO transformation of two-electron integrals
-    !--------------------------------------------------
- 
-    do ixyz=1,ncart
-        call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
-        call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
-    end do
-
-    ! 4-index transform for (aa|aa) block
-
-    call AOtoMO_ERI_UHF(1,1,nBas,c,ERI_AO,ERI_aaaa)
-
-    ! 4-index transform for (aa|bb) block
-
-    call AOtoMO_ERI_UHF(1,2,nBas,c,ERI_AO,ERI_aabb)
-
-    ! 4-index transform for (bb|bb) block
-
-    call AOtoMO_ERI_UHF(2,2,nBas,c,ERI_AO,ERI_bbbb)
-
-    ! Compute linear response
-
-    call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph)
-    if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
-    
-    call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
-
-    !----------------------!
-    ! Excitation densities !
-    !----------------------!
-
-    call UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
-
-    !------------------------------------------------!
-    ! Compute self-energy and renormalization factor !
-    !------------------------------------------------!
-
-    call UGW_SRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z)
-
-    ! Make correlation self-energy Hermitian and transform it back to AO basis
-   
-    do is=1,nspin
-      call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
-    end do
- 
-    ! Solve the quasi-particle equation
-
-    do is=1,nspin
-      F(:,:,is) = Hc(:,:) + J(:,:,is) + J(:,:,mod(is,2)+1) + K(:,:,is) + SigCp(:,:,is)
-    end do
-
-   ! Check convergence 
-
-    do is=1,nspin
-      err(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is))
-    end do
-
-    if(nSCF > 1) Conv = maxval(abs(err))
-
-    ! DIIS extrapolation
-
-    if(max_diis > 1) then
-
-      n_diis = min(n_diis+1,max_diis)
-      do is=1,nspin
-        if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,is), &
-                                               F_diis(:,1:n_diis,is),err(:,:,is),F(:,:,is))
-      end do
-
-    end if
-
-    ! Transform Fock matrix in orthogonal basis
-
-    do is=1,nspin
-      Fp(:,:,is) = matmul(transpose(X(:,:)),matmul(F(:,:,is),X(:,:)))
-    end do
-
-    ! Diagonalize Fock matrix to get eigenvectors and eigenvalues
-
-    cp(:,:,:) = Fp(:,:,:)
-    do is=1,nspin
-      call diagonalize_matrix(nBas,cp(:,:,is),eGW(:,is))
-    end do
-
-    ! Back-transform eigenvectors in non-orthogonal basis
-
-    do is=1,nspin
-      c(:,:,is) = matmul(X(:,:),cp(:,:,is))
-    end do
-
-    ! Back-transform self-energy
-
-    do is=1,nspin
-      call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
-    end do
-
-    ! Compute density matrix 
-
-    do is=1,nspin
-      P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is)))
-    end do
-
-    !------------------------------------------------------------------------
-    !   Compute total energy
-    !------------------------------------------------------------------------
-
-    ! Kinetic energy
-
-    do is=1,nspin
-      ET(is) = trace_matrix(nBas,matmul(P(:,:,is),T(:,:)))
-    end do
-
-    ! Potential energy
-
-    do is=1,nspin
-      EV(is) = trace_matrix(nBas,matmul(P(:,:,is),V(:,:)))
-    end do
-
-    ! Hartree energy
-
-    EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
-    EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) &
-          + 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1)))
-    EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
-
-    ! Exchange energy
-
-    do is=1,nspin
-      EK(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is)))
-    end do
-
-    ! Total energy
-
-    EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(EK(:))
-
-    !------------------------------------------------------------------------
-    ! Print results
-    !------------------------------------------------------------------------
-
-    call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,EK,EcGM,EcRPA(ispin),EqsGW,SigCp,Z,dipole)
-
-  end do
-!------------------------------------------------------------------------
-! End main loop
-!------------------------------------------------------------------------
-
-! Did it actually converge?
-
-  if(nSCF == maxSCF) then
-
-    write(*,*)
-    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
-    write(*,*)'                 Convergence failed                 '
-    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
-    write(*,*)
-
-    stop
-
-  end if
-
-! Deallocate memory
-
-  deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis)
-
-! Perform BSE calculation
-
-  if(BSE) then
-
-    call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
-                   S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE)
-
-    if(exchange_kernel) then
-
-      EcBSE(1) = 0.5d0*EcBSE(1)
-      EcBSE(2) = 0.5d0*EcBSE(2)
-
-    else
-
-      EcBSE(2) = 0.0d0
-
-    end if
-
-    write(*,*)
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-flip)      = ',EcBSE(2),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy                  = ',sum(EcBSE),' au'
-    write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF total       energy                  = ',ENuc + EqsGW + sum(EcBSE),' au'
-    write(*,*)'-------------------------------------------------------------------------------'
-    write(*,*)
-
-!   Compute the BSE correlation energy via the adiabatic connection 
-
-    if(doACFDT) then
-
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*) ' Adiabatic connection version of BSE@qsUGW correlation energy '
-      write(*,*) '--------------------------------------------------------------'
-      write(*,*)
-
-      if(doXBS) then
-
-        write(*,*) '*** scaled screening version (XBS) ***'
-        write(*,*)
-
-      end if
-
-      call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, &
-                       eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA)
-
-      write(*,*)
-      write(*,*)'-------------------------------------------------------------------------------'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcRPA(1),' au'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy (spin-flip)      = ',EcRPA(2),' au'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy                  = ',sum(EcRPA),' au'
-      write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF total       energy                  = ',ENuc + EqsGW + sum(EcRPA),' au'
-      write(*,*)'-------------------------------------------------------------------------------'
-      write(*,*)
-
-    end if
-
-  end if
-
-! Testing zone
-  
-  if(dotest) then
-  
-    call dump_test_value('U','qsGW correlation energy',EcRPA)
-    call dump_test_value('U','qsGW HOMOa energy',eGW(nO(1),1))
-    call dump_test_value('U','qsGW LUMOa energy',eGW(nO(1)+1,1))
-    call dump_test_value('U','qsGW HOMOa energy',eGW(nO(2),2))
-    call dump_test_value('U','qsGW LUMOa energy',eGW(nO(2)+1,2))
-
-  end if
-
-end subroutine 
diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90
index 7f8c8f9..3643c89 100644
--- a/src/GW/UG0W0.f90
+++ b/src/GW/UG0W0.f90
@@ -86,11 +86,13 @@ subroutine UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTD
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Memory allocation
diff --git a/src/GW/UGW.f90 b/src/GW/UGW.f90
index 8cdb3e9..6a95e25 100644
--- a/src/GW/UGW.f90
+++ b/src/GW/UGW.f90
@@ -1,6 +1,6 @@
-subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT,       &
+subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_diis,doACFDT,                 &
               exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
-              linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,            & 
+              linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,                 & 
               ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
 
 ! GW module
@@ -17,7 +17,6 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
   logical,intent(in)            :: doqsGW
   logical,intent(in)            :: doufG0W0
   logical,intent(in)            :: doufGW
-  logical,intent(in)            :: doSRGqsGW
 
   integer,intent(in)            :: maxSCF
   integer,intent(in)            :: max_diis
@@ -36,7 +35,7 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
   logical,intent(in)            :: spin_flip
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -79,7 +78,7 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
     
     call wall_time(start_GW)
     call UG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & 
-               linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
+               linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,            & 
                dipole_int_aa,dipole_int_bb,cHF,eHF)
     call wall_time(end_GW)
   
@@ -97,7 +96,7 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
     call wall_time(start_GW)
     call evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & 
-              spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,  & 
+              spin_conserved,spin_flip,linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,  & 
               ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF)
     call wall_time(end_GW)
 
@@ -115,7 +114,7 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
     call wall_time(start_GW)
     call qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & 
-               spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,    & 
+               spin_conserved,spin_flip,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,         & 
                S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
     call wall_time(end_GW)
 
@@ -125,24 +124,6 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
 
   end if
 
-!------------------------------------------------------------------------
-! Perform SRG-qsGW calculation
-!------------------------------------------------------------------------
-
-  if(doSRGqsGW) then 
-
-    call wall_time(start_GW)
-    call SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & 
-                   spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,    & 
-                   S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF)
-    call wall_time(end_GW)
-
-    t_GW = end_GW - start_GW
-    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
-    write(*,*)
-
-  end if
-
 !------------------------------------------------------------------------
 ! Perform ufG0W0 calculatiom
 !------------------------------------------------------------------------
diff --git a/src/GW/UGW_phBSE.f90 b/src/GW/UGW_phBSE.f90
index 6d14f64..302a19b 100644
--- a/src/GW/UGW_phBSE.f90
+++ b/src/GW/UGW_phBSE.f90
@@ -57,7 +57,7 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO
 
   double precision,intent(out)  :: EcBSE(nspin)
 
-  ! Memory allocation
+! Memory allocation
 
   nS_aa = nS(1)
   nS_bb = nS(2)
@@ -69,6 +69,15 @@ subroutine UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO
   
   allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin))
 
+!-----!
+! TDA !
+!-----!
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
+    write(*,*)
+  end if
+
 !--------------------------!
 ! Spin-conserved screening !
 !--------------------------!
diff --git a/src/GW/evGGW.f90 b/src/GW/evGGW.f90
index 7c862e9..cc9717e 100644
--- a/src/GW/evGGW.f90
+++ b/src/GW/evGGW.f90
@@ -1,5 +1,5 @@
 subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-                 linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
+                 linearize,eta,doSRG,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF)
 
 ! Perform self-consistent eigenvalue-only GW calculation
 
@@ -27,7 +27,7 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   logical,intent(in)            :: doppBSE
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -79,11 +79,13 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Linear mixing
@@ -127,9 +129,11 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     ! Compute correlation part of the self-energy 
 
-    if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
-
-    call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    if(doSRG) then
+      call GGW_SRG_self_energy_diag(nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    else
+      call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    end if
 
     ! Solve the quasi-particle equation
 
diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90
index 5229e96..f7e9258 100644
--- a/src/GW/evRGW.f90
+++ b/src/GW/evRGW.f90
@@ -83,6 +83,15 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     write(*,*)
   end if
 
+! SRG regularization
+
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
+    write(*,*)
+
+  end if
+
 ! Linear mixing
 
   linear_mixing = .false.
diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90
index f073194..b8a9085 100644
--- a/src/GW/evUGW.f90
+++ b/src/GW/evUGW.f90
@@ -90,11 +90,13 @@ subroutine evUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Initialization
diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90
index ce6a64a..0f32c34 100644
--- a/src/GW/qsGGW.f90
+++ b/src/GW/qsGGW.f90
@@ -1,5 +1,5 @@
 subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & 
-                 eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,Ov,Or,T,V,Hc,ERI_AO,  & 
+                 eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nBas2,nC,nO,nV,nR,nS,EGHF,Ov,Or,T,V,Hc,ERI_AO,  & 
                  ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Generalized version of quasiparticle self-consistent GW 
@@ -25,7 +25,7 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   logical,intent(in)            :: dTDA
   logical,intent(in)            :: doppBSE
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -130,11 +130,13 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Memory allocation
@@ -261,9 +263,11 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     call GGW_excitation_density(nBas2,nC,nO,nR,nS,ERI_MO,XpY,rho)
 
-    if(regularize) call GW_regularization(nBas2,nC,nO,nV,nR,nS,eGW,Om,rho)
-
-    call GGW_self_energy(eta,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    if(doSRG) then
+      call GGW_SRG_self_energy(nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    else
+      call GGW_self_energy(eta,nBas2,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+    end if
 
 !   Make correlation self-energy Hermitian and transform it back to AO basis
    
diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90
index 0894911..f9ebb50 100644
--- a/src/GW/qsUGW.f90
+++ b/src/GW/qsUGW.f90
@@ -126,11 +126,13 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
     write(*,*)
   end if
 
-! TDA 
+! SRG regularization
 
-  if(TDA) then 
-    write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGW scheme ***'
     write(*,*)
+
   end if
 
 ! Memory allocation
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index df2b8bc..8745ef3 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -12,7 +12,7 @@ program QuAcK
   logical                       :: doCIS,doCIS_D,doCID,doCISD,doFCI
   logical                       :: dophRPA,dophRPAx,docrRPA,doppRPA
   logical                       :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3
-  logical                       :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
+  logical                       :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical                       :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh
 
   integer                       :: nNuc
@@ -103,7 +103,7 @@ program QuAcK
                     dophRPA,dophRPAx,docrRPA,doppRPA,      &
                     doG0F2,doevGF2,doqsGF2,doufG0F02,      & 
                     doG0F3,doevGF3,                        &
-                    doG0W0,doevGW,doqsGW,doSRGqsGW,        &
+                    doG0W0,doevGW,doqsGW,                  &
                     doufG0W0,doufGW,                       &
                     doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
                     doG0T0eh,doevGTeh,doqsGTeh,            &
@@ -210,7 +210,7 @@ program QuAcK
   if(doRQuAcK) &
     call RQuAcK(doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,              &
                 dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
-                doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,  &
+                doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,            &
                 doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,                                & 
                 nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                                       &
                 S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
@@ -226,7 +226,7 @@ program QuAcK
   if(doUQuAcK) &
     call UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                     &
                 dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
-                doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,  &
+                doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,            &
                 doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,                                & 
                 nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                                            &
                 S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90
index 2573137..57f0fe5 100644
--- a/src/QuAcK/RQuAcK.f90
+++ b/src/QuAcK/RQuAcK.f90
@@ -1,6 +1,6 @@
 subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,               &   
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & 
-                  doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,  &
+                  doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,            &
                   doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,                                & 
                   nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc,                                                       &
                   S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
@@ -25,7 +25,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
   logical,intent(in)            :: doCIS,doCIS_D,doCID,doCISD,doFCI
   logical,intent(in)            :: dophRPA,dophRPAx,docrRPA,doppRPA
   logical,intent(in)            :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3
-  logical,intent(in)            :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
+  logical,intent(in)            :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical,intent(in)            :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
   logical,intent(in)            :: doG0T0eh,doevGTeh,doqsGTeh
 
@@ -309,12 +309,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
 ! GW module !
 !-----------!
 
-  doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. doSRGqsGW
+  doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW 
 
   if(doGW) then
     
     call wall_time(start_GW)
-    call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,      &
+    call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW,                &
              doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
              lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,               &
              V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90
index b5dc405..a4fb603 100644
--- a/src/QuAcK/read_methods.f90
+++ b/src/QuAcK/read_methods.f90
@@ -6,7 +6,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
                         dophRPA,dophRPAx,docrRPA,doppRPA,      & 
                         doG0F2,doevGF2,doqsGF2,doufG0F02,      &
                         doG0F3,doevGF3,                        & 
-                        doG0W0,doevGW,doqsGW,doSRGqsGW,        & 
+                        doG0W0,doevGW,doqsGW,                  & 
                         doufG0W0,doufGW,                       & 
                         doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, &
                         doG0T0eh,doevGTeh,doqsGTeh,            &
@@ -25,7 +25,7 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
   logical,intent(out)           :: doCIS,doCIS_D,doCID,doCISD,doFCI
   logical,intent(out)           :: dophRPA,dophRPAx,docrRPA,doppRPA
   logical,intent(out)           :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3  
-  logical,intent(out)           :: doG0W0,doevGW,doqsGW,doSRGqsGW,doufG0W0,doufGW
+  logical,intent(out)           :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical,intent(out)           :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp 
   logical,intent(out)           :: doG0T0eh,doevGTeh,doqsGTeh
 
@@ -147,18 +147,16 @@ subroutine read_methods(doRHF,doUHF,doGHF,doROHF,              &
   doG0W0    = .false.
   doevGW    = .false.
   doqsGW    = .false.
-  doSRGqsGW = .false.
   doufG0W0  = .false.
   doufGW    = .false.
 
   read(1,*) 
-  read(1,*) ans1,ans2,ans3,ans4,ans5,ans6
+  read(1,*) ans1,ans2,ans3,ans4,ans5
   if(ans1 == 'T') doG0W0    = .true.
   if(ans2 == 'T') doevGW    = .true.
   if(ans3 == 'T') doqsGW    = .true.
-  if(ans4 == 'T') doSRGqsGW = .true.
-  if(ans5 == 'T') doufG0W0  = .true.
-  if(ans6 == 'T') doufGW    = .true.
+  if(ans4 == 'T') doufG0W0  = .true.
+  if(ans5 == 'T') doufGW    = .true.
 
 ! Read GTpp methods