diff --git a/input/methods b/input/methods
index 74f3c34..4fb62d7 100644
--- a/input/methods
+++ b/input/methods
@@ -13,9 +13,9 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F      F      F    F
 # G0W0* evGW* qsGW* ufG0W0 ufGW
-  F     F     F     F      F
+  T     F     F     F      F
 # G0T0 evGT qsGT 
-   F    F    T
+   F    F    F
 # MCMP2
   F
 # * unrestricted version available
diff --git a/input/options b/input/options
index 3a2b74a..1256f09 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess level_shift stability
-      256   0.0000001  T    5     2          1          F         0.0         F
+      256   0.0000001  T    5     2          1          F         1.0         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -9,12 +9,12 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm reg
       256    0.00001 T    5      T   0.0 3      F
 # GW: maxSCF thresh  DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg
-      256    0.00001 T    5      T   0.0 F      F     T     F   F   F
+      256    0.00001 T    5      T   0.0 F      F     F     F   F   F
 # GT: maxSCF thresh  DIIS n_diis lin eta TDA_T reg
       10    0.00001 T    5      T   0.0 F     F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        F    T    T    F
+        T    T    T    T
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90
index bdebd2a..26ba94d 100644
--- a/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90
+++ b/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90
@@ -29,7 +29,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
   integer                       :: ia
 
-  integer,parameter             :: maxS = 10
+  integer                       :: maxS = 10
   double precision              :: gapGW
 
   integer                       :: nSCF
@@ -43,16 +43,23 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
   double precision,allocatable  :: X(:)
   double precision,allocatable  :: Y(:)
 
-  double precision,allocatable  :: Ap_dyn(:,:)
-  double precision,allocatable  :: Am_dyn(:,:)
-  double precision,allocatable  :: Bp_dyn(:,:)
-  double precision,allocatable  :: Bm_dyn(:,:)
+  double precision,allocatable  ::  Ap_dyn(:,:)
+  double precision,allocatable  :: ZAp_dyn(:,:)
+
+  double precision,allocatable  ::  Bp_dyn(:,:)
+
+  double precision,allocatable  ::  Am_dyn(:,:)
+  double precision,allocatable  :: ZAm_dyn(:,:)
+
+  double precision,allocatable  ::  Bm_dyn(:,:)
 
 ! Memory allocation
 
-  allocate(OmDyn(nS),OmOld(nS),X(nS),Y(nS),Ap_dyn(nS,nS))
+  maxS = min(nS,maxS)
+  allocate(OmDyn(maxS),OmOld(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS))
 
-  if(.not.dTDA) allocate(Am_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS))
+  if(.not.dTDA) &
+    allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS))
 
   if(dTDA) then
     write(*,*)
@@ -64,7 +71,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
   Conv = 1d0
   nSCF = 0
-  OmOld(:) = OmBSE(:)
+  OmOld(1:maxS) = OmBSE(1:maxS)
 
   write(*,*) '---------------------------------------------------------------------------------------------------'
   write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies                     '
@@ -83,9 +90,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
     write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Convergence (eV)'
     write(*,*) '---------------------------------------------------------------------------------------------------'
 
-    do ia=1,min(nS,maxS)
+    do ia=1,maxS
 
- 
       X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
       Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
  
@@ -95,7 +101,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
 
        ! Resonant part of the BSE correction
  
-        call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn)
+        call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,ZAp_dyn)
  
         OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:)))
  
@@ -118,8 +124,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,
  
     end do
 
-    Conv = maxval(abs(OmBSE(:) + OmDyn(:) - OmOld(:)))*HaToeV
-    OmOld(:) = OmBSE(:) + OmDyn(:)
+    Conv = maxval(abs(OmBSE(1:maxS) + OmDyn(:) - OmOld(:)))*HaToeV
+    OmOld(:) = OmBSE(1:maxS) + OmDyn(:)
 
     write(*,*) '---------------------------------------------------------------------------------------------------'
     write(*,'(2X,A20,1X,F10.6)') ' Convergence = ',Conv