subroutine actHamiltonian(ilut, csf_i, excitations, nTot, t_singles_only, &
t_print_time, t_full)
! subroutine to calculate the action of the full Hamiltonian on a
! a single CSF given in ilut bit representation and outputs a list
! of excitations also in ilut format, where the exact matrix element
! is stored, where usually the signed walker weight of an occupied
! determinant is stored
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer(n_int), intent(out), allocatable :: excitations(:, :)
integer, intent(out) :: nTot
logical, intent(in), optional :: t_singles_only, t_print_time, t_full
character(*), parameter :: this_routine = "actHamiltonian"
type(timer), save :: proc_timer
integer(n_int), allocatable :: tmp_all_excits(:, :)
integer :: i, j, k, l, nExcits, nMax, ierr
integer(n_int), allocatable :: tempExcits(:,:)
logical :: t_singles_only_, t_print_time_, t_full_
integer :: n
real(dp) :: cmp
def_default(t_singles_only_, t_singles_only, .false.)
def_default(t_print_time_, t_print_time, .false.)
def_default(t_full_, t_full, .true.)
if (.not. present(t_singles_only)) then
t_singles_only_ = .false.
else
t_singles_only_ = t_singles_only
end if
ASSERT(isProperCSF_ilut(ilut))
! finally add some timing routines to this as it gets quite heavy
proc_timer%timer_name = this_routine
call set_timer(proc_timer)
! essentially have to calculate the diagonal matrix element of ilut
! and additionally loop over all possible single excitaiton (i,j)
! and double excitation (i,j,k,l) index combinations and for them
! calculate all possible excitations and matrix elements
! have to somehow store all the already created excitations,...
! so I have to initialize an array somehow... how exactly am i going
! to do that... the ilut list combinig routine says there should be no
! repetitions in the list, but if i initialize a huge array to zero
! and fill that up, does that cause problems? -> ask simon/ali
! and how big does this array have to be?
! for each excitation 2**nOpenOrbs and there are at worst
! (nOrbitals/2)**4 excitation, but not all with the full excitation
! range.. but for now initiate it with this worse than worst case
! hm it seems i can give a number of determinants input to
! routine add_ilut_lists, so maybe i constrain the number of
! processed states with that
nMax = 6 + 4 * (nSpatOrbs)**4 * (count_open_orbs(ilut) + 1)
! todo: get an exact maximum formula for that which holds in all cases!
! because that number above is WAYY too big!
! also for excitations where the indices are known!
! and fot those pescy non-overlap excitations!
allocate(tmp_all_excits(0:nifguga, nMax), stat=ierr)
call LogMemAlloc('tmp_all_excits', (nifguga + 1) * nMax, 8, this_routine, tag_tmp_excits)
! maybe have to set to zero here then, or initialize somehow different
nTot = 0
! single excitations:
! does it help to not calculate stuff if not necessary..
! probably yes..
if (.not. ((thub .and. .not. treal) .or. t_heisenberg_model)) then
do i = 1, nSpatOrbs
do j = 1, nSpatOrbs
if (i == j) cycle ! do not calc. diagonal elements here
! need integral contributions too...
! maybe only loop over non-zero index combinations or
! otherwise a lot of effort wasted ...
! and i might have to loop over created excitations here too
! to multiply with the integrals... better to do that
! within the excitation creation!
! have to think where and when to allocate excitations
! or maybe call the routine below with a dummy variable.
! and store calculated ones in a different variable here..
call calcAllExcitations(ilut, csf_i, i, j, tempExcits, nExcits, t_full_)
#ifdef DEBUG_
do n = 1, nExcits
ASSERT(isProperCSF_ilut(tempExcits(:, n), .true.))
end do
#endif
! merge created lists.. think how that will be implemented
if (nExcits > 0) then
call add_guga_lists(nTot, nExcits, tmp_all_excits, tempExcits)
! nTot gets automatically updated too
end if
end do
end do
end if
! double excitations
! do it really primitive for now. -> make it more elaborate later
if (.not. t_singles_only_) then
if (.not. (thub .and. treal)) then
do i = 1, nSpatOrbs
do j = 1, nSpatOrbs
do k = 1, nSpatOrbs
do l = 1, nSpatOrbs
if (i == j .and. k == l) cycle
! if (t_heisenberg_model) then
! if (.not. (i == l .and. j == k)) cycle
! end if
! not only i=j and k=l index combinations can lead to
! diagonal contributions but also i = l and k = j
! mixed fullstart-> fullstop excitations can have a
! diagonal term -> not only can, but always will have
! so i have to exclude this term in the exact excitation
! calculation when states get added up, since otherwise
! i would count these terms to often, as they are
! already taken into account in the diagonal term
! calculation
! integral contributions
call calcAllExcitations(ilut, csf_i, i, j, k, l, tempExcits, &
nExcits, t_full_)
#ifdef DEBUG_
do n = 1, nExcits
if (.not. isProperCSF_ilut(tempExcits(:, n), .true.)) then
call write_guga_list(stdout, tempExcits(:, 1:nExcits))
print *, "ijkl:", i, j, k, l
end if
ASSERT(isProperCSF_ilut(tempExcits(:, n), .true.))
end do
#endif
if (i == l .and. k == j) then
! exclude the possible diagonal term
! i assume for now that this term is always there
! atleast the x0-matrix element is never zero for
! these diagonal excitations -> so check if it
! stays in the same position: yes it stays in
! first position always..
! for non-open orbitals i
! NO assert since i exclude excitations where
! x1 is always 0
! or otherwise somthings wrong... and then only
! update the list if another excitation is
! encountered
if (nExcits > 1) then
call add_guga_lists(nTot, nExcits - 1, tmp_all_excits, &
tempExcits(:, 2:))
end if
else
if (nExcits > 0) then
call add_guga_lists(nTot, nExcits, tmp_all_excits, &
tempExcits)
end if
end if
end do
end do
end do
end do
end if
end if
! do an additional check if matrix element is 0
j = 1
if (t_matele_cutoff) then
cmp = matele_cutoff
else
cmp = EPS
end if
do i = 1, nTot
if (abs(extract_h_element(tmp_all_excits(:, i))) < cmp) cycle
tmp_all_excits(:, j) = tmp_all_excits(:, i)
j = j + 1
end do
nTot = j - 1
allocate(excitations(0:nifguga, nTot), stat=ierr)
! hm to log that does not make so much sense.. since it gets called
! more than once and is only a temporary array..
call LogMemAlloc('excitations', nTot, 8, this_routine, tag_excitations)
excitations = tmp_all_excits(:, 1:nTot)
deallocate(tmp_all_excits)
call LogMemDealloc(this_routine, tag_tmp_excits)
! allocate the outputted excitations here since otherwise i store
! a way too big array than necessary
! i probably should resize the excitaions array.. since usually nMax
! is way bigger then nTot..
call halt_timer(proc_timer)
if (t_print_time_) then
write(stdout, *) " Exact Hamiltonian application done! "
write(stdout, *) " Elapsed time: ", get_total_time(proc_timer)
call neci_flush(stdout)
end if
end subroutine actHamiltonian