Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|

integer | :: | nI(NEl) |
||||

integer(kind=n_int) | :: | iLut(0:NIfTot) |
||||

integer | :: | Elec1Ind |
||||

integer | :: | Elec2Ind |
||||

integer | :: | SymProduct |
||||

integer | :: | iSpn |
||||

integer | :: | OrbA |
||||

integer | :: | OrbB |
||||

integer | :: | nParts |
||||

integer | :: | iCreate |
||||

real(kind=dp) | :: | Tau |

```
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
```

NECI was developed by Alavi Group

© 2023

Documentation generated by FORD