diff --git a/TODO b/TODO
index 11be081c..7ef32c80 100644
--- a/TODO
+++ b/TODO
@@ -57,4 +57,14 @@ Doc: plugins et qp_plugins
 
 Ajouter les symetries dans devel
 
-Compiler ezfio avec openmp
+
+
+
+IMPORTANT:
+
+Davidson Diagonalization
+------------------------
+
+ Not enough memory: aborting in davidson_diag_hjj_sjj
+
+
diff --git a/ocaml/.gitignore b/ocaml/.gitignore
index 50a9344d..74e83c02 100644
--- a/ocaml/.gitignore
+++ b/ocaml/.gitignore
@@ -9,6 +9,7 @@ Input_ao_two_e_erf_ints.ml
 Input_ao_two_e_ints.ml
 Input_auto_generated.ml
 Input_becke_numerical_grid.ml
+Input_champ.ml
 Input_davidson.ml
 Input_density_for_dft.ml
 Input_determinants.ml
@@ -21,6 +22,7 @@ Input_nuclei.ml
 Input_perturbation.ml
 Input_pseudo.ml
 Input_scf_utils.ml
+Input_variance.ml
 qp_create_ezfio
 qp_create_ezfio.native
 qp_edit
diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index afa7a636..fe99ab5c 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -609,6 +609,12 @@ BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)]
  BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)]
 &BEGIN_PROVIDER[ integer, pt2_R, (N_det_generators)]
   implicit none
+  BEGIN_DOC
+! pt2_J contains the list of generators after ordering them according to the
+! Monte Carlo sampling.
+!
+! pt2_R(i) is the number of combs drawn when determinant i is computed.
+  END_DOC
   integer                :: N_c, N_j
   integer                :: U, t, i
   double precision       :: v
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 69540afc..f1c0ba88 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -286,7 +286,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
     end if
   end do
   deallocate(indices)
-  
+
+!  !$OMP CRITICAL
+!  print *,  'Step1: ', i_generator, preinteresting(0)
+!  !$OMP END CRITICAL
+    
   allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2))
   allocate (mat(N_states, mo_num, mo_num))
   maskInd = -1
@@ -485,6 +489,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
             
             call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
             if(fullMatch) cycle
+! !$OMP CRITICAL
+!  print *,  'Step3: ', i_generator, h1, interesting(0)
+! !$OMP END CRITICAL
             
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
             
@@ -623,6 +630,11 @@ end
 subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
   use bitmasks
   implicit none
+  BEGIN_DOC
+! Computes the contributions A(r,s) by
+! comparing the external determinant to all the internal determinants det(i).
+! an applying two particles (r,s) to the mask.
+  END_DOC
 
   integer, intent(in)            :: sp, i_gen, N_sel
   integer, intent(in)            :: interesting(0:N_sel)
@@ -1198,6 +1210,10 @@ end
 subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
   use bitmasks
   implicit none
+  BEGIN_DOC
+! Identify the determinants in det which are in the internal space. These are
+! the determinants that can be produced by creating two particles on the mask.
+  END_DOC
 
   integer, intent(in) :: i_gen, N
   integer, intent(in) :: interesting(0:N)
@@ -1217,16 +1233,19 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
   end do
 
   genl : do i=1, N
+    ! If det(i) can't be generated by the mask, cycle
     do j=1, N_int
       if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl
       if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
     end do
 
+    ! If det(i) < det(i_gen), it hs already been considered
     if(interesting(i) < i_gen) then
       fullMatch = .true.
       return
     end if
 
+    ! Identify the particles
     do j=1, N_int
       myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1))
       myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f
index b0a419de..ba77f2dd 100644
--- a/src/cipsi/stochastic_cipsi.irp.f
+++ b/src/cipsi/stochastic_cipsi.irp.f
@@ -29,7 +29,7 @@ subroutine run_stochastic_cipsi
   pt2 = -huge(1.e0)
   rpt2 = -huge(1.e0)
   norm = 0.d0
-  variance = 0.d0
+  variance = huge(1.e0)
 
   if (s2_eig) then
     call make_s2_eigenfunction
diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f
index 538085f0..1b9d6011 100644
--- a/src/cipsi/zmq_selection.irp.f
+++ b/src/cipsi/zmq_selection.irp.f
@@ -90,7 +90,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
   mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 +  3.d0 + 5.d0) / (1024.d0**3)
   call write_double(6,mem,'Estimated memory/thread (Gb)')
   if (qp_max_mem > 0) then
-    nproc_target = max(1,int(dble(qp_max_mem)/mem))
+    nproc_target = max(1,int(dble(qp_max_mem)/(0.1d0 + mem)))
     nproc_target = min(nproc_target,nproc)
   endif
 
diff --git a/src/iterations/iterations.irp.f b/src/iterations/iterations.irp.f
index f4c90769..2c9faaf8 100644
--- a/src/iterations/iterations.irp.f
+++ b/src/iterations/iterations.irp.f
@@ -20,9 +20,15 @@ subroutine save_iterations(e_, pt2_,n_)
   END_DOC
   integer, intent(in) :: n_
   double precision, intent(in) :: e_(N_states), pt2_(N_states)
+  integer :: i
 
-  if (N_iter > 100) then
-    return
+  if (N_iter == 101) then
+    do i=2,N_iter-1
+      energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
+      pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) 
+    enddo
+    N_iter = N_iter-1
+    TOUCH N_iter
   endif
 
   energy_iterations(1:N_states,N_iter) = e_(1:N_states)