SUBROUTINE PickBOrb(nI, iSpn, ILUT, ClassCountUnocc2, SpinOrbA, OrbA, SymA, OrbB, SymB, NExcit, MlA, MlB, NExcitOtherWay)
integer, intent(in) :: nI(nel), iSpn, SpinOrbA, OrbA, SymA, SymB
integer, intent(in) :: MlA, MlB
integer, intent(in) :: ClassCountUnocc2(ScratchSize)
integer, intent(out) :: nExcitOtherWay, nExcit, OrbB
integer(n_int), intent(in) :: iLut(0:NIfTot)
integer :: norbs, i, z, ind, ChosenUnocc, attempts, SpinOrbB
real(dp) :: r
!We want to calculate the number of possible B's given the symmetry
! and spin it has to be since we have already picked A.
!We have calculated in NExcit the number of orbitals available for B given A,
! but we also need to know the number of orbitals to choose from for A IF
!we had picked B first.
ind = 0
IF (iSpn == 2) THEN
!If iSpn=2, then we want to find a spinorbital of the opposite spin of SpinOrbA
IF (SpinOrbA == -1) THEN
!We have already picked a beta orbital, so now we want to pick an alpha orbital.
! Find out how many of these there are.
NExcit=ClassCountUnocc2(ClassCountInd(1,SymB,MlB))
NExcitOtherWay=ClassCountUnocc2(ClassCountInd(2,SymA,MlA))
SpinOrbB=1 !This is defined differently to SpinOrbA. 1=Alpha, 2=Beta.
ELSE
!Want to pick an beta orbital.
NExcit = ClassCountUnocc2(ClassCountInd(2, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(1, SymA, MlA))
SpinOrbB = 2
end if
else if (iSpn == 1) THEN
!Definitely want a beta orbital
NExcit = ClassCountUnocc2(ClassCountInd(2, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(2, SymA, MlA))
SpinOrbB = 2
ELSE
!Definitely want an alpha orbital
NExcit = ClassCountUnocc2(ClassCountInd(1, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(1, SymA, MlA))
SpinOrbB = 1
end if
IF ((iSpn /= 2) .and. (SymA == SymB) .and. (MlA == MlB)) THEN
!In this case, we need to check that we do not pick the same orbital as OrbA. If we do this, then we need to redraw.
!Only when SymProduct=0 will the classes of a and b be the same, and the spins will be different if iSpn=2,
! so this is the only possibility of a clash.
NExcit = NExcit - 1 !Subtract 1 from the number of possible orbitals since we cannot choose orbital A.
NExcitOtherWay = NExcitOtherWay - 1 !The same goes for the probabilities the other way round.
end if
!All orbitals with the specified symmetry and spin should be allowed unless it is OrbA.
!There will be NExcit of these. Pick one at random.
!Check that orbital is not in ILUT and is not = OrbA (Although this can only happen
!in the circumstance indicated earlier).
!Now we need to choose the final unoccupied orbital.
!There are two ways to do this.
! We can either choose the orbital we want out of the NExcit possible unoccupied orbitals.
!It would then be necessary to cycle through all orbitals of that symmetry and spin,
! only counting the unoccupied ones to find the correct determinant.
!Alternatively, we could pick orbitals at random and redraw until we find an allowed one.
! This would probably be preferable for larger systems.
IF (tCycleOrbs) THEN
! METHOD 1 (Run though all orbitals in symmetry class with needed spin to find allowed one out of NExcit)
! ==========================
!Choose the unoccupied orbital to exite to
r = genrand_real2_dSFMT()
ChosenUnocc = INT(NExcit * r) + 1
!Now run through all allowed orbitals until we find this one.
IF (tNoSymGenRandExcits) THEN
nOrbs = nBasis / 2 !No symmetry, therefore all orbitals of allowed spin possible to generate.
ELSE
Ind = ClassCountInd(SpinOrbB, SymB, MlB)
nOrbs = SymLabelCounts2(2, Ind)
end if
z = 0 !z is the counter for the number of allowed unoccupied orbitals we have gone through
do i = 0, nOrbs - 1
IF (tNoSymGenRandExcits) THEN
OrbB = (2 * (i + 1)) - (SpinOrbB - 1)
ELSE
!Find the spin orbital index. SymLabelCounts has the index of the state for the given symmetry.
OrbB = SymLabelList2(SymLabelCounts2(1, Ind) + i)
ENDIF
!Find out if the orbital is in the determinant, or is the other unocc picked
IF ((.not. (BTEST(ILUT((OrbB - 1) / bits_n_int), MOD((OrbB - 1), bits_n_int)))) .and. (OrbB /= OrbA)) THEN
!The orbital is not found in the original determinant - increment z
z = z + 1
IF (z == ChosenUnocc) THEN
!We have got to the determinant that we want to pick.
EXIT
end if
end if
end do
!We now have our final orbitals.
IF (z /= ChosenUnocc) THEN
write(stdout, *) "This is a problem, since there should definitely be an allowed beta orbital once alpha is chosen..."
CALL Stop_All("PickBOrb", "Could not find allowed unoccupied orbital to excite to.")
end if
ELSE
! METHOD 2 (Keep drawing orbitals from the desired symmetry and spin until we find one unoccupied)
! =========================
IF (tNoSymGenRandExcits) THEN
nOrbs = nBasis / 2
ELSE
Ind = ClassCountInd(SpinOrbB, SymB, MlB)
nOrbs = SymLabelCounts2(2, Ind)
end if
Attempts = 0
do while (.true.)
!Draw randomly from the set of orbitals
r = genrand_real2_dSFMT()
ChosenUnocc = INT(nOrbs * r)
IF (tNoSymGenRandExcits) THEN
OrbB = (2 * (ChosenUnocc + 1)) - (SpinOrbB - 1)
ELSE
OrbB = SymLabelList2(SymLabelCounts2(1, Ind) + ChosenUnocc)
ENDIF
!Find out if orbital is in nI or not. Accept if it isn't in it.
IF ((.not. (BTEST(ILUT((OrbB - 1) / bits_n_int), MOD((OrbB - 1), bits_n_int)))) .and. (OrbB /= OrbA)) THEN
!Orbital not in nI. Accept.
EXIT
end if
IF (Attempts > 1000) THEN
write(stdout, *) "Cannot find double excitation unoccupied orbital after 1000 attempts..."
write(stdout, *) "This is a problem, since there should definitly be an allowed beta orbital once alpha is chosen..."
write(stdout, *) "nI: "
call write_det(stdout, nI, .TRUE.)
write(stdout, *) "iSpn: ", iSpn
write(stdout, *) "ClassCountUnocc2: ", ClassCountUnocc2(:)
write(stdout, *) "NExcit", NExcit
CALL Stop_All("PickBOrb", "Cannot find double excitation unoccupied orbital after 250 attempts...")
end if
Attempts = Attempts + 1
end do
end if
END SUBROUTINE PickBOrb