program mpn
  implicit none
  BEGIN_DOC
! TODO : Put the documentation of the program here
  END_DOC
  integer :: i, k, l
  double precision, allocatable :: c_pert(:,:)
  double precision, allocatable :: e_pert(:)
  double precision, allocatable :: hc(:), s2(:)

  n_states_diag = 1
  TOUCH n_states_diag
  call generate_fci_space
  allocate(c_pert(N_det,0:mp_order))
  allocate(s2(N_det))
  allocate(e_pert(mp_order+1))
  e_pert = 0.d0
  c_pert(:,:) = 0.d0
  c_pert(1,0) = 1.d0

double precision :: hij

  do k=1,mp_order
    ! H_ij C^(k-1)
    call h_s2_u_0_nstates_zmq(c_pert(1,k),s2,c_pert(1,k-1),1,N_det)
    e_pert(k) = c_pert(1,k)
    print *, k, e_pert(k), sum(e_pert) + nuclear_repulsion

    c_pert(1,k) = 0.d0
    c_pert(:,k) = -c_pert(:,k)
    do l=1,k-1
      do i=2,N_det
        c_pert(i,k) = c_pert(i,k) + e_pert(l) * c_pert(i,k-l)
      enddo
    enddo
    do i=2,N_det
      c_pert(i,k) = c_pert(i,k) + energy_det_i(i) * c_pert(i,k-1)
    enddo
    do i=2,N_det
      c_pert(i,k) = c_pert(i,k) / (energy_det_i(i) - energy_det_i(1))
    enddo

  enddo


end

subroutine generate_fci_space
  use bitmasks
  implicit none
  integer :: i, sze
  integer(bit_kind) :: o(N_int,2)

  if (mo_num > 64) then
     stop 'No more than 64 MOs'
  endif
  o(:,1) = full_ijkl_bitmask(:)
  o(:,2) = 0_bit_kind

  call configuration_to_dets_size(o,n_det_alpha_unique,elec_alpha_num,N_int)
  TOUCH n_det_alpha_unique

  integer :: k,n,m, t, t1, t2
  integer(bit_kind) :: u
  k=0
  n = elec_alpha_num
  m = mo_num - n
  u = shiftl(1_bit_kind,n) -1
  do while (u < shiftl(1_bit_kind,n+m))
    k = k+1
    psi_det_alpha_unique(1,k) = u
    t = ior(u,u-1)
    t1 = t+1
    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
    u = ior(t1,t2)
  enddo


  call configuration_to_dets_size(o,n_det_beta_unique,elec_beta_num,N_int)
  TOUCH n_det_beta_unique

  k=0
  n = elec_beta_num
  m = mo_num - n
  u = shiftl(1_bit_kind,n) -1
  do while (u < shiftl(1_bit_kind,n+m))
    k = k+1
    psi_det_beta_unique(1,k) = u
    t = ior(u,u-1)
    t1 = t+1
    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
    u = ior(t1,t2)
  enddo

  call generate_all_alpha_beta_det_products

  print *, N_det


end