diff --git a/src/tools/NEED b/src/tools/NEED
index 886ee1bb..08682750 100644
--- a/src/tools/NEED
+++ b/src/tools/NEED
@@ -1,2 +1,3 @@
 fci
 mo_two_e_erf_integrals
+aux_quantities
diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f
new file mode 100644
index 00000000..5c0743dd
--- /dev/null
+++ b/src/tools/print_wf.irp.f
@@ -0,0 +1,117 @@
+program print_wf
+ implicit none 
+ BEGIN_DOC
+ ! print the wave function stored in the EZFIO folder in the intermediate normalization 
+ !
+ ! it also prints a lot of information regarding the excitation operators from the reference determinant
+ ! 
+ ! and a first-order perturbative analysis of the wave function. 
+ ! 
+ ! If the wave function strongly deviates from the first-order analysis, something funny is going on :)
+ END_DOC
+
+
+ ! this has to be done in order to be sure that N_det, psi_det and psi_coef are the wave function stored in the EZFIO folder 
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none 
+ integer :: i
+ integer :: degree
+ double precision :: hij,hii,coef_1,h00
+ integer          :: exc(0:2,2,2)
+ double precision :: phase
+ integer :: h1,p1,h2,p2,s1,s2
+ double precision :: get_mo_bielec_integral
+ 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 :: delta_e,coef_2_2
+ norm_mono_a = 0.d0
+ norm_mono_b = 0.d0
+ norm_mono_a_2 = 0.d0
+ norm_mono_b_2 = 0.d0
+ norm_mono_a_pert = 0.d0
+ norm_mono_b_pert = 0.d0
+ norm_mono_a_pert_2 = 0.d0
+ norm_mono_b_pert_2 = 0.d0
+ do i = 1, min(10000,N_det)
+  print*,''
+  print*,'i = ',i
+  call debug_det(psi_det(1,1,i),N_int)
+  call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,1),degree,N_int)
+  print*,'degree = ',degree
+  if(degree == 0)then
+   print*,'Reference determinant '
+   call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,h00)
+  else 
+   call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii)
+   call i_H_j(psi_det(1,1,1),psi_det(1,1,i),N_int,hij)
+   delta_e = hii - h00
+   coef_1 = hij/(h00-hii)
+   if(hij.ne.0.d0)then
+    if (delta_e > 0.d0) then
+      coef_2_2 = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij * hij ))/ hij 
+    else
+      coef_2_2 = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * hij * hij )) /hij 
+    endif
+   endif
+   call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
+   call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+   print*,'phase = ',phase
+   if(degree == 1)then
+    print*,'s1',s1
+    print*,'h1,p1 = ',h1,p1
+    if(s1 == 1)then
+     norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
+     norm_mono_a_2 += dabs(psi_coef(i,1)/psi_coef(1,1))**2
+     norm_mono_a_pert += dabs(coef_1)
+     norm_mono_a_pert_2 += dabs(coef_1)**2
+    else
+     norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
+     norm_mono_b_2 += dabs(psi_coef(i,1)/psi_coef(1,1))**2
+     norm_mono_b_pert += dabs(coef_1)
+     norm_mono_b_pert_2 += dabs(coef_1)**2
+    endif
+    double precision :: hmono,hdouble
+    call  i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble,phase)
+    print*,'hmono         = ',hmono
+    print*,'hdouble       = ',hdouble
+    print*,'hmono+hdouble = ',hmono+hdouble
+    print*,'hij           = ',hij
+   else
+    print*,'s1',s1
+    print*,'h1,p1 = ',h1,p1
+    print*,'s2',s2
+    print*,'h2,p2 = ',h2,p2
+   endif
+   
+   print*,'<Ref| H |D_I> = ',hij
+   print*,'Delta E       = ',h00-hii
+   print*,'coef pert (1) = ',coef_1
+   print*,'coef 2x2      = ',coef_2_2
+   print*,'Delta E_corr  = ',psi_coef(i,1)/psi_coef(1,1) * hij
+  endif
+   print*,'amplitude     = ',psi_coef(i,1)/psi_coef(1,1)
+
+ enddo
+
+
+ print*,''
+ print*,'L1 norm of mono alpha = ',norm_mono_a
+ print*,'L1 norm of mono beta  = ',norm_mono_b
+ print*, '---'
+ print*,'L2 norm of mono alpha = ',norm_mono_a_2
+ print*,'L2 norm of mono beta  = ',norm_mono_b_2
+ print*, '-- perturbative mono'
+ print*,''
+ print*,'L1 norm of pert alpha = ',norm_mono_a_pert
+ print*,'L1 norm of pert beta  = ',norm_mono_b_pert
+ print*,'L2 norm of pert alpha = ',norm_mono_a_pert_2
+ print*,'L2 norm of pert beta  = ',norm_mono_b_pert_2
+
+end
diff --git a/src/aux_quantities/save_one_body_dm.irp.f b/src/tools/save_one_body_dm.irp.f
similarity index 100%
rename from src/aux_quantities/save_one_body_dm.irp.f
rename to src/tools/save_one_body_dm.irp.f