SUBROUTINE PickAOrb(nI, iSpn, ILUT, ClassCountUnocc2, NExcit, Elec1Ind, Elec2Ind, SpinOrbA, OrbA, SymA, SymB, &
& SymProduct, SumMl, MlA, MlB, ForbiddenOrbs, tAOrbFail)
INTEGER, INTENT(IN) :: nI(NEl), iSpn, Elec1Ind, Elec2Ind, ForbiddenOrbs, SymProduct, SumMl, ClassCountUnocc2(ScratchSize)
INTEGER, INTENT(OUT) :: SpinOrbA, SymA, MlA, MlB, NExcit, SymB, OrbA
INTEGER(KIND=n_int), INTENT(IN) :: ILUT(0:NIfTot)
LOGICAL, INTENT(OUT) :: tAOrbFail
INTEGER :: AttemptsOverall, ChosenUnocc, z, i, Attempts
real(dp) :: r
IF (iSpn == 2) THEN
!There is no restriction on whether we choose an alpha or beta spin, so there are nBasis-NEl possible spinorbitals to choose from.
!Therefore the spinOrbA variable has to be set after we have chosen one. For the other iSpn types, we can set it now.
NExcit = nBasis - NEl
else if (iSpn == 1) THEN
!SpinOrbA indicates the spin of the chosen A orb. This is only really useful for iSpn=2
SpinOrbA = -1 !Going to pick a beta orb
NExcit = (nBasis / 2) - nOccBeta !This is the number of unoccupied beta spinorbitals there are.
ELSE !iSpn is 3 - alpha/alpha pair.
SpinOrbA = 1 !Going to pick an alpha orb
NExcit = (nBasis / 2) - nOccAlpha !This is the number of unoccupied alpha spinorbitals there are.
end if
IF (NExcit == ForbiddenOrbs) THEN
tAOrbFail = .true.
RETURN
ELSE
tAOrbFail = .false.
end if
AttemptsOverall = 0
do while (.true.)
!Keep drawing unoccupied orbitals, until we find one which has an allowed partner to form a symmetry-allowed unoccupied pair.
IF (iSpn == 2) THEN
!Electrons chosen were an alpha/beta pair, therefore first randomly chosen orbital can be an alpha OR beta orbital - no restriction.
IF (tCycleOrbs .or. ((NExcit - ForbiddenOrbs) <= 3)) THEN
! METHOD 1 (Run though all orbitals to find desired one out of NExcit-ForbiddenOrbs...now
!only one random number needed, and guarenteed to find excitation.)
! ==========================
!Choose the unoccupied orbital to excite to
r = genrand_real2_dSFMT()
ChosenUnocc = INT((NExcit - ForbiddenOrbs) * r) + 1
z = 1
do i = 0, nBasis - 1
!We need to find if allowed
IF (.not. BTEST(ILUT(i / bits_n_int), MOD(i, bits_n_int))) THEN
!Is not in the original determinant - allow if sym allowed
!Now check that its symmetry allowed
SpinOrbA = G1(i + 1)%Ms
IF (IsAOrbSymAllowed(iSpn, i + 1, SpinOrbA, SymProduct, SumMl, SymA, SymB, MlA, MlB, ClassCountUnocc2)) THEN
IF (z == ChosenUnocc) THEN
!We have found our allowed orbital
OrbA = i + 1
RETURN
end if
z = z + 1
end if
end if
end do
CALL Stop_All("PickAOrb", "Should not get here - have not found allowed A Orb")
ELSE
! METHOD 2 (Keep drawing orbitals randomly until we find one unoccupied). This should
!be more efficient, unless we have v. small basis sets.
! =========================
Attempts = 0
do while (.true.)
!Draw randomly from the set of orbitals
r = genrand_real2_dSFMT()
ChosenUnocc = INT(nBasis * r)
!Find out if orbital is in nI or not. Accept if it isn't in it.
IF (.not. (BTEST(ILUT(ChosenUnocc / bits_n_int), MOD(ChosenUnocc, bits_n_int)))) THEN
!Orbital not in nI. Accept.
EXIT
end if
IF (Attempts > 250) THEN
write(stdout, *) "Cannot find A unoccupied orbital after 250 attempts..."
call write_det(stdout, nI, .TRUE.)
CALL Stop_All("PickAOrb", "Cannot find A unoccupied orbital after 250 attempts...")
end if
Attempts = Attempts + 1
end do
OrbA = ChosenUnocc + 1 !This is the allowed orbital
end if
SpinOrbA = G1(OrbA)%Ms
ELSE
!We are either constrained to choose a beta orbital, or an alpha orbital - we know that there are NExcit of these.
IF (tCycleOrbs .or. ((NExcit - ForbiddenOrbs) <= 3)) THEN
! METHOD 1 (Run though all orbitals to find desired one out of NExcit-ForbiddenOrbs)
! ==========================
r = genrand_real2_dSFMT()
ChosenUnocc = INT((NExcit - ForbiddenOrbs) * r) + 1
SpinOrbA = 0 !The spin of the chosen A is not important, since we have already defined it from iSpn
z = 1
do i = 1, nBasis / 2
IF (iSpn == 1) THEN
!We want to run through all beta orbitals (odd numbered basis function)
OrbA = (2 * i) - 1
ELSE
!We want to run through all alpha orbitals (odd numbered basis functions)
OrbA = (2 * i)
end if
!We need to find if allowed
IF (.not. BTEST(ILUT((OrbA - 1) / bits_n_int), MOD((OrbA - 1), bits_n_int))) THEN
!Is not in the original determinant - allow
!Now check that its symmetry allowed
IF (IsAOrbSymAllowed(iSpn, OrbA, SpinOrbA, SymProduct, SumMl, SymA, SymB, MlA, MlB, ClassCountUnocc2)) THEN
IF (z == ChosenUnocc) THEN
!We have found our allowed orbital
!OrbA is the allowed orbital
RETURN
end if
z = z + 1
end if
end if
end do
CALL Stop_All("PickAOrb", "Should not get here - have not found allowed A Orb")
ELSE
! METHOD 2 (Keep drawing orbitals randomly until we find one unoccupied). This should be more
!efficient, unless we have small basis sets.
! =========================
Attempts = 0
do while (.true.)
!Draw randomly from the set of orbitals
r = genrand_real2_dSFMT()
ChosenUnocc = INT((nBasis / 2) * r) + 1
IF (iSpn == 1) THEN
!We only want to choose beta orbitals(i.e. odd).
ChosenUnocc = (2 * ChosenUnocc) - 1
ELSE
!We only want to choose alpha orbitals(i.e. even).
ChosenUnocc = 2 * ChosenUnocc
end if
!Find out if orbital is in nI or not. Accept if it isn't in it.
IF (.not. (BTEST(ILUT((ChosenUnocc - 1) / bits_n_int), MOD((ChosenUnocc - 1), bits_n_int)))) THEN
!Orbital not in nI. Accept.
EXIT
end if
IF (Attempts > 250) THEN
write(stdout, *) "Cannot find A unoccupied orbital after 250 attempts..."
call write_det(stdout, nI, .TRUE.)
CALL Stop_All("PickAOrb", "Cannot find A unoccupied orbital after 250 attempts...")
end if
Attempts = Attempts + 1
end do
OrbA = ChosenUnocc !This is the allowed orbital
end if
end if
!We now need to test whether this orbital has any symmetry-allowed unoccupied orbitals to form a pair with.
!To test this, we need to find the needed symmetry of B, in order that Sym(A) x Sym(B) = SymProduct
IF (IsAOrbSymAllowed(iSpn, OrbA, SpinOrbA, SymProduct, SumMl, SymA, SymB, MlA, MlB, ClassCountUnocc2)) THEN
!OrbA picked is allowed, exit from loop
EXIT
end if
IF (AttemptsOverall > (nBasis * 20)) THEN
write(stdout, *) "Cannot find first allowed unoccupied orbital for given i,j pair after ", nBasis * 20, " attempts."
write(stdout, *) "It may be that there are no possible excitations from this i,j pair, in which case "
write(stdout, *) "the given algorithm is inadequate to describe excitations from such a small space."
write(stdout, *) "Try reverting to old excitation generators."
call write_det(stdout, nI, .TRUE.)
write(stdout, *) "***", NExcit, ForbiddenOrbs
write(stdout, *) "I,J pair; sym_i, sym_j: ", nI(Elec1Ind), nI(Elec2Ind), G1(nI(Elec1Ind))%Sym%S, G1(nI(Elec2Ind))%Sym%S
CALL Stop_All("PickAOrb", "Cannot find first allowed unocc orb for double excitation")
end if
AttemptsOverall = AttemptsOverall + 1
end do
END SUBROUTINE PickAOrb