diff --git a/config/gfortran_avx.cfg b/config/gfortran_avx.cfg
index 80bbbec9..f065d133 100644
--- a/config/gfortran_avx.cfg
+++ b/config/gfortran_avx.cfg
@@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -ffree-line-length-none -I . -mavx -g
+FC           : gfortran -ffree-line-length-none -I . -mavx -g 
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32
@@ -35,7 +35,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #                           -ffast-math and the Fortran-specific
 #                           -fno-protect-parens and -fstack-arrays.
 [OPT]
-FCFLAGS : -Ofast 
+FCFLAGS : -Ofast -march=native
 
 # Profiling flags
 #################
diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg
index 4b06c5e9..f0c6e320 100644
--- a/config/gfortran_debug.cfg
+++ b/config/gfortran_debug.cfg
@@ -51,7 +51,7 @@ FCFLAGS : -Ofast
 # -g           : Extra debugging information
 #
 [DEBUG]
-FCFLAGS : -g -msse4.2 
+FCFLAGS : -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant 
 
 # OpenMP flags
 #################
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index 76080b02..6cc83745 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@@ -7,6 +7,7 @@ module Determinants_by_hand : sig
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -18,11 +19,14 @@ module Determinants_by_hand : sig
   val to_rst : t -> Rst_string.t
   val of_rst : Rst_string.t -> t option
   val read_n_int : unit -> N_int_number.t
+  val update_ndet : Det_number.t -> unit
+  val extract_state : States_number.t -> unit
 end = struct
   type t = 
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
+      n_states               : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -129,12 +133,12 @@ end = struct
     |> Array.map ~f:Det_coef.of_float
   ;;
 
-  let write_psi_coef ~n_det c =
+  let write_psi_coef ~n_det ~n_states c =
     let n_det = Det_number.to_int n_det
     and c = Array.to_list c
             |> List.map ~f:Det_coef.to_float
     and n_states = 
-      read_n_states () |> States_number.to_int
+      States_number.to_int n_states
     in
     Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
     |> Ezfio.set_determinants_psi_coef 
@@ -200,6 +204,7 @@ end = struct
           expected_s2            = read_expected_s2 ()          ;
           psi_coef               = read_psi_coef ()             ;
           psi_det                = read_psi_det ()              ;
+          n_states               = read_n_states ()             ;
         }
     else
       failwith "No molecular orbitals, so no determinants"
@@ -222,12 +227,14 @@ end = struct
               expected_s2          ;
               psi_coef             ;
               psi_det              ;
+              n_states             ;
             } =
      write_n_int n_int ;
      write_bit_kind bit_kind;
      write_n_det n_det;
+     write_n_states n_states;
      write_expected_s2 expected_s2;
-     write_psi_coef ~n_det:n_det psi_coef ;
+     write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
      write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
   ;;
 
@@ -298,6 +305,7 @@ Determinants ::
 n_int                  = %s
 bit_kind               = %s
 n_det                  = %s
+n_states               = %s
 expected_s2            = %s
 psi_coef               = %s
 psi_det                = %s
@@ -305,6 +313,7 @@ psi_det                = %s
      (b.n_int         |> N_int_number.to_string)
      (b.bit_kind      |> Bit_kind.to_string)
      (b.n_det         |> Det_number.to_string)
+     (b.n_states      |> States_number.to_string)
      (b.expected_s2   |> Positive_float.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")
@@ -433,14 +442,83 @@ psi_det                = %s
       |> Bit_kind.to_int)
     and n_int =
       Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
+    and n_states =
+      Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ())
     in
     let s = 
-       String.concat [ header ; bitkind ; n_int ; psi_coef ; psi_det]
+       String.concat [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det]
     in
 
+
+
+
     Generic_input_of_rst.evaluate_sexp t_of_sexp s
   ;;
 
+  let update_ndet n_det_new =
+    Printf.printf "Reducing n_det to %d\n" (Det_number.to_int n_det_new);
+    let n_det_new =
+      Det_number.to_int n_det_new
+    in
+    let det =
+      read ()
+    in
+    let n_det_old, n_states =
+      Det_number.to_int det.n_det,
+      States_number.to_int det.n_states
+    in
+    if n_det_new = n_det_old then
+      ()
+    ;
+    if n_det_new > n_det_new then
+      failwith @@ Printf.sprintf "Requested n_det should be less than %d" n_det_old
+    ;
+    for j=0 to (n_states-1) do
+      let ishift_old, ishift_new =
+        j*n_det_old,
+        j*n_det_new
+      in
+      for i=0 to (n_det_new-1) do
+        det.psi_coef.(i+ishift_new) <- det.psi_coef.(i+ishift_old)
+      done
+    done
+    ;
+    let new_det =
+      { det with n_det = (Det_number.of_int n_det_new) }
+    in
+    write new_det
+  ;;
+
+  let extract_state istate =
+    Printf.printf "Extracting state %d\n" (States_number.to_int istate);
+    let det =
+      read ()
+    in
+    let n_det, n_states =
+      Det_number.to_int det.n_det,
+      States_number.to_int det.n_states
+    in
+    if (States_number.to_int istate) > n_states then
+      failwith "State to extract should not be greater than n_states"
+    ;
+    let j =
+      (States_number.to_int istate) - 1
+    in
+    begin
+      if (j>0) then
+        let ishift =
+          j*n_det
+        in
+        for i=0 to (n_det-1) do
+          det.psi_coef.(i) <- det.psi_coef.(i+ishift)
+        done
+    end;
+    let new_det =
+      { det with n_states = (States_number.of_int 1) }
+    in
+    write new_det
+  ;;
+
 end
 
 
diff --git a/ocaml/Symmetry.ml b/ocaml/Symmetry.ml
index 5849e116..8647ae99 100644
--- a/ocaml/Symmetry.ml
+++ b/ocaml/Symmetry.ml
@@ -85,7 +85,7 @@ module Xyz = struct
   let of_string s =
     let flush state accu number =
       let n = 
-        if (number = "") then 0
+        if (number = "") then 1
         else (Int.of_string number) 
       in
       match state with
diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml
index a1625719..abc2de1d 100644
--- a/ocaml/TaskServer.ml
+++ b/ocaml/TaskServer.ml
@@ -47,6 +47,14 @@ let debug str =
 let zmq_context =
   ZMQ.Context.create ()
 
+let () =
+  let nproc =
+    match Sys.getenv "OMP_NUM_THREADS" with
+    | Some m -> int_of_string m
+    | None -> 2
+  in
+  ZMQ.Context.set_io_threads zmq_context nproc
+
 
 let bind_socket ~socket_type ~socket ~port =
   let rec loop = function
diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml
index c79bf550..7c07ffe5 100644
--- a/ocaml/qp_create_ezfio_from_xyz.ml
+++ b/ocaml/qp_create_ezfio_from_xyz.ml
@@ -21,6 +21,9 @@ let spec =
      ~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)"
   +> anon ("(xyz_file|zmt_file)" %: file )
 
+type element =
+| Element of Element.t
+| Int_elem of (Nucl_number.t * Element.t)
 
 (** Handle dummy atoms placed on bonds *)
 let dummy_centers ~threshold ~molecule ~nuclei =
@@ -115,17 +118,14 @@ let run ?o b c d m p cart xyz_file =
   (* Open basis set channels *)
   let basis_channel element =
     let key =
-      Element.to_string element
+      match element with
+      | Element e -> Element.to_string e
+      | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i)  (Element.to_string e)
     in
     match Hashtbl.find basis_table key with
     | Some in_channel -> 
         in_channel
-    | None -> 
-         let msg = 
-           Printf.sprintf "%s is not defined in basis %s.%!"
-           (Element.to_long_string element) b ;
-         in
-         failwith msg
+    | None -> raise Not_found
   in
   
   let temp_filename =
@@ -189,12 +189,21 @@ let run ?o b c d m p cart xyz_file =
       | Some (key, basis) -> (*Aux basis *)
         begin
           let elem  = 
-            Element.of_string key
+            try
+              Element (Element.of_string key)
+            with Element.ElementError _ ->
+              let result = 
+                match (String.split ~on:',' key) with
+                | i :: k :: [] -> (Nucl_number.of_int @@ int_of_string i, Element.of_string k)
+                | _ -> failwith "Expected format is int,Element:basis"
+              in Int_elem result
           and basis =
             String.lowercase basis
           in
           let key = 
-             Element.to_string elem
+             match elem with
+             | Element e -> Element.to_string e
+             | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e)
           in
           let new_channel =
             fetch_channel basis
@@ -202,7 +211,13 @@ let run ?o b c d m p cart xyz_file =
           begin
             match Hashtbl.add basis_table ~key:key ~data:new_channel with
             | `Ok -> ()
-            | `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem))
+            | `Duplicate ->
+              let e = 
+                match elem with
+                | Element e -> e
+                | Int_elem (_,e) -> e
+              in
+              failwith ("Duplicate definition of basis for "^(Element.to_long_string e))
           end
        end
     end;
@@ -537,7 +552,20 @@ let run ?o b c d m p cart xyz_file =
               | Element.X -> Element.H
               | e -> e
             in
-            Basis.read_element (basis_channel x.Atom.element) i e
+            let key =
+              Int_elem (i,x.Atom.element)
+            in
+            try
+              Basis.read_element (basis_channel key) i e
+            with Not_found ->
+              let key =
+                Element x.Atom.element
+              in
+              try
+                Basis.read_element (basis_channel key) i e
+              with Not_found ->
+                failwith (Printf.sprintf "Basis not found for atom %d (%s)" (Nucl_number.to_int i)
+                 (Element.to_string x.Atom.element) )
           with
           | End_of_file -> failwith
                 ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
@@ -647,6 +675,7 @@ atoms are taken from the same basis set, otherwise specific elements can be
 defined as follows:
 
  -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\"
+ -b \"cc-pvtz | 1,H:sto-3g | 3,H:6-31g\"
 
 If a file with the same name as the basis set exists, this file will be read.
 Otherwise, the basis set is obtained from the database.
diff --git a/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f b/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
index e200322f..ff5dd509 100644
--- a/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
+++ b/plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
@@ -41,8 +41,8 @@ subroutine run_selection_slave(thread,iproc,energy)
     if (done) then
       ctask = ctask - 1
     else
-      integer :: i_generator, i_generator_start, i_generator_max, step, N
-      read (task,*) i_generator_start, i_generator_max, step, N
+      integer :: i_generator, N
+      read (task,*) i_generator, N
       if(buf%N == 0) then
         ! Only first time 
         call create_selection_buffer(N, N*2, buf)
@@ -50,9 +50,7 @@ subroutine run_selection_slave(thread,iproc,energy)
       else
         if(N /= buf%N) stop "N changed... wtf man??"
       end if
-      do i_generator=i_generator_start,i_generator_max,step
-        call select_connected(i_generator,energy,pt2,buf)
-      enddo
+      call select_connected(i_generator,energy,pt2,buf)
     endif
 
     if(done .or. ctask == size(task_id)) then
diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f
index c5285bdf..ddad71db 100644
--- a/plugins/CAS_SD_ZMQ/selection.irp.f
+++ b/plugins/CAS_SD_ZMQ/selection.irp.f
@@ -671,13 +671,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       if(mat(1, p1, p2) == 0d0) cycle
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
       logical, external              :: is_in_wavefunction
-      if (is_in_wavefunction(det,N_int)) then
-        stop 'is_in_wf'
-        cycle
-      endif
       
       if (do_ddci) then
-        integer, external              :: is_a_two_holes_two_particles
+        logical, external              :: is_a_two_holes_two_particles
         if (is_a_two_holes_two_particles(det)) then
           cycle
         endif
@@ -1219,35 +1215,42 @@ subroutine ZMQ_selection(N_in, pt2)
   
   implicit none
   
-  character*(512)                :: task 
   integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
   integer, intent(in)            :: N_in
   type(selection_buffer)         :: b
   integer                        :: i, N
   integer, external              :: omp_get_thread_num
   double precision, intent(out)  :: pt2(N_states)
+  integer, parameter             :: maxtasks=10000
   
   
+  N = max(N_in,1)
   if (.True.) then
     PROVIDE pt2_e0_denominator
-    N = max(N_in,1)
     provide nproc
     call new_parallel_job(zmq_to_qp_run_socket,"selection")
     call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
-    call zmq_set_running(zmq_to_qp_run_socket)
     call create_selection_buffer(N, N*2, b)
   endif
 
-  integer :: i_generator, i_generator_start, i_generator_max, step
+  character*(20*maxtasks) :: task
+  task = ' '
 
-  step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num ))
-  step = max(1,step)
-  do i= 1, N_det_generators,step
-    i_generator_start = i
-    i_generator_max = min(i+step-1,N_det_generators)
-    write(task,*) i_generator_start, i_generator_max, 1, N
+  integer :: k
+  k=0
+  do i= 1, N_det_generators
+    k = k+1
+    write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
+    k = k+20
+    if (k>20*maxtasks) then
+       k=0
+       call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+    endif
+  enddo
+  if (k > 0) then
     call add_task_to_taskserver(zmq_to_qp_run_socket,task)
-  end do
+  endif
+  call zmq_set_running(zmq_to_qp_run_socket)
 
   !$OMP PARALLEL DEFAULT(shared)  SHARED(b, pt2)  PRIVATE(i) NUM_THREADS(nproc+1)
   i = omp_get_thread_num()
@@ -1264,6 +1267,7 @@ subroutine ZMQ_selection(N_in, pt2)
     if (s2_eig) then
       call make_s2_eigenfunction
     endif
+    call save_wavefunction
   endif
 end subroutine
 
diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f
index e1296ce0..03a0ace5 100644
--- a/plugins/DFT_Utils/grid_density.irp.f
+++ b/plugins/DFT_Utils/grid_density.irp.f
@@ -1,7 +1,7 @@
-BEGIN_PROVIDER [integer, n_points_angular_grid]
- implicit none
- n_points_angular_grid = 50
-END_PROVIDER 
+!BEGIN_PROVIDER [integer, n_points_integration_angular_lebedev]
+!implicit none
+!n_points_integration_angular_lebedev = 50
+!END_PROVIDER 
 
 BEGIN_PROVIDER [integer, n_points_radial_grid]
  implicit none
@@ -9,15 +9,15 @@ BEGIN_PROVIDER [integer, n_points_radial_grid]
 END_PROVIDER 
 
 
- BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_angular_grid,3) ]
-&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_angular_grid)]
+ BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular_lebedev,3) ]
+&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular_lebedev)]
  implicit none
  BEGIN_DOC
 ! weights and grid points for the integration on the angular variables on 
 ! the unit sphere centered on (0,0,0)
 ! According to the LEBEDEV scheme
  END_DOC
- call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points)
+ call cal_quad(n_points_integration_angular_lebedev, angular_quadrature_points,weights_angular_points)
  include 'constants.include.F'
  integer :: i
  double precision :: accu
@@ -63,7 +63,7 @@ END_PROVIDER
 
 END_PROVIDER 
 
-BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid,n_points_radial_grid,nucl_num)]
+BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num)]
  BEGIN_DOC
 ! points for integration over space
  END_DOC
@@ -79,7 +79,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid
    double precision :: x,r
    x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1]
    r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration 
-   do k = 1, n_points_angular_grid  ! explicit values of the grid points centered around each atom 
+   do k = 1, n_points_integration_angular_lebedev  ! explicit values of the grid points centered around each atom 
     grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r
     grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * r
     grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r
@@ -88,7 +88,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid
  enddo
 END_PROVIDER 
 
-BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
+BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
  BEGIN_DOC 
 ! Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
 ! the "n" discrete variable represents the nucleis which in this array is represented by the last dimension 
@@ -102,7 +102,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
  ! run over all points in space
   do j = 1, nucl_num  ! that are referred to each atom 
    do k = 1, n_points_radial_grid -1  !for each radial grid attached to the "jth" atom
-    do l = 1, n_points_angular_grid ! for each angular point attached to the "jth" atom
+    do l = 1, n_points_integration_angular_lebedev ! for each angular point attached to the "jth" atom
      r(1) = grid_points_per_atom(1,l,k,j)
      r(2) = grid_points_per_atom(2,l,k,j)
      r(3) = grid_points_per_atom(3,l,k,j)
@@ -123,8 +123,8 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
 
 END_PROVIDER 
 
- BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
-&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
+ BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
+&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
  implicit none
  integer :: i,j,k,l,m
  double precision :: contrib
@@ -132,7 +132,7 @@ END_PROVIDER
  double precision :: aos_array(ao_num),mos_array(mo_tot_num)
   do j = 1, nucl_num
    do k = 1, n_points_radial_grid -1
-    do l = 1, n_points_angular_grid
+    do l = 1, n_points_integration_angular_lebedev
      one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0
      one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0
      r(1) = grid_points_per_atom(1,l,k,j)
diff --git a/plugins/DFT_Utils/integration_radial.irp.f b/plugins/DFT_Utils/integration_radial.irp.f
index 4943783b..e9500ded 100644
--- a/plugins/DFT_Utils/integration_radial.irp.f
+++ b/plugins/DFT_Utils/integration_radial.irp.f
@@ -4,7 +4,7 @@
  double precision :: accu
  integer :: i,j,k,l
  double precision :: x
- double precision :: integrand(n_points_angular_grid), weights(n_points_angular_grid)
+ double precision :: integrand(n_points_integration_angular_lebedev), weights(n_points_integration_angular_lebedev)
  double precision :: f_average_angular_alpha,f_average_angular_beta
  double precision :: derivative_knowles_function,knowles_function
 
@@ -12,7 +12,7 @@
  ! according ot equation (6) of the paper of Becke (JCP, (88), 1988)
  ! Here the m index is referred to the w_m(r) weight functions of equation (22)
    ! Run over all points of integrations : there are  
-   ! n_points_radial_grid (i) * n_points_angular_grid (k) 
+   ! n_points_radial_grid (i) * n_points_integration_angular_lebedev (k) 
    do j = 1, nucl_num 
     integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
     integral_density_beta_knowles_becke_per_atom(j) = 0.d0
@@ -20,7 +20,7 @@
      ! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
      f_average_angular_alpha = 0.d0
      f_average_angular_beta = 0.d0
-     do k = 1, n_points_angular_grid
+     do k = 1, n_points_integration_angular_lebedev
       f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
       f_average_angular_beta  += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j)  * weight_functions_at_grid_points(k,i,j)
      enddo
diff --git a/plugins/Full_CI_ZMQ/energy.irp.f b/plugins/Full_CI_ZMQ/energy.irp.f
index 7b7cc75b..5f9baf46 100644
--- a/plugins/Full_CI_ZMQ/energy.irp.f
+++ b/plugins/Full_CI_ZMQ/energy.irp.f
@@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ]
  ! E0 in the denominator of the PT2
  END_DOC
  if (initialize_pt2_E0_denominator) then
-  pt2_E0_denominator(1:N_states) = CI_electronic_energy(1:N_states)
+  pt2_E0_denominator(1:N_states) = psi_energy(1:N_states)
 ! pt2_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
 ! pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
   call write_double(6,pt2_E0_denominator(1)+nuclear_repulsion, 'PT2 Energy denominator')
diff --git a/plugins/Full_CI_ZMQ/pt2_stoch.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch.irp.f
index 8ffd8e60..914e7138 100644
--- a/plugins/Full_CI_ZMQ/pt2_stoch.irp.f
+++ b/plugins/Full_CI_ZMQ/pt2_stoch.irp.f
@@ -1,8 +1,7 @@
 program pt2_stoch
   implicit none
-  initialize_pt2_E0_denominator = .False.
   read_wf = .True.
-  SOFT_TOUCH initialize_pt2_E0_denominator read_wf
+  SOFT_TOUCH read_wf
   PROVIDE mo_bielec_integrals_in_map
   call run
 end
@@ -17,31 +16,23 @@ subroutine run
   integer                        :: n_det_before, to_select
   double precision               :: threshold_davidson_in
 
-  double precision               :: E_CI_before(N_states), relative_error
+  double precision               :: E_CI_before, relative_error
 
-  if (.true.) then
-    call ezfio_get_full_ci_zmq_energy(E_CI_before(1))
-    pt2_e0_denominator(:) = E_CI_before(1) - nuclear_repulsion
-    SOFT_TOUCH pt2_e0_denominator read_wf
-  endif
   allocate (pt2(N_states))
   pt2 = 0.d0
   
+  E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion
   threshold_selectors = 1.d0
   threshold_generators = 1d0 
-  relative_error = 1.d-6
+  relative_error = 1.d-3
   call ZMQ_pt2(pt2, relative_error)
   print *,  'Final step'
   print *,  'N_det    = ', N_det
-  print *,  'N_states = ', N_states
-  do k=1,N_states
-    print *, 'State', k
-    print *,  'PT2      = ', pt2
-    print *,  'E        = ', E_CI_before
-    print *,  'E+PT2    = ', E_CI_before+pt2
-    print *,  '-----'
-  enddo
-  call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
+  print *,  'PT2      = ', pt2
+  print *,  'E        = ', E_CI_before
+  print *,  'E+PT2    = ', E_CI_before+pt2
+  print *,  '-----'
+  call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before+pt2(1))
 end
 
 
diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
index eb706ccb..afb1a50c 100644
--- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
+++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
@@ -20,7 +20,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
   double precision, allocatable :: pt2_detail(:,:), comb(:)
   logical, allocatable :: computed(:)
   integer, allocatable :: tbc(:)
-  integer :: i, j, Ncomb, generator_per_task, i_generator_end
+  integer :: i, j, k, Ncomb, generator_per_task, i_generator_end
   integer, external :: pt2_find
   
   double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
@@ -32,7 +32,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
   sum2above = 0d0
   Nabove = 0d0
 
-  provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral
+  provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight
 
   !call random_seed()
   
@@ -69,18 +69,18 @@ subroutine ZMQ_pt2(pt2,relative_error)
 
     do i=1,tbc(0)
       if(tbc(i) > fragment_first) then
-        write(task(ipos:ipos+20),'(I9,X,I9,''|'')') 0, tbc(i)
+        write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, tbc(i)
         ipos += 20
-        if (ipos > 64000) then
+        if (ipos > 63980) then
           call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
           ipos=1
           tasks = .True.
         endif
       else
         do j=1,fragment_count
-          write(task(ipos:ipos+20),'(I9,X,I9,''|'')') j, tbc(i)
+          write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, tbc(i)
           ipos += 20
-          if (ipos > 64000) then
+          if (ipos > 63980) then
             call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
             ipos=1
             tasks = .True.
@@ -108,7 +108,12 @@ subroutine ZMQ_pt2(pt2,relative_error)
       call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
 
     else
-       pt2(1) = sum(pt2_detail(1,:))
+       pt2 = 0.d0
+       do i=1,N_det_generators
+         do k=1,N_states
+            pt2(k) = pt2(k) + pt2_detail(k,i)
+         enddo
+       enddo
     endif
 
     tbc(0) = 0
@@ -117,6 +122,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
     endif
   end do
 
+  deallocate(pt2_detail, comb, computed, tbc)
 
 end subroutine
 
@@ -196,11 +202,15 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
   
   allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
     pt2_mwen(N_states, N_det_generators) )
-  actually_computed(:) = computed(:)
+  do i=1,N_det_generators
+    actually_computed(i) = computed(i)
+  enddo
   
   parts_to_get(:) = 1
   if(fragment_first > 0) then
-    parts_to_get(1:fragment_first) = fragment_count
+    do i=1,fragment_first
+      parts_to_get(i) = fragment_count
+    enddo
   endif
 
   do i=1,tbc(0)
@@ -223,7 +233,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
   pullLoop : do while (more == 1)
     call pull_pt2_results(zmq_socket_pull, Nindex, index, pt2_mwen, task_id, ntask)
     do i=1,Nindex
-      pt2_detail(:, index(i)) += pt2_mwen(:,i)
+      pt2_detail(1:N_states, index(i)) += pt2_mwen(1:N_states,i)
       parts_to_get(index(i)) -= 1
       if(parts_to_get(index(i)) < 0) then 
         print *, i, index(i), parts_to_get(index(i)), Nindex
@@ -273,12 +283,11 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
       if (dabs(eqt/avg) < relative_error) then
         pt2(1) = avg
 !        exit pullLoop
+      else
+        print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
       endif
-      print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
     end if
   end do pullLoop
-  print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
-
 
   call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
   call end_zmq_pull_socket(zmq_socket_pull)
@@ -375,9 +384,9 @@ END_PROVIDER
 
 subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
   implicit none
+  integer, intent(inout) :: Ncomb
   double precision, intent(out) :: comb(Ncomb)
   integer, intent(inout) :: tbc(0:size_tbc)
-  integer, intent(inout) :: Ncomb
   logical, intent(inout) :: computed(N_det_generators)
   integer :: i, j, last_full, dets(comb_teeth), tbc_save
   integer :: icount, n
@@ -515,12 +524,15 @@ end subroutine
     pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
   end do
   
-  pt2_weight = pt2_weight / pt2_cweight(N_det_generators)
-  pt2_cweight = pt2_cweight / pt2_cweight(N_det_generators)
+  do i=1,N_det_generators
+    pt2_weight(i)  = pt2_weight(i) / pt2_cweight(N_det_generators)
+    pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(N_det_generators)
+  enddo
   
   norm_left = 1d0
   
   comb_step = 1d0/dfloat(comb_teeth)
+  first_det_of_comb = 1
   do i=1,N_det_generators
     if(pt2_weight(i)/norm_left < comb_step*.5d0) then
       first_det_of_comb = i
diff --git a/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
index 452b446b..5a246319 100644
--- a/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
@@ -25,7 +25,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
   integer :: index
   integer :: Nindex
 
-  allocate(pt2_detail(N_states, N_det))
+  allocate(pt2_detail(N_states, N_det_generators))
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_push      = new_zmq_push_socket(thread)
   call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
@@ -101,7 +101,7 @@ subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntas
   implicit none
 
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_push
-  double precision, intent(in)   :: pt2_detail(N_states, N_det)
+  double precision, intent(in)   :: pt2_detail(N_states, N_det_generators)
   integer, intent(in) :: ntask, N, index, task_id(*)
   integer :: rc
 
@@ -133,7 +133,7 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
   use selection_types
   implicit none
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
-  double precision, intent(inout) :: pt2_detail(N_states, N_det)
+  double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
   integer, intent(out) :: index
   integer, intent(out) :: N, ntask, task_id(*)
   integer :: rc, rn, i
@@ -150,18 +150,22 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
   rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
   if(rc /= 4) stop "pull"
 
-  rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
+  rc = f77_zmq_recv( zmq_socket_pull, task_id, ntask*4, 0)
   if(rc /= 4*ntask) stop "pull"
 
 ! Activate is zmq_socket_pull is a REP
   rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
+
+  do i=N+1,N_det_generators
+    pt2_detail(1:N_states,i) = 0.d0
+  enddo
 end subroutine
  
  
-BEGIN_PROVIDER [ double precision, pt2_workload, (N_det) ]
+BEGIN_PROVIDER [ double precision, pt2_workload, (N_det_generators) ]
   integer :: i
-  do i=1,N_det
-    pt2_workload(:) = dfloat(N_det - i + 1)**2
+  do i=1,N_det_generators
+    pt2_workload(i) = dfloat(N_det_generators - i + 1)**2
   end do
   pt2_workload = pt2_workload / sum(pt2_workload)
 END_PROVIDER
diff --git a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
index 59b2ba1f..85b52c30 100644
--- a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
@@ -26,7 +26,6 @@ subroutine run_selection_slave(thread,iproc,energy)
   call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
   if(worker_id == -1) then
     print *, "WORKER -1"
-    !call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
     call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
     call end_zmq_push_socket(zmq_socket_push,thread)
     return
diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f
index 43f43211..6fd4fd5e 100644
--- a/plugins/Full_CI_ZMQ/selection.irp.f
+++ b/plugins/Full_CI_ZMQ/selection.irp.f
@@ -1110,3 +1110,1115 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint)
   enddo
   
 end
+=======
+use bitmasks
+
+BEGIN_PROVIDER [ integer, fragment_count ]
+  implicit none
+  BEGIN_DOC
+  ! Number of fragments for the deterministic part
+  END_DOC
+  fragment_count = (elec_alpha_num-n_core_orb)**2
+END_PROVIDER
+
+
+double precision function integral8(i,j,k,l)
+  implicit none
+  
+  integer, intent(in) :: i,j,k,l
+  double precision, external :: get_mo_bielec_integral
+  integer :: ii
+  ii = l-mo_integrals_cache_min
+  ii = ior(ii, k-mo_integrals_cache_min)
+  ii = ior(ii, j-mo_integrals_cache_min)
+  ii = ior(ii, i-mo_integrals_cache_min)
+  if (iand(ii, -64) /= 0) then
+    integral8 = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
+  else
+    ii = l-mo_integrals_cache_min
+    ii = ior( ishft(ii,6), k-mo_integrals_cache_min)
+    ii = ior( ishft(ii,6), j-mo_integrals_cache_min)
+    ii = ior( ishft(ii,6), i-mo_integrals_cache_min)
+    integral8 = mo_integrals_cache(ii)
+  endif
+end function
+
+
+BEGIN_PROVIDER [ integer(1), psi_phasemask, (N_int*bit_kind_size, 2, N_det)]
+  use bitmasks
+  implicit none
+  
+  integer :: i
+  do i=1, N_det
+    call get_mask_phase(psi_det_sorted(1,1,i), psi_phasemask(1,1,i))
+  end do
+END_PROVIDER
+
+
+subroutine assert(cond, msg)
+  character(*), intent(in) :: msg
+  logical, intent(in) :: cond
+  
+  if(.not. cond) then
+    print *, "assert failed: "//msg
+    stop
+  end if
+end 
+
+
+subroutine get_mask_phase(det, phasemask)
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: det(N_int, 2)
+  integer(1), intent(out) :: phasemask(2,N_int*bit_kind_size)
+  integer :: s, ni, i
+  logical :: change
+
+  phasemask = 0_1
+  do s=1,2
+    change = .false.
+    do ni=1,N_int
+      do i=0,bit_kind_size-1
+        if(BTEST(det(ni, s), i)) change = .not. change
+          if(change) phasemask(s, (ni-1)*bit_kind_size + i + 1) = 1_1
+      end do
+    end do
+  end do
+end 
+
+
+subroutine select_connected(i_generator,E0,pt2,b,subset)
+  use bitmasks
+  use selection_types
+  implicit none
+  integer, intent(in)            :: i_generator, subset
+  type(selection_buffer), intent(inout) :: b
+  double precision, intent(inout)  :: pt2(N_states)
+  integer :: k,l
+  double precision, intent(in)   :: E0(N_states)
+
+  integer(bit_kind)              :: hole_mask(N_int,2), particle_mask(N_int,2)
+  double precision               :: fock_diag_tmp(2,mo_tot_num+1)
+  
+  call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
+
+  do l=1,N_generators_bitmask
+    do k=1,N_int
+      hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole,l), psi_det_generators(k,1,i_generator))
+      hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator))
+      particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) )
+      particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) )
+
+    enddo
+    call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset)
+  enddo
+end 
+
+
+double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2)
+  use bitmasks
+  implicit none
+
+  integer(1), intent(in) :: phasemask(2,*)
+  integer, intent(in) :: s1, s2, h1, h2, p1, p2
+  logical :: change
+  integer(1) :: np1
+  integer :: np
+  double precision, save :: res(0:1) = (/1d0, -1d0/)
+
+  np1 = phasemask(s1,h1) + phasemask(s1,p1) + phasemask(s2,h2) + phasemask(s2,p2)
+  np = np1
+  if(p1 < h1) np = np + 1
+  if(p2 < h2) np = np + 1
+  
+  if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1
+  get_phase_bi = res(iand(np,1))
+end 
+
+
+
+subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti
+  double precision :: hij
+  double precision, external :: get_phase_bi, integral8
+  
+  integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
+  integer, parameter :: turn2(2) = (/2,1/) 
+  
+  if(h(0,sp) == 2) then
+    h1 = h(1, sp)
+    h2 = h(2, sp)
+    do i=1,3
+      puti = p(i, sp)
+      if(bannedOrb(puti)) cycle
+      p1 = p(turn3_2(1,i), sp)
+      p2 = p(turn3_2(2,i), sp)
+      hij = integral8(p1, p2, h1, h2) - integral8(p2, p1, h1, h2)
+      hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2)
+      vect(:, puti) += hij * coefs
+    end do
+  else if(h(0,sp) == 1) then
+    sfix = turn2(sp)
+    hfix = h(1,sfix)
+    pfix = p(1,sfix)
+    hmob = h(1,sp)
+    do j=1,2
+      puti = p(j, sp)
+      if(bannedOrb(puti)) cycle
+      pmob = p(turn2(j), sp)
+      hij = integral8(pfix, pmob, hfix, hmob)
+      hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix)
+      vect(:, puti) += hij * coefs
+    end do
+  else
+    puti = p(1,sp)
+    if(.not. bannedOrb(puti)) then
+      sfix = turn2(sp)
+      p1 = p(1,sfix)
+      p2 = p(2,sfix)
+      h1 = h(1,sfix)
+      h2 = h(2,sfix)
+      hij = (integral8(p1,p2,h1,h2) - integral8(p2,p1,h1,h2))
+      hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2)
+      vect(:, puti) += hij * coefs
+    end if
+  end if
+end 
+
+
+
+subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i, hole, p1, p2, sh
+  logical :: ok, lbanned(mo_tot_num)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision :: hij
+  double precision, external :: get_phase_bi, integral8
+  
+  lbanned = bannedOrb
+  sh = 1
+  if(h(0,2) == 1) sh = 2
+  hole = h(1, sh)
+  lbanned(p(1,sp)) = .true.
+  if(p(0,sp) == 2) lbanned(p(2,sp)) = .true.
+  !print *, "SPm1", sp, sh
+  
+  p1 = p(1, sp)
+  
+  if(sp == sh) then
+    p2 = p(2, sp)
+    lbanned(p2) = .true.
+    
+    do i=1,hole-1
+      if(lbanned(i)) cycle
+      hij = (integral8(p1, p2, i, hole) - integral8(p2, p1, i, hole))
+      hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2)
+      vect(:,i) += hij * coefs
+    end do
+    do i=hole+1,mo_tot_num
+      if(lbanned(i)) cycle
+      hij = (integral8(p1, p2, hole, i) - integral8(p2, p1, hole, i))
+      hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
+      vect(:,i) += hij * coefs
+    end do
+
+    call apply_particle(mask, sp, p2, det, ok,  N_int)
+    call i_h_j(gen, det, N_int, hij)
+    vect(:, p2) += hij * coefs
+  else
+    p2 = p(1, sh)
+    do i=1,mo_tot_num
+      if(lbanned(i)) cycle
+      hij = integral8(p1, p2, i, hole)
+      hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
+      vect(:,i) += hij * coefs
+    end do
+  end if
+
+  call apply_particle(mask, sp, p1, det, ok,  N_int)
+  call i_h_j(gen, det, N_int, hij)
+  vect(:, p1) += hij * coefs
+end 
+
+
+subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: vect(N_states, mo_tot_num)
+  integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2)
+  integer :: i
+  logical :: ok, lbanned(mo_tot_num)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision :: hij
+  
+  lbanned = bannedOrb
+  lbanned(p(1,sp)) = .true.
+  do i=1,mo_tot_num
+    if(lbanned(i)) cycle
+    call apply_particle(mask, sp, i, det, ok, N_int)
+    call i_h_j(gen, det, N_int, hij)
+    vect(:, i) += hij * coefs
+  end do
+end 
+
+subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in)            :: i_generator, subset
+  integer(bit_kind), intent(in)  :: hole_mask(N_int,2), particle_mask(N_int,2)
+  double precision, intent(in)   :: fock_diag_tmp(mo_tot_num)
+  double precision, intent(in)   :: E0(N_states)
+  double precision, intent(inout) :: pt2(N_states)
+  type(selection_buffer), intent(inout) :: buf
+  
+  double precision                :: mat(N_states, mo_tot_num, mo_tot_num)
+  integer                         :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii
+  integer(bit_kind)               :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2)
+  logical                         :: fullMatch, ok
+  
+  integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2)
+  integer,allocatable               :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:)
+  integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
+  
+  logical :: monoAdo, monoBdo;
+  integer :: maskInd
+
+  PROVIDE fragment_count
+
+  monoAdo = .true.
+  monoBdo = .true.
+  
+  allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det))
+  allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det))
+  
+  do k=1,N_int
+    hole    (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))
+    hole    (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2))
+    particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1))
+    particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2))
+  enddo
+
+  integer                        :: N_holes(2), N_particles(2)
+  integer                        :: hole_list(N_int*bit_kind_size,2)
+  integer                        :: particle_list(N_int*bit_kind_size,2)
+
+  call bitstring_to_list_ab(hole    , hole_list    , N_holes    , N_int)
+  call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
+
+!  ! ======
+!  ! If the subset doesn't exist, return
+!  logical :: will_compute
+!  will_compute = subset == 0
+!
+!  if (.not.will_compute) then
+!    maskInd = N_holes(1)*N_holes(2) + N_holes(2)*((N_holes(2)-1)/2) + N_holes(1)*((N_holes(1)-1)/2)
+!    will_compute = (maskInd >= subset)
+!    if (.not.will_compute) then
+!      return
+!    endif
+!  endif
+!  ! ======
+
+  
+  integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
+  allocate (preinteresting_det(N_int,2,N_det))
+
+  preinteresting(0) = 0
+  prefullinteresting(0) = 0
+  
+  do i=1,N_int
+    negMask(i,1) = not(psi_det_generators(i,1,i_generator))
+    negMask(i,2) = not(psi_det_generators(i,2,i_generator))
+  end do
+  
+  do i=1,N_det
+    mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
+    mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
+    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) 
+    do j=2,N_int
+      mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i))
+      mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
+      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) 
+    end do
+
+    if(nt <= 4) then
+      if(i <= N_det_selectors) then
+        preinteresting(0) += 1
+        preinteresting(preinteresting(0)) = i
+        do j=1,N_int
+          preinteresting_det(j,1,preinteresting(0)) = psi_det_sorted(j,1,i)
+          preinteresting_det(j,2,preinteresting(0)) = psi_det_sorted(j,2,i)
+        enddo
+      else if(nt <= 2) then
+        prefullinteresting(0) += 1
+        prefullinteresting(prefullinteresting(0)) = i
+      end if
+    end if
+  end do
+  
+
+  maskInd = -1
+  integer :: nb_count
+  do s1=1,2
+    do i1=N_holes(s1),1,-1   ! Generate low excitations first
+      
+      h1 = hole_list(i1,s1)
+      call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int)
+      
+      negMask = not(pmask)
+      
+      interesting(0) = 0
+      fullinteresting(0) = 0
+      
+      do ii=1,preinteresting(0)
+        i = preinteresting(ii)
+        mobMask(1,1) = iand(negMask(1,1), preinteresting_det(1,1,ii))
+        mobMask(1,2) = iand(negMask(1,2), preinteresting_det(1,2,ii))
+        nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+        do j=2,N_int
+          mobMask(j,1) = iand(negMask(j,1), preinteresting_det(j,1,ii))
+          mobMask(j,2) = iand(negMask(j,2), preinteresting_det(j,2,ii))
+          nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+        end do
+        
+        if(nt <= 4) then
+          interesting(0) += 1
+          interesting(interesting(0)) = i
+          minilist(1,1,interesting(0)) = preinteresting_det(1,1,ii)
+          minilist(1,2,interesting(0)) = preinteresting_det(1,2,ii)
+          do j=2,N_int
+            minilist(j,1,interesting(0)) = preinteresting_det(j,1,ii)
+            minilist(j,2,interesting(0)) = preinteresting_det(j,2,ii)
+          enddo
+          if(nt <= 2) then
+            fullinteresting(0) += 1
+            fullinteresting(fullinteresting(0)) = i
+            fullminilist(1,1,fullinteresting(0)) = preinteresting_det(1,1,ii)
+            fullminilist(1,2,fullinteresting(0)) = preinteresting_det(1,2,ii)
+            do j=2,N_int
+              fullminilist(j,1,fullinteresting(0)) = preinteresting_det(j,1,ii)
+              fullminilist(j,2,fullinteresting(0)) = preinteresting_det(j,2,ii)
+            enddo
+          end if
+        end if
+      end do
+      
+      do ii=1,prefullinteresting(0)
+        i = prefullinteresting(ii)
+        nt = 0
+        mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
+        mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i))
+        nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+        do j=2,N_int
+          mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i))
+          mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
+          nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+        end do
+
+        if(nt <= 2) then
+          fullinteresting(0) += 1
+          fullinteresting(fullinteresting(0)) = i
+          fullminilist(1,1,fullinteresting(0)) = psi_det_sorted(1,1,i)
+          fullminilist(1,2,fullinteresting(0)) = psi_det_sorted(1,2,i)
+          do j=2,N_int
+            fullminilist(j,1,fullinteresting(0)) = psi_det_sorted(j,1,i)
+            fullminilist(j,2,fullinteresting(0)) = psi_det_sorted(j,2,i)
+          enddo
+        end if
+      end do
+      
+
+
+      do s2=s1,2
+        sp = s1
+        
+        if(s1 /= s2) sp = 3
+        
+        ib = 1
+        if(s1 == s2) ib = i1+1
+        monoAdo = .true.
+        do i2=N_holes(s2),ib,-1   ! Generate low excitations first
+          logical                        :: banned(mo_tot_num, mo_tot_num,2)
+          logical                        :: bannedOrb(mo_tot_num, 2)
+          
+          h2 = hole_list(i2,s2)
+          call apply_hole(pmask, s2,h2, mask, ok, N_int)
+          banned = .false.
+          do j=1,mo_tot_num
+            bannedOrb(j, 1) = .true.
+            bannedOrb(j, 2) = .true.
+          enddo
+          do s3=1,2
+            do i=1,N_particles(s3)
+              bannedOrb(particle_list(i,s3), s3) = .false.
+            enddo
+          enddo
+          if(s1 /= s2) then
+            if(monoBdo) then
+              bannedOrb(h1,s1) = .false.
+            end if
+            if(monoAdo) then
+              bannedOrb(h2,s2) = .false.
+              monoAdo = .false.
+            end if
+          end if
+
+          maskInd += 1
+          if(subset == 0 .or. mod(maskInd, fragment_count) == (subset-1)) then  
+            
+            call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
+            if(fullMatch) cycle
+          
+            mat = 0d0
+            call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
+
+            call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf)
+          end if
+        enddo
+        if(s1 /= s2) monoBdo = .false.
+      enddo
+    enddo
+  enddo
+end 
+
+
+
+subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in) :: i_generator, sp, h1, h2
+  double precision, intent(in) :: mat(N_states, mo_tot_num, mo_tot_num)
+  logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num)
+  double precision, intent(in)           :: fock_diag_tmp(mo_tot_num)
+  double precision, intent(in)    :: E0(N_states)
+  double precision, intent(inout) :: pt2(N_states) 
+  type(selection_buffer), intent(inout) :: buf
+  logical :: ok
+  integer :: s1, s2, p1, p2, ib, j, istate
+  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
+  double precision :: e_pert, delta_E, val, Hii, max_e_pert,tmp
+  double precision, external :: diag_H_mat_elem_fock
+  
+  logical, external :: detEq
+  
+  
+  if(sp == 3) then
+    s1 = 1
+    s2 = 2
+  else
+    s1 = sp
+    s2 = sp
+  end if
+  
+  call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
+  
+  do p1=1,mo_tot_num
+    if(bannedOrb(p1, s1)) cycle
+    ib = 1
+    if(sp /= 3) ib = p1+1
+    do p2=ib,mo_tot_num
+      if(bannedOrb(p2, s2)) cycle
+      if(banned(p1,p2)) cycle
+      if(mat(1, p1, p2) == 0d0) cycle
+      call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
+      
+      Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
+      max_e_pert = 0d0
+      
+      do istate=1,N_states
+        delta_E = E0(istate) - Hii
+        val = mat(istate, p1, p2) + mat(istate, p1, p2) 
+        tmp = dsqrt(delta_E * delta_E + val * val)
+        if (delta_E < 0.d0) then
+            tmp = -tmp
+        endif
+        e_pert = 0.5d0 * ( tmp - delta_E)
+        pt2(istate) = pt2(istate) + e_pert
+        max_e_pert = min(e_pert,max_e_pert)
+!        ci(istate) = e_pert / mat(istate, p1, p2)
+      end do
+      
+      if(dabs(max_e_pert) > buf%mini) then
+        call add_to_selection_buffer(buf, det, max_e_pert)
+      end if
+    end do
+  end do
+end 
+
+
+subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
+  use bitmasks
+  implicit none
+  
+  integer, intent(in)            :: sp, i_gen, N_sel
+  integer, intent(in)            :: interesting(0:N_sel)
+  integer(bit_kind),intent(in)   :: mask(N_int, 2), det(N_int, 2, N_sel)
+  logical, intent(inout)         :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+
+  integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt
+  integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
+!   logical :: bandon
+!   
+!   bandon = .false.
+  PROVIDE psi_phasemask psi_selectors_coef_transp
+  mat = 0d0
+  
+  do i=1,N_int
+    negMask(i,1) = not(mask(i,1))
+    negMask(i,2) = not(mask(i,2))
+  end do
+
+  do i=1, N_sel ! interesting(0)
+    !i = interesting(ii)
+    if (interesting(i) < 0) then
+      stop 'prefetch interesting(i)'
+    endif
+
+    
+    mobMask(1,1) = iand(negMask(1,1), det(1,1,i))
+    mobMask(1,2) = iand(negMask(1,2), det(1,2,i))
+    nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
+
+    if(nt > 4) cycle
+
+    do j=2,N_int
+      mobMask(j,1) = iand(negMask(j,1), det(j,1,i))
+      mobMask(j,2) = iand(negMask(j,2), det(j,2,i))
+      nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
+    end do
+
+    if(nt > 4) cycle
+    
+    if (interesting(i) == i_gen) then
+        if(sp == 3) then
+          do j=1,mo_tot_num
+            do k=1,mo_tot_num
+              banned(j,k,2) = banned(k,j,1)
+            enddo
+          enddo
+        else
+          do k=1,mo_tot_num
+          do l=k+1,mo_tot_num
+            banned(l,k,1) = banned(k,l,1)
+          end do
+          end do
+        end if
+    end if
+
+    call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
+    call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
+    
+    perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
+    perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
+    do j=2,N_int
+      perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
+      perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
+    end do
+
+    call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
+    call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
+
+    if (interesting(i) >= i_gen) then
+        if(nt == 4) then
+          call get_d2(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        else if(nt == 3) then
+          call get_d1(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        else
+          call get_d0(det(1,1,i), psi_phasemask(1,1,interesting(i)), bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
+        end if
+    else 
+        if(nt == 4) call past_d2(banned, p, sp)
+        if(nt == 3) call past_d1(bannedOrb, p)
+    end if
+  end do
+end 
+
+
+subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  double precision, external :: get_phase_bi, integral8
+  
+  integer :: i, j, tip, ma, mi, puti, putj
+  integer :: h1, h2, p1, p2, i1, i2
+  double precision :: hij, phase
+  
+  integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0,  3,4, 0,0, 0,0,  2,4, 1,4, 0,0,  2,3, 1,3, 1,2 /), (/2,3,4/))
+  integer, parameter :: turn2(2) = (/2, 1/)
+  integer, parameter :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+  
+  integer :: bant
+  bant = 1
+
+  tip = p(0,1) * p(0,2)
+  
+  ma = sp
+  if(p(0,1) > p(0,2)) ma = 1
+  if(p(0,1) < p(0,2)) ma = 2
+  mi = mod(ma, 2) + 1
+  
+  if(sp == 3) then
+    if(ma == 2) bant = 2
+    
+    if(tip == 3) then
+      puti = p(1, mi)
+      do i = 1, 3
+        putj = p(i, ma)
+        if(banned(putj,puti,bant)) cycle
+        i1 = turn3(1,i)
+        i2 = turn3(2,i)
+        p1 = p(i1, ma)
+        p2 = p(i2, ma)
+        h1 = h(1, ma)
+        h2 = h(2, ma)
+        
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
+        if(ma == 1) then
+          mat(:, putj, puti) += coefs * hij
+        else
+          mat(:, puti, putj) += coefs * hij
+        end if
+      end do
+    else
+      h1 = h(1,1)
+      h2 = h(1,2)
+      do j = 1,2
+        putj = p(j, 2)
+        p2 = p(turn2(j), 2)
+        do i = 1,2
+          puti = p(i, 1)
+          
+          if(banned(puti,putj,bant)) cycle
+          p1 = p(turn2(i), 1)
+          
+          hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          mat(:, puti, putj) += coefs * hij
+        end do
+      end do
+    end if
+
+  else
+    if(tip == 0) then
+      h1 = h(1, ma)
+      h2 = h(2, ma)
+      do i=1,3
+      puti = p(i, ma)
+      do j=i+1,4
+        putj = p(j, ma)
+        if(banned(puti,putj,1)) cycle
+        
+        i1 = turn2d(1, i, j)
+        i2 = turn2d(2, i, j)
+        p1 = p(i1, ma)
+        p2 = p(i2, ma)
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2)
+        mat(:, puti, putj) += coefs * hij
+      end do
+      end do
+    else if(tip == 3) then
+      h1 = h(1, mi)
+      h2 = h(1, ma)
+      p1 = p(1, mi)
+      do i=1,3
+        puti = p(turn3(1,i), ma)
+        putj = p(turn3(2,i), ma)
+        if(banned(puti,putj,1)) cycle
+        p2 = p(i, ma)
+        
+        hij = integral8(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2)
+        mat(:, min(puti, putj), max(puti, putj)) += coefs * hij
+      end do
+    else ! tip == 4
+      puti = p(1, sp)
+      putj = p(2, sp)
+      if(.not. banned(puti,putj,1)) then
+        p1 = p(1, mi)
+        p2 = p(2, mi)
+        h1 = h(1, mi)
+        h2 = h(2, mi)
+        hij = (integral8(p1, p2, h1, h2) - integral8(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2)
+        mat(:, puti, putj) += coefs * hij
+      end if
+    end if
+  end if
+end 
+
+
+subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
+  integer(1),intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+  double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num)
+  double precision, external :: get_phase_bi, integral8
+  
+  logical :: lbanned(mo_tot_num, 2), ok
+  integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j, hfix, pfix, h1, h2, p1, p2, ib
+  
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  integer, parameter :: turn2(2) = (/2,1/)
+  integer, parameter :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+  
+  integer :: bant
+  
+  
+  lbanned = bannedOrb
+    
+  do i=1, p(0,1)
+    lbanned(p(i,1), 1) = .true.
+  end do
+  do i=1, p(0,2)
+    lbanned(p(i,2), 2) = .true.
+  end do
+  
+  ma = 1
+  if(p(0,2) >= 2) ma = 2
+  mi = turn2(ma)
+  
+  bant = 1
+
+  if(sp == 3) then
+    !move MA
+    if(ma == 2) bant = 2
+    puti = p(1,mi)
+    hfix = h(1,ma)
+    p1 = p(1,ma)
+    p2 = p(2,ma)
+    if(.not. bannedOrb(puti, mi)) then
+      tmp_row = 0d0
+      do putj=1, hfix-1
+        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
+        hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+      end do
+      do putj=hfix+1, mo_tot_num
+        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
+        hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+      end do
+
+      if(ma == 1) then           
+        mat(1:N_states,1:mo_tot_num,puti) += tmp_row(1:N_states,1:mo_tot_num)
+      else
+        mat(1:N_states,puti,1:mo_tot_num) += tmp_row(1:N_states,1:mo_tot_num)
+      end if
+    end if
+
+    !MOVE MI
+    pfix = p(1,mi)
+    tmp_row = 0d0
+    tmp_row2 = 0d0
+    do puti=1,mo_tot_num
+      if(lbanned(puti,mi)) cycle
+      !p1 fixed
+      putj = p1
+      if(.not. banned(putj,puti,bant)) then
+        hij = integral8(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
+        tmp_row(:,puti) += hij * coefs
+      end if
+      
+      putj = p2
+      if(.not. banned(putj,puti,bant)) then
+        hij = integral8(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
+        tmp_row2(:,puti) += hij * coefs
+      end if
+    end do
+    
+    if(mi == 1) then
+      mat(:,:,p1) += tmp_row(:,:)
+      mat(:,:,p2) += tmp_row2(:,:)
+    else
+      mat(:,p1,:) += tmp_row(:,:)
+      mat(:,p2,:) += tmp_row2(:,:)
+    end if
+  else
+    if(p(0,ma) == 3) then
+      do i=1,3
+        hfix = h(1,ma)
+        puti = p(i, ma)
+        p1 = p(turn3(1,i), ma)
+        p2 = p(turn3(2,i), ma)
+        tmp_row = 0d0
+        do putj=1,hfix-1
+          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
+          hij = (integral8(p1, p2, putj, hfix)-integral8(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+          tmp_row(:,putj) += hij * coefs
+        end do
+        do putj=hfix+1,mo_tot_num
+          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
+          hij = (integral8(p1, p2, hfix, putj)-integral8(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+          tmp_row(:,putj) += hij * coefs
+        end do
+
+        mat(:, :puti-1, puti) += tmp_row(:,:puti-1)
+        mat(:, puti, puti:) += tmp_row(:,puti:)
+      end do
+    else
+      hfix = h(1,mi)
+      pfix = p(1,mi)
+      p1 = p(1,ma)
+      p2 = p(2,ma)
+      tmp_row = 0d0
+      tmp_row2 = 0d0
+      do puti=1,mo_tot_num
+        if(lbanned(puti,ma)) cycle
+        putj = p2
+        if(.not. banned(puti,putj,1)) then
+          hij = integral8(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
+          tmp_row(:,puti) += hij * coefs
+        end if
+        
+        putj = p1
+        if(.not. banned(puti,putj,1)) then
+          hij = integral8(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
+          tmp_row2(:,puti) += hij * coefs
+        end if
+      end do
+      mat(:,:p2-1,p2) += tmp_row(:,:p2-1)
+      mat(:,p2,p2:) += tmp_row(:,p2:)
+      mat(:,:p1-1,p1) += tmp_row2(:,:p1-1)
+      mat(:,p1,p1:) += tmp_row2(:,p1:)
+    end if
+  end if
+
+ !! MONO
+    if(sp == 3) then
+      s1 = 1
+      s2 = 2
+    else
+      s1 = sp
+      s2 = sp
+    end if
+
+    do i1=1,p(0,s1)
+      ib = 1
+      if(s1 == s2) ib = i1+1
+      do i2=ib,p(0,s2)
+        p1 = p(i1,s1)
+        p2 = p(i2,s2)
+        if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
+        call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
+        call i_h_j(gen, det, N_int, hij)
+        mat(:, p1, p2) += coefs * hij
+      end do
+    end do
+end 
+
+
+
+
+subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
+  integer(1), intent(in) :: phasemask(2,N_int*bit_kind_size)
+  logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision, intent(in) :: coefs(N_states)
+  double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+  
+  integer :: i, j, s, h1, h2, p1, p2, puti, putj
+  double precision :: hij, phase
+  double precision, external :: get_phase_bi, integral8
+  logical :: ok
+  
+  integer :: bant
+  bant = 1
+  
+
+  if(sp == 3) then ! AB
+    h1 = p(1,1)
+    h2 = p(1,2)
+    do p1=1, mo_tot_num
+      if(bannedOrb(p1, 1)) cycle
+      do p2=1, mo_tot_num
+        if(bannedOrb(p2,2)) cycle
+        if(banned(p1, p2, bant)) cycle ! rentable?
+        if(p1 == h1 .or. p2 == h2) then
+          call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
+          call i_h_j(gen, det, N_int, hij)
+        else
+          phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
+          hij = integral8(p1, p2, h1, h2) * phase
+        end if
+        mat(:, p1, p2) += coefs(:) * hij
+      end do
+    end do
+  else ! AA BB
+    p1 = p(1,sp)
+    p2 = p(2,sp)
+    do puti=1, mo_tot_num
+      if(bannedOrb(puti, sp)) cycle
+      do putj=puti+1, mo_tot_num
+        if(bannedOrb(putj, sp)) cycle
+        if(banned(puti, putj, bant)) cycle ! rentable?
+        if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
+          call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
+          call i_h_j(gen, det, N_int, hij)
+        else
+          hij = (integral8(p1, p2, puti, putj) -  integral8(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
+        end if
+        mat(:, puti, putj) += coefs(:) * hij
+      end do
+    end do
+  end if
+end 
+ 
+
+subroutine past_d1(bannedOrb, p)
+  use bitmasks
+  implicit none
+
+  logical, intent(inout) :: bannedOrb(mo_tot_num, 2)
+  integer, intent(in) :: p(0:4, 2)
+  integer :: i,s
+
+  do s = 1, 2
+    do i = 1, p(0, s)
+      bannedOrb(p(i, s), s) = .true.
+    end do
+  end do
+end 
+
+
+subroutine past_d2(banned, p, sp)
+  use bitmasks
+  implicit none
+
+  logical, intent(inout) :: banned(mo_tot_num, mo_tot_num)
+  integer, intent(in) :: p(0:4, 2), sp
+  integer :: i,j
+
+  if(sp == 3) then
+    do i=1,p(0,1)
+      do j=1,p(0,2)
+        banned(p(i,1), p(j,2)) = .true.
+      end do
+    end do
+  else
+    do i=1,p(0, sp)
+      do j=1,i-1
+        banned(p(j,sp), p(i,sp)) = .true.
+        banned(p(i,sp), p(j,sp)) = .true.
+      end do
+    end do
+  end if
+end 
+
+
+
+subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
+  use bitmasks
+  implicit none
+  
+  integer, intent(in) :: i_gen, N
+  integer, intent(in) :: interesting(0:N)
+  integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N)
+  logical, intent(inout) :: banned(mo_tot_num, mo_tot_num)
+  logical, intent(out) :: fullMatch
+
+
+  integer :: i, j, na, nb, list(3)
+  integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2)
+
+  fullMatch = .false.
+
+  do i=1,N_int
+    negMask(i,1) = not(mask(i,1))
+    negMask(i,2) = not(mask(i,2))
+  end do
+
+  genl : do i=1, N
+    do j=1, N_int
+      if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl
+      if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
+    end do
+
+    if(interesting(i) < i_gen) then
+      fullMatch = .true.
+      return
+    end if
+
+    do j=1, N_int
+      myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1))
+      myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2))
+    end do
+
+    call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int)
+    call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int)
+    banned(list(1), list(2)) = .true.
+  end do genl
+end 
+
+
+subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Gives the inidices(+1) of the bits set to 1 in the bit string
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: string(Nint)
+  integer, intent(out)           :: list(Nint*bit_kind_size)
+  integer, intent(out)           :: n_elements
+  
+  integer                        :: i, ishift
+  integer(bit_kind)              :: l
+  
+  n_elements = 0
+  ishift = 2
+  do i=1,Nint
+    l = string(i)
+    do while (l /= 0_bit_kind)
+      n_elements = n_elements+1
+      list(n_elements) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
+      l = iand(l,l-1_bit_kind)
+    enddo
+    ishift = ishift + bit_kind_size
+  enddo
+  
+end
diff --git a/plugins/Full_CI_ZMQ/selection_buffer.irp.f b/plugins/Full_CI_ZMQ/selection_buffer.irp.f
index d296b399..8a47cb9d 100644
--- a/plugins/Full_CI_ZMQ/selection_buffer.irp.f
+++ b/plugins/Full_CI_ZMQ/selection_buffer.irp.f
@@ -41,31 +41,33 @@ subroutine sort_selection_buffer(b)
   implicit none
 
   type(selection_buffer), intent(inout) :: b
-  double precision, allocatable :: vals(:), absval(:)
+  double precision, allocatable::  absval(:)
   integer, allocatable :: iorder(:)
-  integer(bit_kind), allocatable :: detmp(:,:,:)
+  double precision, pointer :: vals(:)
+  integer(bit_kind), pointer :: detmp(:,:,:)
   integer :: i, nmwen
   logical, external :: detEq
   nmwen = min(b%N, b%cur)
 
 
-  allocate(iorder(b%cur), detmp(N_int, 2, nmwen), absval(b%cur), vals(nmwen))
+  allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
   absval = -dabs(b%val(:b%cur))
   do i=1,b%cur
     iorder(i) = i
   end do
-  call dsort(absval, iorder, b%cur)
-
+  ! Optimal for almost sorted data
+  call insertion_dsort(absval, iorder, b%cur)
   do i=1, nmwen
     detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
     detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
     vals(i) = b%val(iorder(i))
   end do
-  b%det = 0_bit_kind
-  b%val = 0d0
-  b%det(1:N_int,1,1:nmwen) = detmp(1:N_int,1,1:nmwen)
-  b%det(1:N_int,2,1:nmwen) = detmp(1:N_int,2,1:nmwen)
-  b%val(1:nmwen) = vals(1:nmwen)
+  do i=nmwen+1, size(vals)
+    vals(i) = 0.d0
+  enddo
+  deallocate(b%det, b%val)
+  b%det => detmp
+  b%val => vals
   b%mini = max(b%mini,dabs(b%val(b%N)))
   b%cur = nmwen
 end subroutine
diff --git a/plugins/Full_CI_ZMQ/selection_types.f90 b/plugins/Full_CI_ZMQ/selection_types.f90
index 9506629c..29e48524 100644
--- a/plugins/Full_CI_ZMQ/selection_types.f90
+++ b/plugins/Full_CI_ZMQ/selection_types.f90
@@ -1,9 +1,9 @@
 module selection_types
   type selection_buffer
     integer :: N, cur
-    integer(8), allocatable :: det(:,:,:)
-    double precision, allocatable :: val(:)
-    double precision :: mini
+    integer(8)      , pointer :: det(:,:,:)
+    double precision, pointer :: val(:)
+    double precision          :: mini
   endtype
 end module
 
diff --git a/plugins/Full_CI_ZMQ/zmq_selection.irp.f b/plugins/Full_CI_ZMQ/zmq_selection.irp.f
index 8aaddc19..62703a43 100644
--- a/plugins/Full_CI_ZMQ/zmq_selection.irp.f
+++ b/plugins/Full_CI_ZMQ/zmq_selection.irp.f
@@ -10,26 +10,38 @@ subroutine ZMQ_selection(N_in, pt2)
   integer                        :: i, N
   integer, external              :: omp_get_thread_num
   double precision, intent(out)  :: pt2(N_states)
+  integer, parameter             :: maxtasks=10000
   
   
   PROVIDE fragment_count
 
+  N = max(N_in,1)
   if (.True.) then
     PROVIDE pt2_e0_denominator
-    N = max(N_in,1)
     provide nproc
     call new_parallel_job(zmq_to_qp_run_socket,"selection")
     call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
-    call zmq_set_running(zmq_to_qp_run_socket)
     call create_selection_buffer(N, N*2, b)
   endif
 
-  character(len=:), allocatable  :: task 
-  task = repeat(' ',20*N_det_generators)
+  character*(20*maxtasks) :: task
+  task = ' ' 
+
+  integer :: k
+  k=0
   do i= 1, N_det_generators
-    write(task(20*(i-1)+1:20*i),'(I9,X,I9,''|'')') i, N
+    k = k+1
+    write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
+    k = k+20
+    if (k>20*maxtasks) then
+       k=0
+       call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+    endif
   end do
-  call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+  if (k > 0) then
+    call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+  endif
+  call zmq_set_running(zmq_to_qp_run_socket)
 
   !$OMP PARALLEL DEFAULT(shared)  SHARED(b, pt2)  PRIVATE(i) NUM_THREADS(nproc+1)
   i = omp_get_thread_num()
@@ -48,6 +60,7 @@ subroutine ZMQ_selection(N_in, pt2)
     endif
     call save_wavefunction
   endif
+
 end subroutine
 
 
@@ -83,7 +96,7 @@ subroutine selection_collector(b, pt2)
   real :: time, time0
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_pull = new_zmq_pull_socket()
-  allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
+  allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators))
   done = 0
   more = 1
   pt2(:) = 0d0
diff --git a/plugins/MRCC_Utils/amplitudes.irp.f b/plugins/MRCC_Utils/amplitudes.irp.f
index f9cb51ad..ccbe700d 100644
--- a/plugins/MRCC_Utils/amplitudes.irp.f
+++ b/plugins/MRCC_Utils/amplitudes.irp.f
@@ -23,33 +23,39 @@
  allocate(pathTo(N_det_non_ref))
 
  pathTo(:) = 0
- is_active_exc(:) = .false.
+ is_active_exc(:) = .True. 
  n_exc_active = 0
 
- do hh = 1, hh_shortcut(0)
-   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
-     do II = 1, N_det_ref
+! do hh = 1, hh_shortcut(0)
+!   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
+!     do II = 1, N_det_ref
+!
+!       call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
+!       if(.not. ok) cycle
+!
+!       call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
+!       if(.not. ok) cycle
+!
+!       ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
+!       if(ind == -1) cycle
+!
+!       logical, external :: is_a_two_holes_two_particles
+!       if (is_a_two_holes_two_particles(myDet)) then
+!         is_active_exc(pp) = .False.
+!       endif
 
-       call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
-       if(.not. ok) cycle
+!       ind = psi_non_ref_sorted_idx(ind)
+!       if(pathTo(ind) == 0) then
+!         pathTo(ind) = pp
+!       else
+!         is_active_exc(pp) = .true.
+!         is_active_exc(pathTo(ind)) = .true.
+!       end if
 
-       call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
-       if(.not. ok) cycle
+!     end do
+!   end do
+! end do
 
-       ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
-       if(ind == -1) cycle
-
-       ind = psi_non_ref_sorted_idx(ind)
-       if(pathTo(ind) == 0) then
-         pathTo(ind) = pp
-       else
-         is_active_exc(pp) = .true.
-         is_active_exc(pathTo(ind)) = .true.
-       end if
-     end do
-   end do
- end do
-!is_active_exc=.true.
  do hh = 1, hh_shortcut(0)
   do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
     if(is_active_exc(pp)) then
@@ -66,6 +72,32 @@
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ logical, has_a_unique_parent, (N_det_non_ref) ]
+ implicit none
+ BEGIN_DOC
+ ! True if the determinant in the non-reference has a unique parent
+ END_DOC
+ integer :: i,j,n
+ integer :: degree
+ do j=1,N_det_non_ref
+  has_a_unique_parent(j) = .True. 
+  n=0
+  do i=1,N_det_ref
+   call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,j), degree, N_int)
+   if (degree < 2) then
+     n = n+1
+     if (n > 1) then
+       has_a_unique_parent(j) = .False.
+       exit
+     endif
+   endif
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+
 BEGIN_PROVIDER [ integer, n_exc_active_sze ]
  implicit none
  BEGIN_DOC
@@ -96,7 +128,7 @@ END_PROVIDER
      !$OMP   active_excitation_to_determinants_val, active_excitation_to_determinants_idx)&
      !$OMP shared(hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, &
      !$OMP   psi_non_ref_sorted_idx, psi_ref, N_det_ref, N_states)&
-     !$OMP shared(is_active_exc, active_hh_idx, active_pp_idx, n_exc_active)&
+     !$OMP shared(active_hh_idx, active_pp_idx, n_exc_active)&
      !$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh, s)
  allocate(lref(N_det_non_ref))
  !$OMP DO schedule(dynamic)
diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index 69e5a662..7ba210ca 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -351,11 +351,11 @@ logical function is_generable(det1, det2, Nint)
   integer, intent(in) :: Nint
   integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
   integer :: degree, f, exc(0:2, 2, 2), t
-  integer*2 :: h1, h2, p1, p2, s1, s2
+  integer :: h1, h2, p1, p2, s1, s2
   integer, external :: searchExc
   logical, external :: excEq
   double precision :: phase
-  integer*2 :: tmp_array(4)
+  integer :: tmp_array(4)
   
   is_generable = .false.
   call get_excitation(det1, det2, exc, degree, phase, Nint)
@@ -366,7 +366,7 @@ logical function is_generable(det1, det2, Nint)
   end if
   if(degree > 2) stop "?22??"
   
-  call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
   
   if(degree == 1) then
     h2 = h1
@@ -454,7 +454,7 @@ integer function searchExc(excs, exc, n)
   use bitmasks
   
   integer, intent(in) :: n
-  integer*2,intent(in) :: excs(4,n), exc(4)
+  integer,intent(in) :: excs(4,n), exc(4)
   integer :: l, h, c
   integer, external :: excCmp
   logical, external :: excEq
@@ -519,8 +519,8 @@ subroutine sort_exc(key, N_key)
   
 
   integer, intent(in)                   :: N_key
-  integer*2,intent(inout)       :: key(4,N_key)
-  integer*2                     :: tmp(4)
+  integer,intent(inout)       :: key(4,N_key)
+  integer                     :: tmp(4)
   integer                               :: i,ni
   
   
@@ -542,7 +542,7 @@ end subroutine
 
 logical function exc_inf(exc1, exc2)
   implicit none
-  integer*2,intent(in) :: exc1(4), exc2(4)
+  integer,intent(in) :: exc1(4), exc2(4)
   integer :: i
   exc_inf = .false.
   do i=1,4
@@ -564,9 +564,9 @@ subroutine tamise_exc(key, no, n, N_key)
 ! Uncodumented : TODO
   END_DOC
   integer,intent(in)            :: no, n, N_key
-  integer*2,intent(inout)       :: key(4, N_key)
+  integer,intent(inout)       :: key(4, N_key)
   integer                       :: k,j
-  integer*2                     :: tmp(4)
+  integer                     :: tmp(4)
   logical                       :: exc_inf
   integer                       :: ni
   
@@ -595,8 +595,9 @@ end subroutine
 
 subroutine dec_exc(exc, h1, h2, p1, p2)
   implicit none
-  integer :: exc(0:2,2,2), s1, s2, degree
-  integer*2, intent(out) :: h1, h2, p1, p2
+  integer, intent(in)            :: exc(0:2,2,2)
+  integer, intent(out)           :: h1, h2, p1, p2
+  integer                        :: degree, s1, s2
   
   degree = exc(0,1,1) + exc(0,1,2)
   
@@ -607,7 +608,7 @@ subroutine dec_exc(exc, h1, h2, p1, p2)
     
   if(degree == 0) return
   
-  call decode_exc_int2(exc, degree, h1, p1, h2, p2, s1, s2)
+  call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
   
   h1 += mo_tot_num * (s1-1)
   p1 += mo_tot_num * (s1-1)
@@ -639,7 +640,7 @@ end subroutine
 &BEGIN_PROVIDER [ integer, N_ex_exists ]
   implicit none
   integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
-  integer*2 :: h1, h2, p1, p2
+  integer :: h1, h2, p1, p2
   double precision :: phase
   logical,allocatable :: hh(:,:) , pp(:,:)
   
@@ -739,12 +740,12 @@ END_PROVIDER
   double precision               :: phase
   
   
-  double precision, allocatable :: rho_mrcc_init(:)
+  double precision, allocatable :: rho_mrcc_inact(:)
   integer                        :: a_coll, at_roww
   
   print *, "TI", hh_nex, N_det_non_ref
 
-  allocate(rho_mrcc_init(N_det_non_ref))
+  allocate(rho_mrcc_inact(N_det_non_ref))
   allocate(x_new(hh_nex))
   allocate(x(hh_nex), AtB(hh_nex))
 
@@ -756,7 +757,7 @@ END_PROVIDER
         !$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
         !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx)
     
-    !$OMP DO schedule(dynamic, 100)
+    !$OMP DO schedule(static, 100)
     do at_roww = 1, n_exc_active ! hh_nex
       at_row = active_pp_idx(at_roww)
       do i=1,active_excitation_to_determinants_idx(0,at_roww)
@@ -775,7 +776,7 @@ END_PROVIDER
       X(a_col) = AtB(a_col)
     end do
     
-    rho_mrcc_init = 0d0
+    rho_mrcc_inact(:) = 0d0
     
     allocate(lref(N_det_ref))
     do hh = 1, hh_shortcut(0)
@@ -799,19 +800,15 @@ END_PROVIDER
         X(pp) =  AtB(pp) 
         do II=1,N_det_ref
           if(lref(II) > 0) then
-            rho_mrcc_init(lref(II)) = psi_ref_coef(II,s) * X(pp)
+            rho_mrcc_inact(lref(II)) = psi_ref_coef(II,s) * X(pp)
           else if(lref(II) < 0) then
-            rho_mrcc_init(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
+            rho_mrcc_inact(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
           end if
         end do
       end do
     end do
     deallocate(lref)
 
-    do i=1,N_det_non_ref
-      rho_mrcc(i,s) = rho_mrcc_init(i) 
-    enddo
-
     x_new = x
     
     double precision               :: factor, resold
@@ -839,7 +836,10 @@ END_PROVIDER
         print *, k, res, 1.d0 - res/resold
       endif
       
-      if ( (res < 1d-10).or.(res/resold > 0.99d0) ) then
+      if ( res < 1d-10 ) then
+        exit
+      endif
+      if ( (res/resold > 0.99d0) ) then
         exit
       endif
       resold = res
@@ -848,38 +848,60 @@ END_PROVIDER
     dIj_unique(1:size(X), s) = X(1:size(X))
     print *, k, res, 1.d0 - res/resold
 
-  enddo
 
-  do s=1,N_states
+    do i=1,N_det_non_ref
+      rho_mrcc(i,s) = 0.d0
+    enddo
 
     do a_coll=1,n_exc_active
       a_col = active_pp_idx(a_coll)
       do j=1,N_det_non_ref
         i = active_excitation_to_determinants_idx(j,a_coll)
         if (i==0) exit
+        if (rho_mrcc_inact(i) /= 0.d0) then
+          call debug_det(psi_non_ref(1,1,i),N_int)
+          stop
+        endif
         rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s)
       enddo
     end do
 
-    norm = 0.d0
-    do i=1,N_det_non_ref
-      norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
-    enddo
-    ! Norm now contains the norm of A.X
-    
+    double precision :: norm2_ref, norm2_inact, a, b, c, Delta
+    ! Psi = Psi_ref + Psi_inactive + f*Psi_active
+    ! Find f to normalize Psi
+
+    norm2_ref = 0.d0
     do i=1,N_det_ref
-      norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
+      norm2_ref = norm2_ref + psi_ref_coef(i,s)*psi_ref_coef(i,s)
     enddo
-    ! Norm now contains the norm of Psi + A.X
-    
+
+    a = 0.d0
+    do i=1,N_det_non_ref
+      a = a + rho_mrcc(i,s)*rho_mrcc(i,s)
+    enddo
+
+    norm = a + norm2_ref
     print *, "norm : ", sqrt(norm)
-   enddo
+
+    norm = sqrt((1.d0-norm2_ref)/a)
+
+    ! Renormalize Psi+A.X
+    do i=1,N_det_non_ref
+      rho_mrcc(i,s) = rho_mrcc(i,s) * norm 
+    enddo
+
+!norm = norm2_ref
+!do i=1,N_det_non_ref
+!  norm = norm + rho_mrcc(i,s)**2
+!enddo
+!print *,  'check', norm
+!stop
+
      
         
-   do s=1,N_states
      norm = 0.d0
      double precision               :: f, g, gmax
-     gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s)))
+     gmax = maxval(dabs(psi_non_ref_coef(:,s)))
      do i=1,N_det_non_ref
        if (lambda_type == 2) then
          f = 1.d0
@@ -891,41 +913,22 @@ END_PROVIDER
         f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
 
         ! Avoid numerical instabilities
-!        g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) )
         g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax))
         f = min(f, g)
         f = max(f,-g)
+
       endif
 
-       norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
+       norm = norm + (rho_mrcc(i,s)*f)**2
        rho_mrcc(i,s) = f
      enddo
-     ! norm now contains the norm of |T.Psi_0>
-     ! rho_mrcc now contains the f factors
-
-     f = 1.d0/norm
-     ! f now contains 1/ <T.Psi_0|T.Psi_0>
-
-     norm = 0.d0
-     do i=1,N_det_non_ref
-       norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
-     enddo
-     ! norm now contains <Psi_SD|Psi_SD>
-     f = dsqrt(f*norm)
-     ! f normalises T.Psi_0 such that (1+T)|Psi> is normalized
+     ! rho_mrcc now contains the mu_i factors
 
      print *,  'norm of |T Psi_0> = ', dsqrt(norm)
-     norm = norm*f
-     if (dsqrt(norm) > 1.d0) then
+     if (norm > 1.d0) then
        stop 'Error : Norm of the SD larger than the norm of the reference.'
      endif
 
-     do i=1,N_det_non_ref
-       rho_mrcc(i,s) = rho_mrcc(i,s) * f
-     enddo
-     ! rho_mrcc now contains the product of the scaling factors and the
-     ! normalization constant
-    
   end do
 
 END_PROVIDER
@@ -1028,11 +1031,11 @@ double precision function get_dij(det1, det2, s, Nint)
   integer, intent(in) :: s, Nint
   integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
   integer :: degree, f, exc(0:2, 2, 2), t
-  integer*2 :: h1, h2, p1, p2, s1, s2
+  integer :: h1, h2, p1, p2, s1, s2
   integer, external :: searchExc
   logical, external :: excEq
   double precision :: phase
-  integer*2 :: tmp_array(4)
+  integer :: tmp_array(4)
   
   get_dij = 0d0
   call get_excitation(det1, det2, exc, degree, phase, Nint)
@@ -1041,7 +1044,7 @@ double precision function get_dij(det1, det2, s, Nint)
     stop "get_dij"
   end if
   
-  call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
   
   if(degree == 1) then
     h2 = h1
@@ -1074,8 +1077,8 @@ double precision function get_dij(det1, det2, s, Nint)
 end function
 
 
- BEGIN_PROVIDER [ integer*2, hh_exists, (4, N_hh_exists) ]
-&BEGIN_PROVIDER [ integer*2, pp_exists, (4, N_pp_exists) ]
+ BEGIN_PROVIDER [ integer, hh_exists, (4, N_hh_exists) ]
+&BEGIN_PROVIDER [ integer, pp_exists, (4, N_pp_exists) ]
 &BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_hh_exists + 1) ]
 &BEGIN_PROVIDER [ integer, hh_nex ]
   implicit none
@@ -1090,9 +1093,9 @@ end function
   ! hh_nex : Total number of excitation operators
   !
   END_DOC
-  integer*2,allocatable :: num(:,:)
+  integer,allocatable :: num(:,:)
   integer :: exc(0:2, 2, 2), degree, n, on, s, l, i
-  integer*2 :: h1, h2, p1, p2
+  integer :: h1, h2, p1, p2
   double precision :: phase
   logical, external :: excEq
   
@@ -1118,19 +1121,19 @@ end function
   
   hh_shortcut(0) = 1
   hh_shortcut(1) = 1
-  hh_exists(:,1) = (/1_2, num(1,1), 1_2, num(2,1)/)
-  pp_exists(:,1) = (/1_2, num(3,1), 1_2, num(4,1)/)
+  hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/)
+  pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/)
   s = 1
   do i=2,n
     if(.not. excEq(num(1,i), num(1,s))) then
       s += 1
       num(:, s) = num(:, i)
-      pp_exists(:,s) = (/1_2, num(3,s), 1_2, num(4,s)/)
+      pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/)
       if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. &
             hh_exists(4, hh_shortcut(0)) /= num(2,s)) then
         hh_shortcut(0) += 1
         hh_shortcut(hh_shortcut(0)) = s
-        hh_exists(:,hh_shortcut(0)) = (/1_2, num(1,s), 1_2, num(2,s)/)
+        hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/)
       end if
     end if
   end do
@@ -1178,7 +1181,7 @@ END_PROVIDER
 
 logical function excEq(exc1, exc2)
   implicit none
-  integer*2, intent(in) :: exc1(4), exc2(4)
+  integer, intent(in) :: exc1(4), exc2(4)
   integer :: i
   excEq = .false.
   do i=1, 4
@@ -1190,7 +1193,7 @@ end function
 
 integer function excCmp(exc1, exc2)
   implicit none
-  integer*2, intent(in) :: exc1(4), exc2(4)
+  integer, intent(in) :: exc1(4), exc2(4)
   integer :: i
   excCmp = 0
   do i=1, 4
@@ -1209,8 +1212,8 @@ subroutine apply_hole_local(det, exc, res, ok, Nint)
   use bitmasks
   implicit none
   integer, intent(in) :: Nint
-  integer*2, intent(in) :: exc(4)
-  integer*2 :: s1, s2, h1, h2
+  integer, intent(in) :: exc(4)
+  integer :: s1, s2, h1, h2
   integer(bit_kind),intent(in) :: det(Nint, 2)
   integer(bit_kind),intent(out) :: res(Nint, 2)
   logical, intent(out) :: ok
@@ -1246,8 +1249,8 @@ subroutine apply_particle_local(det, exc, res, ok, Nint)
   use bitmasks
   implicit none
   integer, intent(in) :: Nint
-  integer*2, intent(in) :: exc(4)
-  integer*2 :: s1, s2, p1, p2
+  integer, intent(in) :: exc(4)
+  integer :: s1, s2, p1, p2
   integer(bit_kind),intent(in) :: det(Nint, 2)
   integer(bit_kind),intent(out) :: res(Nint, 2)
   logical, intent(out) :: ok
diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f
index ef026e07..e8d19166 100644
--- a/plugins/MRPT_Utils/energies_cas.irp.f
+++ b/plugins/MRPT_Utils/energies_cas.irp.f
@@ -1167,11 +1167,11 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
        print*, 'e corr perturb EN',accu(state_target)
        print*, ''
        print*, 'coef diagonalized'
-       write(*,'(100(F16.10,X))')psi_in_out_coef(:,state_target)
+       write(*,'(100(F16.10,1X))')psi_in_out_coef(:,state_target)
        print*, 'coef_perturb'
-       write(*,'(100(F16.10,X))')coef_perturb(:)
+       write(*,'(100(F16.10,1X))')coef_perturb(:)
        print*, 'coef_perturb EN'
-       write(*,'(100(F16.10,X))')coef_perturb_bis(:)
+       write(*,'(100(F16.10,1X))')coef_perturb_bis(:)
       endif
       integer :: k
       do k = 1, N_det_ref
diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f
index 4042d90b..9376e0cc 100644
--- a/plugins/MRPT_Utils/excitations_cas.irp.f
+++ b/plugins/MRPT_Utils/excitations_cas.irp.f
@@ -22,7 +22,7 @@ subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, &
 
   integer :: elec_num_tab_local(2)
   integer :: i,j,accu_elec,k
-  integer :: det_tmp(N_int), det_tmp_bis(N_int)
+  integer(bit_kind) :: det_tmp(N_int), det_tmp_bis(N_int)
   double precision :: phase
   double precision :: norm_factor
 ! print*, orb,hole_particle,spin_exc
diff --git a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f
index 676e14e9..b67f7498 100644
--- a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f
+++ b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f
@@ -210,10 +210,6 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
 
                  ! < det_tmp     | H | det_tmp_bis >  = F_{aorb,borb}
                  hab = (fock_operator_local(aorb,borb,kspin) ) * phase  
-            if(isnan(hab))then
-             print*, '1'
-             stop
-            endif
                  ! < jdet | H | det_tmp_bis > = phase * (ir|cv)
                  call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
                  if(ispin == jspin)then
@@ -255,7 +251,8 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
                  call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
 !                ! < det_tmp | H | det_tmp_bis >  = F_{aorb,borb}
                  hab = fock_operator_local(aorb,borb,kspin) * phase  
-            if(isnan(hab))then
+!            if(isnan(hab))then
+            if(hab /= hab)then
              print*, '2'
              stop
             endif
diff --git a/plugins/Properties/delta_rho.irp.f b/plugins/Properties/delta_rho.irp.f
index 7803ba3d..8fd08246 100644
--- a/plugins/Properties/delta_rho.irp.f
+++ b/plugins/Properties/delta_rho.irp.f
@@ -6,7 +6,7 @@
  z_min = 0.d0
  z_max = 10.d0
  delta_z = 0.005d0
- N_z_pts = (z_max - z_min)/delta_z
+ N_z_pts = int( (z_max - z_min)/delta_z )
  print*,'N_z_pts = ',N_z_pts
 
 END_PROVIDER
diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f
index 6fa39278..91b26dc8 100644
--- a/plugins/Properties/hyperfine_constants.irp.f
+++ b/plugins/Properties/hyperfine_constants.irp.f
@@ -151,7 +151,7 @@ subroutine print_hcc
  integer :: i,j
  print*,'Z               AU           GAUSS              MHZ             cm^-1'
  do i = 1, nucl_num
-  write(*,'(I2,X,F4.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
+  write(*,'(I2,1X,F4.1,1X,4(F16.6,1X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
  enddo
 
 end
diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f
index deeb90bf..68b620c5 100644
--- a/plugins/Properties/mulliken.irp.f
+++ b/plugins/Properties/mulliken.irp.f
@@ -126,7 +126,7 @@ subroutine print_mulliken_sd
  accu = 0.d0
  do i = 1, ao_num
   accu += spin_gross_orbital_product(i)
-  write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
+  write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
  enddo
  print*,'sum = ',accu
  accu = 0.d0
@@ -142,7 +142,7 @@ subroutine print_mulliken_sd
   accu = 0.d0
   do i = 0,  ao_l_max
    accu += spin_population_angular_momentum_per_atom(i,j)
-   write(*,'(XX,I3,XX,A4,X,A4,X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
+   write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
    print*,'sum = ',accu
   enddo
  enddo
diff --git a/plugins/analyze_wf/analyze_wf.irp.f b/plugins/analyze_wf/analyze_wf.irp.f
index 6d8bffcf..7d005a05 100644
--- a/plugins/analyze_wf/analyze_wf.irp.f
+++ b/plugins/analyze_wf/analyze_wf.irp.f
@@ -18,7 +18,7 @@ subroutine run
   write(*,'(A)')  '============='
   write(*,'(A)')  ''
   do istate=1,N_states
-    call get_occupation_from_dets(occupation,1)
+    call get_occupation_from_dets(occupation,istate)
     write(*,'(A)')  ''
     write(*,'(A,I3)'),  'State ', istate
     write(*,'(A)')  '---------------'
@@ -30,13 +30,13 @@ subroutine run
       if (occupation(i) > 1.999d0) then
         class(0,1) += 1
         class( class(0,1), 1) = i
-      else if (occupation(i) > 1.95d0) then
+      else if (occupation(i) > 1.97d0) then
         class(0,2) += 1
         class( class(0,2), 2) = i
       else if (occupation(i) < 0.001d0) then
         class(0,5) += 1
         class( class(0,5), 5) = i
-      else if (occupation(i) < 0.01d0) then
+      else if (occupation(i) < 0.03d0) then
         class(0,4) += 1
         class( class(0,4), 4) = i
       else 
diff --git a/plugins/mrcc_selected/ezfio_interface.irp.f b/plugins/mrcc_selected/ezfio_interface.irp.f
index 47e7cea5..54d993fe 100644
--- a/plugins/mrcc_selected/ezfio_interface.irp.f
+++ b/plugins/mrcc_selected/ezfio_interface.irp.f
@@ -1,6 +1,6 @@
 ! DO NOT MODIFY BY HAND
 ! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
-! from file /ccc/work/cont003/gen1738/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
+! from file /panfs/panasas/cnt0024/cpq1738/scemama/workdir/quantum_package/src/mrcc_selected/EZFIO.cfg
 
 
 BEGIN_PROVIDER [ double precision, thresh_dressed_ci  ]
diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg
index a2dc1bb3..53519ec7 100644
--- a/plugins/mrcepa0/EZFIO.cfg
+++ b/plugins/mrcepa0/EZFIO.cfg
@@ -18,7 +18,7 @@ interface: ezfio
 type: logical
 doc: Compute perturbative contribution of the Triples
 interface: ezfio,provider,ocaml
-default: true
+default: false
 
 [energy]
 type: double precision
diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f
index 61443246..d2311676 100644
--- a/plugins/mrcepa0/dressing.irp.f
+++ b/plugins/mrcepa0/dressing.irp.f
@@ -38,7 +38,7 @@ use bitmasks
       do p=hh_shortcut(h), hh_shortcut(h+1)-1
         call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int)
         if(ok) n = n + 1
-        if(n > N_det_non_ref) stop "MRCC..."
+        if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
       end do
       n = n - 1
 
diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
index 6823df81..0c5e1b37 100755
--- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
+++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py
@@ -3,7 +3,7 @@
 convert output of gamess/GAU$$IAN to ezfio
 
 Usage:
-    qp_convert_output_to_ezfio.py <file.out> [--ezfio=<folder.ezfio>]
+    qp_convert_output_to_ezfio.py <file.out> [-o <ezfio_directory>]
 
 Option:
     file.out is the file to check (like gamess.out)
@@ -272,7 +272,7 @@ def write_ezfio(res, filename):
     #
 
     # INPUT
-    # {% for lanel,zcore, l_block in l_atom  $}
+    # {% for label,zcore, l_block in l_atom  $}
     #       #local l_block l=0}
     #       {label} GEN {zcore} {len(l_block)-1 #lmax_block}
     #       {% for l_param in l_block%}
@@ -280,6 +280,7 @@ def write_ezfio(res, filename):
     #                {% for coef,n,zeta for l_param}
     #                    {coef,n, zeta}
 
+
     # OUTPUT
 
     # Local are 1 array padded by max(n_max_block) when l == 0  (output:k_loc_max)
@@ -309,8 +310,16 @@ def write_ezfio(res, filename):
             array_l_max_block.append(l_max_block)
             array_z_remove.append(z_remove)
 
-            matrix.append([[coef_n_zeta.split()[1:] for coef_n_zeta in l.split('\n')] for l in array_party[1:]])
-
+            x = [[coef_n_zeta.split() for coef_n_zeta in l.split('\n')] \
+               for l in array_party[1:] ]
+            x = []
+            for l in array_party[1:]:
+              y = []
+              for coef_n_zeta in l.split('\n'):
+                z = coef_n_zeta.split()
+                if z : y.append(z)
+              x.append(y)
+            matrix.append(x)
         return (matrix, array_l_max_block, array_z_remove)
 
     def get_local_stuff(matrix):
@@ -319,7 +328,6 @@ def write_ezfio(res, filename):
         k_loc_max = max(len(i) for i in matrix_local_unpad)
 
         matrix_local = [ pad(ll, k_loc_max, [0., 2, 0.]) for ll in matrix_local_unpad]
-
         m_coef = [[float(i[0]) for i in atom] for atom in matrix_local]
         m_n = [[int(i[1]) - 2 for i in atom] for atom in matrix_local]
         m_zeta = [[float(i[2]) for i in atom] for atom in matrix_local]
@@ -343,9 +351,20 @@ def write_ezfio(res, filename):
         return (l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc)
 
     try:
-        pseudo_str = res_file.get_pseudo()
+        pseudo_str = []
+        label = ezfio.get_nuclei_nucl_label()
+        for ecp in res.pseudo:
+          pseudo_str += [ "%(label)s GEN %(zcore)d %(lmax)d" % { "label": label[ ecp["atom"]-1 ],
+                "zcore": ecp["zcore"], "lmax": ecp["lmax"] } ]
+          lmax = ecp["lmax"]
+          for l in [lmax] + list(range(0,lmax)):
+            pseudo_str += [ "%d"%len(ecp[str(l)]) ]
+            for t in ecp[str(l)]:
+              pseudo_str += [ "%f  %d  %f"%t ]
+          pseudo_str += [""]
+        pseudo_str = "\n".join(pseudo_str)
+             
         matrix, array_l_max_block, array_z_remove = parse_str(pseudo_str)
-
     except:
         ezfio.set_pseudo_do_pseudo(False)
     else:
@@ -359,10 +378,12 @@ def write_ezfio(res, filename):
         ezfio.nuclei_nucl_charge = [i - j for i, j in zip(ezfio.nuclei_nucl_charge, array_z_remove)]
 
         import math
-        num_elec = sum(ezfio.nuclei_nucl_charge)
+        num_elec_diff = sum(array_z_remove)/2
+        nalpha = ezfio.get_electrons_elec_alpha_num() - num_elec_diff
+        nbeta  = ezfio.get_electrons_elec_beta_num() - num_elec_diff
 
-        ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.))
-        ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.))
+        ezfio.set_electrons_elec_alpha_num(nalpha)
+        ezfio.set_electrons_elec_beta_num( nbeta )
 
         # Change all the array 'cause EZFIO
         #   v_kl (v, l) => v_kl(l,v)
@@ -408,8 +429,8 @@ if __name__ == '__main__':
 
     file_ = get_full_path(arguments['<file.out>'])
 
-    if arguments["--ezfio"]:
-        ezfio_file = get_full_path(arguments["--ezfio"])
+    if arguments["-o"]:
+        ezfio_file = get_full_path(arguments["<ezfio_directory>"])
     else:
         ezfio_file = "{0}.ezfio".format(file_)
 
@@ -421,3 +442,4 @@ if __name__ == '__main__':
         print file_, 'recognized as', str(res_file).split('.')[-1].split()[0]
 
     write_ezfio(res_file, ezfio_file)
+    os.system("qp_run save_ortho_mos "+ezfio_file)
diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template
index 9c7a1386..af9b295c 100644
--- a/scripts/ezfio_interface/qp_edit_template
+++ b/scripts/ezfio_interface/qp_edit_template
@@ -1,6 +1,10 @@
-open Qputils;;
-open Qptypes;;
-open Core.Std;;
+(*
+  vim::syntax=ocaml
+ *)
+
+open Qputils
+open Qptypes
+open Core.Std
 
 (** Interactive editing of the input.
 
@@ -18,7 +22,7 @@ type keyword =
 | Mo_basis
 | Nuclei
 {keywords}
-;;
+
 
 
 let keyword_to_string = function
@@ -28,7 +32,7 @@ let keyword_to_string = function
 | Mo_basis              -> "MO basis"
 | Nuclei                -> "Molecule"
 {keywords_to_string}
-;;
+
 
 
 
@@ -42,7 +46,7 @@ let file_header filename =
 Editing file `%s`
 
 " filename
-;;
+
   
 
 (** Creates the header of a section *)
@@ -50,7 +54,7 @@ let make_header kw =
   let s = keyword_to_string kw in
   let l = String.length s in
   "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
-;;
+
 
 
 (** Returns the rst string of section [s] *)
@@ -82,7 +86,7 @@ let get s =
     | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
   in 
   rst
-;;
+
 
 
 (** Applies the changes from the string [str] corresponding to section [s] *)
@@ -121,7 +125,7 @@ let set str s =
         | Ao_basis         -> () (* TODO *)
         | Mo_basis         -> () (* TODO *)
     end
-;;
+
 
 
 (** Creates the temporary file for interactive editing *)
@@ -135,11 +139,19 @@ let create_temp_file ezfio_filename fields =
       )
   end
   ; temp_filename
-;;
 
 
+                                                 
 
-let run check_only ezfio_filename =
+
+let run check_only ?ndet ?state ezfio_filename =
+
+  (* Set check_only if the arguments are not empty *)
+  let check_only =
+    match ndet, state with
+    | None, None -> check_only
+    | _ -> true
+  in
 
   (* Open EZFIO *)
   if (not (Sys.file_exists_exn ezfio_filename)) then
@@ -147,6 +159,19 @@ let run check_only ezfio_filename =
 
   Ezfio.set_file ezfio_filename;
 
+  begin
+    match ndet with
+    | None -> ()
+    | Some n -> Input.Determinants_by_hand.update_ndet (Det_number.of_int n)
+  end;
+
+  begin
+    match state with
+    | None -> ()
+    | Some n -> Input.Determinants_by_hand.extract_state (States_number.of_int n)
+  end;
+
+
   (*
   let output = (file_header ezfio_filename) :: (
     List.map ~f:get [
@@ -196,7 +221,7 @@ let run check_only ezfio_filename =
 
   (* Remove temp_file *)
   Sys.remove temp_filename
-;;
+
 
 
 (** Create a backup file in case of an exception *)
@@ -207,7 +232,7 @@ let create_backup ezfio_filename =
   "
     ezfio_filename ezfio_filename ezfio_filename
   |> Sys.command_exn
-;;
+
 
 
 (** Restore the backup file when an exception occuprs *)
@@ -215,7 +240,7 @@ let restore_backup ezfio_filename =
   Printf.sprintf "tar -zxf %s/backup.tgz"
     ezfio_filename 
   |> Sys.command_exn
-;;
+
 
 
 let spec =
@@ -223,12 +248,12 @@ let spec =
   empty 
   +> flag "-c" no_arg
      ~doc:"Checks the input data"
-(*
-  +> flag "o" (optional string)
-     ~doc:"Prints output data"
-*)
+  +> flag "ndet" (optional int)
+     ~doc:"int Truncate the wavefunction to the target number of determinants"
+  +> flag "state" (optional int)
+     ~doc:"int Pick the state as a new wavefunction."
   +> anon ("ezfio_file" %: string)
-;;
+
 
 let command = 
     Command.basic 
@@ -245,9 +270,9 @@ Edit input data
        with
        | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
     *)
-    (fun c ezfio_file () ->
+    (fun c ndet state ezfio_file () ->
        try
-          run c ezfio_file ;
+          run c ?ndet ?state ezfio_file ;
           (* create_backup ezfio_file; *)
        with
        | Failure exc 
@@ -268,12 +293,12 @@ Edit input data
            raise e
          end
     )
-;;
+
 
 let () =
   Command.run command;
   exit 0
-;;
+
 
 
 
diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f
index 10ab6f67..e50cf25a 100644
--- a/src/Bitmask/bitmasks.irp.f
+++ b/src/Bitmask/bitmasks.irp.f
@@ -2,11 +2,16 @@ use bitmasks
 
 BEGIN_PROVIDER [ integer, N_int ]
   implicit none
+  include 'Utils/constants.include.F'
   BEGIN_DOC
   ! Number of 64-bit integers needed to represent determinants as binary strings
   END_DOC
   N_int = (mo_tot_num-1)/bit_kind_size + 1
- call write_int(6,N_int, 'N_int')
+  call write_int(6,N_int, 'N_int')
+  if (N_int > N_int_max) then
+    stop 'N_int > N_int_max'
+  endif
+
 END_PROVIDER
 
 
diff --git a/src/Davidson/diagonalization.irp.f b/src/Davidson/diagonalization.irp.f
index 9bbd00f5..fe82a8fb 100644
--- a/src/Davidson/diagonalization.irp.f
+++ b/src/Davidson/diagonalization.irp.f
@@ -355,7 +355,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = ' Iter'
   do i=1,N_st
-    write_buffer = trim(write_buffer)//'           Energy         Residual'
+    write_buffer = trim(write_buffer)//'            Energy          Residual'
   enddo
   write(iunit,'(A)') trim(write_buffer)
   write_buffer = '===== '
@@ -502,7 +502,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
         endif
       enddo
       
-      write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))')  iter, to_print(:,1:N_st)
+      write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,E16.6))')  iter, to_print(:,1:N_st)
       call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
       if (converged) then
         exit
diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f
index f458e32d..bf56855a 100644
--- a/src/Davidson/diagonalization_hs2.irp.f
+++ b/src/Davidson/diagonalization_hs2.irp.f
@@ -413,7 +413,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
         endif
       enddo
       
-      write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))')  iter, to_print(1:3,1:N_st)
+      write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))')  iter, to_print(1:3,1:N_st)
       call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
       do k=1,N_st
         if (residual_norm(k) > 1.e8) then
@@ -838,7 +838,7 @@ subroutine davidson_diag_hjj_sjj_mmap(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz
         endif
       enddo
       
-      write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))')  iter, to_print(1:3,1:N_st)
+      write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))')  iter, to_print(1:3,1:N_st)
       call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
       do k=1,N_st
         if (residual_norm(k) > 1.e8) then
diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f
index 8ef574d4..b096d407 100644
--- a/src/Davidson/u0Hu0.irp.f
+++ b/src/Davidson/u0Hu0.irp.f
@@ -90,7 +90,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
   enddo
   !$OMP END DO
 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,2)
     do i=shortcut(sh,2),shortcut(sh+1,2)-1
       org_i = sort_idx(i,2)
@@ -123,7 +123,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
   enddo
   !$OMP END DO
 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,1)
     do sh2=1,shortcut(0,1)
       if (sh==sh2) cycle
@@ -342,7 +342,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_
 !    istep = 1+ int(workload*target_workload_inv)
     istep = 1
     do blockb2=0, istep-1
-      write(tmp_task,'(3(I9,X),''|'',X)') sh, blockb2, istep
+      write(tmp_task,'(3(I9,1X),''|'',1X)') sh, blockb2, istep
       task = task//tmp_task
       ipos += 32
       if (ipos+32 > iposmax) then
@@ -444,7 +444,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   enddo
   !$OMP END DO
 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,4)
   do sh=1,shortcut(0,2)
     do i=shortcut(sh,2),shortcut(sh+1,2)-1
       org_i = sort_idx(i,2)
@@ -477,7 +477,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   enddo
   !$OMP END DO
 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,4)
   do sh=1,shortcut(0,1)
     do sh2=1,shortcut(0,1)
       if (sh==sh2) cycle
@@ -609,3 +609,352 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   deallocate (shortcut, sort_idx, sorted, version, ut)
 end
 
+
+
+
+
+subroutine H_S2_u_0_nstates_new(v_0,s_0,N_st,sze_8)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
+  !
+  ! n : number of determinants
+  !
+  ! H_jj : array of <j|H|j>
+  !
+  ! S2_jj : array of <j|S^2|j>
+  END_DOC
+  integer, intent(in)            :: N_st,sze_8
+  double precision, intent(out)  :: v_0(sze_8,N_st), s_0(sze_8,N_st)
+
+  
+  PROVIDE ref_bitmask_energy
+
+  double precision               :: hij, s2 
+  integer                        :: i,j
+  integer                        :: k_a, k_b, l_a, l_b, m_a, m_b
+  integer                        :: degree, istate
+  integer                        :: krow, kcol, krow_b, kcol_b
+  integer                        :: lrow, lcol
+  integer                        :: mrow, mcol
+  integer(bit_kind)              :: spindet(N_int)
+  integer(bit_kind)              :: tmp_det(N_int,2)
+  integer(bit_kind)              :: tmp_det2(N_int,2)
+  integer(bit_kind)              :: tmp_det3(N_int,2)
+  integer(bit_kind), allocatable :: buffer(:,:)
+  double precision               :: ck(N_st), cl(N_st), cm(N_st)
+  integer                        :: n_singles, n_doubles
+  integer, allocatable           :: singles(:), doubles(:)
+  integer, allocatable           :: idx(:), idx0(:)
+  logical, allocatable           :: is_single_a(:)
+
+  allocate( buffer(N_int,N_det_alpha_unique),                        &
+      singles(N_det_alpha_unique), doubles(N_det_alpha_unique),      &
+      is_single_a(N_det_alpha_unique),                               &
+      idx(N_det_alpha_unique), idx0(N_det_alpha_unique) )
+  
+  v_0 = 0.d0
+
+  do k_a=1,N_det-1
+    
+    ! Initial determinant is at k_a in alpha-major representation
+    ! -----------------------------------------------------------------------
+    
+    krow = psi_bilinear_matrix_rows(k_a)
+    kcol = psi_bilinear_matrix_columns(k_a)
+    
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
+    
+    ! Initial determinant is at k_b in beta-major representation
+    ! ----------------------------------------------------------------------
+    
+    k_b = psi_bilinear_matrix_order_reverse(k_a)
+
+    ! Diagonal contribution
+    ! ---------------------
+
+    double precision, external :: diag_H_mat_elem
+  
+    v_0(k_a,1:N_st) = v_0(k_a,1:N_st) + diag_H_mat_elem(tmp_det,N_int) * &
+        psi_bilinear_matrix_values(k_a,1:N_st)
+
+    
+    ! Get all single and double alpha excitations
+    ! ===========================================
+    
+    spindet(1:N_int) = tmp_det(1:N_int,1)
+    
+    ! Loop inside the beta column to gather all the connected alphas
+    i=1
+    l_a = k_a+1
+    lcol = psi_bilinear_matrix_columns(l_a)
+    do while ( (lcol == kcol).and.(l_a <= N_det) )
+      lrow = psi_bilinear_matrix_rows(l_a)
+      buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow)
+      idx(i) = lrow
+      i=i+1
+      l_a = l_a + 1
+      lcol = psi_bilinear_matrix_columns(l_a)
+    enddo
+    i = i-1
+    
+    call get_all_spin_singles_and_doubles(                           &
+        buffer, idx, spindet, N_int, i,                              &
+        singles, doubles, n_singles, n_doubles )
+
+    ! Compute Hij for all alpha singles
+    ! ----------------------------------
+
+    l_a = k_a
+    lrow = psi_bilinear_matrix_rows(l_a)
+    tmp_det2(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
+    do i=1,n_singles
+      do while ( lrow < singles(i) )
+        l_a = l_a+1
+        lrow = psi_bilinear_matrix_rows(l_a)
+      enddo
+      tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
+      call i_H_j_mono_spin( tmp_det, tmp_det2, N_int, 1, hij)
+      v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
+      v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
+    enddo
+    
+    ! Compute Hij for all alpha doubles
+    ! ----------------------------------
+    
+    l_a = k_a
+    lrow = psi_bilinear_matrix_rows(l_a)
+    do i=1,n_doubles
+      do while (lrow < doubles(i))
+        l_a = l_a+1
+        lrow = psi_bilinear_matrix_rows(l_a)
+      enddo
+      call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, doubles(i)), N_int, hij)
+      v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
+      v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
+    enddo
+    
+    
+    
+    ! Get all single and double beta excitations
+    ! ===========================================
+    
+    spindet(1:N_int) = tmp_det(1:N_int,2)
+    
+    ! Loop inside the alpha row to gather all the connected betas
+    i=1
+    l_b = k_b+1
+    lrow = psi_bilinear_matrix_transp_rows(l_b)
+    do while ( (lrow == krow).and.(l_b <= N_det) )
+      lcol = psi_bilinear_matrix_transp_columns(l_b)
+      buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol)
+      idx(i) = lcol
+      i=i+1
+      l_b = l_b + 1
+      lrow = psi_bilinear_matrix_transp_rows(l_b)
+    enddo
+    i = i-1
+    
+    call get_all_spin_singles_and_doubles(                           &
+        buffer, idx, spindet, N_int, i,                              &
+        singles, doubles, n_singles, n_doubles )
+    
+    ! Compute Hij for all beta singles
+    ! ----------------------------------
+    
+    l_b = k_b
+    lcol = psi_bilinear_matrix_transp_columns(l_b)
+    tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
+    do i=1,n_singles
+      do while ( lcol < singles(i) )
+        l_b = l_b+1
+        lcol = psi_bilinear_matrix_transp_columns(l_b)
+      enddo
+      tmp_det2(1:N_int,2) = psi_det_beta_unique (1:N_int, lcol)
+      l_a = psi_bilinear_matrix_transp_order(l_b)
+      call i_H_j_mono_spin( tmp_det, tmp_det2, N_int, 2, hij)
+      v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
+      v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
+    enddo
+    
+    ! Compute Hij for all beta doubles
+    ! ----------------------------------
+    
+    l_b = k_b
+    lcol = psi_bilinear_matrix_transp_columns(l_b)
+    do i=1,n_doubles
+      do while (lcol < doubles(i))
+        l_b = l_b+1
+        lcol = psi_bilinear_matrix_transp_columns(l_b)
+      enddo
+      l_a = psi_bilinear_matrix_transp_order(l_b)
+      call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, doubles(i)), N_int, hij)
+      v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
+      v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st)
+    enddo
+    
+  end do
+
+
+  ! Alpha/Beta double excitations
+  ! =============================
+    
+  do i=1,N_det_beta_unique
+    idx0(i) = i
+  enddo
+  is_single_a(:) = .False.
+
+  k_a=1
+  do i=1,N_det_beta_unique
+
+    ! Select a beta determinant
+    ! -------------------------
+
+    spindet(1:N_int) = psi_det_beta_unique(1:N_int, i)
+    tmp_det(1:N_int,2) = spindet(1:N_int)
+
+    call get_all_spin_singles(                                       &
+        psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, &
+        singles, n_singles )
+
+    do j=1,n_singles
+      is_single_a( singles(j) ) = .True.
+    enddo
+
+    ! For all alpha.beta pairs with the selected beta
+    ! -----------------------------------------------
+
+    kcol = psi_bilinear_matrix_columns(k_a)
+    do while (kcol < i)
+      k_a = k_a+1
+      if (k_a > N_det) exit
+      kcol = psi_bilinear_matrix_columns(k_a)
+    enddo
+
+    do while (kcol == i)
+
+      krow = psi_bilinear_matrix_rows(k_a)
+      tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow)
+
+      ! Loop over all alpha.beta pairs with a single exc alpha
+      ! ------------------------------------------------------
+
+      l_a = k_a+1
+      if (l_a > N_det) exit
+      lrow = psi_bilinear_matrix_rows(l_a)
+      lcol = psi_bilinear_matrix_columns(l_a)
+
+      do while (lrow == krow)
+
+        ! Loop over all alpha.beta pairs with a single exc alpha
+        ! ------------------------------------------------------
+        if (is_single_a(lrow)) then
+
+           tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow)
+          
+           ! Build list of singly excited beta
+           ! ---------------------------------
+
+           m_b = psi_bilinear_matrix_order_reverse(l_a)
+           m_b = m_b+1
+           j=1
+           do while ( (mrow == lrow) )
+             mcol = psi_bilinear_matrix_transp_columns(m_b)
+             buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol)
+             idx(j) = mcol
+             j = j+1
+             m_b = m_b+1
+             if (m_b <= N_det) exit
+             mrow = psi_bilinear_matrix_transp_rows(m_b)
+           enddo
+           j=j-1
+
+           call get_all_spin_singles(                                &
+               buffer, idx, tmp_det(1,2), N_int, j,                  &
+               doubles, n_doubles)
+
+           ! Compute Hij for all doubles 
+           ! ---------------------------
+
+           m_b = psi_bilinear_matrix_order(l_a)+1
+           mcol = psi_bilinear_matrix_transp_columns(m_b)
+           do j=1,n_doubles
+             tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) )
+             call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij)
+             do while (mcol /= doubles(j))
+               m_b = m_b+1
+               if (m_b > N_det) exit
+               mcol = psi_bilinear_matrix_transp_columns(m_b)
+             enddo
+             m_a = psi_bilinear_matrix_order_reverse(m_b)
+!             v_0(m_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st)
+!             v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(m_a,1:N_st)
+           enddo
+
+        endif
+        l_a = l_a+1
+        if (l_a > N_det) exit
+        lrow = psi_bilinear_matrix_rows(l_a)
+        lcol = psi_bilinear_matrix_columns(l_a)
+      enddo
+
+      k_b = k_b+1
+      if (k_b > N_det) exit
+      kcol = psi_bilinear_matrix_transp_columns(k_b)
+    enddo
+    
+    do j=1,n_singles
+      is_single_a( singles(j) ) = .False.
+    enddo
+
+  enddo  
+
+    
+end
+
+
+subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: N_st,n,Nint, sze_8
+  integer(bit_kind), intent(in)  :: keys_tmp(Nint,2,n)
+  double precision, intent(out)  :: v_0(sze_8,N_st), s_0(sze_8,N_st)
+  double precision, intent(in)   :: u_0(sze_8,N_st)
+  double precision, intent(in)   :: H_jj(n), S2_jj(n)
+  
+  PROVIDE ref_bitmask_energy
+
+  double precision, allocatable :: vt(:,:)
+  integer, allocatable :: idx(:)
+  integer                        :: i,j, jj
+  double precision               :: hij
+
+  do i=1,n
+    v_0(i,:) = H_jj(i) * u_0(i,:)
+  enddo
+
+  allocate(idx(0:n), vt(N_st,n))
+  Vt = 0.d0
+  do i=2,n
+    idx(0) = i
+    call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
+    do jj=1,idx(0)
+      j = idx(jj)
+      double precision :: phase
+      integer :: degree
+      integer :: exc(0:2,2,2)
+      call get_excitation(keys_tmp(1,1,j),keys_tmp(1,1,i),exc,degree,phase,Nint)
+      if ((degree == 2).and.(exc(0,1,1)==1)) cycle
+!      if (exc(0,1,2) /= 0) cycle
+      call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
+      vt (:,i) = vt (:,i) + hij*u_0(j,:)
+      vt (:,j) = vt (:,j) + hij*u_0(i,:)
+    enddo
+  enddo
+  do i=1,n
+    v_0(i,:) = v_0(i,:) + vt(:,i)
+  enddo
+end
+
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index 41433ed1..7274760e 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -78,96 +78,68 @@ END_PROVIDER
    double precision               :: ck, cl, ckl
    double precision               :: phase
    integer                        :: h1,h2,p1,p2,s1,s2, degree
-   integer(bit_kind)              :: tmp_det(N_int,2), tmp_det2(N_int,2)
-   integer                        :: exc(0:2,2,2),n_occ(2)
+   integer(bit_kind)              :: tmp_det(N_int,2), tmp_det2(N_int)
+   integer                        :: exc(0:2,2),n_occ(2)
    double precision, allocatable  :: tmp_a(:,:,:), tmp_b(:,:,:)
    integer                        :: krow, kcol, lrow, lcol
 
-     one_body_dm_mo_alpha = 0.d0
-     one_body_dm_mo_beta  = 0.d0
-     !$OMP PARALLEL DEFAULT(NONE)                                         &
-        !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
-        !$OMP  tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
-        !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
-        !$OMP  elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
-        !$OMP  mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
-        !$OMP  psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
-        !$OMP  psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
-     allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
-     tmp_a = 0.d0
-     tmp_b = 0.d0
-     !$OMP DO SCHEDULE(guided)
-     do k=1,N_det
-       krow = psi_bilinear_matrix_rows(k) 
-       kcol = psi_bilinear_matrix_columns(k) 
-       tmp_det(:,1) = psi_det(:,1, krow)
-       tmp_det(:,2) = psi_det(:,2, kcol)
-       call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
-       do m=1,N_states
-         ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
-         do l=1,elec_alpha_num
-           j = occ(l,1)
-           tmp_a(j,j,m) += ck
-         enddo
-         do l=1,elec_beta_num
-           j = occ(l,2)
-           tmp_b(j,j,m) += ck
-         enddo
-       enddo
+   PROVIDE psi_det
 
-       l = k+1
-       lrow = psi_bilinear_matrix_rows(l) 
-       lcol = psi_bilinear_matrix_columns(l) 
-       do while ( lcol == kcol )
-         tmp_det2(:,1) = psi_det(:,1, lrow)
-         tmp_det2(:,2) = psi_det(:,2, lcol)
-         call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
-         if (degree == 1) then
-           call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
-           call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-           do m=1,N_states
-             ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
-             if (s1==1) then
-               tmp_a(h1,p1,m) += ckl
-               tmp_a(p1,h1,m) += ckl
-             else
-               tmp_b(h1,p1,m) += ckl
-               tmp_b(p1,h1,m) += ckl
-             endif
-           enddo
-         endif
-         l = l+1
-         if (l>N_det) exit
-         lrow = psi_bilinear_matrix_rows(l) 
-         lcol = psi_bilinear_matrix_columns(l) 
-       enddo
+  one_body_dm_mo_alpha = 0.d0
+  one_body_dm_mo_beta  = 0.d0
+  !$OMP PARALLEL DEFAULT(NONE)                                         &
+    !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
+    !$OMP  tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
+    !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
+    !$OMP  elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
+    !$OMP  mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
+    !$OMP  psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
+    !$OMP  psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique, &
+    !$OMP  psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
+  allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
+  tmp_a = 0.d0
+  tmp_b = 0.d0
+  !$OMP DO SCHEDULE(guided)
+  do k=1,N_det
+    krow = psi_bilinear_matrix_rows(k) 
+    kcol = psi_bilinear_matrix_columns(k) 
+    tmp_det(:,1) = psi_det_alpha_unique(:,krow)
+    tmp_det(:,2) = psi_det_beta_unique (:,kcol)
+    call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
+    do m=1,N_states
+      ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
+      do l=1,elec_alpha_num
+        j = occ(l,1)
+        tmp_a(j,j,m) += ck
+      enddo
+      do l=1,elec_beta_num
+        j = occ(l,2)
+        tmp_b(j,j,m) += ck
+      enddo
+    enddo
 
-       l = k+1
-       lrow = psi_bilinear_matrix_transp_rows(l) 
-       lcol = psi_bilinear_matrix_transp_columns(l) 
-       do while ( lrow == krow )
-         tmp_det2(:,1) = psi_det(:,1, lrow)
-         tmp_det2(:,2) = psi_det(:,2, lcol)
-         call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
-         if (degree == 1) then
-           call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
-           call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-           do m=1,N_states
-             ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase
-             if (s1==1) then
-               tmp_a(h1,p1,m) += ckl
-               tmp_a(p1,h1,m) += ckl
-             else
-               tmp_b(h1,p1,m) += ckl
-               tmp_b(p1,h1,m) += ckl
-             endif
-           enddo
-         endif
-         l = l+1
-         if (l>N_det) exit
-         lrow = psi_bilinear_matrix_transp_rows(l) 
-         lcol = psi_bilinear_matrix_transp_columns(l) 
-       enddo
+    l = k+1
+    lrow = psi_bilinear_matrix_rows(l) 
+    lcol = psi_bilinear_matrix_columns(l) 
+    ! Fix beta determinant, loop over alphas
+    do while ( lcol == kcol )
+      tmp_det2(:) = psi_det_alpha_unique(:, lrow)
+      call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int)
+      if (degree == 1) then
+        exc = 0
+        call get_mono_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int)
+        call decode_exc_spin(exc,h1,p1,h2,p2)
+        do m=1,N_states
+          ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
+          tmp_a(h1,p1,m) += ckl
+          tmp_a(p1,h1,m) += ckl
+        enddo
+      endif
+      l = l+1
+      if (l>N_det) exit
+      lrow = psi_bilinear_matrix_rows(l) 
+      lcol = psi_bilinear_matrix_columns(l) 
+    enddo
 
      enddo
      !$OMP END DO NOWAIT
@@ -364,3 +336,74 @@ END_PROVIDER
 
 END_PROVIDER
 
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_old, (mo_tot_num_align,mo_tot_num,N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_old, (mo_tot_num_align,mo_tot_num,N_states) ]
+   implicit none
+   BEGIN_DOC
+   ! Alpha and beta one-body density matrix for each state
+   END_DOC
+
+   integer                        :: j,k,l,m
+   integer                        :: occ(N_int*bit_kind_size,2)
+   double precision               :: ck, cl, ckl
+   double precision               :: phase
+   integer                        :: h1,h2,p1,p2,s1,s2, degree
+   integer                        :: exc(0:2,2,2),n_occ(2)
+   double precision, allocatable  :: tmp_a(:,:,:), tmp_b(:,:,:)
+
+     one_body_dm_mo_alpha_old = 0.d0
+     one_body_dm_mo_beta_old  = 0.d0
+     !$OMP PARALLEL DEFAULT(NONE)                                         &
+        !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
+        !$OMP  tmp_a, tmp_b, n_occ)&
+        !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
+        !$OMP  elec_beta_num,one_body_dm_mo_alpha_old,one_body_dm_mo_beta_old,N_det,mo_tot_num_align,&
+        !$OMP  mo_tot_num)
+     allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
+     tmp_a = 0.d0
+     tmp_b = 0.d0
+     !$OMP DO SCHEDULE(dynamic)
+     do k=1,N_det
+       call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int)
+       do m=1,N_states
+         ck = psi_coef(k,m)*psi_coef(k,m) 
+         do l=1,elec_alpha_num
+           j = occ(l,1)
+           tmp_a(j,j,m) += ck
+         enddo
+         do l=1,elec_beta_num
+           j = occ(l,2)
+           tmp_b(j,j,m) += ck
+         enddo
+       enddo
+       do l=1,k-1
+         call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
+         if (degree /= 1) then
+           cycle
+         endif
+         call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
+         call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+         do m=1,N_states
+           ckl = psi_coef(k,m) * psi_coef(l,m) * phase 
+           if (s1==1) then
+             tmp_a(h1,p1,m) += ckl
+             tmp_a(p1,h1,m) += ckl
+           else
+             tmp_b(h1,p1,m) += ckl
+             tmp_b(p1,h1,m) += ckl
+           endif
+         enddo
+       enddo
+     enddo
+     !$OMP END DO NOWAIT
+     !$OMP CRITICAL
+     one_body_dm_mo_alpha_old(:,:,:) = one_body_dm_mo_alpha_old(:,:,:) + tmp_a(:,:,:)
+     !$OMP END CRITICAL
+     !$OMP CRITICAL
+     one_body_dm_mo_beta_old(:,:,:)  = one_body_dm_mo_beta_old(:,:,:)  + tmp_b(:,:,:)
+     !$OMP END CRITICAL
+     deallocate(tmp_a,tmp_b)
+     !$OMP END PARALLEL
+
+END_PROVIDER
diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f
index bed3327d..2644801e 100644
--- a/src/Determinants/determinants.irp.f
+++ b/src/Determinants/determinants.irp.f
@@ -23,7 +23,7 @@ BEGIN_PROVIDER [ integer, N_det ]
  ! Number of determinants in the wave function
  END_DOC
  logical                        :: exists
- character*64                   :: label
+ character*(64)                   :: label
  PROVIDE ezfio_filename
  PROVIDE nproc
  if (read_wf) then
@@ -88,7 +88,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
  END_DOC
   integer                        :: i
   logical                        :: exists
-  character*64                   :: label
+  character*(64)                   :: label
   
   psi_det = 0_bit_kind
   if (read_wf) then
diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f
index a807513c..a6e69fb5 100644
--- a/src/Determinants/s2.irp.f
+++ b/src/Determinants/s2.irp.f
@@ -252,8 +252,8 @@ end
 subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates)
   implicit none
   use bitmasks
-  integer(bit_kind), intent(in)  :: psi_keys_tmp(N_int,2,nmax_keys)
   integer, intent(in)            :: n,nmax_coefs,nmax_keys,nstates
+  integer(bit_kind), intent(in)  :: psi_keys_tmp(N_int,2,nmax_keys)
   double precision, intent(in)   :: psi_coefs_tmp(nmax_coefs,nstates)
   double precision, intent(out)  :: s2(nstates,nstates)
   double precision               :: s2_tmp,accu
@@ -344,7 +344,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nsta
   
   print*,'S^2 matrix in the basis of the states considered'
   do i = 1, nstates
-    write(*,'(100(F5.2,X))')s2(i,:)
+    write(*,'(100(F5.2,1X))')s2(i,:)
   enddo
   
   double precision               :: accu_precision_diag,accu_precision_of_diag
@@ -370,7 +370,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nsta
   
   print*,'Modified S^2 matrix that will be diagonalized'
   do i = 1, nstates
-    write(*,'(10(F5.2,X))')s2(i,:)
+    write(*,'(10(F5.2,1X))')s2(i,:)
     s2(i,i) = s2(i,i)
   enddo
   
diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f
index 2a76a079..f4af1b60 100644
--- a/src/Determinants/slater_rules.irp.f
+++ b/src/Determinants/slater_rules.irp.f
@@ -1,32 +1,59 @@
 subroutine get_excitation_degree(key1,key2,degree,Nint)
   use bitmasks
+  include 'Utils/constants.include.F'
   implicit none
   BEGIN_DOC
   ! Returns the excitation degree between two determinants
   END_DOC
   integer, intent(in)            :: Nint
-  integer(bit_kind), intent(in)  :: key1(Nint,2)
-  integer(bit_kind), intent(in)  :: key2(Nint,2)
+  integer(bit_kind), intent(in)  :: key1(Nint*2)
+  integer(bit_kind), intent(in)  :: key2(Nint*2)
   integer, intent(out)           :: degree
   
+  integer(bit_kind)              :: xorvec(2*N_int_max)
   integer                        :: l
   
   ASSERT (Nint > 0)
   
-  degree = popcnt(xor( key1(1,1), key2(1,1))) +                      &
-      popcnt(xor( key1(1,2), key2(1,2)))
-  !DIR$ NOUNROLL
-  do l=2,Nint
-    degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) +            &
-        popcnt(xor( key1(l,2), key2(l,2)))
-  enddo
-  ASSERT (degree >= 0)
+  select case (Nint)
+
+    case (1)
+      xorvec(1) = xor( key1(1), key2(1))
+      xorvec(2) = xor( key1(2), key2(2))
+      degree = sum(popcnt(xorvec(1:2)))
+
+    case (2)
+      xorvec(1) = xor( key1(1), key2(1))
+      xorvec(2) = xor( key1(2), key2(2))
+      xorvec(3) = xor( key1(3), key2(3))
+      xorvec(4) = xor( key1(4), key2(4))
+      degree = sum(popcnt(xorvec(1:4)))
+
+    case (3)
+      do l=1,6
+        xorvec(l) = xor( key1(l), key2(l))
+      enddo
+      degree = sum(popcnt(xorvec(1:6)))
+
+    case (4)
+      do l=1,8
+        xorvec(l) = xor( key1(l), key2(l))
+      enddo
+      degree = sum(popcnt(xorvec(1:8)))
+
+    case default
+      do l=1,ishft(Nint,1)
+        xorvec(l) = xor( key1(l), key2(l))
+      enddo
+      degree = sum(popcnt(xorvec(1:l)))
+  
+  end select
+
   degree = ishft(degree,-1)
   
 end
 
 
-
 subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
   use bitmasks
   implicit none
@@ -139,72 +166,6 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
 end
 
 
-subroutine decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Decodes the exc arrays returned by get_excitation.
-  ! h1,h2 : Holes
-  ! p1,p2 : Particles
-  ! s1,s2 : Spins (1:alpha, 2:beta)
-  ! degree : Degree of excitation
-  END_DOC
-  integer, intent(in)            :: exc(0:2,2,2),degree
-  integer*2, intent(out)           :: h1,h2,p1,p2,s1,s2
-  ASSERT (degree > 0)
-  ASSERT (degree < 3)
-  
-  select case(degree)
-    case(2)
-      if (exc(0,1,1) == 2) then
-        h1 = exc(1,1,1)
-        h2 = exc(2,1,1)
-        p1 = exc(1,2,1)
-        p2 = exc(2,2,1)
-        s1 = 1
-        s2 = 1
-      else if (exc(0,1,2) == 2) then
-        h1 = exc(1,1,2)
-        h2 = exc(2,1,2)
-        p1 = exc(1,2,2)
-        p2 = exc(2,2,2)
-        s1 = 2
-        s2 = 2
-      else
-        h1 = exc(1,1,1)
-        h2 = exc(1,1,2)
-        p1 = exc(1,2,1)
-        p2 = exc(1,2,2)
-        s1 = 1
-        s2 = 2
-      endif
-    case(1)
-      if (exc(0,1,1) == 1) then
-        h1 = exc(1,1,1)
-        h2 = 0
-        p1 = exc(1,2,1)
-        p2 = 0
-        s1 = 1
-        s2 = 0
-      else
-        h1 = exc(1,1,2)
-        h2 = 0
-        p1 = exc(1,2,2)
-        p2 = 0
-        s1 = 2
-        s2 = 0
-      endif
-    case(0)
-      h1 = 0
-      p1 = 0
-      h2 = 0
-      p2 = 0
-      s1 = 0
-      s2 = 0
-  end select
-end
-
-
 subroutine get_double_excitation(det1,det2,exc,phase,Nint)
   use bitmasks
   implicit none
@@ -2154,8 +2115,8 @@ end
 subroutine get_phase(key1,key2,phase,Nint)
   use bitmasks
   implicit none
-  integer(bit_kind), intent(in)  :: key1(Nint,2), key2(Nint,2)
   integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key1(Nint,2), key2(Nint,2)
   double precision, intent(out)  :: phase
   BEGIN_DOC
 ! Returns the phase between key1 and key2
@@ -2224,3 +2185,423 @@ subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze)
   call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
   e_0 =  u_dot_v(v_0,u_0,sze)
 end
+
+
+
+! Spin-determinant routines
+! -------------------------
+
+subroutine get_excitation_degree_spin(key1,key2,degree,Nint)
+  use bitmasks
+  include 'Utils/constants.include.F'
+  implicit none
+  BEGIN_DOC
+  ! Returns the excitation degree between two determinants
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key1(Nint)
+  integer(bit_kind), intent(in)  :: key2(Nint)
+  integer, intent(out)           :: degree
+  
+  integer(bit_kind)              :: xorvec(N_int_max)
+  integer                        :: l
+  
+  ASSERT (Nint > 0)
+  
+  select case (Nint)
+
+    case (1)
+      xorvec(1) = xor( key1(1), key2(1))
+      degree = popcnt(xorvec(1))
+
+    case (2)
+      xorvec(1) = xor( key1(1), key2(1))
+      xorvec(2) = xor( key1(2), key2(2))
+      degree = popcnt(xorvec(1))+popcnt(xorvec(2))
+
+    case (3)
+      xorvec(1) = xor( key1(1), key2(1))
+      xorvec(2) = xor( key1(2), key2(2))
+      xorvec(3) = xor( key1(3), key2(3))
+      degree = sum(popcnt(xorvec(1:3)))
+
+    case (4)
+      xorvec(1) = xor( key1(1), key2(1))
+      xorvec(2) = xor( key1(2), key2(2))
+      xorvec(3) = xor( key1(3), key2(3))
+      xorvec(4) = xor( key1(4), key2(4))
+      degree = sum(popcnt(xorvec(1:4)))
+
+    case default
+      do l=1,N_int
+        xorvec(l) = xor( key1(l), key2(l))
+      enddo
+      degree = sum(popcnt(xorvec(1:Nint)))
+  
+  end select
+
+  degree = ishft(degree,-1)
+  
+end
+
+
+subroutine get_excitation_spin(det1,det2,exc,degree,phase,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns the excitation operators between two determinants and the phase
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: det1(Nint)
+  integer(bit_kind), intent(in)  :: det2(Nint)
+  integer, intent(out)           :: exc(0:2,2)
+  integer, intent(out)           :: degree
+  double precision, intent(out)  :: phase
+  ! exc(number,hole/particle)
+  ! ex :
+  ! exc(0,1) = number of holes
+  ! exc(0,2) = number of particles
+  ! exc(1,2) = first particle 
+  ! exc(1,1) = first hole     
+  
+  ASSERT (Nint > 0)
+  
+  !DIR$ FORCEINLINE
+  call get_excitation_degree_spin(det1,det2,degree,Nint)
+  select case (degree)
+      
+    case (3:)
+      degree = -1
+      return
+      
+    case (2)
+      call get_double_excitation_spin(det1,det2,exc,phase,Nint)
+      return
+      
+    case (1)
+      call get_mono_excitation_spin(det1,det2,exc,phase,Nint)
+      return
+      
+    case(0)
+      return
+      
+  end select
+end
+
+subroutine decode_exc_spin(exc,h1,p1,h2,p2)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Decodes the exc arrays returned by get_excitation.
+  ! h1,h2 : Holes
+  ! p1,p2 : Particles
+  END_DOC
+  integer, intent(in)            :: exc(0:2,2)
+  integer, intent(out)           :: h1,h2,p1,p2
+  
+  select case (exc(0,1))
+    case(2)
+      h1 = exc(1,1)
+      h2 = exc(2,1)
+      p1 = exc(1,2)
+      p2 = exc(2,2)
+    case(1)
+      h1 = exc(1,1)
+      h2 = 0
+      p1 = exc(1,2)
+      p2 = 0
+    case default
+      h1 = 0
+      p1 = 0
+      h2 = 0
+      p2 = 0
+  end select
+end
+
+
+subroutine get_double_excitation_spin(det1,det2,exc,phase,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns the two excitation operators between two doubly excited spin-determinants
+  ! and the phase
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: det1(Nint)
+  integer(bit_kind), intent(in)  :: det2(Nint)
+  integer, intent(out)           :: exc(0:2,2)
+  double precision, intent(out)  :: phase
+  integer                        :: tz
+  integer                        :: l, idx_hole, idx_particle, ishift
+  integer                        :: nperm
+  integer                        :: i,j,k,m,n
+  integer                        :: high, low
+  integer                        :: a,b,c,d
+  integer(bit_kind)              :: hole, particle, tmp
+  double precision, parameter    :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
+  
+  ASSERT (Nint > 0)
+  nperm = 0
+  exc(0,1) = 0
+  exc(0,2) = 0
+  
+  idx_particle = 0
+  idx_hole = 0
+  ishift = 1-bit_kind_size
+  do l=1,Nint
+    ishift = ishift + bit_kind_size
+    if (det1(l) == det2(l)) then
+      cycle
+    endif
+    tmp = xor( det1(l), det2(l) )
+    particle = iand(tmp, det2(l))
+    hole     = iand(tmp, det1(l))
+    do while (particle /= 0_bit_kind)
+      tz = trailz(particle)
+      idx_particle = idx_particle + 1
+      exc(0,2) = exc(0,2) + 1
+      exc(idx_particle,2) = tz+ishift
+      particle = iand(particle,particle-1_bit_kind)
+    enddo
+    if (iand(exc(0,1),exc(0,2))==2) then  ! exc(0,1)==2 or exc(0,2)==2
+      exit
+    endif
+    do while (hole /= 0_bit_kind)
+      tz = trailz(hole)
+      idx_hole = idx_hole + 1
+      exc(0,1) = exc(0,1) + 1
+      exc(idx_hole,1) = tz+ishift
+      hole = iand(hole,hole-1_bit_kind)
+    enddo
+    if (iand(exc(0,1),exc(0,2))==2) then ! exc(0,1)==2 or exc(0,2)==2
+      exit
+    endif
+  enddo
+  
+  select case (exc(0,1))
+      
+    case(1)
+      low  = min(exc(1,1), exc(1,2))
+      high = max(exc(1,1), exc(1,2))
+      
+      ASSERT (low > 0)
+      j = ishft(low-1,-bit_kind_shift)+1   ! Find integer in array(Nint)
+      n = iand(low-1,bit_kind_size-1)+1        ! mod(low,bit_kind_size)
+      ASSERT (high > 0)
+      k = ishft(high-1,-bit_kind_shift)+1
+      m = iand(high-1,bit_kind_size-1)+1
+      
+      if (j==k) then
+        nperm = nperm + popcnt(iand(det1(j),                         &
+            iand( ibset(0_bit_kind,m-1)-1_bit_kind,                  &
+            ibclr(-1_bit_kind,n)+1_bit_kind ) ))
+      else
+        nperm = nperm + popcnt(iand(det1(k),                         &
+            ibset(0_bit_kind,m-1)-1_bit_kind))
+        if (n < bit_kind_size) then
+          nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
+        endif
+        do i=j+1,k-1
+          nperm = nperm + popcnt(det1(i))
+        end do
+      endif
+      
+    case (2)
+      
+      do i=1,2
+        low  = min(exc(i,1), exc(i,2))
+        high = max(exc(i,1), exc(i,2))
+        
+        ASSERT (low > 0)
+        j = ishft(low-1,-bit_kind_shift)+1   ! Find integer in array(Nint)
+        n = iand(low-1,bit_kind_size-1)+1        ! mod(low,bit_kind_size)
+        ASSERT (high > 0)
+        k = ishft(high-1,-bit_kind_shift)+1
+        m = iand(high-1,bit_kind_size-1)+1
+        
+        if (j==k) then
+          nperm = nperm + popcnt(iand(det1(j),                       &
+              iand( ibset(0_bit_kind,m-1)-1_bit_kind,                &
+              ibclr(-1_bit_kind,n)+1_bit_kind ) ))
+        else
+          nperm = nperm + popcnt(iand(det1(k),                       &
+              ibset(0_bit_kind,m-1)-1_bit_kind))
+          if (n < bit_kind_size) then
+            nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
+          endif
+          do l=j+1,k-1
+            nperm = nperm + popcnt(det1(l))
+          end do
+        endif
+        
+      enddo
+      
+      a = min(exc(1,1), exc(1,2))
+      b = max(exc(1,1), exc(1,2))
+      c = min(exc(2,1), exc(2,2))
+      d = max(exc(2,1), exc(2,2))
+      if (c>a .and. c<b .and. d>b) then
+        nperm = nperm + 1
+      endif
+  end select
+  
+  phase = phase_dble(iand(nperm,1))
+  
+end
+
+subroutine get_mono_excitation_spin(det1,det2,exc,phase,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns the excitation operator between two singly excited determinants and the phase
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: det1(Nint)
+  integer(bit_kind), intent(in)  :: det2(Nint)
+  integer, intent(out)           :: exc(0:2,2)
+  double precision, intent(out)  :: phase
+  integer                        :: tz
+  integer                        :: l, idx_hole, idx_particle, ishift
+  integer                        :: nperm
+  integer                        :: i,j,k,m,n
+  integer                        :: high, low
+  integer                        :: a,b,c,d
+  integer(bit_kind)              :: hole, particle, tmp
+  double precision, parameter    :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
+  
+  ASSERT (Nint > 0)
+  nperm = 0
+  exc(0,1) = 0
+  exc(0,2) = 0
+  
+  ishift = 1-bit_kind_size
+  do l=1,Nint
+    ishift = ishift + bit_kind_size
+    if (det1(l) == det2(l)) then
+      cycle
+    endif
+    tmp = xor( det1(l), det2(l) )
+    particle = iand(tmp, det2(l))
+    hole     = iand(tmp, det1(l))
+    if (particle /= 0_bit_kind) then
+      tz = trailz(particle)
+      exc(0,2) = 1
+      exc(1,2) = tz+ishift
+    endif
+    if (hole /= 0_bit_kind) then
+      tz = trailz(hole)
+      exc(0,1) = 1
+      exc(1,1) = tz+ishift
+    endif
+    
+    if ( iand(exc(0,1),exc(0,2)) /= 1) then  ! exc(0,1)/=1 and exc(0,2) /= 1
+      cycle
+    endif
+    
+    low = min(exc(1,1),exc(1,2))
+    high = max(exc(1,1),exc(1,2))
+    
+    ASSERT (low > 0)
+    j = ishft(low-1,-bit_kind_shift)+1   ! Find integer in array(Nint)
+    n = iand(low-1,bit_kind_size-1)+1      ! mod(low,bit_kind_size)
+    ASSERT (high > 0)
+    k = ishft(high-1,-bit_kind_shift)+1
+    m = iand(high-1,bit_kind_size-1)+1
+    if (j==k) then
+      nperm = popcnt(iand(det1(j),                                   &
+          iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind)))
+    else
+      nperm = nperm + popcnt(iand(det1(k),ibset(0_bit_kind,m-1)-1_bit_kind))
+      if (n < bit_kind_size) then
+        nperm = nperm + popcnt(iand(det1(j),ibclr(-1_bit_kind,n)+1_bit_kind))
+      endif
+      do i=j+1,k-1
+        nperm = nperm + popcnt(det1(i))
+      end do
+    endif
+    phase = phase_dble(iand(nperm,1))
+    return
+    
+  enddo
+end
+
+subroutine i_H_j_mono_spin(key_i,key_j,Nint,spin,hij)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants differing by a single excitation
+  END_DOC
+  integer, intent(in)            :: Nint, spin
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij
+  
+  integer                        :: exc(0:2,2)
+  double precision               :: phase
+
+  PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
+
+  call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint)
+  call get_mono_excitation_from_fock(key_i,key_j,exc(1,2),exc(1,1),spin,phase,hij)
+end
+
+subroutine i_H_j_double_spin(key_i,key_j,Nint,hij)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants differing by a same-spin double excitation
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint), key_j(Nint)
+  double precision, intent(out)  :: hij
+  
+  integer                        :: exc(0:2,2)
+  double precision               :: phase
+  double precision, external     :: get_mo_bielec_integral
+
+  PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
+
+  call get_double_excitation_spin(key_i,key_j,exc,phase,Nint)
+  hij = phase*(get_mo_bielec_integral(                               &
+      exc(1,1),                                                      &
+      exc(2,1),                                                      &
+      exc(1,2),                                                      &
+      exc(2,2), mo_integrals_map) -                                  &
+      get_mo_bielec_integral(                                        &
+      exc(1,1),                                                      &
+      exc(2,1),                                                      &
+      exc(2,2),                                                      &
+      exc(1,2), mo_integrals_map) )
+end
+
+subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants differing by an opposite-spin double excitation
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij
+  
+  integer                        :: exc(0:2,2,2)
+  double precision               :: phase
+  double precision, external     :: get_mo_bielec_integral
+
+  PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
+
+  call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint)
+  call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase,Nint)
+  if (exc(1,1,1) == exc(1,2,2)) then
+    hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
+  else if (exc(1,2,1) == exc(1,1,2)) then
+    hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
+  else
+    hij = phase*get_mo_bielec_integral(                              &
+        exc(1,1,1),                                                  &
+        exc(1,1,2),                                                  &
+        exc(1,2,1),                                                  &
+        exc(1,2,2) ,mo_integrals_map)
+  endif
+end
+
+
diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f
index acc49d50..4bb35979 100644
--- a/src/Determinants/spindeterminants.irp.f
+++ b/src/Determinants/spindeterminants.irp.f
@@ -386,8 +386,9 @@ END_PROVIDER
 !==============================================================================!
 
 BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ]
-&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows, (N_det) ]
+&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
+&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order  , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -395,19 +396,20 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
 !  D_a^t C D_b
 !
 ! Rows are alpha determinants and columns are beta.
+!
+! Order refers to psi_det
   END_DOC
   integer                        :: i,j,k, l
   integer(bit_kind)               :: tmp_det(N_int,2)
-  integer                        :: idx
   integer, external              :: get_index_in_psi_det_sorted_bit
 
 
   PROVIDE psi_coef_sorted_bit
 
-  integer, allocatable :: iorder(:), to_sort(:)
+  integer, allocatable :: to_sort(:)
   integer, external :: get_index_in_psi_det_alpha_unique
   integer, external :: get_index_in_psi_det_beta_unique
-  allocate(iorder(N_det), to_sort(N_det))
+  allocate(to_sort(N_det))
   do k=1,N_det
     i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int)
     j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int)
@@ -418,36 +420,41 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
     psi_bilinear_matrix_rows(k) = i
     psi_bilinear_matrix_columns(k) = j
     to_sort(k) = N_det_alpha_unique * (j-1) + i
-    iorder(k) = k
+    psi_bilinear_matrix_order(k) = k
   enddo
-  call isort(to_sort, iorder, N_det)
-  call iset_order(psi_bilinear_matrix_rows,iorder,N_det)
-  call iset_order(psi_bilinear_matrix_columns,iorder,N_det)
+  call isort(to_sort, psi_bilinear_matrix_order, N_det)
+  call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det)
+  call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det)
   do l=1,N_states
-    call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det)
+    call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det)
   enddo
-  deallocate(iorder,to_sort)
+  deallocate(to_sort)
 END_PROVIDER
 
 
 BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
-&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows, (N_det) ]
+&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
+&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order  , (N_det) ]
+&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
 ! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
 !  D_a^t C D_b
 !
-! Rows are Beta determinants and columns are alpha
+! Rows are Alpha determinants and columns are beta, but the matrix is stored in row major
+! format
+!
+! Order refers to psi_bilinear_matrix
   END_DOC
   integer                        :: i,j,k,l
 
 
   PROVIDE psi_coef_sorted_bit
 
-  integer, allocatable :: iorder(:), to_sort(:)
-  allocate(iorder(N_det), to_sort(N_det))
+  integer, allocatable :: to_sort(:)
+  allocate(to_sort(N_det))
   do l=1,N_states
     do k=1,N_det
       psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l)
@@ -459,15 +466,18 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
     i = psi_bilinear_matrix_transp_columns(k) 
     j = psi_bilinear_matrix_transp_rows   (k)
     to_sort(k) = N_det_beta_unique * (j-1) + i
-    iorder(k) = k
+    psi_bilinear_matrix_transp_order(k) = k
   enddo
-  call isort(to_sort, iorder, N_det)
-  call iset_order(psi_bilinear_matrix_transp_rows,iorder,N_det)
-  call iset_order(psi_bilinear_matrix_transp_columns,iorder,N_det)
+  call isort(to_sort, psi_bilinear_matrix_transp_order, N_det)
+  call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det)
+  call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det)
   do l=1,N_states
-    call dset_order(psi_bilinear_matrix_transp_values(1,l),iorder,N_det)
+    call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det)
   enddo
-  deallocate(iorder,to_sort)
+  do k=1,N_det
+    psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_transp_order(k)) = k
+  enddo
+  deallocate(to_sort)
 END_PROVIDER
 
 
@@ -552,7 +562,7 @@ subroutine generate_all_alpha_beta_det_products
 !  Create a wave function from all possible alpha x beta determinants
   END_DOC
   integer                        :: i,j,k,l
-  integer                        :: idx, iproc
+  integer                        :: iproc
   integer, external              :: get_index_in_psi_det_sorted_bit
   integer(bit_kind), allocatable :: tmp_det(:,:,:)
   logical, external              :: is_in_wavefunction
@@ -561,7 +571,7 @@ subroutine generate_all_alpha_beta_det_products
   !$OMP PARALLEL DEFAULT(NONE) SHARED(psi_coef_sorted_bit,N_det_beta_unique,&
       !$OMP N_det_alpha_unique, N_int, psi_det_alpha_unique, psi_det_beta_unique,&
       !$OMP N_det)                                                &
-      !$OMP PRIVATE(i,j,k,l,tmp_det,idx,iproc)
+      !$OMP PRIVATE(i,j,k,l,tmp_det,iproc)
   !$ iproc = omp_get_thread_num()
   allocate (tmp_det(N_int,2,N_det_alpha_unique))
   !$OMP DO
@@ -586,3 +596,782 @@ subroutine generate_all_alpha_beta_det_products
 end
 
 
+
+
+subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buffer, singles, doubles, n_singles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single and double excitations in the list of
+! unique alpha determinants.
+!
+! /!\ : The buffer is transposed !
+!
+  END_DOC
+  integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(Nint)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_singles
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  select case (Nint)
+    case (1)
+      call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles)
+      return
+!    case (2)
+!      call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+!      return
+    case (3)
+      call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+      return
+  end select
+
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
+
+  do k=1,Nint
+    do i=1,size_buffer
+      xorvec(i, k) = xor( spindet(k), buffer(k,i) )
+    enddo
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  do k=2,Nint
+    !DIR$ VECTOR ALIGNED
+    do i=1,size_buffer
+      if ( (degree(i) <= 4).and.(xorvec(i,k) /= 0_8) ) then
+        degree(i) = degree(i) + popcnt(xorvec(i,k))
+      endif
+    enddo
+  enddo
+
+  n_singles = 1
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_singles(buffer, idx, spindet, Nint, size_buffer, singles, n_singles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(Nint)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: n_singles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  select case (Nint)
+    case (1)
+      call get_all_spin_singles_1(buffer, idx, spindet(1), size_buffer, singles, n_singles)
+      return
+    case (2)
+      call get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
+      return
+    case (3)
+      call get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
+      return
+  end select
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
+
+  do k=1,Nint
+    do i=1,size_buffer
+      xorvec(i, k) = xor( spindet(k), buffer(k,i) )
+    enddo
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  do k=2,Nint
+    !DIR$ VECTOR ALIGNED
+    do i=1,size_buffer
+      if ( (degree(i) <= 2).and.(xorvec(i,k) /= 0_8) ) then
+        degree(i) = degree(i) + popcnt(xorvec(i,k))
+      endif
+    enddo
+  enddo
+
+  n_singles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_doubles(buffer, idx, spindet, Nint, size_buffer, doubles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the double excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: Nint, size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(Nint,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(Nint)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  select case (Nint)
+    case (1)
+      call get_all_spin_doubles_1(buffer, idx, spindet(1), size_buffer, doubles, n_doubles)
+      return
+    case (2)
+      call get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
+      return
+    case (3)
+      call get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
+      return
+  end select
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, Nint), degree(size_buffer) )
+
+  do k=1,Nint
+    do i=1,size_buffer
+      xorvec(i, k) = xor( spindet(k), buffer(k,i) )
+    enddo
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  do k=2,Nint
+    !DIR$ VECTOR ALIGNED
+    do i=1,size_buffer
+      if ( (degree(i) <= 4).and.(xorvec(i,k) /= 0_8) ) then
+        degree(i) = degree(i) + popcnt(xorvec(i,k))
+      endif
+    enddo
+  enddo
+
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+  enddo
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single and double excitations in the list of
+! unique alpha determinants.
+!
+! /!\ : The buffer is transposed !
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer
+  integer, intent(in)            :: idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(size_buffer)
+  integer(bit_kind), intent(in)  :: spindet
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_singles
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:)
+  integer                        :: degree
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align) )
+
+  do i=1,size_buffer
+    xorvec(i) = xor( spindet, buffer(i) )
+  enddo
+  
+  n_singles = 1
+  n_doubles = 1
+
+  do i=1,size_buffer
+    degree = popcnt(xorvec(i))
+    if ( degree == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+    if ( degree == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  n_doubles = n_doubles-1
+
+  deallocate(xorvec)
+end
+
+
+subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_singles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(size_buffer)
+  integer(bit_kind), intent(in)  :: spindet
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: n_singles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:)
+
+  allocate( xorvec(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i) = xor( spindet, buffer(i) )
+  enddo
+  
+  n_singles = 1
+  do i=1,size_buffer
+    if ( popcnt(xorvec(i)) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the double excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(size_buffer)
+  integer(bit_kind), intent(in)  :: spindet
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:)
+
+  integer, external              :: align_double
+
+  allocate( xorvec(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i) = xor( spindet, buffer(i) )
+  enddo
+  
+  n_doubles = 1
+
+  do i=1,size_buffer
+    if ( popcnt(xorvec(i)) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+  enddo
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single and double excitations in the list of
+! unique alpha determinants.
+!
+! /!\ : The buffer is transposed !
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(2,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(2)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_singles
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+
+  n_singles = 1
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(2,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(2)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: n_singles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 2).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+
+  n_singles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the double excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(2,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(2)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 2), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+  enddo
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+subroutine get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single and double excitations in the list of
+! unique alpha determinants.
+!
+! /!\ : The buffer is transposed !
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(3,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(3)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_singles
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+    xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,3) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,3))
+    endif
+  enddo
+
+  n_singles = 1
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i) 
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the single excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(3,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(3)
+  integer, intent(out)           :: singles(size_buffer)
+  integer, intent(out)           :: n_singles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+    xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 2).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 2).and.(xorvec(i,3) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,3))
+    endif
+  enddo
+
+  n_singles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 2 ) then
+      singles(n_singles) = idx(i)
+      n_singles = n_singles+1
+    endif
+  enddo
+  n_singles = n_singles-1
+  deallocate(xorvec)
+  
+end
+
+
+subroutine get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+!
+! Returns the indices of all the double excitations in the list of
+! unique alpha determinants.
+!
+  END_DOC
+  integer, intent(in)            :: size_buffer, idx(size_buffer)
+  integer(bit_kind), intent(in)  :: buffer(3,size_buffer)
+  integer(bit_kind), intent(in)  :: spindet(3)
+  integer, intent(out)           :: doubles(size_buffer)
+  integer, intent(out)           :: n_doubles
+
+  integer                        :: i,k
+  integer(bit_kind), allocatable :: xorvec(:,:)
+  integer,           allocatable :: degree(:)
+  integer                        :: size_buffer_align
+
+  integer, external              :: align_double
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
+
+  size_buffer_align = align_double(size_buffer)
+  allocate( xorvec(size_buffer_align, 3), degree(size_buffer) )
+
+  do i=1,size_buffer
+    xorvec(i, 1) = xor( spindet(1), buffer(1,i) )
+    xorvec(i, 2) = xor( spindet(2), buffer(2,i) )
+    xorvec(i, 3) = xor( spindet(3), buffer(3,i) )
+  enddo
+  
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if (xorvec(i,1) /= 0_8) then
+      degree(i) = popcnt(xorvec(i,1))
+    else
+      degree(i) = 0
+    endif
+  enddo
+
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,2) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,2))
+    endif
+  enddo
+  !DIR$ VECTOR ALIGNED
+  do i=1,size_buffer
+    if ( (degree(i) <= 4).and.(xorvec(i,3) /= 0_8) ) then
+      degree(i) = degree(i) + popcnt(xorvec(i,3))
+    endif
+  enddo
+
+  n_doubles = 1
+  do i=1,size_buffer
+    if ( degree(i) == 4 ) then
+      doubles(n_doubles) = idx(i)
+      n_doubles = n_doubles+1
+    endif
+  enddo
+  n_doubles = n_doubles-1
+  deallocate(xorvec)
+  
+end
+
diff --git a/src/Determinants/usefull_for_ovb.irp.f b/src/Determinants/useful_for_ovb.irp.f
similarity index 97%
rename from src/Determinants/usefull_for_ovb.irp.f
rename to src/Determinants/useful_for_ovb.irp.f
index 7b89897b..25bdb03a 100644
--- a/src/Determinants/usefull_for_ovb.irp.f
+++ b/src/Determinants/useful_for_ovb.irp.f
@@ -2,7 +2,8 @@
 integer function n_open_shell(det_in,nint)
  implicit none
  use bitmasks
- integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer, intent(in) :: nint
+ integer(bit_kind), intent(in) :: det_in(nint,2)
  integer :: i
  n_open_shell = 0
  do i=1,Nint
@@ -13,7 +14,8 @@ end
 integer function n_closed_shell(det_in,nint)
  implicit none
  use bitmasks
- integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer, intent(in) :: nint
+ integer(bit_kind), intent(in) :: det_in(nint,2)
  integer :: i
  n_closed_shell = 0
  do i=1,Nint
@@ -24,7 +26,8 @@ end
 integer function n_closed_shell_cas(det_in,nint)
  implicit none
  use bitmasks
- integer(bit_kind), intent(in) :: det_in(nint,2),nint
+ integer, intent(in) :: nint
+ integer(bit_kind), intent(in) :: det_in(nint,2)
  integer(bit_kind) :: det_tmp(nint,2)
  integer :: i
  n_closed_shell_cas = 0
diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f
index 1f2a7a1b..82b89f22 100644
--- a/src/Integrals_Bielec/map_integrals.irp.f
+++ b/src/Integrals_Bielec/map_integrals.irp.f
@@ -44,8 +44,8 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
   l(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i2)+1.d0)-1.d0))
   i3   = i1 - ishft(i2*i2-i2,-1)
   k(1) = ceiling(0.5d0*(dsqrt(8.d0*dble(i3)+1.d0)-1.d0))
-  j(1) = i2 - ishft(l(1)*l(1)-l(1),-1)
-  i(1) = i3 - ishft(k(1)*k(1)-k(1),-1)
+  j(1) = int(i2 - ishft(l(1)*l(1)-l(1),-1),4)
+  i(1) = int(i3 - ishft(k(1)*k(1)-k(1),-1),4)
 
               !ijkl
   i(2) = i(1) !ilkj
diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
index 186938c2..bfe10b91 100644
--- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
+++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
@@ -51,7 +51,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
   print*, 'Providing the nuclear electron pseudo integrals (local)'
   
   call wall_time(wall_1)
-  wall_0 = wall_1
   call cpu_time(cpu_1)
   
   thread_num = 0
@@ -66,6 +65,8 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
       !$OMP         wall_1)
   
   !$ thread_num = omp_get_thread_num()
+
+  wall_0 = wall_1
   !$OMP DO SCHEDULE (guided)
   
   do j = 1, ao_num
@@ -148,7 +149,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
   print*, 'Providing the nuclear electron pseudo integrals (non-local)'
   
   call wall_time(wall_1)
-  wall_0 = wall_1
   call cpu_time(cpu_1)
   thread_num = 0
 
@@ -164,6 +164,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
   
   !$ thread_num = omp_get_thread_num()
   
+  wall_0 = wall_1
   !$OMP DO SCHEDULE (guided)
 !  
   do j = 1, ao_num
diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f
index 95a771b0..48341129 100644
--- a/src/MO_Basis/ao_ortho_canonical.irp.f
+++ b/src/MO_Basis/ao_ortho_canonical.irp.f
@@ -42,7 +42,7 @@
         9;;
       END_TEMPLATE
       case default
-        stop 'Error in ao_cart_to_sphe'
+        stop 'Error in ao_cart_to_sphe : angular momentum too high'
     end select
   enddo
 
diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f
index 0f338877..750e3420 100644
--- a/src/MO_Basis/utils.irp.f
+++ b/src/MO_Basis/utils.irp.f
@@ -88,7 +88,7 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
     enddo
   endif
   do i=1,m
-    write (output_mo_basis,'(I8,X,F16.10)')  i,eigvalues(i)
+    write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
   enddo
   write (output_mo_basis,'(A)') '======== ================'
   write (output_mo_basis,'(A)')  ''
@@ -135,7 +135,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   write (output_mo_basis,'(A)')  '======== ================'
 
   do i=1,m
-    write (output_mo_basis,'(I8,X,F16.10)')  i,D(i)
+    write (output_mo_basis,'(I8,1X,F16.10)')  i,D(i)
   enddo
   write (output_mo_basis,'(A)')  '======== ================'
   write (output_mo_basis,'(A)')  ''
@@ -215,7 +215,7 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,
   write (output_mo_basis,'(A)')  ''
   write (output_mo_basis,'(A)')  '======== ================'
   do i = 1, m
-   write (output_mo_basis,'(I8,X,F16.10)')  i,eigvalues(i)
+   write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
   enddo
   write (output_mo_basis,'(A)')  '======== ================'
   write (output_mo_basis,'(A)')  ''
diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f
index b4da5fb1..34fae989 100644
--- a/src/Nuclei/nuclei.irp.f
+++ b/src/Nuclei/nuclei.irp.f
@@ -37,8 +37,8 @@ BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num_aligned,3) ]
    enddo
    deallocate(buffer)
    
-   character*(64), parameter      :: f = '(A16, 4(X,F12.6))'
-   character*(64), parameter      :: ft= '(A16, 4(X,A12  ))'
+   character*(64), parameter      :: f = '(A16, 4(1X,F12.6))'
+   character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
    double precision, parameter    :: a0= 0.529177249d0
    call write_time(output_Nuclei)
    write(output_Nuclei,'(A)') ''
diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f
index 47fb4355..32090f01 100644
--- a/src/Utils/LinearAlgebra.irp.f
+++ b/src/Utils/LinearAlgebra.irp.f
@@ -30,7 +30,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
   lwork = -1
   call dgesvd('A','A', m, n, A_tmp, LDA,                             &
       D, U, LDU, Vt, LDVt, work, lwork, info)
-  lwork = work(1)
+  lwork = int(work(1))
   deallocate(work)
 
   allocate(work(lwork))
@@ -153,7 +153,7 @@ subroutine ortho_qr(A,LDA,m,n)
   allocate (jpvt(n), tau(n), work(1))
   LWORK=-1
   call  dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
-  LWORK=2*WORK(1)
+  LWORK=2*int(WORK(1))
   deallocate(WORK)
   allocate(WORK(LWORK))
   call  dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
@@ -297,7 +297,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
     print *,  info, ': SVD failed'
     stop
   endif
-  lwork = work(1)
+  lwork = int(work(1))
   deallocate(work)
   allocate(work(lwork))
   call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info)
diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F
index 991ef80a..4974fd8e 100644
--- a/src/Utils/constants.include.F
+++ b/src/Utils/constants.include.F
@@ -1,5 +1,6 @@
 integer, parameter :: max_dim = 511
 integer, parameter :: SIMD_vector = 32
+integer, parameter :: N_int_max = 16
 
 double precision, parameter :: pi =  dacos(-1.d0)
 double precision, parameter :: sqpi =  dsqrt(dacos(-1.d0))
diff --git a/src/Utils/map_functions.irp.f b/src/Utils/map_functions.irp.f
index f5d6f4f8..0378c253 100644
--- a/src/Utils/map_functions.irp.f
+++ b/src/Utils/map_functions.irp.f
@@ -76,8 +76,8 @@ subroutine map_load_from_disk(filename,map)
   double precision               :: x
   type(c_ptr)                    :: c_pointer(3)
   integer                        :: fd(3)
-  integer*8                      :: i,k, l
-  integer                        :: n_elements, j
+  integer*8                      :: i,k,l
+  integer*4                      :: j,n_elements
 
 
 
@@ -105,14 +105,14 @@ subroutine map_load_from_disk(filename,map)
     map % map(i) % value => map % consolidated_value ( map % consolidated_idx (i+1) :)
     map % map(i) % key   => map % consolidated_key   ( map % consolidated_idx (i+1) :)
     map % map(i) % sorted = .True.
-    n_elements = map % consolidated_idx (i+2) - k
+    n_elements = int( map % consolidated_idx (i+2) - k, 4)
     k = map % consolidated_idx (i+2)
     map % map(i) % map_size = n_elements
     map % map(i) % n_elements = n_elements
     ! Load memory from disk
     do j=1,n_elements
       x = x + map % map(i) % value(j)
-      l = iand(l,map % map(i) % key(j))
+      l = iand(l,int(map % map(i) % key(j),8))
       if (map % map(i) % value(j) > 1.e30) then
         stop 'Error in integrals file'
       endif
diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90
index 80260233..ac16f97e 100644
--- a/src/Utils/map_module.f90
+++ b/src/Utils/map_module.f90
@@ -53,17 +53,17 @@ module map_module
 end module map_module
 
 
-real function map_mb(map)
+double precision function map_mb(map)
   use map_module
   use omp_lib
   implicit none
   type (map_type), intent(in)    :: map
   integer(map_size_kind)         :: i
   
-  map_mb = 8+map_size_kind+map_size_kind+omp_lock_kind+4
+  map_mb = dble(8+map_size_kind+map_size_kind+omp_lock_kind+4)
   do i=0,map%map_size
-    map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) +&
-        8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind
+    map_mb = map_mb + dble(map%map(i)%map_size*(cache_key_kind+integral_kind) +&
+        8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind)
   enddo
   map_mb = map_mb / (1024.d0*1024.d0)
 end
@@ -406,8 +406,8 @@ subroutine map_update(map, key, value, sze, thr)
               call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
               call cache_map_shrink(local_map,thr)
             endif
-            cache_key = iand(key(i),map_mask)
-            local_map%n_elements = local_map%n_elements + 1_8
+            cache_key = int(iand(key(i),map_mask),2)
+            local_map%n_elements = local_map%n_elements + 1
             local_map%value(local_map%n_elements) = value(i)
             local_map%key(local_map%n_elements) = cache_key
             local_map%sorted = .False.
@@ -464,7 +464,7 @@ subroutine map_append(map, key, value, sze)
     if (n_elements == map%map(idx_cache)%map_size) then
       call cache_map_reallocate(map%map(idx_cache), n_elements+ ishft(n_elements,-1))
     endif
-    cache_key = iand(key(i),map_mask)
+    cache_key = int(iand(key(i),map_mask),2)
     map%map(idx_cache)%value(n_elements) = value(i)
     map%map(idx_cache)%key(n_elements) = cache_key
     map%map(idx_cache)%n_elements = n_elements
@@ -615,7 +615,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
     idx = -1
     return
   endif
-  cache_key = iand(key,map_mask)
+  cache_key = int(iand(key,map_mask),2)
   ibegin = min(ibegin_in,sze)
   iend   = min(iend_in,sze)
   if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
@@ -723,7 +723,7 @@ subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in
     value = 0.d0
     return
   endif
-  cache_key = iand(key,map_mask)
+  cache_key = int(iand(key,map_mask),2)
   ibegin = min(ibegin_in,sze)
   iend   = min(iend_in,sze)
   if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f
index dd7fbc33..dc91ab3a 100644
--- a/src/Utils/sort.irp.f
+++ b/src/Utils/sort.irp.f
@@ -292,18 +292,17 @@ BEGIN_TEMPLATE
   ! contains the new order of the elements.
   ! iradix should be -1 in input.
   END_DOC
-  $int_type, intent(in)          :: isize
-  $int_type, intent(inout)       :: iorder(isize)
-  $type, intent(inout)           :: x(isize)
+  integer*$int_type, intent(in)  :: isize
+  integer*$int_type, intent(inout) :: iorder(isize)
+  integer*$type, intent(inout)   :: x(isize)
   integer, intent(in)            :: iradix
   integer                        :: iradix_new
-  $type, allocatable             :: x2(:), x1(:)
-  $type                          :: i4
-  $int_type, allocatable         :: iorder1(:),iorder2(:)
-  $int_type                      :: i0, i1, i2, i3, i
+  integer*$type, allocatable     :: x2(:), x1(:)
+  integer*$type                  :: i4
+  integer*$int_type, allocatable :: iorder1(:),iorder2(:)
+  integer*$int_type              :: i0, i1, i2, i3, i
   integer, parameter             :: integer_size=$octets
-  $type, parameter               :: zero=$zero
-  $type                          :: mask
+  integer*$type                  :: mask
   integer                        :: nthreads, omp_get_num_threads
   !DIR$ ATTRIBUTES ALIGN : 128   :: iorder1,iorder2, x2, x1
   
@@ -311,16 +310,16 @@ BEGIN_TEMPLATE
     
     ! Find most significant bit
     
-    i0 = 0_8
-    i4 = -1_8
+    i0 = 0_$int_type
+    i4 = -1_$type
     
     do i=1,isize
       i4 = max(i4,x(i))
     enddo
-    i3 = i4  ! Type conversion
+    i3 = int(i4,$int_type)
     
     iradix_new = integer_size-1-leadz(i3)
-    mask = ibset(zero,iradix_new)
+    mask = ibset(0_$type,iradix_new)
     nthreads = 1
     !   nthreads = 1+ishft(omp_get_num_threads(),-1)
     
@@ -331,22 +330,22 @@ BEGIN_TEMPLATE
       stop
     endif
     
-    i1=1_8
-    i2=1_8
+    i1=1_$int_type
+    i2=1_$int_type
     
     do i=1,isize
-      if (iand(mask,x(i)) == zero) then
+      if (iand(mask,x(i)) == 0_$type) then
         iorder1(i1) = iorder(i)
         x1(i1) = x(i)
-        i1 = i1+1_8
+        i1 = i1+1_$int_type
       else
         iorder2(i2) = iorder(i)
         x2(i2) = x(i)
-        i2 = i2+1_8
+        i2 = i2+1_$int_type
       endif
     enddo
-    i1=i1-1_8
-    i2=i2-1_8
+    i1=i1-1_$int_type
+    i2=i2-1_$int_type
     
     do i=1,i1
       iorder(i0+i) = iorder1(i)
@@ -399,12 +398,12 @@ BEGIN_TEMPLATE
   endif
   
   
-  mask = ibset(zero,iradix)
+  mask = ibset(0_$type,iradix)
   i0=1
   i1=1
   
   do i=1,isize
-    if (iand(mask,x(i)) == zero) then
+    if (iand(mask,x(i)) == 0_$type) then
       iorder(i0) = iorder(i)
       x(i0) = x(i)
       i0 = i0+1
@@ -443,12 +442,12 @@ BEGIN_TEMPLATE
   
  end
 
-SUBST [ X, type, octets, is_big, big, int_type, zero ]
- i  ; integer   ; 32 ; .False. ; ; integer ; 0;;
- i8 ; integer*8 ; 32 ; .False. ; ; integer ; 0_8;;
- i2 ; integer*2 ; 32 ; .False. ; ; integer ; 0;;
- i  ; integer   ; 64 ; .True. ; _big ; integer*8 ; 0 ;;
- i8 ; integer*8 ; 64 ; .True. ; _big ; integer*8 ; 0_8 ;;
+SUBST [ X, type, octets, is_big, big, int_type ]
+ i  ; 4 ; 32 ; .False. ;      ; 4 ;;
+ i8 ; 8 ; 32 ; .False. ;      ; 4 ;;
+ i2 ; 2 ; 32 ; .False. ;      ; 4 ;;
+ i  ; 4 ; 64 ; .True.  ; _big ; 8 ;;
+ i8 ; 8 ; 64 ; .True.  ; _big ; 8 ;;
 END_TEMPLATE
 
 
diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f
index 8e3a94e5..91ed9200 100644
--- a/src/ZMQ/utils.irp.f
+++ b/src/ZMQ/utils.irp.f
@@ -1,11 +1,8 @@
 use f77_zmq
 use omp_lib
 
-integer, pointer :: thread_id 
-integer(omp_lock_kind) :: zmq_lock
-
-
-BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
+ BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
+&BEGIN_PROVIDER [ integer(omp_lock_kind), zmq_lock ]
   use f77_zmq
   implicit none
   BEGIN_DOC
@@ -407,7 +404,9 @@ subroutine end_zmq_sub_socket(zmq_socket_sub)
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_sub
   integer                        :: rc
   
+  call omp_set_lock(zmq_lock)
   rc = f77_zmq_close(zmq_socket_sub)
+  call omp_unset_lock(zmq_lock)
   if (rc /= 0) then
     print *,  'f77_zmq_close(zmq_socket_sub)'
     stop 'error'
@@ -426,7 +425,9 @@ subroutine end_zmq_pair_socket(zmq_socket_pair)
   integer                        :: rc
   character*(8), external        :: zmq_port
   
+  call omp_set_lock(zmq_lock)
   rc = f77_zmq_close(zmq_socket_pair)
+  call omp_unset_lock(zmq_lock)
   if (rc /= 0) then
     print *,  'f77_zmq_close(zmq_socket_pair)'
     stop 'error'
@@ -444,7 +445,9 @@ subroutine end_zmq_pull_socket(zmq_socket_pull)
   integer                        :: rc
   character*(8), external        :: zmq_port
   
+  call omp_set_lock(zmq_lock)
   rc = f77_zmq_close(zmq_socket_pull)
+  call omp_unset_lock(zmq_lock)
   if (rc /= 0) then
     print *,  'f77_zmq_close(zmq_socket_pull)'
     stop 'error'
@@ -469,7 +472,9 @@ subroutine end_zmq_push_socket(zmq_socket_push,thread)
     stop 'Unable to set ZMQ_LINGER on push socket'
   endif
 
+  call omp_set_lock(zmq_lock)
   rc = f77_zmq_close(zmq_socket_push)
+  call omp_unset_lock(zmq_lock)
   if (rc /= 0) then
     print *,  'f77_zmq_close(zmq_socket_push)'
     stop 'error'
@@ -500,10 +505,17 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,name_in)
   integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
   integer(ZMQ_PTR), intent(out)  :: zmq_to_qp_run_socket
 
+  call omp_set_lock(zmq_lock)
   zmq_context = f77_zmq_ctx_new ()
+  call omp_unset_lock(zmq_lock)
   if (zmq_context == 0_ZMQ_PTR) then
      stop 'ZMQ_PTR is null'
   endif
+!  rc = f77_zmq_ctx_set(zmq_context, ZMQ_IO_THREADS, nproc)
+!  if (rc /= 0) then
+!    print *,  'Unable to set the number of ZMQ IO threads to', nproc
+!  endif
+
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   name = name_in
   sze = len(trim(name))
@@ -584,7 +596,10 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,name_in)
   zmq_state = 'No_state'
   call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
 
+  call omp_set_lock(zmq_lock)
   rc = f77_zmq_ctx_term(zmq_context)
+  zmq_context = 0_ZMQ_PTR
+  call omp_unset_lock(zmq_lock)
   if (rc /= 0) then
     print *,  'Unable to terminate ZMQ context'
     stop 'error'
diff --git a/tests/bats/cassd.bats b/tests/bats/cassd.bats
index 1b845e91..67c35235 100644
--- a/tests/bats/cassd.bats
+++ b/tests/bats/cassd.bats
@@ -13,7 +13,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
   qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" 
   qp_run cassd_zmq $INPUT  
   energy="$(ezfio get cas_sd_zmq energy_pt2)"
-  eq $energy -76.231084536315 5.E-5
+  eq $energy -76.231248286858 5.E-5
 
   ezfio set determinants n_det_max 1024 
   ezfio set determinants read_wf True
@@ -21,6 +21,6 @@ source $QP_ROOT/tests/bats/common.bats.sh
   qp_run cassd_zmq $INPUT  
   ezfio set determinants read_wf False
   energy="$(ezfio get cas_sd_zmq energy)"
-  eq $energy -76.2225863580749   2.E-5
+  eq $energy -76.2225678834779   2.E-5
 }
 
diff --git a/tests/bats/fci.bats b/tests/bats/fci.bats
index 79ff91ab..6cded581 100644
--- a/tests/bats/fci.bats
+++ b/tests/bats/fci.bats
@@ -42,11 +42,13 @@ function run_FCI_ZMQ() {
   qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
 }
 @test "FCI H2O cc-pVDZ" {
-  run_FCI h2o.ezfio 2000  -0.761255633582109E+02 -0.761258377850042E+02
+  run_FCI h2o.ezfio 2000  -76.1253758241716     -76.1258130146102     
 }
 
+
+
 @test "FCI-ZMQ H2O cc-pVDZ" {
-  run_FCI_ZMQ h2o.ezfio 2000  -0.761255633582109E+02 -0.761258377850042E+02
+  run_FCI_ZMQ h2o.ezfio 2000 -76.1250552686394     -76.1258817228809  
 }
 
 
diff --git a/tests/bats/mrcepa0.bats b/tests/bats/mrcepa0.bats
index 2420955c..9a62885e 100644
--- a/tests/bats/mrcepa0.bats
+++ b/tests/bats/mrcepa0.bats
@@ -32,7 +32,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
   ezfio set mrcepa0 n_it_max_dressed_ci 3
   qp_run $EXE $INPUT  
   energy="$(ezfio get mrcepa0 energy_pt2)"
-  eq $energy -76.2381673136696 2.e-4
+  eq $energy -76.2381754078899   1.e-4
 }
 
 @test "MRSC2 H2O cc-pVDZ" {