pick_three_elecs Subroutine

    subroutine pick_three_elecs(nI, elecs, src, sym_prod, pgen, ms)
        ! picks three random electrons from nI, biased towards 0, 1, 2 or 3 beta electrons
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: elecs(3), src(3), sym_prod, ms
        real(dp), intent(out) :: pgen
        integer :: ispn, sum_ml
        integer :: nOcc
        real(dp) :: r
        logical :: pickAlpha

        if (tContact) then
            call pick_oppspin_elecs(nI, elecs(1:2), src(1:2), sym_prod, ispn, sum_ml, pgen)
        else
            call pick_biased_elecs(nI, elecs(1:2), src(1:2), sym_prod, ispn, sum_ml, &
                                   pgen, p0B + p0A, p0B)
        end if

        if (ispn == 3 .or. ispn == 1) then
            ! all elecs have the same spin
            pickAlpha = ispn == 3
            if (pickAlpha) then
                nOcc = nOccAlpha
                ms = 3
            else
                nOcc = nOccBeta
                ms = -3
            end if
            call get_missing_elec(nI, elecs, nOcc, 2, pickAlpha, pgen)
            pgen = pgen * 3.0_dp
        else
            ! first picked one alpha and the beta
            r = genrand_real2_dSFMT()
            if (r < p2B) then
                ! we have 2 beta elecs + 1 alpha
                ! then pick the second beta
                call get_missing_elec(nI, elecs, nOccBeta, 1, .false., pgen)
                ms = -1
                pgen = pgen * p2B / (1.0_dp - (p0B + p0A))
            else
                ! we have 2 alpha elecs + 1 beta
                ! then pick the second alpha
                call get_missing_elec(nI, elecs, nOccAlpha, 1, .true., pgen)
                ms = 1
                pgen = pgen * p1B / (1.0_dp - (p0B + p0A))
            end if
            pgen = pgen * 2.0_dp
        end if

        ! sort the generated electrons
        elecs = sort_unique(elecs)
        ! check for invalid excitation
        ! after sorting, invalid electrons are definitly at position 1
        if (elecs(1) == 0) then
            src = 0
        else
            src = nI(elecs)
        end if

    end subroutine pick_three_elecs