SUBROUTINE GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, tij_lt_ab_only)
! This generates one by one, all possible double excitations.
! This involves a way of ordering the electron pairs i,j and a,b so that given an i,j and a,b we can find the next.
! The overall symmetry must also be maintained - i.e. if i and j are alpha and beta, a and b must be alpha and beta
! or vice versa.
USE SystemData, only: ElecPairs, tFixLz, iMaxLz
USE GenRandSymExcitNUMod, only: PickElecPair
use constants, only: bits_n_int
INTEGER :: nI(NEl), Orbj, Orbi, Orba, Orbb, Syma, Symb
INTEGER(KIND=n_int) :: iLut(0:NIfTot)
INTEGER :: Elec1Ind, Elec2Ind, SymProduct, iSpn, Spinb, nJ(NEl), ExcitMat3(2, 2), SumMl
INTEGER, SAVE :: ijInd, OrbaChosen, OrbbIndex, Spina, SymInd
LOGICAL :: tDoubleExcitFound, tFirsta, tFirstb, tNewij, tNewa, tAllExcitFound, tParity, tij_lt_ab_only
INTEGER :: Mla, Mlb, Indij
tDoubleExcitFound = .false.
tFirsta = .false.
tFirstb = .false.
Orbi = ExcitMat3(1, 1)
Orbj = ExcitMat3(1, 2)
Orba = ExcitMat3(2, 1)
Orbb = ExcitMat3(2, 2)
IF (Orbi == 0) THEN
ijInd = 1
! If Orbi, then we are choosing the first double.
! It is therefore also the first set of a and b for this electron pair i,j.
tFirsta = .true.
tFirstb = .true.
end if
lp: do while (.not. tDoubleExcitFound)
! Otherwise we use the previous ijInd and the saved indexes for a and b.
! This routine allows us to pick an electron pair i,j specified by the index ijInd.
! The i and j orbitals are then given by nI(Elec1Ind) and nI(Elec2Ind), and the symmetry product of the two is
! SymProduct and the spin iSpn.
! iSpn=2 for alpha beta pair, ispn=3 for alpha alpha pair and ispn=1 for beta beta pair.
CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, ijInd)
Indij = ((((nI(Elec2Ind) - 2) * (nI(Elec2Ind) - 1)) / 2) + nI(Elec1Ind))
tNewij = .false.
! This becomes true when we can no longer find an allowed orbital a for this ij pair and we need to move onto the next.
do while ((.not. tNewij) .and. (.not. tDoubleExcitFound)) ! This loop runs through the allowed a orbitals
! until a double excitation is found.
IF (tFirsta) THEN
! If this is the first double we are picking with this ij, we start with the alpha spin, unless i and j are both beta.
! There is no restriction on the symmetry for orbital a - although clearly the symmetry we pick determins b.
! write(stdout,*) "iSpn",iSpn
IF (iSpn == 1) THEN
!beta beta pair
Spina = 2
!Want to start at the first beta orbital
OrbaChosen = 1
else if (iSpn == 3) THEN
!alpha alpha pair
Spina = 1
OrbaChosen = 2
ELSE
Spina = 2
OrbaChosen = 1
end if
end if
! If it is not the first, we have stored the previous spina and orba index - need to start with these and see
! if any more double remain.
Orba = OrbaChosen
! write(stdout,*) "Chosen index, orbital for a: ",OrbaChosen,Orba
! The orbital chosen must be unoccupied. This is just a test to make sure this is the case.
do while ((BTEST(iLut((Orba - 1) / bits_n_int), MOD((Orba - 1), bits_n_int))) .or. (abs(SumMl - G1(Orba)%Ml) > iMaxLz))
!We also test that the summl value and ml of Orba is such that it is possible for orbb to conserve ml.
!Will get into this loop if the orbital is occupied, or if the ml is such that no orbb is possible.
! If not, we move onto the next orbital.
IF (iSpn /= 2) THEN
!Increment by two, since we want to look at the same spin state.
OrbaChosen = OrbaChosen + 2
ELSE
!Increment by one, since we want to look at both alpha and beta spins.
OrbaChosen = OrbaChosen + 1
IF (Spina == 2) THEN
Spina = 1
ELSE
Spina = 2
end if
end if
IF (OrbaChosen > nBasis) THEN
!We have reached the end of all allowed symmetries for the a orbital, only taking
!into account spin symmetry. Choose new ij pair now.
tNewij = .true.
EXIT
end if
! Otherwise the new orbital a is the first unoccupied orbital of allowed symmetry etc.
Orba = OrbaChosen
! write(stdout,*) "Chosen index, orbital for a: ",OrbaChosen,Orba
end do
! If we have got to the end of the a orbitals, and need a new i,j pair, we increment ijInd and check
! this hasn't gone beyond the limits - bail out if we have.
IF (tNewij) THEN
ijInd = ijInd + 1
IF (ijInd > ElecPairs) THEN
tDoubleExcitFound = .true.
tAllExcitFound = .true.
! AllExcitFound true indicates there are no more symmetry allowed double excitations.
end if
EXIT
end if
tNewa = .false.
!Find a b
do while ((.not. tNewa) .and. (.not. tDoubleExcitFound))
! We now have i,j,a and we just need to pick b.
! First find the spin of b.
IF (iSpn == 1) THEN
Spinb = 2
else if (iSpn == 3) THEN
Spinb = 1
ELSE
IF (Spina == 1) THEN
Spinb = 2
ELSE
Spinb = 1
end if
end if
! Then find the symmetry of b.
IF (tNoSymGenRandExcits) THEN
Syma = 0
ELSE
Syma = INT(G1(Orba)%Sym%S, 4)
end if
Symb = IEOR(Syma, SymProduct)
! Then find the ml of b.
IF (tFixLz) THEN
Mla = G1(Orba)%Ml
Mlb = SumMl - Mla
ELSE
Mla = 0
Mlb = 0
end if
! If this is the first time we've picked an orbital b for these i,j and a, begin at the start of the symmetry block.
! Otherwise pick up where we left off last time.
IF (tFirstb) THEN
SymInd = ClassCountInd(Spinb, Symb, Mlb)
OrbbIndex = SymLabelCounts2(1, SymInd)
ELSE
!Update new orbital b index
OrbbIndex = OrbbIndex + 1
end if
! If the new b orbital is still within the limits, check it is unoccupied and move onto the next orbital if it is.
IF (OrbbIndex > nBasis) THEN
tNewa = .true.
tFirsta = .false.
end if
IF (.not. tNewa) THEN
IF (OrbbIndex > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN
! If we have already gone beyond the symmetry limits by choosing the next b orbital, pick a new a orbital.
tNewa = .true.
tFirsta = .false.
ELSE
Orbb = SymLabelList2(OrbbIndex)
! Checking the orbital b is unoccupied and > a.
do while (((BTEST(iLut((Orbb - 1) / bits_n_int), MOD((Orbb - 1), bits_n_int))) .or. (Orbb <= Orba)) .or. &
(tij_lt_ab_only .and. (((((Orbb - 2) * (Orbb - 1)) / 2) + Orba) < Indij)))
!Orbital is occupied - try again
OrbbIndex = OrbbIndex + 1
IF (OrbbIndex > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN
!Reached end of symmetry block - need new a
! write(stdout,*) "Reached end of sym block",Orbb,Orba
tNewa = .true.
tFirsta = .false.
EXIT
end if
! write(stdout,*) "Cycling through orbitals: ",OrbbIndex,Orbb
!Update new orbital b index
Orbb = SymLabelList2(OrbbIndex)
! write(stdout,*) "Attempting again with orbital: ",Orbb
end do
end if
end if
! If we are moving onto the next a orbital, check we don't also need a new ij pair.
IF (tNewa) THEN
IF (iSpn /= 2) THEN
!Increment by two, since we want to look at the same spin state.
OrbaChosen = OrbaChosen + 2
! tFirsta=.false.
ELSE
!Increment by one, since we want to look at both alpha and beta spins.
IF (Spina == 1) THEN
Spina = 2
ELSE
Spina = 1
end if
OrbaChosen = OrbaChosen + 1
! tFirsta=.false.
end if
tFirstb = .true.
! write(stdout,*) "New OrbaChosen: ",OrbaChosen
IF (OrbaChosen > nBasis) THEN
!We have reached the end of all allowed symmetries for the a orbital, only taking
!into account spin symmetry. Choose new ij pair now.
tNewij = .true.
ijInd = ijInd + 1
! write(stdout,*) "ijInd: ",ijInd
IF (ijInd > ElecPairs) THEN
tAllExcitFound = .true.
tDoubleExcitFound = .false.
EXIT lp
end if
end if
ELSE
!If we don't need a new a, we have found an excitation ij -> ab that is accepted.
tDoubleExcitFound = .true.
end if
end do
end do
! This is the loop for new ij pairs - if we are choosing a new ij we are automatically choosing a new a and b also.
tFirsta = .true.
tFirstb = .true.
end do lp
if (tDoubleExcitFound .and. (.not. tAllExcitFound)) then
call make_double(nI, nJ, elec1ind, elec2ind, orbA, orbB, &
ExcitMat3, tParity)
end if
ENDSUBROUTINE GenDoubleExcit