get_diag_helement_heisenberg Function

public function get_diag_helement_heisenberg(nI) result(hel)


Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)

Return Value real(kind=dp)


Source Code

    function get_diag_helement_heisenberg(nI) result(hel)
        integer, intent(in) :: nI(nel)
        HElement_t(dp) :: hel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_diag_helement_heisenberg"
        integer :: i, j, src, flip
        integer, allocatable :: spin_neighbors(:)
        integer(n_int) :: ilut(0:NIfTot)

        ! here i have to sum over all the nearest neigbhor pairs and
        ! find the spin alignment, if it is parallel or not..
        ! in the heisenberg case the contribution form n_i n_j form
        ! nearest neighbors is always the same for every determinant nI
        ! but for the tJ model with less than half-filling this quantitiy is
        ! not a constant.. so i should differentiate in here between tJ
        ! and heisenberg models ..
        ! and what about double counting? should i loop over all the
        ! orbitals and their neighbors or must i differentiate between already
        ! counted contributions?

        ! it is easier to check occupancy in the ilut format


        call EncodeBitDet(nI, ilut)
        hel = h_cast(0.0_dp)

        do i = 1, nel
            src = ni(i)
            spin_neighbors = lat%get_spinorb_neighbors(src)
            if (is_beta(src)) then
                flip = +1
                flip = -1
            end if
            do j = 1, size(spin_neighbors)
                if (IsOcc(ilut, spin_neighbors(j))) then
                    ! then it is same spin
                    ! but i really think that we have to take the
                    ! occupancy into account
                    ! and in the tJ model this cancels for parallel spin
                    if (t_heisenberg_model) then
                        hel = hel + h_cast(exchange_j / 4.0_dp)
                    end if
                else if (IsOcc(ilut, spin_neighbors(j) + flip)) then
                    if (t_heisenberg_model) then
                        hel = hel - h_cast(exchange_j / 4.0_dp)
                        hel = hel - h_cast(exchange_j / 2.0_dp)
                    end if
                    ! it can be empty too, then there is no contribution
                end if
            end do

        end do

        ! i think i would double count if i do not divide by 2..
        ! i hope i am right..
        hel = hel / 2.0_dp

    end function get_diag_helement_heisenberg