subroutine calc_hamil_exact(nvecs, krylov_array, array_len, h_matrix, h_diag)
use DetBitOps, only: FindBitExcitLevel
use Determinants, only: get_helement
use FciMCData, only: Hii, exact_subspace_h_time
use global_det_data, only: det_diagH
use hphf_integrals, only: hphf_off_diag_helement
use SystemData, only: tHPHF, nel
use timing_neci, only: set_timer, halt_timer
integer, intent(in) :: nvecs
integer(n_int), intent(in) :: krylov_array(0:, :)
integer, intent(in) :: array_len
real(dp), intent(out) :: h_matrix(:, :)
real(dp), intent(in), optional :: h_diag(:)
integer :: i, j, idet, jdet, ic
integer(n_int) :: ilut_1(0:NIfTot), ilut_2(0:NIfTot)
integer(n_int) :: int_sign(lenof_all_signs)
integer :: nI(nel), nJ(nel)
real(dp) :: real_sign_1(lenof_all_signs), real_sign_2(lenof_all_signs)
real(dp) :: h_elem
logical :: any_occ, occ_1, occ_2
integer(4), allocatable :: occ_flags(:)
call set_timer(exact_subspace_h_time)
h_matrix = 0.0_dp
allocate(occ_flags(array_len))
occ_flags = 0
ilut_1 = 0_n_int
ilut_2 = 0_n_int
! Check to see if there are any replica 1 or 2 walkers on this determinant.
do idet = 1, array_len
int_sign = krylov_array(IlutBits%ind_pop:IlutBits%ind_pop + lenof_all_signs - 1, idet)
any_occ = .false.
do i = 1, nvecs
any_occ = any_occ .or. (int_sign(kp_ind_1(i)) /= 0)
end do
if (any_occ) occ_flags = ibset(occ_flags(idet), 0)
any_occ = .false.
do i = 1, nvecs
any_occ = any_occ .or. (int_sign(kp_ind_2(i)) /= 0)
end do
if (any_occ) occ_flags = ibset(occ_flags(idet), 1)
end do
! Loop over all determinants in krylov_array.
do idet = 1, array_len
ilut_1(0:nifd) = krylov_array(0:nifd, idet)
call decode_bit_det(nI, ilut_1)
int_sign = krylov_array(IlutBits%ind_pop:IlutBits%ind_pop + lenof_all_signs - 1, idet)
real_sign_1 = transfer(int_sign, real_sign_1)
occ_1 = btest(occ_flags(idet), 0)
occ_2 = btest(occ_flags(idet), 1)
do jdet = idet, array_len
! If one of the two determinants are unoccupied in all vectors,
! then cycle.
if (.not. ((occ_1 .and. btest(occ_flags(jdet), 1)) .or. &
(occ_2 .and. btest(occ_flags(jdet), 0)))) cycle
ilut_2(0:nifd) = krylov_array(0:nifd, jdet)
ic = FindBitExcitLevel(ilut_1, ilut_2)
if (ic > 2) cycle
call decode_bit_det(nJ, ilut_2)
int_sign = krylov_array(IlutBits%ind_pop:IlutBits%ind_pop + lenof_all_signs - 1, jdet)
real_sign_2 = transfer(int_sign, real_sign_1)
if (idet == jdet) then
if (present(h_diag)) then
h_elem = h_diag(idet) + Hii
else
h_elem = det_diagH(idet) + Hii
end if
else
if (tHPHF) then
h_elem = hphf_off_diag_helement(nI, nJ, ilut_1, ilut_2)
else
if (tGUGA) then
call stop_all("calc_hamil_exact", "modify for GUGA")
end if
h_elem = get_helement(nI, nJ, ic, ilut_1, ilut_2)
end if
end if
! Finally, add in the contribution to all of the Hamiltonian elements.
do i = 1, nvecs
do j = i, nvecs
if (idet == jdet) then
h_matrix(i, j) = h_matrix(i, j) + &
h_elem * (real_sign_1(kp_ind_1(i)) * real_sign_2(kp_ind_2(j)) + &
real_sign_1(kp_ind_2(i)) * real_sign_2(kp_ind_1(j))) / 2
else
h_matrix(i, j) = h_matrix(i, j) + &
h_elem * (real_sign_1(kp_ind_1(i)) * real_sign_2(kp_ind_2(j)) + &
real_sign_1(kp_ind_2(i)) * real_sign_2(kp_ind_1(j)) + &
real_sign_1(kp_ind_1(j)) * real_sign_2(kp_ind_2(i)) + &
real_sign_1(kp_ind_2(j)) * real_sign_2(kp_ind_1(i))) / 2
end if
end do
end do
end do
end do
do i = 1, nvecs
do j = 1, i - 1
h_matrix(i, j) = h_matrix(j, i)
end do
end do
deallocate(occ_flags)
call halt_timer(exact_subspace_h_time)
end subroutine calc_hamil_exact