diff --git a/INSTALL.rst b/INSTALL.rst
index 7bcc114c..acf918f8 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -43,6 +43,7 @@ Requirements
 - `Zlib`_
 - `GNU Patch`_
 - |ZeroMQ| : networking library
+- `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library
 - |OCaml| compiler with |OPAM| package manager 
 - |Ninja| : a parallel build system
 
@@ -84,6 +85,7 @@ The following packages are supported by the :command:`configure` installer:
 * irpf90 
 * zeromq 
 * f77zmq
+* gmp
 * ocaml  ( :math:`\approx` 10 minutes)
 * ezfio 
 * docopt 
@@ -110,11 +112,11 @@ If the :command:`configure` executable fails to install a specific dependency
 -----------------------------------------------------------------------------
 
 If the :command:`configure` executable does not succeed to install a specific dependency, 
-there are some proposition of how to download and install the minimal dependencies to compile and use the |QP|
+there are some proposition of how to download and install the minimal dependencies to compile and use the |QP|.
 
 
 Before doing anything below, try to install the packages with your package manager
-(:command:`apt`, :command:`yum`, etc)
+(:command:`apt`, :command:`yum`, etc).
 
 
 Ninja
@@ -216,6 +218,29 @@ With Debian or Ubuntu, you can use
 
    sudo apt install zlib1g-dev
 
+GMP
+---
+
+GMP is the GNU Multiple Precision Arithmetic Library.
+
+* Download the latest version of gmp here:
+  `<ftp://ftp.gnu.org/gnu/gmp/gmp-6.1.2.tar.bz2>`_
+  and move it in the :file:`${QP_ROOT}/external` directory
+
+* Extract the archive, go into the :file:`gmp-*` directory and run
+  the following commands
+
+.. code:: bash
+
+   ./configure --prefix=${QP_ROOT}
+   make
+   make install
+
+With Debian or Ubuntu, you can use
+
+.. code:: bash
+
+   sudo apt install libgmp-dev
 
 
 OCaml
diff --git a/TODO b/TODO
index dbe90ecf..7774e5f7 100644
--- a/TODO
+++ b/TODO
@@ -57,6 +57,13 @@ Doc: plugins et qp_plugins
 
 Ajouter les symetries dans devel
 
-Compiler ezfio avec openmp
 
 # Parallelize i_H_psi
+
+
+IMPORTANT:
+
+Davidson Diagonalization
+------------------------
+
+ Not enough memory: aborting in davidson_diag_hjj_sjj
diff --git a/config/cray.cfg b/config/cray.cfg
new file mode 100644
index 00000000..9a4c19cf
--- /dev/null
+++ b/config/cray.cfg
@@ -0,0 +1,67 @@
+# Common flags
+##############
+#
+# -mkl=[parallel|sequential] : Use the MKL library
+# --ninja                 : Allow the utilisation of ninja. It is mandatory !
+# --align=32                 : Align all provided arrays on a 32-byte boundary
+#
+[COMMON]
+FC           : ftn -dynamic -fpic
+LAPACK_LIB   : 
+IRPF90       : irpf90
+IRPF90_FLAGS : --ninja --align=32 -DMPI  
+# For KNL, use:
+#IRPF90_FLAGS : --ninja --align=64 -DMPI  
+
+# Global options
+################
+#
+# 1 : Activate
+# 0 : Deactivate
+# 
+[OPTION]
+MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+CACHE   : 0          ; Enable cache_compile.py
+OPENMP  : 1          ; Append OpenMP flags
+
+# Optimization flags
+####################
+#
+# -xHost                     : Compile a binary optimized for the current architecture
+# -O2                        : O3 not better than O2.
+# -ip                        : Inter-procedural optimizations
+# -ftz                       : Flushes denormal results to zero
+#
+[OPT]
+FCFLAGS  : -xCORE-AVX2 -O2 -ip -ftz -g -traceback -qopt-prefetch=5 -qopt-prefetch-issue-excl-hint -unroll-aggressive 
+#
+#For KNL, use:
+#FCFLAGS : -xMIC-AVX512 -O2 -ip -ftz -g -traceback -qopt-prefetch=5 -qopt-prefetch-issue-excl-hint -unroll-aggressive 
+#
+# Profiling flags
+#################
+#
+[PROFILE]
+FC       : -p -g
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz -qopt-prefetch
+
+# Debugging flags
+#################
+#
+# -traceback   : Activate backtrace on runtime
+# -fpe0        : All floating point exaceptions
+# -C           : Checks uninitialized variables,  array subscripts, etc...
+# -g           : Extra debugging information
+# -xSSE2       : Valgrind needs a very simple x86 executable
+#
+[DEBUG]
+FC      : -g -traceback
+FCFLAGS : -xCORE-AVX2 -C -fpe0 -traceback
+
+# OpenMP flags
+#################
+#
+[OPENMP]
+FC           : -qopenmp
+IRPF90_FLAGS : --openmp
+
diff --git a/configure b/configure
index 5735563a..767900cd 100755
--- a/configure
+++ b/configure
@@ -183,7 +183,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
 fi
 
 if [[ ${PACKAGES} = all ]] ; then
-        PACKAGES="zlib ninja irpf90 zeromq f77zmq ocaml ezfio docopt resultsFile bats"
+        PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp ocaml ezfio docopt resultsFile bats"
 fi
 
 
@@ -200,6 +200,21 @@ for PACKAGE in ${PACKAGES} ; do
 EOF
 
 
+    elif [[ ${PACKAGE} = gmp ]] ; then
+
+            download \
+              "ftp://ftp.gnu.org/gnu/gmp/gmp-6.1.2.tar.bz2" \
+              "${QP_ROOT}"/external/gmp.tar.bz2
+            execute << EOF
+              cd "\${QP_ROOT}"/external
+              tar --bzip2 --extract --file gmp.tar.bz2
+              rm gmp.tar.bz2
+              cd gmp*
+              ./configure --prefix=$QP_ROOT && make -j 8
+              make install
+EOF
+    
+
     elif [[ ${PACKAGE} = irpf90 ]] ; then
 
             # When changing version of irpf90, don't forget to update etc/irpf90.rc
@@ -388,6 +403,12 @@ if [[ ${F77ZMQ} = $(not_found) ]] ; then
         fail
 fi
 
+GMP=$(find_lib -lgmp)
+if [[ ${ZLIB} = $(not_found) ]] ; then
+        error "GMP (gmp) is not installed."
+        fail
+fi
+
 OPAM=$(find_exe opam)
 if [[ ${OPAM} = $(not_found) ]] ; then
         error "OPAM (ocaml) package manager is not installed."
diff --git a/docs/source/research.bib b/docs/source/research.bib
index 1676cf24..cf56f2d1 100644
--- a/docs/source/research.bib
+++ b/docs/source/research.bib
@@ -35,6 +35,17 @@
 
 
 %%%% PUBLISHED PAPERS
+@phdthesis{yann_garniron_2019_2558127,
+  author       = {Yann Garniron},
+  title        = {{Development and parallel implementation of 
+                   selected configuration interaction methods}},
+  school       = {Université de Toulouse},
+  year         = 2019,
+  month        = feb,
+  doi          = {10.5281/zenodo.2558127},
+  url          = {https://doi.org/10.5281/zenodo.2558127}
+}
+
 @article{Giner_2018,
     doi = {10.1063/1.5052714},
     url = {https://doi.org/10.1063%2F1.5052714},
diff --git a/etc/local.rc b/etc/local.rc
index 772914df..83d4b34b 100644
--- a/etc/local.rc
+++ b/etc/local.rc
@@ -10,10 +10,10 @@
 # export QP_NTHREADS_PT2=32
 
 # Number of threads for MKL
-# exort MKL_NUM_THREADS=16
+# export MKL_NUM_THREADS=16
 
 # Number of threads for OpenMP
-# exort OMP_NUM_THREADS=16
+# export OMP_NUM_THREADS=16
 
 # Name of the network interface to be chosen
 # export QP_NIC=ib0
diff --git a/etc/qp.rc b/etc/qp.rc
index 387025c7..fd9b4031 100644
--- a/etc/qp.rc
+++ b/etc/qp.rc
@@ -163,7 +163,7 @@ _qp_Complete()
       create_ezfio)
         case "${prev}" in
           create_ezfio)
-            COMPREPLY=( $(compgen -W "-b -a -c -d -h -m -o -p -x" -- $cur ) )
+            COMPREPLY=( $(compgen -W "-b -a -c -d -h -m -o -p -x $(\ls)" -- $cur ) )
             return 0
             ;;
           -m|-d|-c)
diff --git a/ocaml/.gitignore b/ocaml/.gitignore
index fdb5adf4..10d69488 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
diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index afa7a636..bf78bb21 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -11,7 +11,7 @@ END_PROVIDER
   implicit none
   logical, external :: testTeethBuilding
   integer :: i,j
-  pt2_n_tasks_max = elec_beta_num*elec_beta_num + elec_alpha_num*elec_beta_num  - n_core_orb*2 
+  pt2_n_tasks_max = elec_alpha_num*elec_alpha_num + elec_alpha_num*elec_beta_num  - n_core_orb*2 
   pt2_n_tasks_max = min(pt2_n_tasks_max,1+N_det_generators/10000)
   call write_int(6,pt2_n_tasks_max,'pt2_n_tasks_max')
 
@@ -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 791656ce..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
@@ -64,6 +64,7 @@ subroutine run_stochastic_cipsi
   do while (                                                         &
         (N_det < N_det_max) .and.                                    &
         (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (maxval(abs(variance(1:N_states))) > variance_max) .and.     &
         (correlation_energy_ratio <= correlation_energy_ratio_max)   &
         )
       write(*,'(A)')  '--------------------------------------------------------------------------------'
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/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f
index fdc7a6e4..4a62dcf3 100644
--- a/src/hartree_fock/fock_matrix_hf.irp.f
+++ b/src/hartree_fock/fock_matrix_hf.irp.f
@@ -4,7 +4,7 @@
  use map_module
  implicit none
  BEGIN_DOC
- ! Alpha Fock matrix in AO basis set
+ ! Alpha and Beta Fock matrices in AO basis set
  END_DOC
 
  integer                        :: i,j,k,l,k1,r,s
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)
diff --git a/src/perturbation/EZFIO.cfg b/src/perturbation/EZFIO.cfg
index b90a8961..168c3b0d 100644
--- a/src/perturbation/EZFIO.cfg
+++ b/src/perturbation/EZFIO.cfg
@@ -10,6 +10,12 @@ doc: The selection process stops when the largest |PT2| (for all the state) is l
 interface: ezfio,provider,ocaml
 default: 0.0001
 
+[variance_max]
+type: Positive_float
+doc: The selection process stops when the largest variance (for all the state) is lower than `variance_max` in absolute value
+interface: ezfio,provider,ocaml
+default: 0.0
+
 [pt2_relative_error]
 type: Normalized_float
 doc: Stop stochastic |PT2| when the relative error is smaller than `pT2_relative_error`
diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f
index ed7a17cd..8186037b 100644
--- a/src/scf_utils/diagonalize_fock.irp.f
+++ b/src/scf_utils/diagonalize_fock.irp.f
@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num) ]
    implicit none
    BEGIN_DOC
-   ! Eigenvector of the Fock matrix in the MO basis obtained with level shift.
+   ! Eigenvectors of the Fock matrix in the MO basis obtained with level shift.
    END_DOC
 
    integer                        :: i,j