diff --git a/configure b/configure
index 893c7148..c328c4f0 100755
--- a/configure
+++ b/configure
@@ -231,7 +231,7 @@ EOF
 EOF
     elif [[ ${PACKAGE} = qmckl ]] ; then
 
-            VERSION=0.5.3
+            VERSION=0.5.4
             execute << EOF
               cd "\${QP_ROOT}"/external
               wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
diff --git a/external/ezfio b/external/ezfio
index dba01c4f..d5805497 160000
--- a/external/ezfio
+++ b/external/ezfio
@@ -1 +1 @@
-Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3
+Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f
index 7b559925..1ead9d78 100644
--- a/src/davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/davidson/diagonalization_hs2_dressed.irp.f
@@ -286,7 +286,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
       ! Small
       h(N_st_diag*itermax,N_st_diag*itermax),                        &
-      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
+!      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
       y(N_st_diag*itermax,N_st_diag*itermax),                        &
       s_(N_st_diag*itermax,N_st_diag*itermax),                       &
       s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
@@ -340,7 +340,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       exit
     endif
 
-    do iter=1,itermax-1
+    iter = 0
+    do while (iter < itermax-1)
+      iter += 1
+!    do iter=1,itermax-1
 
       shift  = N_st_diag*(iter-1)
       shift2 = N_st_diag*iter
@@ -430,30 +433,30 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
       call dgemm('T','N', shift2, shift2, sze,                       &
           1.d0, U, size(U,1), W, size(W,1),                          &
-          0.d0, h, size(h_p,1))
+          0.d0, h, size(h,1))
       call dgemm('T','N', shift2, shift2, sze,                       &
           1.d0, U, size(U,1), U, size(U,1),                          &
           0.d0, s_tmp, size(s_tmp,1))
 
-      ! Penalty method
-      ! --------------
-
-      if (s2_eig) then
-        h_p = s_
-        do k=1,shift2
-          h_p(k,k) = h_p(k,k) - expected_s2
-        enddo
-        if (only_expected_s2) then
-          alpha = 0.1d0
-          h_p = h + alpha*h_p
-        else
-          alpha = 0.0001d0
-          h_p = h + alpha*h_p
-        endif
-      else
-        h_p = h
-        alpha = 0.d0
-      endif
+!      ! Penalty method
+!      ! --------------
+!
+!      if (s2_eig) then
+!        h_p = s_
+!        do k=1,shift2
+!          h_p(k,k) = h_p(k,k) - expected_s2
+!        enddo
+!        if (only_expected_s2) then
+!          alpha = 0.1d0
+!          h_p = h + alpha*h_p
+!        else
+!          alpha = 0.0001d0
+!          h_p = h + alpha*h_p
+!        endif
+!      else
+!        h_p = h
+!        alpha = 0.d0
+!      endif
 
       ! Diagonalize h_p
       ! ---------------
@@ -473,8 +476,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        call dsygv(1,'V','U',shift2,y,size(y,1), &
           s_tmp,size(s_tmp,1), lambda, work,lwork,info)
        deallocate(work)
-       if (info /= 0) then
-         stop 'DSYGV Diagonalization failed'
+       if (info > 0) then
+         ! Numerical errors propagate. We need to reduce the number of iterations
+         itermax = iter-1
+         exit
        endif
 
       ! Compute Energy for each eigenvector
diff --git a/src/iterations/print_extrapolation.irp.f b/src/iterations/print_extrapolation.irp.f
index a7f85693..24c9845f 100644
--- a/src/iterations/print_extrapolation.irp.f
+++ b/src/iterations/print_extrapolation.irp.f
@@ -37,7 +37,7 @@ subroutine print_extrapolated_energy
     write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
     do k=2,N_iter_p
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
+      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter_p+1-k), extrapolated_energy(k,i), &
           extrapolated_energy(k,i) - extrapolated_energy(k,1), &
           (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
     enddo