SUBROUTINE FindSymMCSizeofSpace(IUNIT)
use SymData, only: TwoCycleSymGens
use SystemData, only: nEl, G1, nBasis, nOccAlpha, nOccBeta, tMolpro
use SystemData, only: tUEG, tHPHF, tHub, tKPntSym, Symmetry
use SystemData, only: CalcDetCycles, CalcDetPrint, tFixLz
use CalcData, only: tTruncNOpen, trunc_nopen_max, &
t_trunc_nopen_diff, trunc_nopen_diff
use DeterminantData, only: FDet
use DetCalcData, only: ICILevel
use dSFMT_interface
use soft_exit, only: ChangeVars, tSoftExitFound
use Parallel_neci
use DetBitops, only: EncodeBitDet, IsAllowedHPHF, count_open_orbs
use bit_rep_data, only: NIfTot
use sym_mod, only: SymProd
use fcimcdata, only: ilutRef
IMPLICIT NONE
INTEGER :: IUNIT, j, SpatOrbs, FDetMom, ExcitLev
INTEGER(KIND=n_int) :: FDetiLut(0:NIfTot), iLut(0:NIfTot)
INTEGER :: FDetSym, TotalSym, TotalMom, alpha, beta, ierr, Momx, Momy
INTEGER :: Momz, nopenorbs, Space_unit
integer(int64) :: Accept, AcceptAll, i
integer(int64) :: ExcitBin(0:NEl), ExcitBinAll(0:NEl)
real(dp) :: FullSpace, r, Frac
real(dp) :: SizeLevel(0:NEl)
LOGICAL :: truncate_space, tDummy, tDummy2
LOGICAL :: tNotAllowed, tAcc
type(Symmetry) :: FDetKPntMom, KPntMom
write(IUNIT, *) "Calculating exact size of symmetry-allowed " &
//"determinant space using MC..."
write(IUNIT, *) CalcDetCycles, " MC cycles will be used, and " &
//"statistics printed out every ", CalcDetPrint, " cycles."
FDetSym = 0
FDetMom = 0
ExcitBin(:) = 0
ExcitBinAll(:) = 0
Momx = 0
Momy = 0
Momz = 0
FDetKPntMom%S = 0
do i = 1, NEl
IF (tFixLz) FDetMom = FDetMom + G1(FDet(i))%Ml
if (tKPntSym) FDetKPntMom = Symprod(FDetKPntMom, G1(FDet(i))%Sym)
IF (tUEG .or. tHub) THEN
Momx = Momx + G1(FDet(i))%k(1)
Momy = Momy + G1(FDet(i))%k(2)
Momz = Momz + G1(FDet(i))%k(3)
else
FDetSym = IEOR(FDetSym, int(G1(FDet(i))%Sym%S))
end if
end do
IF (.not. ((Momx == 0) .and. (Momy == 0) .and. (Momz == 0))) THEN
write(stdout, *) "Momx: ", Momx, "Momy: ", Momy, "Momz: ", Momz
call stop_all("FindSymMCSizeofSpace", "Cannot calculate MC size of space with non-zero momentum")
end if
IF (ICILevel > 0) THEN
truncate_space = .true.
ELSE
truncate_space = .false.
end if
CALL EncodeBitDet(FDet, FDetiLut)
write(IUNIT, *) "Symmetry of HF determinant is: ", FDetSym
IF (tFixLz) THEN
write(IUNIT, *) "Imposing momentum sym on size calculation"
write(IUNIT, *) "Momentum of HF determinant is: ", FDetMom
end if
if (tKPntSym) then
write(iunit, *) "Using k-point symmetry."
write(iunit, *) "K-point sym of HF determinant is: ", FDetKPntMom%S
end if
IF (tHPHF) THEN
write(stdout, *) "Imposing time-reversal symmetry (HPHF) on " &
//"size of space calculation"
end if
SpatOrbs = nBasis / 2
Accept = 0
FullSpace = choose_i64(SpatOrbs, nOccAlpha)
FullSpace = FullSpace * choose_i64(SpatOrbs, nOccBeta)
write(IUNIT, *) "Size of space neglecting all but Sz symmetry: " &
, FullSpace
CALL neci_flush(IUNIT)
IF (iProcIndex == 0) THEN
Space_unit = get_free_unit()
OPEN (Space_unit, file="SpaceMCStats", status='unknown', &
form='formatted')
end if
! With MerTwistRan the default seed was being used.
! dSFMT does not initialise itself if not already initialised.
call dSFMT_init(5489)
do i = 1, CalcDetCycles
KPntMom%S = 0
TotalSym = 0
TotalMom = 0
ExcitLev = 0
Momx = 0
Momy = 0
Momz = 0
iLut(:) = 0
!Create random determinant (Correct Sz symmetry)
!Loop over alpha electrons
do j = 1, nOccAlpha
tNotAllowed = .true.
do while (tNotAllowed)
r = genrand_real2_dSFMT()
alpha = 2 * (INT(SpatOrbs * r) + 1)
IF (.not. BTEST(iLut((alpha - 1) / bits_n_int), mod((alpha - 1), bits_n_int))) THEN
!Has *not* been picked before
iLut((alpha - 1) / bits_n_int) = &
IBSET(iLut((alpha - 1) / bits_n_int), mod(alpha - 1, bits_n_int))
tNotAllowed = .false.
end if
end do
IF (tFixLz) THEN
TotalMom = TotalMom + G1(alpha)%Ml
end if
IF (tUEG .or. tHub) THEN
Momx = Momx + G1(alpha)%k(1)
Momy = Momy + G1(alpha)%k(2)
Momz = Momz + G1(alpha)%k(3)
else if (tKPntSym) then
KPntMom = Symprod(KPntMom, G1(alpha)%Sym)
else
TotalSym = IEOR(TotalSym, int((G1(alpha)%Sym%S)))
end if
IF (.not. BTEST(FDetiLut((alpha - 1) / bits_n_int), mod((alpha - 1), bits_n_int))) THEN
!orbital chosen is *not* in the reference determinant
ExcitLev = ExcitLev + 1
end if
end do
!Loop over beta electrons
do j = 1, nOccBeta
tNotAllowed = .true.
do while (tNotAllowed)
r = genrand_real2_dSFMT()
beta = 2 * (INT(SpatOrbs * r) + 1) - 1
IF (.not. BTEST(iLut((beta - 1) / bits_n_int), mod((beta - 1), bits_n_int))) THEN
!Has *not* been picked before
iLut((beta - 1) / bits_n_int) = &
IBSET(iLut((beta - 1) / bits_n_int), mod(beta - 1, bits_n_int))
tNotAllowed = .false.
end if
end do
IF (tFixLz) THEN
TotalMom = TotalMom + G1(beta)%Ml
end if
IF (tUEG .or. tHub) THEN
Momx = Momx + G1(beta)%k(1)
Momy = Momy + G1(beta)%k(2)
Momz = Momz + G1(beta)%k(3)
else if (tKPntSym) then
KPntMom = Symprod(KPntMom, G1(beta)%Sym)
ELSE
TotalSym = IEOR(TotalSym, int((G1(beta)%Sym%S)))
end if
IF (.not. BTEST(FDetiLut((beta - 1) / bits_n_int), mod((beta - 1), bits_n_int))) THEN
!orbital chosen is *not* in the reference determinant
ExcitLev = ExcitLev + 1
end if
end do
if (tTruncNOpen .or. t_trunc_nopen_diff) nOpenOrbs = count_open_orbs(iLut)
tAcc = .true.
if (TotalSym /= FDetSym) then
tAcc = .false.
end if
if (tAcc .and. tHPHF) then
if (.not. IsAllowedHPHF(iLut)) then
tAcc = .false.
end if
end if
if (tAcc .and. tFixLz) then
if (TotalMom /= FDetMom) then
tAcc = .false.
end if
end if
if (tAcc .and. truncate_space) then
if (ExcitLev > ICILevel) then
tAcc = .false.
end if
end if
if (tAcc .and. (tUEG .or. tHub)) then
if ((Momx /= 0) .or. (Momy /= 0) .or. (Momz /= 0)) then
tAcc = .false.
end if
end if
if (tAcc .and. (tKPntSym)) then
if (KPntMom%S /= FDetKPntMom%S) then
tAcc = .false.
end if
end if
if (tAcc .and. tTruncNOpen) then
if (nOpenOrbs > trunc_nopen_max) then
tAcc = .false.
end if
end if
if (tAcc .and. t_trunc_nopen_diff) then
if (abs(nOpenOrbs - count_open_orbs(ilutRef(:, 1))) > trunc_nopen_diff) then
tAcc = .false.
end if
end if
if (tAcc) Accept = Accept + 1
IF (tAcc) THEN
!Add to correct bin for the excitation level
ExcitBin(ExcitLev) = ExcitBin(ExcitLev) + 1
end if
IF (mod(i, CalcDetPrint) == 0) THEN
!Write out statistics
call MPIAllReduce(Accept, MPI_SUM, AcceptAll)
call MPIAllReduce(ExcitBin(0:NEl), MPI_SUM, ExcitBinAll(0:NEl))
Frac = REAL(AcceptAll, dp) / REAL(i * nProcessors, dp)
do j = 0, NEl
SizeLevel(j) = (REAL(ExcitBinAll(j), dp) / REAL(AcceptAll, dp)) * Frac * FullSpace
end do
IF (iProcIndex == 0) THEN
if (tMolpro) then
write(Space_unit, "(I16,G35.15)") i, Frac * FullSpace
else
write(Space_unit, "(2I16,2G35.15)", advance='no') i, AcceptAll, Frac, Frac * FullSpace
do j = 0, NEl
write(Space_unit, "(F30.5)", advance='no') SizeLevel(j)
end do
write(Space_unit, "(A)") ""
end if
end if
AcceptAll = 0
ExcitBinAll(0:NEl) = 0
CALL ChangeVars(tDummy, tDummy2)
IF (tSoftExitFound) EXIT
end if
end do
call MPIAllReduce(Accept, MPI_SUM, AcceptAll)
call MPIAllReduce(ExcitBin(0:NEl), MPI_SUM, ExcitBinAll(0:NEl))
Frac = REAL(AcceptAll, dp) / REAL(i * nProcessors, dp)
do j = 0, NEl
SizeLevel(j) = (REAL(ExcitBinAll(j), dp) / REAL(AcceptAll, dp)) * Frac * FullSpace
end do
IF (iProcIndex == 0) THEN
if (tMolpro) then
write(Space_unit, "(I16,G35.15)") i, Frac * FullSpace
else
write(Space_unit, "(2I16,2G35.15)", advance='no') i, AcceptAll, Frac, Frac * FullSpace
do j = 0, NEl
write(Space_unit, "(F30.5)", advance='no') SizeLevel(j)
end do
write(Space_unit, "(A)") ""
end if
close(Space_unit)
end if
write(IUNIT, *) "MC size of space: ", Frac * FullSpace
write(IUNIT, *) "Individual excitation level contributions: "
do j = 0, NEl
write(IUNIT, "(I5,F30.5)") j, SizeLevel(j)
end do
CALL neci_flush(IUNIT)
END SUBROUTINE FindSymMCSizeofSpace