diff --git a/INSTALL.rst b/INSTALL.rst
index 59a680d7..7a665bf9 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -201,7 +201,7 @@ OCaml
 
   .. code:: bash
 
-      opam install cryptokit zmq core sexplib ppx_sexp_conv ppx_deriving
+      opam install ocamlbuild cryptokit zmq core sexplib ppx_sexp_conv ppx_deriving getopt
 
 
 EZFIO
diff --git a/TODO b/TODO
index cfb152b5..624a90cb 100644
--- a/TODO
+++ b/TODO
@@ -46,3 +46,12 @@ Refaire les benchmarks
 # Documentation de qpsh
 
 # Documentation de /etc
+
+# 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
diff --git a/configure b/configure
index 9aec9834..7a9cdad3 100755
--- a/configure
+++ b/configure
@@ -37,6 +37,7 @@ EOF
 }
 
 PACKAGES=""
+OCAML_PACKAGES="ocamlbuild cryptokit zmq core sexplib ppx_sexp_conv ppx_deriving getopt"
 
 while : ; do
     case "$1" in
@@ -233,8 +234,7 @@ EOF
                  --yes --comp=4.07.0
 
             eval $(${QP_ROOT}/bin/opam env)
-            opam install -y ocamlbuild cryptokit zmq core sexplib \
-                 ppx_sexp_conv ppx_deriving
+            opam install -y ${OCAML_PACKAGES}
 
             
     elif [[ ${PACKAGE} = ezfio ]] ; then
diff --git a/docs/source/Makefile b/docs/source/Makefile
index e11deede..bf51e44b 100644
--- a/docs/source/Makefile
+++ b/docs/source/Makefile
@@ -1,6 +1,9 @@
 default:
 	./auto_generate.py
 	make -C ../ html
+	rm -rf ../build/man/
+	make -C ../ man
+	cp -r ../build/man/* $(QP_ROOT)/man/
 clean:
 	make -C ../ clean
 	rm modules/*.rst
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 8aa624c3..edd568a6 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -157,10 +157,14 @@ latex_documents = [
 
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
-man_pages = [
-    (master_doc, 'quantumpackage', 'Quantum Package Documentation',
-     [author], 1)
-]
+man_pages = []
+qpdoc = ' | Quantum Package >'
+import os
+for f in os.listdir("users_guide"):
+  name = f.split('.')[0]
+  if name not in ["index","quickstart"]:
+    filename = os.path.join("users_guide",name)
+    man_pages.append( (filename, name, qpdoc, [author], 1) )
 
 
 # -- Options for Texinfo output ----------------------------------------------
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 0e068bd3..6c1b242d 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -37,6 +37,7 @@
 
    programmers_guide/programming
    programmers_guide/ezfio
+   /programmers_guide/plugins
    programmers_guide/index
 
 
diff --git a/docs/source/modules/cipsi.rst b/docs/source/modules/cipsi.rst
index 935a8933..bf084f35 100644
--- a/docs/source/modules/cipsi.rst
+++ b/docs/source/modules/cipsi.rst
@@ -755,6 +755,20 @@ Subroutines / functions
 
 
 
+.. c:function:: remove_duplicates_in_selection_buffer
+
+    .. code:: text
+
+        subroutine remove_duplicates_in_selection_buffer(b)
+
+    File: :file:`selection_buffer.irp.f`
+
+    
+
+
+
+
+
 .. c:function:: run_cipsi
 
     .. code:: text
diff --git a/docs/source/modules/davidson.rst b/docs/source/modules/davidson.rst
index dbc7f628..f2d9fbe7 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_468`
+    File: :file:`u0_h_u0.irp.f_template_477`
 
     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_468`
+    File: :file:`u0_h_u0.irp.f_template_477`
 
     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_468`
+    File: :file:`u0_h_u0.irp.f_template_477`
 
     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_468`
+    File: :file:`u0_h_u0.irp.f_template_477`
 
     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_468`
+    File: :file:`u0_h_u0.irp.f_template_477`
 
     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 c2d56042..715d065a 100644
--- a/docs/source/modules/determinants.rst
+++ b/docs/source/modules/determinants.rst
@@ -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_141`
+    File: :file:`spindeterminants.irp.f_template_143`
 
     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_141`
+    File: :file:`spindeterminants.irp.f_template_143`
 
     Unique  :math:`\beta`  determinants
 
@@ -1215,7 +1215,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_141`
+    File: :file:`spindeterminants.irp.f_template_143`
 
     Unique  :math:`\alpha`  determinants
 
@@ -1242,7 +1242,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_141`
+    File: :file:`spindeterminants.irp.f_template_143`
 
     Unique  :math:`\beta`  determinants
 
@@ -2337,7 +2337,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2355,7 +2355,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2373,7 +2373,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_4(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2391,7 +2391,7 @@ Subroutines / functions
 
         subroutine get_all_spin_doubles_N_int(buffer, idx, spindet, size_buffer, doubles, n_doubles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2445,7 +2445,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2463,7 +2463,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2481,7 +2481,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_4(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2539,7 +2539,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_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2559,7 +2559,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_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2579,7 +2579,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_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2599,7 +2599,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_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
@@ -2619,7 +2619,7 @@ Subroutines / functions
 
         subroutine get_all_spin_singles_N_int(buffer, idx, spindet, size_buffer, singles, n_singles)
 
-    File: :file:`spindeterminants.irp.f_template_1218`
+    File: :file:`spindeterminants.irp.f_template_1221`
 
     
 
diff --git a/docs/source/modules/dft_utils_one_e.rst b/docs/source/modules/dft_utils_one_e.rst
index 09708468..5ed535d1 100644
--- a/docs/source/modules/dft_utils_one_e.rst
+++ b/docs/source/modules/dft_utils_one_e.rst
@@ -1090,14 +1090,12 @@ Providers
 
     .. code:: text
 
-        double precision, allocatable	:: potential_sr_x_alpha_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_x_beta_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_c_alpha_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_c_beta_ao_lda	(ao_num,ao_num,N_states)
 
     File: :file:`sr_pot_ao.irp.f`
 
-    short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
+    short range correlation alpha/beta potentials with LDA functional on the |AO| basis
 
 
 
@@ -1122,14 +1120,12 @@ Providers
 
     .. code:: text
 
-        double precision, allocatable	:: potential_sr_x_alpha_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_x_beta_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_c_alpha_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_c_beta_ao_lda	(ao_num,ao_num,N_states)
 
     File: :file:`sr_pot_ao.irp.f`
 
-    short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
+    short range correlation alpha/beta potentials with LDA functional on the |AO| basis
 
 
 
@@ -1156,12 +1152,10 @@ Providers
 
         double precision, allocatable	:: potential_sr_x_alpha_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_x_beta_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_c_alpha_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_c_beta_ao_lda	(ao_num,ao_num,N_states)
 
     File: :file:`sr_pot_ao.irp.f`
 
-    short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
+    short range exchange alpha/beta potentials with LDA functional on the |AO| basis
 
 
 
@@ -1188,12 +1182,10 @@ Providers
 
         double precision, allocatable	:: potential_sr_x_alpha_ao_lda	(ao_num,ao_num,N_states)
         double precision, allocatable	:: potential_sr_x_beta_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_c_alpha_ao_lda	(ao_num,ao_num,N_states)
-        double precision, allocatable	:: potential_sr_c_beta_ao_lda	(ao_num,ao_num,N_states)
 
     File: :file:`sr_pot_ao.irp.f`
 
-    short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
+    short range exchange alpha/beta potentials with LDA functional on the |AO| basis
 
 
 
diff --git a/docs/source/modules/fci.rst b/docs/source/modules/fci.rst
index 35f79f49..d4dce0d8 100644
--- a/docs/source/modules/fci.rst
+++ b/docs/source/modules/fci.rst
@@ -36,6 +36,23 @@ EZFIO parameters
 
 
 
+Providers
+---------
+
+
+.. c:var:: nthreads_pt2
+
+    .. code:: text
+
+        integer	:: nthreads_pt2
+
+    File: :file:`environment.irp.f`
+
+    Number of threads for Davidson
+
+
+
+
 Subroutines / functions
 -----------------------
 
diff --git a/docs/source/modules/hartree_fock.rst b/docs/source/modules/hartree_fock.rst
index 2005e9af..7ca55dfe 100644
--- a/docs/source/modules/hartree_fock.rst
+++ b/docs/source/modules/hartree_fock.rst
@@ -63,6 +63,79 @@ Providers
 ---------
 
 
+.. c:var:: ao_two_e_integral_alpha
+
+    .. code:: text
+
+        double precision, allocatable	:: ao_two_e_integral_alpha	(ao_num,ao_num)
+        double precision, allocatable	:: ao_two_e_integral_beta	(ao_num,ao_num)
+
+    File: :file:`fock_matrix_hf.irp.f`
+
+    Alpha Fock matrix in AO basis set
+
+
+
+
+.. c:var:: ao_two_e_integral_beta
+
+    .. code:: text
+
+        double precision, allocatable	:: ao_two_e_integral_alpha	(ao_num,ao_num)
+        double precision, allocatable	:: ao_two_e_integral_beta	(ao_num,ao_num)
+
+    File: :file:`fock_matrix_hf.irp.f`
+
+    Alpha Fock matrix in AO basis set
+
+
+
+
+.. c:var:: extra_e_contrib_density
+
+    .. code:: text
+
+        double precision	:: extra_e_contrib_density
+
+    File: :file:`hf_energy.irp.f`
+
+    Extra contribution to the SCF energy coming from the density. 
+
+    For a Hartree-Fock calculation: extra_e_contrib_density = 0 
+
+    For a Kohn-Sham or Range-separated Kohn-Sham: the exchange/correlation - trace of the V_xc potential
+
+
+
+
+.. c:var:: fock_matrix_ao_alpha
+
+    .. code:: text
+
+        double precision, allocatable	:: fock_matrix_ao_alpha	(ao_num,ao_num)
+        double precision, allocatable	:: fock_matrix_ao_beta	(ao_num,ao_num)
+
+    File: :file:`fock_matrix_hf.irp.f`
+
+    Alpha Fock matrix in AO basis set
+
+
+
+
+.. c:var:: fock_matrix_ao_beta
+
+    .. code:: text
+
+        double precision, allocatable	:: fock_matrix_ao_alpha	(ao_num,ao_num)
+        double precision, allocatable	:: fock_matrix_ao_beta	(ao_num,ao_num)
+
+    File: :file:`fock_matrix_hf.irp.f`
+
+    Alpha Fock matrix in AO basis set
+
+
+
+
 .. c:var:: hf_energy
 
     .. code:: text
@@ -113,6 +186,34 @@ Subroutines / functions
 
 
 
+.. c:function:: create_guess
+
+    .. code:: text
+
+        subroutine create_guess
+
+    File: :file:`scf.irp.f`
+
+    Create a MO guess if no MOs are present in the EZFIO directory
+
+
+
+
+
+.. c:function:: run
+
+    .. code:: text
+
+        subroutine run
+
+    File: :file:`scf.irp.f`
+
+    Run SCF calculation
+
+
+
+
+
 .. c:function:: scf
 
     .. code:: text
diff --git a/docs/source/modules/iterations.rst b/docs/source/modules/iterations.rst
index 81b89008..ed35c921 100644
--- a/docs/source/modules/iterations.rst
+++ b/docs/source/modules/iterations.rst
@@ -91,7 +91,7 @@ Subroutines / functions
 
     .. code:: text
 
-        subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_)
+        subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_st)
 
     File: :file:`print_summary.irp.f`
 
diff --git a/docs/source/modules/kohn_sham.rst b/docs/source/modules/kohn_sham.rst
index 1dac5f0f..4710a3d8 100644
--- a/docs/source/modules/kohn_sham.rst
+++ b/docs/source/modules/kohn_sham.rst
@@ -52,178 +52,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:: ao_two_e_integral_alpha
-
-    .. code:: text
-
-        double precision, allocatable	:: ao_two_e_integral_alpha	(ao_num,ao_num)
-        double precision, allocatable	:: ao_two_e_integral_beta	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_ks.irp.f`
-
-    Alpha Fock matrix in ao basis set
-
-
-
-
-.. c:var:: ao_two_e_integral_beta
-
-    .. code:: text
-
-        double precision, allocatable	:: ao_two_e_integral_alpha	(ao_num,ao_num)
-        double precision, allocatable	:: ao_two_e_integral_beta	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_ks.irp.f`
-
-    Alpha Fock matrix in ao basis set
-
-
-
-
-.. 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:: extra_e_contrib_density
-
-    .. code:: text
-
-        double precision	:: extra_e_contrib_density
-
-    File: :file:`ks_enery.irp.f`
-
-    Extra contribution to the SCF energy coming from the density. 
-
-    For a Hartree-Fock calculation: extra_e_contrib_density = 0 
-
-    For a Kohn-Sham or Range-separated Kohn-Sham: the exchange/correlation - 1/2 trace of the V_xc potential
-
-
-
-
-.. 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_ks.irp.f`
-
-    Mono electronic an Coulomb matrix in ao basis set
-
-
-
-
-.. c:var:: fock_matrix_ao_alpha
-
-    .. code:: text
-
-        double precision, allocatable	:: fock_matrix_ao_alpha	(ao_num,ao_num)
-        double precision, allocatable	:: fock_matrix_ao_beta	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_ks.irp.f`
-
-    Alpha Fock matrix in ao basis set
-
-
-
-
-.. c:var:: fock_matrix_ao_beta
-
-    .. code:: text
-
-        double precision, allocatable	:: fock_matrix_ao_alpha	(ao_num,ao_num)
-        double precision, allocatable	:: fock_matrix_ao_beta	(ao_num,ao_num)
-
-    File: :file:`fock_matrix_ks.irp.f`
-
-    Alpha Fock 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_ks.irp.f`
-
-    Mono electronic an Coulomb matrix in ao basis set
-
-
-
-
-.. c:var:: fock_matrix_energy
-
-    .. code:: text
-
-        double precision	:: ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`ks_enery.irp.f`
-
-    Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
 .. c:var:: ks_energy
 
     .. code:: text
@@ -241,104 +69,11 @@ Providers
 
 
 
-.. c:var:: one_electron_energy
-
-    .. code:: text
-
-        double precision	:: ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`ks_enery.irp.f`
-
-    Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
-.. c:var:: trace_potential_xc
-
-    .. code:: text
-
-        double precision	:: ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`ks_enery.irp.f`
-
-    Kohn-Sham energy containing the nuclear repulsion energy, and the various components of this quantity.
-
-
-
-
-.. c:var:: two_electron_energy
-
-    .. code:: text
-
-        double precision	:: ks_energy
-        double precision	:: two_electron_energy
-        double precision	:: one_electron_energy
-        double precision	:: fock_matrix_energy
-        double precision	:: trace_potential_xc
-
-    File: :file:`ks_enery.irp.f`
-
-    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:`ks_scf.irp.f`
-
-    
-
-
-
-
-
-.. c:function:: create_guess
-
-    .. code:: text
-
-        subroutine create_guess
-
-    File: :file:`ks_scf.irp.f`
-
-    Create a MO guess if no MOs are present in the EZFIO directory
-
-
-
-
-
-.. c:function:: run
-
-    .. code:: text
-
-        subroutine run
-
-    File: :file:`ks_scf.irp.f`
-
-    Run SCF calculation
-
-
-
-
-
 .. c:function:: srs_ks_cf
 
     .. code:: text
diff --git a/docs/source/modules/kohn_sham_rs.rst b/docs/source/modules/kohn_sham_rs.rst
index 9ab0aa5e..ca51ab5b 100644
--- a/docs/source/modules/kohn_sham_rs.rst
+++ b/docs/source/modules/kohn_sham_rs.rst
@@ -60,6 +60,122 @@ 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
@@ -77,11 +193,59 @@ 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 bf647f5c..421fb197 100644
--- a/docs/source/modules/mo_basis.rst
+++ b/docs/source/modules/mo_basis.rst
@@ -343,6 +343,20 @@ 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:: save_mos
 
     .. code:: text
diff --git a/docs/source/modules/scf_utils.rst b/docs/source/modules/scf_utils.rst
index d9ec8f4f..8866bd75 100644
--- a/docs/source/modules/scf_utils.rst
+++ b/docs/source/modules/scf_utils.rst
@@ -97,9 +97,9 @@ EZFIO parameters
     Calculated HF energy
 
 
-.. option:: no_oa_or_av_opt
+.. option:: frozen_orb_scf
 
-    If true, leave the active orbitals untouched in the SCF procedure
+    If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class
 
     Default: False
 
diff --git a/docs/source/programmers_guide/index.rst b/docs/source/programmers_guide/index.rst
index b0049ea4..3b4b7ec8 100644
--- a/docs/source/programmers_guide/index.rst
+++ b/docs/source/programmers_guide/index.rst
@@ -10,7 +10,7 @@ Index of Modules
    :glob:
 
    /modules/*
-   /programmers_guide/plugins
+   /programmers_guide/qp_*
 
 
 .. Auto-generated file
diff --git a/docs/source/programmers_guide/index_providers.rst b/docs/source/programmers_guide/index_providers.rst
index 01b305aa..c95aa43a 100644
--- a/docs/source/programmers_guide/index_providers.rst
+++ b/docs/source/programmers_guide/index_providers.rst
@@ -231,6 +231,7 @@ Index of Providers
 * :c:data:`fock_wee_closed_shell`
 * :c:data:`fps_spf_matrix_ao`
 * :c:data:`fps_spf_matrix_mo`
+* :c:data:`frozen_orb_scf`
 * :c:data:`full_ijkl_bitmask`
 * :c:data:`full_ijkl_bitmask_4`
 * :c:data:`gauleg_t2`
@@ -438,6 +439,7 @@ Index of Providers
 * :c:data:`no_vvvv_integrals`
 * :c:data:`nproc`
 * :c:data:`nthreads_davidson`
+* :c:data:`nthreads_pt2`
 * :c:data:`nucl_aos`
 * :c:data:`nucl_aos_transposed`
 * :c:data:`nucl_charge`
@@ -1213,8 +1215,10 @@ Index of Subroutines/Functions
 * :c:func:`read_dets`
 * :c:func:`recentered_poly2`
 * :c:func:`remove_duplicates_in_psi_det`
+* :c:func:`remove_duplicates_in_selection_buffer`
 * :c:func:`remove_small_contributions`
 * :c:func:`reorder_active_orb`
+* :c:func:`reorder_core_orb`
 * :c:func:`repeat_all_e_corr`
 * :c:func:`reset_zmq_addresses`
 * :c:func:`resident_memory`
diff --git a/docs/source/users_guide/qp_name.rst b/docs/source/programmers_guide/qp_name.rst
similarity index 100%
rename from docs/source/users_guide/qp_name.rst
rename to docs/source/programmers_guide/qp_name.rst
diff --git a/docs/source/users_guide/qp_test.rst b/docs/source/programmers_guide/qp_test.rst
similarity index 100%
rename from docs/source/users_guide/qp_test.rst
rename to docs/source/programmers_guide/qp_test.rst
diff --git a/etc/qp.rc b/etc/qp.rc
index 524b4dca..8270a73e 100644
--- a/etc/qp.rc
+++ b/etc/qp.rc
@@ -29,11 +29,19 @@ Usage:
 "
 }
 
+#function test_ezfio()
+#{
+#      if [[ ! -d ${EZFIO_FILE} ]] ; then
+#        echo "qp: cannot access ${EZFIO_FILE}: No such file or directory"
+#        return 1
+#      fi
+#}
+
 function qp()
 {
   case $1 in
     "has"|"set"|"get"|"set_file"|"unset_file")
-      ezfio $@
+        ezfio $@
       ;;
 
     "set_frozen_core")
@@ -43,9 +51,9 @@ function qp()
 
     "create_ezfio_from_xyz")
       shift
-      [[ -n $EZFIO_FILE ]] && ezfio unset_file
       NAME=$(qp_create_ezfio_from_xyz $@)
       if [[ -d $NAME ]] ; then
+        [[ -d $EZFIO_FILE ]] && ezfio unset_file
         ezfio set_file $NAME
       else
         echo $NAME | more
@@ -77,6 +85,11 @@ function qp()
       qp_mpirun $@ ${EZFIO_FILE}
       ;;
 
+    "man")
+      shift
+      man $QP_ROOT/man/${1}.?
+      ;;
+
     *)
       _qp_usage
       ;;
@@ -125,6 +138,10 @@ _Complete()
         COMPREPLY=( $(compgen -W "-h -core -inact -act -virt -del" -- $cur ) )
         return 0
         ;;
+      man)
+        COMPREPLY=( $(compgen -W "$(cd ${QP_ROOT}/man ; \ls | sed '/\.[1-9] / /')" -- $cur ) )
+        return 0
+        ;;
       set|has|get)
         COMPREPLY=( $(compgen -W "$(cd ${EZFIO_FILE} ; \ls -d */ | sed 's|/||g')" -- $cur ) )
         return 0
@@ -141,6 +158,10 @@ _Complete()
   else
 
     case "${prev}" in
+      man)
+        COMPREPLY=( $(compgen -W "$(cd ${QP_ROOT}/man ; echo * | sed 's|\.[1-9] | |g')" -- $cur ) )
+        return 0
+        ;;
       set_file)
         COMPREPLY=( $(compgen -W "$(\ls -d */ | sed 's|/||g')" -- ${cur} ) )
         return 0
@@ -151,6 +172,7 @@ _Complete()
         ;;
       *)
         COMPREPLY=( $(compgen -W 'set_file \
+                                  man \
                                   create_ezfio_from_xyz \
                                   -h' -- $cur ) )
         return 0
diff --git a/man/.gitignore b/man/.gitignore
new file mode 100644
index 00000000..e69de29b
diff --git a/man/configure.1 b/man/configure.1
new file mode 100644
index 00000000..643db372
--- /dev/null
+++ b/man/configure.1
@@ -0,0 +1,76 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "CONFIGURE" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+configure \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+\&./configure [\-h | \-c <file> | \-i <package>]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-c <file>, \-\-config <file>
+Define a configuration file, in :file\(ga${QP_ROOT}/config/\(ga
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-h, \-\-help
+Print the help message
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-i <package>, \-\-install <package>
+Try to install <package>. Use at your own risk.
+.UNINDENT
+.SH EXAMPLE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+\&./configure
+\&./configure \-c config/gfortran.cfg
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/excited_states.1 b/man/excited_states.1
new file mode 100644
index 00000000..0b76385f
--- /dev/null
+++ b/man/excited_states.1
@@ -0,0 +1,63 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "EXCITED_STATES" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+excited_states \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+It is possible to run excited states calculations with the quantum package.  To
+do this, set \fBdeterminants n_states\fP to the number of requested states.
+The selection criterion will be the maximum of the selection criteria for each
+state. If the Davidson diagonalization has difficulties to converge, increase
+the \fBdavidson n_states_diag\fP value.
+.sp
+When computing multiple states, it is good to have the \fBdeterminants
+s2_eig\fP flag \fBtrue\fP\&. This will force the Davidson algorithm to choose only
+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
+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.
+.sp
+To include excited state of all possible symmetries, a simple trick is to
+run a preliminary multi\-state CIS calculation using the CIS program,
+and then running the selected FCI restarting from the CIS states, setting
+\fBdeterminants read_wf\fP  to \fBtrue\fP\&.
+.sp
+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.
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/interfaces.1 b/man/interfaces.1
new file mode 100644
index 00000000..91e362f8
--- /dev/null
+++ b/man/interfaces.1
@@ -0,0 +1,62 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "INTERFACES" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+interfaces \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+A few interfaces to external codes are available.
+.SH * -> QUANTUM PACKAGE
+.INDENT 0.0
+.TP
+.B \fI\%GAMESS\fP / Gaussian
+Using the \fI\%resultsFile\fP Python library, the geometry and MOs can be read.
+This is useful to make calculations with CAS \- SCF orbitals
+.UNINDENT
+.SH QUANTUM PACKAGE -> *
+.INDENT 0.0
+.TP
+.B \fI\%Molden\fP
+3D plots of Molecular Orbitals
+.TP
+.B FCIDUMP
+Interface with the FCI \- QMC program \fI\%NECI\fP, or the semi\-stochastic
+Heat\-Bath CI program \fI\%Dice\fP\&.
+.UNINDENT
+.sp
+\fI\%QMCPack\fP / \fI\%CHAMP\fP /
+\fI\%QMC=Chem\fP
+Trial wave functions can be used for QMC, with or without pseudo\-potentials.
+These interfaces are provided as \fI\%external plugins\fP\&.
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/natural_orbitals.1 b/man/natural_orbitals.1
new file mode 100644
index 00000000..488a671b
--- /dev/null
+++ b/man/natural_orbitals.1
@@ -0,0 +1,53 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "NATURAL_ORBITALS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+natural_orbitals \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+To produce state\-average natural orbitals, run
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_run save_natorb file.ezfio
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+The MOs will be replaced, so the two\-electron integrals and the wave function
+are invalidated as well.
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/plugins.1 b/man/plugins.1
new file mode 100644
index 00000000..cc3dad9f
--- /dev/null
+++ b/man/plugins.1
@@ -0,0 +1,101 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "PLUGINS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+plugins \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+\fIQuantum Package\fP has very few executables out of the box. Most of the time, external
+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,
+the remote repository needs to be downloaded, and the plugins of the
+repository can be selected for installation.
+.sp
+To download an external repository of plugins, run the following command:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins download http://somewhere/over/the/rainbow/ext_repo
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+This downloads a copy of the repository of external plugins \fBext_repo\fP
+in \fB$QP_ROOT/plugins\fP\&.
+.sp
+The list of available uninstalled plugins can be seen using:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins list \-u
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+Now, the specific plugin \fBext_module\fP contained in the repository
+\fBext_repo\fP can be installed using:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins install ext_module
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+The module is now accessible via a symbolic link in \fB$QP_ROOT/src\fP,
+and can be compiled as any module, running \fI\%Ninja\fP\&.
+.sp
+To remove the module, run
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins uninstall ext_module
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_convert_output_to_ezfio.1 b/man/qp_convert_output_to_ezfio.1
new file mode 100644
index 00000000..dbd78656
--- /dev/null
+++ b/man/qp_convert_output_to_ezfio.1
@@ -0,0 +1,88 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_CONVERT_OUTPUT_TO_EZFIO" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_convert_output_to_ezfio \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+This Python script uses the \fI\%resultsFile\fP Python library to gather the
+geometry, AOs and MOs from output files of \fI\%GAMESS\fP or Gaussian, and
+puts this data in an \fI\%EZFIO\fP database. Some constraints are necessary in the
+output file : the run needs to be a single point HF, DFT or CAS SCF\&.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_convert_output_to_ezfio <FILE.out> [\-o <EZFIO_DIRECTORY>]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-o <EZFIO_DIRECTORY>
+Renames the \fI\%EZFIO\fP directory. If this option is not present, the default
+name fill be \fB<FILE.out>.ezfio\fP
+.UNINDENT
+.sp
+\fBNOTE:\fP
+.INDENT 0.0
+.INDENT 3.5
+The following keywords are necessary for Gaussian
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+GFPRINT pop=Full
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.SH EXAMPLE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_convert_output_to_ezfio h2o.out \-o h2o
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_create_ezfio_from_xyz.1 b/man/qp_create_ezfio_from_xyz.1
new file mode 100644
index 00000000..2b7edc51
--- /dev/null
+++ b/man/qp_create_ezfio_from_xyz.1
@@ -0,0 +1,231 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_CREATE_EZFIO_FROM_XYZ" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_create_ezfio_from_xyz \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+This command creates an \fI\%EZFIO\fP directory from a standard \fIxyz\fP file or from a
+\fIz\-matrix\fP file in Gaussian format.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_create_ezfio_from_xyz [FLAGS] (<xyz_file> | <zmt_file>)
+Flags :
+   \-b [\-au] [\-c int] [\-cart] [\-d float]
+   [\-m int] [\-o file] [\-p string] [\-help]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-b <basis_name>
+Name of basis set.  The basis set is defined as a single string if all the
+atoms are taken from the same basis set, otherwise specific elements can be
+defined as follows:
+.INDENT 7.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+\-b "cc\-pcvdz | H:cc\-pvdz | C:6\-31g"
+\-b "cc\-pvtz | 1,H:sto\-3g | 3,H:6\-31g"
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+By default, the basis set is obtained from the local database of the \fIQuantum Package\fP\&.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-au
+If present, input geometry is in atomic units.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-c <float>
+Total charge of the molecule. Default is 0.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-cart
+Compute AOs in the Cartesian basis set (6d, 10f, ...)
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-d <float>
+Add dummy atoms. x * (covalent radii of the atoms)
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-m <int>
+Spin multiplicity (2S+1) of the molecule. Default is 1.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-o <EZFIO_DIRECTORY>
+Name of the created \fI\%EZFIO\fP directory.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-p <string>
+Name of the pseudo\-potential
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-help, \-?
+Print the help text and exit
+.UNINDENT
+.SH USING CUSTOM ATOMIC BASIS SETS
+.sp
+If a file with the same name as the basis set exists, this file will be read.
+For example, if the file containing the basis set is named \fBcustom.basis\fP,
+and the \fIxyz\fP geometry is in \fBmolecule.xyz\fP, the following should be used:
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_create_ezfio_from_xyz \-b custom.basis molecule.xyz
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+Basis set files should be given in \fI\%GAMESS\fP format, where the full names of the
+atoms are given, and the basis sets for each element are separated by a blank line.
+Here is an example
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+HYDROGEN
+S   3
+1     13.0100000              0.0196850
+2      1.9620000              0.1379770
+3      0.4446000              0.4781480
+S   1
+1      0.1220000              1.0000000
+P   1
+1      0.7270000              1.0000000
+
+BORON
+S   8
+1   4570.0000000              0.0006960
+2    685.9000000              0.0053530
+3    156.5000000              0.0271340
+4     44.4700000              0.1013800
+5     14.4800000              0.2720550
+6      5.1310000              0.4484030
+7      1.8980000              0.2901230
+8      0.3329000              0.0143220
+S   8
+1   4570.0000000             \-0.0001390
+2    685.9000000             \-0.0010970
+3    156.5000000             \-0.0054440
+4     44.4700000             \-0.0219160
+5     14.4800000             \-0.0597510
+6      5.1310000             \-0.1387320
+7      1.8980000             \-0.1314820
+8      0.3329000              0.5395260
+S   1
+1      0.1043000              1.0000000
+P   3
+1      6.0010000              0.0354810
+2      1.2410000              0.1980720
+3      0.3364000              0.5052300
+P   1
+1      0.0953800              1.0000000
+D   1
+1      0.3430000              1.0000000
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH USING CUSTOM PSEUDO-POTENTIALS
+.sp
+As for the basis set, if a file with the same name as the pseudo\-potential
+exists, this file will be read.
+For example, if the file containing the custom pseudo\-potential is named
+\fBcustom.pseudo\fP, the basis set is named \fBcustom.basis\fP, and the \fIxyz\fP
+geometry is in \fBmolecule.xyz\fP, the following command should be used
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_create_ezfio_from_xyz \-b custom.basis \-p custom.pseudo molecule.xyz
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+Pseudo\-potential files should be given in a format very close to \fI\%GAMESS\fP
+format. The first line should be formatted as \fB%s GEN %d %d\fP where the
+first string is the chemical symbol, the first integer is the number of
+core electrons to be removed and the second integer is LMAX+1 as in \fI\%GAMESS\fP
+format.
+The pseudo\-potential for each element are separated by a blank line.
+Here is an example
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+Ne GEN 2 1
+3
+8.00000000 1 10.74945199
+85.99561593 3 10.19801460
+\-56.79004456 2 10.18694048
+1
+55.11144535 2 12.85042963
+
+F GEN 2 1
+3
+7.00000000 1 11.39210685
+79.74474797 3 10.74911370
+\-49.45159098 2 10.45120693
+1
+50.25646328 2 11.30345826
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_edit.1 b/man/qp_edit.1
new file mode 100644
index 00000000..b1eeff37
--- /dev/null
+++ b/man/qp_edit.1
@@ -0,0 +1,103 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_EDIT" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_edit \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+This command reads the content of the \fI\%EZFIO\fP directory and creates a temporary
+file containing the data. The data is presented as a \fIReStructured Text\fP (rst)
+document, where each section corresponds to the corresponding \fIQuantum Package\fP module.
+The content of the file can be modified to change the input parameters. When
+the text editor is closed, the updated data is saved into the \fI\%EZFIO\fP directory.
+.sp
+\fBNOTE:\fP
+.INDENT 0.0
+.INDENT 3.5
+The text editor which will be opened is defined by the \fBEDITOR\fP
+environment variable.
+.UNINDENT
+.UNINDENT
+.sp
+\fBWARNING:\fP
+.INDENT 0.0
+.INDENT 3.5
+When the wave function is too large (more than 10 000 determinants), the
+determinants are not displayed.
+.UNINDENT
+.UNINDENT
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_edit [FLAGS] <EZFIO_DIRECTORY>
+Flags :
+  [\-c]          [\-ndet int]
+  [\-state int]  [\-help]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-c
+Checks the input data
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-ndet <int>
+Truncates the wavefunction to the target number of determinants
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-state <int>
+Pick the target state as a new wavefunction.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-help
+Print the help text and exits
+.UNINDENT
+.sp
+Here is a short list of important control parameters :
+.INDENT 0.0
+.TP
+.B read_wf
+If \fBfalse\fP, initialize the calculation with a single\-determinant wave
+function. If \fBtrue\fP, initialize the calculation with the wave function stored
+in the \fI\%EZFIO\fP directory.
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_export_as_tgz.1 b/man/qp_export_as_tgz.1
new file mode 100644
index 00000000..897bcae9
--- /dev/null
+++ b/man/qp_export_as_tgz.1
@@ -0,0 +1,64 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_EXPORT_AS_TGZ" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_export_as_tgz \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.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.
+.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.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_export_as_tgz
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.sp
+\fBNOTE:\fP
+.INDENT 0.0
+.INDENT 3.5
+There can be conflicts due to the version of Glibc. The machine on which \fIQuantum Package\fP is
+compiled should be the oldest one.
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_plugins.1 b/man/qp_plugins.1
new file mode 100644
index 00000000..7f7f5a4f
--- /dev/null
+++ b/man/qp_plugins.1
@@ -0,0 +1,109 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_PLUGINS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_plugins \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\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
+be downloaded, and the plugins in these repositories can be
+installed/uninstalled of created.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_plugins list [ \-i | \-u | \-q ]
+qp_plugins download <url>
+qp_plugins install <name>...
+qp_plugins uninstall <name>
+qp_plugins create \-n <name> [\-r <repository>] [<needed_modules>...]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B list
+List all the available plugins.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-i
+List all the \fIinstalled\fP plugins.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-u
+List all the \fIuninstalled\fP plugins.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-q
+List all the downloaded repositories.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B download <url>
+Download an external repository. The URL points to a tar.gz file or a
+git repository, for example:
+.INDENT 7.0
+.IP \(bu 2
+\fI\%http://example.com/site/example.tar.gz\fP
+.IP \(bu 2
+\fI\%git@gitlab.com\fP:user/example_repository
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B install <plugin_name>
+Install the plugin \fBplugin_name\fP\&.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B uninstall <plugin_name>
+Uninstall the plugin \fBplugin_name\fP\&.
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-n <plugin_name>
+Create a new plugin named \fBplugin_name\fP (in local repository by default).
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-r <repository>
+Specify in which repository the new plugin will be created.
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_run.1 b/man/qp_run.1
new file mode 100644
index 00000000..1a747f00
--- /dev/null
+++ b/man/qp_run.1
@@ -0,0 +1,73 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_RUN" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_run \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+Command used to run a calculation.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_run [\-slave] [\-help] <PROGRAM> <EZFIO_DIRECTORY>
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-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
+another running instance of the \fIQuantum Package\fP\&.
+.UNINDENT
+.SH EXAMPLE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_run FCI h2o.ezfio
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_set_frozen_core.1 b/man/qp_set_frozen_core.1
new file mode 100644
index 00000000..0d72b392
--- /dev/null
+++ b/man/qp_set_frozen_core.1
@@ -0,0 +1,64 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_SET_FROZEN_CORE" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_set_frozen_core \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+Automatically finds n, the number of core electrons. Calls
+qp_set_mo_class setting all MOs as \fBActive\fP, except the n/2
+first ones which are set as \fBCore\fP\&.  If pseudo\-potentials are used, all the
+MOs are set as \fBActive\fP\&.
+.sp
+For elements on the right of the periodic table, \fIqp_set_frozen_core\fP will
+work as expected. But for elements on the left, a small core will be chosen. For
+example, a Carbon atom will have 2 core electrons, but a Lithium atom will have
+zero.
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_set_frozen_core [\-q]  <EZFIO_DIRECTORY>
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-q
+Prints in the standard output the number of core electrons.
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/qp_set_mo_class.1 b/man/qp_set_mo_class.1
new file mode 100644
index 00000000..205e1656
--- /dev/null
+++ b/man/qp_set_mo_class.1
@@ -0,0 +1,118 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "QP_SET_MO_CLASS" "1" "Jan 11, 2019" "2.0" "Quantum Package"
+.SH NAME
+qp_set_mo_class \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.sp
+This command sets the orbital classes in an \fI\%EZFIO\fP directory.
+.INDENT 0.0
+.TP
+.B Core
+MOs which are always doubly occupied
+.TP
+.B Deleted
+MOs which are never occupied
+.TP
+.B Active
+MOs in which any number of holes/particles can be made
+.TP
+.B Inactive
+MOs in which only holes can be made
+.TP
+.B Virtual
+MOs in which only particles can be made
+.UNINDENT
+.sp
+To avoid errors, all the MOs should be given a class.
+The range of MOs are given like the ranges in \fI\%SLURM\fP commands. For example,
+\fB"[36\-53,72\-107,126\-131]"\fP\&.
+.sp
+\fBTIP:\fP
+.INDENT 0.0
+.INDENT 3.5
+To quickly setup a frozen core calculation, the script qp_set_frozen_core
+can be used.
+.UNINDENT
+.UNINDENT
+.SH USAGE
+.INDENT 0.0
+.INDENT 3.5
+.sp
+.nf
+.ft C
+qp_set_mo_class [FLAGS] <EZFIO_DIRECTORY>
+Flags:
+   [\-act range]    [\-core range]   [\-del range]
+   [\-inact range]  [\-q]            [\-virt range]
+   [\-help]
+.ft P
+.fi
+.UNINDENT
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-act <range>
+Range of active orbitals
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-core <range>
+Range of core orbitals
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-del <range>
+Range of deleted orbitals
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-inact <range>
+Range of inactive orbitals
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-q
+Query: print the current masks
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-virt <range>
+Range of virtual orbitals
+.UNINDENT
+.INDENT 0.0
+.TP
+.B \-help, \-?
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2018, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/ocaml/Command_line.ml b/ocaml/Command_line.ml
new file mode 100644
index 00000000..8e7b407c
--- /dev/null
+++ b/ocaml/Command_line.ml
@@ -0,0 +1,161 @@
+type argument = With_arg | Without_arg | With_opt_arg 
+
+let anon_args = ref []
+and header_doc = ref ""
+and footer_doc = ref ""
+and specs = ref []
+
+let set_header_doc s = header_doc := s
+let set_footer_doc s = footer_doc := s
+
+let dict = 
+  let d = Hashtbl.create 67 in
+  d
+
+let get x = 
+  try Some (Hashtbl.find dict x)
+  with Not_found -> None
+
+let get_bool x = 
+  Hashtbl.mem dict x
+
+let show_help () =
+  get_bool "help"
+
+let anonymous ?(optional=false) name doc =
+  ( ' ', name, doc, 
+    if optional then With_opt_arg else Without_arg
+  )
+
+let help () =
+
+  Format.printf "@[%s@]@.@." !header_doc;
+  let get_param_from_doc doc =
+      match Str.split (Str.regexp "[ \n\r\x0c\t]+") doc with
+      | param :: doc -> param, (String.concat " " doc)
+      | _ -> failwith "Bad format for documentation"
+  in
+
+  let anon = 
+    List.filter (fun (x,_,_,_) -> x = ' ') !specs
+    |> List.map (fun x ->
+      match x with 
+      | (_,name,doc,Without_arg) -> (name,doc,false)
+      | (_,name,doc,_)  -> (name,doc,true)
+      )
+  in
+
+  let options = 
+    List.filter (fun (x,_,_,_) -> x <> ' ') !specs
+    |> List.sort (fun (x,_,_,_) (y,_,_,_) -> Char.compare x y)
+    |> List.map (fun x ->
+      match x with
+      | (short,long,doc,With_arg) -> (* with arg *)
+            let param, doc = get_param_from_doc doc in
+            (Printf.sprintf "-%c %s" short param,
+             Printf.sprintf "--%s=%s" long  param,
+             doc) 
+      | (short,long,doc,Without_arg) ->  (* without arg *)
+            (Printf.sprintf "-%c" short,
+             Printf.sprintf "--%s" long,
+             doc) 
+      | (short,long,doc,With_opt_arg) -> (* with or without arg *)
+            let param, doc = get_param_from_doc doc in
+            (Printf.sprintf "-%c [%s]" short param,
+             Printf.sprintf "--%s[=%s]" long  param,
+             doc) 
+      ) 
+  in
+        
+  let max_short =
+    List.map (fun (x,_,_) -> String.length x) options
+    |> List.fold_left max 0
+  in
+
+  let max_long  =
+    List.map (fun (_,x,_) -> String.length x) options
+    |> List.fold_left max 0
+  in
+
+  let fmt_opt max_w o = 
+    let l = String.length o in
+    o^(String.make (max_w-l) ' ')
+  in
+
+
+  let output_option ?(fixed_width=false) (short, long, doc) = 
+      if fixed_width then
+        Format.printf "@[%s  %s@]" 
+          (fmt_opt max_short short) (fmt_opt max_long long)
+      else
+        Format.printf "@[%s|%s@]" short long
+  in
+
+  let output_anon ?(fixed_width=false) (name, doc, optional) = 
+      if optional then
+        Format.printf "@[[%s]@]" name
+      else
+        Format.printf "@[%s@]" name
+  in
+
+  Format.printf "@[<v>@[<v 2>Usage:@,@,@[<hov 4>@[%s@]" Sys.argv.(0);
+  List.iter (fun x -> 
+    Format.printf "@ @[[";
+    output_option ~fixed_width:false x;
+    Format.printf "]@]"
+    ) options;
+  Format.printf "@ @[[--]@]";
+  List.iter (fun x -> 
+    Format.printf "@ @[";
+    output_anon ~fixed_width:false x;
+    Format.printf "@]"
+    ) anon;
+  Format.printf "@]@]@,@,";
+    
+  Format.printf "@[<v>Arguments:@,";
+
+  Format.printf "@[<v 2>@," ;
+  List.iter (fun (name,doc,optional) ->
+    Format.printf "@[<h>";
+    output_anon ~fixed_width:true (name,doc,optional);
+    Format.printf "@   @[<v 0>%s@]@]@," doc
+    ) anon;
+  Format.printf "@]@;";
+
+  Format.printf "@[<v>Options:@,";
+
+  Format.printf "@[<v 2>@," ;
+  List.iter (fun (short,long,doc) ->
+    Format.printf "@[<h>";
+    output_option ~fixed_width:true (short,long,doc);
+    Format.printf "@   @[<v 0>%s@]@]@," doc
+    ) options;
+  Format.printf "@]@;";
+
+  Format.printf "@[%s@]@." !footer_doc
+
+let set_specs specs_in =
+  specs := ( 'h', "help", "Prints the help message", Without_arg) :: specs_in;
+
+  let specs = 
+    List.filter (fun (x,_,_,_) -> x != ' ') !specs
+    |> List.map (fun x ->
+     match x with
+     | (short, long, doc, With_arg) ->
+         (short, long, None, Some (fun x -> Hashtbl.replace dict long x) )
+     | (short, long, doc, Without_arg) ->
+         (short, long, Some (fun () -> Hashtbl.replace dict long ""), None)
+     | (short, long, doc, With_opt_arg) ->
+         (short, long, Some (fun () -> Hashtbl.replace dict long ""),
+         Some (fun x -> Hashtbl.replace dict long x) )
+    ) 
+  in
+
+  Getopt.parse_cmdline specs (fun x -> anon_args := !anon_args @ [x]);
+
+  if show_help () then
+      (help () ; exit 0)
+
+let anon_args () = !anon_args
+
+
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index 25b983de..80a81503 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@@ -80,7 +80,10 @@ end = struct
   ;;
 
   let write_n_det n =
-    Det_number.to_int n
+    let n_det_old =
+      Ezfio.get_determinants_n_det ()
+    in
+    min n_det_old (Det_number.to_int n)
     |> Ezfio.set_determinants_n_det
   ;;
 
diff --git a/ocaml/_tags b/ocaml/_tags
index bf267415..42843d25 100644
--- a/ocaml/_tags
+++ b/ocaml/_tags
@@ -1,3 +1,3 @@
-true: package(core,cryptokit,zmq,str,ppx_sexp_conv,ppx_deriving)
+true: package(core,cryptokit,zmq,str,ppx_sexp_conv,ppx_deriving,getopt)
 true: thread
 false: profile
diff --git a/ocaml/myocamlbuild.ml b/ocaml/myocamlbuild.ml
index 2980af57..25b40953 100644
--- a/ocaml/myocamlbuild.ml
+++ b/ocaml/myocamlbuild.ml
@@ -1,5 +1,4 @@
 open Ocamlbuild_plugin;;
-open Command;;
 
 dispatch begin function
   | Before_rules -> 
diff --git a/ocaml/qp_basis_clean.ml b/ocaml/qp_basis_clean.ml
deleted file mode 100644
index 2cc6d349..00000000
--- a/ocaml/qp_basis_clean.ml
+++ /dev/null
@@ -1,72 +0,0 @@
-open Core
-
-let filenames =
-  let dir_name = Qpackage.root^"/data/basis/"
-  in
-  Sys.readdir dir_name
-  |> Array.map ~f:(fun x -> dir_name^x)
-  |> Array.to_list
-;;
-
-let clean_file filename =
-  let command =
-    Printf.sprintf "cp -f %s %s.old" filename filename
-  in 
-  let () = 
-    match Sys.command command with
-    | 0 -> ()
-    | i -> failwith (Printf.sprintf "Command %s exited with code %d\n" command i)
-  in
-
-  let lines =
-    In_channel.with_file filename ~f:In_channel.input_lines
-  in
-
-  Out_channel.with_file filename ~f:(fun out_channel ->
-
-    let rec loop ~do_s = function
-    | [] -> ()
-    | line :: tail ->
-      begin
-        let buffer = String.strip line
-            |> String.split ~on:' ' 
-            |> List.filter ~f:(fun x -> x <> "")
-        in
-        let () = 
-          match buffer with
-          | []       -> Printf.fprintf out_channel "\n"
-          | [ atom ] -> Printf.fprintf out_channel "%s\n" atom
-          | [ i ; expo ; coef ]  ->
-              Printf.fprintf out_channel "%3s  %14s  %14s\n" i expo coef
-          | [ i ; expo ; coef ; coef2 ]  ->
-              if (do_s) then
-                Printf.fprintf out_channel "%3s  %14s  %14s\n" i expo coef
-              else
-                Printf.fprintf out_channel "%3s  %14s  %14s\n" i expo coef2
-          | [ sym ; n ] -> 
-              if (sym = "L") then
-                 let () = 
-                   Printf.fprintf out_channel "S %3s\n" n
-                 in
-                 let rec build_newlist accu = function
-                   | (0, _)  
-                   | (_,[])  -> List.rev ((Printf.sprintf "P %3s\n" n)::accu)
-                   | (i,head::tail) ->
-                      build_newlist (head::accu) ( i-1, tail )
-                 in
-                 let newlist = build_newlist [] ((Int.of_string n),tail)
-                 in
-                 loop ~do_s:true newlist
-              else
-                 Printf.fprintf out_channel "%s %3s\n" sym n
-          | _ -> ()
-        in
-        loop ~do_s:do_s tail
-      end
-    in loop ~do_s:false lines
-  )
-;;
-
-List.iter ~f:clean_file filenames
-;;
-
diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml
index 4873c95b..aba08688 100644
--- a/ocaml/qp_create_ezfio_from_xyz.ml
+++ b/ocaml/qp_create_ezfio_from_xyz.ml
@@ -2,27 +2,6 @@ open Qputils
 open Qptypes
 open Core
 
-let spec =
-  let open Command.Spec in
-  empty 
-  +> flag "o" (optional string)
-     ~doc:"file Name of the created EZFIO file."
-  +> flag "b" (required string)
-     ~doc:"string Name of basis set."
-  +> flag "au" no_arg
-     ~doc:"Input geometry is in atomic units."
-  +> flag "c" (optional_with_default 0 int)
-     ~doc:"int Total charge of the molecule. Default is 0."
-  +> flag "d" (optional_with_default 0. float)
-     ~doc:"float Add dummy atoms. x * (covalent radii of the atoms)"
-  +> flag "m" (optional_with_default 1 int)
-     ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
-  +> flag "p" (optional string)
-     ~doc:"string Name of the pseudopotential"
-  +> flag "cart" no_arg
-     ~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)"
-  +> anon ("(xyz_file|zmt_file)" %: file )
-
 type element =
 | Element of Element.t
 | Int_elem of (Nucl_number.t * Element.t)
@@ -660,16 +639,16 @@ let run ?o b au c d m p cart xyz_file =
         end; 
         raise ex;
       end
-  in print_endline ezfio_file
+  in
+  ignore @@ Sys.command ("qp_edit -c "^ezfio_file);
+  print_endline ezfio_file
  
 
 
-let command = 
-    Command.basic_spec
-    ~summary: "Quantum Package command"
-    ~readme:(fun () -> "
 
-=== Available basis sets ===
+let () =
+
+  "=== Available basis sets ===
 
 " ^ (list_basis ()) ^ "
 
@@ -686,15 +665,68 @@ defined as follows:
 If a file with the same name as the basis set exists, this file will be read.
 Otherwise, the basis set is obtained from the database.
 
-" )
-    spec
-    (fun o b au c d m p cart xyz_file () ->
-       run ?o b au c d m p cart xyz_file )
-
-
-let () =
-    Command.run command
-
-
+" |> Command_line.set_header_doc ;
+
+  [ ( 'o', "output", "file Name of the created EZFIO file.", Command_line.With_arg) ;
+    ( 'b', "basis", "string Name of basis set.", Command_line.With_arg) ;
+    ( 'a', "au", "Input geometry is in atomic units.", Command_line.Without_arg) ;
+    ( 'c', "charge", "int Total charge of the molecule. Default is 0.", Command_line.With_arg) ;
+    ( 'd', "dummy", "float Add dummy atoms. x * (covalent radii of the atoms)",  Command_line.With_arg);
+    ( 'm', "multiplicity", "int Spin multiplicity (2S+1) of the molecule. Default is 1.",  Command_line.With_arg);
+    ( 'p', "pseudo", "string Name of the pseudopotential.", Command_line.With_arg);
+    ( 'x', "cartesian", "Compute AOs in the Cartesian basis set (6d, 10f, ...)", Command_line.Without_arg);
+    Command_line.anonymous "(xyz_file|zmt_file)" "input file in xyz format or z-matrix."
+  ]
+  |> Command_line.set_specs ;
+  
+
+  (*  Handle options *)
+  let output =
+    Command_line.get "output" 
+  in
+
+  let basis =
+    match Command_line.get "basis" with
+    | None -> (Command_line.help () ; failwith "Error: [-b|--basis] option is missing.")
+    | Some x -> x
+  in
+
+  let au =
+    Command_line.get_bool "au"
+  in
+
+  let charge =
+    match Command_line.get "charge" with
+    | None -> 0
+    | Some x -> int_of_string x
+  in
+
+  let dummy =
+    match Command_line.get "dummy" with
+    | None -> 0.
+    | Some x -> float_of_string x
+  in
+
+  let multiplicity =
+    match Command_line.get "multiplicity" with
+    | None -> 1
+    | Some n -> int_of_string n
+  in
+
+  let pseudo =
+    Command_line.get "pseudo" 
+  in
+
+  let cart =
+    Command_line.get_bool "cartesian" 
+  in
+
+  let xyz_filename =
+    match Command_line.anon_args () with
+    | [x] -> x
+    | _ -> (Command_line.help () ; failwith "input file is missing")
+  in
+
+  run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename
 
 
diff --git a/ocaml/qp_create_guess.ml b/ocaml/qp_create_guess.ml
deleted file mode 100644
index 6e0d6cfe..00000000
--- a/ocaml/qp_create_guess.ml
+++ /dev/null
@@ -1,142 +0,0 @@
-open Qputils
-open Qptypes
-open Core
-
-let run ~multiplicity ezfio_file =
-  if (not (Sys.file_exists_exn ezfio_file)) then
-    failwith ("EZFIO directory "^ezfio_file^" not found");
-  Ezfio.set_file ezfio_file;
-  let d =
-    Input.Determinants_by_hand.read ()
-  in
-  let m =
-    Multiplicity.of_int multiplicity
-  in
-  let ne = 
-    Ezfio.get_electrons_elec_alpha_num () +
-    Ezfio.get_electrons_elec_beta_num () 
-    |> Elec_number.of_int
-  in
-  let alpha, beta = 
-    let (a,b) = 
-      Multiplicity.to_alpha_beta ne m 
-    in
-    (Elec_alpha_number.to_int a, Elec_beta_number.to_int b)
-  in
-  let n_open_shells =
-    alpha - beta
-  in
-  let mo_num = 
-    Ezfio.get_mo_basis_mo_num ()
-  in
-  let build_list_of_dets ne n_closed n_open =
-    let init = 
-      Array.create ~len:n_closed Bit.One 
-      |> Array.to_list
-    in
-    let rec set_electron accu = function
-    | 1 ->  [ Bit.One :: accu ]
-    | i -> 
-          assert (i>1);
-          let rest =
-            set_electron (Bit.Zero :: accu) (i-1)
-          in
-          (Bit.One::accu) :: rest
-    in
-    let rec extend accu = function 
-    | 0 -> List.rev accu
-    | i -> extend (Bit.Zero::accu) (i-1)
-    in
-    let rec set_n_electrons accu imax = function
-    | 0 -> []
-    | 1 -> set_electron accu imax
-    | i ->
-          assert (i>1);
-          let l =
-            set_electron accu (imax-1)
-          in
-          List.map ~f:(fun x -> set_n_electrons x (imax-1) (i-1)) l
-          |> List.concat
-    in
-    set_n_electrons init n_open ne
-          |> List.filter ~f:(fun x -> List.length x <= n_closed+n_open) 
-          |> List.map ~f:(fun x -> extend x (((mo_num-1)/64+1)*64 - List.length x))
-  in
-
-  let alpha_new = 
-    (Elec_number.to_int ne + 1)/2
-  and beta_new = 
-    Elec_number.to_int ne/2
-  in
-  let l_alpha = 
-    build_list_of_dets ((alpha-beta+1)/2) beta n_open_shells
-  in 
-  let l_beta = 
-    if alpha_new = beta_new then
-      l_alpha
-    else
-      build_list_of_dets ((alpha-beta)/2)beta n_open_shells 
-  in
-
-  let n_int = 
-    Bitlist.n_int_of_mo_num mo_num
-  in
-  let determinants = 
-    List.map l_alpha ~f:(fun x -> List.map l_beta ~f:(fun y -> (x,y) ))
-    |> List.concat
-    |> List.map ~f:(fun pair -> Determinant.of_bitlist_couple ~n_int
-      ~alpha:(Elec_alpha_number.of_int alpha_new)
-      ~beta:(Elec_beta_number.of_int beta_new) pair )
-  in
-
-  let c = 
-    Array.init (List.length determinants) (fun _ -> Det_coef.of_float ((Random.float 2.)-.1.))
-  in
- 
-  determinants 
-  |> List.map ~f:(fun x -> Determinant.to_string ~mo_num:(MO_number.of_int mo_num) x) 
-  |> List.iter ~f:(fun x -> Printf.printf "%s\n\n%!" x);
-
-  let l =
-    List.length determinants
-  in
-  if l > 0 then
-    begin
-      let d = 
-        let s = (Float.of_int (alpha - beta)) *. 0.5 in
-        let open Input.Determinants_by_hand in
-        { d with n_int ; 
-          n_det = Det_number.of_int ~min:1 ~max:l  l;
-          expected_s2 = Positive_float.of_float (s *. (s +. 1.)) ;
-          psi_coef = c;
-          psi_det = Array.of_list determinants;
-        }
-      in
-      Input.Determinants_by_hand.write d;
-      Ezfio.set_determinants_read_wf true
-    end
-  else
-    Ezfio.set_determinants_read_wf false
-  
-
-
-let spec =
-  let open Command.Spec in
-  empty
-  +> flag "m" (required int)
-     ~doc:"int Spin multiplicity"
-  +> anon ("ezfio_file" %: string)
-
-let () = 
-  Command.basic_spec
-  ~summary: "Quantum Package command"
-  ~readme:( fun () -> "
-Creates an open-shell multiplet initial guess\n\n" )
-  spec
-  (fun multiplicity ezfio_file () ->
-    run ~multiplicity ezfio_file
-  )
-  |> Command.run   ~version: Git.sha1   ~build_info: Git.message
-
-
-
diff --git a/ocaml/qp_export.ml b/ocaml/qp_export.ml
deleted file mode 100644
index d3d7153e..00000000
--- a/ocaml/qp_export.ml
+++ /dev/null
@@ -1,176 +0,0 @@
-open Core
-open Qptypes
-
-let get_data () =
-  let mos = 
-    Input.Mo_basis.read ()
-  in
-  let molecule =
-    let e = 
-      match Input.Electrons.read () with
-      | Some x -> x
-      | None -> failwith "No electrons"
-    in
-    { Molecule.nuclei =
-        begin
-          match Input.Nuclei_by_hand.read () with
-          | None -> failwith "No nuclei"
-          | Some l -> Input.Nuclei_by_hand.to_atom_list l
-        end ;
-      Molecule.elec_alpha = e.Input.Electrons.elec_alpha_num ;
-      Molecule.elec_beta  = e.Input.Electrons.elec_beta_num  ;
-    }
-  in
-  let charge = 
-    Molecule.get_charge molecule
-    |> Charge.to_int
-  and mult = 
-    Molecule.get_multiplicity molecule
-    |> Multiplicity.to_int
-  in 
-  let hf =
-    if mult = 0 then "RHF" else "ROHF"
-  in
-  let guess = 
-    if mos = None then `Huckel else `Read
-  in
-  (mos, molecule, charge, mult, hf, guess)
-
-
-
-
-let run_g09 () =
-  let (mos, molecule, charge, mult, hf, guess) =
-    get_data ()
-  in
-  let guess = 
-    match guess with
-    | `Read -> "cards"
-    | `Huckel -> "huckel"
-  in
-  let result = [
-    Printf.sprintf "# %s/Gen guess=%s" hf guess;
-    "" ;
-    Molecule.name molecule ;
-    "" ;
-    Printf.sprintf "%d %d" charge mult ;
-    Printf.sprintf "%s" (
-      match String.split ~on:'\n' @@ Molecule.to_xyz molecule  with
-      | _ :: _ :: rest -> String.concat ~sep:"\n" rest
-      | _ -> assert false
-    );
-    "";
-    begin
-      match Input.Ao_basis.read () with
-      | Some x -> Input.Ao_basis.to_basis x
-      | None -> failwith "No basis"
-    end 
-    |> Basis.to_string ~fmt:Gto.Gaussian ~ele_array:(List.map ~f:(fun x -> x.Atom.element)
-                                                       molecule.Molecule.nuclei |> Array.of_list) ;
-    "";
-    begin
-      match mos with
-      | None -> ""
-      | Some mos' -> 
-        begin
-          Printf.sprintf "(E20.12)\n%s\n    0" 
-            ( 
-              Array.map mos'.Input.Mo_basis.mo_coef ~f:(fun x -> 
-                  Array.map x ~f:(fun y -> Printf.sprintf "%20.12E" (MO_coef.to_float y))
-                  |> Array.to_list
-                  |> String.concat ~sep:"\n")
-              |> Array.mapi ~f:(fun i -> Printf.sprintf "%5d\n%s" (i+1))
-              |> Array.to_list
-              |> String.concat ~sep:"\n"
-            )
-        end
-    end;
-    "" ; "" ;
-    "" ; "" ;
-  ]
-  in
-  List.iter ~f:(fun x -> print_endline x) result
-
-
-
-
-let run_gamess () =
-  failwith "Not yet implemented"
-(* TODO
-  let (mos, molecule, charge, mult, hf, guess) =
-    get_data ()
-  in
-  let guess = 
-    match guess with
-    | `Read -> "MOREAD"
-    | `Huckel -> "HUCKEL"
-  in
-  and coord  =
-    Printf.sprintf "%s" (
-      match String.split ~on:'\n' @@ Molecule.to_xyz molecule  with
-      | _ :: _ :: rest -> String.concat ~sep:"\n" rest
-      | _ -> assert false
-    );
-  and basis  =
-    begin
-      match !read_data.basis with
-      | None |> failwith "No basis set defined in command line"
-      | Some b |> b
-    end
-  and typ    = !read_data.typ
-  and vecfile =
-    begin
-      match !read_data.filename with
-      | None |> ""
-      | Some filename |> filename
-    end
-  and nstate = !read_data.nstate
-  and guess  = !read_data.guess
-  in
-
-  let system =
-    Gamess.{ mult ; charge ; basis ; coord }
-  in
-  Gamess.create_input ~vecfile ~system ~guess ~nstate typ
-  |> print_endline 
-*)
-
-
-let spec =
-  let open Command.Spec in
-  empty
-  +> flag "-gamess" no_arg ~doc:"GAMESS(US) input"
-  +> flag "-g09" no_arg ~doc:"Gaussian 09 input"
-  +> anon ("ezfio_file" %: string)
-
-let command = 
-  Command.basic_spec
-  ~summary: "Quantum Package command"
-  ~readme:(fun () -> 
-      "
-Exports the computed data for other programs.
-      ")
-  spec
-  (fun gamess g09 ezfio_file () -> 
-    let count_flags l =
-        List.fold_left ~f:(fun a x -> 
-            if x then a+1 else a) ~init:0 l
-    in
-    let () =
-      match count_flags [ g09 ; gamess ] with
-      | 0  -> raise (Invalid_argument "No file format specified")
-      | 1  -> ()
-      | _  -> raise (Invalid_argument "Too many file formats specified")
-    in
-    Ezfio.set_file ezfio_file;
-    if g09 then
-      run_g09 ()
-    else if gamess then
-      run_gamess ()
-  )
-
-    
-let () = 
-  Command.run command;
-  exit 0
-   
diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml
deleted file mode 100644
index 97360270..00000000
--- a/ocaml/qp_find_pi_space.ml
+++ /dev/null
@@ -1,164 +0,0 @@
-open Core
-open Qputils
-open Qptypes
-
-let run ?(sym="None") ?(apply="no") fc ezfio_filename =
-  Ezfio.set_file ezfio_filename ;
-
-  let aos = 
-    match Input.Ao_basis.read () with
-    | Some aos -> aos
-    | None -> failwith "Unable to read AOs"
-  and mos = 
-    match Input.Mo_basis.read () with
-    | Some mos -> mos
-    | None -> failwith "Unable to read MOs"
-  in
-  let rec find_power symmetry accu = function
-    | -1 -> accu
-    | i  -> let new_accu = 
-              if aos.Input.Ao_basis.ao_power.(i) = symmetry then (i::accu) 
-              else accu
-           in
-           find_power symmetry new_accu (i-1)
-  and n =
-    (AO_number.to_int aos.Input.Ao_basis.ao_num) 
-  and m =
-    (MO_number.to_int mos.Input.Mo_basis.mo_num) 
-  and one  = Positive_int.of_int 1
-  and zero =  Positive_int.of_int 0
-  in
-
-  (* Indices of the px, py and pz AOs *)
-  let x_indices =
-    find_power Symmetry.Xyz.{x=one  ; y=zero ; z=zero} [] (n-1)
-  and y_indices =
-    find_power Symmetry.Xyz.{x=zero ; y=one  ; z=zero} [] (n-1)
-  and z_indices =
-    find_power Symmetry.Xyz.{x=zero ; y=zero ; z=one } [] (n-1)
-  in
-
-  (* Compute the relative weight of each MO on the px, py, pz spaces *)
-  let compute_weight mo_i list_aos =
-    let num =
-       List.fold_left ~f:(fun s i -> s +. (MO_coef.to_float @@ mos.Input.Mo_basis.mo_coef.(mo_i).(i)) ** 2.) ~init:0. list_aos
-    and denom =
-       Array.fold ~f:(fun s x -> s +. (MO_coef.to_float x) ** 2.) ~init:0. mos.Input.Mo_basis.mo_coef.(mo_i) 
-    in
-    num /. denom
-  in
-  let result = 
-    Array.init ~f:(fun mo_i -> 
-      (mo_i+1,
-      compute_weight mo_i x_indices,
-      compute_weight mo_i y_indices,
-      compute_weight mo_i z_indices) ) m
-    |> Array.to_list
-  in
-
-  let pi, sigma =
-    let rec aux test_xyz (accu_pi: int list) (accu_sigma: int list) = function
-    | [] -> (List.rev accu_pi, List.rev accu_sigma)
-    | (mo_i,x,y,z)::rest ->
-        if test_xyz (x,y,z) then
-           aux test_xyz (mo_i::accu_pi) accu_sigma rest
-         else
-           aux test_xyz accu_pi (mo_i::accu_sigma) rest
-    in
-    match sym with
-    | "x" | "X" -> aux (fun (x,y,z) -> (x>y && x>z)) [] [] result
-    | "y" | "Y" -> aux (fun (x,y,z) -> (y>x && y>z)) [] [] result
-    | "z" | "Z" -> aux (fun (x,y,z) -> (z>x && z>y)) [] [] result
-    | _         -> ([],[])
-  in
-
-  begin
-    match sym with
-    | "x" | "X"  | "y" | "Y"  | "z" | "Z" -> 
-        if apply = "cas" then
-          begin
-            let ne = Ezfio.get_electrons_elec_alpha_num () in
-            let command = 
-              "qp_set_mo_class -act \"[" ^
-              (String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^
-              "]\" -core \"[" ^
-              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i <= ne) sigma) ^
-              "]\" -del \"[" ^
-              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > ne) sigma) ^
-              "]\" " ^ ezfio_filename
-            in
-            print_endline command;
-            if Sys.command command <> 0 then
-              failwith "Command failed"
-          end
-        else if apply = "cassd" then
-          begin
-            let ne = Ezfio.get_electrons_elec_alpha_num () in
-            let nfc = 
-              if fc then
-                let ic =
-                  Unix.open_process_in ("qp_set_frozen_core.py -q "^ezfio_filename)
-                in
-                let result = int_of_string @@
-                match In_channel.input_line ic with 
-                | Some line -> line
-                | None -> failwith "Failed in qp_set_frozen_core"
-                in
-                In_channel.close ic; result
-              else 0
-            in
-            let command = 
-              "qp_set_mo_class" ^
-              ( if nfc > 0 then
-                Printf.sprintf " -core \"[1-%d]\"" nfc else "") ^
-              " -inact \"[" ^
-              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > nfc && mo_i <= ne) sigma) ^ "]\"" ^
-              " -act \"[" ^ (String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^ "]\"" ^
-              " -virt \"[" ^
-              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > ne) sigma) ^ "]\" " ^ ezfio_filename
-            in
-            print_endline command;
-            if Sys.command command <> 0 then
-              failwith "Command failed"
-          end
-        else
-          begin
-            Printf.printf "Pi: [";
-            List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi;
-            Printf.printf "\b]\n\nSigma: [";
-            List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma;
-            Printf.printf "\b]\n"
-          end
-    | _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result
-  end
-
-  
-
-
-
-let spec =
-  let open Command.Spec in
-  empty
-  +> flag "sym" (optional string)
-    ~doc:"{x,y,z} Axis perpendicular to the plane"
-  +> flag "apply" (optional string)
-    ~doc:"[cas|cassd] Set the pi space as the active space, all other MOs frozen. If [cassd] is chosen, set inactive and virtual mos."
-  +> flag "fc" no_arg
-    ~doc:"freeze core MOs with -apply cassd"
-  +> anon ("ezfio_filename" %: string)
-
-
-
-let command =
-    Command.basic_spec
-    ~summary: "Quantum Package command"
-    ~readme:(fun () ->
-     "Find all the pi molecular orbitals to create a pi space.
-     ")
-    spec
-    (fun sym apply fc ezfio_filename () -> run ?sym ?apply fc ezfio_filename)
-
-
-let () =
-    Command.run command
-
diff --git a/ocaml/qp_overlap_of_wf.ml b/ocaml/qp_overlap_of_wf.ml
deleted file mode 100644
index 3197bae5..00000000
--- a/ocaml/qp_overlap_of_wf.ml
+++ /dev/null
@@ -1,93 +0,0 @@
-(** 
- * Computes the overlap <Psi_0 | Psi_1> where both Psi_0 and Psi_1 are truncated in the set
- * of common determinants and normalized
- *)
-
-open Input_determinants_by_hand
-open Qptypes
-
-
-let () =
-  let ezfio, ezfio' =
-    try
-      Sys.argv.(1), Sys.argv.(2)
-    with Invalid_argument _ ->
-      raise (Invalid_argument (Printf.sprintf
-         "Syntax : %s EZFIO1 EZFIO2" Sys.argv.(0)))
-  in
-
-  let fetch_wf ~state filename =
-    (* State 0 is the ground state *)
-    Ezfio.set_file filename;
-    let mo_num =
-      Ezfio.get_mo_basis_mo_num ()
-      |> MO_number.of_int
-    in
-    let d =
-      Determinants_by_hand.read ()
-    in
-    let n_det =
-      Det_number.to_int d.Determinants_by_hand.n_det
-    in
-    let state_shift = 
-      state*n_det
-    in
-    let keys = 
-      Array.map (Determinant.to_string ~mo_num) 
-        d.Determinants_by_hand.psi_det
-    and values =
-      Array.map Det_coef.to_float
-        d.Determinants_by_hand.psi_coef
-    in
-    let hash = 
-      Hashtbl.create n_det
-    in
-    for i=0 to n_det-1
-    do
-      Hashtbl.add hash keys.(i) values.(state_shift+i);
-    done;
-    hash
-  in
-
-  let overlap wf wf' =
-    let result, norm, norm' = 
-      Hashtbl.fold (fun k c (accu,norm,norm') ->
-        let (c',c) =
-          try  (Hashtbl.find wf' k, c)
-          with Caml.Not_found -> (0.,0.)
-        in
-        (accu +. c *. c' , 
-        norm +. c *. c  , 
-        norm'+. c'*. c' ) 
-      ) wf (0.,0.,0.)
-    in 
-    result /. (sqrt (norm *. norm'))
-  in
-
-  let n_st1 =
-      Ezfio.set_file ezfio;
-      Ezfio.get_determinants_n_states ()
-  and n_st2 = 
-      Ezfio.set_file ezfio';
-      Ezfio.get_determinants_n_states ()
-  in
-  Array.init n_st2 (fun i -> i)
-  |> Array.iter (fun state_j -> 
-      Printf.printf "%d  " (state_j+1);
-      let wf'  = 
-           fetch_wf ~state:state_j ezfio'
-      in
-      Array.init n_st1 (fun i -> i)
-      |> Array.iter (fun state_i ->
-        let wf = 
-           fetch_wf ~state:state_i ezfio
-        in
-        let o = 
-          overlap wf wf'
-        in
-        Printf.printf "%f %!" (abs_float o)
-      );
-      Printf.printf "\n%!"
-  )
-
-
diff --git a/ocaml/qp_print.ml b/ocaml/qp_print.ml
deleted file mode 100644
index efbdc01e..00000000
--- a/ocaml/qp_print.ml
+++ /dev/null
@@ -1,183 +0,0 @@
-open Qputils;;
-open Qptypes;;
-open Core;;
-
-type input_action =
-  | Basis
-  | Nuclei
-  | Charge
-  | Multiplicity
-  | Electrons
-;;
-
-let create_i_action = function
-  | "basis" -> Basis
-  | "nucl" -> Nuclei
-  | "charge" -> Charge
-  | "mult" -> Multiplicity
-  | "elec" -> Electrons
-  | _ -> raise (Failure "Action should be:
-  * basis
-  * nucl
-  * charge
-  * mult
-  * elec
-")
-
-;;
-
-let spec =
-  let open Command.Spec in
-  empty 
-  +> flag "i"  (optional string)
-     ~doc:"Prints input data"
-  +> flag "o" (optional string)
-     ~doc:"Prints output data"
-  +> anon ("ezfio_file" %: string)
-;;
-
-let run_i ~action ezfio_filename =
-
-  let action = create_i_action action in
-
-  if (not (Sys.file_exists_exn ezfio_filename)) then
-    failwith (ezfio_filename^" does not exist");
-
-  Ezfio.set_file ezfio_filename;
-
-  let print_basis () =
-    Printf.printf "%s\n" (Ezfio.get_ao_basis_ao_basis ())
-  in
-
-
-  let compute_charge () =
-    let input = match Input.Electrons.read () with
-    | Some x -> x
-    | None   -> assert false
-    in
-    let nucl_charge = Ezfio.get_nuclei_nucl_charge ()
-      |> Ezfio.flattened_ezfio |> Array.map ~f:(Float.to_int) 
-    and n_alpha = input.Input.Electrons.elec_alpha_num
-      |> Elec_alpha_number.to_int
-    and n_beta  = input.Input.Electrons.elec_beta_num 
-      |> Elec_beta_number.to_int
-    in Array.fold ~init:(-n_alpha-n_beta) ~f:(fun x y -> x+y) nucl_charge
-    |> Charge.of_int
-  in
-
-  let compute_multiplicity () =
-    let input = match Input.Electrons.read () with
-    | Some x -> x
-    | None   -> assert false
-    in
-    let n_alpha = input.Input.Electrons.elec_alpha_num
-    and n_beta  = input.Input.Electrons.elec_beta_num 
-    in Multiplicity.of_alpha_beta n_alpha n_beta
-  in
-
-  let create_molecule () =
-    let nucl_num = Ezfio.get_nuclei_nucl_num () 
-    and nucl_charge = Ezfio.get_nuclei_nucl_charge ()
-       |> Ezfio.flattened_ezfio
-    and nucl_coord = Ezfio.get_nuclei_nucl_coord ()
-       |> Ezfio.flattened_ezfio
-    in
-    let nucl_label =
-      if (Ezfio.has_nuclei_nucl_label ()) then
-       Ezfio.get_nuclei_nucl_label () |> Ezfio.flattened_ezfio
-      else
-        Array.map ~f:(fun x-> x
-          |> Charge.of_float
-          |> Element.of_charge 
-          |> Element.to_string ) nucl_charge
-    in
-    let buffer = ref "" in
-    for i=0 to (nucl_num-1) do
-      buffer := !buffer ^ (Printf.sprintf "%s %f %f %f %f\n"
-        nucl_label.(i)
-        nucl_charge.(i)
-        nucl_coord.(i)
-        nucl_coord.(i+nucl_num)
-        nucl_coord.(i+nucl_num+nucl_num)
-      )
-    done ;
-    let charge = compute_charge () in
-    let mult = compute_multiplicity () in
-    Molecule.of_xyz_string ~charge:charge ~multiplicity:mult !buffer
-  in
-
-  let print_nuclei () =
-    let molecule = create_molecule () in
-    print_endline (Molecule.to_string molecule)
-
-  and print_charge () = 
-    let molecule = create_molecule () in
-    Printf.printf "%s" (Charge.to_string (Molecule.get_charge molecule))
-
-  and print_multiplicity () = 
-    let molecule = create_molecule () in
-    Printf.printf "%s" (Multiplicity.to_string (Molecule.get_multiplicity
-    molecule))
-
-  and print_electrons () = ()
-
-  in
-  match action with
-  | Basis        -> print_basis ()
-  | Nuclei       -> print_nuclei ()
-  | Charge       -> print_charge ()
-  | Multiplicity -> print_multiplicity ()
-  | Electrons    -> print_electrons ()
-;;
-
-let run_o ~action ezfio_filename =
-
-  if (not (Sys.file_exists_exn ezfio_filename)) then
-    failwith (ezfio_filename^" does not exist");
-
-  (* Open EZFIO *)
-  Ezfio.set_file ezfio_filename;
-
-;;
-
-let command = 
-    Command.basic_spec
-    ~summary: "Quantum Package command"
-    ~readme:(fun () ->
-      "
-Prints data contained into the EZFIO file.
-
-Input data :
-
-  * basis
-  * nucl
-  * charge
-  * mult
-  * elec
-
-Output data :
-
-  * 
-      ")
-    spec
-    (fun i o ezfio_file () ->
-  try 
-      match (i,o) with
-      | (Some i, None) -> run_i ~action:i ezfio_file
-      | (None, Some o) -> run_o ~action:o ezfio_file
-      | (Some _, Some _) -> 
-          raise (Failure "Error : please specify -i or -o but not both.")
-      | (None, None) -> 
-          raise (Failure "Error : please specify -i or -o.")
-  with
-  | Failure msg -> print_string ("Error\n"^msg)
-  | _ -> ()
-    )
-;;
-
-let () =
-    Command.run command
-;;
-
-
-
diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml
index 81d7721e..d41e6081 100644
--- a/ocaml/qp_run.ml
+++ b/ocaml/qp_run.ml
@@ -1,4 +1,3 @@
-open Core
 open Qputils
 
 (* Environment variables :
@@ -10,7 +9,7 @@ open Qputils
 
 let print_list () =
   Lazy.force Qpackage.executables 
-  |> List.iter ~f:(fun (x,_) -> Printf.printf " * %s\n" x) 
+  |> List.iter (fun (x,_) -> Printf.printf " * %s\n" x) 
 
 let () = 
   Random.self_init ()
@@ -27,13 +26,13 @@ let run slave exe ezfio_file =
     in
     let rec try_new_port port_number =
       try 
-        List.iter [ 0;1;2;3;4;5;6;7;8;9 ] ~f:(fun i ->
+        List.iter (fun i ->
             let address = 
               Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) (port_number+i)
             in
             Zmq.Socket.bind dummy_socket address;
             Zmq.Socket.unbind dummy_socket address;
-        );
+        ) [ 0;1;2;3;4;5;6;7;8;9 ] ;
         port_number
       with
       | Unix.Unix_error _ -> try_new_port (port_number+100)
@@ -47,21 +46,21 @@ let run slave exe ezfio_file =
   in
 
   let time_start = 
-    Time.now ()
+    Core.Time.now ()
   in
 
-  if (not (Sys.file_exists_exn ezfio_file)) then
+  if (not (Sys.file_exists ezfio_file)) then
     failwith ("EZFIO directory "^ezfio_file^" not found");
 
   let executables = Lazy.force Qpackage.executables in
-  if (not (List.exists ~f:(fun (x,_) -> x = exe) executables)) then
+  if (not (List.exists (fun (x,_) -> x = exe) executables)) then
     begin
         Printf.printf "\nPossible choices:\n";
-        List.iter executables ~f:(fun (x,_) -> Printf.printf "* %s\n%!" x);
+        List.iter (fun (x,_) -> Printf.printf "* %s\n%!" x) executables;
         failwith ("Executable "^exe^" not found")
     end;
 
-  Printf.printf "%s\n" (Time.to_string time_start);
+  Printf.printf "%s\n" (Core.Time.to_string time_start);
   Printf.printf "===============\nQuantum Package\n===============\n\n";
   Printf.printf "Git Commit: %s\n" Git.message;
   Printf.printf "Git Date  : %s\n" Git.date;
@@ -86,10 +85,10 @@ let run slave exe ezfio_file =
     if slave then
       try
         let address = 
-          In_channel.read_all qp_run_address_filename
-          |> String.strip
+          Core.In_channel.read_all qp_run_address_filename
+          |> String.trim
         in
-        Unix.putenv ~key:"QP_RUN_ADDRESS_MASTER" ~data:address
+        Unix.putenv "QP_RUN_ADDRESS_MASTER" address
       with Sys_error _ -> failwith "No master is not running"
   in
        
@@ -104,23 +103,21 @@ let run slave exe ezfio_file =
   let address =
     Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) port_number
   in
-  Unix.putenv ~key:"QP_RUN_ADDRESS" ~data:address;
+  Unix.putenv "QP_RUN_ADDRESS" address;
   let () = 
     if (not slave) then
-      Out_channel.with_file qp_run_address_filename  ~f:(
-        fun oc -> Out_channel.output_lines oc [address]) 
+      Core.Out_channel.with_file qp_run_address_filename  ~f:(
+        fun oc -> Core.Out_channel.output_lines oc [address]) 
   in
 
 
   (** Run executable *)
   let prefix = 
-    match Sys.getenv "QP_PREFIX" with
-    | Some x -> x^" "
-    | None -> ""
+    try (Sys.getenv "QP_PREFIX")^" " with
+    | Not_found -> ""
   and exe =
-    match (List.find ~f:(fun (x,_) -> x = exe) executables) with
-    | Some (_,x) -> x^" "
-    | None -> assert false
+    match (List.find (fun (x,_) -> x = exe) executables) with
+    | (_,exe) -> exe^" "
   in
   let exit_code = 
     match (Sys.command (prefix^exe^ezfio_file)) with
@@ -133,37 +130,46 @@ let run slave exe ezfio_file =
   if (not slave) then
     Sys.remove qp_run_address_filename;
 
-  let duration = Time.diff (Time.now()) time_start 
-  |> Time.Span.to_string in
+  let duration = Core.Time.diff (Core.Time.now()) time_start 
+  |> Core.Time.Span.to_string in
   Printf.printf "Wall time : %s\n\n" duration;
   if (exit_code <> 0) then
     exit exit_code
 
-let spec = 
-  let open Command.Spec in
-  empty
-  +> flag "slave" no_arg
-     ~doc:(" Required for slave tasks")
-  +> anon ("executable" %: string)
-  +> anon ("ezfio_file" %: string)
-;;
 
 
 
 let () =
-  Command.basic_spec
-  ~summary: "Quantum Package command"
-  ~readme:( fun () -> "
-Executes a Quantum Package binary file among these:\n\n"
+  (* Command-line specs *)
+  "Executes a Quantum Package binary file among these:\n\n"
 ^ (Lazy.force Qpackage.executables
-    |> List.map ~f:(fun (x,_) -> Printf.sprintf " * %s" x )
-    |> String.concat ~sep:"\n" 
-    )
-  )
-  spec
-  (fun slave exe ezfio_file () ->
-    run slave exe ezfio_file
-  )
-  |> Command.run   ~version: Git.sha1   ~build_info: Git.message
+    |> List.map (fun (x,_) -> Printf.sprintf " * %s" x )
+    |> String.concat "\n"
+  ) |> Command_line.set_header_doc;
+
+  [ (
+      '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";
+    ]
+  |> Command_line.set_specs ;
+
+
+  (*  Handle options *)
+  if Command_line.show_help () then
+    exit 0;
+
+  let slave = 
+    match Command_line.get "slave" with
+    | None -> false
+    | _ -> true
+  in
+
+  (* Run the program *)
+  match Command_line.anon_args () with
+  | exe :: ezfio_file :: [] -> run slave exe ezfio_file
+  | _ ->  (Command_line.help () ; failwith "Inconsistent command line")
+  
 
 
diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml
index 819b8ef7..6f235d85 100644
--- a/ocaml/qp_set_mo_class.ml
+++ b/ocaml/qp_set_mo_class.ml
@@ -312,30 +312,40 @@ let default range =
   )
 
 
-let spec =
-  let open Command.Spec in
-  empty 
-  +> flag "core"   (optional string) ~doc:"range Range of core orbitals"
-  +> flag "inact"  (optional string) ~doc:"range Range of inactive orbitals"
-  +> flag "act"    (optional string) ~doc:"range Range of active orbitals"
-  +> flag "virt"   (optional string) ~doc:"range Range of virtual orbitals"
-  +> flag "del"    (optional string) ~doc:"range Range of deleted orbitals"
-  +> flag "q"       no_arg ~doc:" Query: print the current masks"
-  +> anon ("ezfio_filename" %: ezfio_file)
-
-
-let command = 
-    Command.basic_spec
-    ~summary: "Quantum Package command"
-    ~readme:(fun () ->
-     "Set the orbital classes in an EZFIO directory
-      The range of MOs has the form : \"[36-53,72-107,126-131]\"
-        ")
-    spec
-    (fun core inact act virt del q ezfio_filename () -> run ~q ?core ?inact ?act ?virt ?del ezfio_filename )
-
 
 let () =
-    Command.run command
+  "Set the orbital classes in an EZFIO directory.
+The range of MOs has the form : \"[36-53,72-107,126-131]\"." |> Command_line.set_header_doc ;
+
+  [ ( 'c',  "core", "range Range of     core MOs", Command_line.With_opt_arg);
+    ( 'i', "inact", "range Range of inactive MOs", Command_line.With_opt_arg);
+    ( 'a',   "act", "range Range of   active MOs", Command_line.With_opt_arg);
+    ( 'v',  "virt", "range Range of  virtual MOs", Command_line.With_opt_arg);
+    ( 'd',   "del", "range Range of  deleted MOs", Command_line.With_opt_arg);
+    ( 'q', "query", "Print the current MOs classes", Command_line.Without_arg);
+    Command_line.anonymous "<EZFIO_FILE>" "EZFIO directory";
+  ] |> Command_line.set_specs ;
+
+  (*  Handle options *)
+  if Command_line.show_help () then
+    exit 0;
+
+  let core  = Command_line.get   "core" 
+  and inact = Command_line.get  "inact" 
+  and   act = Command_line.get    "act" 
+  and  virt = Command_line.get   "virt" 
+  and   del = Command_line.get    "del" 
+  and     q = Command_line.get_bool "q"
+  in
+
+  let ezfio_filename = 
+    match Command_line.anon_args () with
+    | [x] -> x
+    | _ -> (Command_line.help () ; failwith "EZFIO_FILE is missing")
+  in
+  run ~q ?core ?inact ?act ?virt ?del ezfio_filename
+
+
+
 
 
diff --git a/qpsh b/qpsh
index bea64a09..967f9429 100755
--- a/qpsh
+++ b/qpsh
@@ -3,8 +3,21 @@
 export QP_ROOT=$(dirname $0)
 
 exec bash --init-file <(cat << EOF
-PS1="[\u@\h \W] |\\\${EZFIO_FILE}> "
-source $QP_ROOT/quantum_package.rc
+
+[[ -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_convert_output_to_ezfio b/scripts/ezfio_interface/qp_convert_output_to_ezfio
index 3805323b..b9d5c704 100755
--- a/scripts/ezfio_interface/qp_convert_output_to_ezfio
+++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio
@@ -3,11 +3,11 @@
 convert output of gamess/GAU$$IAN to ezfio
 
 Usage:
-    qp_convert_output_to_ezfio.py <file.out> [-o <ezfio_directory>]
+    qp_convert_output_to_ezfio.py <file.out> [-o|--output <EZFIO_DIRECTORY>]
 
 Option:
     file.out is the file to check (like gamess.out)
-    folder.ezfio is the name you whant for the ezfio
+    EZFIO_DIRECTORY is the name of the produced directory
         (by default is file.out.ezfio)
 
 """
diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template
index a79cea25..b3fdde79 100644
--- a/scripts/ezfio_interface/qp_edit_template
+++ b/scripts/ezfio_interface/qp_edit_template
@@ -132,13 +132,14 @@ let set str s =
 let create_temp_file ezfio_filename fields =
   let temp_filename  = Filename.temp_file "qp_edit_" ".rst" in
   begin
-      Out_channel.with_file temp_filename ~f:(fun out_channel ->
+    Out_channel.with_file temp_filename ~f:(fun out_channel ->
         (file_header ezfio_filename) :: (List.map ~f:get fields) 
         |> String.concat ~sep:"\n" 
         |> Out_channel.output_string out_channel 
-      )
+      );
+    at_exit (fun () -> Sys.remove temp_filename);
+    temp_filename
   end
-  ; temp_filename
 
 
                                                  
@@ -220,88 +221,116 @@ let run check_only ?ndet ?state ezfio_filename =
     In_channel.with_file temp_filename ~f:(fun in_channel ->
       In_channel.input_all in_channel) 
   in
-  List.iter ~f:(fun x -> set temp_string x) tasks;
-
-  (* Remove temp_file *)
-  Sys.remove temp_filename
+  List.iter ~f:(fun x -> set temp_string x) tasks
 
 
 
+
+
+(** Remove the backup file *)
+let remove_backup ezfio_filename =
+  let backup_filename = 
+    Printf.sprintf "%s/work/backup.tar" ezfio_filename
+  in
+  try Sys.remove backup_filename
+  with _ -> ()
+
 (** Create a backup file in case of an exception *)
 let create_backup ezfio_filename =
-  Printf.sprintf "
-    rm -f %s/backup.tgz ;
-    tar -zcf .backup.tgz %s && mv .backup.tgz %s/backup.tgz
-  "
-    ezfio_filename ezfio_filename ezfio_filename
-  |> Sys.command_exn
-
+  remove_backup ezfio_filename;
+  let backup_filename = 
+    Printf.sprintf "%s/work/backup.tar" ezfio_filename
+  in
+  try
+    Printf.sprintf "
+      mkdir -p %s/work ;
+      tar -cf .backup.tar --exclude=\"work/*\" %s && (mv .backup.tar %s || rm .backup.tar)
+    "
+      ezfio_filename ezfio_filename backup_filename
+    |> Sys.command_exn
+  with _ -> ()
 
 
 (** Restore the backup file when an exception occuprs *)
 let restore_backup ezfio_filename =
-  Printf.sprintf "tar -zxf %s/backup.tgz"
-    ezfio_filename 
-  |> Sys.command_exn
+  let filename = 
+    Printf.sprintf "%s/work/backup.tar" ezfio_filename
+  in
+  if Sys.file_exists_exn filename then
+    begin
+      Printf.sprintf "tar -xf %s" filename
+      |> Sys.command_exn;
+      remove_backup ezfio_filename
+    end
 
 
 
-let spec =
-  let open Command.Spec in
-  empty 
-  +> flag "-c" no_arg
-     ~doc:"Checks the input data"
-  +> flag "ndet" (optional int)
-     ~doc:"int Truncate the wavefunction to the target number of determinants"
-  +> flag "state" (optional int)
-     ~doc:"int Pick the state as a new wavefunction."
-  +> anon ("ezfio_file" %: string)
-
-
-let command = 
-    Command.basic_spec
-    ~summary: "Quantum Package command"
-    ~readme:(fun () ->
-      "
-Edit input data
-      ")
-    spec
-(*    (fun i o ezfio_file () -> *)
-    (*fun ezfio_file () -> 
-       try 
-           run ezfio_file
-       with
-       | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
-    *)
-    (fun c ndet state ezfio_file () ->
-       try
-          run c ?ndet ?state ezfio_file ;
-          (* create_backup ezfio_file; *)
-       with
-       | Failure exc 
-       | Invalid_argument exc as e -> 
-         begin
-           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
-           Printf.eprintf "%s\n\n" exc;
-           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
-           (* restore_backup ezfio_file; *) 
-           raise e
-         end
-       | Assert_failure (file, line, ch) as e ->
-         begin
-           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
-           Printf.eprintf "Assert error in file $QP_ROOT/ocaml/%s, line %d, character %d\n\n" file line ch;
-           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
-           (* restore_backup ezfio_file; *)
-           raise e
-         end
-    )
-
 
 let () =
-  Command.run command;
-  exit 0
+  (* Command-line specs *)
+  Command_line.set_header_doc "Opens a text editor to edit the parameters of a EZFIO directory.";
+
+  [ ( '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);
+    Command_line.anonymous "<EZFIO_FILE>" "EZFIO directory";
+  ]
+  |> Command_line.set_specs ;
 
 
+  (*  Handle options *)
+  if Command_line.show_help () then
+    exit 0;
 
+  let ndet =
+    match Command_line.get "ndet" with
+    | None -> None
+    | Some s -> (try Some (int_of_string s)
+            with _ -> failwith "[-n|--ndet] expects an integer")
+  in
+  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")
+  in
+
+  let c =
+    match Command_line.get "check" with
+    | None -> false
+    | _ -> true
+  in
+
+  let ezfio_filename =
+    match Command_line.anon_args () with
+    | [x] -> x
+    | _ -> (Command_line.help () ; failwith "EZFIO_FILE is missing")
+  in
+
+  at_exit (fun () -> remove_backup ezfio_filename);
+  (* Run the program *)
+  try
+    if (not c) then create_backup ezfio_filename;
+    run c ?ndet ?state ezfio_filename
+  with
+  | Failure exc
+  | Invalid_argument exc ->
+    begin
+      Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+      Printf.eprintf "%s\n\n" exc;
+      Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+      restore_backup ezfio_filename;
+      ignore @@ exit 1
+    end
+  | Assert_failure (file, line, ch) ->
+    begin
+      Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+      Printf.eprintf "Assert error in file $QP_ROOT/ocaml/%s, line %d, character %d\n\n" file line ch;
+      Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+      restore_backup ezfio_filename;
+      ignore @@ exit 2
+    end;
+
+    exit 0
 
diff --git a/scripts/qp_set_frozen_core b/scripts/qp_set_frozen_core
index 16e5820f..4f9c6921 100755
--- a/scripts/qp_set_frozen_core
+++ b/scripts/qp_set_frozen_core
@@ -13,10 +13,10 @@ zero.
 
 
 Usage:
-      qp_set_frozen_core [-q] <EZFIO_DIRECTORY>
+      qp_set_frozen_core [-q|--query] <EZFIO_DIRECTORY>
 
 Options:
-    -q        Prints in the standard output the number of core electrons
+    -q --query   Prints in the standard output the number of core electrons
 
 """
 
@@ -58,15 +58,15 @@ def main(arguments):
 
     mo_num = ezfio.mo_basis_mo_num
 
-    if arguments["-q"]:
+    if arguments["--query"]:
         print nb
         sys.exit(0)
 
     if nb == 0:
-        os.system( """qp_set_mo_class -act "[1-%d]" %s""" % 
+        os.system( """qp_set_mo_class -a "[1-%d]" %s""" % 
             (mo_num, sys.argv[1]) )
     else:
-        os.system( """qp_set_mo_class -core "[1-%d]" -act "[%d-%d]" %s""" % 
+        os.system( """qp_set_mo_class -c "[1-%d]" -a "[%d-%d]" %s""" % 
             (nb, nb+1, mo_num, sys.argv[1]) )
 
 
diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f
index 95f613d6..d8f97154 100644
--- a/src/determinants/spindeterminants.irp.f
+++ b/src/determinants/spindeterminants.irp.f
@@ -10,7 +10,7 @@ integer*8 function spin_det_search_key(det,Nint)
   use bitmasks
   implicit none
   BEGIN_DOC
-! Return an integer(8) corresponding to a determinant index for searching
+! Returns an integer(8) corresponding to a determinant index for searching
   END_DOC
   integer, intent(in) :: Nint
   integer(bit_kind), intent(in) :: det(Nint)
@@ -25,32 +25,32 @@ end
 
 
 BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha, (N_int,psi_det_size) ]
- implicit none
- BEGIN_DOC
-! List of $\alpha$ determinants of psi_det
- END_DOC
- integer :: i,k
-
- do i=1,N_det
-   do k=1,N_int
-     psi_det_alpha(k,i) = psi_det(k,1,i)
-   enddo
- enddo
-
+  implicit none
+  BEGIN_DOC
+  ! List of $\alpha$ determinants of psi_det
+  END_DOC
+  integer                        :: i,k
+  
+  do i=1,N_det
+    do k=1,N_int
+      psi_det_alpha(k,i) = psi_det(k,1,i)
+    enddo
+  enddo
+  
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta, (N_int,psi_det_size) ]
- implicit none
- BEGIN_DOC
-! List of $\beta$ determinants of psi_det
- END_DOC
- integer :: i,k
-
- do i=1,N_det
-   do k=1,N_int
-     psi_det_beta(k,i) = psi_det(k,2,i)
-   enddo
- enddo
+  implicit none
+  BEGIN_DOC
+  ! List of $\beta$ determinants of psi_det
+  END_DOC
+  integer                        :: i,k
+  
+  do i=1,N_det
+    do k=1,N_int
+      psi_det_beta(k,i) = psi_det(k,2,i)
+    enddo
+  enddo
 END_PROVIDER
 
 
@@ -58,80 +58,108 @@ BEGIN_TEMPLATE
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_$alpha_unique, (N_int,psi_det_size) ]
 &BEGIN_PROVIDER [ integer, N_det_$alpha_unique ]
- implicit none
- BEGIN_DOC
- ! Unique $\\$alpha$ determinants
- END_DOC
+  implicit none
+  BEGIN_DOC
+  ! Unique $\\$alpha$ determinants
+  END_DOC
+  
+  integer                        :: i,j,k
+  integer, allocatable           :: iorder(:)
+  integer*8, allocatable         :: bit_tmp(:)
+  integer*8                      :: last_key
+  integer*8, external            :: spin_det_search_key
+  logical,allocatable            :: duplicate(:)
+  
+  allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
+  
+  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
+  do i=1,N_det
+    iorder(i) = i
+    bit_tmp(i) = spin_det_search_key(psi_det_$alpha(1,i),N_int)
+  enddo
+  !$OMP END PARALLEL DO
+  
+  
+  call i8sort(bit_tmp,iorder,N_det)
+  
+  N_det_$alpha_unique = 0
+  last_key = 0_8
+  do i=1,N_det
+    last_key = bit_tmp(i)
+    N_det_$alpha_unique += 1
+    do k=1,N_int
+      psi_det_$alpha_unique(k,N_det_$alpha_unique) = psi_det_$alpha(k,iorder(i))
+    enddo
+    duplicate(i) = .False.
+  enddo
+  
+  j=1
+  do i=1,N_det_$alpha_unique-1
+    if (duplicate(i)) then
+      cycle
+    endif
+    j = i+1
+    do while (bit_tmp(j)==bit_tmp(i))
+      if (duplicate(j)) then
+        j += 1
+        cycle
+      endif
+      duplicate(j) = .True.
+      do k=1,N_int
+        if (psi_det_$alpha_unique(k,i) /= psi_det_$alpha_unique(k,j)) then
+          duplicate(j) = .False.
+          exit
+        endif
+      enddo
+      j+=1
+      if (j > N_det_$alpha_unique) then
+        exit
+      endif
+    enddo
+  enddo
+  
+  j=1
+  do i=2,N_det_$alpha_unique
+    if (duplicate(i)) then
+      cycle
+    else
+      j += 1
+      psi_det_$alpha_unique(:,j) = psi_det_$alpha_unique(:,i)
+    endif
+  enddo
+  N_det_$alpha_unique = j
+  
+  call write_int(6,N_det_$alpha_unique,'Number of unique $alpha determinants')
+  
+  deallocate (iorder, bit_tmp, duplicate)
+  
+END_PROVIDER
 
- integer                        :: i,j,k
- integer, allocatable           :: iorder(:)
- integer*8, allocatable         :: bit_tmp(:)
- integer*8                      :: last_key
- integer*8, external            :: spin_det_search_key
- logical,allocatable            :: duplicate(:)
+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
 
- allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
- 
- !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
- do i=1,N_det
-   iorder(i) = i
-   bit_tmp(i) = spin_det_search_key(psi_det_$alpha(1,i),N_int)
- enddo
- !$OMP END PARALLEL DO
-
-
- call i8sort(bit_tmp,iorder,N_det)
-
- N_det_$alpha_unique = 0
- last_key = 0_8
- do i=1,N_det
-   last_key = bit_tmp(i)
-   N_det_$alpha_unique += 1
-   do k=1,N_int
-     psi_det_$alpha_unique(k,N_det_$alpha_unique) = psi_det_$alpha(k,iorder(i))
-   enddo
-   duplicate(i) = .False.
- enddo
-
- j=1
- do i=1,N_det_$alpha_unique-1
-   if (duplicate(i)) then
-     cycle
-   endif
-   j = i+1
-   do while (bit_tmp(j)==bit_tmp(i))
-     if (duplicate(j)) then
-       j += 1
-       cycle
-     endif
-     duplicate(j) = .True.
-     do k=1,N_int
-       if (psi_det_$alpha_unique(k,i) /= psi_det_$alpha_unique(k,j)) then
-         duplicate(j) = .False.
-         exit
-       endif
-     enddo
-     j+=1
-     if (j > N_det_$alpha_unique) then
-       exit
-     endif
-   enddo
- enddo
-
- j=1
- do i=2,N_det_$alpha_unique
-   if (duplicate(i)) then
-     cycle
-  else
-    j += 1
-    psi_det_$alpha_unique(:,j) = psi_det_$alpha_unique(:,i)
-  endif
- enddo
- N_det_$alpha_unique = j
-
- call write_int(6,N_det_$alpha_unique,'Number of unique $alpha determinants')
-
- deallocate (iorder, bit_tmp, duplicate)
+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
 
@@ -148,29 +176,29 @@ 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 ``psi_det_alpha_unique`` array
   END_DOC
   implicit none
-
+  
   integer, intent(in)            :: Nint
   integer(bit_kind), intent(in)  :: key(Nint)
-
+  
   integer                        :: i, ibegin, iend, istep, l
-  integer(8)                      :: det_ref, det_search
-  integer(8), external            :: spin_det_search_key
+  integer(8)                     :: det_ref, det_search
+  integer(8), external           :: spin_det_search_key
   logical                        :: in_wavefunction
-
+  
   in_wavefunction = .False.
   get_index_in_psi_det_alpha_unique = 0
   ibegin = 1
   iend   = N_det_alpha_unique + 1
-
+  
   !DIR$ FORCEINLINE
   det_ref = spin_det_search_key(key,Nint)
-
+  
   !DIR$ FORCEINLINE
   det_search = spin_det_search_key(psi_det_alpha_unique(1,1),Nint)
-
+  
   istep = shiftr(iend-ibegin,1)
   i=ibegin+istep
   do while (istep > 0)
@@ -186,7 +214,7 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint)
     istep = shiftr(iend-ibegin,1)
     i = ibegin + istep
   end do
-
+  
   !DIR$ FORCEINLINE
   do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref)
     i = i-1
@@ -195,9 +223,9 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint)
     endif
   enddo
   i += 1
-
+  
   ASSERT (i <= N_det_alpha_unique)
-
+  
   !DIR$ FORCEINLINE
   do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref)
     if (key(1) /= psi_det_alpha_unique(1,i)) then
@@ -221,36 +249,36 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint)
       ASSERT (get_index_in_psi_det_alpha_unique > 0)
       return
     endif
-
+    
   enddo
-
+  
 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 ``psi_det_beta_unique`` array
   END_DOC
   implicit none
-
+  
   integer, intent(in)            :: Nint
   integer(bit_kind), intent(in)  :: key(Nint)
-
+  
   integer                        :: i, ibegin, iend, istep, l
-  integer(8)                      :: det_ref, det_search
-  integer(8), external            :: spin_det_search_key
+  integer(8)                     :: det_ref, det_search
+  integer(8), external           :: spin_det_search_key
   logical                        :: in_wavefunction
-
+  
   in_wavefunction = .False.
   get_index_in_psi_det_beta_unique = 0
   ibegin = 1
   iend   = N_det_beta_unique + 1
-
+  
   !DIR$ FORCEINLINE
   det_ref = spin_det_search_key(key,Nint)
   !DIR$ FORCEINLINE
   det_search = spin_det_search_key(psi_det_beta_unique(1,1),Nint)
-
+  
   istep = shiftr(iend-ibegin,1)
   i=ibegin+istep
   do while (istep > 0)
@@ -266,7 +294,7 @@ integer function get_index_in_psi_det_beta_unique(key,Nint)
     istep = shiftr(iend-ibegin,1)
     i = ibegin + istep
   end do
-
+  
   !DIR$ FORCEINLINE
   do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref)
     i = i-1
@@ -275,9 +303,9 @@ integer function get_index_in_psi_det_beta_unique(key,Nint)
     endif
   enddo
   i += 1
-
+  
   ASSERT (i <= N_det_beta_unique)
-
+  
   !DIR$ FORCEINLINE
   do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref)
     if (key(1) /= psi_det_beta_unique(1,i)) then
@@ -301,22 +329,22 @@ integer function get_index_in_psi_det_beta_unique(key,Nint)
       ASSERT (get_index_in_psi_det_beta_unique > 0)
       return
     endif
-
+    
   enddo
-
+  
 end
 
 
 subroutine write_spindeterminants
   use bitmasks
   implicit none
-  integer(8), allocatable         :: tmpdet(:,:)
+  integer(8), allocatable        :: tmpdet(:,:)
   integer                        :: N_int2
   integer                        :: i,j,k
-  integer(8)                      :: det_8(100)
+  integer(8)                     :: det_8(100)
   integer(bit_kind)              :: det_bk((100*8)/bit_kind)
   equivalence (det_8, det_bk)
-
+  
   N_int2 = (N_int*bit_kind)/8
   call ezfio_set_spindeterminants_n_det_alpha(N_det_alpha_unique)
   call ezfio_set_spindeterminants_n_det_beta(N_det_beta_unique)
@@ -324,7 +352,7 @@ subroutine write_spindeterminants
   call ezfio_set_spindeterminants_n_int(N_int)
   call ezfio_set_spindeterminants_bit_kind(bit_kind)
   call ezfio_set_spindeterminants_n_states(N_states)
-
+  
   allocate(tmpdet(N_int2,N_det_alpha_unique))
   do i=1,N_det_alpha_unique
     do k=1,N_int
@@ -357,31 +385,31 @@ end
 
  BEGIN_PROVIDER [ double precision, det_alpha_norm, (N_det_alpha_unique) ]
 &BEGIN_PROVIDER [ double precision, det_beta_norm, (N_det_beta_unique) ]
- implicit none
- BEGIN_DOC
- ! Norm of the $\alpha$ and $\beta$ spin determinants in the wave function:
- !
- ! ||Da||_i \sum_j C_{ij}**2
- END_DOC
-
- integer :: i,j,k,l
- double precision :: f
-
- det_alpha_norm = 0.d0
- det_beta_norm  = 0.d0
- do k=1,N_det
-   i = psi_bilinear_matrix_rows(k)
-   j = psi_bilinear_matrix_columns(k)
-   f = 0.d0
-   do l=1,N_states
-    f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l) * state_average_weight(l)
-   enddo
-   det_alpha_norm(i) += f
-   det_beta_norm(j)  += f
- enddo
- det_alpha_norm = det_alpha_norm 
- det_beta_norm = det_beta_norm 
-
+  implicit none
+  BEGIN_DOC
+  ! Norm of the $\alpha$ and $\beta$ spin determinants in the wave function:
+  !
+  ! $||D_\alpha||_i = \sum_j C_{ij}^2$
+  END_DOC
+  
+  integer                        :: i,j,k,l
+  double precision               :: f
+  
+  det_alpha_norm = 0.d0
+  det_beta_norm  = 0.d0
+  do k=1,N_det
+    i = psi_bilinear_matrix_rows(k)
+    j = psi_bilinear_matrix_columns(k)
+    f = 0.d0
+    do l=1,N_states
+      f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l) * state_average_weight(l)
+    enddo
+    det_alpha_norm(i) += f
+    det_beta_norm(j)  += f
+  enddo
+  det_alpha_norm = det_alpha_norm
+  det_beta_norm = det_beta_norm
+  
 END_PROVIDER
 
 
@@ -391,41 +419,41 @@ END_PROVIDER
 !                                                                              !
 !==============================================================================!
 
-BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ]
+ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order  , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
-!  D_a^t C D_b
-!
-! Rows are $\alpha$ determinants and columns are $\beta.$
-!
-! Order refers to psi_det
+  ! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
+  !  $D_\alpha^\dagger.C.D_\beta$
+  !
+  ! Rows are $\alpha$ determinants and columns are $\beta$.
+  !
+  ! Order refers to psi_det
   END_DOC
   integer                        :: i,j,k, l
-  integer(bit_kind)               :: tmp_det(N_int,2)
+  integer(bit_kind)              :: tmp_det(N_int,2)
   integer, external              :: get_index_in_psi_det_sorted_bit
-
-
+  
+  
   PROVIDE psi_coef_sorted_bit
-
-  integer*8, allocatable :: to_sort(:)
-  integer, external :: get_index_in_psi_det_alpha_unique
-  integer, external :: get_index_in_psi_det_beta_unique
+  
+  integer*8, allocatable         :: to_sort(:)
+  integer, external              :: get_index_in_psi_det_alpha_unique
+  integer, external              :: get_index_in_psi_det_beta_unique
   allocate(to_sort(N_det))
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l)
   do k=1,N_det
     i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int)
     ASSERT (i>0)
     ASSERT (i<=N_det_alpha_unique)
-
+    
     j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int)
     ASSERT (j>0)
     ASSERT (j<=N_det_beta_unique)
-
+    
     do l=1,N_states
       psi_bilinear_matrix_values(k,l) = psi_coef(k,l)
     enddo
@@ -437,7 +465,7 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
   enddo
   !$OMP END PARALLEL DO
   call i8sort(to_sort, psi_bilinear_matrix_order, N_det)
-  !$OMP PARALLEL 
+  !$OMP PARALLEL
   !$OMP SINGLE
   call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det)
   !$OMP END SINGLE
@@ -457,7 +485,7 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
   ASSERT (maxval(psi_bilinear_matrix_rows) == N_det_alpha_unique)
   ASSERT (maxval(psi_bilinear_matrix_columns) == N_det_beta_unique)
   ASSERT (maxval(psi_bilinear_matrix_order) == N_det)
-
+  
 END_PROVIDER
 
 
@@ -465,7 +493,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 ``psi_bilinear_matrix`` to ``psi_det``
   END_DOC
   integer                        :: k
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k)
@@ -482,16 +510,16 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1)
   use bitmasks
   implicit none
   BEGIN_DOC
-! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
-!
-! $D_\alpha^\dagger.C.D_\beta$
-!
-! Rows are $\alpha$ determinants and columns are $\beta.$
-!
-! Order refers to psi_det
+  ! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
+  !
+  ! $D_\alpha^\dagger.C.D_\beta$
+  !
+  ! Rows are $\alpha$ determinants and columns are $\beta$.
+  !
+  ! Order refers to ``psi_det``
   END_DOC
   integer                        :: i,j,k, l
-
+  
   l = psi_bilinear_matrix_columns(1)
   psi_bilinear_matrix_columns_loc(l) = 1
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k,l)
@@ -503,35 +531,35 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1)
       psi_bilinear_matrix_columns_loc(l) = k
     endif
     if (psi_bilinear_matrix_columns(k) < 1) then
-     stop '(psi_bilinear_matrix_columns(k) < 1)'
+      stop '(psi_bilinear_matrix_columns(k) < 1)'
     endif
   enddo
   !$OMP END PARALLEL DO
   psi_bilinear_matrix_columns_loc(N_det_beta_unique+1) = N_det+1
   ASSERT (minval(psi_bilinear_matrix_columns_loc) == 1)
   ASSERT (maxval(psi_bilinear_matrix_columns_loc) == N_det+1)
-
+  
 END_PROVIDER
 
-BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
+ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order  , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-! Transpose of psi_bilinear_matrix
-!
-! $D_\beta^\dagger.C^\dagger.D_\alpha$
-!
-! Rows are $\alpha$ determinants and columns are $\beta$, but the matrix is stored in row major
-! format
+  ! Transpose of ``psi_bilinear_matrix``
+  !
+  ! $D_\beta^\dagger.C^\dagger.D_\alpha$
+  !
+  ! Rows are $\alpha$ determinants and columns are $\beta$, but the matrix is stored in row major
+  ! format.
   END_DOC
   integer                        :: i,j,k,l
-
+  
   PROVIDE psi_coef_sorted_bit
-
-  integer*8, allocatable :: to_sort(:)
+  
+  integer*8, allocatable         :: to_sort(:)
   allocate(to_sort(N_det))
   !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l)
   do l=1,N_states
@@ -541,17 +569,17 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
     enddo
     !$OMP ENDDO NOWAIT
   enddo
-  !$OMP DO 
+  !$OMP DO
   do k=1,N_det
     psi_bilinear_matrix_transp_columns(k) = psi_bilinear_matrix_columns(k)
     ASSERT (psi_bilinear_matrix_transp_columns(k) > 0)
     ASSERT (psi_bilinear_matrix_transp_columns(k) <= N_det)
-
+    
     psi_bilinear_matrix_transp_rows   (k) = psi_bilinear_matrix_rows   (k)
     ASSERT (psi_bilinear_matrix_transp_rows(k) > 0)
     ASSERT (psi_bilinear_matrix_transp_rows(k) <= N_det)
-
-    i = psi_bilinear_matrix_transp_columns(k) 
+    
+    i = psi_bilinear_matrix_transp_columns(k)
     j = psi_bilinear_matrix_transp_rows   (k)
     to_sort(k) = int(N_det_beta_unique,8) * int(j-1,8) + int(i,8)
     ASSERT (to_sort(k) > 0)
@@ -574,17 +602,17 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
   ASSERT (maxval(psi_bilinear_matrix_transp_columns) == N_det_beta_unique)
   ASSERT (maxval(psi_bilinear_matrix_transp_rows) == N_det_alpha_unique)
   ASSERT (maxval(psi_bilinear_matrix_transp_order) == N_det)
-
+  
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows_loc, (N_det_alpha_unique+1) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-! Location of the columns in the psi_bilinear_matrix
+  ! Location of the columns in the ``psi_bilinear_matrix``
   END_DOC
   integer                        :: i,j,k, l
-
+  
   l = psi_bilinear_matrix_transp_rows(1)
   psi_bilinear_matrix_transp_rows_loc(l) = 1
   !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k,l)
@@ -600,14 +628,15 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows_loc, (N_det_alpha_uniq
   psi_bilinear_matrix_transp_rows_loc(N_det_alpha_unique+1) = N_det+1
   ASSERT (minval(psi_bilinear_matrix_transp_rows_loc) == 1)
   ASSERT (maxval(psi_bilinear_matrix_transp_rows_loc) == N_det+1)
-
+  
 END_PROVIDER
 
 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 ``psi_bilinear_matrix_order_transp`` to
+  ! ``psi_bilinear_matrix``
   END_DOC
   integer                        :: k
   psi_bilinear_matrix_order_transp_reverse = -1
@@ -624,11 +653,11 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_det_beta_unique,N_states) ]
   implicit none
   BEGIN_DOC
-! Coefficient matrix if the wave function is expressed in a bilinear form :
-!
-! $D_\alpha^\dagger.C.D_\beta$
+  ! Coefficient matrix if the wave function is expressed in a bilinear form :
+  !
+  ! $D_\alpha^\dagger.C.D_\beta$
   END_DOC
-  integer :: i,j,k,istate
+  integer                        :: i,j,k,istate
   psi_bilinear_matrix = 0.d0
   do k=1,N_det
     i = psi_bilinear_matrix_rows(k)
@@ -643,8 +672,8 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
   use bitmasks
   implicit none
   BEGIN_DOC
-! Generate a wave function containing all possible products 
-! of $\alpha$ and $\beta$ determinants
+  ! Generates a wave function containing all possible products
+  ! of $\alpha$ and $\beta$ determinants
   END_DOC
   logical, intent(in)            :: truncate
   integer                        :: i,j,k
@@ -653,7 +682,7 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
   integer, external              :: get_index_in_psi_det_sorted_bit
   double precision               :: norm(N_states)
   PROVIDE psi_bilinear_matrix
-
+  
   call generate_all_alpha_beta_det_products
   norm = 0.d0
   !$OMP PARALLEL DO DEFAULT(NONE)                                    &
@@ -673,7 +702,7 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
       idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int)
       if (idx > 0) then
         do k=1,N_states
-          psi_coef_sorted_bit(idx,k) = psi_bilinear_matrix(i,j,k) 
+          psi_coef_sorted_bit(idx,k) = psi_bilinear_matrix(i,j,k)
           !$OMP ATOMIC
           norm(k) += psi_bilinear_matrix(i,j,k)*psi_bilinear_matrix(i,j,k)
         enddo
@@ -681,7 +710,7 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
     enddo
   enddo
   !$OMP END PARALLEL DO
-
+  
   do k=1,N_states
     norm(k) = 1.d0/dsqrt(norm(k))
     do i=1,N_det
@@ -707,13 +736,13 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
   enddo
   N_det = min(i,N_det)
   SOFT_TOUCH psi_det psi_coef N_det
-
+  
 end
 
 subroutine generate_all_alpha_beta_det_products
   implicit none
   BEGIN_DOC
-!  Create a wave function from all possible $\alpha \times \beta$ determinants
+  ! Creates a wave function from all possible $\alpha \times \beta$ determinants
   END_DOC
   integer                        :: i,j,k,l
   integer                        :: iproc
@@ -721,10 +750,10 @@ subroutine generate_all_alpha_beta_det_products
   integer(bit_kind), allocatable :: tmp_det(:,:,:)
   logical, external              :: is_in_wavefunction
   PROVIDE H_apply_buffer_allocated
-
+  
   !$OMP PARALLEL DEFAULT(NONE) SHARED(psi_coef_sorted_bit,N_det_beta_unique,&
       !$OMP N_det_alpha_unique, N_int, psi_det_alpha_unique, psi_det_beta_unique,&
-      !$OMP N_det)                                                &
+      !$OMP N_det)                                                   &
       !$OMP PRIVATE(i,j,k,l,tmp_det,iproc)
   !$ iproc = omp_get_thread_num()
   allocate (tmp_det(N_int,2,N_det_alpha_unique))
@@ -742,7 +771,7 @@ subroutine generate_all_alpha_beta_det_products
     enddo
     call fill_H_apply_buffer_no_selection(l-1, tmp_det, N_int, iproc)
   enddo
-  !$OMP END DO 
+  !$OMP END DO
   deallocate(tmp_det)
   !$OMP END PARALLEL
   call copy_H_apply_buffer_to_wf
@@ -756,12 +785,14 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single and double excitations in the list of
-! unique $\alpha$ determinants.
-!
-! /!\ : The buffer is transposed !
-!
+  !
+  ! Returns the indices of all the single and double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
+  ! ..warning::
+  ! 
+  !     The buffer is transposed.
+  !
   END_DOC
   integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
@@ -770,7 +801,7 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf
   integer, intent(out)           :: doubles(size_buffer)
   integer, intent(out)           :: n_singles
   integer, intent(out)           :: n_doubles
-
+  
   select case (Nint)
     case (1)
       call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles)
@@ -780,10 +811,10 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf
       call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
     case (4)
       call get_all_spin_singles_and_doubles_4(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
-    case default
+      case default
       call get_all_spin_singles_and_doubles_N_int(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
   end select
-
+  
 end
 
 
@@ -791,17 +822,17 @@ subroutine get_all_spin_singles(buffer, idx, spindet, Nint, size_buffer, singles
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the single excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
   integer(bit_kind), intent(in)  :: spindet(Nint)
   integer, intent(out)           :: singles(size_buffer)
   integer, intent(out)           :: n_singles
-
+  
   select case (N_int)
     case (1)
       call get_all_spin_singles_1(buffer, idx, spindet(1), size_buffer, singles, n_singles)
@@ -812,7 +843,7 @@ subroutine get_all_spin_singles(buffer, idx, spindet, Nint, size_buffer, singles
       call get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
     case (4)
       call get_all_spin_singles_4(buffer, idx, spindet, size_buffer, singles, n_singles)
-    case default
+      case default
       call get_all_spin_singles_N_int(buffer, idx, spindet, size_buffer, singles, n_singles)
   end select
   
@@ -823,17 +854,17 @@ subroutine get_all_spin_doubles(buffer, idx, spindet, Nint, size_buffer, doubles
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the double excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
   integer(bit_kind), intent(in)  :: spindet(Nint)
   integer, intent(out)           :: doubles(size_buffer)
   integer, intent(out)           :: n_doubles
-
+  
   select case (N_int)
     case (1)
       call get_all_spin_doubles_1(buffer, idx, spindet(1), size_buffer, doubles, n_doubles)
@@ -843,10 +874,10 @@ subroutine get_all_spin_doubles(buffer, idx, spindet, Nint, size_buffer, doubles
       call get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
     case (4)
       call get_all_spin_doubles_4(buffer, idx, spindet, size_buffer, doubles, n_doubles)
-    case default
+      case default
       call get_all_spin_doubles_N_int(buffer, idx, spindet, size_buffer, doubles, n_doubles)
   end select
-
+  
 end
 
 
@@ -856,7 +887,7 @@ end
 subroutine copy_psi_bilinear_to_psi(psi, isize)
   implicit none
   BEGIN_DOC
-! Overwrites psi_det and psi_coef with the wf in bilinear order
+  ! Overwrites ``psi_det`` and ``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)
@@ -870,77 +901,77 @@ subroutine copy_psi_bilinear_to_psi(psi, isize)
 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)
+  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
- END_DOC
- integer                        :: i,j
- integer, allocatable           :: idx0(:), s(:)
- allocate (idx0(N_det_alpha_unique))
- do i=1, N_det_alpha_unique
-   idx0(i) = i
- enddo
-
- !$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   PRIVATE(i,s,j)
- allocate (s(singles_alpha_size))
- !$OMP DO SCHEDULE(static,1)
- 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, &
-       N_det_alpha_unique, s, j)
-   singles_alpha_csc_idx(i+1) = int(j,8)
- enddo
- !$OMP END DO
- deallocate(s)
- !$OMP END PARALLEL 
- deallocate(idx0)
-
- singles_alpha_csc_idx(1) = 1_8
- do i=2, N_det_alpha_unique+1
-   singles_alpha_csc_idx(i) = singles_alpha_csc_idx(i) + singles_alpha_csc_idx(i-1)
- enddo
- singles_alpha_csc_size = singles_alpha_csc_idx(N_det_alpha_unique+1)
+  implicit none
+  BEGIN_DOC
+  ! Dimension of the ``singles_alpha`` array
+  END_DOC
+  integer                        :: i,j
+  integer, allocatable           :: idx0(:), s(:)
+  allocate (idx0(N_det_alpha_unique))
+  do i=1, N_det_alpha_unique
+    idx0(i) = i
+  enddo
+  
+  !$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   PRIVATE(i,s,j)
+  allocate (s(singles_alpha_size))
+  !$OMP DO SCHEDULE(static,1)
+  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,&
+        N_det_alpha_unique, s, j)
+    singles_alpha_csc_idx(i+1) = int(j,8)
+  enddo
+  !$OMP END DO
+  deallocate(s)
+  !$OMP END PARALLEL
+  deallocate(idx0)
+  
+  singles_alpha_csc_idx(1) = 1_8
+  do i=2, N_det_alpha_unique+1
+    singles_alpha_csc_idx(i) = singles_alpha_csc_idx(i) + singles_alpha_csc_idx(i-1)
+  enddo
+  singles_alpha_csc_size = singles_alpha_csc_idx(N_det_alpha_unique+1)
 END_PROVIDER
 
 
 BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ]
- implicit none
- BEGIN_DOC
- ! Dimension of the singles_alpha array
- END_DOC
- integer                        :: i, k
- integer, allocatable           :: idx0(:)
- allocate (idx0(N_det_alpha_unique))
- do i=1, N_det_alpha_unique
-   idx0(i) = i
- enddo
-
- !$OMP PARALLEL DO DEFAULT(NONE) &
- !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique, &
- !$OMP          idx0, N_int, singles_alpha_csc, singles_alpha_csc_idx) &
- !$OMP   PRIVATE(i,k) SCHEDULE(static,1)
- 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, &
-       N_det_alpha_unique, singles_alpha_csc(singles_alpha_csc_idx(i)), &
-       k)
- enddo
- !$OMP END PARALLEL DO
- deallocate(idx0)
-
+  implicit none
+  BEGIN_DOC
+  ! Dimension of the singles_alpha array
+  END_DOC
+  integer                        :: i, k
+  integer, allocatable           :: idx0(:)
+  allocate (idx0(N_det_alpha_unique))
+  do i=1, N_det_alpha_unique
+    idx0(i) = i
+  enddo
+  
+  !$OMP PARALLEL DO DEFAULT(NONE)                                    &
+      !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique,       &
+      !$OMP          idx0, N_int, singles_alpha_csc, singles_alpha_csc_idx)&
+      !$OMP   PRIVATE(i,k) SCHEDULE(static,1)
+  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,&
+        N_det_alpha_unique, singles_alpha_csc(singles_alpha_csc_idx(i)),&
+        k)
+  enddo
+  !$OMP END PARALLEL DO
+  deallocate(idx0)
+  
 END_PROVIDER
 
 
@@ -950,12 +981,12 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer,
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single and double excitations in the list of
-! unique $\alpha$ determinants.
-!
-! /!\ : The buffer is transposed !
-!
+  !
+  ! Returns the indices of all the single and double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
+  ! /!\ : The buffer is transposed !
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(size_buffer)
@@ -964,12 +995,12 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer,
   integer, intent(out)           :: doubles(size_buffer)
   integer, intent(out)           :: n_singles
   integer, intent(out)           :: n_doubles
-
+  
   integer                        :: i
   include 'utils/constants.include.F'
   integer                        :: degree
-
-
+  
+  
   n_singles = 1
   n_doubles = 1
   do i=1,size_buffer
@@ -993,10 +1024,10 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the single excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(size_buffer)
@@ -1007,10 +1038,10 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_
   integer(bit_kind)              :: v
   integer                        :: degree
   include 'utils/constants.include.F'
-
+  
   n_singles = 1
   do i=1,size_buffer
-    degree = popcnt(xor( spindet, buffer(i) )) 
+    degree = popcnt(xor( spindet, buffer(i) ))
     if (degree == 2) then
       singles(n_singles) = idx(i)
       n_singles = n_singles+1
@@ -1025,10 +1056,10 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the double excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer(size_buffer)
@@ -1038,7 +1069,7 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_
   integer                        :: i
   include 'utils/constants.include.F'
   integer                        :: degree
-
+  
   n_doubles = 1
   do i=1,size_buffer
     degree = popcnt(xor( spindet, buffer(i) ))
@@ -1059,12 +1090,12 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single and double excitations in the list of
-! unique $\alpha$ determinants.
-!
-! /!\ : The buffer is transposed !
-!
+  !
+  ! Returns the indices of all the single and double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
+  ! /!\ : The buffer is transposed !
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer($N_int,size_buffer)
@@ -1073,17 +1104,17 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
   integer, intent(out)           :: doubles(size_buffer)
   integer, intent(out)           :: n_singles
   integer, intent(out)           :: n_doubles
-
+  
   integer                        :: i,k
   integer(bit_kind)              :: xorvec($N_int)
   integer                        :: degree
-
+  
   n_singles = 1
   n_doubles = 1
   do i=1,size_buffer
-
+    
     do k=1,$N_int
-        xorvec(k) = xor( spindet(k), buffer(k,i) )
+      xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
     
     if (xorvec(1) /= 0_8) then
@@ -1091,13 +1122,13 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
     else
       degree = 0
     endif
-  
+    
     do k=2,$N_int
       if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
       endif
     enddo
-
+    
     if ( degree == 4 ) then
       doubles(n_doubles) = idx(i)
       n_doubles = n_doubles+1
@@ -1105,7 +1136,7 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
       singles(n_singles) = idx(i)
       n_singles = n_singles+1
     endif
-
+    
   enddo
   n_singles = n_singles-1
   n_doubles = n_doubles-1
@@ -1117,27 +1148,27 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the single excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the single excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer($N_int,size_buffer)
   integer(bit_kind), intent(in)  :: spindet($N_int)
   integer, intent(out)           :: singles(size_buffer)
   integer, intent(out)           :: n_singles
-
+  
   integer                        :: i,k
   include 'utils/constants.include.F'
   integer(bit_kind)              :: xorvec($N_int)
   integer                        :: degree
-
+  
   n_singles = 1
   do i=1,size_buffer
-
+    
     do k=1,$N_int
-       xorvec(k) = xor( spindet(k), buffer(k,i) )
+      xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
     
     if (xorvec(1) /= 0_8) then
@@ -1145,18 +1176,18 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single
     else
       degree = 0
     endif
-
+    
     do k=2,$N_int
       if ( (degree <= 2).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
       endif
     enddo
-
+    
     if ( degree == 2 ) then
-        singles(n_singles) = idx(i)
-        n_singles = n_singles+1
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
     endif
-
+    
   enddo
   n_singles = n_singles-1
   
@@ -1167,24 +1198,24 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double
   use bitmasks
   implicit none
   BEGIN_DOC
-!
-! Returns the indices of all the double excitations in the list of
-! unique $\alpha$ determinants.
-!
+  !
+  ! Returns the indices of all the double excitations in the list of
+  ! unique $\alpha$ determinants.
+  !
   END_DOC
   integer, intent(in)            :: size_buffer, idx(size_buffer)
   integer(bit_kind), intent(in)  :: buffer($N_int,size_buffer)
   integer(bit_kind), intent(in)  :: spindet($N_int)
   integer, intent(out)           :: doubles(size_buffer)
   integer, intent(out)           :: n_doubles
-
+  
   integer                        :: i,k, degree
   include 'utils/constants.include.F'
   integer(bit_kind)              :: xorvec($N_int)
-
+  
   n_doubles = 1
   do i=1,size_buffer
-
+    
     do k=1,$N_int
       xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
@@ -1194,20 +1225,20 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double
     else
       degree = 0
     endif
-  
+    
     do k=2,$N_int
       if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
       endif
     enddo
-
+    
     if ( degree == 4 ) then
       doubles(n_doubles) = idx(i)
       n_doubles = n_doubles+1
     endif
-
+    
   enddo
-
+  
   n_doubles = n_doubles-1
   
 end
@@ -1225,8 +1256,8 @@ subroutine wf_of_psi_bilinear_matrix(truncate)
   use bitmasks
   implicit none
   BEGIN_DOC
-! Generate a wave function containing all possible products 
-! of $\alpha$ and $\beta$ determinants
+  ! Generate a wave function containing all possible products
+  ! of $\alpha$ and $\beta$ determinants
   END_DOC
   logical, intent(in)            :: truncate
   integer                        :: i,j,k
@@ -1234,22 +1265,22 @@ subroutine wf_of_psi_bilinear_matrix(truncate)
   integer                        :: idx
   integer, external              :: get_index_in_psi_det_sorted_bit
   double precision               :: norm(N_states)
-
+  
   do k=1,N_det
-   i = psi_bilinear_matrix_rows(k)
-   j = psi_bilinear_matrix_columns(k)
-   psi_det(1:N_int,1,k) = psi_det_alpha_unique(1:N_int,i)
-   psi_det(1:N_int,2,k) = psi_det_beta_unique (1:N_int,j)
+    i = psi_bilinear_matrix_rows(k)
+    j = psi_bilinear_matrix_columns(k)
+    psi_det(1:N_int,1,k) = psi_det_alpha_unique(1:N_int,i)
+    psi_det(1:N_int,2,k) = psi_det_beta_unique (1:N_int,j)
   enddo
   psi_coef(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
-  TOUCH psi_det psi_coef 
-
+  TOUCH psi_det psi_coef
+  
   psi_det   = psi_det_sorted
   psi_coef  = psi_coef_sorted
   do while (sum( dabs(psi_coef(N_det,1:N_states)) ) == 0.d0)
     N_det -= 1
   enddo
   SOFT_TOUCH psi_det psi_coef N_det
-
+  
 end
 
diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats
index 4d6983e1..eed81922 100644
--- a/src/hartree_fock/10.hf.bats
+++ b/src/hartree_fock/10.hf.bats
@@ -20,7 +20,7 @@ function run() {
 }
 
 @test "SO" { # 0.539000
-  run so.ezfio -25.71752633718843
+  run so.ezfio -25.7175126082701
 }
 
 @test "HCO" { # 0.636700