diff --git a/TODO b/TODO
index 645de479..4853fdac 100644
--- a/TODO
+++ b/TODO
@@ -55,12 +55,17 @@ Refaire les benchmarks
 # Documentation de qpsh
 
 # Documentation de /etc
+# Extrapolation qui prend aussi en compe la variance
+# Parler dans le papier de rPT2
 
 # Toto
-Selection d'etats avec qp_edit --state
-singles_alpha_csc_idx, singles_alpha_size  | git diff 
-Environment variables dans qp_run (prefix,etc)
-
 Si un provider est un programme, generer une page a lui tout seul avec le man
 
 Options obligatoires dans Command_line.ml
+Documentation des commandes a updater
+
+# Completion foireuse:
+  qp create_ezfio_from_xyz -b 6-31g ./<Tab>
+  qp create_ezfio_from_xyz -b <Tab>
+
+# Verifier la doc des commandes avec les options - et --
diff --git a/bin/qpsh b/bin/qpsh
new file mode 100755
index 00000000..f7a00546
--- /dev/null
+++ b/bin/qpsh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+export QP_ROOT=$(dirname $0)/..
+
+exec bash --init-file <(cat << EOF
+[[ -f \${HOME}/.bashrc ]] && source \${HOME}/.bashrc 
+source \${QP_ROOT}/quantum_package.rc
+qp prompt
+EOF
+)
+
+
+
+
diff --git a/configure b/configure
index 7a9cdad3..47dce84c 100755
--- a/configure
+++ b/configure
@@ -3,7 +3,7 @@
 # Quantum Package configuration script
 #
 
-export QP_ROOT="$( cd "$(dirname "$0")" ; pwd --physical )"
+export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )"
 echo "QP_ROOT="$QP_ROOT
 
 
@@ -61,7 +61,7 @@ done
 # Trim leading and trailing spaces
 PACKAGES=$(echo $PACKAGES | xargs)
 
-echo "export QP_ROOT=\"$QP_ROOT\"" > ${QP_ROOT}/etc/00.qp_root
+echo "export QP_ROOT=\"$QP_ROOT\"" > ${QP_ROOT}/etc/00.qp_root.rc
 
 source quantum_package.rc
 
@@ -149,6 +149,7 @@ for PACKAGE in ${PACKAGES} ; do
 
     elif [[ ${PACKAGE} = irpf90 ]] ; then
 
+            # When changing version of irpf90, don't forget to update etc/irpf90.rc
             download \
               "https://gitlab.com/scemama/irpf90/-/archive/v1.7.4/irpf90-v1.7.4.tar.gz" \
               "${QP_ROOT}"/external/irpf90.tar.gz
@@ -158,15 +159,6 @@ for PACKAGE in ${PACKAGES} ; do
               rm irpf90.tar.gz
               cd irpf90-*
               make
-              for i in irpf90 irpman irpf90_indent
-              do
-                      cat << EOF > "${QP_ROOT}"/bin/$i
-#!/bin/sh
-exec $PWD/bin/$i \$@
-EOF
-
-                      chmod +x "${QP_ROOT}"/bin/$i
-              done
             )
     
 
@@ -234,7 +226,7 @@ EOF
                  --yes --comp=4.07.0
 
             eval $(${QP_ROOT}/bin/opam env)
-            opam install -y ${OCAML_PACKAGES}
+            opam install -y ${OCAML_PACKAGES}  || exit 1
 
             
     elif [[ ${PACKAGE} = ezfio ]] ; then
diff --git a/docs/source/modules/cipsi.rst b/docs/source/modules/cipsi.rst
index bf084f35..6dc9f0af 100644
--- a/docs/source/modules/cipsi.rst
+++ b/docs/source/modules/cipsi.rst
@@ -517,48 +517,6 @@ Subroutines / functions
 
 
 
-.. c:function:: get_m0
-
-    .. code:: text
-
-        subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-
-    File: :file:`selection.irp.f`
-
-    
-
-
-
-
-
-.. c:function:: get_m1
-
-    .. code:: text
-
-        subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-
-    File: :file:`selection.irp.f`
-
-    
-
-
-
-
-
-.. c:function:: get_m2
-
-    .. code:: text
-
-        subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-
-    File: :file:`selection.irp.f`
-
-    
-
-
-
-
-
 .. c:function:: get_mask_phase
 
     .. code:: text
diff --git a/docs/source/modules/davidson.rst b/docs/source/modules/davidson.rst
index f2d9fbe7..0b8c30bd 100644
--- a/docs/source/modules/davidson.rst
+++ b/docs/source/modules/davidson.rst
@@ -449,7 +449,7 @@ Subroutines / functions
 
         subroutine H_S2_u_0_nstates_openmp_work_1(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
 
-    File: :file:`u0_h_u0.irp.f_template_477`
+    File: :file:`u0_h_u0.irp.f_template_576`
 
     Computes  :math:`v_t = H|u_t angle`  and  :math:`s_t = S^2 |u_t angle` 
 
@@ -465,7 +465,7 @@ Subroutines / functions
 
         subroutine H_S2_u_0_nstates_openmp_work_2(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
 
-    File: :file:`u0_h_u0.irp.f_template_477`
+    File: :file:`u0_h_u0.irp.f_template_576`
 
     Computes  :math:`v_t = H|u_t angle`  and  :math:`s_t = S^2 |u_t angle` 
 
@@ -481,7 +481,7 @@ Subroutines / functions
 
         subroutine H_S2_u_0_nstates_openmp_work_3(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
 
-    File: :file:`u0_h_u0.irp.f_template_477`
+    File: :file:`u0_h_u0.irp.f_template_576`
 
     Computes  :math:`v_t = H|u_t angle`  and  :math:`s_t = S^2 |u_t angle` 
 
@@ -497,7 +497,7 @@ Subroutines / functions
 
         subroutine H_S2_u_0_nstates_openmp_work_4(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
 
-    File: :file:`u0_h_u0.irp.f_template_477`
+    File: :file:`u0_h_u0.irp.f_template_576`
 
     Computes  :math:`v_t = H|u_t angle`  and  :math:`s_t = S^2 |u_t angle` 
 
@@ -513,7 +513,7 @@ Subroutines / functions
 
         subroutine H_S2_u_0_nstates_openmp_work_N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
 
-    File: :file:`u0_h_u0.irp.f_template_477`
+    File: :file:`u0_h_u0.irp.f_template_576`
 
     Computes  :math:`v_t = H|u_t angle`  and  :math:`s_t = S^2 |u_t angle` 
 
diff --git a/docs/source/modules/determinants.rst b/docs/source/modules/determinants.rst
index 3878e945..c71fd51a 100644
--- a/docs/source/modules/determinants.rst
+++ b/docs/source/modules/determinants.rst
@@ -199,7 +199,7 @@ Providers
 
     Norm of the  :math:`\alpha`  and  :math:`\beta`  spin determinants in the wave function: 
 
-    ||Da||_i \sum_j C_{ij}**2
+     :math:`||D_\alpha||_i = \sum_j C_{ij}^2`
 
 
 
@@ -215,7 +215,7 @@ Providers
 
     Norm of the  :math:`\alpha`  and  :math:`\beta`  spin determinants in the wave function: 
 
-    ||Da||_i \sum_j C_{ij}**2
+     :math:`||D_\alpha||_i = \sum_j C_{ij}^2`
 
 
 
@@ -451,7 +451,7 @@ Providers
         integer(bit_kind), allocatable	:: psi_det_alpha_unique	(N_int,psi_det_size)
         integer	:: n_det_alpha_unique
 
-    File: :file:`spindeterminants.irp.f_template_143`
+    File: :file:`spindeterminants.irp.f_template_144`
 
     Unique  :math:`\alpha`  determinants
 
@@ -465,7 +465,7 @@ Providers
         integer(bit_kind), allocatable	:: psi_det_beta_unique	(N_int,psi_det_size)
         integer	:: n_det_beta_unique
 
-    File: :file:`spindeterminants.irp.f_template_143`
+    File: :file:`spindeterminants.irp.f_template_144`
 
     Unique  :math:`\beta`  determinants
 
@@ -769,9 +769,9 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Sparse coefficient matrix if the wave function is expressed in a bilinear form : D_a^t C D_b 
+    Sparse coefficient matrix if the wave function is expressed in a bilinear form :  :math:`D_\alpha^\dagger.C.D_\beta` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta.` 
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` . 
 
     Order refers to psi_det
 
@@ -790,9 +790,9 @@ Providers
 
      :math:`D_\alpha^\dagger.C.D_\beta` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta.` 
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` . 
 
-    Order refers to psi_det
+    Order refers to :c:data:`psi_det`
 
 
 
@@ -808,9 +808,9 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Sparse coefficient matrix if the wave function is expressed in a bilinear form : D_a^t C D_b 
+    Sparse coefficient matrix if the wave function is expressed in a bilinear form :  :math:`D_\alpha^\dagger.C.D_\beta` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta.` 
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` . 
 
     Order refers to psi_det
 
@@ -825,7 +825,7 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Order which allows to go from psi_bilinear_matrix to psi_det
+    Order which allows to go from :c:data:`psi_bilinear_matrix` to :c:data:`psi_det`
 
 
 
@@ -838,7 +838,7 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Order which allows to go from psi_bilinear_matrix_order_transp to psi_bilinear_matrix
+    Order which allows to go from :c:data:`psi_bilinear_matrix_order_transp` to :c:data:`psi_bilinear_matrix`
 
 
 
@@ -854,9 +854,9 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Sparse coefficient matrix if the wave function is expressed in a bilinear form : D_a^t C D_b 
+    Sparse coefficient matrix if the wave function is expressed in a bilinear form :  :math:`D_\alpha^\dagger.C.D_\beta` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta.` 
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` . 
 
     Order refers to psi_det
 
@@ -874,11 +874,11 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Transpose of psi_bilinear_matrix 
+    Transpose of :c:data:`psi_bilinear_matrix` 
 
      :math:`D_\beta^\dagger.C^\dagger.D_\alpha` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format.
 
 
 
@@ -894,11 +894,11 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Transpose of psi_bilinear_matrix 
+    Transpose of :c:data:`psi_bilinear_matrix` 
 
      :math:`D_\beta^\dagger.C^\dagger.D_\alpha` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format.
 
 
 
@@ -914,11 +914,11 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Transpose of psi_bilinear_matrix 
+    Transpose of :c:data:`psi_bilinear_matrix` 
 
      :math:`D_\beta^\dagger.C^\dagger.D_\alpha` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format.
 
 
 
@@ -931,7 +931,7 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Location of the columns in the psi_bilinear_matrix
+    Location of the columns in the :c:data:`psi_bilinear_matrix`
 
 
 
@@ -947,11 +947,11 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Transpose of psi_bilinear_matrix 
+    Transpose of :c:data:`psi_bilinear_matrix` 
 
      :math:`D_\beta^\dagger.C^\dagger.D_\alpha` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` , but the matrix is stored in row major format.
 
 
 
@@ -967,9 +967,9 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Sparse coefficient matrix if the wave function is expressed in a bilinear form : D_a^t C D_b 
+    Sparse coefficient matrix if the wave function is expressed in a bilinear form :  :math:`D_\alpha^\dagger.C.D_\beta` 
 
-    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta.` 
+    Rows are  :math:`\alpha`  determinants and columns are  :math:`\beta` . 
 
     Order refers to psi_det
 
@@ -1187,7 +1187,7 @@ Providers
         integer(bit_kind), allocatable	:: psi_det_alpha_unique	(N_int,psi_det_size)
         integer	:: n_det_alpha_unique
 
-    File: :file:`spindeterminants.irp.f_template_143`
+    File: :file:`spindeterminants.irp.f_template_144`
 
     Unique  :math:`\alpha`  determinants
 
@@ -1214,7 +1214,7 @@ Providers
         integer(bit_kind), allocatable	:: psi_det_beta_unique	(N_int,psi_det_size)
         integer	:: n_det_beta_unique
 
-    File: :file:`spindeterminants.irp.f_template_143`
+    File: :file:`spindeterminants.irp.f_template_144`
 
     Unique  :math:`\beta`  determinants
 
@@ -1582,7 +1582,7 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Dimension of the singles_alpha array
+    Indices of all single excitations
 
 
 
@@ -1596,7 +1596,9 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Dimension of the singles_alpha array
+    singles_alpha_csc_size : Dimension of the :c:data:`singles_alpha_csc` array 
+
+    singles_alpha_csc_idx  : Index where the single excitations of determinant i start
 
 
 
@@ -1610,20 +1612,54 @@ Providers
 
     File: :file:`spindeterminants.irp.f`
 
-    Dimension of the singles_alpha array
+    singles_alpha_csc_size : Dimension of the :c:data:`singles_alpha_csc` array 
+
+    singles_alpha_csc_idx  : Index where the single excitations of determinant i start
 
 
 
 
-.. c:var:: singles_alpha_size
+.. c:var:: singles_beta_csc
 
     .. code:: text
 
-        integer	:: singles_alpha_size
+        integer, allocatable	:: singles_beta_csc	(singles_beta_csc_size)
 
     File: :file:`spindeterminants.irp.f`
 
-    Dimension of the singles_alpha array
+    Indices of all single excitations
+
+
+
+
+.. c:var:: singles_beta_csc_idx
+
+    .. code:: text
+
+        integer*8, allocatable	:: singles_beta_csc_idx	(N_det_beta_unique+1)
+        integer*8	:: singles_beta_csc_size
+
+    File: :file:`spindeterminants.irp.f`
+
+    singles_beta_csc_size : Dimension of the :c:data:`singles_beta_csc` array 
+
+    singles_beta_csc_idx  : Index where the single excitations of determinant i start
+
+
+
+
+.. c:var:: singles_beta_csc_size
+
+    .. code:: text
+
+        integer*8, allocatable	:: singles_beta_csc_idx	(N_det_beta_unique+1)
+        integer*8	:: singles_beta_csc_size
+
+    File: :file:`spindeterminants.irp.f`
+
+    singles_beta_csc_size : Dimension of the :c:data:`singles_beta_csc` array 
+
+    singles_beta_csc_idx  : Index where the single excitations of determinant i start
 
 
 
@@ -1901,7 +1937,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Overwrites psi_det and psi_coef with the wf in bilinear order
+    Overwrites :c:data:`psi_det` and :c:data:`psi_coef` with the wave function in bilinear order
 
 
 
@@ -1957,7 +1993,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Generate a wave function containing all possible products of  :math:`\alpha`  and  :math:`\beta`  determinants
+    Generates a wave function containing all possible products of  :math:`\alpha`  and  :math:`\beta`  determinants
 
 
 
@@ -2261,7 +2297,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Create a wave function from all possible  :math:`\alpha \times \beta`  determinants
+    Creates a wave function from all possible  :math:`\alpha \times \beta`  determinants
 
 
 
@@ -2309,7 +2345,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2327,7 +2363,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2345,7 +2381,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_4(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2363,7 +2399,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_N_int(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2417,7 +2453,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2435,7 +2471,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2453,7 +2489,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_4(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2477,7 +2513,7 @@ Subroutines / functions
 
     Returns the indices of all the single and double excitations in the list of unique  :math:`\alpha`  determinants. 
 
-    /!\ : The buffer is transposed ! 
+    Warning: The buffer is transposed. 
 
 
 
@@ -2511,7 +2547,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2531,7 +2567,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2551,7 +2587,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_and_doubles_4(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2571,7 +2607,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_and_doubles_N_int(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2591,7 +2627,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_N_int(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1221`
+    File: :file:`spindeterminants.irp.f_template_1291`
 
     
 
@@ -2765,7 +2801,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Returns the index of the determinant in the ``psi_det_alpha_unique`` array
+    Returns the index of the determinant in the :c:data:`psi_det_alpha_unique` array
 
 
 
@@ -2779,7 +2815,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Returns the index of the determinant in the ``psi_det_beta_unique`` array
+    Returns the index of the determinant in the :c:data:`psi_det_beta_unique` array
 
 
 
@@ -3545,7 +3581,7 @@ Subroutines / functions
 
     File: :file:`spindeterminants.irp.f`
 
-    Return an integer(8) corresponding to a determinant index for searching
+    Returns an integer(8) corresponding to a determinant index for searching
 
 
 
diff --git a/docs/source/modules/kohn_sham.rst b/docs/source/modules/kohn_sham.rst
index 032455e4..21ad409d 100644
--- a/docs/source/modules/kohn_sham.rst
+++ b/docs/source/modules/kohn_sham.rst
@@ -52,7 +52,6 @@ Providers
 ---------
 
 
-<<<<<<< HEAD
 .. c:var:: ao_potential_alpha_xc
 
     .. code:: text
diff --git a/docs/source/modules/kohn_sham_rs.rst b/docs/source/modules/kohn_sham_rs.rst
index 4976fc11..544bb544 100644
--- a/docs/source/modules/kohn_sham_rs.rst
+++ b/docs/source/modules/kohn_sham_rs.rst
@@ -60,122 +60,6 @@ Providers
 ---------
 
 
-.. c:var:: ao_potential_alpha_xc
-
-    .. code:: text
-
-        double precision, allocatable	:: ao_potential_alpha_xc	(ao_num,ao_num)
-        double precision, allocatable	:: ao_potential_beta_xc	(ao_num,ao_num)
-
-    File: :file:`pot_functionals.irp.f`
-
-    
-
-
-
-
-.. c:var:: ao_potential_beta_xc
-
-    .. code:: text
-
-        double precision, allocatable	:: ao_potential_alpha_xc	(ao_num,ao_num)
-        double precision, allocatable	:: ao_potential_beta_xc	(ao_num,ao_num)
-
-    File: :file:`pot_functionals.irp.f`
-
-    
-
-
-
-
-.. c:var:: e_correlation_dft
-
-    .. code:: text
-
-        double precision	:: e_correlation_dft
-
-    File: :file:`pot_functionals.irp.f`
-
-    
-
-
-
-
-.. c:var:: e_exchange_dft
-
-    .. code:: text
-
-        double precision	:: e_exchange_dft
-
-    File: :file:`pot_functionals.irp.f`
-
-    
-
-
-
-
-.. c:var:: fock_matrix_alpha_no_xc_ao
-
-    .. code:: text
-
-        double precision, allocatable	:: fock_matrix_alpha_no_xc_ao	(ao_num,ao_num)
-        double precision, allocatable	:: fock_matrix_beta_no_xc_ao	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_rs_ks.irp.f`
-
-    Mono electronic an Coulomb matrix in AO basis set
-
-
-
-
-.. c:var:: fock_matrix_beta_no_xc_ao
-
-    .. code:: text
-
-        double precision, allocatable	:: fock_matrix_alpha_no_xc_ao	(ao_num,ao_num)
-        double precision, allocatable	:: fock_matrix_beta_no_xc_ao	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_rs_ks.irp.f`
-
-    Mono electronic an Coulomb matrix in AO basis set
-
-
-
-
-.. c:var:: fock_matrix_energy
-
-    .. code:: text
-
-        double precision	:: rs_ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`rs_ks_energy.irp.f`
-
-    Range-separated Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
-.. c:var:: one_electron_energy
-
-    .. code:: text
-
-        double precision	:: rs_ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`rs_ks_energy.irp.f`
-
-    Range-separated Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
 .. c:var:: rs_ks_energy
 
     .. code:: text
@@ -193,59 +77,11 @@ Providers
 
 
 
-.. c:var:: trace_potential_xc
-
-    .. code:: text
-
-        double precision	:: rs_ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`rs_ks_energy.irp.f`
-
-    Range-separated Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
-.. c:var:: two_electron_energy
-
-    .. code:: text
-
-        double precision	:: rs_ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`rs_ks_energy.irp.f`
-
-    Range-separated Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
 Subroutines / functions
 -----------------------
 
 
 
-.. c:function:: check_coherence_functional
-
-    .. code:: text
-
-        subroutine check_coherence_functional
-
-    File: :file:`rs_ks_scf.irp.f`
-
-    
-
-
-
-
-
 .. c:function:: rs_ks_scf
 
     .. code:: text
diff --git a/docs/source/modules/mo_basis.rst b/docs/source/modules/mo_basis.rst
index f7bf538c..044702fe 100644
--- a/docs/source/modules/mo_basis.rst
+++ b/docs/source/modules/mo_basis.rst
@@ -329,20 +329,6 @@ Subroutines / functions
 
 
 
-.. c:function:: reorder_core_orb
-
-    .. code:: text
-
-        subroutine reorder_core_orb
-
-    File: :file:`track_orb.irp.f`
-
-    routines that takes the current :c:data:`mo_coef` and reorder the core orbitals (see :c:data:`list_core` and :c:data:`n_core_orb`) according to the overlap with :c:data:`mo_coef_begin_iteration`
-
-
-
-
-
 .. c:function:: reorder_core_orb
 
     .. code:: text
diff --git a/docs/source/modules/tools.rst b/docs/source/modules/tools.rst
index 9305a423..af5ebbfe 100644
--- a/docs/source/modules/tools.rst
+++ b/docs/source/modules/tools.rst
@@ -25,7 +25,11 @@ Subroutines / functions
 
     File: :file:`diagonalize_h.irp.f`
 
-    program that extracts the `n_states` lowest states of the Hamiltonian within the set of Slater determinants stored in the EZFIO folder
+    Program that extracts the :option:`determinants n_states` lowest states of the Hamiltonian within the set of Slater determinants stored in the EZFIO folder. 
+
+    If :option:`determinants s2_eig` = True, it will retain only states 
+
+    which corresponds to the desired value of :option:`determinants expected_s2`.
 
 
 
@@ -39,9 +43,9 @@ Subroutines / functions
 
     File: :file:`fcidump.irp.f`
 
-    Produce a regular FCIDUMP file from the MOs stored in the EZFIO folder. 
+    Produce a regular FCIDUMP file from the |MOs| stored in the |EZFIO| folder. 
 
-    To specify an active space, the class of the mos have to set in the EZFIO folder (see set_mo_class/qp_set_mo_class). 
+    To specify an active space, the class of the mos have to set in the |EZFIO| folder (see :ref:`qp_set_mo_class`). 
 
     The fcidump program supports 3 types of MO_class : 
 
@@ -63,9 +67,9 @@ Subroutines / functions
 
     File: :file:`four_idx_transform.irp.f`
 
-    4-index transformation of two-electron integrals from AO to MO integrals. 
+    4-index transformation of two-electron integrals from |AO| to |MO| integrals. 
 
-    This program will compute the two-electron integrals on the MO basis and store it into the EZFIO folder. 
+    This program will compute the two-electron integrals on the |MO| basis and store it into the |EZFIO| folder. 
 
     This program can be useful if the AO --> MO transformation is an expensive step by itself.
 
@@ -99,9 +103,9 @@ Subroutines / functions
 
     for all istate, this program produces 
 
-    * a file "EZFIO. :math:`istate.conv" containing the variational and var+PT2 energies as a function of N_det 
+    * a file "EZFIO.istate.conv" containing the variational and var+PT2 energies as a function of N_det 
 
-    * for istate > 1, a file EZFIO.` istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
+    * for istate > 1, a file EZFIO.istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
 
 
 
@@ -115,9 +119,9 @@ Subroutines / functions
 
     File: :file:`print_wf.irp.f`
 
-    print the wave function stored in the EZFIO folder in the intermediate normalization. 
+    Print the ground state wave function stored in the |EZFIO| folder in the intermediate normalization. 
 
-    it also prints a lot of information regarding the excitation operators from the reference determinant 
+    It also prints a lot of information regarding the excitation operators from the reference determinant 
 
     and a first-order perturbative analysis of the wave function. 
 
@@ -141,6 +145,34 @@ Subroutines / functions
 
 
 
+.. c:function:: routine_e_conv
+
+    .. code:: text
+
+        subroutine routine_e_conv
+
+    File: :file:`print_e_conv.irp.f`
+
+    routine called by :c:func:`print_e_conv`
+
+
+
+
+
+.. c:function:: routine_save_one_e_dm
+
+    .. code:: text
+
+        subroutine routine_save_one_e_dm
+
+    File: :file:`save_one_e_dm.irp.f`
+
+    routine called by :c:func:`save_one_e_dm`
+
+
+
+
+
 .. c:function:: save_natorb
 
     .. code:: text
diff --git a/docs/source/programmers_guide/index_providers.rst b/docs/source/programmers_guide/index_providers.rst
index b787d0bd..4cadc9ac 100644
--- a/docs/source/programmers_guide/index_providers.rst
+++ b/docs/source/programmers_guide/index_providers.rst
@@ -672,7 +672,9 @@ Index of Providers
 * :c:data:`singles_alpha_csc`
 * :c:data:`singles_alpha_csc_idx`
 * :c:data:`singles_alpha_csc_size`
-* :c:data:`singles_alpha_size`
+* :c:data:`singles_beta_csc`
+* :c:data:`singles_beta_csc_idx`
+* :c:data:`singles_beta_csc_size`
 * :c:data:`size_select_max`
 * :c:data:`slater_bragg_radii`
 * :c:data:`slater_bragg_radii_per_atom`
@@ -950,9 +952,6 @@ Index of Subroutines/Functions
 * :c:func:`get_index_in_psi_det_beta_unique`
 * :c:func:`get_index_in_psi_det_sorted_bit`
 * :c:func:`get_inverse`
-* :c:func:`get_m0`
-* :c:func:`get_m1`
-* :c:func:`get_m2`
 * :c:func:`get_mask_phase`
 * :c:func:`get_mo_erf_map_size`
 * :c:func:`get_mo_map_size`
@@ -1233,7 +1232,9 @@ Index of Subroutines/Functions
 * :c:func:`rintgauss`
 * :c:func:`roothaan_hall_scf`
 * :c:func:`routine`
+* :c:func:`routine_e_conv`
 * :c:func:`routine_example_psi_det`
+* :c:func:`routine_save_one_e_dm`
 * :c:func:`rs_ks_scf`
 * :c:func:`run`
 * :c:func:`run_cipsi`
@@ -1391,4 +1392,4 @@ Index of Subroutines/Functions
 * :c:func:`zmq_put_psi_det_beta_unique`
 * :c:func:`zmq_put_psi_det_size`
 * :c:func:`zmq_selection`
-* :c:func:`zmq_set_running`
+* :c:func:`zmq_set_running`
\ No newline at end of file
diff --git a/docs/source/users_guide/natural_orbitals.rst b/docs/source/users_guide/natural_orbitals.rst
index 99ef2c1c..12c4fd03 100644
--- a/docs/source/users_guide/natural_orbitals.rst
+++ b/docs/source/users_guide/natural_orbitals.rst
@@ -34,7 +34,7 @@ These orbitals are in general known to be better than the usual |HF| |MOs| as th
 
 .. code::
 
-    qp_run save_natorb hcn
+    qp_run save_natorb file.ezfio
 
 
 Hands on
@@ -50,5 +50,5 @@ Hands on
 
 .. seealso:: 
 
-    The documentation of the :ref:`save_natorb` program. 
+    The documentation of the :c:func:`save_natorb` program. 
  
diff --git a/docs/source/users_guide/printing.rst b/docs/source/users_guide/printing.rst
index 386f7e22..4f94334d 100644
--- a/docs/source/users_guide/printing.rst
+++ b/docs/source/users_guide/printing.rst
@@ -7,7 +7,7 @@ If the wave function contains less than :math:`10^4` determinants, you can direc
 
 .. code::
 
-    qp_edit hcn
+    qp_edit file.ezfio
 
 
 .. important::
@@ -22,7 +22,7 @@ interactively in :ref:`qp_edit` mode. An alternative is to use the :command:`pri
 
 .. code::
 
-    qp_run print_wf hcn | tee hcn.fci_natorb.wf
+    qp_run print_wf file.ezfio | tee file.ezfio.fci_natorb.wf
 
 This program will, by default, print out the first :math:`10^4` determinants whatever the size of the wave function stored in the |EZFIO| folder. If you want to change the number of printed Slater determinants, just change the :option:`determinants n_det_print_wf` keyword using the :ref:`qp_edit` tool. 
 
diff --git a/docs/source/users_guide/qp_run.rst b/docs/source/users_guide/qp_run.rst
index 31d010cc..832616c6 100644
--- a/docs/source/users_guide/qp_run.rst
+++ b/docs/source/users_guide/qp_run.rst
@@ -14,17 +14,27 @@ Usage
 
 .. code:: bash
 
-   qp_run [-slave] [-help] <PROGRAM> <EZFIO_DIRECTORY>
+  qp_run [-h|--help] [-p <string>|--prefix=<string>] [-s|--slave] [--]
+      PROGRAM EZFIO_DIR
+
+``PROGRAM`` is the name of the |QP| program to be run, and ``EZFIO_DIR`` is
+the name of the |EZFIO| directory containing the data.
 
 
-.. option:: -help
+.. option:: -h, --help
 
    Displays the list of available |qp| programs. 
 
 
-.. option:: -slave
+.. option:: -p <string>, --prefix=<string>
 
-   This option needs to be set when ``PROGRAM`` is a slave program, to accelerate
+   Prefix before running the program. This option is used to run programs like
+   like gdb or valgrind.
+
+
+.. option:: -s, --slave
+
+   This option needs to be set to run a slave job for ``PROGRAM``, to accelerate
    another running instance of the |qp|.
 
 
@@ -33,6 +43,8 @@ Example
 
 .. code:: bash
 
-   qp_run FCI h2o.ezfio
+   qp_run fci h2o.ezfio &
+   srun qp_run --slave fci h2o.ezfio
+   wait
 
 
diff --git a/etc/ezfio.rc b/etc/ezfio.rc
index 67a8a008..ec557139 100644
--- a/etc/ezfio.rc
+++ b/etc/ezfio.rc
@@ -2,11 +2,15 @@
 
 export QP_EZFIO=${QP_ROOT}/external/ezfio
 
-if [[ -f ${QP_EZFIO}/Bash/ezfio.sh ]]; then
-      if [[ "$(ps -p $$ -ocomm=)" == "zsh" ]] ; then
-        autoload bashcompinit
-        bashcompinit
-      fi
-      source ${QP_EZFIO}/Bash/ezfio.sh
-fi
+function source_if_exists() {
+  if [[ -f $1 ]]; then
+        if [[ "$(ps -p $$ -ocomm=)" == "zsh" ]] ; then
+          autoload bashcompinit
+          bashcompinit
+        fi
+        source $1
+  fi
+} 
+
+source_if_exists "${QP_EZFIO}/Bash/ezfio.sh"
 
diff --git a/etc/irpf90.rc b/etc/irpf90.rc
index e3a57d1d..cbf5552b 100644
--- a/etc/irpf90.rc
+++ b/etc/irpf90.rc
@@ -1,9 +1,21 @@
 # Configuration of IRPF90 package
 
-export IRPF90=${QP_ROOT}/bin/irpf90
+# Set the path of IRPF90 here:
+export IRPF90_PATH=${QP_ROOT}/external/irpf90-v1.7.4
+export PATH=${PATH}:${IRPF90_PATH}/bin
 
-# Load irpman shell completion
-irpman=$(tail -1 "${QP_ROOT}/bin/irpman" | cut --delimiter " " --field=2)
-source $(dirname ${irpman})/../irpman-completions.bash
+export IRPF90=${IRPF90_PATH}/bin/irpf90
+
+function source_if_exists() {
+  if [[ -f $1 ]]; then
+        if [[ "$(ps -p $$ -ocomm=)" == "zsh" ]] ; then
+          autoload bashcompinit
+          bashcompinit
+        fi
+        source $1
+  fi
+} 
+
+source_if_exists "${IRPF90_PATH}/irpman-completions.bash"
   
 
diff --git a/etc/qp.rc b/etc/qp.rc
index 8270a73e..ad34b786 100644
--- a/etc/qp.rc
+++ b/etc/qp.rc
@@ -24,8 +24,9 @@ Usage:
   qp mpirun      PROGRAM
 
   qp set_frozen_core
-  qp create_ezfio_from_xyz -help
-  qp set_mo_class          -help
+  qp create_ezfio_from_xyz    --help
+  qp convert_output_to_ezfio  --help
+  qp set_mo_class             --help
 "
 }
 
@@ -41,17 +42,17 @@ function qp()
 {
   case $1 in
     "has"|"set"|"get"|"set_file"|"unset_file")
-        ezfio $@
+        ezfio "$@"
       ;;
 
     "set_frozen_core")
       shift
-      qp_set_frozen_core ${EZFIO_FILE}
+      qp_set_frozen_core "$@" ${EZFIO_FILE}
       ;;
 
     "create_ezfio_from_xyz")
       shift
-      NAME=$(qp_create_ezfio_from_xyz $@)
+      NAME=$(qp_create_ezfio_from_xyz "$@")
       if [[ -d $NAME ]] ; then
         [[ -d $EZFIO_FILE ]] && ezfio unset_file
         ezfio set_file $NAME
@@ -60,29 +61,34 @@ function qp()
       fi
       ;;
 
+    "convert_output_to_ezfio")
+      shift
+      qp_convert_output_to_ezfio "$@"
+      ;;
+
     "set_mo_class")
       shift
-      qp_set_mo_class ${EZFIO_FILE} $@
+      qp_set_mo_class "$@" -- ${EZFIO_FILE} 
       ;;
 
     "edit")
       shift
-      qp_edit ${EZFIO_FILE}
+      qp_edit "$@" -- ${EZFIO_FILE} 
       ;;
 
     "run")
       shift
-      qp_run $@ ${EZFIO_FILE}
+      qp_run "$@" -- ${EZFIO_FILE}
       ;;
 
     "srun")
       shift
-      qp_srun $@ ${EZFIO_FILE}
+      qp_srun "$@" ${EZFIO_FILE}
       ;;
 
     "mpirun")
       shift
-      qp_mpirun $@ ${EZFIO_FILE}
+      qp_mpirun "$@" ${EZFIO_FILE}
       ;;
 
     "man")
@@ -90,6 +96,19 @@ function qp()
       man $QP_ROOT/man/${1}.?
       ;;
 
+    "prompt")
+      shift
+      function _check_ezfio() {
+        if [[ -d ${EZFIO_FILE} ]] ; then
+          printf "\e[0;32m|${EZFIO_FILE}>\e[m"
+        else
+          printf "\e[0;31m|${EZFIO_FILE}>\e[m"
+        fi
+      }
+#      PS1="${PS1%\\\$ } \$(_check_ezfio) $ "
+      PS1="\$(_check_ezfio)\n$PS1"
+      ;;
+      
     *)
       _qp_usage
       ;;
@@ -150,6 +169,7 @@ _Complete()
         COMPREPLY=( $(compgen -W 'has get set unset_file edit \
                                   run srun mpirun set_frozen_core \
                                   set_mo_class create_ezfio_from_xyz \
+                                  convert_output_to_ezfio \
                                   -h' -- $cur ) )
         return 0
         ;;
@@ -166,7 +186,7 @@ _Complete()
         COMPREPLY=( $(compgen -W "$(\ls -d */ | sed 's|/||g')" -- ${cur} ) )
         return 0
         ;;
-      create_ezfio_from_xyz)
+      convert_output_to_ezfio|create_ezfio_from_xyz)
         COMPREPLY=( $(compgen -W "$(\ls)" -- ${cur} ) )
         return 0
         ;;
@@ -174,6 +194,7 @@ _Complete()
         COMPREPLY=( $(compgen -W 'set_file \
                                   man \
                                   create_ezfio_from_xyz \
+                                  convert_output_to_ezfio \
                                   -h' -- $cur ) )
         return 0
         ;;
diff --git a/man/configure.1 b/man/configure.1
index 643db372..60341d38 100644
--- a/man/configure.1
+++ b/man/configure.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "CONFIGURE" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "CONFIGURE" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 configure \-  | Quantum Package >
 .
@@ -30,6 +30,8 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .in \\n[rst2man-indent\\n[rst2man-indent-level]]u
 ..
+.sp
+Program that can either configure the compilations options and download/install external dependencies (see the installation description).
 .SH USAGE
 .INDENT 0.0
 .INDENT 3.5
diff --git a/man/excited_states.1 b/man/excited_states.1
index 0b76385f..c0414666 100644
--- a/man/excited_states.1
+++ b/man/excited_states.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "EXCITED_STATES" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "EXCITED_STATES" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 excited_states \-  | Quantum Package >
 .
@@ -42,7 +42,7 @@ s2_eig\fP flag \fBtrue\fP\&. This will force the Davidson algorithm to choose on
 vectors with a value of \ewidehat{S^2} equal to \fBdeterminants expected_s2\fP\&.
 Otherwise, different spin states will come out in the diagonalization.
 .sp
-The \fIQuantum Package\fP doesn\(aqt take account of the symmetry. Due to numerical noise, excited
+The \fIQuantum Package\fP doesn’t take account of the symmetry. Due to numerical noise, excited
 states of different symmetries may enter in the calculation. Note that it is
 possible to make state\-average calculation of states with different symmetries
 and/or different spin multiplicities.
@@ -55,6 +55,13 @@ and then running the selected FCI restarting from the CIS states, setting
 Usually, it is good practice to use state\-averaged natural MOs so that all
 states have MOs of comparable quality. This allows for a faster convergence
 of excitation energies.
+.sp
+\fBSEE ALSO:\fP
+.INDENT 0.0
+.INDENT 3.5
+The documentation of the \fBscf()\fP, \fBcis()\fP and  \fBfci()\fP programs.
+.UNINDENT
+.UNINDENT
 .SH AUTHOR
 A. Scemama, E. Giner
 .SH COPYRIGHT
diff --git a/man/interfaces.1 b/man/interfaces.1
index 91e362f8..6708c93c 100644
--- a/man/interfaces.1
+++ b/man/interfaces.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "INTERFACES" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "INTERFACES" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 interfaces \-  | Quantum Package >
 .
diff --git a/man/natural_orbitals.1 b/man/natural_orbitals.1
index 488a671b..c3f5fcaf 100644
--- a/man/natural_orbitals.1
+++ b/man/natural_orbitals.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "NATURAL_ORBITALS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "NATURAL_ORBITALS" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 natural_orbitals \-  | Quantum Package >
 .
@@ -30,6 +30,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .in \\n[rst2man-indent\\n[rst2man-indent-level]]u
 ..
+.SH SUMMARY
 .sp
 To produce state\-average natural orbitals, run
 .INDENT 0.0
@@ -45,6 +46,59 @@ qp_run save_natorb file.ezfio
 .sp
 The MOs will be replaced, so the two\-electron integrals and the wave function
 are invalidated as well.
+.SH EXTRACTING NATURAL ORBITALS
+.sp
+Once obtained the near FCI wave function, one can obtain many quantities related to it.
+One of these quantities are the natural orbitals which have the properties of diagonalizing the one\-body density matrix:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.ce
+\erho_{ij} = \edelta_{ij}
+.ce 0
+.UNINDENT
+.UNINDENT
+.sp
+where the element of the one\-body density matrix \erho_{ij} is define as:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.ce
+\erho_{ij} = \elangle \ePsi | \eleft( a^{\edagger}_{j,\ealpha} a_{i,\ealpha} + a^{\edagger}_{j,\ebeta} a_{i,\ebeta} \eright) |\ \ePsi\ \erangle
+.ce 0
+.UNINDENT
+.UNINDENT
+.sp
+These orbitals are in general known to be better than the usual HF MOs as they are obtained from a correlated wave function. To use these orbitals for future calculations, one has to replace the current MOs by the natural orbitals. To do so, just run:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_run save_natorb file.ezfio
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH HANDS ON
+.sp
+\fBIMPORTANT:\fP
+.INDENT 0.0
+.INDENT 3.5
+As the MOs are changed, for the sake of coherence of future calculations, the save_natorb program
+\fIautomatically removes the current wave function\fP stored in the \fI\%EZFIO\fP database and replace
+it by a single Slater determinant corresponding to a HF occupation of the new spin orbitals.
+Also, all the keywords to read the one\- and two\-electron integrals on the MO basis are set to “None”
+in order to be sure not to read integrals not corresponding to the current set of MOs\&.
+.UNINDENT
+.UNINDENT
+.sp
+\fBSEE ALSO:\fP
+.INDENT 0.0
+.INDENT 3.5
+The documentation of the \fBsave_natorb()\fP program.
+.UNINDENT
+.UNINDENT
 .SH AUTHOR
 A. Scemama, E. Giner
 .SH COPYRIGHT
diff --git a/man/plugins.1 b/man/plugins.1
index cc3dad9f..389f0bd8 100644
--- a/man/plugins.1
+++ b/man/plugins.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "PLUGINS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "PLUGINS" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 plugins \-  | Quantum Package >
 .
@@ -35,7 +35,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 plugins need to be downloaded and installed in the \fB$QP_ROOT/plugins\fP
 directory.
 .sp
-Plugins are usually hosted in external repositories. To dowload a plugin,
+Plugins are usually hosted in external repositories. To download a plugin,
 the remote repository needs to be downloaded, and the plugins of the
 repository can be selected for installation.
 .sp
@@ -93,6 +93,13 @@ qp_plugins uninstall ext_module
 .fi
 .UNINDENT
 .UNINDENT
+.sp
+\fBSEE ALSO:\fP
+.INDENT 0.0
+.INDENT 3.5
+For a more detailed explanation and an example, see qp_plugins\&.
+.UNINDENT
+.UNINDENT
 .SH AUTHOR
 A. Scemama, E. Giner
 .SH COPYRIGHT
diff --git a/man/qp_convert_output_to_ezfio.1 b/man/qp_convert_output_to_ezfio.1
index dbd78656..99f75dab 100644
--- a/man/qp_convert_output_to_ezfio.1
+++ b/man/qp_convert_output_to_ezfio.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_CONVERT_OUTPUT_TO_EZFIO" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_CONVERT_OUTPUT_TO_EZFIO" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_convert_output_to_ezfio \-  | Quantum Package >
 .
diff --git a/man/qp_create_ezfio_from_xyz.1 b/man/qp_create_ezfio_from_xyz.1
index 2b7edc51..1d33db40 100644
--- a/man/qp_create_ezfio_from_xyz.1
+++ b/man/qp_create_ezfio_from_xyz.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_CREATE_EZFIO_FROM_XYZ" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_CREATE_EZFIO_FROM_XYZ" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_create_ezfio_from_xyz \-  | Quantum Package >
 .
@@ -80,7 +80,7 @@ Total charge of the molecule. Default is 0.
 .INDENT 0.0
 .TP
 .B \-cart
-Compute AOs in the Cartesian basis set (6d, 10f, ...)
+Compute AOs in the Cartesian basis set (6d, 10f, …)
 .UNINDENT
 .INDENT 0.0
 .TP
diff --git a/man/qp_edit.1 b/man/qp_edit.1
index b1eeff37..295cea94 100644
--- a/man/qp_edit.1
+++ b/man/qp_edit.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_EDIT" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_EDIT" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_edit \-  | Quantum Package >
 .
diff --git a/man/qp_export_as_tgz.1 b/man/qp_export_as_tgz.1
index 897bcae9..e494c16c 100644
--- a/man/qp_export_as_tgz.1
+++ b/man/qp_export_as_tgz.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_EXPORT_AS_TGZ" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_EXPORT_AS_TGZ" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_export_as_tgz \-  | Quantum Package >
 .
@@ -33,7 +33,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .sp
 In some HPC facilities, the access to the internet is limited for security reasons.
 In such an environment, the installation of \fIQuantum Package\fP is sometimes very painful because
-the OCaml compiler and the libraries can\(aqt be installed by a non\-root user.
+the OCaml compiler and the libraries can’t be installed by a non\-root user.
 .sp
 This command creates a self\-contained binary distribution in the form of a \fItar.gz\fP file
 that can be copied on another machine.
diff --git a/man/qp_plugins.1 b/man/qp_plugins.1
index 7f7f5a4f..019d7c62 100644
--- a/man/qp_plugins.1
+++ b/man/qp_plugins.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_PLUGINS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_PLUGINS" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_plugins \-  | Quantum Package >
 .
@@ -31,7 +31,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
 .in \\n[rst2man-indent\\n[rst2man-indent-level]]u
 ..
 .sp
-This command deals with all external plugings of \fIQuantum Package\fP\&. Plugin repositories can
+This command deals with all external plugins of \fIQuantum Package\fP\&. Plugin repositories can
 be downloaded, and the plugins in these repositories can be
 installed/uninstalled of created.
 .SH USAGE
@@ -101,6 +101,69 @@ Create a new plugin named \fBplugin_name\fP (in local repository by default).
 .B \-r <repository>
 Specify in which repository the new plugin will be created.
 .UNINDENT
+.SH EXAMPLE
+.sp
+Let us download, install and compile some specific external plugins from \fI\%https://gitlab.com/eginer/qp_plugins_eginer\fP .
+.sp
+First, download the git repo associated to these plugins. To do so, first go to the \fIplugins\fP directory in the \fIQuantum Package\fP and execute:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins download https://gitlab.com/eginer/qp_plugins_eginer
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+This will create in the directory \fIplugins\fP a local copy of the git repo located at the URL you indicated.
+Then, go in \fIqp_plugins_eginer/stable/\fP
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+cd qp_plugins_eginer/stable/
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+In the directory \fIstable\fP, there are many directories which all correspond to a specific plugin that have been developed
+by the person in charge of the git repo. All these plugins might use some global variables and routines contained
+in the core modules of the \fIQuantum Package\fP\&.
+.sp
+Now let us install the plugin \fIrsdft_cipsi\fP:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins install rsdft_cipsi
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+This will link this directory to the \fIQuantum Package\fP which means that when the code will be compiled, this plugin will be compiled to and therefore
+all the executables/scripts/input keywords contained in this module will be available as if there were part of the core of the \fIQuantum Package\fP\&.
+.sp
+Then, to compile the new plugin, just recompile the \fIQuantum Package\fP as usual by going at the root of the \fIQuantum Package\fP directory:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+cd $QP_ROOT
+ninja
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+Then, if you go back to the plugin directory you just installed, you should see all the executables/scripts which have been created and
+which are now available with the \fIqp_run\fP command.
 .SH AUTHOR
 A. Scemama, E. Giner
 .SH COPYRIGHT
diff --git a/man/qp_run.1 b/man/qp_run.1
index 1a747f00..95c070c4 100644
--- a/man/qp_run.1
+++ b/man/qp_run.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_RUN" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_RUN" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_run \-  | Quantum Package >
 .
@@ -38,20 +38,30 @@ Command used to run a calculation.
 .sp
 .nf
 .ft C
-qp_run [\-slave] [\-help] <PROGRAM> <EZFIO_DIRECTORY>
+qp_run [\-h|\-\-help] [\-p <string>|\-\-prefix=<string>] [\-s|\-\-slave] [\-\-]
+    PROGRAM EZFIO_DIR
 .ft P
 .fi
 .UNINDENT
 .UNINDENT
+.sp
+\fBPROGRAM\fP is the name of the \fIQuantum Package\fP program to be run, and \fBEZFIO_DIR\fP is
+the name of the \fI\%EZFIO\fP directory containing the data.
 .INDENT 0.0
 .TP
-.B \-help
+.B \-h, \-\-help
 Displays the list of available \fIQuantum Package\fP programs.
 .UNINDENT
 .INDENT 0.0
 .TP
-.B \-slave
-This option needs to be set when \fBPROGRAM\fP is a slave program, to accelerate
+.B \-p <string>, \-\-prefix=<string>
+Prefix before running the program. This option is used to run programs like
+like gdb or valgrind.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-s, \-\-slave
+This option needs to be set to run a slave job for \fBPROGRAM\fP, to accelerate
 another running instance of the \fIQuantum Package\fP\&.
 .UNINDENT
 .SH EXAMPLE
@@ -60,7 +70,9 @@ another running instance of the \fIQuantum Package\fP\&.
 .sp
 .nf
 .ft C
-qp_run FCI h2o.ezfio
+qp_run fci h2o.ezfio &
+srun qp_run \-\-slave fci h2o.ezfio
+wait
 .ft P
 .fi
 .UNINDENT
diff --git a/man/qp_set_frozen_core.1 b/man/qp_set_frozen_core.1
index 0d72b392..41eea988 100644
--- a/man/qp_set_frozen_core.1
+++ b/man/qp_set_frozen_core.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_SET_FROZEN_CORE" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_SET_FROZEN_CORE" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_set_frozen_core \-  | Quantum Package >
 .
diff --git a/man/qp_set_mo_class.1 b/man/qp_set_mo_class.1
index 205e1656..7c8f1245 100644
--- a/man/qp_set_mo_class.1
+++ b/man/qp_set_mo_class.1
@@ -1,6 +1,6 @@
 .\" Man page generated from reStructuredText.
 .
-.TH "QP_SET_MO_CLASS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.TH "QP_SET_MO_CLASS" "1" "Jan 14, 2019" "2.0" "Quantum Package"
 .SH NAME
 qp_set_mo_class \-  | Quantum Package >
 .
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index 80a81503..7fbaff09 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@@ -22,6 +22,7 @@ module Determinants_by_hand : sig
   val read_n_int : unit -> N_int_number.t
   val update_ndet : Det_number.t -> unit
   val extract_state : States_number.t -> unit
+  val extract_states : Range.t -> unit
 end = struct
   type t = 
     { n_int                  : N_int_number.t;
@@ -592,6 +593,47 @@ psi_det                = %s
     write new_det
   ;;
 
+  let extract_states range =
+    Printf.printf "Extracting states %s\n" (Range.to_string range);
+    let det =
+      read ()
+    in
+    let n_det, n_states =
+      Det_number.to_int det.n_det,
+      States_number.to_int det.n_states
+    in
+    Range.to_int_list range
+    |> List.iter ~f:(fun istate ->
+      if istate > n_states then
+        failwith "State to extract should not be greater than n_states") 
+    ;
+    let sorted_list = 
+      Range.to_int_list range
+      |> List.sort ~compare 
+    in
+    let state_shift = ref 0 in
+    List.iter ~f:(fun istate ->
+      let j =
+        istate - 1
+      in
+      begin
+        if (j>0) then
+          let ishift =
+            j*n_det
+          in
+          for i=0 to (n_det-1) do
+            det.psi_coef.(!state_shift+i) <- det.psi_coef.(i+ishift)
+          done
+      end;
+      state_shift := !state_shift + n_det
+    ) sorted_list
+    ;
+    let new_det =
+      { det with n_states = (States_number.of_int @@ List.length sorted_list) }
+    in
+    write new_det
+  ;;
+
 end
 
 
diff --git a/ocaml/Range.ml b/ocaml/Range.ml
index 91fcbcce..d11952b4 100644
--- a/ocaml/Range.ml
+++ b/ocaml/Range.ml
@@ -14,6 +14,8 @@ open Sexplib.Std
 
 type t = int list [@@deriving sexp] 
 
+let to_int_list r = r
+
 let expand_range r =
   match String_ext.lsplit2 ~on:'-' r with
   | Some (s, f) ->
@@ -50,6 +52,10 @@ let of_string s =
 
 
 let to_string l =
+  "[" ^ 
+  (List.map string_of_int l
+   |> String.concat ",") ^ "]"
+(*
   let rec do_work buf symbol = function
     | [] -> buf
     | a::([] as t) -> 
@@ -58,15 +64,15 @@ let to_string l =
         if (b-a = 1) then
           do_work buf "-" t
         else
-          do_work (buf^symbol^(string_of_int a)^","^(string_of_int b)) "" t
+          do_work (buf^symbol^","^(string_of_int b)) "" t
   in
   let result = 
     match l with
-    | [] ->
-        "[]"
+    | [] -> "[]"
     | h::t  ->
         do_work ("["^(string_of_int h)) "" l in
   (String.sub result 0 ((String.length result)))^"]"
+  *)
 
 
 let test_module () =
diff --git a/ocaml/Range.mli b/ocaml/Range.mli
index e186ccf9..2a930d80 100644
--- a/ocaml/Range.mli
+++ b/ocaml/Range.mli
@@ -8,3 +8,4 @@ type t = int list [@@deriving sexp]
  *)
 val of_string : string -> t
 val to_string : t -> string
+val to_int_list : t -> int list
diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml
index d41e6081..a94c16db 100644
--- a/ocaml/qp_run.ml
+++ b/ocaml/qp_run.ml
@@ -2,7 +2,6 @@ open Qputils
 
 (* Environment variables :
 
-   QP_PREFIX=gdb   : to run gdb (or valgrind, or whatever)
    QP_TASK_DEBUG=1 : debug task server
 
 *)
@@ -14,7 +13,7 @@ let print_list () =
 let () = 
   Random.self_init ()
 
-let run slave exe ezfio_file =
+let run slave ?prefix exe ezfio_file =
 
   (** Check availability of the ports *)
   let port_number = 
@@ -113,8 +112,9 @@ let run slave exe ezfio_file =
 
   (** Run executable *)
   let prefix = 
-    try (Sys.getenv "QP_PREFIX")^" " with
-    | Not_found -> ""
+    match prefix with
+    | Some x -> x^" "
+    | None -> ""
   and exe =
     match (List.find (fun (x,_) -> x = exe) executables) with
     | (_,exe) -> exe^" "
@@ -147,11 +147,12 @@ let () =
     |> String.concat "\n"
   ) |> Command_line.set_header_doc;
 
-  [ (
-      's', "slave", "Required to run slave tasks in distributed environments",
+  [ ( 's', "slave", "Required to run slave tasks in distributed environments",
         Command_line.Without_arg);
-      Command_line.anonymous "<EXECUTABLE>" "Name of the QP program to be run";
-      Command_line.anonymous "<EZFIO_FILE>" "EZFIO directory";
+    ( 'p', "prefix", "<string> Prefix before running the program, like gdb or valgrind",
+        Command_line.With_arg );
+      Command_line.anonymous "PROGRAM" "Name of the QP program to be run";
+      Command_line.anonymous "EZFIO_DIR" "EZFIO directory";
     ]
   |> Command_line.set_specs ;
 
@@ -166,9 +167,13 @@ let () =
     | _ -> true
   in
 
+  let prefix = 
+    Command_line.get "prefix"
+  in
+
   (* Run the program *)
   match Command_line.anon_args () with
-  | exe :: ezfio_file :: [] -> run slave exe ezfio_file
+  | exe :: ezfio_file :: [] -> run slave ?prefix exe ezfio_file
   | _ ->  (Command_line.help () ; failwith "Inconsistent command line")
   
 
diff --git a/qpsh b/qpsh
deleted file mode 100755
index 967f9429..00000000
--- a/qpsh
+++ /dev/null
@@ -1,24 +0,0 @@
-#!/bin/bash
-
-export QP_ROOT=$(dirname $0)
-
-exec bash --init-file <(cat << EOF
-
-[[ -f \${HOME}/.bashrc ]] && source \${HOME}/.bashrc 
-
-ESC=\$(printf "\\e")
-
-function check_ezfio() {
-  if [[ -d \${EZFIO_FILE} ]] ; then
-    printf "\\e[0;32m|\${EZFIO_FILE}>\\e[m"
-  else
-    printf "\\e[0;31m|\${EZFIO_FILE}>\\e[m"
-  fi
-}
-
-PS1="\${PS1%$ }\\\$(check_ezfio) \$ "
-source \${QP_ROOT}/quantum_package.rc
-EOF
-)
-
-
diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template
index b3fdde79..c832d1b6 100644
--- a/scripts/ezfio_interface/qp_edit_template
+++ b/scripts/ezfio_interface/qp_edit_template
@@ -169,9 +169,9 @@ let run check_only ?ndet ?state ezfio_filename =
   begin
     match state with
     | None -> ()
-    | Some n -> 
+    | Some range -> 
        begin
-         Input.Determinants_by_hand.extract_state (States_number.of_int n)
+         Input.Determinants_by_hand.extract_states range
        end
   end;
 
@@ -273,7 +273,7 @@ let () =
   [ ( 'c', "check", "Checks the input data", Command_line.Without_arg);
     ( 'n', "ndet", "<int> Truncate the wavefunction to the target number of determinants",
       Command_line.With_arg);
-    ( 's', "state", "<int> Pick the state as a new wavefunction", Command_line.With_arg);
+    ( 's', "state", "<range> Extract selected states, for example \"[1,3-5]\"", Command_line.With_arg);
     Command_line.anonymous "<EZFIO_FILE>" "EZFIO directory";
   ]
   |> Command_line.set_specs ;
@@ -292,8 +292,8 @@ let () =
   let state =
     match Command_line.get "state" with
     | None -> None
-    | Some s -> (try Some (int_of_string s)
-            with _ -> failwith "[-s|--state] expects an integer")
+    | Some s -> (try Some (Range.of_string s)
+            with _ -> failwith "[-s|--state] expects a range")
   in
 
   let c =
diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio b/scripts/qp_convert_output_to_ezfio
similarity index 95%
rename from scripts/ezfio_interface/qp_convert_output_to_ezfio
rename to scripts/qp_convert_output_to_ezfio
index b9d5c704..e1c2add5 100755
--- a/scripts/ezfio_interface/qp_convert_output_to_ezfio
+++ b/scripts/qp_convert_output_to_ezfio
@@ -1,14 +1,13 @@
 #!/usr/bin/env python2
 """
-convert output of gamess/GAU$$IAN to ezfio
+convert output of GAMESS/GAU$$IAN to ezfio
 
 Usage:
-    qp_convert_output_to_ezfio.py <file.out> [-o|--output <EZFIO_DIRECTORY>]
+    qp_convert_output_to_ezfio [-o EZFIO_FILE|--output=EZFIO_FILE] <file.out>
 
-Option:
-    file.out is the file to check (like gamess.out)
-    EZFIO_DIRECTORY is the name of the produced directory
-        (by default is file.out.ezfio)
+Options:
+    -o EZFIO_FILE  --output=EZFIO_FILE    Produced directory
+                                          by default is file.out.ezfio
 
 """
 
@@ -358,10 +357,13 @@ if __name__ == '__main__':
 
     file_ = get_full_path(arguments['<file.out>'])
 
-    if arguments["-o"]:
-        ezfio_file = get_full_path(arguments["<ezfio_directory>"])
+    if arguments["--output"]:
+        ezfio_file = get_full_path(arguments["--output"])
+    elif arguments["-o"]:
+        ezfio_file = get_full_path(arguments["-o"])
     else:
         ezfio_file = "{0}.ezfio".format(file_)
+    print ezfio_file
 
     try:
         res_file = getFile(file_)
@@ -371,10 +373,11 @@ if __name__ == '__main__':
         print file_, 'recognized as', str(res_file).split('.')[-1].split()[0]
 
     write_ezfio(res_file, ezfio_file)
+    sys.stdout.flush()
     if os.system("qp_run save_ortho_mos "+ezfio_file) != 0:
       print """Warning: You need to run 
 
-         qp_run save_ortho_mos """+ezfio_file+"""
+         qp run save_ortho_mos 
 
 to be sure your MOs will be orthogonal, which is not the case when
 the MOs are read from output files (not enough precision in output)."""
diff --git a/scripts/qp_export_as_tgz b/scripts/qp_export_as_tgz
index f440cd1f..695840ef 100755
--- a/scripts/qp_export_as_tgz
+++ b/scripts/qp_export_as_tgz
@@ -18,6 +18,8 @@ cd ${QP_ROOT}
 if [[ -f quantum_package.rc \
    && -f README.md \
    && -d src  \
+   && -d etc  \
+   && -d man  \
    && -d bin  \
    && -d ocaml \
    && -d external \
@@ -43,7 +45,7 @@ QPACKAGE_STATIC=${QP_ROOT}/quantum_package_static
 
 function find_libs ()
 {
-  for i in $@
+  for i in "$@"
   do
       ldd $i
   done | sort | grep '/' | cut --delimiter=' ' --fields=3 | uniq
@@ -61,10 +63,10 @@ echo "Creating root of static directory"
 #     ---------------------------------
 
 rm --recursive --force -- "${QPACKAGE_STATIC}"
-mkdir --parents -- ${QPACKAGE_STATIC}/{bin,lib,extra_lib,external}
+mkdir --parents -- ${QPACKAGE_STATIC}/{bin,etc,man,lib,extra_lib,external}
 if [[ $? -ne 0 ]] ;
 then
-  echo "Error creating ${QPACKAGE_STATIC}/{bin,lib,extra_lib,external}"
+  echo "Error creating ${QPACKAGE_STATIC}/{bin,lib,etc,man,extra_lib,external}"
   exit 1
 fi
 
@@ -102,6 +104,8 @@ cd ${QPACKAGE_STATIC}/bin
   done
 )
 
+cp ${QP_ROOT}/bin/qpsh ${QPACKAGE_STATIC}/bin
+
 cp --recursive -- ${QP_ROOT}/data ${QPACKAGE_STATIC}/data
 for i in ${FORTRAN_EXEC}
 do
@@ -163,13 +167,21 @@ cp --recursive -- ${QP_ROOT}/external/Python ${QPACKAGE_STATIC}/external/
 mkdir ${QPACKAGE_STATIC}/external/ezfio
 cp --recursive -- ${QP_ROOT}/external/ezfio/Python ${QPACKAGE_STATIC}/external/ezfio/
 cp --recursive -- ${QP_ROOT}/external/ezfio/Bash   ${QPACKAGE_STATIC}/external/ezfio/
+cp --recursive -- ${QP_ROOT}/man   ${QPACKAGE_STATIC}/man
 
 
 #
 echo "Creating quantum_package.rc"
 #     ---------------------------
 
-sed "s!^export QP_ROOT=.*\$!export QP_ROOT=\$( cd \$(dirname "\${BASH_SOURCE}")  ; pwd -P )!" ${QP_ROOT}/quantum_package.rc.default > ${QPACKAGE_STATIC}/quantum_package.rc
+cp ${QP_ROOT}/quantum_package.rc ${QPACKAGE_STATIC}/
+cp ${QP_ROOT}/etc/* ${QPACKAGE_STATIC}/etc/
+cat << EOF > ${QPACKAGE_STATIC}/etc/00.qp_root.rc
+export QP_ROOT="\$( cd \$(dirname \\\${BASH_SOURCE}) ; cd .. ; pwd -P )"
+EOF
+
+
+
 
 #exit 0
 #
diff --git a/scripts/qp_plugins b/scripts/qp_plugins
index 0baf377f..54cd5eb2 100755
--- a/scripts/qp_plugins
+++ b/scripts/qp_plugins
@@ -127,7 +127,6 @@ def main(arguments):
     if arguments["create"]:                        
         m_instance = ModuleHandler([QP_SRC])       
                                                    
-        print arguments
         l_children = arguments["<needed_modules>"] 
                                                    
         name = arguments["<name>"][0]              
@@ -226,7 +225,6 @@ def main(arguments):
 
         for name in l_name:
 
-            print name
             if name in d_local:
                 print "{0} Is already installed".format(name)
 
@@ -285,7 +283,6 @@ def main(arguments):
         for module in set(l_name_to_remove):
 
             uninstall = os.path.join(QP_SRC,module,"uninstall")
-            print uninstall
             if os.path.isfile(uninstall):
                 subprocess.check_call([uninstall])
 
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 6aa3769b..db221f45 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -134,190 +134,6 @@ end
 
 
 
-subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-  use bitmasks
-  implicit none
-  
-  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
-  integer(bit_kind), intent(in) :: phasemask(N_int,2)
-  logical, intent(in) :: bannedOrb(mo_num)
-  double precision, intent(in) :: coefs(N_states)
-  double precision, intent(inout) :: vect(N_states, mo_num)
-  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
-  integer :: i, j, k, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti
-  double precision :: hij
-  double precision, external :: get_phase_bi, mo_two_e_integral
-  
-  integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
-  integer, parameter :: turn2(2) = (/2,1/) 
-  
-  if(h(0,sp) == 2) then
-    h1 = h(1, sp)
-    h2 = h(2, sp)
-    do i=1,3
-      puti = p(i, sp)
-      if(bannedOrb(puti)) cycle
-      p1 = p(turn3_2(1,i), sp)
-      p2 = p(turn3_2(2,i), sp)
-      hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2, p1, h1, h2)
-      hij = hij * get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2, N_int)
-      do k=1,N_states
-        vect(k,puti) = vect(k,puti) + hij * coefs(k)
-      enddo
-    end do
-  else if(h(0,sp) == 1) then
-    sfix = turn2(sp)
-    hfix = h(1,sfix)
-    pfix = p(1,sfix)
-    hmob = h(1,sp)
-    do j=1,2
-      puti = p(j, sp)
-      if(bannedOrb(puti)) cycle
-      pmob = p(turn2(j), sp)
-      hij = mo_two_e_integral(pmob, pfix, hmob, hfix)
-      hij = hij * get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix, N_int)
-      do k=1,N_states
-        vect(k,puti) = vect(k,puti) + hij * coefs(k)
-      enddo
-    end do
-  else
-    puti = p(1,sp)
-    if(.not. bannedOrb(puti)) then
-      sfix = turn2(sp)
-      p1 = p(1,sfix)
-      p2 = p(2,sfix)
-      h1 = h(1,sfix)
-      h2 = h(2,sfix)
-      hij = (mo_two_e_integral(p1,p2,h1,h2) - mo_two_e_integral(p2,p1,h1,h2))
-      hij = hij * get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2, N_int)
-      do k=1,N_states
-        vect(k,puti) = vect(k,puti) + hij * coefs(k)
-      enddo
-    end if
-  end if
-end 
-
-
-
-subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-  use bitmasks
-  implicit none
-  
-  integer(bit_kind), intent(in)  :: gen(N_int, 2), mask(N_int, 2)
-  integer(bit_kind), intent(in) :: phasemask(N_int,2)
-  logical, intent(in)            :: bannedOrb(mo_num)
-  double precision, intent(in)   :: coefs(N_states)
-  double precision, intent(inout) :: vect(N_states, mo_num)
-  integer, intent(in)            :: sp, h(0:2, 2), p(0:3, 2)
-  integer                        :: i, hole, p1, p2, sh, k
-  logical                        :: ok
-
-  logical, allocatable           :: lbanned(:)
-  integer(bit_kind)              :: det(N_int, 2)
-  double precision               :: hij
-  double precision, external     :: get_phase_bi, mo_two_e_integral
-  
-  allocate (lbanned(mo_num))
-  lbanned = bannedOrb
-  sh = 1
-  if(h(0,2) == 1) sh = 2
-  hole = h(1, sh)
-  lbanned(p(1,sp)) = .true.
-  if(p(0,sp) == 2) lbanned(p(2,sp)) = .true.
-  !print *, "SPm1", sp, sh
-  
-  p1 = p(1, sp)
-  
-  if(sp == sh) then
-    p2 = p(2, sp)
-    lbanned(p2) = .true.
-    
-
-    double precision :: hij_cache(mo_num,2)
-    call get_mo_two_e_integrals(hole,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map)
-    call get_mo_two_e_integrals(hole,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map)
-
-    do i=1,hole-1
-      if(lbanned(i)) cycle
-      hij = hij_cache(i,1)-hij_cache(i,2)
-      if (hij /= 0.d0) then
-        hij = hij * get_phase_bi(phasemask, sp, sp, i, p1, hole, p2, N_int)
-        do k=1,N_states
-          vect(k,i) = vect(k,i) + hij * coefs(k)
-        enddo
-      endif
-    end do
-    do i=hole+1,mo_num
-      if(lbanned(i)) cycle
-      hij = hij_cache(i,2)-hij_cache(i,1)
-      if (hij /= 0.d0) then
-        hij = hij * get_phase_bi(phasemask, sp, sp, hole, p1, i, p2, N_int)
-        do k=1,N_states
-          vect(k,i) = vect(k,i) + hij * coefs(k)
-        enddo
-      endif
-    end do
-
-    call apply_particle(mask, sp, p2, det, ok,  N_int)
-    call i_h_j(gen, det, N_int, hij)
-    do k=1,N_states
-      vect(k,p2) = vect(k,p2) + hij * coefs(k)
-    enddo
-  else
-    p2 = p(1, sh)
-    call get_mo_two_e_integrals(hole,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map)
-    do i=1,mo_num
-      if(lbanned(i)) cycle
-      hij = hij_cache(i,1)
-      if (hij /= 0.d0) then
-        hij = hij * get_phase_bi(phasemask, sp, sh, i, p1, hole, p2, N_int)
-        do k=1,N_states
-          vect(k,i) = vect(k,i) + hij * coefs(k)
-        enddo
-      endif
-    end do
-  end if
-  deallocate(lbanned)
-
-  call apply_particle(mask, sp, p1, det, ok,  N_int)
-  call i_h_j(gen, det, N_int, hij)
-  do k=1,N_states
-    vect(k,p1) = vect(k,p1) + hij * coefs(k)
-  enddo
-end 
-
-
-subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
-  use bitmasks
-  implicit none
-  
-  integer(bit_kind), intent(in)  :: gen(N_int, 2), mask(N_int, 2)
-  integer(bit_kind), intent(in) :: phasemask(N_int,2)
-  logical, intent(in)            :: bannedOrb(mo_num)
-  double precision, intent(in)   :: coefs(N_states)
-  double precision, intent(inout) :: vect(N_states, mo_num)
-  integer, intent(in)            :: sp, h(0:2, 2), p(0:3, 2)
-  integer                        :: i,k
-  logical                        :: ok
-
-  logical, allocatable           :: lbanned(:)
-  integer(bit_kind)              :: det(N_int, 2)
-  double precision               :: hij
-  
-  allocate(lbanned(mo_num))
-  lbanned = bannedOrb
-  lbanned(p(1,sp)) = .true.
-  do i=1,mo_num
-    if(lbanned(i)) cycle
-    call apply_particle(mask, sp, i, det, ok, N_int)
-    call i_h_j(gen, det, N_int, hij)
-    do k=1,N_states
-      vect(k,i) = vect(k,i) + hij * coefs(k)
-    enddo
-  end do
-  deallocate(lbanned)
-end 
-
 
 subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,variance,norm,buf,subset,csubset)
   use bitmasks
@@ -354,6 +170,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
   double precision :: rss
   double precision, external :: memory_of_double, memory_of_int
+
   rss = memory_of_int( (8*N_int+5)*N_det + N_det_alpha_unique + 4*N_int*N_det_selectors)
   rss += memory_of_double(mo_num*mo_num*(N_states+1))
   call check_mem(rss,irp_here)
@@ -656,7 +473,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
           maskInd = maskInd + 1
           if(mod(maskInd, csubset) == (subset-1)) then  
-            
+
             call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
             if(fullMatch) cycle
           
diff --git a/src/cis/20.cis.bats b/src/cis/20.cis.bats
index 0a2ce66d..918a41f4 100644
--- a/src/cis/20.cis.bats
+++ b/src/cis/20.cis.bats
@@ -1,20 +1,21 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 function run() {
   thresh=1.e-5
   test_exe cis || skip
-  qp_edit -c $1
-  ezfio set_file $1
-  ezfio set determinants n_states  3
-  ezfio set davidson threshold_davidson 1.e-12
-  ezfio set mo_two_e_ints io_mo_two_e_integrals Write
-  qp_set_frozen_core $1
-  qp_run cis $1
-  energy1="$(ezfio get cis energy | tr '[]' ' ' | cut -d ',' -f 1)"
-  energy2="$(ezfio get cis energy | tr '[]' ' ' | cut -d ',' -f 2)"
-  energy3="$(ezfio get cis energy | tr '[]' ' ' | cut -d ',' -f 3)"
+  qp set_file $1
+  qp edit --check
+  qp set determinants n_states  3
+  qp set davidson threshold_davidson 1.e-12
+  qp set mo_two_e_ints io_mo_two_e_integrals Write
+  qp set_frozen_core 
+  qp run cis 
+  energy1="$(qp get cis energy | tr '[]' ' ' | cut -d ',' -f 1)"
+  energy2="$(qp get cis energy | tr '[]' ' ' | cut -d ',' -f 2)"
+  energy3="$(qp get cis energy | tr '[]' ' ' | cut -d ',' -f 3)"
   eq $energy1 $2 $thresh
   eq $energy2 $3 $thresh
   eq $energy3 $4 $thresh
diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats
index fdc959b2..5a828064 100644
--- a/src/cisd/30.cisd.bats
+++ b/src/cisd/30.cisd.bats
@@ -1,55 +1,63 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 function run() {
   thresh=1.e-5
   test_exe cisd || skip
-  qp_edit -c $1
-  ezfio set_file $1
-  ezfio set determinants n_states  2
-  ezfio set davidson threshold_davidson 1.e-12
-  ezfio set davidson n_states_diag 24
-  qp_run cisd $1
-  energy1="$(ezfio get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
-  energy2="$(ezfio get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
-  eq $energy1 $2 $thresh
-  eq $energy2 $3 $thresh
+  qp edit --check
+  qp set determinants n_states  2
+  qp set davidson threshold_davidson 1.e-12
+  qp set davidson n_states_diag 24
+  qp run cisd 
+  energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
+  energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
+  eq $energy1 $1 $thresh
+  eq $energy2 $2 $thresh
 }
 
 
 @test "SiH2_3B1" { # 1.53842s
-  run sih2_3b1.ezfio -290.015949171697 -289.805036176618
+  qp set_file sih2_3b1.ezfio
+  run -290.015949171697 -289.805036176618
 }
 
 @test "HBO" { # 4.42968s
-  run  hbo.ezfio  -100.2019254455993 -99.79484127741013 
+  qp set_file hbo.ezfio
+  run -100.2019254455993 -99.79484127741013 
 }
 
 @test "HCO" { # 6.6077s
-  run hco.ezfio -113.288687359997  -113.122945162967
+  qp set_file hco.ezfio
+  run -113.288687359997  -113.122945162967
 }
 
 @test "H2O" { # 7.0651s
-  run  h2o.ezfio  -76.22975602077072 -75.80609108747208 
+  qp set_file h2o.ezfio
+  run -76.22975602077072 -75.80609108747208 
 }
 
 @test "H2S" { # 7.42152s
-  run h2s.ezfio -398.853701416768 -398.519020035337
+  qp set_file h2s.ezfio
+  run -398.853701416768 -398.519020035337
 }
 
 @test "N2H4" { # 15.8394s
-  qp_set_mo_class n2h4.ezfio -core "[1-2]" -act "[3-24]" -del "[25-48]"
-  run n2h4.ezfio -111.366247464687 -110.990795989548
+  qp set_file n2h4.ezfio
+  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
+  run -111.366247464687 -110.990795989548
 }
 
 @test "H2O2" { # 16.3164s
-  qp_set_mo_class h2o2.ezfio -core "[1-2]" -act "[3-24]" -del "[25-38]"
-  run h2o2.ezfio -151.003775695363 -150.650247854914
+  qp set_file h2o2.ezfio
+  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]"
+  run -151.003775695363 -150.650247854914
 }
 
 @test "OH" { # 18.2159s
-  run oh.ezfio -75.6087472926588 -75.5370393736601
+  qp set_file oh.ezfio
+  run -75.6087472926588 -75.5370393736601
 }
 
 
@@ -58,86 +66,101 @@ function run() {
 
 @test "CH4" { # 19.821s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class ch4.ezfio -core "[1]" -act "[2-30]" -del "[31-59]"
-  run ch4.ezfio -40.2403962667047 -39.8433221754964
+  qp set_file ch4.ezfio
+  qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
+  run -40.2403962667047 -39.8433221754964
 }
 
 @test "SiH3" { # 20.2202s
   [[ -n $TRAVIS ]] && skip
-  run sih3.ezfio -5.57096611856522  -5.30950347928823
+  qp set_file sih3.ezfio
+  run -5.57096611856522  -5.30950347928823
 }
 
 @test "NH3" { # 20.6771s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class nh3.ezfio -core "[1-4]" -act "[5-72]"
-  run nh3.ezfio -56.2447484835843 -55.9521689975716
+  qp set_file nh3.ezfio
+  qp set_mo_class --core="[1-4]" --act="[5-72]"
+  run -56.2447484835843 -55.9521689975716
 }
 
 @test "DHNO" { # 24.7077s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class dhno.ezfio -core "[1-7]" -act "[8-64]" 
-  run dhno.ezfio -130.458814562403 -130.356308303681
+  qp set_file dhno.ezfio
+  qp set_mo_class --core="[1-7]" --act="[8-64]" 
+  run -130.458814562403 -130.356308303681
 }
 
 @test "H3COH" { # 24.7248s
   [[ -n $TRAVIS ]] && skip
-  run h3coh.ezfio  -115.204958752377 -114.755913828245
+  qp set_file h3coh.ezfio
+  run -115.204958752377 -114.755913828245
 }
 
 @test "[Cu(NH3)4]2+" { # 29.9956s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class cu_nh3_4_2plus.ezfio -core "[1-24]" -act "[25-45]" -del "[46-87]"
-  run  cu_nh3_4_2plus.ezfio -1862.98659549315  -1862.68813764356
+  qp set_file cu_nh3_4_2plus.ezfio
+  qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
+  run -1862.98684406958  -1862.68818035746
 }
 
 @test "ClF" { # 30.3225s
   [[ -n $TRAVIS ]] && skip
-  run clf.ezfio  -559.162476603880  -558.792395927088
+  qp set_file clf.ezfio
+  run -559.162476603880  -558.792395927088
 }
 
 @test "C2H2" { # 35.3324s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class c2h2.ezfio -act "[1-30]" -del "[31-36]"
-  run c2h2.ezfio -12.3566731164213 -11.9495394759914 
+  qp set_file c2h2.ezfio
+  qp set_mo_class --act="[1-30]" --del="[31-36]"
+  run -12.3566731164213 -11.9495394759914 
 }
 
 @test "ClO" { # 37.6949s
   [[ -n $TRAVIS ]] && skip
-  run clo.ezfio -534.5404021326773 -534.3818725793897 
+  qp set_file clo.ezfio
+  run -534.5404021326773 -534.3818725793897 
 }
 
 @test "F2" { # 45.2078s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class f2.ezfio -core "[1,2]" -act "[3-30]" -del "[31-62]"
-  run f2.ezfio  -199.056829527539 -198.731828008346
+  qp set_file f2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-62]"
+  run -199.056829527539 -198.731828008346
 }
 
 @test "SO2" { # 47.6922s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class so2.ezfio -core "[1-8]" -act "[9-87]" 
-  run so2.ezfio -41.5746738710350 -41.3800467740750
+  qp set_file so2.ezfio
+  qp set_mo_class --core="[1-8]" --act="[9-87]" 
+  run -41.5746738710350 -41.3800467740750
 }
 
 @test "SO" { # 51.2476s
   [[ -n $TRAVIS ]] && skip
-  run so.ezfio -26.0131812819785 -25.7053111980226
+  qp set_file so.ezfio
+  run -26.0131812819785 -25.7053111980226
 }
 
 @test "CO2" { # 95.3736s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class co2.ezfio -core "[1,2]" -act "[3-30]" -del "[31-42]"
-  run co2.ezfio -187.959378390998  -187.432502050556
+  qp set_file co2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
+  run -187.959378390998  -187.432502050556
 }
 
 @test "N2" { # 133.1814
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class n2.ezfio -core "[1,2]" -act "[3-40]" -del "[41-60]"
-  run n2.ezfio -109.275693633982 -108.757794570948 
+  qp set_file n2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]"
+  run -109.275693633982 -108.757794570948 
 }
 
 @test "HCN" { # 133.8696s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class hcn.ezfio -core "[1,2]" -act "[3-40]" -del "[41-55]"
-  run hcn.ezfio -93.0776334511721 -92.6684633795506
+  qp set_file hcn.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
+  run -93.0776334511721 -92.6684633795506
 }
 
diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f
index 3638d5b5..ae808d12 100644
--- a/src/davidson/u0_h_u0.irp.f
+++ b/src/davidson/u0_h_u0.irp.f
@@ -138,10 +138,18 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   integer                        :: maxab, n_singles_a, n_singles_b, kcol_prev
   integer*8                      :: k8
   logical                        :: compute_singles
+  integer*8                      :: last_found, left, right, right_max
+  double precision               :: rss, mem, ratio
 
-  !TODO 
-  compute_singles = .True.
-
+!  call resident_memory(rss)
+!  mem = dble(singles_beta_csc_size) / 1024.d0**3
+!
+!  compute_singles = (mem+rss > qp_max_mem)
+!
+!  if (.not.compute_singles) then
+!    provide singles_beta_csc
+!  endif
+compute_singles=.True.
 
   maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
   allocate(idx0(maxab))
@@ -154,7 +162,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   ! -------------------------------------------------
 
   PROVIDE N_int nthreads_davidson
-  !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)        &
+  !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads_davidson)        &
       !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,                &
       !$OMP          psi_bilinear_matrix_columns,                    &
       !$OMP          psi_det_alpha_unique, psi_det_beta_unique,      &
@@ -166,13 +174,15 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
       !$OMP          psi_bilinear_matrix_columns_loc,                &
       !$OMP          psi_bilinear_matrix_transp_rows_loc,            &
       !$OMP          istart, iend, istep, irp_here, v_t, s_t,        &
-      !$OMP          ishift, idx0, u_t, maxab, compute_singles)      &
+      !$OMP          ishift, idx0, u_t, maxab, compute_singles,      &
+      !$OMP          singles_alpha_csc,singles_alpha_csc_idx,        &
+      !$OMP          singles_beta_csc,singles_beta_csc_idx)          &        
       !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
       !$OMP          lcol, lrow, l_a, l_b,                           &
       !$OMP          buffer, doubles, n_doubles,                     &
       !$OMP          tmp_det2, hij, sij, idx, l, kcol_prev,          &
-      !$OMP          singles_a, n_singles_a, singles_b,              &
-      !$OMP          n_singles_b, k8)
+      !$OMP          singles_a, n_singles_a, singles_b, ratio,       &
+      !$OMP          n_singles_b, k8, last_found,left,right,right_max)
   
   ! Alpha/Beta double excitations
   ! =============================
@@ -201,12 +211,19 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
 
     if (kcol /= kcol_prev) then
+      tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
       if (compute_singles) then
-        tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
         call get_all_spin_singles_$N_int(                              &
             psi_det_beta_unique, idx0,                                 &
             tmp_det(1,2), N_det_beta_unique,                           &
             singles_b, n_singles_b)
+      else
+        n_singles_b = 0
+        !DIR$ LOOP COUNT avg(1000)
+        do k8=singles_beta_csc_idx(kcol),singles_beta_csc_idx(kcol+1)-1
+          n_singles_b = n_singles_b+1
+          singles_b(n_singles_b) = singles_beta_csc(k8)
+        enddo
       endif
     endif
     kcol_prev = kcol
@@ -214,15 +231,20 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! Loop over singly excited beta columns
     ! -------------------------------------
 
+    !DIR$ LOOP COUNT avg(1000)
     do i=1,n_singles_b
       lcol = singles_b(i)
 
       tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
 
-      l_a = psi_bilinear_matrix_columns_loc(lcol)
-      ASSERT (l_a <= N_det)
+!---
+!      if (compute_singles) then
 
-      if (compute_singles) then
+        l_a = psi_bilinear_matrix_columns_loc(lcol)
+        ASSERT (l_a <= N_det)
+
+        !DIR$ UNROLL(8)
+        !DIR$ LOOP COUNT avg(50000)
         do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
           lrow = psi_bilinear_matrix_rows(l_a)
           ASSERT (lrow <= N_det_alpha_unique)
@@ -238,11 +260,75 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
         call get_all_spin_singles_$N_int(                              &
             buffer, idx, tmp_det(1,1), j,                              &
             singles_a, n_singles_a )
-      endif
+
+!-----
+!      else
+!
+! ! Search for singles
+!
+!call cpu_time(time0)
+! ! Right boundary
+!        l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1
+!        ASSERT (l_a <= N_det)
+!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+!          lrow = psi_bilinear_matrix_rows(l_a)
+!          ASSERT (lrow <= N_det_alpha_unique)
+!
+!          left = singles_alpha_csc_idx(krow)
+!          right_max = -1_8
+!          right = singles_alpha_csc_idx(krow+1)
+!          do while (right-left>0_8)
+!            k8 = shiftr(right+left,1)
+!            if (singles_alpha_csc(k8) > lrow) then
+!              right = k8 
+!            else if (singles_alpha_csc(k8) < lrow) then
+!              left = k8 + 1_8
+!            else
+!              right_max = k8+1_8
+!              exit
+!            endif
+!          enddo
+!          if (right_max > 0_8) exit
+!          l_a = l_a-1
+!        enddo
+!        if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow)
+!
+! ! Search
+!        n_singles_a = 0
+!        l_a = psi_bilinear_matrix_columns_loc(lcol)
+!        ASSERT (l_a <= N_det)
+!
+!        last_found = singles_alpha_csc_idx(krow)
+!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+!          lrow = psi_bilinear_matrix_rows(l_a)
+!          ASSERT (lrow <= N_det_alpha_unique)
+!
+!          left = last_found
+!          right = right_max
+!          do while (right-left>0_8)
+!            k8 = shiftr(right+left,1)
+!            if (singles_alpha_csc(k8) > lrow) then
+!              right = k8 
+!            else if (singles_alpha_csc(k8) < lrow) then
+!              left = k8 + 1_8
+!            else
+!              n_singles_a += 1
+!              singles_a(n_singles_a) = l_a
+!              last_found = k8+1_8
+!              exit
+!            endif
+!          enddo
+!          l_a = l_a+1
+!        enddo
+!        j = j-1
+!
+!      endif
+!-----
 
       ! Loop over alpha singles
       ! -----------------------
 
+      !DIR$ LOOP COUNT avg(1000)
       do k = 1,n_singles_a
         l_a = singles_a(k)
         ASSERT (l_a <= N_det)
@@ -253,9 +339,10 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
         tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
         call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij)
         call get_s2(tmp_det,tmp_det2,$N_int,sij)
+        !DIR$ LOOP COUNT AVG(4)
         do l=1,N_st
           v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
-          s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
+          s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a) 
         enddo
       enddo
 
@@ -295,6 +382,8 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! Loop inside the beta column to gather all the connected alphas
     lcol = psi_bilinear_matrix_columns(k_a)
     l_a = psi_bilinear_matrix_columns_loc(lcol)
+
+    !DIR$ LOOP COUNT avg(200000)
     do i=1,N_det_alpha_unique
       if (l_a > N_det) exit
       lcol = psi_bilinear_matrix_columns(l_a)
@@ -316,6 +405,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! ----------------------------------
 
     tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
+    !DIR$ LOOP COUNT avg(1000)
     do i=1,n_singles_a
       l_a = singles_a(i)
       ASSERT (l_a <= N_det)
@@ -326,6 +416,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
       tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
       call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 1, hij)
 
+      !DIR$ LOOP COUNT AVG(4)
       do l=1,N_st
         v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
         ! single => sij = 0 
@@ -336,6 +427,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! Compute Hij for all alpha doubles
     ! ----------------------------------
     
+    !DIR$ LOOP COUNT avg(50000)
     do i=1,n_doubles
       l_a = doubles(i)
       ASSERT (l_a <= N_det)
@@ -344,6 +436,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
       ASSERT (lrow <= N_det_alpha_unique)
 
       call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
+      !DIR$ LOOP COUNT AVG(4)
       do l=1,N_st
         v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
         ! same spin => sij = 0
@@ -375,6 +468,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! Loop inside the alpha row to gather all the connected betas
     lrow = psi_bilinear_matrix_transp_rows(k_b)
     l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
+    !DIR$ LOOP COUNT avg(200000)
     do i=1,N_det_beta_unique
       if (l_b > N_det) exit
       lrow = psi_bilinear_matrix_transp_rows(l_b)
@@ -396,6 +490,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! ----------------------------------
     
     tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
+    !DIR$ LOOP COUNT avg(1000)
     do i=1,n_singles_b
       l_b = singles_b(i)
       ASSERT (l_b <= N_det)
@@ -407,6 +502,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
       call i_H_j_mono_spin( tmp_det, tmp_det2, $N_int, 2, hij)
       l_a = psi_bilinear_matrix_transp_order(l_b)
       ASSERT (l_a <= N_det)
+      !DIR$ LOOP COUNT AVG(4)
       do l=1,N_st
         v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
         ! single => sij = 0 
@@ -416,6 +512,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
     ! Compute Hij for all beta doubles
     ! ----------------------------------
     
+    !DIR$ LOOP COUNT avg(50000)
     do i=1,n_doubles
       l_b = doubles(i)
       ASSERT (l_b <= N_det)
@@ -427,6 +524,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
       l_a = psi_bilinear_matrix_transp_order(l_b)
       ASSERT (l_a <= N_det)
 
+      !DIR$ LOOP COUNT AVG(4)
       do l=1,N_st
         v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
         ! same spin => sij = 0 
@@ -454,6 +552,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   
     hij = diag_H_mat_elem(tmp_det,$N_int) 
     sij = diag_S_mat_elem(tmp_det,$N_int)
+    !DIR$ LOOP COUNT AVG(4)
     do l=1,N_st
       v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,k_a)
       s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,k_a)
diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f
index d8f97154..0ea4dbe0 100644
--- a/src/determinants/spindeterminants.irp.f
+++ b/src/determinants/spindeterminants.irp.f
@@ -135,33 +135,6 @@ BEGIN_TEMPLATE
   
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer, n_singles_max_$alpha ]
-  implicit none
-  BEGIN_DOC
-  ! Maximum number of single excitations with the $\\$alpha$ spin
-  END_DOC
-  n_singles_max_$alpha = elec_$alpha_num * (mo_num - elec_$alpha_num)
-  
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, psi_det_$alpha_unique_singles, (n_singles_max_$alpha,N_det_$alpha_unique) ]
-  implicit none
-  BEGIN_DOC
-  ! For each unique spin-determinant, the list of indices of spin-determinant that
-  ! are singly excited.
-  END_DOC
-  integer                        :: i
-  integer, allocatable           :: idx0(:)
-  
-  allocate( idx0(N_det_$alpha_unique) )
-  
-  call get_all_spin_singles_$N_int(                                  &
-      psi_det_beta_unique, idx0,                                     &
-      tmp_det(1,2), N_det_beta_unique,                               &
-      singles_b, n_singles_b)
-endif
-
-END_PROVIDER
 
 SUBST [ alpha ]
 
@@ -176,7 +149,7 @@ END_TEMPLATE
 integer function get_index_in_psi_det_alpha_unique(key,Nint)
   use bitmasks
   BEGIN_DOC
-  ! Returns the index of the determinant in the ``psi_det_alpha_unique`` array
+  ! Returns the index of the determinant in the :c:data:`psi_det_alpha_unique` array
   END_DOC
   implicit none
   
@@ -257,7 +230,7 @@ end
 integer function get_index_in_psi_det_beta_unique(key,Nint)
   use bitmasks
   BEGIN_DOC
-  ! Returns the index of the determinant in the ``psi_det_beta_unique`` array
+  ! Returns the index of the determinant in the :c:data:`psi_det_beta_unique` array
   END_DOC
   implicit none
   
@@ -493,7 +466,7 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Order which allows to go from ``psi_bilinear_matrix`` to ``psi_det``
+  ! Order which allows to go from :c:data:`psi_bilinear_matrix` to :c:data:`psi_det`
   END_DOC
   integer                        :: k
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k)
@@ -516,7 +489,7 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1)
   !
   ! Rows are $\alpha$ determinants and columns are $\beta$.
   !
-  ! Order refers to ``psi_det``
+  ! Order refers to :c:data:`psi_det`
   END_DOC
   integer                        :: i,j,k, l
   
@@ -548,7 +521,7 @@ END_PROVIDER
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Transpose of ``psi_bilinear_matrix``
+  ! Transpose of :c:data:`psi_bilinear_matrix`
   !
   ! $D_\beta^\dagger.C^\dagger.D_\alpha$
   !
@@ -609,7 +582,7 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows_loc, (N_det_alpha_uniq
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Location of the columns in the ``psi_bilinear_matrix``
+  ! Location of the columns in the :c:data:`psi_bilinear_matrix`
   END_DOC
   integer                        :: i,j,k, l
   
@@ -635,8 +608,8 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_transp_reverse , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Order which allows to go from ``psi_bilinear_matrix_order_transp`` to
-  ! ``psi_bilinear_matrix``
+  ! Order which allows to go from :c:data:`psi_bilinear_matrix_order_transp` to
+  ! :c:data:`psi_bilinear_matrix`
   END_DOC
   integer                        :: k
   psi_bilinear_matrix_order_transp_reverse = -1
@@ -757,7 +730,7 @@ subroutine generate_all_alpha_beta_det_products
       !$OMP PRIVATE(i,j,k,l,tmp_det,iproc)
   !$ iproc = omp_get_thread_num()
   allocate (tmp_det(N_int,2,N_det_alpha_unique))
-  !$OMP DO SCHEDULE(static,1)
+  !$OMP DO SCHEDULE(static,8)
   do j=1,N_det_beta_unique
     l = 1
     do i=1,N_det_alpha_unique
@@ -789,9 +762,7 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf
   ! Returns the indices of all the single and double excitations in the list of
   ! unique $\alpha$ determinants.
   !
-  ! ..warning::
-  ! 
-  !     The buffer is transposed.
+  ! Warning: The buffer is transposed.
   !
   END_DOC
   integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
@@ -887,7 +858,8 @@ end
 subroutine copy_psi_bilinear_to_psi(psi, isize)
   implicit none
   BEGIN_DOC
-  ! Overwrites ``psi_det`` and ``psi_coef`` with the wave function in bilinear order
+  ! Overwrites :c:data:`psi_det` and :c:data:`psi_coef` with the wave function
+  ! in bilinear order
   END_DOC
   integer, intent(in)            :: isize
   integer(bit_kind), intent(out) :: psi(N_int,2,isize)
@@ -900,19 +872,14 @@ subroutine copy_psi_bilinear_to_psi(psi, isize)
   enddo
 end
 
-BEGIN_PROVIDER [ integer, singles_alpha_size ]
-  implicit none
-  BEGIN_DOC
-  ! Dimension of the ``singles_alpha`` array
-  END_DOC
-  singles_alpha_size = elec_alpha_num * (mo_num - elec_alpha_num)
-END_PROVIDER
 
  BEGIN_PROVIDER [ integer*8, singles_alpha_csc_idx, (N_det_alpha_unique+1) ]
 &BEGIN_PROVIDER [ integer*8, singles_alpha_csc_size ]
   implicit none
   BEGIN_DOC
-  ! Dimension of the ``singles_alpha`` array
+  ! singles_alpha_csc_size : Dimension of the :c:data:`singles_alpha_csc` array
+  !
+  ! singles_alpha_csc_idx  : Index where the single excitations of determinant i start
   END_DOC
   integer                        :: i,j
   integer, allocatable           :: idx0(:), s(:)
@@ -924,10 +891,10 @@ END_PROVIDER
   !$OMP PARALLEL DEFAULT(NONE)                                       &
       !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique,       &
       !$OMP          idx0, N_int, singles_alpha_csc,                 &
-      !$OMP          singles_alpha_size, singles_alpha_csc_idx)      &
+      !$OMP          elec_alpha_num, mo_num, singles_alpha_csc_idx)      &
       !$OMP   PRIVATE(i,s,j)
-  allocate (s(singles_alpha_size))
-  !$OMP DO SCHEDULE(static,1)
+  allocate (s(elec_alpha_num * (mo_num-elec_alpha_num) ))
+  !$OMP DO SCHEDULE(static,64)
   do i=1, N_det_alpha_unique
     call get_all_spin_singles(                                       &
         psi_det_alpha_unique, idx0, psi_det_alpha_unique(1,i), N_int,&
@@ -950,7 +917,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ]
   implicit none
   BEGIN_DOC
-  ! Dimension of the singles_alpha array
+  ! Indices of all single excitations
   END_DOC
   integer                        :: i, k
   integer, allocatable           :: idx0(:)
@@ -977,6 +944,78 @@ END_PROVIDER
 
 
 
+ BEGIN_PROVIDER [ integer*8, singles_beta_csc_idx, (N_det_beta_unique+1) ]
+&BEGIN_PROVIDER [ integer*8, singles_beta_csc_size ]
+  implicit none
+  BEGIN_DOC
+  ! singles_beta_csc_size : Dimension of the :c:data:`singles_beta_csc` array
+  !
+  ! singles_beta_csc_idx  : Index where the single excitations of determinant i start
+  END_DOC
+  integer                        :: i,j
+  integer, allocatable           :: idx0(:), s(:)
+  allocate (idx0(N_det_beta_unique))
+  do i=1, N_det_beta_unique
+    idx0(i) = i
+  enddo
+  
+  !$OMP PARALLEL DEFAULT(NONE)                                       &
+      !$OMP   SHARED(N_det_beta_unique, psi_det_beta_unique,       &
+      !$OMP          idx0, N_int, singles_beta_csc,                 &
+      !$OMP          elec_beta_num, mo_num, singles_beta_csc_idx)      &
+      !$OMP   PRIVATE(i,s,j)
+  allocate (s(elec_beta_num*(mo_num-elec_beta_num)))
+  !$OMP DO SCHEDULE(static,1)
+  do i=1, N_det_beta_unique
+    call get_all_spin_singles(                                       &
+        psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,&
+        N_det_beta_unique, s, j)
+    singles_beta_csc_idx(i+1) = int(j,8)
+  enddo
+  !$OMP END DO
+  deallocate(s)
+  !$OMP END PARALLEL
+  deallocate(idx0)
+  
+  singles_beta_csc_idx(1) = 1_8
+  do i=2, N_det_beta_unique+1
+    singles_beta_csc_idx(i) = singles_beta_csc_idx(i) + singles_beta_csc_idx(i-1)
+  enddo
+  singles_beta_csc_size = singles_beta_csc_idx(N_det_beta_unique+1)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, singles_beta_csc, (singles_beta_csc_size) ]
+  implicit none
+  BEGIN_DOC
+  ! Indices of all single excitations
+  END_DOC
+  integer                        :: i, k
+  integer, allocatable           :: idx0(:)
+  allocate (idx0(N_det_beta_unique))
+  do i=1, N_det_beta_unique
+    idx0(i) = i
+  enddo
+  
+  !$OMP PARALLEL DO DEFAULT(NONE)                                    &
+      !$OMP   SHARED(N_det_beta_unique, psi_det_beta_unique,       &
+      !$OMP          idx0, N_int, singles_beta_csc, singles_beta_csc_idx)&
+      !$OMP   PRIVATE(i,k) SCHEDULE(static,64)
+  do i=1, N_det_beta_unique
+    call get_all_spin_singles(                                       &
+        psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,&
+        N_det_beta_unique, singles_beta_csc(singles_beta_csc_idx(i)),&
+        k)
+  enddo
+  !$OMP END PARALLEL DO
+  deallocate(idx0)
+  
+END_PROVIDER
+
+
+
+
+
 subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
   use bitmasks
   implicit none
diff --git a/src/ezfio_files/00.create.bats b/src/ezfio_files/00.create.bats
index 3aa1a7b5..0034e981 100644
--- a/src/ezfio_files/00.create.bats
+++ b/src/ezfio_files/00.create.bats
@@ -1,6 +1,7 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 function run {
   local INPUT=$1
@@ -14,12 +15,12 @@ function run {
   fi
   cp ${QP_ROOT}/tests/input/$INPUT .
   rm -rf $EZ
-  qp_create_ezfio_from_xyz \
-     $INPUT -b "$BASIS" -m $MULT -c $CHARGE $PSEUDO -o $EZ
-  qp_edit -c $EZ
-  ezfio set_file $EZ
-  ezfio set scf_utils thresh_scf 1.e-12
-  ezfio set ao_two_e_ints io_ao_two_e_integrals "Write" 
+  qp create_ezfio_from_xyz \
+     $INPUT --basis="$BASIS" -m $MULT -c $CHARGE $PSEUDO -o $EZ
+  qp edit --check
+  qp set scf_utils thresh_scf 1.e-12
+  qp set ao_two_e_ints io_ao_two_e_integrals "Write" 
+  qp set mo_two_e_ints io_mo_two_e_integrals "Write" 
 }
 
 
diff --git a/src/ezfio_files/01.convert.bats b/src/ezfio_files/01.convert.bats
index 6df43d1e..a977a67f 100644
--- a/src/ezfio_files/01.convert.bats
+++ b/src/ezfio_files/01.convert.bats
@@ -1,16 +1,16 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 function run {
   local INPUT=$1
   local EZ=$2
   cp ${QP_ROOT}/tests/input/$INPUT .
-  qp_convert_output_to_ezfio $INPUT -o $EZ
-  qp_edit -c $EZ
-  ezfio set_file $EZ
-  ezfio set scf_utils thresh_scf 1.e-12
-  echo "Write" > ${EZ}/ao_two_e_ints/io_ao_two_e_integrals
+  qp convert_output_to_ezfio $INPUT -o $EZ
+  qp set_file $EZ
+  qp edit --check
+  qp set scf_utils thresh_scf 1.e-12
 }
 
 @test "HBO GAMESS" {
@@ -23,5 +23,5 @@ function run {
 
 @test "[Cu(NH3)4]2+ GAMESS" {
   run cu_nh3_4_2plus.gms.out  cu_nh3_4_2plus.ezfio
-  ezfio set scf_utils thresh_scf 1.e-10
+  qp set scf_utils thresh_scf 1.e-10
 }
diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats
index 243c542b..d9b9198c 100644
--- a/src/fci/40.fci.bats
+++ b/src/fci/40.fci.bats
@@ -1,140 +1,163 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 function run() {
-  thresh=$3
+  thresh=$2
   test_exe fci || skip
-  qp_edit -c $1
-  ezfio set_file $1
-  ezfio set determinants n_det_max 8000
-  ezfio set determinants n_states  1
-  ezfio set davidson threshold_davidson 1.e-10
-  ezfio set davidson n_states_diag 8
-  qp_run fci $1
+  qp edit --check
+  qp set determinants n_det_max 8000
+  qp set determinants n_states  1
+  qp set davidson threshold_davidson 1.e-10
+  qp set davidson n_states_diag 8
+  qp run fci 
   energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
-  eq $energy1 $2 $thresh
+  eq $energy1 $1 $thresh
 }
 
 
 
 
 @test "NH3" { # 10.6657s
-  qp_set_mo_class nh3.ezfio -core "[1-4]" -act "[5-72]"
-  run nh3.ezfio -56.2447484821590 1.e-5
+  qp set_file nh3.ezfio
+  qp set_mo_class --core="[1-4]" --act="[5-72]"
+  run -56.2447484821590 1.e-5
 }
 
 @test "DHNO" { # 11.4721s
-  qp_set_mo_class dhno.ezfio -core "[1-7]" -act "[8-64]" 
-  run dhno.ezfio -130.45902272485 1.e-5
+  qp set_file dhno.ezfio
+  qp set_mo_class --core="[1-7]" --act="[8-64]" 
+  run -130.45902272485 1.e-5
 }
 
 @test "HCO" { # 12.2868s
-  run hco.ezfio -113.297580169167 1.444e-05 
+  qp set_file hco.ezfio
+  run -113.297580169167 1.444e-05 
 }
 
 @test "H2O2" { # 12.9214s
-  qp_set_mo_class h2o2.ezfio -core "[1-2]" -act "[3-24]" -del "[25-38]"
-  run h2o2.ezfio -151.004593814816 0.00011948 
+  qp set_file h2o2.ezfio
+  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]"
+  run -151.004593814816 0.00011948 
 }
 
 @test "HBO" { # 13.3144s
-  run  hbo.ezfio  -100.213113590746 1.36e-05
+  qp set_file hbo.ezfio
+  run -100.213113590746 1.36e-05
 }
 
 @test "H2O" { # 11.3727s
-  run  h2o.ezfio  -76.2358876720796  1.988e-05 
+  qp set_file h2o.ezfio
+  run -76.2358876720796  1.988e-05 
 }
 
 @test "ClO" { # 13.3755s
-  run clo.ezfio -534.545851735243 0.00019344
+  qp set_file clo.ezfio
+  run -534.545851735243 0.00019344
 }
 
 @test "SO" { # 13.4952s
-  run so.ezfio -26.0118045926651 0.00014494
+  qp set_file so.ezfio
+  run -26.0118045926651 0.00014494
 }
 
 @test "H2S" { # 13.6745s
   [[ -n $TRAVIS ]] && skip
-  run h2s.ezfio -398.859198067009 8.46e-06
+  qp set_file h2s.ezfio
+  run -398.859198067009 8.46e-06
 }
 
 @test "OH" { # 13.865s 
   [[ -n $TRAVIS ]] && skip
-  run oh.ezfio -75.6120973654659 1.744e-05 
+  qp set_file oh.ezfio
+  run -75.6120973654659 1.744e-05 
 }
 
 @test "SiH2_3B1" { # 13.938ss
   [[ -n $TRAVIS ]] && skip
-  run sih2_3b1.ezfio -290.017547995946 1.e-5
+  qp set_file sih2_3b1.ezfio
+  run -290.017547995946 1.e-5
 }
 
 @test "H3COH" { # 14.7299s
   [[ -n $TRAVIS ]] && skip
-  run h3coh.ezfio  -115.200348295051 0.00018132
+  qp set_file h3coh.ezfio
+  run -115.200348295051 0.00018132
 }
 
 @test "SiH3" { # 15.99s
   [[ -n $TRAVIS ]] && skip
-  run sih3.ezfio -5.57259338826877 1.116e-05 
+  qp set_file sih3.ezfio
+  run -5.57259338826877 1.116e-05 
 }
 
 @test "CH4" { # 16.1612s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class ch4.ezfio -core "[1]" -act "[2-30]" -del "[31-59]"
-  run ch4.ezfio -40.2410273920655 3.02e-06
+  qp set_file ch4.ezfio
+  qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]"
+  run -40.2410273920655 3.02e-06
 }
 
 @test "ClF" { # 16.8864s
   [[ -n $TRAVIS ]] && skip
-  run clf.ezfio  -559.171627972338 0.00021062
+  qp set_file clf.ezfio
+  run -559.171627972338 0.00021062
 }
 
 @test "SO2" { # 17.5645s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class so2.ezfio -core "[1-8]" -act "[9-87]" 
-  run so2.ezfio -41.5746738710646 1.e-5
+  qp set_file so2.ezfio
+  qp set_mo_class --core="[1-8]" --act="[9-87]" 
+  run -41.5746738710646 1.e-5
 }
 
 @test "C2H2" { # 17.6827s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class c2h2.ezfio -act "[1-30]" -del "[31-36]"
-  run c2h2.ezfio -12.3681909988587 9.402e-05 
+  qp set_file c2h2.ezfio
+  qp set_mo_class --act="[1-30]" --del="[31-36]"
+  run -12.3681909988587 9.402e-05 
 }
 
 @test "N2" { # 18.0198s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class n2.ezfio -core "[1,2]" -act "[3-40]" -del "[41-60]"
-  run n2.ezfio -109.291407960731 0.00010052
+  qp set_file n2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]"
+  run -109.291407960731 0.00010052
 }
 
 @test "N2H4" { # 18.5006s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class n2h4.ezfio -core "[1-2]" -act "[3-24]" -del "[25-48]"
-  run n2h4.ezfio -111.367266319251 0.00010255
+  qp set_file n2h4.ezfio
+  qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]"
+  run -111.367266319251 0.00010255
 }
 
 @test "CO2" { # 21.1748s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class co2.ezfio -core "[1,2]" -act "[3-30]" -del "[31-42]"
-  run co2.ezfio -187.968251806361 0.00028902
+  qp set_file co2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]"
+  run -187.968251806361 0.00028902
 }
 
 @test "F2" { # 21.331s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class f2.ezfio -core "[1,2]" -act "[3-30]" -del "[31-62]"
-  run f2.ezfio  -199.068698950474 0.00014534
+  qp set_file f2.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-62]"
+  run -199.068698950474 0.00014534
 }
 
 @test "[Cu(NH3)4]2+" { # 25.0417s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class cu_nh3_4_2plus.ezfio -core "[1-24]" -act "[25-45]" -del "[46-87]"
-  run  cu_nh3_4_2plus.ezfio -1862.98632761077  5.e-07
+  qp set_file cu_nh3_4_2plus.ezfio
+  qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
+  run -1862.98632761077  5.e-07
 }
 
 @test "HCN" { # 20.3273s
   [[ -n $TRAVIS ]] && skip
-  qp_set_mo_class hcn.ezfio -core "[1,2]" -act "[3-40]" -del "[41-55]"
-  run hcn.ezfio -93.0774580352237 0.00016522
+  qp set_file hcn.ezfio
+  qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]"
+  run -93.0774580352237 0.00016522
 }
 
diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats
index eed81922..7d352041 100644
--- a/src/hartree_fock/10.hf.bats
+++ b/src/hartree_fock/10.hf.bats
@@ -1,15 +1,16 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 
 function run() {
   thresh=1.e-8
   test_exe scf || skip
-  qp_edit -c $1
-  ezfio set_file $1
-  qp_run scf $1
-  qp_set_frozen_core $1
+  qp set_file $1
+  qp edit --check
+  qp run scf 
+  qp set_frozen_core 
   energy="$(ezfio get hartree_fock energy)"
   eq $energy $2 $thresh
 }
@@ -20,7 +21,7 @@ function run() {
 }
 
 @test "SO" { # 0.539000
-  run so.ezfio -25.7175126082701
+  run so.ezfio -25.7175263371942
 }
 
 @test "HCO" { # 0.636700
@@ -102,8 +103,8 @@ function run() {
 
 @test "[Cu(NH3)4]2+" { # 59.610100
   [[ -n $TRAVIS ]] && skip
-  ezfio set_file cu_nh3_4_2plus.ezfio
-  ezfio set scf_utils thresh_scf 1.e-10
+  qp set_file cu_nh3_4_2plus.ezfio
+  qp set scf_utils thresh_scf 1.e-10
   run  cu_nh3_4_2plus.ezfio -1862.97590388214
 }
 
diff --git a/src/kohn_sham_rs/61.rsks.bats b/src/kohn_sham_rs/61.rsks.bats
index 6ee103a5..2705181f 100644
--- a/src/kohn_sham_rs/61.rsks.bats
+++ b/src/kohn_sham_rs/61.rsks.bats
@@ -1,19 +1,20 @@
 #!/usr/bin/env bats
 
 source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
 
 
 function run() {
   thresh=1.e-8
-  qp_edit -c $1
   functional=$2
-  ezfio set_file $1
-  ezfio set scf_utils thresh_scf 1.e-10
-  ezfio set dft_keywords exchange_functional $functional
-  ezfio set dft_keywords correlation_functional $functional
-  ezfio set ao_two_e_erf_ints mu_erf 0.5 
-  ezfio set becke_numerical_grid grid_type_sgn 1 
-  qp_run rs_ks_scf $1
+  qp set_file $1
+  qp edit --check
+  qp set scf_utils thresh_scf 1.e-10
+  qp set dft_keywords exchange_functional $functional
+  qp set dft_keywords correlation_functional $functional
+  qp set ao_two_e_erf_ints mu_erf 0.5 
+  qp set becke_numerical_grid grid_type_sgn 1 
+  qp run rs_ks_scf 
   energy="$(ezfio get kohn_sham_rs energy)"
   eq $energy $3 $thresh
 }
diff --git a/src/tools/print_e_conv.irp.f b/src/tools/print_e_conv.irp.f
index 8781ee22..d2ca4108 100644
--- a/src/tools/print_e_conv.irp.f
+++ b/src/tools/print_e_conv.irp.f
@@ -7,11 +7,11 @@ program print_e_conv
 !
 ! * a file "EZFIO.istate.conv" containing the variational and var+PT2 energies as a function of N_det
 !
-! * for istate > 1, a file EZFIO.$istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
+! * for istate > 1, a file EZFIO.istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
  END_DOC
 
   provide ezfio_filename
-  call routine
+  call routine_e_conv
  end
 
 subroutine routine_e_conv
diff --git a/src/tools/save_one_e_dm.irp.f b/src/tools/save_one_e_dm.irp.f
index 8c6c4772..7bcc5773 100644
--- a/src/tools/save_one_e_dm.irp.f
+++ b/src/tools/save_one_e_dm.irp.f
@@ -11,7 +11,7 @@ program save_one_e_dm
  END_DOC
   read_wf = .True.
   touch read_wf
-  call routine
+  call routine_save_one_e_dm
 
 end
 
diff --git a/src/utils/map_module.f90 b/src/utils/map_module.f90
index 9600cddf..e7bd7fdd 100644
--- a/src/utils/map_module.f90
+++ b/src/utils/map_module.f90
@@ -549,13 +549,13 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
   double precision, pointer :: v(:)
   integer :: i
   
-!  call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
-  call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend)
-!  if (idx > 0) then
-!    value = v(idx)
-!  else
-!    value = 0._integral_kind
-!  endif
+  call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
+  if (idx > 0) then
+    value = map%value(idx)
+  else
+    value = 0._integral_kind
+  endif
+!  call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend)
 end
 
 
@@ -665,7 +665,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
     
     istep = shiftr(iend-ibegin,1)
     idx = ibegin + istep
-    do while (istep > 64)
+    do while (istep > 4)
       idx = ibegin + istep
       ! TODO : Cache misses 
       if (cache_key < X(idx)) then
@@ -703,17 +703,17 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
       endif
     enddo
     idx = ibegin
-    if (min(iend_in,sze) > ibegin+64) then
-      iend = ibegin+64
+    if (min(iend_in,sze) > ibegin+4) then
+      iend = ibegin+4
+      !DIR$ LOOP COUNT MAX(4)
       do while (cache_key > X(idx))
         idx = idx+1
       end do
     else
+      !DIR$ LOOP COUNT MAX(4)
       do while (cache_key > X(idx))
         idx = idx+1
-        if (idx /= iend) then
-          cycle
-        else
+        if (idx == iend) then
           exit
         endif
       end do
@@ -771,10 +771,11 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
   iend   = min(iend_in,sze)
   if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
     
-    istep = shiftr(iend-ibegin,1)
+    istep = shiftr(iend+ibegin,1)
     idx = ibegin + istep
-    do while (istep > 64)
+    do while (istep > 4)
       idx = ibegin + istep
+      ! TODO : Cache misses 
       if (cache_key < X(idx)) then
         iend = idx
         istep = shiftr(idx-ibegin,1)
@@ -813,20 +814,17 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
       endif
     enddo
     idx = ibegin
-    value = Y(idx)
-    if (min(iend_in,sze) > ibegin+64) then
-      iend = ibegin+64
+    if (min(iend_in,sze) > ibegin+4) then
+      iend = ibegin+4
+      !DIR$ LOOP COUNT MAX(4)
       do while (cache_key > X(idx))
         idx = idx+1
-        value = Y(idx)
       end do
     else
+      !DIR$ LOOP COUNT MAX(4)
       do while (cache_key > X(idx))
         idx = idx+1
-        value = Y(idx)
-        if (idx /= iend) then
-          cycle
-        else
+        if (idx == iend) then
           exit
         endif
       end do
@@ -834,6 +832,8 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
     if (cache_key /= X(idx)) then
       idx = 1-idx
       value = 0.d0
+    else
+      value = Y(idx)
     endif
     return