diff --git a/plugins/DavidsonDressed/diagonalize_CI.irp.f b/plugins/DavidsonDressed/diagonalize_CI.irp.f
index e01dddb3..fb80bb9a 100644
--- a/plugins/DavidsonDressed/diagonalize_CI.irp.f
+++ b/plugins/DavidsonDressed/diagonalize_CI.irp.f
@@ -176,29 +176,27 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
 
   h_matrix_dressed(1:N_det,1:N_det) = h_matrix_all_dets(1:N_det,1:N_det)
 
-!  do k=1,N_states
-!    do j=1,N_det
-!      do i=1,N_det
-!        h_matrix_dressed(i,j) = h_matrix_dressed(i,j) +       &
-!            0.5d0 * (dressing_column_h(i,k) * psi_coef(j,k) + &
-!                     dressing_column_h(j,k) * psi_coef(i,k))
-!      enddo
-!    enddo
-!  enddo
-
-  integer :: l,jj
-  double precision :: f
-  do k=1,N_states
-    l = dressed_column_idx(k)
-    f = 1.d0/psi_coef(l,k)
+  if (N_states == 1) then
+    integer                        :: l,jj
+    double precision               :: f
+    l = dressed_column_idx(1)
+    f = 0.5d0/psi_coef(l,1)
     do i=1,N_det
-       if (i==l) cycle
-       h_matrix_dressed(i,l) = h_matrix_dressed(i,l) + dressing_column_h(i,k) *f
-       h_matrix_dressed(l,i) = h_matrix_dressed(l,i) + dressing_column_h(i,k) *f
-       h_matrix_dressed(l,l) = h_matrix_dressed(l,l) - psi_coef(i,k) * dressing_column_h(i,k) *f*f
+      h_matrix_dressed(i,l) = h_matrix_dressed(i,l) + dressing_column_h(i,1) *f
+      h_matrix_dressed(l,i) = h_matrix_dressed(l,i) + dressing_column_h(i,1) *f
     enddo
+  else
+    do k=1,N_states
+      do j=1,N_det
+        do i=1,N_det
+          h_matrix_dressed(i,j) = h_matrix_dressed(i,j) +            &
+              0.5d0 *(dressing_column_h(i,k) * psi_coef(j,k) +       &
+              dressing_column_h(j,k) * psi_coef(i,k))
+        enddo
+      enddo
+    enddo
+  endif
 
-  enddo
 
 END_PROVIDER
 
diff --git a/plugins/mrcepa0/dressing_vector.irp.f b/plugins/mrcepa0/dressing_vector.irp.f
index 2f2a535c..b233d941 100644
--- a/plugins/mrcepa0/dressing_vector.irp.f
+++ b/plugins/mrcepa0/dressing_vector.irp.f
@@ -16,12 +16,15 @@
 
  do k=1,N_states
   integer :: jj
+  l = dressed_column_idx(k)
    do jj=1,N_det_non_ref
      j = idx_non_ref(jj)
-     dressing_column_h(j,k) = delta_ij   (k,jj) 
-     dressing_column_s(j,k) = delta_ij_s2(k,jj)
+     dressing_column_h(j,k) = 2.d0 * delta_ij   (k,jj) 
+     dressing_column_s(j,k) = 2.d0 * delta_ij_s2(k,jj)
+     dressing_column_h(l,k) -= psi_coef(j,k) * delta_ij(k,jj) /psi_coef(l,k)
   enddo
  enddo
+
 END_PROVIDER
 
 
diff --git a/src/Davidson/diagonalization_hs2_dressed.irp.f b/src/Davidson/diagonalization_hs2_dressed.irp.f
index f629c0c9..7e4dc047 100644
--- a/src/Davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/Davidson/diagonalization_hs2_dressed.irp.f
@@ -281,35 +281,53 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       
       if (dressing_state > 0) then
 
-        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
-          psi_coef, size(psi_coef,1), &
-          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+        if (N_st == 1) then
 
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
-          1.d0, W(1,shift+1), size(W,1))
+          l = dressed_column_idx(1)
+          double precision :: f
+          f = 0.5d0/psi_coef(l,1)
+          do istate=1,N_st_diag
+            do i=1,sze
+              W(i,shift+istate) += dressing_column_h(i,1) *f * U(l,shift+istate)
+              W(l,shift+istate) += dressing_column_h(i,1) *f * U(i,shift+istate)
+              S(i,shift+istate) += dressing_column_s(i,1) *f * U(l,shift+istate)
+              S(l,shift+istate) += dressing_column_s(i,1) *f * U(i,shift+istate)
+            enddo
 
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
-          1.d0, S(1,shift+1), size(S,1))
+          enddo
+
+        else 
+
+          call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+            psi_coef, size(psi_coef,1), &
+            U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+
+          call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+            dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
+            1.d0, W(1,shift+1), size(W,1))
+
+          call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+            dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
+            1.d0, S(1,shift+1), size(S,1))
 
 
-        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
-          dressing_column_h, size(dressing_column_h,1), &
-          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+          call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+            dressing_column_h, size(dressing_column_h,1), &
+            U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
 
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
-          1.d0, W(1,shift+1), size(W,1))
+          call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+            psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+            1.d0, W(1,shift+1), size(W,1))
 
-        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
-          dressing_column_s, size(dressing_column_s,1), &
-          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+          call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+            dressing_column_s, size(dressing_column_s,1), &
+            U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
 
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
-          1.d0, S(1,shift+1), size(S,1))
+          call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+            psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+            1.d0, S(1,shift+1), size(S,1))
 
+        endif
       endif
 
       ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>