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