SUBROUTINE GenSingleExcit(nI, iLut, nJ, exflag, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only)
! Despite being fed four indices, this routine finds single excitations. Orbi -> Orba. (Orbj and Orbb remain 0).
! Feeding in 0 indices indicates it is the first excitation that needs to be found.
! The single excitation goes from orbital i to a, from determinant nI to nJ.
! When the last single is found it then finds the first double excitation, unless exflag=1 in which tAllExcitFound
! becomes true and no more excitations are generated.
use SystemData, only: tFixLz
use constants, only: bits_n_int
INTEGER :: nI(NEl), Orbi, Orba, Symi, nJ(NEl)
INTEGER(KIND=n_int) :: iLut(0:NIfTot)
INTEGER :: NoOcc, ExcitMat3(2, 2), exflag, SymInd, Spina, Mla
LOGICAL :: tInitOrbsFound, tParity, tAllExcitFound, tEndaOrbs, ti_lt_a_only, tAux
INTEGER, SAVE :: OrbiIndex, OrbaIndex, Spini, NewSym, Mli
tInitOrbsFound = .false.
Orbi = ExcitMat3(1, 1)
Orba = ExcitMat3(2, 1)
IF ((Orbi == 0) .or. (Orba == 0)) THEN ! Want to find the first excitation.
OrbiIndex = 1
Orbi = nI(OrbiIndex) ! Take the first occupied orbital
Symi = SpinOrbSymLabel(Orbi) ! and find its spin and spat symmetries.
IF ((G1(Orbi)%Ms) == -1) Spini = 2
IF ((G1(Orbi)%Ms) == 1) Spini = 1
IF (tFixLz) THEN
Mli = G1(Orbi)%Ml
ELSE
Mli = 0
end if
OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli)) ! Start considering a at the first allowed symmetry.
ELSE
Orbi = nI(OrbiIndex) ! Begin by using the same i as last time - check if there are any
! more possible excitations from this.
! At this stage, OrbaIndex is the a from the previous excitation.
SymInd = ClassCountInd(Spini, SpinOrbSymLabel(Orbi), Mli)
IF (OrbaIndex == (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN
!Orba was the last in the symmetry block. Do not allow OrbaIndex+1
! Either we're got to the final spin symmetry, or the next orbital after
!Orba does not have the same symmetry as Orbi.
! Need to move onto the next i, and find a new a to match.
OrbiIndex = OrbiIndex + 1
IF (OrbiIndex <= NEl) THEN
Orbi = nI(OrbiIndex)
Symi = SpinOrbSymLabel(Orbi)
IF ((G1(Orbi)%Ms) == -1) Spini = 2
IF ((G1(Orbi)%Ms) == 1) Spini = 1
IF (tFixLz) THEN
Mli = G1(Orbi)%Ml
ELSE
Mli = 0
end if
OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli))
ELSE
IF (exflag /= 1) THEN
ExcitMat3(:, :) = 0
CALL GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only)
exflag = 2
ELSE
tAllExcitFound = .true.
tInitOrbsFound = .true.
end if
end if
ELSE
! There are more possible excitations from orbital a, simply check the next orbital after the current a.
OrbaIndex = OrbaIndex + 1
Symi = SpinOrbSymLabel(Orbi)
end if
end if
do while (.not. tInitOrbsFound)
tEndaOrbs = .false.
IF (OrbiIndex > NEl) THEN
! If we've read in the last single, set orbi, orbj, orba, and orbb to 0 and call gendoubleexcit.
IF (exflag /= 1) THEN
ExcitMat3(:, :) = 0
CALL GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only)
exflag = 2
ELSE
tAllExcitFound = .true.
end if
EXIT
end if
! To find Orba, take the first in SymLabelList2 with the same symmetry and spin.
! SymLabelCounts2(spin,1,symmetry) gives the index in SymLabelList2 where that spin and symmetry starts.
IF (OrbaIndex > nBasis) THEN
tEndaOrbs = .true.
ELSE
tEndaOrbs = .false.
Orba = SymLabelList2(OrbaIndex)
end if
SymInd = ClassCountInd(Spini, SpinOrbSymLabel(Orbi), Mli)
! Need to also make sure orbital a is unoccupied, so make sure the orbital is not in nI.
NoOcc = 0
IF (.not. tEndaOrbs) THEN
do while ((BTEST(iLut((Orba - 1) / bits_n_int), MOD((Orba - 1), bits_n_int))) .or. &
(ti_lt_a_only .and. (Orba < Orbi)))
! While this is true, Orba is occupied, so keep incrementing Orba until it is not.
NoOcc = NoOcc + 1
IF (OrbaIndex + NoOcc > nBasis) THEN
!We have reached the end of all a orbitals. Now we need to pick a new i
tEndaOrbs = .true.
EXIT
ELSE
Orba = SymLabelList2(OrbaIndex + NoOcc)
IF ((OrbaIndex + NoOcc) > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) EXIT
end if
end do
end if
tAux = .false.
IF (.not. tEndaOrbs) THEN
! Then check we have not overrun the symmetry block while skipping the occupied orbitals.
NewSym = SpinOrbSymLabel(Orba)
IF ((G1(Orba)%Ms) == -1) Spina = 2
IF ((G1(Orba)%Ms) == 1) Spina = 1
IF (tFixLz) THEN
Mla = G1(Orba)%Ml
ELSE
Mla = 0
end if
IF (NewSym == Symi .and. (Spina == Spini) .and. (Mli == Mla)) THEN
! If not, then these are the new Orbi and Orba.
tInitOrbsFound = .true.
OrbaIndex = OrbaIndex + NoOcc
ELSE
tAux = .true.
end if
end if
! If we have, move onto the next occupied orbital i, no symmetry allowed single excitations exist from the first.
IF (tAux .or. tEndaOrbs) THEN
OrbiIndex = OrbiIndex + 1
IF (OrbiIndex <= NEl) THEN
Orbi = nI(OrbiIndex)
Symi = SpinOrbSymLabel(Orbi) ! and find its spin and spat symmetries.
IF ((G1(Orbi)%Ms) == -1) Spini = 2
IF ((G1(Orbi)%Ms) == 1) Spini = 1
IF (tFixLz) THEN
Mli = G1(Orbi)%Ml
ELSE
Mli = 0
end if
OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli))
end if
end if
end do
IF ((ExcitMat3(1, 2) == 0) .and. (.not. tAllExcitFound)) then
CALL FindNewSingDet(nI, nJ, OrbiIndex, OrbA, ExcitMat3, tParity)
end if
ENDSUBROUTINE GenSingleExcit