subroutine calculate_sparse_hamiltonian(num_states, ilut_list)
integer, intent(in) :: num_states
integer(n_int), intent(in) :: ilut_list(0:NIfTot, num_states)
integer :: i, j, counter, ierr
integer :: nI(nel), nJ(nel)
HElement_t(dp), allocatable, dimension(:) :: hamiltonian_row
integer, allocatable, dimension(:) :: sparse_diag_positions, sparse_row_sizes, indices
integer(TagIntType) :: HRTag, SRTag, SDTag, ITag
character(len=*), parameter :: t_r = "calculate_sparse_hamiltonian"
integer :: pos, nexcits
integer(n_int) :: ilutG(0:nifguga)
integer(n_int), pointer :: excitations(:, :)
type(ExcitationInformation_t) :: excitInfo
type(CSF_Info_t) :: csf_i, csf_j
allocate(sparse_ham(num_states))
allocate(SparseHamilTags(2, num_states))
allocate(hamiltonian_row(num_states), stat=ierr)
call LogMemAlloc('hamiltonian_row', num_states, 8, t_r, HRTag, ierr)
allocate(hamil_diag(num_states), stat=ierr)
call LogMemAlloc('hamil_diag', num_states, 8, t_r, HDiagTag, ierr)
allocate(sparse_row_sizes(num_states), stat=ierr)
call LogMemAlloc('sparse_row_sizes', num_states, sizeof_int, t_r, SRTag, ierr)
allocate(sparse_diag_positions(num_states), stat=ierr)
call LogMemAlloc('sparse_diag_positions', num_states, sizeof_int, t_r, SDTag, ierr)
allocate(indices(num_states), stat=ierr)
call LogMemAlloc('indices', num_states, sizeof_int, t_r, ITag, ierr)
! Set each element to one to count the diagonal elements straight away.
sparse_row_sizes = 1
! also need to change this routine here for the guga implementation
do i = 1, num_states
hamiltonian_row = 0.0_dp
call decode_bit_det(nI, ilut_list(:, i))
! sparse_diag_positions(i) stores the number of non-zero elements in row i of the
! Hamiltonian, up to and including the diagonal element.
! sparse_row_sizes(i) stores this number currently as all non-zero elements before
! the diagonal have been counted (as the Hamiltonian is symmetric).
sparse_diag_positions(i) = sparse_row_sizes(i)
if (tGUGA) csf_i = CSF_Info_t(ilut_list(0:nifd, i))
do j = i, num_states
call decode_bit_det(nJ, ilut_list(:, j))
if (tGUGA) csf_j = CSF_Info_t(ilut_list(:, j))
! If on the diagonal of the Hamiltonian.
if (i == j) then
if (tHPHF) then
hamiltonian_row(j) = hphf_diag_helement(nI, ilut_list(:, i))
else if (tGUGA) then
hamiltonian_row(j) = calcDiagMatEleGuga_nI(nI)
else
hamiltonian_row(j) = get_helement(nI, nJ, 0)
end if
hamil_diag(j) = hamiltonian_row(j)
else
if (tHPHF) then
hamiltonian_row(j) = hphf_off_diag_helement(nI, nJ, ilut_list(:, i), &
ilut_list(:, j))
else if (tGUGA) then
call calc_guga_matrix_element(ilut_list(:, i), csf_i, ilut_list(:, j), csf_j, &
excitInfo, hamiltonian_row(j), .true.)
#ifdef CMPLX_
hamiltonian_row(j) = conjg(hamiltonian_row(j))
#endif
! call calc_guga_matrix_element(ilut_list(:,j), ilut_list(:,i), &
! excitInfo, hamiltonian_row(j), .true., 2)
else
hamiltonian_row(j) = get_helement(nI, nJ, ilut_list(:, i), &
ilut_list(:, j))
end if
if (abs(hamiltonian_row(j)) > 0.0_dp) then
! If element is nonzero, update the following sizes.
sparse_row_sizes(i) = sparse_row_sizes(i) + 1
sparse_row_sizes(j) = sparse_row_sizes(j) + 1
end if
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(sparse_ham, i, sparse_row_sizes(i), "sparse_ham", SparseHamilTags(:, i))
sparse_ham(i)%elements = 0.0_dp
sparse_ham(i)%positions = 0
sparse_ham(i)%num_elements = sparse_row_sizes(i)
! Now fill in all matrix elements beyond and including the diagonal, as these are
! stored in hamiltonian_row.
counter = sparse_diag_positions(i)
do j = i, num_states
if (abs(hamiltonian_row(j)) > 0.0_dp) then
sparse_ham(i)%positions(counter) = j
sparse_ham(i)%elements(counter) = hamiltonian_row(j)
counter = counter + 1
end if
end do
end do
! At this point, sparse_ham has been allocated with the correct sizes, but only the
! matrix elements above and including the diagonal have been filled in. Now we must
! fill in the other elements. To do this, cycle through every element above the
! diagonal and fill in every corresponding element below it:
indices = 1
do i = 1, num_states
do j = sparse_diag_positions(i) + 1, sparse_row_sizes(i)
sparse_ham(sparse_ham(i)%positions(j))%&
&positions(indices(sparse_ham(i)%positions(j))) = i
sparse_ham(sparse_ham(i)%positions(j))%&
&elements(indices(sparse_ham(i)%positions(j))) = h_conjg(sparse_ham(i)%elements(j))
indices(sparse_ham(i)%positions(j)) = indices(sparse_ham(i)%positions(j)) + 1
end do
end do
deallocate(hamiltonian_row, stat=ierr)
call LogMemDealloc(t_r, HRTag, ierr)
deallocate(sparse_row_sizes, stat=ierr)
call LogMemDealloc(t_r, SRTag, ierr)
deallocate(sparse_diag_positions, stat=ierr)
call LogMemDealloc(t_r, SDTag, ierr)
deallocate(indices, stat=ierr)
call LogMemDealloc(t_r, ITag, ierr)
end subroutine calculate_sparse_hamiltonian