SUBROUTINE CreateSingleExcit(nI, nJ, ClassCount2, ClassCountUnocc2, ILUT, ExcitMat, tParity, pGen)
INTEGER :: ElecsWNoExcits, i, Attempts, nOrbs, z, Orb
INTEGER :: Eleci, ElecSym, nI(NEl), nJ(NEl), NExcit, iSpn, ChosenUnocc
INTEGER :: ExcitMat(2, maxExcit)
INTEGER :: ClassCount2(ScratchSize)
INTEGER :: ClassCountUnocc2(ScratchSize), ElecK, SymIndex
INTEGER(KIND=n_int) :: ILUT(0:NIfTot)
real(dp) :: r, pGen
LOGICAL :: tParity
character(*), parameter :: t_r = 'CreateSingleExcit'
CALL CheckIfSingleExcits(ElecsWNoExcits, ClassCount2, ClassCountUnocc2, nI)
IF (ElecsWNoExcits == NEl) THEN
!There are no single excitations from this determinant - return a null excitation
nJ(1) = 0
pgen = 0.0_dp
RETURN
end if
Attempts = 0
elecK = 0
orb = 0
do while (.true.)
!Choose an electron randomly...
r = genrand_real2_dSFMT()
Eleci = INT(NEl * r) + 1
!Find symmetry of chosen electron
IF (tNoSymGenRandExcits) THEN
ElecSym = 0
ELSE
!For abelian symmetry, the irrep of i and a must be the same.
!For solids, this means that the excitation must be within the same k-point
ElecSym = SpinOrbSymLabel(nI(Eleci))
ElecK = G1(nI(Eleci))%Ml
end if
IF (G1(nI(Eleci))%Ms == 1) THEN
!Alpha orbital - see how many single excitations there are from this electron...
iSpn = 1
ELSE
!Beta orbital
iSpn = 2
end if
SymIndex = ClassCountInd(iSpn, ElecSym, ElecK)
NExcit = ClassCountUnocc2(SymIndex)
IF (NExcit /= 0) EXIT !Have found electron with allowed excitations
IF (Attempts > 250) THEN
write(stdout, *) "Cannot find single excitation from electrons after 250 attempts..."
call write_det(stdout, nI, .true.)
write(stdout, *) "***"
call stop_all(t_r, "Cannot find single excitation from &
&electrons after 250 attempts...")
end if
Attempts = Attempts + 1
end do
!Now we need to choose the unoccupied orbital for the chosen electron.
!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
ELSE
nOrbs = OrbClassCount(SymIndex)
end if
z = 0 !z is the counter for the number of allowed unoccupied orbitals we have gone through
do i = 0, nOrbs - 1
!Find the spin orbital index. SymLabelCounts has the index of the state for the given symmetry.
IF (tNoSymGenRandExcits) THEN
Orb = (2 * (i + 1)) - (iSpn - 1)
ELSE
Orb = SymLabelList2(SymLabelCounts2(1, SymIndex) + i)
end if
!Find out if the orbital is in the determinant.
IF (.not. (BTEST(ILUT((Orb - 1) / bits_n_int), MOD(Orb - 1, bits_n_int)))) 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. i=nI(Eleci). a=Orb.
IF (z /= ChosenUnocc) THEN
call stop_all(t_r, "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
nOrbs = OrbClassCount(SymIndex)
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
Orb = (2 * (ChosenUnocc + 1)) - (iSpn - 1)
ELSE
Orb = SymLabelList2(SymLabelCounts2(1, SymIndex) + ChosenUnocc)
end if
!Find out if orbital is in nI or not. Accept if it isn't in it.
IF (.not. (BTEST(ILUT((Orb - 1) / bits_n_int), MOD(Orb - 1, bits_n_int)))) THEN
!Orbital not in nI. Accept.
EXIT
end if
IF (Attempts > 250) THEN
write(stdout, *) "Cannot find single excitation unoccupied orbital after 250 attempts..."
write(stdout, *) "Desired symmetry of unoccupied orbital = ", ElecSym
write(stdout, *) "Number of orbitals (of correct spin) in symmetry = ", nOrbs
write(stdout, *) "Number of orbitals to legitimatly pick = ", NExcit
call write_det(stdout, nI, .true.)
call stop_all(t_r, "Cannot find single excitation &
&unoccupied orbital after 250 attempts...")
end if
Attempts = Attempts + 1
end do
end if
! Construct the new determinant, excitation matrix and parity
call make_single(nI, nJ, eleci, orb, ExcitMat, tParity)
#ifdef DEBUG_
! These are useful (but O[N]) operations to test the determinant
! generated. If there are any problems with then excitations, I
! recommend uncommenting these tests to check the results.
if (.not. SymAllowedExcit(nI, nJ, 1, ExcitMat)) &
call stop_all(t_r, 'Generated excitation invalid')
#endif
!Now we need to find the probability of creating this excitation.
!This is: P_single x P(i) x P(a|i) x N/(N-ElecsWNoExcits)
pGen = (1 - pDoubNew) / (REAL((NExcit * (NEl - ElecsWNoExcits)), dp))
END SUBROUTINE CreateSingleExcit