Merge branch 'LCPQ-master'

2024-06-29 16:34:50 +02:00 · 2015-11-07 01:38:37 +01:00 · 2015-11-07 01:38:37 +01:00 · 8b5ad98793
commit 8b5ad98793
parent f82605fb8a ca780fa6bd
22 changed files with 1308 additions and 863 deletions
--- a/.travis.yml
+++ b/.travis.yml
@ -27,7 +27,7 @@ python:
 script: 
  - ./configure  --production ./config/gfortran.cfg
  - source ./quantum_package.rc
-  - qp_install_module.py install Full_CI Hartree_Fock
+  - qp_install_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD
  - ninja
  - cd ocaml ; make ; cd -
  - cd testing_no_regression ; ./unit_test.py
--- a/README.md
+++ b/README.md
@ -30,7 +30,7 @@ For example you can type `./configure config/gfortran.cfg --production`
 This command have to purpose :
 - Download and install all the requirements.
-   Installing OCaml and the Core library may take somme time (up to 20min on an old machine).
+   Installing OCaml and the Core library may take some time (up to 20min on an old machine).
 - And create the file which contains all the tree dependencies for the binaries.  
   It's not a Makefile, but a Ninja file (so don't type `make` is hopeless, type `ninja` instead)
@ -173,4 +173,4 @@ You have two or more ezfio configuration file for the same variable. Check in `$
 #### Fix
    - rm $QP_ROOT/install/EZFIO/config/*
-    - ninja 
+    - ninja 
--- a/config/gfortran.cfg
+++ b/config/gfortran.cfg
@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -ffree-line-length-none -I .
+FC           : gfortran -g -ffree-line-length-none -I .
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32
--- a/ocaml/.gitignore
+++ b/ocaml/.gitignore
@ -1,50 +1,47 @@
 .gitignore
 ezfio.ml
 Qptypes.ml
 qptypes_generator.byte
 _build
-qp_basis_clean.native
+ezfio.ml
-qp_create_ezfio_from_xyz.native
+.gitignore
-qp_edit.native
+Input_auto_generated.ml
 qp_print.native
 qp_run.native
 qp_set_ddci.native
 qp_set_mo_class.native
 qp_edit.native
 test_atom.byte
 test_basis.byte
 test_bitlist.byte
 test_determinants.byte
 test_elements.byte
 test_excitation.byte
 test_gto.byte
 test_mo_label.byte
 test_molecule.byte
 test_point3d.byte
 test_atom
 test_basis
 test_bitlist
 test_determinants
 test_elements
 test_excitation
 test_gto
 test_mo_label
 test_molecule
 test_point3d
 qp_basis_clean
 qp_create_ezfio_from_xyz
 qp_edit
 qp_print
 qp_run
 qp_set_ddci
 qp_set_mo_class
 Input_determinants.ml
 Input_hartree_fock.ml
 Input_integrals_bielec.ml
 Input_perturbation.ml
 Input_properties.ml
 Input_pseudo.ml
-qp_edit.ml
+qp_basis_clean
 qp_basis_clean.native
 qp_create_ezfio_from_xyz
 qp_create_ezfio_from_xyz.native
 qp_edit
 qp_edit.ml
 qp_edit.native
-Input_auto_generated.ml
+qp_print
 qp_print.native
 qp_run
 qp_run.native
 qp_set_ddci
 qp_set_ddci.native
 qp_set_mo_class
 qp_set_mo_class.native
 qptypes_generator.byte
 Qptypes.ml
 test_atom
 test_atom.byte
 test_basis
 test_basis.byte
 test_bitlist
 test_bitlist.byte
 test_determinants
 test_determinants.byte
 test_elements
 test_elements.byte
 test_excitation
 test_excitation.byte
 test_gto
 test_gto.byte
 test_mo_label
 test_mo_label.byte
 test_molecule
 test_molecule.byte
 test_point3d
 test_point3d.byte
--- a/ocaml/Makefile
+++ b/ocaml/Makefile
@ -30,7 +30,7 @@ default: $(ALL_TESTS) $(ALL_EXE) .gitignore
 	qp_edit.ml qp_edit qp_edit.native Input_auto_generated.ml;\
 	do \
 	   echo $$i ; \
-	done > .gitignore
+	done |sort | uniq > .gitignore
 executables: $(QP_ROOT)/data/executables
--- a/plugins/CAS_SD/README.rst
+++ b/plugins/CAS_SD/README.rst
@ -118,69 +118,101 @@ Documentation
  Undocumented
-`h_apply_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L414>`_
+`h_apply_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L537>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_cas_sd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1>`_
+`h_apply_cas_sd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L3>`_
  Undocumented
 `h_apply_cas_sd_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L119>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L269>`_
+`h_apply_cas_sd_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L96>`_
  Undocumented
 `h_apply_cas_sd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L387>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2610>`_
+`h_apply_cas_sd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L3175>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_cas_sd_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2118>`_
+`h_apply_cas_sd_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2540>`_
  Undocumented
 `h_apply_cas_sd_pt2_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2678>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2427>`_
+`h_apply_cas_sd_pt2_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2644>`_
  Undocumented
 `h_apply_cas_sd_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2987>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1872>`_
+`h_apply_cas_sd_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2291>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_cas_sd_selected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1346>`_
+`h_apply_cas_sd_selected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1620>`_
  Undocumented
 `h_apply_cas_sd_selected_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1760>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_selected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1675>`_
+`h_apply_cas_sd_selected_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1725>`_
  Undocumented
 `h_apply_cas_sd_selected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2089>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_selected_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1128>`_
+`h_apply_cas_sd_selected_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1399>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_cas_sd_selected_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L602>`_
+`h_apply_cas_sd_selected_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L728>`_
  Undocumented
 `h_apply_cas_sd_selected_no_skip_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L868>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_cas_sd_selected_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L931>`_
+`h_apply_cas_sd_selected_no_skip_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L833>`_
  Undocumented
 `h_apply_cas_sd_selected_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1197>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
--- a/plugins/Full_CI/README.rst
+++ b/plugins/Full_CI/README.rst
@ -31,137 +31,201 @@ Documentation
  Undocumented
-`h_apply_fci <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L527>`_
+`h_apply_fci <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L674>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_fci_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1>`_
+`h_apply_fci_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3>`_
  Undocumented
 `h_apply_fci_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L143>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_mono <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2744>`_
+`h_apply_fci_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L108>`_
  Undocumented
 `h_apply_fci_mono <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3333>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_fci_mono_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2216>`_
+`h_apply_fci_mono_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2660>`_
  Undocumented
 `h_apply_fci_mono_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2800>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_mono_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2545>`_
+`h_apply_fci_mono_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2765>`_
  Undocumented
 `h_apply_fci_mono_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3129>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L330>`_
+`h_apply_fci_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L472>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1998>`_
+`h_apply_fci_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2439>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_fci_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1472>`_
+`h_apply_fci_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1768>`_
  Undocumented
 `h_apply_fci_no_skip_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1908>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1801>`_
+`h_apply_fci_no_skip_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1873>`_
  Undocumented
 `h_apply_fci_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2237>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1265>`_
+`h_apply_fci_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1558>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_fci_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L773>`_
+`h_apply_fci_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L923>`_
  Undocumented
 `h_apply_fci_pt2_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1061>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_fci_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1082>`_
+`h_apply_fci_pt2_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1027>`_
  Undocumented
 `h_apply_fci_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1370>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_pt2_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4258>`_
+`h_apply_pt2_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5141>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_pt2_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3764>`_
+`h_apply_pt2_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4504>`_
  Undocumented
 `h_apply_pt2_mono_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4642>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_pt2_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4073>`_
+`h_apply_pt2_mono_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4608>`_
  Undocumented
 `h_apply_pt2_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4951>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_pt2_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5729>`_
+`h_apply_pt2_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6906>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_pt2_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5237>`_
+`h_apply_pt2_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6271>`_
  Undocumented
 `h_apply_pt2_mono_di_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6409>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_pt2_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5546>`_
+`h_apply_pt2_mono_di_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6375>`_
  Undocumented
 `h_apply_pt2_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6718>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_select_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3518>`_
+`h_apply_select_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4255>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_select_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2990>`_
+`h_apply_select_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3582>`_
  Undocumented
 `h_apply_select_mono_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3722>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_select_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3319>`_
+`h_apply_select_mono_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3687>`_
  Undocumented
 `h_apply_select_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4051>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_select_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4991>`_
+`h_apply_select_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6022>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_select_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4465>`_
+`h_apply_select_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5351>`_
  Undocumented
 `h_apply_select_mono_di_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5491>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_select_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4794>`_
+`h_apply_select_mono_di_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5456>`_
  Undocumented
 `h_apply_select_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5820>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
--- a/plugins/Full_CI/full_ci.irp.f
+++ b/plugins/Full_CI/full_ci.irp.f
@ -11,6 +11,7 @@ program full_ci
  pt2 = 1.d0
  diag_algorithm = "Lapack"
  if (N_det > N_det_max) then
    call diagonalize_CI
    call save_wavefunction
--- a/plugins/MRCC_Utils/H_apply.irp.f
+++ b/plugins/MRCC_Utils/H_apply.irp.f
@ -9,7 +9,7 @@ s.data["declarations"] += """
    double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
    double precision, intent(in) :: delta_ii_(Ndet_ref,*)
 """
-s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc)"
+s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
 s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
 s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
 s.data["decls_main"] += """
--- a/plugins/MRCC_Utils/README.rst
+++ b/plugins/MRCC_Utils/README.rst
@ -1,171 +1,3 @@
 ===========
 MRCC Module
 ===========
 Needed Modules
 ==============
 .. Do not edit this section. It was auto-generated from the
 .. by the `update_README.py` script.
 .. image:: tree_dependency.png
 * `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
 * `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
 Documentation
 =============
 .. Do not edit this section. It was auto-generated from the
 .. by the `update_README.py` script.
 `ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L166>`_
  Eigenvectors/values of the CI matrix
 `ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L167>`_
  Eigenvectors/values of the CI matrix
 `ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L165>`_
  Eigenvectors/values of the CI matrix
 `ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L232>`_
  N_states lowest eigenvalues of the dressed CI matrix
 `davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L56>`_
  Davidson diagonalization with specific diagonal elements of the H matrix
  .br
  H_jj : specific diagonal H matrix elements to diagonalize de Davidson
  .br
  dets_in : bitmasks corresponding to determinants
  .br
  u_in : guess coefficients on the various states. Overwritten
  on exit
  .br
  dim_in : leftmost dimension of u_in
  .br
  sze : Number of determinants
  .br
  N_st : Number of eigenstates
  .br
  iunit : Unit for the I/O
  .br
  Initial guess vectors are not necessarily orthonormal
 `davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L1>`_
  Davidson diagonalization.
  .br
  dets_in : bitmasks corresponding to determinants
  .br
  u_in : guess coefficients on the various states. Overwritten
  on exit
  .br
  dim_in : leftmost dimension of u_in
  .br
  sze : Number of determinants
  .br
  N_st : Number of eigenstates
  .br
  iunit : Unit number for the I/O
  .br
  Initial guess vectors are not necessarily orthonormal
 `delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L104>`_
  Dressing matrix in N_det basis
 `delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L103>`_
  Dressing matrix in N_det basis
 `diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L247>`_
  Replace the coefficients of the CI states by the coefficients of the
  eigenstates of the CI matrix
 `find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L206>`_
  Undocumented
 `h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L422>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 `h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L1>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
 `h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L273>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
 `h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L140>`_
  Dressed H with Delta_ij
 `h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L360>`_
  Computes v_0 = H|u_0>
  .br
  n : number of determinants
  .br
  H_jj : array of <j|H|j>
 `lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L5>`_
  cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 `lambda_mrcc_tmp <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L81>`_
  Undocumented
 `lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L6>`_
  cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 `mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L15>`_
  Undocumented
 `mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L160>`_
  Undocumented
 `mrcc_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L7>`_
  Undocumented
 `oscillations <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L86>`_
  Undocumented
 `pert_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L1>`_
  Undocumented
 `psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L3>`_
  Locks on ref determinants to fill delta_ij
 `run_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L1>`_
  Undocumented
 `set_generators_bitmasks_as_holes_and_particles <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L69>`_
  Undocumented
 Needed Modules
 ==============
 .. Do not edit this section It was auto-generated
@ -258,7 +90,11 @@ Documentation
  N_states lowest eigenvalues of the dressed CI matrix
-`davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L56>`_
+`create_minilist <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L17>`_
  Undocumented
 `davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L59>`_
  Davidson diagonalization with specific diagonal elements of the H matrix
  .br
  H_jj : specific diagonal H matrix elements to diagonalize de Davidson
@ -279,7 +115,7 @@ Documentation
  Initial guess vectors are not necessarily orthonormal
-`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L1>`_
+`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L4>`_
  Davidson diagonalization.
  .br
  dets_in : bitmasks corresponding to determinants
@ -370,7 +206,7 @@ Documentation
  Find A.C = B
-`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L206>`_
+`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L266>`_
  Undocumented
@ -396,6 +232,22 @@ Documentation
  Undocumented
 `gen_det_idx <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L258>`_
  Undocumented
 `gen_det_shortcut <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L256>`_
  Undocumented
 `gen_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L255>`_
  Undocumented
 `gen_det_version <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L257>`_
  Undocumented
 `get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/LinearAlgebra.irp.f#L95>`_
  Find C = A^-1
@ -431,18 +283,26 @@ Documentation
  Undocumented
-`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L422>`_
+`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L553>`_
  Calls H_apply on the HF determinant and selects all connected single and double
  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L1>`_
+`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L3>`_
  Undocumented
 `h_apply_mrcc_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L127>`_
  Generate all double excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
-`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L273>`_
+`h_apply_mrcc_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L100>`_
  Undocumented
 `h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L399>`_
  Generate all single excitations of key_in using the bit masks of holes and
  particles.
  Assume N_int is already provided.
@ -452,7 +312,7 @@ Documentation
  Dressed H with Delta_ij
-`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L360>`_
+`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L374>`_
  Computes v_0 = H|u_0>
  .br
  n : number of determinants
@ -757,11 +617,7 @@ Documentation
  n!
-`mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L15>`_
+`mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L65>`_
  Undocumented
 `mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L160>`_
  Undocumented
@ -855,7 +711,7 @@ Documentation
  Current status for displaying progress bars. Global variable.
-`psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L3>`_
+`psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L4>`_
  Locks on ref determinants to fill delta_ij
--- a/plugins/MRCC_Utils/davidson.irp.f
+++ b/plugins/MRCC_Utils/davidson.irp.f
@ -1,3 +1,6 @@
 subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,istate)
  use bitmasks
  implicit none
@ -102,7 +105,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
  double precision               :: to_print(2,N_st)
  double precision               :: cpu, wall
-  PROVIDE det_connections
+  integer(bit_kind)              :: dets_in_sorted(Nint,2,sze)
  integer                        :: idx(sze), shortcut(0:sze+1),sh,ii,tmp
  !PROVIDE det_connections
  call write_time(iunit)
  call wall_time(wall)
@ -147,6 +153,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
  ! Initialization
  ! ==============
  dets_in_sorted(:,:,:) = dets_in(:,:,:)
  call sort_dets_ab(dets_in_sorted, idx, shortcut, sze, Nint)
  k_pairs=0
  do l=1,N_st
    do k=1,l
@ -205,6 +215,8 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
      ! ----------------------
      do k=1,N_st
        !call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in_sorted,shortcut,idx,Nint,istate)
        !call H_u_0_mrcc_org(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate)
        call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate)
      enddo
@ -357,6 +369,8 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
  abort_here = abort_all
 end
 subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
  use bitmasks
  implicit none
@ -377,37 +391,104 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
  double precision, allocatable  :: vt(:)
  integer                        :: i,j,k,l, jj,ii
  integer                        :: i0, j0
  integer                        :: shortcut(0:n+1), sort_idx(n)
  integer(bit_kind)              :: sorted(Nint,n), version(Nint,n)
  integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
 !   
  ASSERT (Nint > 0)
  ASSERT (Nint == N_int)
  ASSERT (n>0)
-  PROVIDE ref_bitmask_energy delta_ij 
+  PROVIDE ref_bitmask_energy delta_ij
  integer, parameter             :: block_size = 157
-  !$OMP PARALLEL DEFAULT(NONE)                                       &
+  
-      !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt)                             &
+ 
-      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
+   !$OMP PARALLEL DEFAULT(NONE)                                       &
-  !$OMP DO SCHEDULE(static)
+      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)                             &
      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij,sorted,shortcut,sort_idx,version)
 !$OMP DO SCHEDULE(static)
  do i=1,n  
    v_0(i) = H_jj(i) * u_0(i)
  enddo 
  !$OMP END DO
  allocate(idx(0:n), vt(n))
  Vt = 0.d0
-  !$OMP DO SCHEDULE(guided)
+  
-  do i=1,n
+  
-    idx(0) = i
+  !$OMP SINGLE
-    call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
+  call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
-    do jj=1,idx(0)
+  !$OMP END SINGLE
-      j = idx(jj)
+  
-!     if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
+
-        call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
+  !$OMP DO SCHEDULE(dynamic)
-        hij = hij 
+  do sh=1,shortcut(0)
-        vt (i) = vt (i) + hij*u_0(j)
+  do sh2=1,sh
-        vt (j) = vt (j) + hij*u_0(i)
+    exa = 0
-!     endif
+    do ni=1,Nint
-    enddo
+      exa += popcnt(xor(version(ni,sh), version(ni,sh2)))
    end do
    if(exa > 2) then
      cycle
    end if
    do i=shortcut(sh),shortcut(sh+1)-1
      if(sh==sh2) then
        endi = i-1
      else
        endi = shortcut(sh2+1)-1
      end if
      do j=shortcut(sh2),endi
        ext = exa
        do ni=1,Nint
          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
        end do
        if(ext <= 4) then
          org_i = sort_idx(i)
          org_j = sort_idx(j)
          call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
          vt (org_i) = vt (org_i) + hij*u_0(org_j)
          vt (org_j) = vt (org_j) + hij*u_0(org_i)
        end if
      end do
    end do
  end do
  enddo
  !$OMP END DO
-
+  
  !$OMP SINGLE
  call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
  !$OMP END SINGLE
  !$OMP DO SCHEDULE(dynamic)
  do sh=1,shortcut(0)
    do i=shortcut(sh),shortcut(sh+1)-1
      do j=shortcut(sh),i-1
        ext = 0
        do ni=1,Nint
          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
        end do
        if(ext == 4) then
          org_i = sort_idx(i)
          org_j = sort_idx(j)
          call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
          vt (org_i) = vt (org_i) + hij*u_0(org_j)
          vt (org_j) = vt (org_j) + hij*u_0(org_i)
        end if
      end do
    end do
  enddo
  !$OMP END DO
  !$OMP DO SCHEDULE(guided)
  do ii=1,n_det_ref
    i = idx_ref(ii)
@ -427,4 +508,3 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
  !$OMP END PARALLEL
 end
--- a/plugins/MRCC_Utils/mrcc_dress.irp.f
+++ b/plugins/MRCC_Utils/mrcc_dress.irp.f
@ -1,4 +1,5 @@
 use omp_lib
 use bitmasks
 BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
 implicit none
@ -12,12 +13,61 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
 END_PROVIDER
-subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
+
 subroutine create_minilist(key_mask, fullList, miniList, idx_miniList, N_fullList, N_miniList, Nint)
  use bitmasks
  implicit none
  integer(bit_kind), intent(in)            :: fullList(Nint, 2, N_fullList)
  integer, intent(in)                      :: N_fullList
  integer(bit_kind),intent(out)            :: miniList(Nint, 2, N_fullList)
  integer,intent(out)                      :: idx_miniList(N_fullList), N_miniList
  integer, intent(in)                      :: Nint
  integer(bit_kind)                        :: key_mask(Nint, 2)
  integer                                  :: ni, i, n_a, n_b, e_a, e_b
  n_a = 0
  n_b = 0
  do ni=1,nint
    n_a = n_a + popcnt(key_mask(ni,1))
    n_b = n_b + popcnt(key_mask(ni,2))
  end do
  if(n_a == 0) then
    N_miniList = N_fullList
    miniList(:,:,:) = fullList(:,:,:)
    do i=1,N_fullList
      idx_miniList(i) = i
    end do
    return
  end if
  N_miniList = 0
  do i=1,N_fullList
    e_a = n_a
    e_b = n_b
    do ni=1,nint
      e_a -= popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1)))
      e_b -= popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2)))
    end do
    if(e_a + e_b <= 2) then
      N_miniList = N_miniList + 1
      miniList(:,:,N_miniList) = fullList(:,:,i)
      idx_miniList(N_miniList) = i
    end if
  end do
 end subroutine
 subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
 use bitmasks
 implicit none
  integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet_ref, Ndet_non_ref
+  integer, intent(in)            :: Ndet_ref, Ndet_non_ref
  double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
  double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
@ -40,17 +90,65 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
  integer(bit_kind)              :: tmp_det(Nint,2)
  integer                        :: iint, ipos
  integer                        :: i_state, k_sd, l_sd, i_I, i_alpha
-
+  
-  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
+  integer(bit_kind),allocatable  :: miniList(:,:,:)
  integer(bit_kind),intent(in)   :: key_mask(Nint, 2)
  integer,allocatable            :: idx_miniList(:)
  integer                        :: N_miniList, ni
  allocate(miniList(Nint, 2, max(N_det_generators, N_det_non_ref)), idx_miniList(max(N_det_generators, N_det_non_ref)))
  l = 0
  do ni = 1,Nint
    l += popcnt(key_mask(ni,1)) + popcnt(key_mask(ni,2))
  end do
  if(l == 0) then
    N_miniList = i_generator-1
    miniList(:,:,:N_miniList) = psi_det_generators(:,:,:N_minilist)
  else
    N_miniList = 0
    do i=i_generator-1,1,-1
      k = l
      do ni=1,nint
        k -= popcnt(iand(key_mask(ni,1), psi_det_generators(ni,1,i))) + popcnt(iand(key_mask(ni,2), psi_det_generators(ni,2,i)))
      end do
      if(k == 0) then
        deallocate(miniList, idx_miniList)
        return
      end if
      if(k <= 2) then
        N_minilist += 1
        miniList(:,:,N_minilist) = psi_det_generators(:,:,i)
      end if
    end do                                                                                                                           
  end if
  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
  allocate (dIa_hla(N_states,Ndet_non_ref))
  ! |I>
  ! |alpha>
   if(N_tq > 0) then
     call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
   end if
  do i_alpha=1,N_tq
-    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
+!    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
-
+    call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
    do j=1,idx_alpha(0)
      idx_alpha(j) = idx_miniList(idx_alpha(j))
    end do
    ! |I>
    do i_I=1,N_det_ref
       ! Find triples and quadruple grand parents
@ -149,61 +247,24 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
    enddo
  enddo
  deallocate (dIa_hla)
  deallocate(miniList, idx_miniList)
 end
 BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted,  (N_int,2,N_det_generators,2) ]
 &BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
 &BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
 &BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
  gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
  gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
  call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
  call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
 END_PROVIDER
 subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
 subroutine mrcc_dress_simple(delta_ij_non_ref_,Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
 use bitmasks
 implicit none
  integer, intent(in)            :: i_generator,n_selected, Nint, iproc
  integer, intent(in) :: Ndet_non_ref
  double precision, intent(inout) :: delta_ij_non_ref_(Ndet_non_ref,Ndet_non_ref,*)
  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
  integer                        :: i,j,k,m
  integer                        :: new_size
  integer                        :: degree(psi_det_size)
  integer                        :: idx(0:psi_det_size)
  logical                        :: good
  integer(bit_kind)              :: tq(Nint,2,n_selected)
  integer                        :: N_tq, c_ref
  integer                        :: connected_to_ref
  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
  ! Compute <k|H|a><a|H|j> / (E0 - Haa)
  double precision :: hka, haa
  double precision :: haj
  double precision :: f(N_states)
  do i=1,N_tq
    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i),degree,Nint,Ndet_non_ref,idx)
    call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
    do m=1,N_states
      f(m) = 1.d0/(ci_electronic_energy(m)-haa)
    enddo
    do k=1,idx(0)
      call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(k)),Nint,hka)
      do j=k,idx(0)
        call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(j)),Nint,haj)
        do m=1,N_states
          delta_ij_non_ref_(idx(k), idx(j),m) += haj*hka* f(m)
          delta_ij_non_ref_(idx(j), idx(k),m) += haj*hka* f(m)
        enddo
      enddo 
    enddo
  enddo
 end
 subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
 use bitmasks
 implicit none
@ -218,22 +279,37 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
  integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
  integer, intent(out)           :: N_tq
-  integer                        :: c_ref
+  
-  integer                        :: connected_to_ref
+  
-
+  integer                        :: nt,ni
  integer(bit_kind),intent(in)  :: miniList(Nint,2,N_det_generators)
  integer,intent(in)            :: N_miniList
  N_tq = 0
-  do i=1,N_selected
+  i_loop : do i=1,N_selected
-    c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, &
+    do j=1,N_miniList
-       i_generator,N_det_generators)
+      nt = 0
-
+      do ni=1,Nint
-    if (c_ref /= 0) then
+        nt += popcnt(xor(miniList(ni,1,j), det_buffer(ni,1,i))) + popcnt(xor(miniList(ni,2,j), det_buffer(ni,2,i)))
-      cycle
+      end do
-    endif
+      if(nt <= 4) then
        cycle i_loop
      end if
    end do
 !     if(connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, &
 !        i_generator,N_det_generators) /= 0) then
 !         cycle i_loop
 !     end if
    ! Select determinants that are triple or quadruple excitations
    ! from the ref
    good = .True.
-    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) 
    !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
    do k=1,idx(0)
      if (degree(k) < 3) then
        good = .False.
@ -249,8 +325,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
        enddo
      endif
    endif
-  enddo
+  enddo i_loop
 end
--- a/plugins/Perturbation/README.rst
+++ b/plugins/Perturbation/README.rst
@ -112,42 +112,74 @@ Documentation
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L686>`_
 =======
 `perturb_buffer_by_mono_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L266>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L476>`_
 =======
 `perturb_buffer_by_mono_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L371>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L581>`_
 =======
 `perturb_buffer_by_mono_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L476>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L371>`_
 =======
 `perturb_buffer_by_mono_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L791>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L266>`_
 =======
 `perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L686>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L161>`_
 =======
 `perturb_buffer_by_mono_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L581>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L56>`_
 =======
 `perturb_buffer_by_mono_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L896>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_by_mono_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L896>`_
 =======
 `perturb_buffer_by_mono_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L56>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
  routine.
@ -157,42 +189,74 @@ Documentation
  routine.
 <<<<<<< HEAD
 `perturb_buffer_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L635>`_
 =======
 `perturb_buffer_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L215>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L425>`_
 =======
 `perturb_buffer_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L320>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L530>`_
 =======
 `perturb_buffer_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L425>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L320>`_
 =======
 `perturb_buffer_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L740>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L215>`_
 =======
 `perturb_buffer_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L635>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L110>`_
 =======
 `perturb_buffer_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L530>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L5>`_
 =======
 `perturb_buffer_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L845>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
  routine.
 <<<<<<< HEAD
 `perturb_buffer_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L845>`_
 =======
 `perturb_buffer_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L5>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
  routine.
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@ -99,7 +99,7 @@ class H_apply(object):
  deallocate(H_jj,iorder)
    """
-    s["size_max"] = str(1024*128) 
+    s["size_max"] = "256"
    s["copy_buffer"] = """call copy_H_apply_buffer_to_wf
  if (s2_eig) then
    call make_s2_eigenfunction
@ -266,7 +266,7 @@ class H_apply(object):
      double precision, intent(inout) :: select_max_out"""
      self.data["params_post"] += ", select_max(min(i_generator,size(select_max,1)))"
-      self.data["size_max"] = str(1024*128) 
+      self.data["size_max"] = "256"
      self.data["copy_buffer"] = """
      call copy_H_apply_buffer_to_wf
      if (s2_eig) then
--- a/src/Determinants/H_apply.template.f
+++ b/src/Determinants/H_apply.template.f
@ -1,4 +1,122 @@
-subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters )
+
 subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters )
  integer(bit_kind), intent(in)         :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2)
  integer(bit_kind), intent(in)         :: particl_1(N_int, 2), particl_2(N_int, 2)
  integer(bit_kind)                     :: p1_mask(N_int, 2), p2_mask(N_int, 2), tmp
  integer,intent(in)                    :: i_generator,iproc_in
  integer(bit_kind)                     :: status(N_int*bit_kind_size, 2)
  integer                               :: highest, p1,p2,sp,ni,i,mi,nt,ns
  integer(bit_kind), intent(in)          :: key_prev(N_int, 2, *)
  $declarations
  highest = 0
  status(:,:) = 0
  do sp=1,2
    do ni=1,N_int
      do i=1,bit_kind_size
        if(iand(1,ishft(key_in(ni, sp), -(i-1))) == 0) then
          cycle
        end if
        mi = (ni-1)*bit_kind_size+i
        status(mi, sp) = iand(1,ishft(hole_1(ni, sp), -(i-1)))
        status(mi, sp) = status(mi, sp) + 2*iand(1,ishft(hole_2(ni, sp), -(i-1)))
        if(status(mi, sp) /= 0 .and. mi > highest) then
          highest = mi
        end if
      end do
    end do
  end do
 ! ! GEL D'ELECTRONS
 ! !  nt = 0
 !   do i=1,i_generator-1
 !     if(key_in(1,1) == key_prev(1,1,i)) then
 !       tmp = xor(key_in(1,2), key_prev(1,2,i))
 !       if(popcnt(tmp) == 2) then
 !         ns = 1+trailz(iand(tmp, key_in(1,2)))
 ! !         if(status(ns, 2) /= 0) then
 ! !           nt += 1
 ! !         end if
 !         status(ns, 2) = 0
 !       end if
 !     else if(key_in(1,2) == key_prev(1,2,i)) then
 !       tmp = xor(key_in(1,1), key_prev(1,1,i))
 !       if(popcnt(tmp) == 2) then
 !         ns = 1+trailz(iand(tmp, key_in(1,1)))
 ! !         if(status(ns, 1) /= 0) then
 ! !           nt += 1
 ! !         end if
 !         status(ns, 1) = 0
 !       end if
 !     end if
 !   end do
 ! !  print *, "nt", nt, i_generator
  do sp=1,2
    do p1=1,highest
      if(status(p1, sp) == 0) then
        cycle
      end if
      do p2=1,highest
        if(status(p2, sp) == 0) then
          cycle
        end if
        if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. &
            (status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. &
            (status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then
          call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters )
        end if
      end do
    end do
  end do
  do p1=1,highest
    if(status(p1, 1) == 0) then
      cycle
    end if
    do p2=1,highest
      if(status(p2, 2) == 0) then
        cycle
      end if
      if((status(p1, 1) == 3) .or. &
          (status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. &
          (status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then
          call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters )
      end if
    end do
  end do
 end subroutine
 subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters )
  integer(bit_kind), intent(in)         :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2)
  integer(bit_kind)                     :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2)
  integer,intent(in)                    :: fh1,fh2,fs1,fs2,i_generator,iproc_in
  integer(bit_kind)                     :: miniList(N_int, 2, N_det)
  integer                               :: n_minilist, n_alpha, n_beta, deg(2), i, ni
  $declarations
  p1_mask(:,:) = 0_bit_kind
  p2_mask(:,:) = 0_bit_kind
  p1_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) = ishft(1,iand(fh1-1,bit_kind_size-1))
  p2_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) = ishft(1,iand(fh2-1,bit_kind_size-1))
  key_mask(:,:) = key_in(:,:)
  key_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) -= ishft(1,iand(fh1-1,bit_kind_size-1))
  key_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) -= ishft(1,iand(fh2-1,bit_kind_size-1))
  call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters )
 end subroutine
 subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters )
  use omp_lib
  use bitmasks
  implicit none
@ -10,7 +128,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
  integer,parameter              :: size_max = $size_max
  $declarations
  integer          ,intent(in)   :: i_generator
-  integer(bit_kind),intent(in)   :: key_in(N_int,2)
+  integer(bit_kind),intent(in)   :: key_in(N_int,2), key_mask(N_int, 2)
  integer(bit_kind),allocatable  :: keys_out(:,:,:)
  integer(bit_kind), intent(in)  :: hole_1(N_int,2), particl_1(N_int,2)
  integer(bit_kind), intent(in)  :: hole_2(N_int,2), particl_2(N_int,2)
@ -290,13 +408,18 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
  double precision               :: diag_H_mat_elem
  integer(omp_lock_kind), save   :: lck, ifirst=0
  integer                        :: iproc
-
+  
  integer(bit_kind)              :: key_mask(N_int, 2)
  logical :: check_double_excitation 
  key_mask(:,:) = 0_bit_kind
  iproc = iproc_in
  check_double_excitation = .True.
  $check_double_excitation
-
+  
  if (ifirst == 0) then
    ifirst=1
@ -412,7 +535,6 @@ subroutine $subroutine($params_main)
  nmax = mod( N_det_generators,nproc )
  !$ call omp_init_lock(lck)
  call start_progress(N_det_generators,'Selection (norm)',0.d0)
@ -455,6 +577,7 @@ subroutine $subroutine($params_main)
    enddo
    if($do_double_excitations)then
     call $subroutine_diexc(psi_det_generators(1,1,i_generator),      &
         psi_det_generators(1,1,1),                                   &
         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
         i_generator, iproc $params_post)
@ -515,6 +638,7 @@ subroutine $subroutine($params_main)
    if($do_double_excitations)then
      call $subroutine_diexc(psi_det_generators(1,1,i_generator),    &
        psi_det_generators(1,1,1),                                   &
        mask(1,1,d_hole1), mask(1,1,d_part1),                        &
        mask(1,1,d_hole2), mask(1,1,d_part2),                        &
        i_generator, iproc $params_post)
--- a/src/Determinants/README.rst
+++ b/src/Determinants/README.rst
@ -1,42 +1,3 @@
 ===========
 Dets Module
 ===========
 This module contains the determinants of the CI wave function.
 H is applied on the list of generator determinants. Selected determinants
 are added into the *H_apply buffer*. Then the new wave function is
 constructred as the concatenation of the odl wave function and
 some determinants of the H_apply buffer. Generator determinants are built
 as a subset of the determinants of the wave function.
 Assumptions
 ===========
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES_CHILDREN file by the `update_README.py` script.
 * The MOs are orthonormal
 * All the determinants have the same number of electrons
 * The determinants are orthonormal
 * The number of generator determinants <= the number of determinants
 * All the determinants in the H_apply buffer are supposed to be different from the 
  wave function determinants
 * All the determinants in the H_apply buffer are supposed to be unique
 Needed Modules
 ==============
 .. Do not edit this section It was auto-generated
 .. by the `update_README.py` script.
 .. image:: tree_dependency.png
 * `Integrals_Monoelec <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec>`_
 * `Integrals_Bielec <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Bielec>`_
 Needed Modules
 ==============
 .. Do not edit this section It was auto-generated
@ -54,7 +15,7 @@ Documentation
 .. by the `update_README.py` script.
-`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1108>`_
+`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1110>`_
  Needed for diag_H_mat_elem
@ -66,7 +27,7 @@ Documentation
  Max and min values of the coefficients
-`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1153>`_
+`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1155>`_
  Needed for diag_H_mat_elem
@ -157,11 +118,11 @@ Documentation
  of alpha and beta determinants
-`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L382>`_
+`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L565>`_
  True if the Davidson algorithm is converged
-`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L372>`_
+`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L555>`_
  Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
@ -184,7 +145,7 @@ Documentation
  Initial guess vectors are not necessarily orthonormal
-`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L68>`_
+`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L241>`_
  Davidson diagonalization with specific diagonal elements of the H matrix
  .br
  H_jj : specific diagonal H matrix elements to diagonalize de Davidson
@ -213,7 +174,7 @@ Documentation
  Max number of Davidson sizes
-`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L373>`_
+`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L556>`_
  Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
@ -241,8 +202,8 @@ Documentation
  det_coef
-`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1283>`_
+`det_inf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L69>`_
-  Build connection proxy between determinants
+  Undocumented
 `det_occ <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L226>`_
@ -261,7 +222,7 @@ Documentation
  Diagonalization algorithm (Davidson or Lapack)
-`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1046>`_
+`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1048>`_
  Computes <i|H|i>
@ -321,7 +282,7 @@ Documentation
  idx(0) is the number of determinants that interact with key1
-`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L163>`_
+`filter_connected_davidson_warp <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L102>`_
  Filters out the determinants that are not connected by H
  returns the array idx which contains the index of the
  determinants in the array key1 that interact
@ -331,7 +292,7 @@ Documentation
  key1 should come from psi_det_sorted_ab.
-`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L293>`_
+`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L390>`_
  returns the array idx which contains the index of the
  .br
  determinants in the array key1 that interact
@ -341,7 +302,7 @@ Documentation
  idx(0) is the number of determinants that interact with key1
-`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L392>`_
+`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L489>`_
  standard filter_connected_i_H_psi but returns in addition
  .br
  the array of the index of the non connected determinants to key1
@ -353,6 +314,7 @@ Documentation
  to repeat the excitations
 <<<<<<< HEAD
 `filter_connected_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L101>`_
  Filters out the determinants that are not connected by H
  returns the array idx which contains the index of the
@ -364,10 +326,13 @@ Documentation
 `generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L500>`_
 =======
 `generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L472>`_
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
  Create a wave function from all possible alpha x beta determinants
-`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L141>`_
+`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L142>`_
  Returns the two excitation operators between two doubly excited determinants and the phase
@ -379,7 +344,7 @@ Documentation
  Returns the excitation degree between two determinants
-`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L951>`_
+`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L953>`_
  Applies get_excitation_degree to an array of determinants
@ -395,11 +360,11 @@ Documentation
  Returns the index of the determinant in the ``psi_det_sorted_bit`` array
-`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L274>`_
+`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L275>`_
  Returns the excitation operator between two singly excited determinants and the phase
-`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1201>`_
+`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1203>`_
  Returns a list of occupation numbers from a bitstring
@ -433,7 +398,7 @@ Documentation
  Undocumented
-`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1217>`_
+`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1219>`_
  Computes v_0 = H|u_0>
  .br
  n : number of determinants
@ -441,23 +406,23 @@ Documentation
  H_jj : array of <j|H|j>
-`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L355>`_
+`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L356>`_
  Returns <i|H|j> where i and j are determinants
-`i_h_j_phase_out <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L491>`_
+`i_h_j_phase_out <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L492>`_
  Returns <i|H|j> where i and j are determinants
-`i_h_j_verbose <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L627>`_
+`i_h_j_verbose <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L628>`_
  Returns <i|H|j> where i and j are determinants
-`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L766>`_
+`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L767>`_
  <key|H|psi> for the various Nstates
-`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L848>`_
+`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L850>`_
  <key|H|psi> for the various Nstate
  .br
  returns in addition
@ -471,7 +436,7 @@ Documentation
  to repeat the excitations
-`i_h_psi_sc2_verbose <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L895>`_
+`i_h_psi_sc2_verbose <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L897>`_
  <key|H|psi> for the various Nstate
  .br
  returns in addition
@ -485,7 +450,7 @@ Documentation
  to repeat the excitations
-`i_h_psi_sec_ord <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L801>`_
+`i_h_psi_sec_ord <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L803>`_
  <key|H|psi> for the various Nstates
@ -532,10 +497,6 @@ Documentation
  Energy of the reference bitmask used in Slater rules
 `n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1275>`_
  Number of integers to represent the connections between determinants
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L3>`_
  Number of determinants in the wave function
@ -897,6 +858,18 @@ Documentation
  for a given couple of hole/particle excitations i.
 `sort_dets_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L197>`_
  Undocumented
 `sort_dets_ab_v <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L145>`_
  Undocumented
 `sort_dets_ba_v <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L124>`_
  Undocumented
 `sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L489>`_
  Determinants on which we apply <i|H|j>.
  They are sorted by the 3 highest electrons in the alpha part,
@ -918,8 +891,13 @@ Documentation
  Weights in the state-average calculation of the density matrix
 <<<<<<< HEAD
 `target_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L160>`_
  Energy that should be obtained when truncating the wave function (optional)
 =======
 `tamiser <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L91>`_
  Undocumented
 >>>>>>> 2d3ba8003b05cb406674cc93a9bfc917d94db7fc
 `threshold_convergence_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L18>`_
--- a/src/Determinants/connected_to_ref.irp.f
+++ b/src/Determinants/connected_to_ref.irp.f
@ -169,7 +169,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
  ! output :   0 : not connected
  !            i : connected to determinant i of the past
-  !           -i : is the ith determinant of the refernce wf keys
+  !           -i : is the ith determinant of the reference wf keys
  ASSERT (Nint > 0)
  ASSERT (Nint == N_int)
--- a/src/Determinants/davidson.irp.f
+++ b/src/Determinants/davidson.irp.f
@ -65,6 +65,179 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
  deallocate (H_jj)
 end
 logical function det_inf(key1, key2, Nint)
  use bitmasks
  implicit none
  integer(bit_kind),intent(in)       :: key1(Nint, 2), key2(Nint, 2)
  integer,intent(in)                 :: Nint
  integer                            :: i,j
  det_inf = .false.
  do i=1,2
    do j=Nint,1,-1
      if(key1(j,i) < key2(j,i)) then
        det_inf = .true.
        return
      else if(key1(j,i) > key2(j,i)) then
        return
      end if
    end do
  end do
 end function
 subroutine tamiser(key, idx, no, n, Nint, N_key)
  use bitmasks
  implicit none
  integer(bit_kind),intent(inout)       :: key(Nint, 2, N_key)
  integer,intent(in)                    :: no, n, Nint, N_key
  integer,intent(inout)                 :: idx(N_key)
  integer                               :: k,j,tmpidx
  integer(bit_kind)                     :: tmp(Nint, 2)
  logical                               :: det_inf
  k = no
  j = 2*k
  do while(j <= n)
    if(j < n .and. det_inf(key(:,:,j), key(:,:,j+1), Nint)) then
      j = j+1
    end if
    if(det_inf(key(:,:,k), key(:,:,j), Nint)) then
      tmp(:,:) = key(:,:,k)
      key(:,:,k) = key(:,:,j)
      key(:,:,j) = tmp(:,:)
      tmpidx = idx(k)
      idx(k) = idx(j)
      idx(j) = tmpidx
      k = j
      j = 2*k
    else
      return
    end if
  end do
 end subroutine
 subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
  use bitmasks
  implicit none
  integer(bit_kind),intent(in)          :: key_in(Nint,2,N_key)
  integer(bit_kind)                     :: key(Nint,2,N_key)
  integer(bit_kind),intent(out)         :: key_out(Nint,N_key)
  integer,intent(out)                   :: idx(N_key)
  integer,intent(out)                   :: shortcut(0:N_key+1)
    integer(bit_kind),intent(out)       :: version(Nint,N_key+1)
  integer, intent(in)                   :: Nint, N_key
  integer(bit_kind)                     :: tmp(Nint, 2,N_key)
  key(:,1,:N_key) = key_in(:,2,:N_key)
  key(:,2,:N_key) = key_in(:,1,:N_key)
  call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint)
 end subroutine
 subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
  use bitmasks
  implicit none
  integer(bit_kind),intent(in)          :: key_in(Nint,2,N_key)
  integer(bit_kind)                     :: key(Nint,2,N_key)
  integer(bit_kind),intent(out)         :: key_out(Nint,N_key)
  integer,intent(out)                   :: idx(N_key)
  integer,intent(out)                   :: shortcut(0:N_key+1)
  integer(bit_kind),intent(out)         :: version(Nint,N_key+1)
  integer, intent(in)                   :: Nint, N_key
  integer(bit_kind)                     :: tmp(Nint, 2)
  integer                               :: tmpidx,i,ni
  key(:,:,:) = key_in(:,:,:)
  do i=1,N_key
    idx(i) = i
  end do
  do i=N_key/2,1,-1
    call tamiser(key, idx, i, N_key, Nint, N_key)
  end do
  do i=N_key,2,-1
    tmp(:,:) = key(:,:,i)
    key(:,:,i) = key(:,:,1)
    key(:,:,1) = tmp(:,:)
    tmpidx = idx(i)
    idx(i) = idx(1)
    idx(1) = tmpidx
    call tamiser(key, idx, 1, i-1, Nint, N_key)
  end do
  shortcut(0) = 1
  shortcut(1) = 1
  version(:,1) = key(:,1,1)
  do i=2,N_key
    do ni=1,nint
      if(key(ni,1,i) /= key(ni,1,i-1)) then
        shortcut(0) = shortcut(0) + 1
        shortcut(shortcut(0)) = i
        version(:,shortcut(0)) = key(:,1,i)
        exit
      end if
    end do
  end do
  shortcut(shortcut(0)+1) = N_key+1
  key_out(:,:) = key(:,2,:)
 end subroutine
 c
 subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
  use bitmasks
  implicit none
  integer(bit_kind),intent(inout)       :: key(Nint,2,N_key)
  integer,intent(out)                   :: idx(N_key)
  integer,intent(out)                   :: shortcut(0:N_key+1)
  integer, intent(in)                   :: Nint, N_key
  integer(bit_kind)                     :: tmp(Nint, 2)
  integer                               :: tmpidx,i,ni
  do i=1,N_key
    idx(i) = i
  end do
  do i=N_key/2,1,-1
    call tamiser(key, idx, i, N_key, Nint, N_key)
  end do
  do i=N_key,2,-1
    tmp(:,:) = key(:,:,i)
    key(:,:,i) = key(:,:,1)
    key(:,:,1) = tmp(:,:)
    tmpidx = idx(i)
    idx(i) = idx(1)
    idx(1) = tmpidx
    call tamiser(key, idx, 1, i-1, Nint, N_key)
  end do
  shortcut(0) = 1
  shortcut(1) = 1
  do i=2,N_key
    do ni=1,nint
      if(key(ni,1,i) /= key(ni,1,i-1)) then
        shortcut(0) = shortcut(0) + 1
        shortcut(shortcut(0)) = i
        exit
      end if
    end do
  end do
  shortcut(shortcut(0)+1) = N_key+1
 end subroutine
 subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit)
  use bitmasks
  implicit none
@ -106,7 +279,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  integer                        :: k_pairs, kl
  integer                        :: iter2
-  double precision, allocatable  :: W(:,:,:),  U(:,:,:), R(:,:)
+  double precision, allocatable  :: W(:,:,:),  U(:,:,:), R(:,:), Wt(:)
  double precision, allocatable  :: y(:,:,:,:), h(:,:,:,:), lambda(:)
  double precision               :: diag_h_mat_elem
  double precision               :: residual_norm(N_st)
@ -114,7 +287,10 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  double precision               :: to_print(2,N_st)
  double precision               :: cpu, wall
-  PROVIDE det_connections
+  integer(bit_kind)              :: dets_in_sorted(Nint, 2, sze)
  integer                        :: idx(sze), shortcut(0:sze+1)
  !PROVIDE det_connections
  call write_time(iunit)
  call wall_time(wall)
@ -145,6 +321,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  allocate(                                                          &
      kl_pairs(2,N_st*(N_st+1)/2),                                   &
      W(sze,N_st,davidson_sze_max),                                                   &
      Wt(sze),                                                        &
      U(sze,N_st,davidson_sze_max),                                  &
      R(sze,N_st),                                                   &
      h(N_st,davidson_sze_max,N_st,davidson_sze_max),                &
@ -159,6 +336,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  ! Initialization
  ! ==============
  dets_in_sorted(:,:,:) = dets_in(:,:,:)
  call sort_dets_ab(dets_in_sorted, idx, shortcut, sze, Nint)
  k_pairs=0
  do l=1,N_st
    do k=1,l
@ -170,7 +350,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  !$OMP PARALLEL DEFAULT(NONE)                                       &
      !$OMP  SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs,             &
-      !$OMP  Nint,dets_in,u_in)                                 &
+      !$OMP  Nint,dets_in_sorted,dets_in,u_in)                                 &
      !$OMP  PRIVATE(k,l,kl,i)
@ -217,9 +397,11 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
      ! ----------------------
      do k=1,N_st
-        call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
+!          call H_u_0_org(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
         call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in_sorted,shortcut,idx,Nint)
      enddo
      ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
      ! -------------------------------------------
@ -274,7 +456,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
        to_print(1,k) = lambda(k) + nuclear_repulsion
        to_print(2,k) = residual_norm(k)
      enddo
-
+      
      write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st)
      call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
      if (converged) then
@ -360,6 +542,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
  deallocate (                                                       &
      kl_pairs,                                                      &
      W,                                                             &
      Wt,                                                            &
      U,                                                             &
      R,                                                             &
      h,                                                             &
--- a/src/Determinants/filter_connected.irp.f
+++ b/src/Determinants/filter_connected.irp.f
@ -98,69 +98,8 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
 end
 subroutine filter_connected_sorted_ab(key1,key2,next,Nint,sze,idx)
  use bitmasks
  implicit none
  BEGIN_DOC
  ! Filters out the determinants that are not connected by H
  ! returns the array idx which contains the index of the 
  ! determinants in the array key1 that interact 
  ! via the H operator with key2.
  ! idx(0) is the number of determinants that interact with key1
  !
  ! Determinants are taken from the psi_det_sorted_ab array
  END_DOC
  integer, intent(in)            :: Nint, sze
  integer, intent(in)            :: next(2,N_det)
  integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
  integer(bit_kind), intent(in)  :: key2(Nint,2)
  integer, intent(out)           :: idx(0:sze)
  integer                        :: i,j,l
  integer                        :: degree_x2
  integer(bit_kind)              :: det3_1(Nint,2), det3_2(Nint,2)
  ASSERT (Nint > 0)
  ASSERT (sze >= 0)
-  l=1
+subroutine filter_connected_davidson_warp(key1,warp,key2,Nint,sze,idx)
  call filter_3_highest_electrons( key2(1,1), det3_2(1,1), Nint)
  if (Nint==1) then
    i = 1
    do while ( i<= sze )
      call filter_3_highest_electrons( key1(1,1,i), det3_1(1,1), Nint)
      degree_x2 = popcnt( xor( det3_1(1,1), det3_2(1,1))) 
      if (degree_x2 > 4) then
        i = next(1,i)
        cycle
      else 
        degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1)) ) 
        if (degree_x2 <= 4) then
          degree_x2 += popcnt( xor( key1(1,2,i), key2(1,2)) )
          if (degree_x2 <= 4) then
            idx(l) = i
            l += 1
          endif
        endif
        i += 1
      endif
    enddo
  else
    print *,  'Not implemented', irp_here
    stop 1
  endif
  idx(0) = l-1
 end
 subroutine filter_connected_davidson(key1,key2,Nint,sze,idx)
  use bitmasks
  implicit none
  BEGIN_DOC
@ -177,119 +116,277 @@ subroutine filter_connected_davidson(key1,key2,Nint,sze,idx)
  integer(bit_kind), intent(in)  :: key2(Nint,2)
  integer, intent(out)           :: idx(0:sze)
  integer,intent(in)             :: warp(2,0:sze+1)
  integer                        :: i,j,k,l
  integer                        :: degree_x2
-  integer :: j_int, j_start
+  integer :: i_alpha, i_beta, exc_a, exc_b, endloop, ni
-  integer*8 :: itmp
+  integer(bit_kind) :: tmp1, tmp2
  PROVIDE N_con_int det_connections
  ASSERT (Nint > 0)
  ASSERT (sze >= 0)
  l=1
  i_alpha = 0
-  if (Nint==1) then
+  
-
+  if (Nint /= 1) then
-    i = idx(0)
+    do while(i_alpha < warp(1,0) .and. warp(1,i_alpha+1) <= sze)
-    do j_int=1,N_con_int
+      i_alpha = i_alpha + 1
-      itmp = det_connections(j_int,i)
+      exc_a = 0
-      do while (itmp /= 0_8)
+      do ni=1,Nint
-        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+        exc_a += popcnt(xor(key1(ni,1,warp(1,i_alpha)), key2(ni,1)))
-        do j = j_start+1, min(j_start+32,i-1)
+      end do
-          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +             &
+      endloop = min(warp(2,i_alpha), sze)
-              popcnt(xor( key1(1,2,j), key2(1,2)))
+      if(exc_a == 4) then
-          if (degree_x2 > 4) then
+        beta_loop : do i_beta=warp(1,i_alpha),endloop
-            cycle
+          do ni=1,Nint
-          else
+            if(key1(ni,2,i_beta) /= key2(ni,2)) then
-            idx(l) = j
+              cycle beta_loop
-            l = l+1
+            end if
-          endif
+          end do
-        enddo
+          idx(l) = i_beta
-        itmp = iand(itmp-1_8,itmp)
+          l = l + 1
-      enddo
+          exit beta_loop
-    enddo
+        end do beta_loop
-    
+      else
-  else if (Nint==2) then
+        do i_beta=warp(1,i_alpha),endloop
-    
+          exc_b = 0
-    
+          do ni=1,Nint
-    i = idx(0)
+            exc_b += popcnt(xor(key1(ni,2,i_beta), key2(ni,2)))
-    do j_int=1,N_con_int
+          end do
-      itmp = det_connections(j_int,i)
+          if(exc_b + exc_a <= 4) then
-      do while (itmp /= 0_8)
+            idx(l) = i_beta
-        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+            l = l + 1
-        do j = j_start+1, min(j_start+32,i-1)
+          end if                          
-          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
+        end do
-              popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
+      end if
-              popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
+    end do
              popcnt(xor( key1(2,2,j), key2(2,2)))
          if (degree_x2 > 4) then
            cycle
          else
            idx(l) = j
            l = l+1
          endif
        enddo
        itmp = iand(itmp-1_8,itmp)
      enddo
    enddo
  else if (Nint==3) then
    i = idx(0)
    !DIR$ LOOP COUNT (1000)
    do j_int=1,N_con_int
      itmp = det_connections(j_int,i)
      do while (itmp /= 0_8)
        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
        do j = j_start+1, min(j_start+32,i-1)
          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
              popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
              popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
              popcnt(xor( key1(2,2,j), key2(2,2))) +                 &
              popcnt(xor( key1(3,1,j), key2(3,1))) +                 &
              popcnt(xor( key1(3,2,j), key2(3,2)))
          if (degree_x2 > 4) then
            cycle
          else
            idx(l) = j
            l = l+1
          endif
        enddo
        itmp = iand(itmp-1_8,itmp)
      enddo
    enddo
  else
    do while(i_alpha < warp(1,0) .and. warp(1,i_alpha+1) <= sze)
      i_alpha = i_alpha + 1
      exc_a = popcnt(xor(key1(1,1,warp(1,i_alpha)), key2(1,1)))
      endloop = min(warp(2,i_alpha), sze)
      if(exc_a == 4) then
        do i_beta=warp(1,i_alpha),endloop
          if(key1(1,2,i_beta) == key2(1,2)) then
            idx(l) = i_beta
            l = l + 1
            exit
          end if
        end do
      else
        do i_beta=warp(1,i_alpha),endloop
          exc_b = popcnt(xor(key1(1,2,i_beta), key2(1,2)))
          if(exc_b + exc_a <= 4) then
            idx(l) = i_beta
            l = l + 1
          end if                          
        end do
      end if
    end do
  end if
    i = idx(0)
    !DIR$ LOOP COUNT (1000)
    do j_int=1,N_con_int
      itmp = det_connections(j_int,i)
      do while (itmp /= 0_8)
        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
        do j = j_start+1, min(j_start+32,i-1)
          degree_x2 = 0
          !DEC$ LOOP COUNT MIN(4)
          do k=1,Nint
            degree_x2 = degree_x2+ popcnt(xor( key1(k,1,j), key2(k,1))) +&
                popcnt(xor( key1(k,2,j), key2(k,2)))
            if (degree_x2 > 4) then
              exit
            endif
          enddo
          if (degree_x2 <= 5) then
            idx(l) = j
            l = l+1
          endif
        enddo
        itmp = iand(itmp-1_8,itmp)
      enddo
    enddo
  endif
  idx(0) = l-1
 end
 ! subroutine filter_connected_davidson_shortcut(key1,shortcut,key2,Nint,sze,idx)
 !   use bitmasks
 !   implicit none
 !   BEGIN_DOC
 !   ! Filters out the determinants that are not connected by H
 !   ! returns the array idx which contains the index of the 
 !   ! determinants in the array key1 that interact 
 !   ! via the H operator with key2.
 !   !
 !   ! idx(0) is the number of determinants that interact with key1
 !   ! key1 should come from psi_det_sorted_ab.
 !   END_DOC
 !   integer, intent(in)            :: Nint, sze
 !   integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
 !   integer(bit_kind), intent(in)  :: key2(Nint,2)
 !   integer, intent(out)           :: idx(0:sze)
 !   
 !   integer,intent(in)             :: shortcut(0:sze+1)
 !   
 !   integer                        :: i,j,k,l
 !   integer                        :: degree_x2
 !   integer :: i_alpha, i_beta, exc_a, exc_b, endloop
 !   integer(bit_kind) :: tmp1, tmp2
 !   
 !   ASSERT (Nint > 0)
 !   ASSERT (sze >= 0)
 ! 
 !   l=1
 !   i_alpha = 0
 !   
 !   if (Nint==1) then
 !     do while(shortcut(i_alpha+1) < sze)
 !       i_alpha = i_alpha + 1
 !       exc_a = popcnt(xor(key1(1,1,shortcut(i_alpha)), key2(1,1)))
 !       if(exc_a > 4) then
 !         cycle
 !       end if
 !       endloop = min(shortcut(i_alpha+1)-1, sze)
 !       if(exc_a == 4) then
 !         do i_beta = shortcut(i_alpha), endloop
 !           if(key1(1,2,i_beta) == key2(1,2)) then
 !             idx(l) = i_beta
 !             l = l + 1
 !             exit
 !           end if
 !         end do
 !       else
 !         do i_beta = shortcut(i_alpha), endloop
 !           exc_b = popcnt(xor(key1(1,2,i_beta), key2(1,2)))
 !           if(exc_b + exc_a <= 4) then
 !             idx(l) = i_beta
 !             l = l + 1
 !           end if
 !         end do
 !       end if
 !     end do
 !   else
 !     print *, "TBD : filter_connected_davidson_shortcut Nint>1"
 !     stop
 !   end if
 !     
 !   idx(0) = l-1
 ! end
 ! 
 ! subroutine filter_connected_davidson(key1,key2,Nint,sze,idx)
 !   use bitmasks
 !   implicit none
 !   BEGIN_DOC
 !   ! Filters out the determinants that are not connected by H
 !   ! returns the array idx which contains the index of the 
 !   ! determinants in the array key1 that interact 
 !   ! via the H operator with key2.
 !   !
 !   ! idx(0) is the number of determinants that interact with key1
 !   ! key1 should come from psi_det_sorted_ab.
 !   END_DOC
 !   integer, intent(in)            :: Nint, sze
 !   integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
 !   integer(bit_kind), intent(in)  :: key2(Nint,2)
 !   integer, intent(inout)           :: idx(0:sze)
 !   
 !   integer                        :: i,j,k,l
 !   integer                        :: degree_x2
 !   integer :: j_int, j_start
 !   integer*8 :: itmp
 !   
 !   PROVIDE N_con_int det_connections
 !   
 !   
 !   ASSERT (Nint > 0)
 !   ASSERT (sze >= 0)
 ! 
 !   l=1
 !   
 !   if (Nint==1) then
 ! 
 !     i = idx(0) ! lecture dans un intent(out) ?
 !     do j_int=1,N_con_int
 !       itmp = det_connections(j_int,i)
 !       do while (itmp /= 0_8)
 !         j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
 !         do j = j_start+1, min(j_start+32,i-1)
 !           degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +             &
 !               popcnt(xor( key1(1,2,j), key2(1,2)))
 !           if (degree_x2 > 4) then
 !             cycle
 !           else
 !             idx(l) = j
 !             l = l+1
 !           endif
 !         enddo
 !         itmp = iand(itmp-1_8,itmp)
 !       enddo
 !     enddo
 !     
 !   else if (Nint==2) then
 !     
 !     
 !     i = idx(0)
 !     do j_int=1,N_con_int
 !       itmp = det_connections(j_int,i)
 !       do while (itmp /= 0_8)
 !         j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
 !         do j = j_start+1, min(j_start+32,i-1)
 !           degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
 !               popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
 !               popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
 !               popcnt(xor( key1(2,2,j), key2(2,2)))
 !           if (degree_x2 > 4) then
 !             cycle
 !           else
 !             idx(l) = j
 !             l = l+1
 !           endif
 !         enddo
 !         itmp = iand(itmp-1_8,itmp)
 !       enddo
 !     enddo
 !     
 !   else if (Nint==3) then
 !     
 !     i = idx(0)
 !     !DIR$ LOOP COUNT (1000)
 !     do j_int=1,N_con_int
 !       itmp = det_connections(j_int,i)
 !       do while (itmp /= 0_8)
 !         j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
 !         do j = j_start+1, min(j_start+32,i-1)
 !           degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
 !               popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
 !               popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
 !               popcnt(xor( key1(2,2,j), key2(2,2))) +                 &
 !               popcnt(xor( key1(3,1,j), key2(3,1))) +                 &
 !               popcnt(xor( key1(3,2,j), key2(3,2)))
 !           if (degree_x2 > 4) then
 !             cycle
 !           else
 !             idx(l) = j
 !             l = l+1
 !           endif
 !         enddo
 !         itmp = iand(itmp-1_8,itmp)
 !       enddo
 !     enddo
 !     
 !   else
 !     
 !     i = idx(0)
 !     !DIR$ LOOP COUNT (1000)
 !     do j_int=1,N_con_int
 !       itmp = det_connections(j_int,i)
 !       do while (itmp /= 0_8)
 !         j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
 !         do j = j_start+1, min(j_start+32,i-1)
 !           degree_x2 = 0
 !           !DEC$ LOOP COUNT MIN(4)
 !           do k=1,Nint
 !             degree_x2 = degree_x2+ popcnt(xor( key1(k,1,j), key2(k,1))) +&
 !                 popcnt(xor( key1(k,2,j), key2(k,2)))
 !             if (degree_x2 > 4) then
 !               exit
 !             endif
 !           enddo
 !           if (degree_x2 <= 5) then
 !             idx(l) = j
 !             l = l+1
 !           endif
 !         enddo
 !         itmp = iand(itmp-1_8,itmp)
 !       enddo
 !     enddo
 ! 
 !   endif
 !   idx(0) = l-1
 ! end
 subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
  use bitmasks
  BEGIN_DOC
--- a/src/Determinants/s2.irp.f
+++ b/src/Determinants/s2.irp.f
@ -107,35 +107,103 @@ subroutine get_s2_u0_old(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
 end
 subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
- implicit none
+  implicit none
- use bitmasks
+  use bitmasks
- integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
+  integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
- integer, intent(in) :: n,nmax
+  integer, intent(in) :: n,nmax
- double precision, intent(in) :: psi_coefs_tmp(nmax)
+  double precision, intent(in) :: psi_coefs_tmp(nmax)
- double precision, intent(out) :: s2
+  double precision, intent(out) :: s2
- double precision :: s2_tmp
+  double precision :: s2_tmp
- integer :: i,j,l,jj
+  integer :: i,j,l,jj,ii
- integer, allocatable           :: idx(:)
+  integer, allocatable           :: idx(:)
- s2 = 0.d0
+ 
- !$OMP PARALLEL DEFAULT(NONE)                                        &
+  integer                        :: shortcut(0:n+1), sort_idx(n)
-     !$OMP PRIVATE(i,j,s2_tmp,idx)                                   &
+  integer(bit_kind)              :: sorted(N_int,n), version(N_int,n)
-     !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold)&
+  integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
-     !$OMP REDUCTION(+:s2)
+  double precision               :: davidson_threshold_bis
- allocate(idx(0:n))
+
- !$OMP DO SCHEDULE(dynamic)
+  !PROVIDE davidson_threshold
- do i=1,n
+  
-   idx(0) = i
+  s2 = 0.d0
-   call filter_connected_davidson(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx)
+  davidson_threshold_bis = davidson_threshold
-   do jj=1,idx(0)
+  !$OMP PARALLEL DEFAULT(NONE)                                        &
-     j = idx(jj)
+      !$OMP PRIVATE(i,j,s2_tmp,idx,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)                                   &
-     if ( dabs(psi_coefs_tmp(j)) + dabs(psi_coefs_tmp(i))            &
+      !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)&
-           > davidson_threshold ) then
+      !$OMP REDUCTION(+:s2)
-       call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
+  allocate(idx(0:n))
-       s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
+  
-     endif
+  
-   enddo
+  !$OMP SINGLE
- enddo
+  call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
- !$OMP END DO
+  !$OMP END SINGLE
  !$OMP DO SCHEDULE(dynamic)
  do sh=1,shortcut(0)
  do sh2=1,sh
    exa = 0
    do ni=1,N_int
      exa += popcnt(xor(version(ni,sh), version(ni,sh2)))
    end do
    if(exa > 2) then
      cycle
    end if
    do i=shortcut(sh),shortcut(sh+1)-1
      if(sh==sh2) then
        endi = i-1
      else
        endi = shortcut(sh2+1)-1
      end if
      do j=shortcut(sh2),endi
        ext = exa
        do ni=1,N_int
          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
        end do
        if(ext <= 4) then
          org_i = sort_idx(i)
          org_j = sort_idx(j)
          if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))            &
                > davidson_threshold ) then
            call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int)
            s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp
          endif
        end if
      end do
    end do
  end do
  enddo
  !$OMP END DO
  !$OMP SINGLE
  call sort_dets_ba_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
  !$OMP END SINGLE
  !$OMP DO SCHEDULE(dynamic)
  do sh=1,shortcut(0)
    do i=shortcut(sh),shortcut(sh+1)-1    
      do j=shortcut(sh),i-1
        ext = 0
        do ni=1,N_int
          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
        end do
        if(ext == 4) then
          org_i = sort_idx(i)
          org_j = sort_idx(j)
          if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))            &
                > davidson_threshold ) then
            call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int)
            s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp
          endif
        end if
      end do
    end do          
  enddo
  !$OMP END DO
 deallocate(idx)
 !$OMP END PARALLEL
 s2 = s2+s2
@ -146,3 +214,4 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
 s2 = s2 + S_z2_Sz
 end
--- a/src/Determinants/slater_rules.irp.f
+++ b/src/Determinants/slater_rules.irp.f
@ -138,6 +138,7 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
  end select
 end
 subroutine get_double_excitation(det1,det2,exc,phase,Nint)
  use bitmasks
  implicit none
@ -787,6 +788,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
  ASSERT (Ndet > 0)
  ASSERT (Ndet_max >= Ndet)
  i_H_psi_array = 0.d0
  call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
  do ii=1,idx(0)
    i = idx(ii)
@ -1214,7 +1216,7 @@ subroutine get_occ_from_key(key,occ,Nint)
 end
-subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
+subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint)
  use bitmasks
  implicit none
  BEGIN_DOC
@ -1232,29 +1234,56 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
  integer, allocatable           :: idx(:)
  double precision               :: hij
  double precision, allocatable  :: vt(:)
-  integer                        :: i,j,k,l, jj
+  integer                        :: i,j,k,l, jj,ii,sh
  integer                        :: i0, j0
  integer,intent(in)             :: shortcut(0:n+1), sort_idx(n)
  integer                        :: tmp, warp(2,0:n+1), ni
  ASSERT (Nint > 0)
  ASSERT (Nint == N_int)
  ASSERT (n>0)
  PROVIDE ref_bitmask_energy
  !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt)                             &
+      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,warp,tmp,sh)                             &
-      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold)
+      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,shortcut,sort_idx)
  allocate(idx(0:n), vt(n))
  Vt = 0.d0
  v_0 = 0.d0
-  !$OMP DO SCHEDULE(guided)
+
-  do i=1,n
+  !$OMP DO SCHEDULE(dynamic)
-    idx(0) = i
+  
-    call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
+
-    do jj=1,idx(0)
+  do sh=1,shortcut(0)
-      j = idx(jj)
+    warp(1,0) = 0
-      if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then
+    do ii=1,sh!shortcut(0)
-        call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
+      tmp = 0
-        vt (i) = vt (i) + hij*u_0(j)
+      do ni=1,Nint
-        vt (j) = vt (j) + hij*u_0(i)
+        tmp = popcnt(xor(keys_tmp(ni,1, shortcut(ii)), keys_tmp(ni,1,shortcut(sh))))
-      endif
+      end do
      if(tmp <= 4) then
        warp(1,0) = warp(1,0) + 1
        warp(1,warp(1,0)) = shortcut(ii)
        warp(2,warp(1,0)) = shortcut(ii+1)-1
      end if
    end do
    do ii=shortcut(sh),shortcut(sh+1)-1
      idx(0) = ii
      call filter_connected_davidson_warp(keys_tmp,warp,keys_tmp(1,1,ii),Nint,ii-1,idx)
      i = sort_idx(ii)
      do jj=1,idx(0)
        j = sort_idx(idx(jj))
        if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then
          call i_H_j(keys_tmp(1,1,idx(jj)),keys_tmp(1,1,ii),Nint,hij)
          vt (i) = vt (i) + hij*u_0(j)
          vt (j) = vt (j) + hij*u_0(i)
        endif
      enddo
    enddo
  enddo
  !$OMP END DO
@ -1266,180 +1295,8 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
  deallocate(idx,vt)
  !$OMP END PARALLEL
  do i=1,n
-   v_0(i) += H_jj(i) * u_0(i)
+    v_0(i) += H_jj(i) * u_0(i)
  enddo
 end
 BEGIN_PROVIDER [ integer, N_con_int ]
  implicit none
  BEGIN_DOC
  ! Number of integers to represent the connections between determinants
  END_DOC
  N_con_int = 1 + ishft(N_det-1,-11)
 END_PROVIDER
 BEGIN_PROVIDER [ integer*8, det_connections, (N_con_int,N_det) ]
  implicit none
  BEGIN_DOC
  ! Build connection proxy between determinants
  END_DOC
  integer                        :: i,j
  integer                        :: degree
  integer                        :: j_int, j_k, j_l
  integer, allocatable           :: idx(:)
  integer                        :: thread_num
  integer                        :: omp_get_thread_num
  PROVIDE progress_bar
  call start_progress(N_det,'Det connections',0.d0)
  select case(N_int)
    case(1)
      !$OMP PARALLEL DEFAULT (NONE)                                  &
          !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections, &
          !$OMP progress_bar,progress_value)&
          !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num)
      !$  thread_num = omp_get_thread_num()
      allocate (idx(0:N_det))
      !$OMP DO SCHEDULE(guided)
      do i=1,N_det
        if (thread_num == 0) then
          progress_bar(1) = i
          progress_value = dble(i)
        endif
        do j_int=1,N_con_int
          det_connections(j_int,i) = 0_8
          j_k = ishft(j_int-1,11)
          do j_l = j_k,min(j_k+2047,N_det), 32
            do j = j_l+1,min(j_l+32,i)
              degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + &
                  popcnt(xor( psi_det(1,2,i),psi_det(1,2,j)))
              if (degree < 5) then
                det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) )
                exit
              endif
            enddo
          enddo
        enddo
      enddo
      !$OMP ENDDO
      deallocate(idx)
      !$OMP END PARALLEL
    case(2)
      !$OMP PARALLEL DEFAULT (NONE)                                  &
          !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,&
          !$OMP progress_bar,progress_value)&
          !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num)
      !$  thread_num = omp_get_thread_num()
      allocate (idx(0:N_det))
      !$OMP DO SCHEDULE(guided)
      do i=1,N_det
        if (thread_num == 0) then
          progress_bar(1) = i
          progress_value = dble(i)
        endif
        do j_int=1,N_con_int
          det_connections(j_int,i) = 0_8
          j_k = ishft(j_int-1,11)
          do j_l = j_k,min(j_k+2047,N_det), 32
            do j = j_l+1,min(j_l+32,i)
              degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + &
                  popcnt(xor( psi_det(1,2,i),psi_det(1,2,j))) +      &
                  popcnt(xor( psi_det(2,1,i),psi_det(2,1,j))) +      &
                  popcnt(xor( psi_det(2,2,i),psi_det(2,2,j)))
              if (degree < 5) then
                det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) )
                exit
              endif
            enddo
          enddo
        enddo
      enddo
      !$OMP ENDDO
      deallocate(idx)
      !$OMP END PARALLEL
    case(3)
      !$OMP PARALLEL DEFAULT (NONE)                                  &
          !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,&
          !$OMP progress_bar,progress_value)&
          !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num)
      !$  thread_num = omp_get_thread_num()
      allocate (idx(0:N_det))
      !$OMP DO SCHEDULE(guided)
      do i=1,N_det
        if (thread_num == 0) then
          progress_bar(1) = i
          progress_value = dble(i)
        endif
        do j_int=1,N_con_int
          det_connections(j_int,i) = 0_8
          j_k = ishft(j_int-1,11)
          do j_l = j_k,min(j_k+2047,N_det), 32
            do j = j_l+1,min(j_l+32,i)
              degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + &
                  popcnt(xor( psi_det(1,2,i),psi_det(1,2,j))) +      &
                  popcnt(xor( psi_det(2,1,i),psi_det(2,1,j))) +      &
                  popcnt(xor( psi_det(2,2,i),psi_det(2,2,j))) +      &
                  popcnt(xor( psi_det(3,1,i),psi_det(3,1,j))) +      &
                  popcnt(xor( psi_det(3,2,i),psi_det(3,2,j)))
              if (degree < 5) then
                det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) )
                exit
              endif
            enddo
          enddo
        enddo
      enddo
      !$OMP ENDDO
      deallocate(idx)
      !$OMP END PARALLEL
    case default
      !$OMP PARALLEL DEFAULT (NONE)                                  &
          !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,&
          !$OMP progress_bar,progress_value)&
          !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num)
      !$  thread_num = omp_get_thread_num()
      allocate (idx(0:N_det))
      !$OMP DO SCHEDULE(guided)
      do i=1,N_det
        if (thread_num == 0) then
          progress_bar(1) = i
          progress_value = dble(i)
        endif
        do j_int=1,N_con_int
          det_connections(j_int,i) = 0_8
          j_k = ishft(j_int-1,11)
          do j_l = j_k,min(j_k+2047,N_det), 32
            do j = j_l+1,min(j_l+32,i)
              !DIR$ FORCEINLINE
              call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
              if (degree < 3) then
                det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) )
                exit
              endif
            enddo
          enddo
        enddo
      enddo
      !$OMP ENDDO
      deallocate(idx)
      !$OMP END PARALLEL
  end select
  call stop_progress
 END_PROVIDER
--- a/src/Ezfio_files/README.rst
+++ b/src/Ezfio_files/README.rst
@ -203,119 +203,87 @@ Documentation
  Output file for Bitmask
 `output_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L41>`_
  Output file for CAS_SD
 `output_cis <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L61>`_
  Output file for CIS
 `output_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L81>`_
  Output file for CISD
 `output_cisd_selected <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L101>`_
  Output file for CISD_selected
 `output_cpu_time_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f#L2>`_
  Initial CPU and wall times when printing in the output files
-`output_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L121>`_
+`output_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L41>`_
  Output file for Determinants
-`output_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L141>`_
+`output_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L61>`_
  Output file for Electrons
-`output_ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L161>`_
+`output_ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L81>`_
  Output file for Ezfio_files
-`output_fcidump <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L181>`_
+`output_full_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L101>`_
  Output file for FCIdump
 `output_full_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L201>`_
  Output file for Full_CI
-`output_generators_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L221>`_
+`output_generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L121>`_
  Output file for Generators_CAS
 `output_generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L241>`_
  Output file for Generators_full
-`output_hartree_fock <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L261>`_
+`output_hartree_fock <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L141>`_
  Output file for Hartree_Fock
-`output_integrals_bielec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L281>`_
+`output_integrals_bielec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L161>`_
  Output file for Integrals_Bielec
-`output_integrals_monoelec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L301>`_
+`output_integrals_monoelec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L181>`_
  Output file for Integrals_Monoelec
-`output_mo_basis <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L321>`_
+`output_mo_basis <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L201>`_
  Output file for MO_Basis
-`output_moguess <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L341>`_
+`output_moguess <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L221>`_
  Output file for MOGuess
-`output_mrcc_cassd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L361>`_
+`output_mrcc_cassd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L241>`_
  Output file for MRCC_CASSD
-`output_mrcc_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L381>`_
+`output_mrcc_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L261>`_
  Output file for MRCC_Utils
-`output_nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L401>`_
+`output_nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L281>`_
  Output file for Nuclei
-`output_perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L421>`_
+`output_perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L301>`_
  Output file for Perturbation
-`output_properties <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L441>`_
+`output_properties <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L321>`_
  Output file for Properties
-`output_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L461>`_
+`output_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L341>`_
  Output file for Pseudo
-`output_psiref_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L481>`_
+`output_psiref_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L361>`_
  Output file for Psiref_CAS
-`output_psiref_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L501>`_
+`output_psiref_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L381>`_
  Output file for Psiref_Utils
-`output_qmcchem <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L521>`_
+`output_selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L401>`_
  Output file for QmcChem
 `output_selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L541>`_
  Output file for Selectors_full
-`output_singlerefmethod <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L561>`_
+`output_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L421>`_
  Output file for SingleRefMethod
 `output_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L581>`_
  Output file for Utils