diff --git a/.gitignore b/.gitignore
index 096e385b..085e4f74 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,7 +5,10 @@ build.ninja
 .ninja_deps
 bin/
 lib/
+lib64/
+libexec/
 config/qp_create_ninja.pickle
 src/*/.gitignore
 ezfio_interface.irp.f
 share
+*.swp
diff --git a/external/.gitignore b/external/.gitignore
index 676c79b7..241e560d 100644
--- a/external/.gitignore
+++ b/external/.gitignore
@@ -1,2 +1,2 @@
-#*
+*
 
diff --git a/include/.gitignore b/include/.gitignore
index afcb2de6..72e8ffc0 100644
--- a/include/.gitignore
+++ b/include/.gitignore
@@ -1,8 +1 @@
-zmq.h
-gmp.h
-zconf.h
-zconf.h
-zlib.h
-zmq_utils.h
-f77_zmq_free.h
-f77_zmq.h
+*
diff --git a/ocaml/.gitignore b/ocaml/.gitignore
deleted file mode 100644
index 888fa10e..00000000
--- a/ocaml/.gitignore
+++ /dev/null
@@ -1,57 +0,0 @@
-_build
-element_create_db
-element_create_db.byte
-ezfio.ml
-.gitignore
-Git.ml
-Input_ao_one_e_ints.ml
-Input_ao_two_e_erf_ints.ml
-Input_ao_two_e_ints.ml
-Input_auto_generated.ml
-Input_becke_numerical_grid.ml
-Input_cassd.ml
-Input_champ.ml
-Input_cipsi.ml
-Input_cipsi_tc_bi_ortho.ml
-Input_cosgtos_ao_int.ml
-Input_davidson_keywords.ml
-Input_davidson.ml
-Input_density_for_dft.ml
-Input_determinants.ml
-Input_dft_keywords.ml
-Input_dressing.ml
-Input_fci_tc_bi.ml
-Input_general_mrci.ml
-Input_mo_localization.ml
-Input_mo_one_e_ints.ml
-Input_mo_two_e_erf_ints.ml
-Input_mo_two_e_ints.ml
-Input_mpn.ml
-Input_mu_of_r.ml
-Input_nuclei.ml
-Input_perturbation.ml
-Input_pseudo.ml
-Input_qmcchem.ml
-Input_scf_utils.ml
-Input_some_mu_of_r.ml
-Input_tc_keywords.ml
-Input_three_body_ints.ml
-Input_trust_region.ml
-Input_two_body_rdm.ml
-Input_utils.ml
-Input_utils_trust_region.ml
-qp_create_ezfio
-qp_create_ezfio.native
-qp_edit
-qp_edit.ml
-qp_edit.native
-qp_print_basis
-qp_print_basis.native
-qp_run
-qp_run.native
-qp_set_mo_class
-qp_set_mo_class.native
-qp_tunnel
-qp_tunnel.native
-qptypes_generator.byte
-Qptypes.ml
diff --git a/ocaml/Makefile b/ocaml/Makefile
index 40d292fe..8853a991 100644
--- a/ocaml/Makefile
+++ b/ocaml/Makefile
@@ -29,7 +29,7 @@ tests: $(ALL_TESTS)
 .gitignore: $(MLFILES) $(MLIFILES)
 	@for i in .gitignore ezfio.ml element_create_db Qptypes.ml Git.ml *.byte *.native _build $(ALL_EXE) $(ALL_TESTS) \
 	$(patsubst %.ml,%,$(wildcard test_*.ml)) $(patsubst %.ml,%,$(wildcard qp_*.ml)) \
-	$(shell grep Input Input_auto_generated.ml | awk '{print $$2 ".ml"}') \
+	Input_*.ml \
 	qp_edit.ml qp_edit qp_edit.native Input_auto_generated.ml;\
 	do \
 	   echo $$i ; \
diff --git a/scripts/.gitignore b/scripts/.gitignore
index b44ac5a2..103b3ae9 100644
--- a/scripts/.gitignore
+++ b/scripts/.gitignore
@@ -2,3 +2,6 @@
 *.pyo
 docopt.py
 resultsFile/
+verif_omp/a.out
+src/*/Makefile
+src/*/*/
diff --git a/src/.gitignore b/src/.gitignore
new file mode 100644
index 00000000..6353c21a
--- /dev/null
+++ b/src/.gitignore
@@ -0,0 +1,11 @@
+*
+!README.rst
+!*/
+*/*
+!*/*.*
+*/*.o
+*/build.ninja
+*/ezfio_interface.irp.f
+*/.gitignore
+*/*.swp
+
diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f
index 6949fc48..dc19f6c7 100644
--- a/src/ao_one_e_ints/pot_ao_ints.irp.f
+++ b/src/ao_one_e_ints/pot_ao_ints.irp.f
@@ -38,11 +38,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
 
       ao_integrals_n_e = 0.d0
 
-      !        _
-      ! /|  / |_)
-      !  | /  | \
-      !
-
       !$OMP PARALLEL                                                   &
           !$OMP DEFAULT (NONE)                                         &
           !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
@@ -106,7 +101,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
     endif
 
 
-    IF(DO_PSEUDO) THEN
+    IF(do_pseudo) THEN
        ao_integrals_n_e += ao_pseudo_integrals
     ENDIF
 
diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 7e8020f5..8d436c19 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -292,9 +292,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
 !      call omp_set_max_active_levels(1)
 
 
-      print '(A)', '========== ======================= ===================== ===================== ==========='
-      print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+      print '(A)', ' Samples           Energy             PT2            Variance          Norm^2       Convergence    Seconds'
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
 
       PROVIDE global_selection_buffer
 
@@ -319,7 +319,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call set_multiple_levels_omp(.True.)
 !      call omp_set_max_active_levels(8)
 
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+
 
       do k=1,N_states
         pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
@@ -417,6 +418,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
   double precision :: rss
   double precision, external :: memory_of_double, memory_of_int
 
+  character(len=20) :: format_str1, str_error1, format_str2, str_error2
+  character(len=20) :: format_str3, str_error3, format_str4, str_error4
+  character(len=20) :: format_value1, format_value2, format_value3, format_value4
+  character(len=20) :: str_value1, str_value2, str_value3, str_value4
+  character(len=20) :: str_conv
+  double precision  :: value1, value2, value3, value4
+  double precision  :: error1, error2, error3, error4
+  integer           :: size1,size2,size3,size4
+
+  double precision :: conv_crit
+
   sending =.False.
 
   rss  = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2)
@@ -540,14 +552,60 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
 
           if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
             time1 = time
-            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, &
-              pt2_data     % pt2(pt2_stoch_istate) +E, &
-              pt2_data_err % pt2(pt2_stoch_istate), &
-              pt2_data     % variance(pt2_stoch_istate), &
-              pt2_data_err % variance(pt2_stoch_istate), &
-              pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
-              pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
-              time-time0
+
+            value1 =  pt2_data     % pt2(pt2_stoch_istate) + E
+            error1 =  pt2_data_err % pt2(pt2_stoch_istate)           
+            value2 =  pt2_data     % pt2(pt2_stoch_istate)
+            error2 =  pt2_data_err % pt2(pt2_stoch_istate)           
+            value3 =  pt2_data     % variance(pt2_stoch_istate)           
+            error3 =  pt2_data_err % variance(pt2_stoch_istate)           
+            value4 =  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            error4 =  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            
+            ! Max size of the values (FX.Y) with X=size
+            size1 = 15
+            size2 = 12
+            size3 = 12
+            size4 = 12
+
+            ! To generate the format: number(error)
+            call format_w_error(value1,error1,size1,8,format_value1,str_error1)
+            call format_w_error(value2,error2,size2,8,format_value2,str_error2)
+            call format_w_error(value3,error3,size3,8,format_value3,str_error3)
+            call format_w_error(value4,error4,size4,8,format_value4,str_error4)
+
+            ! value > string with the right format
+            write(str_value1,'('//format_value1//')') value1
+            write(str_value2,'('//format_value2//')') value2
+            write(str_value3,'('//format_value3//')') value3
+            write(str_value4,'('//format_value4//')') value4
+
+            ! Convergence criterion
+            conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+                  (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+            write(str_conv,'(G10.3)') conv_crit
+
+            write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,&
+            adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),&
+            adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),&
+            adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),&
+            adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),&
+            adjustl(str_conv),& 
+            time-time0
+
+            ! Old print
+            !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, &
+            !  pt2_data     % pt2(pt2_stoch_istate) +E, &
+            !  pt2_data_err % pt2(pt2_stoch_istate), &
+            !  pt2_data     % variance(pt2_stoch_istate), &
+            !  pt2_data_err % variance(pt2_stoch_istate), &
+            !  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  time-time0, &
+            !  pt2_data     % pt2(pt2_stoch_istate), &
+            !  dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+            !      (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+
             if (stop_now .or. (                                      &
                   (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
                   (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f
index 50104d69..ddfc050e 100644
--- a/src/cipsi/slave_cipsi.irp.f
+++ b/src/cipsi/slave_cipsi.irp.f
@@ -5,7 +5,6 @@ subroutine run_slave_cipsi
   END_DOC
 
   call set_multiple_levels_omp(.False.)
-!  call omp_set_max_active_levels(1)
   distributed_davidson = .False.
   read_wf = .False.
   SOFT_TOUCH read_wf distributed_davidson
@@ -173,10 +172,8 @@ subroutine run_slave_main
 
       !---
       call set_multiple_levels_omp(.True.)
-!      call omp_set_max_active_levels(8)
       call davidson_slave_tcp(0)
       call set_multiple_levels_omp(.False.)
-!      call omp_set_max_active_levels(1)
       print *,  mpi_rank, ': Davidson done'
       !---
 
diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f
index f72197c2..2b16a5f7 100644
--- a/src/cis/cis.irp.f
+++ b/src/cis/cis.irp.f
@@ -74,7 +74,7 @@ subroutine run
     print*,'******************************************************'
     print*,'Excitation energies (au)                     (eV)'
     do i = 2, N_states
-      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0
+      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev
     enddo
     print*,''
   endif
diff --git a/src/csf/configurations.irp.f b/src/csf/configurations.irp.f
index c11a49a4..aebf53d9 100644
--- a/src/csf/configurations.irp.f
+++ b/src/csf/configurations.irp.f
@@ -856,6 +856,7 @@ end subroutine
 !end subroutine
 !
  BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
+&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
 &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]
 
  implicit none
@@ -944,6 +945,29 @@ end subroutine
  enddo
 
  deallocate(dets, old_order)
+ integer :: ndet_conf
+ do i = 1, N_configuration
+  ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1
+  psi_configuration_n_det(i) = ndet_conf
+ enddo
 
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)]
+ implicit none
+ integer :: i,j,k,l
+ integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int)
+ n_elec_alpha_for_psi_configuration = 0
+ do i = 1, N_configuration
+  j = psi_configuration_to_psi_det(2,i) 
+  det_tmp(:,:) = psi_det(:,:,j)
+  k = 0
+  do l = 1, N_int
+   det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i))
+   k += popcnt(det_alpha(l))
+  enddo
+  n_elec_alpha_for_psi_configuration(i) = k
+ enddo
+
+END_PROVIDER 
diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f
index bdd8b327..494c3bfa 100644
--- a/src/csf/conversion.irp.f
+++ b/src/csf/conversion.irp.f
@@ -10,7 +10,6 @@ BEGIN_PROVIDER [ double precision, psi_csf_coef, (N_csf, N_states) ]
  call convertWFfromDETtoCSF(N_states, buffer, psi_csf_coef)
 END_PROVIDER
 
-
 subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
   use cfunctions
   use bitmasks
diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f
index 4bf154a9..4f88506a 100644
--- a/src/davidson/davidson_parallel.irp.f
+++ b/src/davidson/davidson_parallel.irp.f
@@ -509,7 +509,6 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
 
 
   call set_multiple_levels_omp(.True.)
-!  call omp_set_max_active_levels(5)
 
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
diff --git a/src/davidson/davidson_parallel_csf.irp.f b/src/davidson/davidson_parallel_csf.irp.f
index 4af779ab..d8e9bffa 100644
--- a/src/davidson/davidson_parallel_csf.irp.f
+++ b/src/davidson/davidson_parallel_csf.irp.f
@@ -465,7 +465,7 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
   endif
 
   call set_multiple_levels_omp(.True.)
-!  call omp_set_max_active_levels(4)
+
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
   if (ithread == 0 ) then
diff --git a/src/davidson/davidson_parallel_nos2.irp.f b/src/davidson/davidson_parallel_nos2.irp.f
index d408971c..597b001f 100644
--- a/src/davidson/davidson_parallel_nos2.irp.f
+++ b/src/davidson/davidson_parallel_nos2.irp.f
@@ -465,7 +465,7 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze)
   endif
 
   call set_multiple_levels_omp(.True.)
-!  call omp_set_max_active_levels(4)
+
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
   if (ithread == 0 ) then
diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f
index 3020ecd8..f531fb3a 100644
--- a/src/davidson/diagonalization_hcsf_dressed.irp.f
+++ b/src/davidson/diagonalization_hcsf_dressed.irp.f
@@ -300,7 +300,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
       shift  = N_st_diag*(iter-1)
       shift2 = N_st_diag*iter
 
-      if ((iter > 1).or.(itertot == 1)) then
+!      if ((iter > 1).or.(itertot == 1)) then
         ! Compute |W_k> = \sum_i |i><i|H|u_k>
         ! -----------------------------------
 
@@ -310,10 +310,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
         else
             call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
         endif
-      else
-         ! Already computed in update below
-         continue
-      endif
+!      else
+!         ! Already computed in update below
+!         continue
+!      endif
 
       if (dressing_state > 0) then
 
diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f
index 64a68482..f4c05595 100644
--- a/src/davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/davidson/diagonalization_hs2_dressed.irp.f
@@ -77,9 +77,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   !$OMP END PARALLEL
 
   if (dressing_state > 0) then
-    do k = 1, N_st
+    do k=1,N_st
 
-      do i = 1, sze
+      do i=1,sze
         H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
       enddo
 
diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f
index 7c4a7fec..fb461848 100644
--- a/src/determinants/density_matrix.irp.f
+++ b/src/determinants/density_matrix.irp.f
@@ -304,7 +304,7 @@ BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ]
        c0_weight(i) = c0_weight(i) * c
      enddo
    else
-     c0_weight = 1.d0
+     c0_weight(:) = 1.d0
    endif
 
 END_PROVIDER
@@ -321,7 +321,7 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
    if (weight_one_e_dm == 0) then
      state_average_weight(:) = c0_weight(:)
    else if (weight_one_e_dm == 1) then
-     state_average_weight(:) = 1./N_states
+     state_average_weight(:) = 1.d0/N_states
    else
      call ezfio_has_determinants_state_average_weight(exists)
      if (exists) then
@@ -384,6 +384,14 @@ END_PROVIDER
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)]
+ implicit none
+   BEGIN_DOC
+   !  one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta 
+   END_DOC
+    one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta 
+END_PROVIDER 
+
 
 subroutine get_occupation_from_dets(istate,occupation)
   implicit none
diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f
index c1774a26..a34608f9 100644
--- a/src/determinants/slater_rules.irp.f
+++ b/src/determinants/slater_rules.irp.f
@@ -652,7 +652,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
     case (1)
       call get_single_excitation(key_i,key_j,exc,phase,Nint)
       !DIR$ FORCEINLINE
-      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
       if (exc(0,1,1) == 1) then
         ! Single alpha
         m = exc(1,1,1)
@@ -671,10 +670,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
   end select
 end
 
-
-
-
-
 subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase)
   use bitmasks
   implicit none
@@ -1009,7 +1004,6 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
 
 end
 
-
 subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
   use bitmasks
   implicit none
diff --git a/src/determinants/slater_rules_wee_mono.irp.f b/src/determinants/slater_rules_wee_mono.irp.f
index 4c1c9330..7c2ad148 100644
--- a/src/determinants/slater_rules_wee_mono.irp.f
+++ b/src/determinants/slater_rules_wee_mono.irp.f
@@ -282,9 +282,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
   double precision               :: get_two_e_integral
   integer                        :: m,n,p,q
   integer                        :: i,j,k
-  integer                        :: occ(Nint*bit_kind_size,2)
   double precision               :: diag_H_mat_elem, phase,phase_2
-  integer                        :: n_occ_ab(2)
   PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy
 
   ASSERT (Nint > 0)
@@ -342,7 +340,6 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij)
     case (1)
       call get_single_excitation(key_i,key_j,exc,phase,Nint)
       !DIR$ FORCEINLINE
-      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
       if (exc(0,1,1) == 1) then
         ! Mono alpha
         m = exc(1,1,1)
diff --git a/src/dressing/run_dress_slave.irp.f b/src/dressing/run_dress_slave.irp.f
index 16feb798..08b654c9 100644
--- a/src/dressing/run_dress_slave.irp.f
+++ b/src/dressing/run_dress_slave.irp.f
@@ -73,7 +73,6 @@ subroutine run_dress_slave(thread,iproce,energy)
   ending = dress_N_cp+1
   ntask_tbd = 0
   call set_multiple_levels_omp(.True.)
-!  call omp_set_max_active_levels(8)
 
   !$OMP PARALLEL DEFAULT(SHARED) &
   !$OMP PRIVATE(interesting, breve_delta_m, task_id) &
@@ -86,7 +85,6 @@ subroutine run_dress_slave(thread,iproce,energy)
   integer, external :: connect_to_taskserver
   !$OMP CRITICAL
   call set_multiple_levels_omp(.False.)
-!  call omp_set_max_active_levels(1)
   if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
     print *,  irp_here, ': Unable to connect to task server'
     stop -1
@@ -299,7 +297,6 @@ subroutine run_dress_slave(thread,iproce,energy)
 
   !$OMP END PARALLEL
   call set_multiple_levels_omp(.False.)
-!  call omp_set_max_active_levels(1)
 !  do i=0,dress_N_cp+1
 !    call omp_destroy_lock(lck_sto(i))
 !  end do
diff --git a/src/ezfio_files/trexio.irp.f b/src/ezfio_files/trexio.irp.f
new file mode 100644
index 00000000..07eb6c0b
--- /dev/null
+++ b/src/ezfio_files/trexio.irp.f
@@ -0,0 +1,51 @@
+use trexio
+
+BEGIN_PROVIDER [ logical, use_trexio ]
+ implicit none
+ BEGIN_DOC
+ ! Returns the content of the QP_USE_TREXIO variable
+ END_DOC
+ character(32) :: buffer
+
+ call getenv('QP_USE_TREXIO', buffer)
+ if (trim(buffer) == '0') then
+   print *, 'Using EZFIO storage'
+   use_trexio = .False.
+ else
+   print *, 'Using TREXIO storage'
+   use_trexio = .True.
+ endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ character*(1024), trexio_filename ]
+  implicit none
+  BEGIN_DOC
+  ! Name of the TREXIO file.
+  END_DOC
+  trexio_filename = trim(ezfio_work_dir)//'/trexio.h5'
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer(trexio_back_end_t), trexio_backend ]
+  implicit none
+  BEGIN_DOC
+  ! Name of the TREXIO file.
+  END_DOC
+  trexio_backend = TREXIO_HDF5
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer(trexio_t), trexio_file ]
+  implicit none
+  BEGIN_DOC
+  ! Name of the TREXIO file.
+  END_DOC
+  integer (trexio_exit_code) :: rc
+
+  trexio_file = 0_trexio_t
+  if (use_trexio) then
+    trexio_file = trexio_open(trexio_filename, 'u', trexio_backend, rc)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+  endif
+
+END_PROVIDER
+
diff --git a/src/iterations/print_extrapolation.irp.f b/src/iterations/print_extrapolation.irp.f
index cb46fb67..7c6dbb9b 100644
--- a/src/iterations/print_extrapolation.irp.f
+++ b/src/iterations/print_extrapolation.irp.f
@@ -35,12 +35,13 @@ subroutine print_extrapolated_energy
     do k=2,min(N_iter,8)
       write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
           extrapolated_energy(k,i) - extrapolated_energy(k,1), &
-          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
+          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * ha_to_ev
     enddo
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
   enddo
 
   print *,  ''
+  call ezfio_set_fci_energy_extrapolated(extrapolated_energy(min(N_iter,3),1:N_states))
 
 end subroutine
 
diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f
index 8e6285e2..a0db3534 100644
--- a/src/iterations/print_summary.irp.f
+++ b/src/iterations/print_summary.irp.f
@@ -36,7 +36,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
   write(*,fmt) '# E          ', e_(1:N_states_p)
   if (N_states_p > 1) then
     write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
-    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
+    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*ha_to_ev
   endif
   write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
   write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
@@ -47,8 +47,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
   if (N_states_p > 1) then
     write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
       dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
-    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
+    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*ha_to_ev, &
+      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*ha_to_ev, k=1,N_states_p)
   endif
   write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
   write(*,fmt)
@@ -82,19 +82,19 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
     print *, 'Variational Energy difference (au | eV)'
     do i=2, N_states_p
       print*,'Delta E = ', (e_(i) - e_(1)), &
-        (e_(i) - e_(1)) * 27.211396641308d0
+        (e_(i) - e_(1)) * ha_to_ev
     enddo
     print *,  '-----'
     print*, 'Variational + perturbative Energy difference (au | eV)'
     do i=2, N_states_p
       print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
-        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
+        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * ha_to_ev
     enddo
     print *,  '-----'
     print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
     do i=2, N_states_p
       print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
-        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
+        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * ha_to_ev
     enddo
   endif
 
diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f
index d58932ce..6f4c5c17 100644
--- a/src/mo_two_e_ints/mo_bi_integrals.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f
@@ -53,7 +53,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
 !    call four_idx_novvvv
     call four_idx_novvvv_old
   else
-    call add_integrals_to_map(full_ijkl_bitmask_4)
+    if (32.d-9*dble(ao_num)**4 < dble(qp_max_mem)) then
+      call four_idx_dgemm
+    else
+      call add_integrals_to_map(full_ijkl_bitmask_4)
+    endif
   endif
 
   call wall_time(wall_2)
@@ -77,6 +81,94 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
 
 END_PROVIDER
 
+subroutine four_idx_dgemm
+  implicit none
+  integer :: p,q,r,s,i,j,k,l
+  double precision, allocatable :: a1(:,:,:,:)
+  double precision, allocatable :: a2(:,:,:,:)
+
+  allocate (a1(ao_num,ao_num,ao_num,ao_num))
+
+  print *, 'Getting AOs'
+  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s)
+  do s=1,ao_num
+    do r=1,ao_num
+      do q=1,ao_num
+        call get_ao_two_e_integrals(q,r,s,ao_num,a1(1,q,r,s))
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+  print *, '1st transformation'
+  ! 1st transformation
+  allocate (a2(ao_num,ao_num,ao_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num))
+
+  ! 2nd transformation
+  print *, '2nd transformation'
+  deallocate (a1)
+  allocate (a1(ao_num,ao_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num))
+
+  ! 3rd transformation
+  print *, '3rd transformation'
+  deallocate (a2)
+  allocate (a2(ao_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num))
+
+  ! 4th transformation
+  print *, '4th transformation'
+  deallocate (a1)
+  allocate (a1(mo_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num))
+
+  deallocate (a2)
+
+  integer :: n_integrals, size_buffer
+  integer(key_kind)  , allocatable :: buffer_i(:)
+  real(integral_kind), allocatable :: buffer_value(:)
+  size_buffer = min(ao_num*ao_num*ao_num,16000000)
+
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
+  allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
+
+  n_integrals = 0
+  !$OMP DO
+  do l=1,mo_num
+    do k=1,mo_num
+      do j=1,l
+        do i=1,k
+            if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then
+              cycle
+            endif
+            n_integrals += 1
+            buffer_value(n_integrals) = a1(i,j,k,l)
+            !DIR$ FORCEINLINE
+            call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals))
+            if (n_integrals == size_buffer) then
+              call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+              n_integrals = 0
+            endif
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+
+  call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+
+  deallocate(buffer_i, buffer_value)
+  !$OMP END PARALLEL
+
+  deallocate (a1)
+
+  call map_unique(mo_integrals_map)
+
+  integer*8                      :: get_mo_map_size, mo_map_size
+  mo_map_size = get_mo_map_size()
+
+end subroutine
 
 subroutine add_integrals_to_map(mask_ijkl)
   use bitmasks
diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f
index a567b9c7..5188581a 100644
--- a/src/scf_utils/diagonalize_fock.irp.f
+++ b/src/scf_utils/diagonalize_fock.irp.f
@@ -73,11 +73,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
    liwork = -1
 
    F_save = F
-   !print *, ' Fock matrix'
-   !do i = 1, mo_num
-   !  write(*, '(1000(F16.10,X))') F_save(:,i)
-   !enddo
-
    call dsyevd( 'V', 'U', mo_num, F,                             &
        size(F,1), diag, work, lwork, iwork, liwork, info)
 
@@ -108,16 +103,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
      endif
    endif
 
-  !print *, ' eigenvalues'
-  !do i = 1, mo_num
-  !  write(*, '(1000(F16.10,X))') diag(i)
-  !enddo
-  !print *, ' eigenvectors'
-  !do i = 1, mo_num
-  !  write(*, '(1000(F16.10,X))') F(:,i)
-  !enddo
-
-
    call dgemm('N','N',ao_num,mo_num,mo_num, 1.d0,            &
        mo_coef, size(mo_coef,1), F, size(F,1),                       &
        0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
diff --git a/src/scf_utils/scf_density_matrix_ao.irp.f b/src/scf_utils/scf_density_matrix_ao.irp.f
index 34ac2bec..7759431d 100644
--- a/src/scf_utils/scf_density_matrix_ao.irp.f
+++ b/src/scf_utils/scf_density_matrix_ao.irp.f
@@ -12,16 +12,6 @@ BEGIN_PROVIDER [double precision, SCF_density_matrix_ao_alpha, (ao_num,ao_num) ]
      SCF_density_matrix_ao_alpha = 0.d0
    endif
 
-!  integer          :: i, j
-!  double precision :: trace_density
-!  trace_density = 0.d0
-!  do i = 1, ao_num !elec_alpha_num
-!    do j = 1, ao_num !elec_alpha_num
-!      trace_density = trace_density &
-!                    + SCF_density_matrix_ao_alpha(j,i) * ao_overlap(j,i)
-!    enddo
-!  enddo
-!  print *, ' trace of SCF_density_matrix_ao_alpha =', trace_density
 
 END_PROVIDER
 
diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f
index 417b25ad..830a141e 100644
--- a/src/tools/molden.irp.f
+++ b/src/tools/molden.irp.f
@@ -52,8 +52,8 @@ program molden
         l += 1
         if (l > ao_num) exit
       enddo
-      write(i_unit_output,*)''
     enddo
+    write(i_unit_output,*)''
   enddo
 
 
diff --git a/src/tools/print_cfg.irp.f b/src/tools/print_cfg.irp.f
new file mode 100644
index 00000000..9e9be7cf
--- /dev/null
+++ b/src/tools/print_cfg.irp.f
@@ -0,0 +1,27 @@
+program print_energy
+ implicit none
+ BEGIN_DOC
+ ! Prints the energy of the wave function stored in the |EZFIO| directory.
+ END_DOC
+
+ ! this has to be done in order to be sure that N_det, psi_det and
+ ! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
+ read_wf = .True.
+ touch read_wf
+ PROVIDE N_states
+ call run
+end
+
+subroutine run
+ implicit none
+ integer :: i,j
+
+ do i=1,N_configuration
+   print *, i, sum(popcnt(psi_configuration(:,1,i)))
+ enddo
+
+ print *, ''
+ do i=0,elec_num
+   print *, i, cfg_seniority_index(i)
+ enddo
+end
diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f
index 3b9f6978..46cae2dc 100644
--- a/src/tools/print_wf.irp.f
+++ b/src/tools/print_wf.irp.f
@@ -33,8 +33,9 @@ subroutine routine
  double precision :: norm_mono_a,norm_mono_b
  double precision :: norm_mono_a_2,norm_mono_b_2
  double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2
- double precision :: norm_mono_a_pert,norm_mono_b_pert
+ double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1
  double precision :: delta_e,coef_2_2
+
  norm_mono_a = 0.d0
  norm_mono_b = 0.d0
  norm_mono_a_2 = 0.d0
@@ -43,6 +44,7 @@ subroutine routine
  norm_mono_b_pert = 0.d0
  norm_mono_a_pert_2 = 0.d0
  norm_mono_b_pert_2 = 0.d0
+ norm_double_1 = 0.d0
  do i = 1, min(N_det_print_wf,N_det)
   print*,''
   print*,'i = ',i
@@ -94,6 +96,7 @@ subroutine routine
     print*,'h1,p1 = ',h1,p1
     print*,'s2',s2
     print*,'h2,p2 = ',h2,p2
+    norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))
    endif
 
    print*,'<Ref| H |D_I> = ',hij
@@ -110,6 +113,7 @@ subroutine routine
  print*,''
  print*,'L1 norm of mono alpha = ',norm_mono_a
  print*,'L1 norm of mono beta  = ',norm_mono_b
+ print*,'L1 norm of double exc = ',norm_double_1
  print*, '---'
  print*,'L2 norm of mono alpha = ',norm_mono_a_2
  print*,'L2 norm of mono beta  = ',norm_mono_b_2
diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f
index 6c66c8ec..8ed97391 100644
--- a/src/tools/truncate_wf.irp.f
+++ b/src/tools/truncate_wf.irp.f
@@ -55,10 +55,23 @@ subroutine routine_s2
  integer :: i,j,k
  double precision :: accu(N_states)
 
- print *,  'Weights of the CFG'
+ integer :: weights(0:16), ix
+ double precision :: x
+
+ weights(:) = 0
  do i=1,N_det
-   print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
+   x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0)
+   ix = min(int(x), 16)
+   weights(ix) += 1
  enddo
+
+ print *,  'Histogram of the weights of the CFG'
+ do i=0,15
+   print *, '  10^{-', i, '}   ', weights(i) 
+ end do
+ print *, '< 10^{-', 15, '}   ', weights(16) 
+
+
  print*, 'Min weight of the configuration?'
  read(5,*) wmin
 
diff --git a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f
index eb247dea..26ed5ae6 100644
--- a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f
+++ b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f
@@ -529,10 +529,14 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_
        c_average += c_1(l) * c_1(l) * state_weights(l) 
      enddo
      
-     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     if (nkeys > 0) then
+       call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     endif
      nkeys = 0
      call orb_range_diag_to_all_2_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     if (nkeys > 0) then
+       call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+     endif
      nkeys = 0
      
    end do
diff --git a/src/utils/format_w_error.irp.f b/src/utils/format_w_error.irp.f
new file mode 100644
index 00000000..1378d367
--- /dev/null
+++ b/src/utils/format_w_error.irp.f
@@ -0,0 +1,71 @@
+subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error)
+
+  implicit none
+ 
+  BEGIN_DOC
+  ! Format for double precision, value(error)
+  END_DOC
+
+  ! in
+  ! | value         | double precision | value...      |
+  ! | error         | double precision | error...      |
+  ! | size_nb       | integer          | X in FX.Y     |
+  ! | max_nb_digits | integer          | Max Y in FX.Y |
+
+  ! out
+  ! | format_value  | character | string FX.Y for the format |
+  ! | str_error     | character | string of the error |
+
+  ! internal
+  ! | str_size      | character | size in string format                             |
+  ! | nb_digits     | integer   | number of digits Y in FX.Y depending of the error |
+  ! | str_nb_digits | character | nb_digits in string format                        |
+  ! | str_exp       | character | string of the value in exponential format         |
+
+  ! in
+  double precision, intent(in)   :: error, value
+  integer, intent(in)            :: size_nb, max_nb_digits
+
+  ! out
+  character(len=20), intent(out) :: str_error, format_value
+
+  ! internal
+  character(len=20)              :: str_size, str_nb_digits, str_exp
+  integer                        :: nb_digits
+
+  ! max_nb_digit: Y max
+  ! size_nb = Size of the double: X (FX.Y)
+  write(str_size,'(I3)') size_nb
+
+  ! Error
+  write(str_exp,'(1pE20.0)') error
+  str_error = trim(adjustl(str_exp))
+  
+  ! Number of digit: Y (FX.Y) from the exponent
+  str_nb_digits = str_exp(19:20)
+  read(str_nb_digits,*) nb_digits
+ 
+  ! If the error is 0d0
+  if (error <= 1d-16) then 
+    write(str_nb_digits,*) max_nb_digits
+  endif
+
+  ! If the error is too small 
+  if (nb_digits > max_nb_digits) then
+      write(str_nb_digits,*) max_nb_digits
+      str_error(1:1) = '0'
+  endif
+
+  ! If the error is too big (>= 0.5)
+  if (error >= 0.5d0) then
+    str_nb_digits = '1'
+    str_error(1:1) = '*'
+  endif
+
+  ! FX.Y,A1,A1,A1 for value(str_error)
+  !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1'
+
+  ! FX.Y just for the value 
+  format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
+
+end
diff --git a/src/utils/units.irp.f b/src/utils/units.irp.f
new file mode 100644
index 00000000..1850b28b
--- /dev/null
+++ b/src/utils/units.irp.f
@@ -0,0 +1,22 @@
+BEGIN_PROVIDER [double precision, ha_to_ev]
+
+  implicit none
+  BEGIN_DOC
+  ! Converstion from Hartree to eV
+  END_DOC
+
+  ha_to_ev = 27.211396641308d0
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, au_to_D]
+
+  implicit none
+  BEGIN_DOC
+  ! Converstion from au to Debye
+  END_DOC
+
+  au_to_D = 2.5415802529d0
+
+END_PROVIDER
+
diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f
index 5279fedd..e7f00ce2 100644
--- a/src/utils/util.irp.f
+++ b/src/utils/util.irp.f
@@ -37,6 +37,10 @@ double precision function binom_func(i,j)
   else
     binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) )
   endif
+
+  ! To avoid .999999 numbers
+  binom_func = floor(binom_func + 0.5d0)
+
 end