gen_all_doubles_k_space Subroutine

public subroutine gen_all_doubles_k_space(ni, n_excits, det_list, sign_list)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: ni(nel)
integer, intent(out) :: n_excits
integer(kind=n_int), intent(out), allocatable :: det_list(:,:)
real(kind=dp), intent(out), optional, allocatable :: sign_list(:)

Contents


Source Code

    subroutine gen_all_doubles_k_space(nI, n_excits, det_list, sign_list)
        integer, intent(in) :: ni(nel)
        integer, intent(out) :: n_excits
        integer(n_int), intent(out), allocatable :: det_list(:, :)
        real(dp), intent(out), allocatable, optional :: sign_list(:)
        character(*), parameter :: this_routine = "gen_all_doubles_k_space"

        integer(n_int) :: ilutJ(0:niftot), ilut(0:niftot)
        integer :: n_bound, i, j, a, b, src(2), ex(2, 2), pos, n_par, n_opp, nj(nel)
        integer(n_int), allocatable :: temp_list(:, :)
        HElement_t(dp) :: elem
        logical :: t_sign, tpar
        real(dp), allocatable :: temp_sign(:)

        ! fuck it! these annoying old routines break my balls!
        ! just write a new one to test my excitations with!
        call EncodeBitDet(nI, ilut)

        n_excits = 1

        ! i think a more correct estimat is:
        n_bound = max(int(nel * (nel - 1) * (nBasis - nel) / 4), 10)

        allocate(temp_list(0:niftot, n_bound))

        n_par = 0
        n_opp = 0

        temp_list = 0_n_int

        if (present(sign_list)) then
            t_sign = .true.
            allocate(temp_sign(n_bound), source=0.0_dp)
        else
            t_sign = .false.
        end if

        do i = 1, nel - 1
            do j = i + 1, nel
                src = nI([i, j])

                ex(1, :) = src

                ! if we do not have transcorrelation cyclce for same spins
                if (same_spin(src(1), src(2)) .and. .not. t_trans_corr_2body) cycle

                do a = 1, nbasis
                    if (IsNotOcc(ilut, a)) then
                        b = get_orb_from_kpoints(src(1), src(2), a)

                        if (IsNotOcc(ilut, b) .and. a /= b) then

                            ! do i need to sort ex?? try..
                            ex(2, :) = [a, b]

!                             if (abs(get_offdiag_helement_k_sp_hub(nI, ex, .false.)) > EPS) then
                            ! use the get_helement_lattice to avoid circular
                            ! dependencies
                            call make_double(nI, nJ, i, j, a, b, ex, tpar)
                            elem = get_helement_lattice(nI, nJ)
                            !elem = get_helement_lattice(nI, 2, ex, .false.)
                            if (abs(elem) > EPS) then

                                ilutJ = make_ilutJ(ilut, ex, 2)

                                pos = binary_search_ilut(temp_list(:, 1:n_excits), ilutJ, nifd + 1)

                                if (pos < 0) then

                                    temp_list(:, n_excits) = ilutJ

                                    if (t_sign) then
#ifdef CMPLX_
                                        call stop_all(this_routine, "not implemented for complex")
#else
                                        temp_sign(n_excits) = sign(1.0_dp, elem)
                                        call sort(temp_list(:, 1:n_excits), temp_sign(1:n_excits))
#endif
                                    else
                                        call sort(temp_list(:, 1:n_excits), ilut_lt, ilut_gt)
                                    end if

                                    n_excits = n_excits + 1
                                    ! damn.. i have to sort everytime i guess..

                                    ! to be sure also count seperately if it is
                                    ! a parallel or opposite excitation
                                    if (same_spin(src(1), src(2))) then
                                        ASSERT(same_spin(src(1), a))
                                        ASSERT(same_spin(a, b))
                                        n_par = n_par + 1
                                    else
                                        ASSERT(.not. same_spin(a, b))
                                        n_opp = n_opp + 1
                                    end if
                                end if
                            end if
                        end if
                    end if
                end do
            end do
        end do

        n_excits = n_excits - 1
        allocate(det_list(0:NIfTot, n_excits), source=temp_list(:, 1:n_excits))

        if (t_sign) then
            allocate(sign_list(n_excits), source=temp_sign(1:n_excits))
            call sort(det_list, sign_list)!, ilut_lt, ilut_gt)
        else
            call sort(det_list, ilut_lt, ilut_gt)
        end if

    end subroutine gen_all_doubles_k_space