PickElecPair Subroutine

public subroutine PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, IndInp)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(NEl)
integer, intent(out) :: Elec1Ind
integer, intent(out) :: Elec2Ind
integer, intent(out) :: SymProduct
integer, intent(out) :: iSpn
integer, intent(out) :: SumMl
integer, intent(in) :: IndInp

Source Code

    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