diff --git a/stable/utilities/truncate_nonzero.irp.f b/stable/utilities/truncate_nonzero.irp.f
new file mode 100644
index 0000000..0829d45
--- /dev/null
+++ b/stable/utilities/truncate_nonzero.irp.f
@@ -0,0 +1,54 @@
+program truncate_wf
+ implicit none
+ BEGIN_DOC
+! Truncate the wave function
+ END_DOC
+ read_wf = .True.
+ call routine
+end
+
+subroutine routine
+ implicit none
+ integer :: ndet_max
+ integer(bit_kind), allocatable :: psi_det_tmp(:,:,:)
+ double precision, allocatable :: psi_coef_tmp(:,:)
+
+ integer :: i,j
+ double precision :: accu(N_states)
+
+ ndet_max=N_det
+
+ do i = 1, n_det
+   if (psi_average_norm_contrib_sorted(i) < 1.e-24) then
+     ndet_max = i-1
+     exit
+   endif
+ enddo
+ print *,  "N_det = ", N_det, " -> ", ndet_max
+
+ allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states))
+ 
+ accu = 0.d0
+ do i = 1, ndet_max
+  do j = 1, N_int
+   psi_det_tmp(j,1,i) = psi_det_sorted(j,1,i)
+   psi_det_tmp(j,2,i) = psi_det_sorted(j,2,i)
+  enddo
+  do j = 1, N_states
+   psi_coef_tmp(i,j) = psi_coef_sorted(i,j)
+   accu(j) += psi_coef_tmp(i,j) **2
+  enddo
+ enddo
+ do j = 1, N_states
+  accu(j) = 1.d0/dsqrt(accu(j))
+ enddo
+ do j = 1, N_states
+  do i = 1, ndet_max
+   psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j)
+  enddo
+ enddo
+
+ call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,size(psi_coef_tmp,1),psi_coef_tmp)
+ 
+end
+