Hartree-Fock works well

2024-06-25 22:52:15 +02:00 · 2014-06-25 14:58:58 +02:00 · 2014-06-25 14:58:58 +02:00 · ff86d51c5f
commit ff86d51c5f
parent f3fc0fdb8a
39 changed files with 751 additions and 296 deletions
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@ -216,8 +216,9 @@ class H_apply(object):
      self.data["size_max"] = str(1024*128) 
      self.data["copy_buffer"] = """
      call copy_h_apply_buffer_to_wf
-      selection_criterion_min = selection_criterion_min*0.1d0
+      selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0
      selection_criterion = selection_criterion_min
+      call write_double(output_Dets,selection_criterion,'Selection criterion')
      """
      self.data["keys_work"] = """
      e_2_pert_buffer = 0.d0
@ -230,16 +231,17 @@ class H_apply(object):
      self.data["omp_parallel"]    += """&
 !$OMP REDUCTION (max:select_max_out)"""
      self.data["skip"] = """
-      if ((i_generator < size(select_max)).and. &
-        (select_max(i_generator) < selection_criterion_min*selection_criterion_factor)) then
-        !$ call omp_set_lock(lck)
-        do k=1,N_st
-          norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
-          delta_pt2(k) = 0.d0
-          pt2_old(k) = 0.d0
-        enddo
-        !$ call omp_unset_lock(lck)
-        cycle
+      if (i_generator < size_select_max) then
+        if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then
+          !$ call omp_set_lock(lck)
+          do k=1,N_st
+            norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+            delta_pt2(k) = 0.d0
+            pt2_old(k) = 0.d0
+          enddo
+          !$ call omp_unset_lock(lck)
+          cycle
+        endif
      endif
      select_max(i_generator) = 0.d0
      """
--- a/setup_environment.sh
+++ b/setup_environment.sh
@ -23,9 +23,9 @@ export QPACKAGE_ROOT=${QPACKAGE_ROOT}
 export PYTHONPATH+=:\${QPACKAGE_ROOT}/scripts
 export PATH+=:\${QPACKAGE_ROOT}/scripts
 export PATH+=:\${QPACKAGE_ROOT}/bin
-export QPACKAGE_CACHE_URL="http://qmcchem.ups-tlse.fr/files/scemama/quantum_package/cache
+export QPACKAGE_CACHE_URL="http://qmcchem.ups-tlse.fr/files/scemama/quantum_package/cache"
 export PATH+=:${QPACKAGE_ROOT}/irpf90/bin/
-"
+source ${QPACKAGE_ROOT}/irpf90/bin/irpman
 EOF

 source quantum_package.rc
--- a/src/AOs/ASSUMPTIONS.rst
+++ b/src/AOs/ASSUMPTIONS.rst
@ -4,5 +4,5 @@

   \int \left(\chi_i({\bf r}) \right)^2 dr = 1

-* The AO coefficients in the EZFIO files are not normalized
+* The AO coefficients in the EZFIO files are not necessarily normalized and are normalized after reading
 * The AO coefficients and exponents are ordered in increasing order of exponents
--- a/src/AOs/README.rst
+++ b/src/AOs/README.rst
@ -54,17 +54,17 @@ Documentation
  Coefficients, exponents and powers of x,y and z
  ao_coef(i,j) = coefficient of the jth primitive on the ith ao

-`ao_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L136>`_
+`ao_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L96>`_
  Transposed ao_coef and ao_expo

 `ao_expo <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L20>`_
  Coefficients, exponents and powers of x,y and z
  ao_coef(i,j) = coefficient of the jth primitive on the ith ao

-`ao_expo_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L137>`_
+`ao_expo_transp <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L97>`_
  Transposed ao_coef and ao_expo

-`ao_nucl <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L186>`_
+`ao_nucl <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L146>`_
  Index of the nuclei on which the ao is centered

 `ao_num <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L1>`_
@ -73,21 +73,17 @@ Documentation
 `ao_num_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L2>`_
  Number of atomic orbitals

-`ao_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L96>`_
-  matrix of the overlap for tha AOs
-  S(i,j) = overlap between the ith and the jth atomic basis function
-
 `ao_power <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L19>`_
  Coefficients, exponents and powers of x,y and z
  ao_coef(i,j) = coefficient of the jth primitive on the ith ao

-`ao_prim_num <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L154>`_
+`ao_prim_num <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L114>`_
  Number of primitives per atomic orbital

-`ao_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L176>`_
+`ao_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L136>`_
  Undocumented

-`ao_prim_num_max_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L177>`_
+`ao_prim_num_max_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L137>`_
  Undocumented


--- a/src/BiInts/README.rst
+++ b/src/BiInts/README.rst
@ -36,58 +36,58 @@ Documentation
  integral of the AO basis <ik|jl> or (ij|kl)
  i(r1) j(r1) 1/r12 k(r2) l(r2)

-`ao_bielec_integral_schwartz <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L326>`_
+`ao_bielec_integral_schwartz <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L341>`_
  Needed to compuet Schwartz inequalities

-`ao_bielec_integrals_in_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L189>`_
+`ao_bielec_integrals_in_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L188>`_
  Map of Atomic integrals
  i(r1) j(r2) 1/r12 k(r1) l(r2)

-`compute_ao_bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L148>`_
+`compute_ao_bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L147>`_
  Compute AO 1/r12 integrals for all i and fixed j,k,l

-`eri <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L487>`_
+`eri <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L502>`_
  ATOMIC PRIMTIVE bielectronic integral between the 4 primitives ::
  primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2)
  primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2)
  primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2)
  primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2)

-`general_primitive_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L352>`_
+`general_primitive_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L367>`_
  Computes the integral <pq|rs> where p,q,r,s are Gaussian primitives

-`give_polynom_mult_center_x <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L632>`_
+`give_polynom_mult_center_x <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L647>`_
  subroutine that returns the explicit polynom in term of the "t"
  variable of the following polynomw :
  I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q)

-`i_x1_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L576>`_
+`i_x1_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L591>`_
  recursive function involved in the bielectronic integral

-`i_x1_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L695>`_
+`i_x1_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L710>`_
  recursive function involved in the bielectronic integral

-`i_x1_pol_mult_a1 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L815>`_
+`i_x1_pol_mult_a1 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L830>`_
  recursive function involved in the bielectronic integral

-`i_x1_pol_mult_a2 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L869>`_
+`i_x1_pol_mult_a2 <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L884>`_
  recursive function involved in the bielectronic integral

-`i_x1_pol_mult_recurs <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L729>`_
+`i_x1_pol_mult_recurs <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L744>`_
  recursive function involved in the bielectronic integral

-`i_x2_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L599>`_
+`i_x2_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L614>`_
  recursive function involved in the bielectronic integral

-`i_x2_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L931>`_
+`i_x2_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L946>`_
  recursive function involved in the bielectronic integral

-`integrale_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L531>`_
+`integrale_new <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L546>`_
  calculate the integral of the polynom ::
  I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q)
  between ( 0 ; 1)

-`n_pt_sup <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L618>`_
+`n_pt_sup <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/ao_bi_integrals.irp.f#L633>`_
  Returns the upper boundary of the degree of the polynom involved in the
  bielctronic integral :
  Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z)
@ -112,7 +112,7 @@ Documentation
 `clear_ao_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L128>`_
  Frees the memory of the AO map

-`clear_mo_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L243>`_
+`clear_mo_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L285>`_
  Frees the memory of the MO map

 `get_ao_bielec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L33>`_
@ -136,7 +136,11 @@ Documentation
  Returns multiple integrals <ij|kl> in the MO basis, all
  i for j,k,l fixed.

-`get_mo_map_size <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L235>`_
+`get_mo_bielec_integrals_existing_ik <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L235>`_
+  Returns multiple integrals <ij|kl> in the MO basis, all
+  i for j,k,l fixed.
+
+`get_mo_map_size <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L277>`_
  Return the number of elements in the MO map

 `insert_into_ao_integrals_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L153>`_
@ -154,13 +158,13 @@ Documentation
 `add_integrals_to_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L40>`_
  Adds integrals to tha MO map according to some bitmask

-`mo_bielec_integral_jj <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L296>`_
+`mo_bielec_integral_jj <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L305>`_
  Transform Bi-electronic integrals <ij|ij> and <ij|ji>

-`mo_bielec_integral_jj_anti <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L298>`_
+`mo_bielec_integral_jj_anti <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L307>`_
  Transform Bi-electronic integrals <ij|ij> and <ij|ji>

-`mo_bielec_integral_jj_exchange <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L297>`_
+`mo_bielec_integral_jj_exchange <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L306>`_
  Transform Bi-electronic integrals <ij|ij> and <ij|ji>

 `mo_bielec_integrals_in_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L21>`_
@ -169,22 +173,22 @@ Documentation
 `mo_bielec_integrals_index <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/mo_bi_integrals.irp.f#L1>`_
  Computes an unique index for i,j,k,l integrals

-`ao_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L73>`_
+`ao_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L69>`_
  If <pq|rs> < ao_integrals_threshold, <pq|rs> = 0

-`do_direct_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L111>`_
+`do_direct_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L105>`_
  If True, compute integrals on the fly

-`mo_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L92>`_
+`mo_integrals_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L87>`_
  If <ij|kl> < mo_integrals_threshold, <ij|kl> = 0

-`read_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L55>`_
+`read_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L52>`_
  If true, read AO integrals in EZFIO

-`read_mo_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L37>`_
+`read_mo_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L35>`_
  If true, read MO integrals in EZFIO

-`write_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L19>`_
+`write_ao_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L18>`_
  If true, write AO integrals in EZFIO

 `write_mo_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/options.irp.f#L1>`_
--- a/src/Bitmask/README.rst
+++ b/src/Bitmask/README.rst
@ -57,7 +57,7 @@ Documentation
 `full_ijkl_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L12>`_
  Bitmask to include all possible MOs

-`generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L91>`_
+`generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L88>`_
  Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator).
  3rd index is :
  * 1 : hole     for single exc
@ -70,10 +70,10 @@ Documentation
 `hf_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L32>`_
  Hartree Fock bit mask

-`i_bitmask_gen <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L182>`_
+`i_bitmask_gen <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L174>`_
  Current bitmask for the generators

-`i_bitmask_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L190>`_
+`i_bitmask_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L182>`_
  Current bitmask for the reference

 `n_generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L58>`_
@ -82,13 +82,13 @@ Documentation
 `n_int <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L3>`_
  Number of 64-bit integers needed to represent determinants as binary strings

-`n_reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L126>`_
+`n_reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L122>`_
  Number of bitmasks for reference

 `ref_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L50>`_
  Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask

-`reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L159>`_
+`reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L152>`_
  Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference)

 `bitstring_to_hexa <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L95>`_
--- a/src/Bitmask/bitmasks.irp.f
+++ b/src/Bitmask/bitmasks.irp.f
@ -28,6 +28,18 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int,4) ]
  enddo
 END_PROVIDER

+ 
+BEGIN_PROVIDER [ integer(bit_kind), cis_ijkl_bitmask, (N_int,4) ]
+  implicit none
+  BEGIN_DOC
+  ! Bitmask to include all possible single excitations from Hartree-Fock
+  END_DOC
+  
+  integer                        :: i,j,n
+  cis_ijkl_bitmask = full_ijkl_bitmask
+  cis_ijkl_bitmask(:,1) = HF_bitmask(:,1)
+END_PROVIDER
+

 BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)]
  implicit none
--- a/src/CIS/ASSUMPTIONS.rst
+++ b/src/CIS/ASSUMPTIONS.rst
@ -0,0 +1 @@
+* The molecular orbitals are assumed orthonormal
--- a/src/CIS/H_apply.irp.f
+++ b/src/CIS/H_apply.irp.f
@ -0,0 +1,9 @@
+! Generates subroutine H_apply_cisd
+! ----------------------------------
+
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import H_apply
+H = H_apply("cis",do_double_exc=False)
+print H
+END_SHELL
+
--- a/src/CIS/Makefile
+++ b/src/CIS/Makefile
@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common
--- a/src/CIS/README.rst
+++ b/src/CIS/README.rst
@ -0,0 +1,42 @@
+==========
+CIS Module
+==========
+
+Assumptions
+===========
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* The molecular orbitals are assumed orthonormal
+
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+
--- a/src/CIS/cis.irp.f
+++ b/src/CIS/cis.irp.f
@ -0,0 +1,17 @@
+program cis
+  implicit none
+  integer :: i
+
+  print *,   'HF      = ', HF_energy
+  print *,  'N_states = ', N_states
+  call H_apply_cis
+  print *,  'N_det = ', N_det
+  do i = 1,N_states
+   print *,  'energy  = ',CI_energy(i) 
+   print *,  'E_corr  = ',CI_electronic_energy(i) - ref_bitmask_energy
+  enddo
+  psi_coef = ci_eigenvectors
+  SOFT_TOUCH psi_coef
+  call save_wavefunction
+
+end
--- a/src/CIS/super_ci.irp.f
+++ b/src/CIS/super_ci.irp.f
@ -0,0 +1,61 @@
+program cis
+  implicit none
+  integer                        :: i
+  
+  call super_CI
+  
+end
+
+subroutine super_CI
+  implicit none
+  double precision               :: E, delta_E, delta_D, E_min
+  integer                        :: k
+  character                      :: save_char
+
+  call write_time(output_hartree_fock)
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'),          &
+      '====','================','================','================'
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'),          &
+      '  N ', 'Energy  ', 'Energy diff  ', 'Save  '
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'),          &
+      '====','================','================','================'
+  
+  E = HF_energy + 1.d0
+  delta_D = 0.d0
+  E_min = HF_energy
+  FREE psi_det psi_coef
+  call clear_mo_map
+  N_det = 1
+  SOFT_TOUCH N_det
+  mo_coef = eigenvectors_fock_matrix_mo
+  TOUCH mo_coef
+  do k=1,n_it_scf_max
+    delta_E = HF_energy - E
+    E = HF_energy
+    if (E < E_min) then
+      call save_mos
+      save_char = 'X'
+    else
+      save_char = ' '
+    endif
+    write(output_hartree_fock,'(I4,X,F16.10, X, F16.10, X, A8 )'),&
+        k, E, delta_E, save_char
+    if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then
+      exit
+    endif
+    call H_apply_cis
+    call diagonalize_CI
+    call set_natural_mos
+    FREE psi_det psi_coef
+    call clear_mo_map
+    N_det = 1
+    SOFT_TOUCH N_det
+    mo_coef = eigenvectors_fock_matrix_mo
+    TOUCH mo_coef
+  enddo
+  
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'),          &
+      '====','================','================','================'
+  call write_time(output_hartree_fock)
+end
+    
--- a/src/DensityMatrix/README.rst
+++ b/src/DensityMatrix/README.rst
@ -17,12 +17,6 @@ Documentation
 `iunit_two_body_dm_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L4>`_
  Temporary files for 2-body dm calculation

-`one_body_dm_a <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L213>`_
-  Alpha and beta one-body density matrix
-
-`one_body_dm_b <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L214>`_
-  Alpha and beta one-body density matrix
-
 `two_body_dm_diag_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L170>`_
  diagonal part of the two body density matrix

--- a/src/Dets/README.rst
+++ b/src/Dets/README.rst
@ -152,6 +152,21 @@ Documentation
 `davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L374>`_
  Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]

+`one_body_dm_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/density_matrix.irp.f#L73>`_
+  One-body density matrix
+
+`one_body_dm_mo_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/density_matrix.irp.f#L1>`_
+  Alpha and beta one-body density matrix for each state
+
+`one_body_dm_mo_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/density_matrix.irp.f#L2>`_
+  Alpha and beta one-body density matrix for each state
+
+`save_natural_mos <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/density_matrix.irp.f#L81>`_
+  Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
+
+`state_average_weight <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/density_matrix.irp.f#L99>`_
+  Weights in the state-average calculation of the density matrix
+
 `det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L303>`_
  Return an integer*8 corresponding to a determinant index for searching

@ -307,6 +322,9 @@ Documentation
 `s_z2_sz <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L37>`_
  Undocumented

+`save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/save_natorb.irp.f#L1>`_
+  Undocumented
+
 `a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L721>`_
  Needed for diag_H_mat_elem

--- a/src/Dets/density_matrix.irp.f
+++ b/src/Dets/density_matrix.irp.f
@ -78,10 +78,10 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num)
 one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta
 END_PROVIDER

-subroutine save_natural_mos
+subroutine set_natural_mos
 implicit none
 BEGIN_DOC
- ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
+ ! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
 END_DOC
 character*(64) :: label
 double precision, allocatable :: tmp(:,:)
@ -92,9 +92,18 @@ subroutine save_natural_mos
 label = "Natural"
 call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label)
 deallocate(tmp)
- call save_mos

 end
+subroutine save_natural_mos
+ implicit none
+ BEGIN_DOC
+ ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
+ END_DOC
+ call set_natural_mos
+  call save_mos
+ 
+end
+

 BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
 implicit none
--- a/src/Dets/determinants.ezfio_config
+++ b/src/Dets/determinants.ezfio_config
@ -10,3 +10,7 @@ determinants
  n_det_max_jacobi integer
  threshold_generators  double precision
  threshold_selectors   double precision
+  det_num                integer
+  det_occ                integer           (electrons_elec_alpha_num,determinants_det_num,2)
+  det_coef               double precision  (determinants_det_num)
+
--- a/src/Dets/determinants.irp.f
+++ b/src/Dets/determinants.irp.f
@ -56,7 +56,7 @@ BEGIN_PROVIDER [ integer, N_det_max_jacobi ]
 if (exists) then
   call ezfio_get_determinants_n_det_max_jacobi(N_det_max_jacobi)
 else
-   N_det_max_jacobi = 2000
+   N_det_max_jacobi = 10000
 endif
 call write_int(output_dets,N_det_max_jacobi,'Lapack diagonalization up to')
 ASSERT (N_det_max_jacobi > 0)
--- a/src/Dets/save_for_qmcchem.irp.f
+++ b/src/Dets/save_for_qmcchem.irp.f
@ -0,0 +1,37 @@
+subroutine save_dets_qmcchem
+ use bitmasks
+ implicit none
+ character :: c(mo_tot_num)
+ integer :: i,k
+
+ integer, allocatable :: occ(:,:,:), occ_tmp(:,:)
+ !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp
+
+ call ezfio_set_determinants_det_num(N_det)
+ call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1))
+
+ allocate (occ(elec_alpha_num,N_det,2))
+ ! OMP PARALLEL DEFAULT(NONE) &
+ ! OMP PRIVATE(occ_tmp,i,k)&
+ ! OMP SHARED(N_det,psi_det_sorted,elec_alpha_num, &
+ ! OMP   occ,elec_beta_num,N_int)
+ allocate (occ_tmp(N_int*bit_kind_size,2))
+ occ_tmp = 0
+ ! OMP DO 
+ do i=1,N_det
+  call bitstring_to_list(psi_det_sorted(1,1,i), occ_tmp(1,1), elec_alpha_num, N_int )
+  call bitstring_to_list(psi_det_sorted(1,2,i), occ_tmp(1,2), elec_beta_num, N_int )
+  do k=1,elec_alpha_num
+    occ(k,i,1) = occ_tmp(k,1)
+    occ(k,i,2) = occ_tmp(k,2)
+  enddo
+ enddo
+ ! OMP END DO
+ deallocate(occ_tmp)
+ ! OMP END PARALLEL
+ call ezfio_set_determinants_det_occ(occ)
+ call write_int(output_dets,N_det,'Determinants saved for QMC')
+ deallocate(occ)
+end
+
+
--- a/src/Generators_full/README.rst
+++ b/src/Generators_full/README.rst
@ -5,3 +5,45 @@ Generators_full Module
 All the determinants of the wave function are generators. In this way, the Full CI
 space is explored.

+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`n_det_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full/generators.irp.f#L22>`_
+  For Single reference wave functions, the number of generators is 1 : the
+  Hartree-Fock determinant
+
+`psi_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full/generators.irp.f#L44>`_
+  For Single reference wave functions, the generator is the
+  Hartree-Fock determinant
+
+`select_max <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full/generators.irp.f#L60>`_
+  Memo to skip useless selectors
+
+`threshold_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full/generators.irp.f#L3>`_
+  Percentage of the norm of the state-averaged wave function to
+  consider for the generators
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+
--- a/src/Generators_full/generators.irp.f
+++ b/src/Generators_full/generators.irp.f
@ -57,7 +57,15 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]

 END_PROVIDER

- BEGIN_PROVIDER [ double precision, select_max, (3000) ]
+BEGIN_PROVIDER [ integer, size_select_max]
+ implicit none
+ BEGIN_DOC
+ ! Size of the select_max array
+ END_DOC
+ size_select_max = 10000
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision precision, select_max, (size_select_max) ]
 implicit none
 BEGIN_DOC
 ! Memo to skip useless selectors
--- a/src/Hartree_Fock/Fock_matrix.irp.f
+++ b/src/Hartree_Fock/Fock_matrix.irp.f
@ -260,3 +260,36 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
 endif
 END_PROVIDER

+subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
+  implicit none
+  integer, intent(in)            :: LDFMO ! size(FMO,1)
+  integer, intent(in)            :: LDFAO ! size(FAO,1)
+  double precision, intent(in)   :: FMO(LDFMO,*)
+  double precision, intent(out)  :: FAO(LDFAO,*)
+  
+  double precision, allocatable  :: T(:,:), M(:,:)
+  ! F_ao = S C F_mo C^t S
+  allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      ao_overlap, size(ao_overlap,1),                                &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+      M, size(M,1),                                                  &
+      FMO, size(FMO,1),                                              &
+      0.d0,                                                          &
+      T, size(T,1))
+  call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,            &
+      T, size(T,1),                                                  &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      M, size(M,1),                                                  &
+      ao_overlap, size(ao_overlap,1),                                &
+      0.d0,                                                          &
+      FAO, size(FAO,1))
+  deallocate(T,M)
+end
+
--- a/src/Hartree_Fock/README.rst
+++ b/src/Hartree_Fock/README.rst
@ -71,6 +71,9 @@ Documentation
  K = Fb - Fa
  .br

+`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L263>`_
+  Undocumented
+
 `hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L211>`_
  Hartree-Fock energy

@ -83,25 +86,10 @@ Documentation
 `hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L14>`_
  Beta density matrix in the AO basis

-`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L31>`_
+`run <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L7>`_
  Undocumented

-`insert_new_scf_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L19>`_
-  Undocumented
-
-`it_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L1>`_
-  Number of the current SCF iteration
-
-`scf_density_matrices <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L9>`_
-  Density matrices at every SCF iteration
-
-`scf_energies <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L10>`_
-  Density matrices at every SCF iteration
-
-`scf_interpolation_step <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L64>`_
-  Undocumented
-
-`scf_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L89>`_
+`scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L2>`_
  Undocumented

 `diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L1>`_
@ -110,12 +98,6 @@ Documentation
 `eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L2>`_
  Diagonal Fock matrix in the MO basis

-`run <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L7>`_
-  Undocumented
-
-`scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L2>`_
-  Undocumented
-
 `do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L39>`_
  If True, compute integrals on the fly

--- a/src/Hartree_Fock/SCF.irp.f
+++ b/src/Hartree_Fock/SCF.irp.f
@ -1,108 +1,23 @@
-BEGIN_PROVIDER [ integer, it_scf ]
-  implicit none
-  BEGIN_DOC
-  ! Number of the current SCF iteration
-  END_DOC
-  it_scf = 0
-END_PROVIDER

- BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]
-&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]
-  implicit none
-  BEGIN_DOC
-  ! Density matrices at every SCF iteration
-  END_DOC
-  SCF_density_matrices = 0.d0
-  SCF_energies = 0.d0
-END_PROVIDER
-
-subroutine insert_new_SCF_density_matrix
-  implicit none
-  integer                        :: i,j
-  do j=1,ao_num
-    do i=1,ao_num
-      SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j)
-      SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j)
-    enddo
-  enddo
-  SCF_energies(it_scf) = HF_energy
+program scf
+  call orthonormalize_mos
+  call run
 end

-subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
+subroutine run
+
+  use bitmasks
  implicit none
-  integer, intent(in)            :: LDFMO ! size(FMO,1)
-  integer, intent(in)            :: LDFAO ! size(FAO,1)
-  double precision, intent(in)   :: FMO(LDFMO,*)
-  double precision, intent(out)  :: FAO(LDFAO,*)
+  double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
+  double precision               :: E0
+  integer                        :: i_it, i, j, k
+   
+  E0 = HF_energy 
  
-  double precision, allocatable  :: T(:,:), M(:,:)
-  ! F_ao = S C F_mo C^t S
-  allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
-      ao_overlap, size(ao_overlap,1),                                &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
-      M, size(M,1),                                                  &
-      FMO, size(FMO,1),                                              &
-      0.d0,                                                          &
-      T, size(T,1))
-  call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,            &
-      T, size(T,1),                                                  &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0,                                                          &
-      M, size(M,1))
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
-      M, size(M,1),                                                  &
-      ao_overlap, size(ao_overlap,1),                                &
-      0.d0,                                                          &
-      FAO, size(FAO,1))
-  deallocate(T,M)
-end
-
-subroutine SCF_interpolation_step
-  implicit none
-  integer                        :: i,j
-  double precision               :: c
+  thresh_SCF = 1.d-10
+  call damping_SCF
+  mo_label = "Canonical"
+  TOUCH mo_label mo_coef
+  call save_mos
  
-  if (it_scf == 1) then
-    return
-  endif
-  call random_number(c)
-  c = c*0.5d0
-  do j=1,ao_num
-    do i=1,ao_num
-      HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf-1)+SCF_density_matrices(i,j,1,it_scf) * (1.d0 - c)
-      HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf-1)+SCF_density_matrices(i,j,2,it_scf) * (1.d0 - c)
-    enddo
-  enddo
-  TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
-  
-  ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),&
-      !      Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) )
-  ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),&
-      !      Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) )
-  ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta
-end
-
-subroutine scf_iterations
-  implicit none
-  integer                        :: i,j
-  do i=1,n_it_scf_max
-    it_scf += 1
-    SOFT_TOUCH it_scf
-    mo_coef = eigenvectors_Fock_matrix_mo
-    TOUCH mo_coef
-    call insert_new_SCF_density_matrix
-    print *,  HF_energy
-    if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then
-      call SCF_interpolation_step
-    endif
-    if (it_scf>1 ) then
-      if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then
-        exit
-      endif
-    endif
-  enddo
 end
--- a/src/Hartree_Fock/damping_SCF.irp.f
+++ b/src/Hartree_Fock/damping_SCF.irp.f
@ -0,0 +1,123 @@
+subroutine damping_SCF
+  implicit none
+  double precision               :: E
+  double precision, allocatable  :: D_alpha(:,:), D_beta(:,:)
+  double precision               :: E_new
+  double precision, allocatable  :: D_new_alpha(:,:), D_new_beta(:,:), F_new(:,:)
+  double precision, allocatable  :: delta_alpha(:,:), delta_beta(:,:)
+  double precision               :: lambda, E_half, a, b, delta_D, delta_E, E_min
+  
+  integer                        :: i,j,k
+  logical                        :: saving
+  character                      :: save_char
+
+  allocate(                                                          &
+      D_alpha( ao_num_align, ao_num ),                               &
+      D_beta( ao_num_align, ao_num ),                                &
+      F_new( ao_num_align, ao_num ),                                 &
+      D_new_alpha( ao_num_align, ao_num ),                           &
+      D_new_beta( ao_num_align, ao_num ),                            &
+      delta_alpha( ao_num_align, ao_num ),                           &
+      delta_beta( ao_num_align, ao_num ))
+  
+  do j=1,ao_num
+    do i=1,ao_num
+      D_alpha(i,j) = HF_density_matrix_ao_alpha(i,j)
+      D_beta (i,j) = HF_density_matrix_ao_beta (i,j)
+    enddo
+  enddo
+  
+  
+  call write_time(output_hartree_fock)
+
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '===='
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '  N ', 'Energy  ', 'Energy diff  ', 'Density diff  ', 'Save'
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '===='
+
+  E = HF_energy + 1.d0
+  E_min = HF_energy
+  delta_D = 0.d0
+  do k=1,n_it_scf_max
+    
+    delta_E = HF_energy - E
+    E = HF_energy
+    
+    if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then
+      exit
+    endif
+
+    saving = E < E_min
+    if (saving) then
+      call save_mos
+      save_char = 'X'
+      E_min = E
+    else
+      save_char = ' '
+    endif
+
+    write(output_hartree_fock,'(I4,X,F16.10, X, F16.10, X, F16.10, 3X, A )'), &
+      k, E, delta_E, delta_D, save_char
+    
+    D_alpha = HF_density_matrix_ao_alpha
+    D_beta  = HF_density_matrix_ao_beta 
+    mo_coef = eigenvectors_fock_matrix_mo
+    TOUCH mo_coef
+    
+    D_new_alpha = HF_density_matrix_ao_alpha
+    D_new_beta  = HF_density_matrix_ao_beta 
+    F_new = Fock_matrix_ao
+    E_new = HF_energy
+
+    delta_alpha = D_new_alpha - D_alpha
+    delta_beta  = D_new_beta  - D_beta 
+    
+    lambda = 0.5d0
+    E_half = 0.d0
+    do while (E_half > E)
+      HF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha
+      HF_density_matrix_ao_beta  = D_beta  + lambda * delta_beta
+      TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
+      mo_coef = eigenvectors_fock_matrix_mo
+      TOUCH mo_coef
+      E_half = HF_energy
+      if ((E_half > E).and.(E_new < E)) then
+        lambda = 1.d0
+        exit
+      else if ((E_half > E).and.(lambda > 1.d-3)) then
+        lambda = 0.5d0 * lambda
+        E_new = E_half
+      else
+        exit
+      endif
+    enddo
+
+    a = (E_new + E - 2.d0*E_half)*2.d0
+    b = -E_new - 3.d0*E + 4.d0*E_half
+    lambda = -lambda*b/a
+    D_alpha = (1.d0-lambda) * D_alpha + lambda * D_new_alpha
+    D_beta  = (1.d0-lambda) * D_beta  + lambda * D_new_beta 
+    delta_E = HF_energy - E
+    do j=1,ao_num
+      do i=1,ao_num
+        delta_D = delta_D + &
+        (D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j))*(D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j)) + &
+        (D_beta (i,j) - HF_density_matrix_ao_beta (i,j))*(D_beta (i,j) - HF_density_matrix_ao_beta (i,j))
+      enddo
+    enddo
+    delta_D = dsqrt(delta_D/dble(ao_num)**2)
+    HF_density_matrix_ao_alpha = D_alpha
+    HF_density_matrix_ao_beta  = D_beta
+    TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
+    mo_coef = eigenvectors_fock_matrix_mo
+    TOUCH mo_coef
+
+
+  enddo
+  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '===='
+  write(output_hartree_fock,*)
+  
+  call write_double(output_hartree_fock, HF_energy, 'Hartree-Fock energy')
+  call write_time(output_hartree_fock)
+
+  deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta)
+end
--- a/src/Hartree_Fock/mo_SCF_iterations.irp.f
+++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f
@ -1,23 +0,0 @@
-
-program scf
-  call orthonormalize_mos
-  call run
-end
-
-subroutine run
-
-  use bitmasks
-  implicit none
-  double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
-  double precision               :: E0
-  integer                        :: i_it, i, j, k
-   
-  E0 = HF_energy 
-  
-  thresh_SCF = 1.d-10
-  call damping_SCF
-  mo_label = "Canonical"
-  TOUCH mo_label mo_coef
-  call save_mos
-  
-end
--- a/src/MOGuess/NEEDED_MODULES
+++ b/src/MOGuess/NEEDED_MODULES
@ -1,2 +1,3 @@
-AOs Ezfio_files MOs Nuclei Output Utils MonoInts
+AOs Electrons Ezfio_files MonoInts MOs Nuclei Output Utils
+

--- a/src/MOGuess/README.rst
+++ b/src/MOGuess/README.rst
@ -9,12 +9,13 @@ Needed Modules
 .. NEEDED_MODULES file.

 * `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
 * `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
 * `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
 * `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
 * `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
 * `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
-* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_

 Documentation
 =============
@ -22,19 +23,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.

-`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L/subroutine h_core_guess/;">`_
+`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L1>`_
  Undocumented

-`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]/;">`_
+`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L2>`_
  matrix of the coefficients of the mos generated by the
  orthonormalization by the S^{-1/2} canonical transformation of the aos
  ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital

-`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]/;">`_
+`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L26>`_
  overlap matrix of the ao_ortho_lowdin
  supposed to be the Identity

-`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]/;">`_
+`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L1>`_
  Undocumented


--- a/src/MOs/README.rst
+++ b/src/MOs/README.rst
@ -36,12 +36,15 @@ Documentation
 .. NEEDED_MODULES file.

 `cholesky_mo <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L1>`_
-  Cholesky decomposition of MO Density matrix to
+  Cholesky decomposition of AO Density matrix to
  generate MOs

 `mo_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L44>`_
  Density matrix in MO basis

+`mo_density_matrix_virtual <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L64>`_
+  Density matrix in MO basis (virtual MOs)
+
 `mo_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L22>`_
  Molecular orbital coefficients on AO basis set
  mo_coef(i,j) = coefficient of the ith ao on the jth mo
--- a/src/MOs/mos.ezfio_config
+++ b/src/MOs/mos.ezfio_config
@ -1,5 +1,4 @@
 mo_basis
-  ao_num             integer
  mo_tot_num         integer
  mo_coef            double precision  (ao_basis_ao_num,mo_basis_mo_tot_num)
  mo_label           character*(64)
--- a/src/Makefile.config.example
+++ b/src/Makefile.config.example
@ -1,6 +1,6 @@
 OPENMP  =1
 PROFILE =0
-DEBUG  = 1
+DEBUG  = 0

 IRPF90_FLAGS+= --align=32
 FC     = ifort -g 
--- a/src/Output/output.irp.f
+++ b/src/Output/output.irp.f
@ -55,7 +55,7 @@ subroutine write_double(iunit,value,label)
  integer, intent(in)            :: iunit
  double precision               :: value
  character*(*)                  :: label
-  character*(64), parameter      :: f = '(A50,G16.8)'
+  character*(64), parameter      :: f = '(A50,G24.16)'
  character*(50)                 :: newlabel
  write(newlabel,'(A,A)') '* ',trim(label)
  write(iunit,f) newlabel, value
--- a/src/Perturbation/README.rst
+++ b/src/Perturbation/README.rst
@ -82,7 +82,7 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.

-`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L/subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)/;">`_
+`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L1>`_
  compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
  .br
  for the various n_st states.
@ -92,7 +92,7 @@ Documentation
  e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
  .br

-`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
+`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
  compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
  .br
  for the various N_st states.
@ -102,7 +102,7 @@ Documentation
  e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
  .br

-`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
+`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L43>`_
  compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
  .br
  for the various N_st states.
@ -112,7 +112,7 @@ Documentation
  c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
  .br

-`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
+`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L2>`_
  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
  .br
  for the various N_st states,
@ -135,23 +135,23 @@ Documentation
  .br
  H_pert_diag = <HF|H|det_pert> c_pert

-`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/double precision function repeat_all_e_corr(key_in)/;">`_
+`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L90>`_
  Undocumented

-`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,coef_pert_buffer, &>`_
+`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
  Fill the H_apply buffer with determiants for the selection

-`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine remove_small_contributions/;">`_
+`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L84>`_
  Remove determinants with small contributions. N_states is assumed to be
  provided.

-`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/BEGIN_PROVIDER [ double precision, selection_criterion ]/;">`_
+`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L71>`_
  Threshold to select determinants. Set by selection routines.

-`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_factor ]/;">`_
+`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L73>`_
  Threshold to select determinants. Set by selection routines.

-`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_min ]/;">`_
+`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L72>`_
  Threshold to select determinants. Set by selection routines.


--- a/src/SC2/README.rst
+++ b/src/SC2/README.rst
@ -8,67 +8,9 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.

-`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/SC2.irp.f#L1>`_
-  CISD+SC2 method              :: take off all the disconnected terms of a CISD (selected or not)
-  .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
-  Initial guess vectors are not necessarily orthonormal
-
-`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/SC2.irp.f#L169>`_
-  Undocumented
-
 `cisd <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/cisd_SC2.irp.f#L1>`_
  Undocumented

-`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L19>`_
-  Eigenvectors/values of the CI matrix
-
-`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L18>`_
-  Eigenvectors/values of the CI matrix
-
-`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L1>`_
-  N_states lowest eigenvalues of the CI matrix
-
-`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L38>`_
-  Replace the coefficients of the CI states by the coefficients of the
-  eigenstates of the CI matrix
-
-`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/pert_sc2.irp.f#L2>`_
-  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  .br
-  for the various N_st states,
-  .br
-  but  with the correction in the denominator
-  .br
-  comming from the interaction of that determinant with all the others determinants
-  .br
-  that can be repeated by repeating all the double excitations
-  .br
-  : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  .br
-  that could be repeated to this determinant.
-  .br
-  In addition, for the perturbative energetic contribution you have the standard second order
-  .br
-  e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
-  .br
-  and also the purely projected contribution
-  .br
-  H_pert_diag = <HF|H|det_pert> c_pert
-
-`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/pert_sc2.irp.f#L82>`_
-  Undocumented
-


 Needed Modules
--- a/src/Selectors_full/README.rst
+++ b/src/Selectors_full/README.rst
@ -2,3 +2,96 @@
 Selectors_full Module
 =====================

+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`coef_hf_selector <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/e_corr_selectors.irp.f#L27>`_
+  energy of correlation per determinant respect to the Hartree Fock determinant
+  .br
+  for the all the double excitations in the selectors determinants
+  .br
+  E_corr_per_selectors(i) = <D_i|H|HF> * c(D_i)/c(HF) if |D_i> is a double excitation
+  .br
+  E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation
+  .br
+  coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants
+
+`double_index_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/e_corr_selectors.irp.f#L4>`_
+  degree of excitation respect to Hartree Fock for the wave function
+  .br
+  for the all the selectors determinants
+  .br
+  double_index_selectors = list of the index of the double excitations
+  .br
+  n_double_selectors = number of double excitations in the selectors determinants
+
+`e_corr_per_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/e_corr_selectors.irp.f#L28>`_
+  energy of correlation per determinant respect to the Hartree Fock determinant
+  .br
+  for the all the double excitations in the selectors determinants
+  .br
+  E_corr_per_selectors(i) = <D_i|H|HF> * c(D_i)/c(HF) if |D_i> is a double excitation
+  .br
+  E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation
+  .br
+  coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants
+
+`exc_degree_per_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/e_corr_selectors.irp.f#L3>`_
+  degree of excitation respect to Hartree Fock for the wave function
+  .br
+  for the all the selectors determinants
+  .br
+  double_index_selectors = list of the index of the double excitations
+  .br
+  n_double_selectors = number of double excitations in the selectors determinants
+
+`n_double_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/e_corr_selectors.irp.f#L5>`_
+  degree of excitation respect to Hartree Fock for the wave function
+  .br
+  for the all the selectors determinants
+  .br
+  double_index_selectors = list of the index of the double excitations
+  .br
+  n_double_selectors = number of double excitations in the selectors determinants
+
+`n_det_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/selectors.irp.f#L26>`_
+  For Single reference wave functions, the number of selectors is 1 : the
+  Hartree-Fock determinant
+
+`psi_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/selectors.irp.f#L48>`_
+  Determinants on which we apply <i|H|psi> for perturbation.
+
+`psi_selectors_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/selectors.irp.f#L49>`_
+  Determinants on which we apply <i|H|psi> for perturbation.
+
+`psi_selectors_size <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/selectors.irp.f#L21>`_
+  Undocumented
+
+`threshold_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full/selectors.irp.f#L3>`_
+  Percentage of the norm of the state-averaged wave function to
+  consider for the selectors
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+
--- a/src/SingleRefMethod/README.rst
+++ b/src/SingleRefMethod/README.rst
@ -19,6 +19,9 @@ Documentation
  For Single reference wave functions, the generator is the
  Hartree-Fock determinant

+`select_max <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod/generators.irp.f#L43>`_
+  Memo to skip useless selectors
+


 Needed Modules
--- a/src/SingleRefMethod/generators.irp.f
+++ b/src/SingleRefMethod/generators.irp.f
@ -28,6 +28,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
   call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int)
   if (degree == 0) then
     k = j
+     exit
   endif
 end do

--- a/src/Utils/LinearAlgebra.irp.f
+++ b/src/Utils/LinearAlgebra.irp.f
@ -48,7 +48,6 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
    stop
  endif
  
-  
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
  !$OMP PRIVATE(i,j,k)
@ -178,7 +177,7 @@ subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n)
  call dgemm('N','N',m,n,n,1.d0,A,LDA,R,LDR,0.d0,B,LDB)
 end

-subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
+subroutine lapack_diagd(eigvalues,eigvectors,H,nmax,n)
  implicit none
  BEGIN_DOC
  ! Diagonalize matrix H
@ -244,3 +243,122 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
  enddo
  deallocate(A,eigenvalues)
 end
+
+subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
+  implicit none
+  BEGIN_DOC
+  ! Diagonalize matrix H
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+  double precision, intent(out)  :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+  double precision, intent(in)   :: H(nmax,n)
+  double precision,allocatable   :: eigenvalues(:)
+  double precision,allocatable   :: work(:)
+  double precision,allocatable   :: A(:,:)
+  integer                        :: lwork, info, i,j,l,k, liwork
+
+  allocate(A(nmax,n),eigenvalues(n))
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  A=H
+  lwork = 2*n*n + 6*n+ 1
+  allocate (work(lwork))
+
+  lwork = -1
+  call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    info )
+  if (info < 0) then
+    print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( work( 1 ) )
+  deallocate (work)
+
+  allocate (work(lwork))
+  call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    info )
+  deallocate(work)
+
+  if (info < 0) then
+    print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'DSYEV Failed'
+     stop 1
+  end if
+
+  eigvectors = 0.d0
+  eigvalues = 0.d0
+  do j = 1, n
+    eigvalues(j) = eigenvalues(j)
+    do i = 1, n
+      eigvectors(i,j) = A(i,j)
+    enddo
+  enddo
+  deallocate(A,eigenvalues)
+end
+
+subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st)
+  implicit none
+  BEGIN_DOC
+  ! Diagonalize matrix H
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax,n_st
+  double precision, intent(out)  :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+  double precision, intent(in)   :: H(nmax,n)
+  double precision,allocatable   :: work(:)
+  integer         ,allocatable   :: iwork(:), isuppz(:)
+  double precision,allocatable   :: A(:,:)
+  integer                        :: lwork, info, i,j,l,k,m, liwork
+
+  allocate(A(nmax,n))
+
+  A=H
+  lwork = 2*n*n + 6*n+ 1
+  liwork = 5*n + 3
+  allocate (work(lwork),iwork(liwork),isuppz(2*N_st))
+
+  lwork = -1
+  liwork = -1
+  call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, &
+    iwork, liwork, info )
+  if (info < 0) then
+    print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( work( 1 ) )
+  liwork = iwork(1)
+  deallocate (work,iwork)
+
+  allocate (work(lwork),iwork(liwork))
+  call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, &
+    iwork, liwork, info )
+  deallocate(work,iwork)
+
+  if (info < 0) then
+    print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'DSYEVR Failed'
+     stop 1
+  end if
+
+  deallocate(A)
+end
--- a/src/Utils/README.rst
+++ b/src/Utils/README.rst
@ -92,7 +92,7 @@ Documentation
  into
  fact_k  (x-x_P)^iorder(1)  (y-y_P)^iorder(2)  (z-z_P)^iorder(3) exp(-p(r-P)^2)

-`hermite <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L472>`_
+`hermite <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L477>`_
  Hermite polynomial

 `multiply_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L201>`_
@ -108,10 +108,10 @@ Documentation
  \int_0^1 dx \exp(-p x^2) x^n
  .br

-`rint1 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L528>`_
+`rint1 <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L533>`_
  Standard version of rint

-`rint_large_n <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L497>`_
+`rint_large_n <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L502>`_
  Version of rint for large values of n

 `rint_sum <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L421>`_
				`@ -0,0 +1 @@`
				`* The molecular orbitals are assumed orthonormal`