SUBROUTINE CreateSingleExcitBiased(nI, nJ, iLut, ExcitMat, tParity, ElecsWNoExcits, nParts, WSign, Tau, iCreate)
INTEGER :: ClassCount2(ScratchSize), i, Attempts, OrbA
INTEGER :: ClassCountUnocc2(ScratchSize), Ind
INTEGER :: ElecsWNoExcits, nParts, iCreate, nI(NEl), nJ(NEl)
REAL(dp) :: WSign
INTEGER(KIND=n_int) :: iLut(0:NIfTot)
INTEGER :: ExcitMat(2, maxExcit), SpawnOrb(nBasis), Eleci, ElecSym, NExcit, VecInd, ispn, EndSymState, j
real(dp) :: Tau, SpawnProb(nBasis), NormProb, r, rat
LOGICAL :: tParity
HElement_t(dp) :: rh
debug_function_name("CreateSingleExcitBiased")
!First, we need to do an O[N] operation to find the number of occupied alpha electrons, number of occupied beta electrons
!and number of occupied electrons of each symmetry class and spin. This is similar to the ClassCount array.
!This has the format (Spn,sym), where Spin=1,2 corresponding to alpha and beta.
!For molecular systems, sym runs from 0 to 7. This is NOT general and should be made so using SymLabels.
!This could be stored to save doing this multiple times, but shouldn't be too costly an operation.
CALL construct_class_counts(nI, ClassCount2, ClassCountUnocc2)
!We need to find out if there are any electrons which have no possible excitations.
!This is because these will need to be redrawn and so will affect the probabilities.
ElecsWNoExcits = 0
!Need to look for forbidden electrons through all the irreps.
do i = 0, nSymLabels - 1
!Run through all labels
IF ((ClassCount2(ClassCountInd(1, i, 0)) /= 0) .and. (ClassCountUnocc2(ClassCountInd(1, i, 0)) == 0)) THEN
!If there are alpha electrons in this class with no possible unoccupied alpha orbitals
!in the same class, these alpha electrons have no single excitations.
ElecsWNoExcits = ElecsWNoExcits + ClassCount2(ClassCountInd(1, i, 0))
end if
IF ((ClassCount2(ClassCountInd(2, i, 0)) /= 0) .and. (ClassCountUnocc2(ClassCountInd(2, i, 0)) == 0)) THEN
ElecsWNoExcits = ElecsWNoExcits + ClassCount2(ClassCountInd(2, i, 0))
end if
end do
IF (ElecsWNoExcits == NEl) THEN
!There are no single excitations from this determinant at all.
!Then we will create a double excitation instead.
RETURN
end if
!We want to pick an occupied electron - Prob = 1/(N-ElecsWNoExcits)
Attempts = 0
do while (.true.)
Attempts = Attempts + 1
!Choose an electron randomly...
r = genrand_real2_dSFMT()
Eleci = INT(NEl * r) + 1
!Find symmetry of chosen electron
ElecSym = INT((G1(nI(Eleci))%Sym%S), 4)
IF (G1(nI(Eleci))%Ms == 1) THEN
!Alpha orbital - see how many single excitations there are from this electron...
NExcit = ClassCountUnocc2(ClassCountInd(1, ElecSym, 0))
ispn = 1 !Alpha ispn=1, Beta ispn=2
ELSE
!Beta orbital
NExcit = ClassCountUnocc2(ClassCountInd(2, ElecSym, 0))
ispn = 2
end if
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("CreateSingleExcit", "Cannot find single excitation from electrons after 250 attempts...")
end if
end do
ExcitMat(1, 1) = nI(Eleci)
!Now we want to run through all sym+spin allowed excitations of the chosen electron, and determine their matrix elements
!To run just through the states of the required symmetry we want to use SymLabelCounts.
!We also want to take into account spin. We want the spin of the chosen unoccupied
!orbital to be the same as the chosen occupied orbital.
!Run over all possible a orbitals
Ind = ClassCountInd(ispn, ElecSym, 0)
EndSymState = SymLabelCounts2(1, Ind) + SymLabelCounts2(2, Ind) - 1
VecInd = 1
NormProb = 0.0_dp
do j = SymLabelCounts2(1, Ind), EndSymState
!We want to look through all beta orbitals
OrbA = SymLabelList2(j) !This is the spin orbital chosen for a
IF (BTEST(ILUT((OrbA - 1) / bits_n_int), MOD((OrbA - 1), bits_n_int))) THEN
!Orbital is in nI...not an unoccupied orbital
CYCLE
end if
!Now we want to find the information about this excitation
ExcitMat(2, 1) = OrbA
rh = sltcnd_excit(nI, Excite_1_t(ExcitMat(:, :1)), .false.)
SpawnProb(VecInd) = abs(REAL(rh, dp))
SpawnOrb(VecInd) = OrbA
NormProb = NormProb + SpawnProb(VecInd)
VecInd = VecInd + 1
IF (VecInd >= nBasis) THEN
CALL Stop_All("CreateSingleExcitBiased", "Finding too many virtual pairs...")
end if
end do
IF ((VecInd - 1) /= NExcit) THEN
CALL Stop_All("CreateSingleExcitBiased", "Wrong number of excitations found from chosen electron")
end if
IF (VecInd == 1) THEN
CALL Stop_All("CreateSingleExcitBiased", "No excitations found from electron, but some should exist")
end if
!We now have to find out how many children to spawn, based on the value of normprob.
rat = Tau * NormProb * REAL((NEl - ElecsWNoExcits) * nParts, dp) / (1.0_dp - PDoubNew)
iCreate = INT(rat)
rat = rat - REAL(iCreate, dp)
r = genrand_real2_dSFMT()
IF (rat > r) THEN
!Child is created
iCreate = iCreate + 1
end if
if (tGUGA) then
call stop_all("CreateSingleExcitBiased", &
"modify get_helement for GUGA")
end if
IF (iCreate > 0) THEN
!We want to spawn particles. This only question now is where. Run
!through the ab pairs again and choose based on the SpawnProb element.
r = genrand_real2_dSFMT()
r = r * NormProb
i = 0
do while (r > 0.0_dp)
i = i + 1
r = r - SpawnProb(i)
end do
IF (i > VecInd - 1) THEN
CALL Stop_All("CreateSingleExcitBiased", "Chosen virtual does not correspond to allowed orbital")
end if
OrbA = SpawnOrb(i)
! Construct the new determinant, excitation matrix and parity
call make_single(nI, nJ, eleci, orbA, ExcitMat, tParity)
!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.
ASSERT(SymAllowedExcit(nI, nJ, 1, ExcitMat))
!Once we have the definitive determinant, we also want to find out what sign the particles we want to create are.
!iCreate is initially positive, so its sign can change depending on the sign of the connection and of the parent particle(s)
rh = get_helement(nI, nJ, 1, ExcitMat, tParity)
IF (WSign > 0.0_dp) THEN
!Parent particle is positive
IF (real(rh, dp) > 0.0_dp) THEN
iCreate = -iCreate !-ve walker created
end if
ELSE
IF (real(rh, dp) < 0.0_dp) THEN
iCreate = -iCreate !-ve walkers created
end if
end if
end if
END SUBROUTINE CreateSingleExcitBiased