diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f
index 1279b5cd..cf770049 100644
--- a/src/cipsi/cipsi.irp.f
+++ b/src/cipsi/cipsi.irp.f
@@ -128,7 +128,13 @@ subroutine run_cipsi
     if (qp_stop()) exit
   enddo
 
-  if ((.not.qp_stop()).and.(N_det > N_det_max)) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
+  if ((.not.qp_stop()).and.                                          &
+        (N_det > N_det_max) .and.                                    &
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
+        ) then
     if (do_pt2) then
       call pt2_dealloc(pt2_data)
       call pt2_dealloc(pt2_data_err)
diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f
index deffaeb7..339f7084 100644
--- a/src/cipsi/stochastic_cipsi.irp.f
+++ b/src/cipsi/stochastic_cipsi.irp.f
@@ -118,7 +118,13 @@ subroutine run_stochastic_cipsi
     if (qp_stop()) exit
   enddo
 
-  if ((.not.qp_stop()).and.(N_det > N_det_max)) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
+  if ((.not.qp_stop()).and.                                          &
+        (N_det > N_det_max) .and.                                    &
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
+        ) then
     call pt2_dealloc(pt2_data)
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f
index dc42b9a8..346f1cf9 100644
--- a/src/davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/davidson/diagonalization_hs2_dressed.irp.f
@@ -465,9 +465,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        integer :: lwork, info
        double precision, allocatable :: work(:)
 
-
        y = h
-!       y = h_p
+!       y = h_p   ! Doesn't work for non-singlets
        lwork = -1
        allocate(work(1))
        call dsygv(1,'V','U',shift2,y,size(y,1), &