gen_uniform_double_anti Subroutine

private subroutine gen_uniform_double_anti(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel)
integer(kind=n_int), intent(out) :: ilutJ(0:niftot)
integer, intent(out) :: ex(2,3)
logical, intent(out) :: tParity
real(kind=dp), intent(out) :: pgen

Contents


Source Code

    subroutine gen_uniform_double_anti(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)
        ! routine to do a uniform anti-parallel spin double excitation
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: nJ(nel), ex(2, 3)
        integer(n_int), intent(out) :: ilutJ(0:niftot)
        logical, intent(out) :: tParity
        real(dp), intent(out) :: pgen
        integer :: a, b, i, elecs(2)
        integer, parameter :: max_trials = 1000
        real(dp) :: p_elec, p_orb

        call pick_spin_opp_elecs(nI, elecs, p_elec)

        p_orb = 2.0_dp / real(nbasis - nel, dp)

        nJ(1) = 0

        ! pick 2 holes now
        do i = 1, max_trials

            a = int(nbasis * genrand_real2_dsfmt()) + 1

            if (IsOcc(ilutI, a)) cycle

            b = get_orb_from_kpoints(nI(elecs(1)), nI(elecs(2)), a)

            ! do we have to reject or can we cycle if not fitting?
            ! a == b test has to be here for the spin-parallel
            ! excitations!
            if (IsOcc(ilutI, b) .or. a == b) then
                nJ(1) = 0
                return
            end if

            call make_double(nI, nJ, elecs(1), elecs(2), a, b, ex, tParity)

            ilutJ = make_ilutJ(ilutI, ex, 2)
            exit
        end do

        pgen = p_elec * p_orb

    end subroutine gen_uniform_double_anti