CreateSingleExcitBiased Subroutine

public subroutine CreateSingleExcitBiased(nI, nJ, iLut, ExcitMat, tParity, ElecsWNoExcits, nParts, WSign, Tau, iCreate)


Type IntentOptional Attributes Name
integer :: nI(NEl)
integer :: nJ(NEl)
integer(kind=n_int) :: iLut(0:NIfTot)
integer :: ExcitMat(2,maxExcit)
logical :: tParity
integer :: ElecsWNoExcits
integer :: nParts
real(kind=dp) :: WSign
real(kind=dp) :: Tau
integer :: iCreate


Source Code

    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

!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.
        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
!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
            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
                IF (real(rh, dp) < 0.0_dp) THEN
                    iCreate = -iCreate    !-ve walkers created
                end if
            end if

        end if

    END SUBROUTINE CreateSingleExcitBiased