program cisd
  implicit none
  BEGIN_DOC
  ! Multi-Reference Configuration Interaction with Single and Double excitations.
  !
  ! This program takes a reference wave function as input and
  ! and performs all single and double excitations on top of it, disregarding
  ! spatial symmetry and compute the "n_states" lowest eigenstates of that CI
  ! matrix (see :option:`determinants n_states`).
  !
  END_DOC
  PROVIDE N_states
  read_wf = .True.
  SOFT_TOUCH read_wf
  call run
end

subroutine run
  implicit none
  integer                        :: i,k
  double precision               :: cisdq(N_states), delta_e
  double precision,external      :: diag_h_mat_elem

  if(do_singles)then
    if(do_doubles)then
      call H_apply_cisd
    else
      call H_apply_cis
    endif
  else
    if(do_doubles)then
      call H_apply_cid
    else
      stop "do_singles and do_doubles are both False"
    endif
  endif

  psi_coef = ci_eigenvectors
  SOFT_TOUCH psi_coef
  call save_wavefunction
  call ezfio_set_general_mrci_energy(CI_energy)

  do i = 1,N_states
    k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
    delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
    cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
  enddo
  print *,  'N_det = ', N_det
  print*,''
  print*,'******************************'
  print *,  'MR-CI Energies'
  do i = 1,N_states
    print *,  i, CI_energy(i)
  enddo
  print*,''
  print*,'******************************'
  print *,  'MR-CI+Q Energies'
  do i = 1,N_states
    print *,  i, cisdq(i)
  enddo
  if (N_states > 1) then
    print*,''
    print*,'******************************'
    print*,'Excitation energies (au)    (CISD+Q)'
    do i = 2, N_states
      print*, i ,CI_energy(i) - CI_energy(1), cisdq(i) - cisdq(1)
    enddo
    print*,''
    print*,'******************************'
    print*,'Excitation energies (eV)    (CISD+Q)'
    do i = 2, N_states
      print*, i ,(CI_energy(i) - CI_energy(1))/0.0367502d0, &
        (cisdq(i) - cisdq(1)) / 0.0367502d0
    enddo
  endif

end