diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f
index b8a551b5..3a4922fb 100644
--- a/src/Dets/H_apply.irp.f
+++ b/src/Dets/H_apply.irp.f
@@ -113,6 +113,10 @@ end
 subroutine copy_H_apply_buffer_to_wf
   use omp_lib
   implicit none
+  BEGIN_DOC
+! Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
+! after calling this function.
+  END_DOC
   integer(bit_kind), allocatable :: buffer_det(:,:,:)
   double precision, allocatable  :: buffer_coef(:,:)
   integer                        :: i,j,k
@@ -144,7 +148,6 @@ subroutine copy_H_apply_buffer_to_wf
   do j=0,nproc-1
     N_det = N_det + H_apply_buffer(j)%N_det
   enddo
-  TOUCH N_det
   
   if (psi_det_size < N_det) then
     psi_det_size = N_det
@@ -188,8 +191,7 @@ subroutine copy_H_apply_buffer_to_wf
   H_apply_buffer(j)%N_det = 0
   !$OMP END PARALLEL
   call normalize(psi_coef,N_det)
-  
-  SOFT_TOUCH psi_det psi_coef
+  SOFT_TOUCH psi_det psi_coef N_det
   
 end
 
diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f
index ccab5080..855e526e 100644
--- a/src/Dets/H_apply_template.f
+++ b/src/Dets/H_apply_template.f
@@ -334,7 +334,6 @@ subroutine $subroutine($params_main)
   
   $copy_buffer
   $generate_psi_guess
-  SOFT_TOUCH psi_det psi_coef
   
 end
 
diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f
index cd7e5444..1f1ad1cf 100644
--- a/src/Dets/diagonalize_CI.irp.f
+++ b/src/Dets/diagonalize_CI.irp.f
@@ -66,3 +66,17 @@ END_PROVIDER
   
 END_PROVIDER
 
+subroutine diagonalize_CI
+  implicit none
+  BEGIN_DOC
+!  Replace the coefficients of the CI states by the coefficients of the 
+!  eigenstates of the CI matrix
+  END_DOC
+  integer :: i,j
+  do j=1,N_states
+    do i=1,N_det
+      psi_coef(i,j) = CI_eigenvectors(i,j)
+    enddo
+  enddo
+  SOFT_TOUCH psi_coef psi_det CI_electronic_energy CI_energy CI_eigenvectors
+end
diff --git a/src/Perturbation/cisd_selection.irp.f b/src/Perturbation/cisd_selection.irp.f
index 77adaa0a..fa1891e5 100644
--- a/src/Perturbation/cisd_selection.irp.f
+++ b/src/Perturbation/cisd_selection.irp.f
@@ -1,7 +1,6 @@
 program cisd
   implicit none
   integer                        :: i,k
-  double precision, allocatable  :: eigvalues(:),eigvectors(:,:)
 
   
   double precision, allocatable  :: pt2(:), norm_pert(:)
@@ -9,16 +8,17 @@ program cisd
   integer                        :: N_st, iter
   N_st = N_states
   allocate (pt2(N_st), norm_pert(N_st))
-  allocate(eigvalues(n_states),eigvectors(n_det,n_states))
   
   pt2 = 1.d0
   do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
     E_old = CI_energy(1)
     call H_apply_cisd_selection(pt2, norm_pert, H_pert_diag,  N_st)
+    call diagonalize_CI
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
     print *,  'PT2      = ', pt2
     print *,  'E        = ', E_old
     print *,  'E+PT2    = ', E_old+pt2
   enddo
+  deallocate(pt2,norm_pert)
 end