program e_curve
  use bitmasks
  implicit none
  integer :: i,j,k, nab, m, l
  double precision :: norm 
  integer, allocatable :: iorder(:)
  double precision , allocatable :: norm_sort(:)
  PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals


  if (.not.read_wf) then
    stop 'Please set read_wf to true'
  endif

  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,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

  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)
      call dset_order(psi_coef(1,k),psi_bilinear_matrix_order_reverse,N_det)
    enddo
    TOUCH psi_det psi_coef

    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

    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
  enddo
  print *,  '=========================================================='

  deallocate (iorder, norm_sort)
end