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