actHamiltonian Subroutine

public subroutine actHamiltonian(ilut, csf_i, excitations, nTot, t_singles_only, t_print_time, t_full)


Type IntentOptional Attributes Name
integer(kind=n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer(kind=n_int), intent(out), allocatable :: excitations(:,:)
integer, intent(out) :: nTot
logical, intent(in), optional :: t_singles_only
logical, intent(in), optional :: t_print_time
logical, intent(in), optional :: t_full


Source Code

Source Code

    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.
            t_singles_only_ = t_singles_only
        end if


        ! 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

                ! 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

                        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
                            if (nExcits > 0) then
                                call add_guga_lists(nTot, nExcits, tmp_all_excits, &
                            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
            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)

        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