subroutine calc_determ_hamil_sparse(rep)
type(core_space_t), intent(inout) :: rep
integer :: i, j, row_size, counter, ierr
integer :: nI(nel), nJ(nel)
integer(n_int), allocatable, dimension(:, :) :: temp_store
integer(TagIntType) :: HRTag, TempStoreTag
HElement_t(dp), allocatable, dimension(:) :: hamiltonian_row
integer :: pos, nExcit
integer(n_int) :: ilutG(0:nifguga)
integer(n_int), pointer :: excitations(:, :)
type(ExcitationInformation_t) :: excitInfo
character(len=*), parameter :: this_routine = "calc_determ_hamil_sparse"
integer(n_int) :: tmp(0:NIfD), ilutI_tmp(0:NIfTot), ilutJ_tmp(0:niftot)
integer :: IC, nI_tmp(nel), nJ_tmp(nel)
integer(n_int) :: ilutI(0:niftot), ilutJ(0:niftot)
type(CSF_Info_t) :: csf_i, csf_j
HElement_t(dp) :: tmp_mat, tmp_mat_2, HOffDiag
allocate(rep%sparse_core_ham(rep%determ_sizes(iProcIndex)), stat=ierr)
allocate(rep%SparseCoreHamilTags(2, rep%determ_sizes(iProcIndex)))
allocate(hamiltonian_row(rep%determ_space_size), stat=ierr)
call LogMemAlloc('hamiltonian_row', int(rep%determ_space_size), 8, this_routine, HRTag, ierr)
allocate(rep%core_ham_diag(rep%determ_sizes(iProcIndex)), stat=ierr)
allocate(temp_store(0:NIfTot, rep%determ_space_size), stat=ierr)
call LogMemAlloc('temp_store', rep%determ_space_size * (NIfTot + 1), 8, this_routine, TempStoreTag, ierr)
! Stick together the deterministic states from all processors, on
! all processors.
! n.b. Explicitly use 0:NIfTot, as NIfTot may not equal NIfBCast
call MPIAllGatherV(SpawnedParts(0:NIfTot, 1:rep%determ_sizes(iProcIndex)), &
temp_store(0:niftot, 1:), rep%determ_sizes, rep%determ_displs)
! Loop over all deterministic states on this processor.
do i = 1, rep%determ_sizes(iProcIndex)
ilutI = SpawnedParts(0:niftot, i)
call decode_bit_det(nI, IlutI)
if (tGUGA) csf_i = CSF_Info_t(ilutI)
row_size = 0
hamiltonian_row = 0.0_dp
! Loop over all deterministic states.
do j = 1, rep%determ_space_size
ilutJ = temp_store(:, j)
call decode_bit_det(nJ, ilutJ)
if (tGUGA) csf_j = CSF_Info_t(ilutJ)
if(t_evolve_adjoint(rep%first_run())) then
nI_tmp = nJ
nJ_tmp = nI
ilutI_tmp = ilutJ
ilutJ_tmp = ilutI
else
nI_tmp = nI
nJ_tmp = nJ
ilutI_tmp = ilutI
ilutJ_tmp = ilutJ
end if
! If on the diagonal of the Hamiltonian.
if (DetBitEq(IlutI, ilutJ_tmp, nifd)) then
hamiltonian_row(j) = get_diagonal_matel(nI_tmp, IlutI_tmp) - Hii
rep%core_ham_diag(i) = hamiltonian_row(j)
! We calculate and store the diagonal matrix element at
! this point for later access.
if (.not. tReadPops) then
call set_det_diagH(i, Real(hamiltonian_row(j), dp))
HOffDiag = get_off_diagonal_matel(nI_tmp, ilutI_tmp)
call set_det_offdiagH(i, HOffDiag)
end if
! Always include the diagonal elements.
row_size = row_size + 1
else
if (tHPHF) then
hamiltonian_row(j) = hphf_off_diag_helement(nI_tmp, nJ_tmp, IlutI_tmp, ilutJ_tmp)
else if (tGUGA) then
call calc_guga_matrix_element(&
ilutI_tmp, csf_i, ilutJ_tmp, csf_j, excitInfo, tmp_mat, .true.)
#ifdef DEBUG_
call calc_guga_matrix_element(&
ilutI_tmp, csf_i, ilutJ_tmp, csf_j, excitInfo, tmp_mat_2, .true.)
if (.not. near_zero(tmp_mat - tmp_mat_2)) then
call stop_all(this_routine, "type 1 and 2 do not agree!")
end if
call calc_guga_matrix_element(&
ilutJ_tmp, csf_j, ilutI_tmp, csf_i, excitInfo, tmp_mat_2, .true.)
if (.not. near_zero(tmp_mat - conjgt(tmp_mat_2))) then
call stop_all(this_routine, "not hermititan!")
end if
#endif
#ifdef CMPLX_
hamiltonian_row(j) = conjg(tmp_mat)
#else
hamiltonian_row(j) = tmp_mat
#endif
else
tmp = ieor(IlutI_tmp(0:NIfD), ilutJ_tmp(0:NIfD))
tmp = iand(IlutI_tmp(0:NIfD), tmp)
IC = CountBits(tmp, NIfD)
if (IC <= maxExcit) then
! hamiltonian_row(j) = get_helement(nI_tmp, nJ_tmp, IC, ilutI_tmp, ilutJ_tmp)
hamiltonian_row(j) = get_helement(nJ_tmp, nI_tmp, IC, ilutJ_tmp, ilutI_tmp)
end if
end if
if (abs(hamiltonian_row(j)) > 0.0_dp) row_size = row_size + 1
end if
end do
! Now we know the number of non-zero elements in this row of the Hamiltonian, so allocate it.
call allocate_sparse_ham_row(rep%sparse_core_ham, i, row_size, "sparse_core_ham", rep%SparseCoreHamilTags(:, i))
rep%sparse_core_ham(i)%elements = 0.0_dp
rep%sparse_core_ham(i)%positions = 0
rep%sparse_core_ham(i)%num_elements = row_size
counter = 1
do j = 1, rep%determ_space_size
! If non-zero or a diagonal element.
if (abs(hamiltonian_row(j)) > 0.0_dp .or. (j == i + rep%determ_displs(iProcIndex))) then
rep%sparse_core_ham(i)%positions(counter) = j
rep%sparse_core_ham(i)%elements(counter) = hamiltonian_row(j)
counter = counter + 1
end if
if (counter == row_size + 1) exit
end do
end do
! Don't time the mpi_barrier call, because if we did then we wouldn't
! be able separate out some of the core Hamiltonian creation time from
! the MPIBarrier calls in the main loop.
call MPIBarrier(ierr, tTimeIn=.false.)
deallocate(temp_store, stat=ierr)
call LogMemDealloc(this_routine, TempStoreTag, ierr)
deallocate(hamiltonian_row, stat=ierr)
call LogMemDealloc(this_routine, HRTag, ierr)
end subroutine calc_determ_hamil_sparse