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