#include "macros.h"
MODULE GenRandSymExcitNUMod
!This is a new random excitation generator. It still creates excitations with a normalised and calculable probability,
!but these probabilities are not uniform. They should also be quicker to generate, and have small or no excitation generators
!to store. nI is the root determinant and nJ is the excitation.
!This currently is only available for abelian symmetry. This is because I label the symmetries as integers, but this
!algorithm could theoretically be applied relatively easily to non-abelian symmetry. The symmetries would have to be stored
!in a better way. For molecular systems, symmetry irrep is just symmetry_label-1. Symlabels will need to be used for more
!complicated symmetry tables.
!THEORY
! P[i,j,a,b] = P_doub x ( P[i,j] x P[a,b|i,j] ) x (M-N)/(M-N-Q)
! P[a,b|i,j] = P[b|i,j] x P[a|b,i,j] + P[a|i,j] x P[b|a,i,j]
! P[i,j] = 2/N(N-1) always - others are more difficult to calculate
! Difficulties arise if there is ever no allowed final unoccupied orbital to pick.
! Spin and spatial symmetry needs to be considered.
! The fourth orbital chosen (i.e. one of the unoccupied ones) will have its
! symmetry class specified by the classes of the other ones.
! i and j are picked randomly. They can be alpha,beta | beta,alpha | alpha,alpha | beta,beta
! a is then picked from all orbitals with the allowed spin. This means that a beta,beta i,j pair
! will force a to also be beta. There are also a number of forbidden a orbitals which are counted.
! These are forbidden since they have no possible b orbital which will give rise to a symmetry and
! spin allowed unoccupied a,b pair. The number of these orbitals, Q, is needed to calculate the
! normalised probability of generating the excitation.
use SystemData, only: ALAT, iSpinSkip, tFixLz, iMaxLz, tUEG, tNoFailAb, &
tLatticeGens, tHub, nEl, G1, nBasis, nBasisMax, &
tNoSymGenRandExcits, Arr, nMax, tCycleOrbs, &
nOccAlpha, nOccBeta, ElecPairs, MaxABPairs, &
tKPntSym, lzTot, tNoBrillouin, tUseBrillouin, &
tUEG2, kvec, tAllSymSectors, NMAXX, NMAXY, NMAXZ, &
tOrbECutoff, OrbECutoff, tGUGA, t_pcpp_excitgen, &
t_new_real_space_hubbard
use CalcData, only: tSpinProject
use FciMCData, only: pDoubles, psingles, iter, excit_gen_store_type, iluthf
use Parallel_neci
use IntegralsData, only: UMat
use Determinants, only: get_helement
use DeterminantData, only: write_det, FDet
use SymData, only: nSymLabels, TwoCycleSymGens, SymLabelList, &
SymLabelCounts
use dSFMT_interface, only: genrand_real2_dSFMT
use SymExcitDataMod
use DetBitOps, only: FindExcitBitDet, EncodeBitDet, detbiteq
use excitation_types, only: Excite_1_t
use sltcnd_mod, only: sltcnd_excit
use constants, only: dp, n_int, bits_n_int, maxExcit, stdout
use bit_rep_data, only: NIfTot
use sym_mod, only: mompbcsym, GetLz
use timing_neci
use sym_general_mod
use get_excit, only: make_single, make_double
use procedure_pointers, only: get_umat_el
use back_spawn, only: pick_virtual_electrons_double_hubbard, &
pick_occupied_orbital_hubbard, check_electron_location
use lattice_models_utils, only: get_ispn
use util_mod, only: stop_all
use neci_intfce
use bit_rep_data, only: test_flag
use sym_mod
use SystemData, only: G1, NEl, Symmetry, nBasisMax, BasisFN
use FciMCData, only: HFSym
implicit none
contains
subroutine gen_rand_excit(nI, ilut, nJ, ilutnJ, exFlag, IC, ExcitMat, &
tParity, pGen, HElGen, store, part_type)
! This routine is the same as GenRandSymexcitNu, but you can pass in
! the class count arrays so that they do not have to be recalculated
! each time for the same excitation. If tFilled is false, it will
! assume they are unfilled and calculate them, returning the arrays
! and tFilled=T.
! The two arrays want to be integers, both of size (1, 1:nSymLabels)
integer, intent(in) :: nI(nel), exFlag
integer(n_int), intent(in) :: iLut(0:niftot)
integer, intent(out) :: nJ(nel), IC, ExcitMat(2, maxExcit)
logical, intent(out) :: tParity
real(dp), intent(out) :: pgen
type(excit_gen_store_type), intent(inout), target :: store
! Not used
integer(n_int), intent(out) :: ilutnJ(0:niftot)
HElement_t(dp), intent(out) :: HElGen
integer, intent(in), optional :: part_type
real(dp) :: r
character(*), parameter :: this_routine = 'gen_rand_excit'
! Just in case
ilutnJ(0) = -1
HElGen = 0.0_dp
IF ((tUEG .and. tLatticeGens) .or. (tHub .and. tLatticeGens)) THEN
call CreateExcitLattice(nI, iLut, nJ, tParity, ExcitMat, pGen, part_type)
IC = 2
RETURN
end if
!TODO: Not quite sure what conditions we need to check for now...
if (.not. store%tFilled) then
!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 stored to save doing this multiple times, but shouldn't be too costly an operation.
CALL construct_class_counts(nI, store%ClassCountOcc, &
store%ClassCountUnocc)
store%tFilled = .true.
end if
!ExFlag is 1 for singles, 2 for just doubles, and 3 for both.
IF (ExFlag == 3) THEN
!Choose whether to generate a double or single excitation. Prob of generating a double is given by pDoub.
pDoubNew = pDoubles
r = genrand_real2_dSFMT()
IF (r < pDoubNew) THEN
!A double excitation has been chosen to be created.
IC = 2
ELSE
IC = 1
end if
else if (ExFlag == 2) THEN
IC = 2
pDoubNew = 1.0_dp
else if (ExFlag == 1) THEN
IC = 1
pDoubNew = 0.0_dp
ELSE
CALL Stop_All(this_routine, "Error in choosing excitations to create.")
end if
IF (IC == 2) THEN
CALL CreateDoubExcit(nI, nJ, store%ClassCountUnocc, ILUT, &
ExcitMat, tParity, pGen)
ELSE
CALL CreateSingleExcit(nI, nJ, store%ClassCountOcc, &
store%ClassCountUnocc, ILUT, ExcitMat, &
tParity, pGen)
end if
end subroutine
SUBROUTINE CreateDoubExcit(nI, nJ, ClassCountUnocc2, ILUT, ExcitMat, tParity, pGen)
integer, intent(in) :: nI(nel)
integer, intent(out) :: nJ(nel), ExcitMat(2, maxExcit)
integer, intent(in) :: ClassCountUnocc2(ScratchSize)
integer(n_int), intent(in) :: iLut(0:NIfTot)
real(dp), intent(out) :: pGen
logical, intent(out) :: tParity
integer :: nExcitOtherWay, orbB, nExcitB, SpinOrbA, OrbA, SymA, SymB
integer :: nExcitA, sumMl, mlA, mlB, iSpn, Elec1Ind, Elec2Ind
integer :: SymProduct, ForbiddenOrbs
logical :: tAOrbFail
!First, we need to pick an unbiased distinct electron pair.
!These have symmetry product SymProduct, and spin pair iSpn = 1=beta/beta; 2=alpha/beta; 3=alpha/alpha
CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, -1)
!This routine finds the number of orbitals which are allowed by spin,
! but not part of any spatial symmetry allowed unoccupied pairs.
!This number is needed for the correct normalisation of the probability of drawing
! any given A orbital since these can be chucked and redrawn.
IF (tNoSymGenRandExcits) THEN
CALL FindNumForbiddenOrbsNoSym(ForbiddenOrbs, ClassCountUnocc2, iSpn)
ELSE
CALL FindNumForbiddenOrbs(ForbiddenOrbs, ClassCountUnocc2, SymProduct, iSpn, SumMl)
end if
!Now we have to pick the first unoccupied orbital. If an orbital is not present in any allowed pairs,
! it is chucked and a new one drawn.
!The number NExcit is the number of unoccupied orbitals that the orbital was chosen from
! (including symmetry-forbidden orbital pairs) Arguments:
CALL PickAOrb(nI,iSpn,ILUT,ClassCountUnocc2,NExcitA,Elec1Ind,Elec2Ind,SpinOrbA,OrbA,SymA,SymB, &
& SymProduct,SumMl,MlA,MlB,ForbiddenOrbs,tAOrbFail)
IF(tAOrbFail) THEN
nJ(1)=0
pGen=HUGE(0.0_dp)
RETURN
end if
!This routine will pick an unoccupied orbital at random from a specified spin and symmetry class.
!There should definitely be a possible spinorbital, since A was chosen so that there was one.
!We have to make sure with alpha/alpha or beta/beta pairs and when SymProduct=0,
! that we don't choose the same unoccupied orbital.
!If we do this, then we should chuck and redraw,
! since there should definitely be another allowed spinorbital in the class.
!We return the number of allowed B's for the A we picked in NExcitB,
! however we also need to know the number of allowed A's if we
!had picked B first. This will be returned in NExcitOtherWay.
CALL PickBOrb(nI, iSpn, ILUT, ClassCountUnocc2, SpinOrbA, OrbA, SymA, OrbB, SymB, NExcitB, MlA, MlB, NExcitOtherWay)
call make_double(nI, nJ, elec1ind, elec2ind, orbA, orbB, &
ExcitMat, tParity)
CALL FindDoubleProb(ForbiddenOrbs, NExcitA, NExcitB, NExcitOtherWay, pGen)
END SUBROUTINE CreateDoubExcit
!This routine finds the probability of creating the excitation.
! See the header of the file for more information on how this works.
SUBROUTINE FindDoubleProb(ForbiddenOrbs, NExcitA, NExcitB, NExcitOtherWay, pGen)
INTEGER, INTENT(IN) :: ForbiddenOrbs, NExcitA, NExcitB, NExcitOtherWay
real(dp), INTENT(OUT) :: pGen!,PabGivenij
pGen = pDoubles * ((1.0_dp / real(NExcitB, dp)) + (1.0_dp / real(NExcitOtherWay, dp))) / (REAL((ElecPairs * (NExcitA - ForbiddenOrbs)), dp))
END SUBROUTINE FindDoubleProb
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
!This routine does the same as the FindNumForbiddenOrbs routine,
! but is optimised for when there are no spatial symmetry considerations.
SUBROUTINE FindNumForbiddenOrbsNoSym(ForbiddenOrbs, ClassCountUnocc2, iSpn)
integer, intent(in) :: ClassCountUnocc2(ScratchSize), iSpn
integer, intent(out) :: ForbiddenOrbs
!We know that all orbitals are totally symmetric, and that the symproduct=0
ForbiddenOrbs = 0
IF (iSpn == 2) THEN
IF (ClassCountUnocc2(1) == 0) THEN
!There are no unoccupied alpha orbitals - are there any beta spins which are now forbidden?
ForbiddenOrbs = ClassCountUnocc2(2)
end if
IF (ClassCountUnocc2(2) == 0) THEN
!There are no unoccupied beta orbitals - are there any alpha spins which are now forbidden?
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(1)
end if
else if (iSpn == 1) THEN
!If the symmetry product of the occupied orbitals is 0, then the a,b pair want to be taken from the same class.
!This means that if there is only one spin-allowed orbital in that class, it has no symmetry-allowed pairs, and so is forbidden.
IF (ClassCountUnocc2(2) == 1) THEN
!The one beta orbital in this class is forbidden, since it cannot form a pair.
ForbiddenOrbs = 1
end if
else if (iSpn == 3) THEN
IF (ClassCountUnocc2(1) == 1) THEN
ForbiddenOrbs = 1
end if
end if
END SUBROUTINE FindNumForbiddenOrbsNoSym
!This routine finds the number of orbitals which are allowed by spin,
! but not part of any spatial symmetry allowed unoccupied pairs.
!This number is needed for the correct normalisation of the probability of
! drawing any given A orbital since these can be chucked and redrawn.
!For Lz symmetry, it is generally quicker to count the allowed orbitals, and subtract from all possible ones,
! rather than directly counting
!the forbidden ones.
SUBROUTINE FindNumForbiddenOrbs(ForbiddenOrbs, ClassCountUnocc2, SymProduct, iSpn, SumMl)
integer, intent(in) :: ClassCountUnocc2(ScratchSize)
integer, intent(in) :: SumMl, iSpn, SymProduct
integer, intent(out) :: ForbiddenOrbs
integer :: i, k, Ind, AllowedOrbs, SymInd, OrbAMl, SymOrbs, ConjSym
ForbiddenOrbs = 0
IF (tFixLz) THEN
AllowedOrbs = 0
IF (iSpn == 2) THEN
Ind = 1
do k = -iMaxLz, iMaxLz
OrbAMl = SumMl - k
IF ((k <= OrbAMl) .and. (abs(OrbAMl) <= iMaxLz)) THEN
do i = 0, nSymLabels - 1
SymInd = ClassCountInd(1, IEOR(SymProduct, i), OrbAMl) !Alpha of the corresponding a orbital
IF (ClassCountUnocc2(Ind) /= 0) THEN
!Check the beta conjugate orbital
SymOrbs = ClassCountUnocc2(SymInd + 1)
IF (SymOrbs /= 0) THEN
IF (k == OrbAMl) THEN
!The Ml symmetries are the same! Don't double count.
AllowedOrbs = AllowedOrbs + SymOrbs
ELSE
AllowedOrbs = AllowedOrbs + SymOrbs + ClassCountUnocc2(Ind)
end if
end if
end if
IF (ClassCountUnocc2(Ind + 1) /= 0) THEN
SymOrbs = ClassCountUnocc2(SymInd)
IF (SymOrbs /= 0) THEN
IF (k == OrbAMl) THEN
!The Ml symmetries are the same! Don't double count.
AllowedOrbs = AllowedOrbs + SymOrbs
ELSE
AllowedOrbs = AllowedOrbs + SymOrbs + ClassCountUnocc2(Ind + 1)
end if
end if
end if
Ind = Ind + 2
end do
ELSE
!Move onto the next k-block of B orbitals.
Ind = Ind + nSymLabels * 2
end if
end do
ForbiddenOrbs = nBasis - NEl - AllowedOrbs
else if (iSpn == 3) THEN !alpha/alpha - run through all alpha orbitals
Ind = 1
do k = -iMaxLz, iMaxLz
OrbAMl = SumMl - k
IF ((k <= OrbAMl) .and. (abs(OrbAMl) <= iMaxLz)) THEN
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) /= 0) THEN
IF ((SymProduct == 0) .and. (OrbAMl == k)) THEN
IF (ClassCountUnocc2(Ind) > 1) THEN
AllowedOrbs = AllowedOrbs + ClassCountUnocc2(Ind)
end if
ELSE
SymOrbs = ClassCountUnocc2(ClassCountInd(1, IEOR(SymProduct, i), OrbAMl))
IF (SymOrbs /= 0) THEN
IF (k == OrbAMl) THEN
!The Ml symmetries are the same! Don't double count.
AllowedOrbs = AllowedOrbs + SymOrbs
ELSE
AllowedOrbs = AllowedOrbs + SymOrbs + ClassCountUnocc2(Ind)
end if
end if
end if
end if
Ind = Ind + 2
end do
ELSE
Ind = Ind + nSymLabels * 2
end if
end do
ForbiddenOrbs = (nBasis / 2) - nOccAlpha - AllowedOrbs
ELSE
Ind = 2
do k = -iMaxLz, iMaxLz
OrbAMl = SumMl - k
IF ((k <= OrbAMl) .and. (abs(OrbAMl) <= iMaxLz)) THEN
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) /= 0) THEN
IF ((SymProduct == 0) .and. (OrbAMl == k)) THEN
IF (ClassCountUnocc2(Ind) > 1) THEN
AllowedOrbs = AllowedOrbs + ClassCountUnocc2(Ind)
end if
ELSE
SymOrbs = ClassCountUnocc2(ClassCountInd(2, IEOR(SymProduct, i), OrbAMl))
IF (SymOrbs /= 0) THEN
IF (k == OrbAMl) THEN
!The Ml symmetries are the same! Don't double count.
AllowedOrbs = AllowedOrbs + SymOrbs
ELSE
AllowedOrbs = AllowedOrbs + SymOrbs + ClassCountUnocc2(Ind)
end if
end if
end if
end if
Ind = Ind + 2
end do
ELSE
Ind = Ind + nSymLabels * 2
end if
end do
ForbiddenOrbs = (nBasis / 2) - nOccBeta - AllowedOrbs
end if
ELSE
!Not Lz symmetry...
IF (iSpn == 2) THEN
! write(stdout,*) "Alpha/Beta"
!i,j are an alpha/beta pair. The number of forbidden orbitals includes all alphas and betas.
Ind = 1
do i = 0, nSymLabels - 1
!Run though all symmetries of possible "a" orbital syms. If there aren't any, then we
!know the corresponding "b" orbitals are excluded.
IF (ClassCountUnocc2(Ind) == 0) THEN
!This symmetry has no unoccupied alpha orbitals - does its symmetry conjugate have any
!unoccupied beta orbitals which are now forbidden?
!If there are no unoccupied orbitals in this conjugate symmetry, then it won't increase
!the forbidden orbital number, since it can never be chosen.
! ConjSym=IEOR(SymProduct,i)
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(2, &
& RandExcitSymLabelProd(SymProduct, SymInvLabel(i)), 0))
!No unocc alphas in i, therefore all betas in ConjSym are forbidden
! write(stdout,*) ClassCountUnocc2(2,ConjSym),i,ConjSym
end if
IF (ClassCountUnocc2(Ind + 1) == 0) THEN
!This symmetry has no unoccupied beta orbitals - does its symmetry conjugate have any
!unoccupied alpha orbitals which are now forbidden?
!If there are no unoccupied orbitals in this conjugate symmetry, then it won't increase
!the forbidden orbital number, since it can never be chosen.
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(1, &
& RandExcitSymLabelProd(SymProduct, SymInvLabel(i)), 0))
end if
Ind = Ind + 2
end do
else if (iSpn == 1) THEN
Ind = 2
IF (.not. tKPntSym) THEN
!With molecular systems, the irreps are their own inverses, so it is a little simpler to do the
!two cases seperately.
IF (SymProduct /= 0) THEN
!i,j are a beta/beta pair. The number of forbidden orbitals is just betas
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) == 0) THEN
! ConjSym=IEOR(SymProduct,i)
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(2, RandExcitSymLabelProd(SymProduct, i), 0))
end if
Ind = Ind + 2
end do
ELSE
!There is a subtle point here, which could change the probabilities.
!If the symmetry product of the occupied orbitals is 0, then the a,b pair want to be taken from the same class.
!This means that if there is only one spin-allowed orbital in that class, it has no symmetry-allowed pairs, and so is forbidden.
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) == 1) THEN
!The one beta orbital in this class is forbidden, since it cannot form a pair.
ForbiddenOrbs = ForbiddenOrbs + 1
end if
Ind = Ind + 2
end do
end if
ELSE
!With KPntSym, we have to work out if we are in the case that sym_a^* = sym_b
!Unfortunately, I don't think you can tell from SymProduct when this case is going to be satisfied.
do i = 0, nSymLabels - 1 !Run over symmetries of the orbitals
IF (ClassCountUnocc2(Ind) <= 1) THEN
ConjSym = RandExcitSymLabelProd(SymProduct, SymInvLabel(i))
IF (ConjSym == i) THEN
!A and B come from the same symmetry, so we must have more than one
!orbitals available from there..
IF (ClassCountUnocc2(Ind) == 1) THEN
!The one beta orbital in this class is forbidden, since it cannot form a pair.
ForbiddenOrbs = ForbiddenOrbs + 1
end if
ELSE
IF (ClassCountUnocc2(Ind) == 0) THEN
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(2, ConjSym, 0))
end if
end if
end if
Ind = Ind + 2
end do
end if
else if (iSpn == 3) THEN
Ind = 1
IF (.not. tKPntSym) THEN
IF (SymProduct /= 0) THEN
!i,j are a alpha/alpha pair. The number of forbidden orbitals is just alphas
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) == 0) THEN
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(1, IEOR(SymProduct, i), 0))
end if
Ind = Ind + 2
end do
ELSE
!There is a subtle point here, which could change the probabilities.
!If the symmetry product of the occupied orbitals is 0, then the a,b pair want to be taken from the same class.
!This means that if there is only one spin-allowed orbital in that class, it has no symmetry-allowed pairs, and so is forbidden.
do i = 0, nSymLabels - 1
IF (ClassCountUnocc2(Ind) == 1) THEN
!The one alpha orbital in this class is forbidden, since it cannot form a pair.
ForbiddenOrbs = ForbiddenOrbs + 1
end if
Ind = Ind + 2
end do
end if
ELSE
!With KPntSym, we have to work out if we are in the case that sym_a^* = sym_b
!Unfortunately, I don't think you can tell from SymProduct when this case is going to be satisfied.
do i = 0, nSymLabels - 1 !Run over symmetries of the orbitals
IF (ClassCountUnocc2(Ind) <= 1) THEN
ConjSym = RandExcitSymLabelProd(SymProduct, SymInvLabel(i))
IF (ConjSym == i) THEN
!A and B come from the same symmetry, so we must have more than one
!orbitals available from there..
IF (ClassCountUnocc2(Ind) == 1) THEN
!The one beta orbital in this class is forbidden, since it cannot form a pair.
ForbiddenOrbs = ForbiddenOrbs + 1
end if
ELSE
IF (ClassCountUnocc2(Ind) == 0) THEN
!There aren't any a's in this pair of syms, so the b's are forbidden
ForbiddenOrbs = ForbiddenOrbs + ClassCountUnocc2(ClassCountInd(1, ConjSym, 0))
end if
end if
end if
Ind = Ind + 2
end do
end if
end if
end if
END SUBROUTINE FindNumForbiddenOrbs
!Choose the first unoccupied orbital to excite to. Spatial symmetry does not need to be considered, but spin symmetry does.
!After one has been chosen, we have to make sure that an allowed symmetry pair can be formed, otherwise the orbital is
!chucked and a new one drawn.
!Arguments: NExcit = Number of possible spin-allowed unoccupied spinorbitals, including forbidden orbs (these will be chucked)
! SpinOrbA = Spin of the chosen spin-orbital. 1 is alpha, -1 is beta.
! OrbA = Index of the chosen spinorbital
! SymB = Symmetry required of the second unoccupied spinorbital, so that sym(A) x sym(B) = SymProduct
! SymProduct = (intent in), the symmetry of electron pair chosen
! SumMl = Sum of the Ml values for the two picked electrons
! MlA = Ml of the a orbital chosen
! MlB = Required Ml of the b orbital still to be chosen
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
!This routine will look at an orbital (OrbA) and check whether it is an allowed A orbital to
!pick, i.e. it has allowed B orbitals, given the i,js.
LOGICAL FUNCTION IsAOrbSymAllowed(iSpn, OrbA, SpinOrbA, SymProduct, SumMl, SymA, SymB, MlA, MlB, ClassCountUnocc2)
INTEGER, INTENT(IN) :: iSpn, OrbA, SpinOrbA, SymProduct, SumMl, ClassCountUnocc2(ScratchSize)
INTEGER, INTENT(OUT) :: SymA, SymB, MlA, MlB
IsAOrbSymAllowed = .false.
IF (tNoSymGenRandExcits) THEN
SymA = 0
SymB = 0
MlA = 0
MlB = 0
else if (tKPntSym) THEN
SymA = SpinOrbSymLabel(OrbA)
SymB = RandExcitSymLabelProd(SymInvLabel(SymA), SymProduct)
MlB = 0
MlA = 0
ELSE
SymA = INT(G1(OrbA)%Sym%S, 4)
SymB = IEOR(SymA, SymProduct)
MlB = 0
MlA = 0
IF (tFixLz) THEN
MlA = G1(OrbA)%Ml
MlB = SumMl - MlA
end if
end if
IF (abs(MlB) <= iMaxLz) THEN
!Make sure that the B orbital that we would need to pick to conserve momentum is actually in the available range of Ml values.
IF (iSpn == 2) THEN
!We want an alpha/beta unocc pair.
IF (SpinOrbA == 1) THEN
!We have picked an alpha orbital - check to see if there are allowed beta unoccupied orbitals from the SymB Class.
IF (ClassCountUnocc2(ClassCountInd(2, SymB, MlB)) /= 0) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
ELSE
!We have picked a beta orbital - check to see if there are allowed alpha unoccupied orbitals from the SymB Class.
IF (ClassCountUnocc2(ClassCountInd(1, SymB, MlB)) /= 0) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
end if
else if (iSpn == 1) THEN
!We want a beta/beta pair.
IF ((SymB /= SymA) .or. (MlA /= MlB)) THEN
!Check to see if there are any unoccupied beta orbitals in the SymB Class.
IF (ClassCountUnocc2(ClassCountInd(2, SymB, MlB)) /= 0) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
ELSE
!We want an orbital from the same class. Check that this isn't the only unoccupied beta orbital in the class.
IF (ClassCountUnocc2(ClassCountInd(2, SymB, MlB)) /= 1) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
end if
ELSE
!We want an alpha/alpha pair.
IF ((SymA /= SymB) .or. (MlA /= MlB)) THEN
!Check to see if there are any unoccupied alpha orbitals in the SymB Class.
IF (ClassCountUnocc2(ClassCountInd(1, SymB, MlB)) /= 0) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
ELSE
!We want an orbital from the same class. Check that this isn't the only unoccupied alpha orbital in the class.
IF (ClassCountUnocc2(ClassCountInd(1, SymB, MlB)) /= 1) THEN
!Success! We have found an allowed A orbital!
IsAOrbSymAllowed = .true.
end if
end if
end if
end if
END FUNCTION IsAOrbSymAllowed
!This routine takes determinant nI and returns two randomly chosen electrons, whose index in nI is Elec1Ind and Elec2Ind.
!These electrons have symmetry product SymProduct and spin pairing iSpn, where iSpn = 1=beta/beta; 2=alpha/beta; 3=alpha/alpha.
!If IndInp = -1, the pair is picked randomly with prob = 1/ElecPairs. Otherwise, it will choose electron pair given by index IndInp.
SUBROUTINE PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, IndInp)
INTEGER, INTENT(IN) :: nI(NEl), IndInp
INTEGER, INTENT(OUT) :: Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl
INTEGER :: Ind, X, K, Orb1, Orb2
real(dp) :: r
!Triangular indexing system.
!This is used for picking two distinct electrons out of all N(N-1)/2 pairs.
!
! 12 13 14 15 1 2 3 4
! 23 24 25 5 6 7
! 34 35 => 8 9
! 45 10
!For a given index, ind, there are [N(N-1)/2 - ind] elements at positions larger
!than ind. Call this number X.
!We want to find out how many rows there are after the row containing the element
!ind. t rows has t(t+1)/2 elements in it.
!Therefore, to find out the number of rows after ind, we want to find the largest K, such that K(K+1)/2 =< X
!The solution to this is that K =< (SQRT(8*X+1)-1)/2, therefore we can find K (the
!largest integer which satisfies this).
!We then know the number of rows after the element ind. Therefore, since there are N-1
!rows in total, we know we are on row N-1-K.
!This gives us the index of the first electron.
!To find the second electron (i.e. the column), we know that out of the X elements in
!positions larger than ind, K(K+1)/2 are in the next rows.
!This means that X - K(K+1)/2 are in the same row. There are N-(N-1-K) = 1+K elements
!in the row chosen, and so the number of elements into the
!row it is, is given by (1+K) - (X-K(K+1)/2). However, in row z, the first column
!index is z+1. Therefore the index of the second electron is
!(1+K) - (X-K(K+1)/2) + N-K-1 = N-X+(K(K+1)/2).
! ElecPairs=(NEl*(NEl-1))/2
! If we want to find an index randomly, IndInp will be -1.
IF (IndInp == -1) THEN
!Find an index randomly.
r = genrand_real2_dSFMT()
Ind = INT(ElecPairs * r) + 1
ELSE
!Otherwise we are looking for a specific electron pair specified by IndInp
Ind = IndInp
end if
!X is number of elements at positions larger than ind
X = ElecPairs - Ind
!K is the number of complete rows after the element ind
K = INT((SQRT(8.0_dp * REAL(X, dp) + 1.0_dp) - 1.0_dp) / 2.0_dp)
Elec1Ind = NEl - 1 - K
Elec2Ind = NEl - X + ((K * (K + 1)) / 2)
Orb1 = nI(Elec1Ind)
Orb2 = nI(Elec2Ind)
!We now want to find the symmetry product label of the two electrons, and the spin product of the two electrons.
SymProduct = RandExcitSymLabelProd(SpinOrbSymLabel(Orb1), SpinOrbSymLabel(Orb2))
IF ((G1(Orb1)%Ms) * (G1(Orb2)%Ms) == -1) THEN
!We have an alpha beta pair of electrons.
iSpn = 2
ELSE
IF (G1(Orb1)%Ms == 1) THEN
!We have an alpha alpha pair of electrons.
iSpn = 3
ELSE
!We have a beta beta pair of electrons.
iSpn = 1
end if
end if
SumMl = G1(Orb1)%Ml + G1(Orb2)%Ml
END SUBROUTINE PickElecPair
SUBROUTINE CheckIfSingleExcits(ElecsWNoExcits, ClassCount2, ClassCountUnocc2, nI)
INTEGER, intent(out) :: ElecsWNoExcits
integer, intent(in) :: ClassCount2(ScratchSize), ClassCountUnocc2(ScratchSize)
INTEGER, intent(in) :: nI(NEl)
integer :: i
!First, 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
IF (tFixLz .or. tKPntSym) THEN
!Here, we also have to check that the electron is momentum allowed.
!Since there are many more irreps, it will be quicker here to check all electrons, rather than all the symmetries.
do i = 1, NEl
IF (G1(nI(i))%Ms == 1) THEN
IF (ClassCountUnocc2(ClassCountInd(1, SpinOrbSymLabel(nI(i)), G1(nI(i))%Ml)) == 0) THEN
ElecsWNoExcits = ElecsWNoExcits + 1
end if
ELSE
IF (ClassCountUnocc2(ClassCountInd(2, SpinOrbSymLabel(nI(i)), G1(nI(i))%Ml)) == 0) THEN
ElecsWNoExcits = ElecsWNoExcits + 1
end if
end if
end do
! do i=1,ScratchSize
!!Run through all labels
! IF((ClassCount2(i).ne.0).and.(ClassCountUnocc2(i).eq.0)) THEN
!!If there are electrons in this class with no possible unoccupied orbitals in the
!same class, these electrons have no single excitations.
! ElecsWNoExcits=ElecsWNoExcits+ClassCount2(i)
! end if
! end do
ELSE
do i = 1, ScratchSize
!Run through all labels
IF ((ClassCount2(i) /= 0) .and. (ClassCountUnocc2(i) == 0)) THEN
!If there are electrons in this class with no possible unoccupied orbitals in the same
!class, these electrons have no single excitations.
ElecsWNoExcits = ElecsWNoExcits + ClassCount2(i)
end if
end do
end if
!Rather than choosing a double now if there are no singles, just return a null det.
! IF(ElecsWNoExcits.eq.NEl) THEN
!!There are no single excitations from this determinant at all. This means the
!probability to create a double excitation = 1
!!Then we will create a double excitation instead.
! pDoubNew=1.0_dp
! RETURN
! end if
!
END SUBROUTINE CheckIfSingleExcits
!> Wrapper function for creating a uniform single excitation
subroutine uniform_single_excit_wrapper(nI, ilutI, nJ, ilutJ, ex, tpar, store, pgen)
implicit none
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: nJ(nel), ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:NIfTot)
logical, intent(out) :: tpar
real(dp), intent(out) :: pGen
type(excit_gen_store_type), intent(inout), target :: store
! First, count the number of available orbitals per irrep
if (.not. store%tFilled) then
call construct_class_counts(nI, store%ClassCountOcc, &
store%ClassCountUnocc)
store%tFilled = .true.
end if
pDoubNew = 1.0 - pSingles
! Then, chose uniformly from them
call createSingleExcit(nI, nJ, store%ClassCountOcc, store%classCountUnocc, ilutI, &
ex, tPar, pGen)
if (nJ(1) /= 0) then
! set the output ilut
ilutJ = ilutI
clr_orb(ilutJ, ex(1, 1))
set_orb(ilutJ, ex(2, 1))
else
ilutJ = 0_n_int
end if
end subroutine uniform_single_excit_wrapper
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
subroutine construct_class_counts(nI, CCOcc, CCUnocc)
! Return two arrays of length ScratchSize, containing information on
! the number of orbitals - occupied and unoccupied - in each symmetry.
!
! The arrays are indexed via the indices returned by ClassCountInd
! n.b. this is O[nel], so we should store this if we can.
integer, intent(in) :: nI(nel)
integer, intent(out) :: CCOcc(ScratchSize), CCUnocc(ScratchSize)
integer :: ind_alpha, ind_beta, i, ind
CCOcc = 0
CCUnocc = OrbClassCount
if (tNoSymGenRandExcits .or. t_new_real_space_hubbard) then
ind_alpha = ClassCountInd(1, 0, 0)
ind_beta = ClassCountInd(2, 0, 0)
CCOcc(ind_alpha) = nOccAlpha
CCOcc(ind_beta) = nOccBeta
CCUnocc(ind_alpha) = CCUnocc(ind_alpha) - nOccAlpha
CCUnocc(ind_beta) = CCUnocc(ind_beta) - nOccBeta
else
do i = 1, nel
ind = ClasscountInd(nI(i))
CCOcc(ind) = CCOcc(ind) + 1
CCUnocc(ind) = CCUnocc(ind) - 1
end do
end if
end subroutine
subroutine calc_pgen_symrandexcit2(nI, ex, ic, ClassCount2, &
ClassCountUnocc2, pDoub, pGen)
! This routine will calculate the PGen between two connected
! determinants, nI and nJ which are IC excitations of each other, using
! the unbiased scheme.
!
! Only the excitation matrix is needed (1,*) are the i,j orbs, and
! (2,*) are the a,b orbs. This is the prob of generating nJ FROM nI,
! not the other way round. Passed in is also the ClassCount2 arrays for
! nI, and the probability of picking a double.
!
! A word of warning: The routine does not check that the determinants
! are indeed connected, and may well return a non-zero probability even
! if they arent. Therefore, make sure that they are at most double
! excitations of each other.
!
! nI is the determinant from which the excitation comes from.
integer, intent(in) :: nI(nel), ex(2, 2), ic
integer, intent(in) :: ClassCount2(ScratchSize)
integer, intent(in) :: ClassCountUnocc2(ScratchSize)
real(dp), intent(in) :: pDoub
real(dp), intent(out) :: pGen
integer :: ForbiddenOrbs, symA, symB, sumMl, MlA, MlB, Elec1Ml
integer :: ElecsWNoExcits, NExcitOtherWay, OrbI, OrbJ, iSpn, NExcitA
integer :: NExcitB, ElecSym, orbA, orbB
IF (IC == 1) THEN
!First, 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.
CALL CheckIfSingleExcits(ElecsWNoExcits, ClassCount2, ClassCountUnocc2, nI)
IF (tNoSymGenRandExcits) THEN
!Find symmetry of chosen electron
ElecSym = 0
ELSE
ElecSym = SpinOrbSymLabel((Ex(1, 1)))
IF (tFixLz) THEN
Elec1Ml = G1(Ex(1, 1))%Ml
end if
end if
IF (G1(Ex(1, 1))%Ms == 1) THEN
!Alpha orbital - see how many single excitations there are from this electron...
NExcitA = ClassCountUnocc2(ClassCountInd(1, ElecSym, Elec1Ml))
ELSE
!Beta orbital
NExcitA = ClassCountUnocc2(ClassCountInd(2, ElecSym, Elec1Ml))
end if
!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)
!Prob of generating a single is 1-pDoub
! pGen=(1.0_dp-pDoub)*(1.0_dp/real(NEl,dp))*(1.0_dp/real(NExcitA,dp))*((real(NEl,dp))/(real((NEl-ElecsWNoExcits),dp)))
pGen = (1 - pDoub) / (REAL((NExcitA * (NEl - ElecsWNoExcits)), dp))
ELSE
!Prob of generating a double excitation.
!Find the I and J orbitals
OrbI = Ex(1, 1)
OrbJ = Ex(1, 2)
OrbA = Ex(2, 1)
OrbB = Ex(2, 2)
! write(stdout,*) OrbI,OrbJ,OrbA,OrbB
!Find the spin-product of the occupied pair
IF ((G1(OrbI)%Ms) * (G1(OrbJ)%Ms) == -1) THEN
!We have an alpha beta pair of electrons.
iSpn = 2
!NExcit is the number of allowed unoccupied orbitals to choose
NExcitA = nBasis - NEl
ELSE
IF (G1(OrbI)%Ms == 1) THEN
!We have an alpha alpha pair of electrons.
iSpn = 3
NExcitA = (nBasis / 2) - nOccAlpha !This is the number of unocc alpha spinorbs
ELSE
!We have a beta beta pair of electrons.
iSpn = 1
NExcitA = (nBasis / 2) - nOccBeta !This is the number of unocc beta spinorbs
end if
end if
IF (tFixLz) THEN
SumMl = G1(OrbI)%Ml + G1(OrbJ)%Ml
MlA = G1(OrbA)%Ml
MlB = SumMl - MlA
ELSE
SumMl = 0
MlA = 0
MlB = 0
end if
IF (tNoSymGenRandExcits) THEN
ElecSym = 0
CALL FindNumForbiddenOrbsNoSym(ForbiddenOrbs, ClassCountUnocc2, iSpn)
SymA = 0
SymB = 0
ELSE
!Calculate the symmetry product of the occupied orbital pair
ElecSym = RandExcitSymLabelProd(SpinOrbSymLabel(OrbI), SpinOrbSymLabel(OrbJ))
! ElecSym=INT(IEOR(G1(OrbI)%Sym%S,G1(OrbJ)%Sym%S),4)
!This will calculate the A orbitals which will have no B pair
CALL FindNumForbiddenOrbs(ForbiddenOrbs, ClassCountUnocc2, ElecSym, iSpn, SumMl)
!Need to find the symmetries of the unoccupied A and B orbitals.
SymA = SpinOrbSymLabel(OrbA)
SymB = RandExcitSymLabelProd(SymInvLabel(SymA), ElecSym)
! SymA=INT(G1(OrbA)%Sym%S,4)
! SymB=IEOR(SymA,ElecSym)
end if
! write(stdout,*) "Check: ",ForbiddenOrbs,ElecSym,iSpn,SymA,SymB,OrbA,OrbB
!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.
IF (iSpn == 2) THEN
!If iSpn=2, then we want to find a spinorbital of the opposite spin of SpinOrbA
IF ((G1(OrbA)%Ms) == -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.
NExcitB = ClassCountUnocc2(ClassCountInd(1, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(2, SymA, MlA))
ELSE
!Want to pick an beta orbital.
NExcitB = ClassCountUnocc2(ClassCountInd(2, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(1, SymA, MlA))
end if
else if (iSpn == 1) THEN
!Definitely want a beta orbital
NExcitB = ClassCountUnocc2(ClassCountInd(2, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(2, SymA, MlA))
ELSE
!Definitely want an alpha orbital
NExcitB = ClassCountUnocc2(ClassCountInd(1, SymB, MlB))
NExcitOtherWay = ClassCountUnocc2(ClassCountInd(1, SymA, MlA))
end if
IF ((iSpn /= 2) .and. (SymA == SymB) .and. (MlB == MlA)) 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 ElecSym=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.
NExcitB = NExcitB - 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
pGen = pDoub * ((1.0_dp / real(NExcitB, dp)) + &
(1.0_dp / real(NExcitOtherWay, dp))) / (REAL((ElecPairs * (NExcitA - ForbiddenOrbs)), dp))
end if
end subroutine
!*********************** BIASED EXCITATION GENERATION ROUTINES *************************!
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
SUBROUTINE CreateDoubExcitBiased(nI, nJ, iLut, ExcitMat, tParity, nParts, WSign, Tau, iCreate)
INTEGER :: nI(NEl), nJ(NEl), ExcitMat(2, maxExcit), iCreate, iSpn, OrbA, OrbB, SymProduct
INTEGER(KIND=n_int) :: iLut(0:NIfTot)
INTEGER :: Elec1Ind, Elec2Ind, nParts, SumMl
REAL(dp) :: WSign
HElement_t(dp) :: rh
LOGICAL :: tParity
real(dp) :: Tau
!First, we need to pick an unbiased distinct electron pair.
!These have symmetry product SymProduct, and spin pair iSpn = 1=beta/beta; 2=alpha/beta; 3=alpha/alpha
!The probability for doing this is 1/ElecPairs.
CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, -1)
!This routine runs through all distinct ab pairs for the chosen ij and stochastically chooses how many particles to create.
!If spawning wants to occur, then it runs through the list again and chooses a pair, which it returns.
CALL CalcAllab(nI, iLut, Elec1Ind, Elec2Ind, SymProduct, iSpn, OrbA, OrbB, nParts, iCreate, Tau)
if (tGUGA) then
call stop_all("CreateDoubExcitBiased", &
"modify get_helement for GUGA")
end if
!We now know that we want to create iCreate particles, from orbitals nI(Elec1/2Ind) -> OrbA + OrbB.
IF (iCreate > 0) THEN
call make_double(nI, nJ, elec1ind, elec2ind, orbA, orbB, &
ExcitMat, tParity)
!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, 2, ExcitMat, tParity)
IF (WSign > 0) 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 CreateDoubExcitBiased
SUBROUTINE CalcAllab(nI, ILUT, Elec1Ind, Elec2Ind, SymProduct, iSpn, OrbA, OrbB, nParts, iCreate, Tau)
INTEGER :: nI(NEl), Elec1Ind, Elec2Ind, SymProduct, iSpn, OrbA, OrbB, iCreate
INTEGER(KIND=n_int) :: iLut(0:NIfTot)
INTEGER :: SpatOrbi, SpatOrbj, Spini, Spinj, i, aspn, bspn, SymA, SymB, SpatOrba, EndSymState, VecInd
real(dp) :: Tau, SpawnProb(MaxABPairs), NormProb, rat, r
INTEGER :: SpawnOrbs(2, MaxABPairs), j, nParts, SpinIndex, Ind
HElement_t(dp) :: HEl
!We want the spatial orbital number for the ij pair (Elec1Ind is the index in nI).
!Later, we'll have to use GTID for UHF.
SpatOrbi = ((nI(Elec1Ind) - 1) / 2) + 1
SpatOrbj = ((nI(Elec2Ind) - 1) / 2) + 1
Spini = G1(nI(Elec1Ind))%Ms
Spinj = G1(nI(Elec2Ind))%Ms
VecInd = 1
NormProb = 0.0_dp
do i = 1, nBasis
!Run through all a orbitals
IF (mod(i, 2) == 0) THEN
!We have an alpha spin...
aspn = 1
IF (iSpn == 1) THEN
!ij is an beta/beta pair, so we only want to pick beta a orbitals.
CYCLE
end if
ELSE
!We have a beta spin...
aspn = -1
IF (iSpn == 3) THEN
!ij is a alpha/alpha pair, so we only want to pick alpha a orbitals.
CYCLE
end if
end if
!We also want to check that the a orbital we have picked is not already in the determinant we are exciting from.
IF (BTEST(ILUT((i - 1) / bits_n_int), MOD((i - 1), bits_n_int))) THEN
!Orbital is in nI...do not make a pair from this.
CYCLE
end if
!Now we have to run over all b's which can go with a. We only want unique pairs, so we constrain b to be at least a+1.
!We only want to run over symmetry and spin allowed b's though.
!Find the required symmetry of b.
SymA = INT(G1(i)%Sym%S, 4)
SymB = IEOR(SymA, SymProduct)
SpatOrba = ((i - 1) / 2) + 1
!We also want to take into account spin.
IF (ispn == 1) THEN
bspn = -1 !Want beta spin b orbitals
SpinIndex = 2
else if (ispn == 3) THEN
bspn = 1 !Want alpha spin b orbitals
SpinIndex = 1
ELSE
!ij pair is an alpha/beta spin pair, therefore b wants to be of opposite spin to a.
IF (aspn == -1) THEN
!a orbital is a beta orbital, therefore we want b to be an alpha orbital.
bspn = 1
SpinIndex = 1
ELSE
bspn = -1
SpinIndex = 2
end if
end if
!To run just through the states of the required symmetry we want to use SymLabelCounts.
! StartSymState=SymLabelCounts(1,SymB+1)
Ind = ClassCountInd(SpinIndex, SymB, 0)
EndSymState = SymLabelCounts2(1, Ind) + SymLabelCounts2(2, Ind) - 1
!Run over all possible b orbitals
do j = SymLabelCounts2(1, Ind), EndSymState
! IF(bspn.eq.-1) THEN
! OrbB=(2*SymLabelList(j))-1 !This is the spin orbital chosen for b
! ELSE
! OrbB=(2*SymLabelList(j))
! end if
OrbB = SymLabelList2(j) !This is the spin orbital chosen for b
IF (OrbB <= i) THEN
!Since we only want unique ab pairs, ensure that b > a.
CYCLE
end if
IF (BTEST(ILUT((OrbB - 1) / bits_n_int), MOD((OrbB - 1), bits_n_int))) THEN
!Orbital is in nI...do not make a pair from this.
CYCLE
end if
!We have now found an allowed ab pair to go with the ij pair chosen previously - record its stats.
IF (Spini == aspn .and. Spinj == bspn) THEN
Hel = get_umat_el(SpatOrbi, SpatOrbj, SpatOrba, j)
ELSE
Hel = (0.0_dp)
end if
IF (Spini == bspn .and. Spinj == aspn) THEN
Hel = Hel - get_umat_el(SpatOrbi, SpatOrbj, j, SpatOrba)
end if
SpawnProb(VecInd) = abs(REAL(Hel, dp))
SpawnOrbs(1, VecInd) = i
SpawnOrbs(2, VecInd) = OrbB
NormProb = NormProb + SpawnProb(VecInd)
VecInd = VecInd + 1
IF (VecInd >= MaxABPairs) THEN
CALL Stop_All("CalcAllab", "Finding too many ab pairs...")
end if
end do
end do
VecInd = VecInd - 1 !This now indicates the total number of ab pairs we have found for the chosen ij.
IF (VecInd == 0) THEN
CALL Stop_All("CalcAllab", "No ab pairs found for the chosen ij")
end if
!We now have to find out how many children to spawn, based on the value of normprob.
rat = Tau * NormProb * REAL(ElecPairs * nParts, 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 (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) THEN
CALL Stop_All("CalcAllab", "Chosen ab pair does not correspond to allowed pair")
end if
OrbA = SpawnOrbs(1, i)
OrbB = SpawnOrbs(2, i)
end if
END SUBROUTINE CalcAllab
! For a given ij pair in the UEG or Hubbard model, this generates ab as a double excitation efficiently.
! This takes into account the momentum conservation rule, i.e. that kb=ki+ki-ka(+G).
SUBROUTINE CreateDoubExcitLattice(nI, iLutnI, nJ, tParity, ExcitMat, pGen, Elec1Ind, Elec2Ind, iSpn, part_type)
INTEGER :: i, nI(NEl), nJ(NEl), Elec1Ind, Elec2Ind, iSpn, kb_ms
INTEGER(KIND=n_int) :: iLutnI(0:NIfTot)
INTEGER :: ChosenUnocc, Hole1BasisNum, Hole2BasisNum, ki(3), kj(3), ka(3), kb(3), ExcitMat(2, maxExcit), iSpinIndex, TestEnergyB
LOGICAL :: tAllowedExcit, tParity
real(dp) :: r, pGen, pAIJ
integer, intent(in), optional :: part_type
integer :: elecs(2), loc, temp_part_type
! [W.D.] just a workaround for now, after i rehaul all the lattice
! excitation generators, this will be more optimized
if (present(part_type)) then
temp_part_type = part_type
else
temp_part_type = 1
end if
! This chooses an a of the correct spin, excluding occupied orbitals
! This currently allows b orbitals to be created that are disallowed
hole2basisnum = 0
DO
r = genrand_real2_dSFMT()
! Choose a
IF (iSpn == 2) THEN ! alpha/beta combination
ChosenUnocc = INT(nBasis * r) + 1
ELSE ! alpha/alpha, beta/beta
ChosenUnocc = 2 * (INT(nBasis / 2 * r) + 1) & ! 2*(a number between 1 and nbasis/2) gives the alpha spin
& - (1 - (iSpn / 3)) ! alpha numbered even, iSpn/3 returns 1 for alpha/alpha, 0 for beta/beta
end if
! Check a isn't occupied
IF (.not. (BTEST(iLutnI((ChosenUnocc - 1) / bits_n_int), MOD(ChosenUnocc - 1, bits_n_int)))) THEN
!Orbital not in nI. Accept.
EXIT
end if
end do
Hole1BasisNum = ChosenUnocc
IF ((Hole1BasisNum <= 0) .or. (Hole1BasisNum > nBasis)) THEN
CALL Stop_All("CreateDoubExcitLattice", "Incorrect basis function generated")
end if
!=============================================
if (tUEG2) then
! kb is now uniquely defined
ki = kvec(nI(Elec1Ind), 1:3)
kj = kvec(nI(Elec2Ind), 1:3)
ka = kvec(Hole1BasisNum, 1:3)
kb = ki + kj - ka
! Find the spin of b
IF (iSpn == 2) THEN ! alpha/beta required, therefore b has to be opposite spin to a
kb_ms = 1 - ((G1(Hole1BasisNum)%Ms + 1) / 2) * 2
ELSE ! b is the same spin as a
kb_ms = G1(Hole1BasisNum)%Ms
end if
! Is kb allowed by the size of the space?
tAllowedExcit = .true.
TestEnergyB = kb(1)**2 + kb(2)**2 + kb(3)**2
IF (tOrbECutoff .and. (TestEnergyB > OrbECutoff)) tAllowedExcit = .false.
IF (.not. tAllowedExcit) THEN
nJ(1) = 0
RETURN
end if
iSpinIndex = (kb_ms + 1) / 2 + 1
Hole2BasisNum = kPointToBasisFn(kb(1), kb(2), kb(3), iSpinIndex)
IF (Hole2BasisNum == -1 .or. Hole1BasisNum == Hole2BasisNum) THEN
nJ(1) = 0
RETURN
end if
! Is b occupied?
IF (BTEST(iLutnI((Hole2BasisNum - 1) / bits_n_int), MOD(Hole2BasisNum - 1, bits_n_int))) THEN
!Orbital is in nI. Reject.
tAllowedExcit = .false.
end if
IF (.not. tAllowedExcit) THEN
nJ(1) = 0
RETURN
end if
! Check that the correct kb has been found -- can be commented out later
DO i = 1, 3
IF ((kvec(nI(Elec2Ind), i) + kvec(nI(Elec1Ind), i) &
- kvec(Hole1BasisNum, i) - kvec(Hole2BasisNum, i)) /= 0) THEN
write(stdout, *) "Tried to excite "
write(stdout, *) "ki ", ki
write(stdout, *) "kj ", kj
write(stdout, *) "ka ", ka
write(stdout, *) "kb should be ", kb
write(stdout, *) "but found as ", kvec(Hole2BasisNum, 1:3)
CALL Stop_All("CreateDoubExcitLattice", "Wrong b found")
end if
end do
! Find the new determinant
call make_double(nI, nJ, elec1ind, elec2ind, Hole1BasisNum, &
Hole2BasisNum, ExcitMat, tParity)
!Calculate generation probabilities
IF (iSpn == 2) THEN
pAIJ = 1.0_dp / (nBasis - Nel)
else if (iSpn == 1) THEN
pAIJ = 1.0_dp / (nBasis / 2 - nOccBeta)
ELSE
!iSpn = 3
pAIJ = 1.0_dp / (nBasis / 2 - nOccAlpha)
end if
! Note, p(b|ij)=p(a|ij) for this system
pGen = 2.0_dp / (NEl * (NEl - 1)) * 2.0_dp * pAIJ
return
end if
!=============================================
! kb is now uniquely defined
ki = G1(nI(Elec1Ind))%k
kj = G1(nI(Elec2Ind))%k
ka = G1(Hole1BasisNum)%k
kb = ki + kj - ka
! Find the spin of b
IF (iSpn == 2) THEN ! alpha/beta required, therefore b has to be opposite spin to a
kb_ms = 1 - ((G1(Hole1BasisNum)%Ms + 1) / 2) * 2
ELSE ! b is the same spin as a
kb_ms = G1(Hole1BasisNum)%Ms
end if
IF (tHub) THEN
! Re-map back into cell
CALL MomPbcSym(kb, nBasisMax)
! This is the look-up table method of finding the kb orbital
iSpinIndex = (kb_ms + 1) / 2 + 1
! [W.D.]
! make the access to kpoint same in hubbard and ueg
! Hole2BasisNum=kPointToBasisFn(kb(1),kb(2),1,iSpinIndex)
Hole2BasisNum = kPointToBasisFn(kb(1), kb(2), kb(3), iSpinIndex)
!Hole2BasisNum will be -1 if that orb is frozen
! Is b occupied?
IF (Hole2BasisNum == -1 .or. BTEST(iLutnI((Hole2BasisNum - 1) / bits_n_int), MOD(Hole2BasisNum - 1, bits_n_int))) THEN
!Orbital is in nI. Reject.
nJ(1) = 0
RETURN
end if
! Reject if a=b.
IF (Hole1BasisNum == Hole2BasisNum) THEN
nJ(1) = 0
RETURN
end if
end if
IF (tUEG) THEN
! Is kb allowed by the size of the space?
! Currently only applies when NMAXX etc. are set by the CELL keyword
! Not sure what happens when an energy cutoff is set
tAllowedExcit = .true.
IF (ABS(kb(1)) > NMAXX) tAllowedExcit = .false.
IF (ABS(kb(2)) > NMAXY) tAllowedExcit = .false.
IF (ABS(kb(3)) > NMAXZ) tAllowedExcit = .false.
TestEnergyB = kb(1)**2 + kb(2)**2 + kb(3)**2
IF (tOrbECutoff .and. (TestEnergyB > OrbECutoff)) tAllowedExcit = .false.
IF (.not. tAllowedExcit) THEN
nJ(1) = 0
RETURN
end if
! Which orbital has momentum kb?
!IF (iSpn.eq.2) THEN
! Hole2BasisNum=2*((NMAXZ*2+1)*(NMAXY*2+1)*(kb(1)+NMAXX)+(NMAXZ*2+1)*
!(kb(2)+NMAXY)+(kb(3)+NMAXZ)+1)-(1+G1(Hole1BasisNum)%Ms)/2
!ELSE
! Hole2BasisNum=2*((NMAXZ*2+1)*(NMAXY*2+1)*(kb(1)+NMAXX)+(NMAXZ*2+1)*
!(kb(2)+NMAXY)+(kb(3)+NMAXZ)+1)-1+iSpn/3
! end if
iSpinIndex = (kb_ms + 1) / 2 + 1
Hole2BasisNum = kPointToBasisFn(kb(1), kb(2), kb(3), iSpinIndex)
IF (Hole2BasisNum == -1 .or. Hole1BasisNum == Hole2BasisNum) THEN
nJ(1) = 0
RETURN
end if
! Is b occupied?
IF (BTEST(iLutnI((Hole2BasisNum - 1) / bits_n_int), MOD(Hole2BasisNum - 1, bits_n_int))) THEN
!Orbital is in nI. Reject.
tAllowedExcit = .false.
end if
IF (.not. tAllowedExcit) THEN
nJ(1) = 0
RETURN
end if
! Check that the correct kb has been found -- can be commented out later
DO i = 1, 3
IF ((G1(nI(Elec2Ind))%k(i) + G1(nI(Elec1Ind))%k(i) - G1(Hole1BasisNum)%k(i) - G1(Hole2BasisNum)%k(i)) /= 0) THEN
write(stdout, *) "Tried to excite "
write(stdout, *) "ki ", ki
write(stdout, *) "kj ", kj
write(stdout, *) "ka ", ka
write(stdout, *) "kb should be ", kb
write(stdout, *) "but found as ", G1(Hole2BasisNum)%k
CALL Stop_All("CreateDoubExcitLattice", "Wrong b found")
end if
end do
end if
! Find the new determinant
call make_double(nI, nJ, elec1ind, elec2ind, Hole1BasisNum, &
Hole2BasisNum, ExcitMat, tParity)
IF (tHub) THEN
! Debug to test the resultant determinant
! i think i also have to check here if kpoint symmetry is turned
! on... or is it generally turned on in the hubbard case??
IF (.not. (IsMomentumAllowed(nJ)) .and. tKPntSym) THEN
CALL Stop_All("CreateDoubExcitLattice", "Incorrect kb generated -- momentum not conserved")
end if
end if
!Calculate generation probabilities
IF (iSpn == 2) THEN
! otherwise take the pAIJ set above
pAIJ = 1.0_dp / (nBasis - Nel)
else if (iSpn == 1) THEN
pAIJ = 1.0_dp / (nBasis / 2 - nOccBeta)
ELSE
!iSpn = 3
pAIJ = 1.0_dp / (nBasis / 2 - nOccAlpha)
end if
! Note, p(b|ij)=p(a|ij) for this system
if (tUEG) then
pGen = 2.0_dp / (NEl * (NEl - 1)) * 2.0_dp * pAIJ
else ! i.e. if hubbard model, use modified probabilities
! hubbard model can't spawn alpha/alpha and beta/beta type excitations
pGen = 1.0_dp / (nOccAlpha * nOccBeta) * 2.0_dp * pAIJ
end if
END SUBROUTINE CreateDoubExcitLattice
! This creates a general excitation for the UEG and Hubbard model. Each ij pair is
!picked with equal probability, and then
! each ab is picked for that ij pair with equal probability.
! For UEG there is a more sophisticated algorithm that allows more ijab choices to be
!rejected before going back to the main
! code.
SUBROUTINE CreateExcitLattice(nI, iLutnI, nJ, tParity, ExcitMat, pGen, part_type)
INTEGER :: i, j ! Loop variables
INTEGER :: Elec1, Elec2
INTEGER :: nI(NEl), nJ(NEl), Elec1Ind, Elec2Ind, ExcitMat(2, maxExcit), iSpn
INTEGER(KIND=n_int) :: iLutnI(0:NIfTot)
INTEGER :: ki(3), kj(3), kTrial(3), iElecInExcitRange, iExcludedKFromElec1, iAllowedExcites
INTEGER :: KaXLowerLimit, KaXUpperLimit, KaXRange, KaYLowerLimit, KaYUpperLimit, KaYRange, KaZLowerLimit, KaZUpperLimit, KaZRange
LOGICAL :: tParity, tDoubleCount, tExtraPoint
real(dp) :: r(2), pGen, pAIJ
integer, intent(in), optional :: part_type
INTEGER, ALLOCATABLE :: Excludedk(:, :)
character(*), parameter :: this_routine = "CreateExcitLattice"
integer :: elecs(2), src(2), sym_prod, sum_ml
real(dp) :: pgen_back
integer :: temp_part_type
! CALL PickElecPair(nI,Elec1Ind,Elec2Ind,SymProduct,iSpn,SumMl,-1)
! Completely random ordering of electrons is important when considering ij->ab ij/->ba.
!This affects pgens for alpha/beta pairs.
kaxrange = 0
ielecinexcitrange = 0
kazrange = 0
kayrange = 0
if (present(part_type)) then
temp_part_type = part_type
else
temp_part_type = 1
end if
DO
r(1) = genrand_real2_dSFMT()
Elec1 = INT(r(1) * NEl + 1)
DO
r(2) = genrand_real2_dSFMT()
Elec2 = INT(r(2) * NEl + 1)
IF (Elec2 /= Elec1) EXIT
end do
Elec1Ind = Elec1
Elec2Ind = Elec2
IF ((G1(nI(Elec1Ind))%Ms == -1) .and. (G1(nI(Elec2Ind))%Ms == -1)) THEN
iSpn = 1
ELSE
IF ((G1(nI(Elec1Ind))%Ms == 1) .and. (G1(nI(Elec2Ind))%Ms == 1)) THEN
iSpn = 3
ELSE
iSpn = 2
end if
end if
IF ((tHub .and. iSpn == 2) .or. (tUEG)) EXIT ! alpha/beta pairs are the only pairs generated for the hubbard model
end do
IF (tNoFailAb) THEN ! pGen is calculated first because there might be no excitations available for this ij pair
IF (tHub) CALL Stop_All("CreateExcitLattice", "This doesn't work with Hubbard Model")
! Find the upper and lower ranges of kx, ky and kz from the point of view of electron i
ki = G1(nI(Elec1Ind))%k
kj = G1(nI(Elec2Ind))%k
!================================================
if (tUEG2) then
ki = kvec(nI(Elec1Ind), 1:3)
kj = kvec(nI(Elec2Ind), 1:3)
end if
!================================================
KaXLowerLimit = MAX(-NMAXX, ki(1) - (NMAXX - kj(1)))
KaXUpperLimit = MIN(NMAXX, ki(1) + (NMAXX + kj(1)))
KaXRange = KaXUpperLimit - KaXLowerLimit + 1
KaYLowerLimit = MAX(-NMAXY, ki(2) - (NMAXY - kj(2)))
KaYUpperLimit = MIN(NMAXY, ki(2) + (NMAXY + kj(2)))
KaYRange = KaYUpperLimit - KaYLowerLimit + 1
KaZLowerLimit = MAX(-NMAXZ, ki(3) - (NMAXZ - kj(3)))
KaZUpperLimit = MIN(NMAXZ, ki(3) + (NMAXZ + kj(3)))
KaZRange = KaZUpperLimit - KaZLowerLimit + 1
! How many disallowed excitations are there due to electrons 'blocking' i->a or j->b excitations
iElecInExcitRange = 0
tExtraPoint = .false.
allocate(Excludedk(NEl, 3))
! Is a = b allowed by momentum and disallowed by spin symmetry?
IF ((G1(nI(Elec1Ind))%Ms == G1(nI(Elec2Ind))%Ms) .and. &
& (MOD(ki(1) - kj(1), 2) == 0) .and. &
& (MOD(ki(2) - kj(2), 2) == 0) .and. &
& (MOD(ki(3) - kj(3), 2) == 0)) THEN ! This is the disallowed double-excitation to the same orbital
iElecInExcitRange = iElecInExcitRange + 1
Excludedk(iElecInExcitRange, :) = (ki + kj) / 2 ! Integer division okay because we're already checked for mod2=0
tExtraPoint = .true.
end if
! How many disallowed i->a are there?
DO i = 1, NEl
IF (G1(nI(i))%Ms /= G1(nI(Elec1Ind))%Ms) CYCLE
kTrial = G1(nI(i))%k
!================================================
if (tUEG2) then
kTrial = kvec(nI(i), 1:3)
end if
!================================================
IF (kTrial(1) < KaXLowerLimit) CYCLE
IF (kTrial(1) > KaXUpperLimit) CYCLE
IF (kTrial(2) < KaYLowerLimit) CYCLE
IF (kTrial(2) > KaYUpperLimit) CYCLE
IF (kTrial(3) < KaZLowerLimit) CYCLE
IF (kTrial(3) > KaZUpperLimit) CYCLE
IF (tExtraPoint) THEN
IF ((kTrial(1) == Excludedk(1, 1)) &
& .and. (kTrial(2) == Excludedk(1, 2)) &
& .and. (kTrial(3) == Excludedk(1, 3))) CYCLE
end if
iElecInExcitRange = iElecInExcitRange + 1
Excludedk(iElecInExcitRange, :) = kTrial
end do
iExcludedKFromElec1 = iElecInExcitRange
! How many disallowed j->b are there, given that some k-points have already been elimated by i->a being disallowed?
DO i = 1, NEl
IF (G1(nI(i))%Ms /= G1(nI(Elec2Ind))%Ms) CYCLE
kTrial = ki + kj - G1(nI(i))%k
!================================================
if (tUEG2) then
kTrial = ki + kj - kvec(nI(i), 1:3)
end if
!================================================
IF (kTrial(1) < KaXLowerLimit) CYCLE
IF (kTrial(1) > KaXUpperLimit) CYCLE
IF (kTrial(2) < KaYLowerLimit) CYCLE
IF (kTrial(2) > KaYUpperLimit) CYCLE
IF (kTrial(3) < KaZLowerLimit) CYCLE
IF (kTrial(3) > KaZUpperLimit) CYCLE
! Need to check for this k-point already having been eliminated
! by the previous loop over electrons
tDoubleCount = .false.
DO j = 1, iExcludedKFromElec1
IF ((Excludedk(j, 1) == kTrial(1)) .and. (Excludedk(j, 2) == kTrial(2)) .and. Excludedk(j, 3) == kTrial(3)) then
tDoubleCount = .true.
end if
end do
IF (.not. tDoubleCount) iElecInExcitRange = iElecInExcitRange + 1
end do
DEallocate(Excludedk)
! If there are no excitations for this ij pair then this subroutine must exit
iAllowedExcites = KaXRange * KaYRange * KaZRange - iElecInExcitRange
IF (iAllowedExcites == 0) THEN
! write(stdout,*) "No allowed excitations from this ij pair"
nJ(1) = 0
pGen = 1.0_dp
RETURN
end if
end if
DO i = 1, 10000
IF (i == 10000) THEN ! Arbitrary termination of the loop to prevent hanging
! due to no excitations being allowed from a certain ij pair
write(stdout, *) "nI:", nI
write(stdout, *) "i & j", nI(Elec1Ind), nI(Elec2Ind)
write(stdout, *) "Allowed Excitations", iAllowedExcites
CALL Stop_All("CreateExcitLattice", "Failure to generate a valid excitation from an electron pair combination")
end if
CALL CreateDoubExcitLattice(nI, iLutnI, nJ, tParity, ExcitMat, pGen, Elec1Ind, Elec2Ind, iSpn, temp_part_type)
IF (.not. tNoFailAb) then
! For an invalid excitation we still need to have a valid value for pGen
if (nJ(1) == 0) pGen = 1.0_dp
RETURN
end if
IF (nJ(1) /= 0) EXIT ! i.e. if we are using the NoFail algorithm only exit on successful nJ(1)!=0
end do
! ***This part of the code is only used if tNoFailAb is ON***
! Else the pgen used is from CreateDoubExcitLattice
! Now calculate pgen
pAIJ = 1.0_dp / (KaXRange * KaYRange * KaZRange - iElecInExcitRange)
! pBIJ is zero for this kind of excitation generator for antiparallel spins
! but is equal to pAIJ for parallel spins.
IF (G1(nI(Elec1Ind))%Ms /= G1(nI(Elec2Ind))%Ms) THEN
pGen = 2.0_dp / (NEl * (NEl - 1)) * pAIJ ! Spins not equal
ELSE
pGen = 2.0_dp / (NEl * (NEl - 1)) * 2.0 * pAIJ ! Spins equal
end if
IF (pAIJ <= 0.0_dp) CALL Stop_All("CreateExcitLattice", "pAIJ is less than 0")
END SUBROUTINE CreateExcitLattice
!As with CalcNonUniPgens:
!This routine will calculate the PGen between two connected determinants, nI and nJ which
!are IC excitations of each other, using the unbiased scheme.
!Only the excitation matrix is needed (1,*) are the i,j orbs, and (2,*) are the a,b orbs
!This routine does it for the lattice models: UEG and hubbard model
SUBROUTINE CalcPGenLattice(Ex, pGen)
INTEGER :: Ex(2, 2), iSpin, jSpin
real(dp) :: pGen, pAIJ
character(*), parameter :: this_routine = "CalcPGenLattice"
ASSERT(.not. tNoFailAb)
! can i make this easier?
iSpin = get_ispn(get_src(ex))
if (iSpin == 1) then
pAIJ = 1.0_dp / (nBasis / 2 - nOccBeta)
else if (iSpin == 2) then
pAIJ = 1.0_dp / (nBasis - Nel)
else if (iSpin == 3) then
pAIJ = 1.0_dp / (nBasis / 2 - nOccAlpha)
end if
! Note, p(b|ij)=p(a|ij) for this system
! there is a small difference here and in the excitation generator
! part. is is only multiplied by 2 if it is a parallel excitation
! in the excitation generator.. but not here anymore..
if (tUEG) then
pGen = 2.0_dp / (NEl * (NEl - 1)) * 2.0_dp * pAIJ
else ! i.e. if hubbard model, use modified probabilities
! hubbard model can't spawn alpha/alpha and beta/beta type excitations
pGen = 1.0_dp / (nOccAlpha * nOccBeta) * 2.0_dp * pAIJ
end if
END SUBROUTINE CalcPGenLattice
FUNCTION IsMomentumAllowed(nJ)
use sym_mod, only: mompbcsym
use SystemData, only: kvec, tUEG2, tAllSymSectors
LOGICAL :: IsMomentumAllowed ! Returns whether the determinant is momentum allowed for
! UEG and Hubbard models
! Compares the total k from a determinant nI with kTotal
INTEGER :: nJ(NEl), kx, ky, kz, ktrial(3), i
if (tAllSymSectors) then
IsMomentumAllowed = .true.
return
end if
IsMomentumAllowed = .false.
!====================================================
if (tUEG2) then
! The momentum constraint for UEG: every determinant must have a total momentum
! which is equal.
kx = 0
ky = 0
kz = 0
do i = 1, NEl
kx = kx + kvec(nJ(i), 1)
ky = ky + kvec(nJ(i), 2)
kz = kz + kvec(nJ(i), 3)
end do
IF (kx == kTotal(1) .and. ky == kTotal(2) .and. kz == kTotal(3)) THEN
IsMomentumAllowed = .true.
end if
return
end if
!====================================================
! The momentum constraint for UEG: every determinant must have a total momentum
! which is equal.
IF (tUEG) THEN
kx = 0
ky = 0
kz = 0
do i = 1, NEl
kx = kx + G1(nJ(i))%k(1)
ky = ky + G1(nJ(i))%k(2)
kz = kz + G1(nJ(i))%k(3)
end do
IF (kx == kTotal(1) .and. ky == kTotal(2) .and. kz == kTotal(3)) THEN
IsMomentumAllowed = .true.
end if
end if
! The momentum constraint from Hubbard model: every determinant must have a total momentum
! which is equal to within a reciprocal lattice vector.
IF (tHub) THEN
kx = 0
ky = 0
kz = 0
do i = 1, NEl
kx = kx + G1(nJ(i))%k(1)
ky = ky + G1(nJ(i))%k(2)
kz = kz + G1(nJ(i))%k(3)
end do
ktrial = (/kx, ky, 0/)
CALL MomPbcSym(ktrial, nBasisMax) ! This re-maps the total momentum under PBCs: equivalent to this being equal to
! a value to within a reciproval lattice vector.
IF (ktrial(1) == kTotal(1) .and. ktrial(2) == kTotal(2)) THEN
IsMomentumAllowed = .true.
else
print *, "ERRONEOUS MOMENTUM: ", ktrial, ktotal
end if
end if
END FUNCTION IsMomentumAllowed
! Initialise or clean the excitation generation storage.
subroutine init_excit_gen_store(store)
type(excit_gen_store_type), intent(inout) :: store
call clean_excit_gen_store(store)
store%tFilled = .false.
allocate(store%ClassCountOcc(ScratchSize1))
allocate(store%ClassCountUnocc(ScratchSize2))
allocate(store%scratch3(ScratchSize3))
if (tBuildOccVirtList) then
allocate(store%occ_list(nel, ScratchSize1))
allocate(store%virt_list(maxval(OrbClassCount), Scratchsize1))
end if
if (tBuildSpinSepLists) then
allocate(store%nI_alpha(nel))
allocate(store%nI_beta(nel))
allocate(store%nI_alpha_inds(nel))
allocate(store%nI_beta_inds(nel))
!store%nel_alpha = 0
end if
if (tSpinProject) then
allocate(store%dorder_i(nel))
allocate(store%dorder_j(nel))
end if
if (t_pcpp_excitgen) allocate(store%elec_map(nel))
end subroutine
subroutine clean_excit_gen_store(store)
type(excit_gen_store_type), intent(inout) :: store
if (associated(store%ClassCountOcc)) then
deallocate(store%ClassCountOcc)
nullify (store%ClassCountOcc)
end if
if (associated(store%ClassCountUnocc)) then
deallocate(store%ClassCountUnocc)
nullify (store%ClassCountUnocc)
end if
if (associated(store%scratch3)) then
deallocate(store%scratch3)
nullify (store%scratch3)
end if
if (associated(store%occ_list)) then
deallocate(store%occ_list)
nullify (store%occ_list)
end if
if (associated(store%virt_list)) then
deallocate(store%virt_list)
nullify (store%virt_list)
end if
if (associated(store%dorder_i)) then
deallocate(store%dorder_i)
nullify (store%dorder_i)
end if
if (associated(store%dorder_j)) then
deallocate(store%dorder_j)
nullify (store%dorder_j)
end if
if (associated(store%elec_map)) then
deallocate(store%elec_map)
nullify (store%elec_map)
end if
end subroutine
!Sometimes (especially UHF orbitals), the symmetry routines will not set up the orbitals correctly.
!Therefore, this routine will set up symlabellist and symlabelcounts
!to be cast in terms of spin orbitals, and the symrandexcit2 routines will use these arrays.
SUBROUTINE SpinOrbSymSetup()
use SymExcitDataMod, only: ScratchSize, ScratchSize1, ScratchSize2, &
SymLabelList2, SymLabelCounts2, kTotal, &
OrbClassCount, kPointToBasisFn, sym_label_list_spat
use SymExcitDataMod, only: SpinOrbSymLabel, SymInvLabel, &
SymTableLabels, KPntInvSymOrb
use SymData, only: nSymLabels, SymClasses, SymLabels
use SystemData, only: NMAXZ, tFixLz, iMaxLz, nBasis, tUEG, tKPntSym, G1, tHub, nBasisMax, NEl
use SystemData, only: Symmetry, tHPHF, tSpn, tISKFuncs, Arr, tNoSymGenRandExcits, Elecpairs
use SystemData, only: MaxABPairs, tUEG2, kvec, t_k_space_hubbard
use umatcache, only: gtID
use sym_mod, only: mompbcsym, SymProd, writesym
use util_mod, only: stop_all
use constants, only: dp, stdout
use lattice_mod, only: lat
IMPLICIT NONE
INTEGER :: i, j, k, SymInd, Lab
INTEGER :: Spin, ierr, OrbSym, InvSym
real(dp) :: OrbEnergy
INTEGER, ALLOCATABLE :: Temp(:)
! These are for the hubbard and UEG model look-up table
INTEGER :: kmaxX, kmaxY, kminX, kminY, kminZ, kmaxz, iSpinIndex, ktrial(3)
type(Symmetry) :: SymProduct, SymI, SymJ
character(len=*), parameter :: this_routine = 'SpinOrbSymSetup'
integer :: sym0
ElecPairs = (NEl * (NEl - 1)) / 2
MaxABPairs = (nBasis * (nBasis - 1) / 2)
! write(stdout,*) "SETTING UP SYMMETRY!!",nBasis,elecpairs
IF (tFixLz) THEN
!Calculate the upper array bound for the ClassCount2 arrays. This will be dependant on the number of symmetries needed.
ScratchSize = 2 * nSymLabels * (2 * iMaxLz + 1)
ELSE
ScratchSize = 2 * nSymLabels !For k-point sym, this can be large
end if
IF (tNoSymGenRandExcits .or. tUEG) THEN
ScratchSize = 2
end if
! We only need the first 2 scratch arrays
ScratchSize1 = ScratchSize
ScratchSize2 = ScratchSize
!Create SpinOrbSymLabel array.
!This array will return a number between 0 and nSymLabels-1.
!For molecular systems, this IS the character of the irrep
!For k-point systems, this is an arbitrary label, and is equal to the standard label - 1.
!This is chosen so that the indexing works with the rest of the excitation generators.
if (allocated(SpinOrbSymLabel)) deallocate(SpinOrbSymLabel)
allocate(SpinOrbSymLabel(nBasis))
do i = 1, nBasis
if (tNoSymGenRandExcits .or. tUEG) then
SpinOrbSymLabel(i) = 0
else if (tKPntSym .or. tHub) then
SpinOrbSymLabel(i) = SymClasses(((i + 1) / 2)) - 1 !This ensures that the symmetry labels go from 0 -> nSymLabels-1
else
SpinOrbSymLabel(i) = INT(G1(i)%Sym%S, 4)
end if
end do
#ifdef DEBUG_
write(stdout, *) "SpinOrbSymLabel: "
do i = 1, nBasis
write(stdout, *) i, SpinOrbSymLabel(i)
end do
#endif
if (tKPntSym) then
allocate(SymTableLabels(0:nSymLabels - 1, 0:nSymLabels - 1))
SymTableLabels(:, :) = -9000 !To make it easier to track bugs
do i = 0, nSymLabels - 1
do j = 0, i
SymI = SymLabels(i + 1) !Convert to the other symlabel convention to use SymLabels -
!TODO: I will fix this to make them consistent when working (ghb24)!
SymJ = SymLabels(j + 1)
SymProduct = SymProd(SymI, SymJ)
!Now, we need to find the label according to this symmetry!
!Run through all symmetries to make working (this could be far more efficient, but its only once, so sod it...
do Lab = 1, nSymLabels
if (SymLabels(Lab)%S == SymProduct%S) then
EXIT
end if
end do
if (Lab == nSymLabels + 1) then
call stop_all("SpinOrbSymSetup", "Cannot find symmetry label")
end if
SymTableLabels(i, j) = Lab - 1
SymTableLabels(j, i) = Lab - 1
end do
end do
#ifdef DEBUG_
write(stdout, *) "SymTable:"
do i = 0, nSymLabels - 1
do j = 0, nSymLabels - 1
write(stdout, "(I6)", advance='no') SymTableLabels(i, j)
end do
write(stdout, *) ""
end do
#endif
end if
!SymInvLabel takes the label (0 -> nSymLabels-1) of a spin orbital, and returns the inverse symmetry label, suitable for
!use in ClassCountInd.
if (allocated(SymInvLabel)) deallocate(SymInvLabel)
allocate(SymInvLabel(0:nSymLabels - 1))
SymInvLabel = -999
! Dongxia changes the gamma point away from center.
! SDS: Provide a default sym0 for cases where this doesn't apply
sym0 = 0
do i = 1, nsymlabels
if (symlabels(i)%s == 0) sym0 = i - 1
end do
do i = 0, nSymLabels - 1
if (tKPntSym) then
! Change the sym label back to the representation used by the rest
! of the code, use SymConjTab, then change back to other rep of
! labels SymConjTab only works when all irreps are self-inverse.
! Therefore, instead, we will calculate the inverses by just
! finding the symmetry which will give A1.
do j = 0, nSymLabels - 1
! Run through all labels to find what gives totally symmetric
! rep
if (SymTableLabels(i, j) == sym0) then
if (SymInvLabel(i) /= -999) then
write(stdout, *) "SymLabel: ", i
call stop_all(this_routine, &
"Multiple inverse irreps found - error")
end if
! This is the inverse
SymInvLabel(i) = j
end if
end do
if (SymInvLabel(i) == -999) then
write(stdout, *) "SymLabel: ", i
call stop_all(this_routine, "No inverse symmetry found - error")
end if
else
! If not using k-point sym, then we are self-inverse
SymInvLabel(i) = i
end if
end do
#ifdef DEBUG_
write(stdout, *) "SymInvLabel: "
do i = 0, nSymLabels - 1
write(stdout, *) i, SymInvLabel(i)
end do
#endif
if (tISKFuncs) then
write(stdout, *) "Setting up inverse orbital lookup for use with ISK functions..."
if (.not. tKPntSym) call stop_all(this_routine, "Cannot use ISK funcs without KPoint symmetry")
if (tSpn) call stop_all(this_routine, "Cannot use ISK on open shell systems")
if (tHPHF) call stop_all(this_routine, "HPHF is not yet working with ISK")
allocate(KPntInvSymOrb(1:nBasis), stat=ierr)
if (ierr /= 0) call stop_all(this_routine, "Allocation error")
do i = 1, nBasis
!Find k-inverse orbital for each spin orbital.
!If the orbital is a self inverse, then just lookup itself.
OrbSym = SpinOrbSymLabel(i)
InvSym = SymInvLabel(OrbSym)
if (InvSym == OrbSym) then
!Orbital is self-inverse
KPntInvSymOrb(i) = i
cycle
end if
!Run through all orbitals looking for inverse
OrbEnergy = Arr(i, 2) !Fock energy of orbital
!Ensure that we get a same-spin orbital
do j = 1, nBasis
if ((SpinOrbSymLabel(j) == InvSym) .and. (mod(i, 2) == mod(j, 2))) then
!This orbital is the right symmetry - is it the inverse orbital? Check Energy.
if ((abs(OrbEnergy - Arr(j, 2))) < 1.0e-7_dp) then
!Assume that this is the inverse orbital.
KPntInvSymOrb(i) = j
exit
end if
end if
end do
if (j > nBasis) then
write(stdout, *) "Orbital: ", i
write(stdout, *) "Symmetry label: ", OrbSym
write(stdout, *) "Inverse Symmetry label: ", InvSym
write(stdout, *) "Fock energy: ", Arr(i, 2)
write(stdout, *) "SpinOrbSymLabel: "
do k = 1, nBasis
write(stdout, *) k, SpinOrbSymLabel(k)
end do
write(stdout, *) "SymInvLavel: "
do k = 0, nSymLabels - 1
write(stdout, *) k, SymInvLabel(k)
end do
call stop_all(this_routine, "Could not find inverse orbital pair for ISK setup.")
end if
end do
do i = 1, nBasis
if (KPntInvSymOrb(i) /= i) exit
end do
if (i > nBasis) then
write(stdout, *) "!! All kpoints are self-inverse, i.e. at the Gamma point or BZ boundary !!"
write(stdout, *) "This means that ISK functions cannot be constructed."
write(stdout, *) "However, through correct rotation of orbitals, all orbitals should be " &
& //"made real, and the code run in real mode (with tRotatedOrbsReal set)."
write(stdout, *) "Alternatively, run again in complex mode without ISK functions."
write(stdout, *) "If ISK functions are desired, the kpoint lattice must be shifted from this position."
! call stop_all(this_routine,"All kpoints are self-inverse")
end if
write(stdout, *) "All inverse kpoint orbitals correctly assigned."
write(stdout, *) "Orbital Inverse Orbital"
do i = 1, nBasis
write(stdout, *) i, KPntInvSymOrb(i)
end do
end if
!SymLabelList2 and SymLabelCounts2 are now organised differently, so that it is more efficient, and easier to add new symmetries.
!SymLabelCounts is of size (2,ScratchSize), where 1,x gives the index in SymlabelList2 where the orbitals of symmetry x start.
!SymLabelCounts(2,x) tells us the number of orbitals of spin & symmetry x there are.
!Therefore, if you want to run over all orbitals of a specific symmetry, you want to run over
!SymLabelList from SymLabelCounts(1,sym) to SymLabelCounts(1,sym)+SymLabelCounts(2,sym)-1
allocate(SymLabelList2(nBasis))
if (allocated(sym_label_list_spat)) deallocate(sym_label_list_spat)
allocate(sym_label_list_spat(nBasis))
allocate(SymLabelCounts2(2, ScratchSize))
SymLabelList2(:) = 0 !Indices: spin-orbital number
SymLabelCounts2(:, :) = 0 !Indices: index/Number , symmetry(inc. spin)
allocate(Temp(ScratchSize))
do j = 1, nBasis
IF (G1(j)%Ms == 1) THEN
Spin = 1
ELSE
Spin = 2
end if
! write(stdout,*) "BASIS FN ",j,G1(j)%Sym,SymClasses((j+1)/2)
SymInd = ClassCountInd(Spin, SpinOrbSymLabel(j), G1(j)%Ml)
SymLabelCounts2(2, SymInd) = SymLabelCounts2(2, SymInd) + 1
end do
SymLabelCounts2(1, 1) = 1
do j = 2, ScratchSize
SymLabelCounts2(1, j) = SymLabelCounts2(1, j - 1) + SymLabelCounts2(2, j - 1)
end do
Temp(:) = SymLabelCounts2(1, :)
do j = 1, nBasis
IF (G1(j)%Ms == 1) THEN
Spin = 1
ELSE
Spin = 2
end if
SymInd = ClassCountInd(Spin, SpinOrbSymLabel(j), G1(j)%Ml)
SymLabelList2(Temp(SymInd)) = j
Temp(SymInd) = Temp(SymInd) + 1
end do
! get also only the spatial orbitals
sym_label_list_spat = gtID(SymLabelList2)
! write(stdout,*) "SymLabelCounts2: ",SymLabelCounts2(1,:)
! write(stdout,*) "SymLabelCounts2: ",SymLabelCounts2(2,:)
Deallocate(Temp)
allocate(OrbClassCount(ScratchSize))
OrbClassCount(:) = 0
IF (tNoSymGenRandExcits .or. tUEG) THEN
!All orbitals are in irrep 0
OrbClassCount(ClassCountInd(1, 0, 0)) = (nBasis / 2)
OrbClassCount(ClassCountInd(2, 0, 0)) = (nBasis / 2)
ELSE
do i = 1, nBasis
IF (G1(i)%Ms == 1) THEN
! write(stdout,*) "Index: ",ClassCountInd(1,SpinOrbSymLabel(i),G1(i)%Ml)
! write(stdout,*) i,"SpinOrbSymLabel: ",SpinOrbSymLabel(i)
OrbClassCount(ClassCountInd(1, SpinOrbSymLabel(i), G1(i)%Ml)) = &
& OrbClassCount(ClassCountInd(1, SpinOrbSymLabel(i), G1(i)%Ml)) + 1
ELSE
! write(stdout,*) "Index: ",ClassCountInd(1,SpinOrbSymLabel(i),G1(i)%Ml)
! write(stdout,*) i,"SpinOrbSymLabel: ",SpinOrbSymLabel(i)
OrbClassCount(ClassCountInd(2, SpinOrbSymLabel(i), G1(i)%Ml)) = &
& OrbClassCount(ClassCountInd(2, SpinOrbSymLabel(i), G1(i)%Ml)) + 1
end if
end do
end if
! write(stdout,*) "*******",OrbClassCount(:)
! ELSE
!!SymLabelCounts(2,1:nSymLabels) gives the number of *states* in each symmetry class.
!!There are therefore equal number of alpha and beta orbitals in each state from which to calculate the unoccupied classcount.
! do i=1,nSymLabels
! OrbClassCount(ClassCountInd(1,i-1,0))=SymLabelCounts2(1,2,i)
! OrbClassCount(ClassCountInd(2,i-1,0))=SymLabelCounts2(2,2,i)
! end do
! end if
! This makes a 3D lookup table kPointToBasisFn(kx,ky,1,ms_index) which gives the orbital number for a given kx, ky and ms_index
IF (tHub .and. .not. (NMAXZ /= 0 .and. NMAXZ /= 1)) THEN
! IF(NMAXZ.ne.0.and.NMAXZ.ne.1) CALL Stop_All("SpinOrbSymSetup","This routine doesn't work with non-2D Hubbard model")
kmaxX = 0
kminX = 0
kmaxY = 0
kminY = 0
! [W.D:] can we make this the same as in the UEG:
kminZ = 0
kmaxZ = 0
do i = 1, nBasis ! In the hubbard model with tilted lattice boundary conditions, it's unobvious what the maximum values of
! kx and ky are, so this should be found
IF (G1(i)%k(1) > kmaxX) kmaxX = G1(i)%k(1)
IF (G1(i)%k(1) < kminX) kminX = G1(i)%k(1)
IF (G1(i)%k(2) > kmaxY) kmaxY = G1(i)%k(2)
IF (G1(i)%k(2) < kminY) kminY = G1(i)%k(2)
if (G1(i)%k(3) > kmaxz) kmaxz = g1(i)%k(3)
if (G1(i)%k(3) < kminz) kminz = g1(i)%k(3)
end do
! allocate(kPointToBasisFn(kminX:kmaxX,kminY:kmaxY,1,2))
allocate(kPointToBasisFn(kminX:kmaxX, kminY:kmaxY, kminz:kmaxz, 2))
kPointToBasisFn = -1 !Init to invalid
do i = 1, nBasis
iSpinIndex = (G1(i)%Ms + 1) / 2 + 1 ! iSpinIndex equals 1 for a beta spin (ms=-1), and 2 for an alpha spin (ms=1)
kPointToBasisFn(G1(i)%k(1), G1(i)%k(2), G1(i)%k(3), iSpinIndex) = i
end do
end if
!======================================================
if (tUEG2) then
kmaxX = 0
kminX = 0
kmaxY = 0
kminY = 0
kminZ = 0
kmaxZ = 0
do i = 1, nBasis
IF (kvec(i, 1) > kmaxX) kmaxX = kvec(i, 1)
IF (kvec(i, 1) < kminX) kminX = kvec(i, 1)
IF (kvec(i, 2) > kmaxY) kmaxY = kvec(i, 2)
IF (kvec(i, 2) < kminY) kminY = kvec(i, 2)
IF (kvec(i, 3) > kmaxZ) kmaxZ = kvec(i, 3)
IF (kvec(i, 3) < kminZ) kminZ = kvec(i, 3)
end do
allocate(kPointToBasisFn(kminX:kmaxX, kminY:kmaxY, kminZ:kmaxZ, 2))
kPointToBasisFn = -1 !Init to invalid
do i = 1, nBasis
iSpinIndex = (G1(i)%Ms + 1) / 2 + 1 ! iSpinIndex equals 1 for a beta spin (ms=-1), and 2 for an alpha spin (ms=1)
kPointToBasisFn(kvec(i, 1), kvec(i, 2), kvec(i, 3), iSpinIndex) = i
end do
kTotal(1) = 0
kTotal(2) = 0
kTotal(3) = 0
do j = 1, NEl
kTotal(1) = kTotal(1) + kvec(FDet(j), 1)
kTotal(2) = kTotal(2) + kvec(FDet(j), 2)
kTotal(3) = kTotal(3) + kvec(FDet(j), 3)
end do
write(stdout, *) "Total momentum is", kTotal
return
end if
!======================================================
! This makes a 3D lookup table kPointToBasisFn(kx,ky,kz,ms_index) which
!gives the orbital number for a given kx, ky, kz and ms_index
IF (tUEG) THEN
kmaxX = 0
kminX = 0
kmaxY = 0
kminY = 0
kminZ = 0
kmaxZ = 0
do i = 1, nBasis
IF (G1(i)%k(1) > kmaxX) kmaxX = G1(i)%k(1)
IF (G1(i)%k(1) < kminX) kminX = G1(i)%k(1)
IF (G1(i)%k(2) > kmaxY) kmaxY = G1(i)%k(2)
IF (G1(i)%k(2) < kminY) kminY = G1(i)%k(2)
IF (G1(i)%k(3) > kmaxZ) kmaxZ = G1(i)%k(3)
IF (G1(i)%k(3) < kminZ) kminZ = G1(i)%k(3)
end do
allocate(kPointToBasisFn(kminX:kmaxX, kminY:kmaxY, kminZ:kmaxZ, 2))
kPointToBasisFn = -1 !Init to invalid
do i = 1, nBasis
iSpinIndex = (G1(i)%Ms + 1) / 2 + 1 ! iSpinIndex equals 1 for a beta spin (ms=-1), and 2 for an alpha spin (ms=1)
kPointToBasisFn(G1(i)%k(1), G1(i)%k(2), G1(i)%k(3), iSpinIndex) = i
end do
end if
IF (tUEG .or. tHUB) THEN
kTotal = 0
do j = 1, NEl
if (t_k_space_hubbard) then
kTotal = lat%add_k_vec(kTotal, G1(FDet(j))%k)
else
kTotal = kTotal + G1(FDet(j))%k
end if
! just to be sure..
! not necessary anymore with new add_k_vec
! if (t_k_space_hubbard) then
! ktotal = lat%map_k_vec(ktotal)
! end if
end do
if (tHub .and. .not. t_k_space_hubbard) then
! is this turned off correctly?! check:
ktrial = (/kTotal(1), kTotal(2), 0/)
CALL MomPbcSym(ktrial, nBasisMax)
kTotal(1) = ktrial(1)
kTotal(2) = ktrial(2)
end if
write(stdout, *) "Total momentum is", kTotal
end if
END SUBROUTINE SpinOrbSymSetup
logical function IsMomAllowedDet(nJ)
integer :: nJ(nEl)
IsMomAllowedDet = IsMomAllowedDetAnyParent(nJ, HFSym%Sym)
end function IsMomAllowedDet
LOGICAL FUNCTION IsMomAllowedDetAnyParent(nJ, parentSym)
Type(Symmetry) :: SYM1
Type(Symmetry) :: parentSym
Type(BasisFN) :: iSym
INTEGER :: i, nJ(NEl), KPnt(3)
IsMomAllowedDetAnyParent = .false.
SYM1%S = 0
do i = 1, NEl
SYM1 = SYMPROD(SYM1, G1(nJ(i))%Sym)
end do
IF (SYM1%S /= parentSym%S) THEN
write(stdout, *) "nJ: ", nJ(:)
write(stdout, *) "parentSym,SYM1: ", parentSym%S, SYM1%S
! write(stdout,*) "Counter: ",Counter
CALL DecomposeAbelianSym(SYM1%S, KPnt)
write(stdout, "(A,3I5)") "KPnt for nJ: ", KPnt(1), KPnt(2), KPnt(3)
CALL Stop_All("IsMomAllowedDet", "Momentum forbidden excitation created1.")
ELSE
IsMomAllowedDetAnyParent = .true.
end if
CALL GETSYM(nJ, NEl, G1, nBasisMax, iSym)
IF (iSym%Sym%S /= parentSym%S) THEN
write(stdout, *) "nJ: ", nJ(:)
write(stdout, *) "parentSym,SYM1: ", parentSym%S, iSym%Sym%S
! write(stdout,*) "Counter: ",Counter
CALL DecomposeAbelianSym(iSym%Sym%S, KPnt)
write(stdout, "(A,3I5)") "KPnt for nJ: ", KPnt(1), KPnt(2), KPnt(3)
CALL Stop_All("IsMomAllowedDet", "Momentum forbidden excitation created2.")
ELSE
IsMomAllowedDetAnyParent = .true.
end if
END FUNCTION IsMomAllowedDetAnyParent
END MODULE GenRandSymExcitNUMod