diff --git a/devel/qmcchem/qmc_e_curve.irp.f b/devel/qmcchem/qmc_e_curve.irp.f
index f9bc210..7b26a9f 100644
--- a/devel/qmcchem/qmc_e_curve.irp.f
+++ b/devel/qmcchem/qmc_e_curve.irp.f
@@ -1,33 +1,12 @@
 program e_curve
   use bitmasks
   implicit none
-  integer :: i,j,k, kk, nab, m, l
-  double precision :: norm, E, hij, num, ci, cj
+  integer :: i,j,k, nab, m, l
+  double precision :: norm 
   integer, allocatable :: iorder(:)
   double precision , allocatable :: norm_sort(:)
-  double precision :: e_0(N_states)
   PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
 
-  nab = n_det_alpha_unique+n_det_beta_unique
-  allocate ( norm_sort(0:nab), iorder(0:nab) )
-
-  double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
-  double precision, allocatable  :: u_0(:,:), v_0(:,:)
-
-
-  norm_sort(0) = 0.d0
-  iorder(0) = 0
-  do i=1,n_det_alpha_unique
-   norm_sort(i) = det_alpha_norm(i)
-   iorder(i) = i
-  enddo
-
-  do i=1,n_det_beta_unique
-   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
-   iorder(i+n_det_alpha_unique) = -i
-  enddo
-
-  call dsort(norm_sort(1),iorder(1),nab)
 
   if (.not.read_wf) then
     stop 'Please set read_wf to true'
@@ -35,40 +14,67 @@ program e_curve
 
   PROVIDE psi_bilinear_matrix_values nuclear_repulsion
 
+  if (.True.) then
+    print *,  ''
+    print *,  'Energy: ', psi_energy(1)+nuclear_repulsion
+  endif
   print *,  ''
   print *,  '=============================='
   print *,  'Energies at different cut-offs'
   print *,  '=============================='
   print *,  ''
-  print *,  '=========================================================='
-  print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
-  print *,  '=========================================================='
+  print *,  '=============================================================================='
+
+  print '(A8,2X,A8,A8,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Na', 'Nb', 'Cost', 'Norm', 'E'
+  print *,  '=============================================================================='
   double precision :: thresh
   integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
   integer :: na, nb
-  thresh = 1.d-10
-  na = n_det_alpha_unique
-  nb = n_det_beta_unique
-  do j=0,nab
-    i = iorder(j)
-    if (i<0) then
-      nb -= 1
-      do k=1,n_det
-        if (psi_bilinear_matrix_columns(k) == -i) then
-          psi_bilinear_matrix_values(k,1) = 0.d0
-        endif
-      enddo
-    else
-      na -= 1
-      do k=1,n_det
-        if (psi_bilinear_matrix_rows(k) ==  i) then
-          psi_bilinear_matrix_values(k,1) = 0.d0
-        endif
-      enddo
-    endif
-    if (thresh > norm_sort(j)) then
-      cycle
-    endif
+
+  nab = n_det_alpha_unique+n_det_beta_unique
+  allocate ( norm_sort(0:nab), iorder(0:nab) )
+  thresh = 1.d-8
+
+  do while (N_det > 1)
+    nab = n_det_alpha_unique+n_det_beta_unique
+
+    norm_sort(0) = 0.d0
+    iorder(0) = 0
+    do i=1,n_det_alpha_unique
+     norm_sort(i) = det_alpha_norm(i)
+     iorder(i) = i
+    enddo
+
+    do i=1,n_det_beta_unique
+     norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+     iorder(i+n_det_alpha_unique) = -i
+    enddo
+
+    call dsort(norm_sort(1),iorder(1),nab)
+
+    na = n_det_alpha_unique
+    nb = n_det_beta_unique
+    do j=1,nab
+      i = iorder(j)
+      if ((i<0).and.(nb>1)) then
+        nb -= 1
+        do k=1,n_det
+          if (psi_bilinear_matrix_columns(k) == -i) then
+            psi_bilinear_matrix_values(k,1) = 0.d0
+          endif
+        enddo
+      else if ((i>0).and.(na>1)) then
+        na -= 1
+        do k=1,n_det
+          if (psi_bilinear_matrix_rows(k) ==  i) then
+            psi_bilinear_matrix_values(k,1) = 0.d0
+          endif
+        enddo
+      endif
+      if (thresh < norm_sort(j)) then
+        exit
+      endif
+    enddo
 
     do k=1,N_states
       psi_coef(1:N_det,k) = psi_bilinear_matrix_values(1:N_det,k)
@@ -76,30 +82,31 @@ program e_curve
     enddo
     TOUCH psi_det psi_coef
 
-    m = 0
-    do k=1,n_det
-     if (psi_bilinear_matrix_values(k,1) /= 0.d0) then
-      m = m+1
-     endif
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+
+    do m=1,n_det
+     if (psi_coef_sorted(m,1) == 0.d0) exit
     enddo
 
-    if (m == 0) then
-       exit
+    N_det = m-1
+    TOUCH psi_det psi_coef N_det
+
+    ! Avoid providing psi_energy 
+    if (.True.) then
+      double precision :: cost0, cost
+      cost0 = elec_alpha_num**3 + elec_beta_num**3
+      cost = (na-1) * elec_alpha_num**2 + &
+             (nb-1) * elec_beta_num**2 + &
+              elec_alpha_num**3 + elec_beta_num**3
+      cost = cost/cost0
+
+      double precision :: u_dot_u
+      norm = dsqrt(u_dot_u(psi_coef(1,1),N_det))
+      print '(E9.1,2X,I8,I8,I8,2X,F10.2,2X,F10.8,2X,F15.10)',  thresh, N_det, &
+        na, nb, cost,  norm, psi_energy(1) + nuclear_repulsion
+      thresh = thresh * dsqrt(10.d0)
     endif
-    E = E_0(1) + nuclear_repulsion
-
-    double precision :: cost0, cost
-    cost0 = elec_alpha_num**3 + elec_beta_num**3
-    cost = (na-1) * elec_alpha_num**2 + &
-           (nb-1) * elec_beta_num**2 + &
-            elec_alpha_num**3 + elec_beta_num**3
-    cost = cost/cost0
-
-    double precision :: u_dot_u
-    norm = dsqrt(u_dot_u(psi_coef(1,1),N_det))
-    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)',  thresh, m, &
-      cost,  norm, psi_energy(1)
-    thresh = thresh * dsqrt(10.d0)
   enddo
   print *,  '=========================================================='