PickElecPair Subroutine

    SUBROUTINE PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, IndInp)
        INTEGER, INTENT(IN) :: nI(NEl), IndInp
        INTEGER, INTENT(OUT) :: Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl
        INTEGER :: Ind, X, K, Orb1, Orb2
        real(dp) :: r
!Triangular indexing system.
!This is used for picking two distinct electrons out of all N(N-1)/2 pairs.
!
!   12  13  14  15          1   2   3   4
!       23  24  25              5   6   7
!           34  35      =>          8   9
!               45                      10

!For a given index, ind, there are [N(N-1)/2 - ind] elements at positions larger
!than ind. Call this number X.
!We want to find out how many rows there are after the row containing the element
!ind. t rows has t(t+1)/2 elements in it.
!Therefore, to find out the number of rows after ind, we want to find the largest K, such that K(K+1)/2 =< X
!The solution to this is that K =< (SQRT(8*X+1)-1)/2, therefore we can find K (the
!largest integer which satisfies this).
!We then know the number of rows after the element ind. Therefore, since there are N-1
!rows in total, we know we are on row N-1-K.
!This gives us the index of the first electron.
!To find the second electron (i.e. the column), we know that out of the X elements in
!positions larger than ind, K(K+1)/2 are in the next rows.
!This means that X - K(K+1)/2 are in the same row. There are N-(N-1-K) = 1+K elements
!in the row chosen, and so the number of elements into the
!row it is, is given by (1+K) - (X-K(K+1)/2). However, in row z, the first column
!index is z+1. Therefore the index of the second electron is
!(1+K) - (X-K(K+1)/2) + N-K-1 = N-X+(K(K+1)/2).

!        ElecPairs=(NEl*(NEl-1))/2

! If we want to find an index randomly, IndInp will be -1.
        IF (IndInp == -1) THEN
!Find an index randomly.
            r = genrand_real2_dSFMT()
            Ind = INT(ElecPairs * r) + 1
        ELSE
!Otherwise we are looking for a specific electron pair specified by IndInp
            Ind = IndInp
        end if

!X is number of elements at positions larger than ind
        X = ElecPairs - Ind
!K is the number of complete rows after the element ind
        K = INT((SQRT(8.0_dp * REAL(X, dp) + 1.0_dp) - 1.0_dp) / 2.0_dp)
        Elec1Ind = NEl - 1 - K
        Elec2Ind = NEl - X + ((K * (K + 1)) / 2)

        Orb1 = nI(Elec1Ind)
        Orb2 = nI(Elec2Ind)

!We now want to find the symmetry product label of the two electrons, and the spin product of the two electrons.
        SymProduct = RandExcitSymLabelProd(SpinOrbSymLabel(Orb1), SpinOrbSymLabel(Orb2))

        IF ((G1(Orb1)%Ms) * (G1(Orb2)%Ms) == -1) THEN
!We have an alpha beta pair of electrons.
            iSpn = 2
        ELSE
            IF (G1(Orb1)%Ms == 1) THEN
!We have an alpha alpha pair of electrons.
                iSpn = 3
            ELSE
!We have a beta beta pair of electrons.
                iSpn = 1
            end if
        end if

        SumMl = G1(Orb1)%Ml + G1(Orb2)%Ml

    END SUBROUTINE PickElecPair