SUBROUTINE FindSymMCSizeExcitLevel(IUNIT)
use SymData, only: TwoCycleSymGens
use SystemData, only: nEl, G1, nBasis, iMCCalcTruncLev
use SystemData, only: tUEG, tHPHF, tHub
use SystemData, only: CalcDetCycles, CalcDetPrint, tFixLz
use DeterminantData, only: FDet
use dSFMT_interface
use soft_exit, only: ChangeVars, tSoftExitFound
use Parallel_neci
use DetBitops, only: EncodeBitDet
use bit_rep_data, only: NIfTot
IMPLICIT NONE
INTEGER, intent(in) :: IUNIT
INTEGER :: j, SpatOrbs, FDetMom, Attempts, iExcitLevTest, i, ExcitLev
INTEGER(KIND=n_int) :: FDetiLut(0:NIfTot), iLut(0:NIfTot)
INTEGER :: FDetSym, TotalSym, TotalMom, alpha, beta, ierr, Momx, Momy
integer(int64) :: Accept, AcceptAll, TotalAttempts, TotalAttemptsAll
integer(int64) :: ExcitBin(0:iMCCalcTruncLev), ExcitBinAll(0:iMCCalcTruncLev)
real(dp) :: ExcitLevBias(0:iMCCalcTruncLev)
real(dp) :: FullSpace, r, Frac, SymSpace
real(dp) :: SizeLevel(0:iMCCalcTruncLev)
LOGICAL :: tDummy, tDummy2
iExcitLevTest = iMCCalcTruncLev
write(IUNIT, "(A,I6)") "Calculating MC size of symmetry-allowed " &
//"space of excitation levels up to level: ", iExcitLevTest
write(IUNIT, "(I18,A,I18,A)") CalcDetCycles, " MC cycles will be used, and " &
//"statistics printed out every ", CalcDetPrint, " cycles."
FDetSym = 0
FDetMom = 0
ExcitBin(:) = 0
ExcitBinAll(:) = 0
do i = 1, NEl
FDetSym = IEOR(FDetSym, int(G1(FDet(i))%Sym%S))
IF (tFixLz) FDetMom = FDetMom + G1(FDet(i))%Ml
end do
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 (tHPHF) THEN
write(stdout, *) "Imposing time-reversal symmetry (HPHF) on " &
//"size of space calculation"
end if
Accept = 0
AcceptAll = 0
TotalAttemptsAll = 0
TotalAttempts = 0 !This is the total number of attempts at creating a sym-allowed det (including successful ones)
!Sz symmetry could be put in here to make it more efficient
!(would be a little fiddly for OS systems though)
FullSpace = choose_i64(NEl, iExcitLevTest) !Pick 4 holes
FullSpace = FullSpace * choose_i64(nBasis - NEl + iExcitLevTest, iExcitLevTest) !Pick total unoccupied space
!Calculate excitation level bias due to the way the determinants are constructed.
do i = 0, iExcitLevTest
ExcitLevBias(i) = choose_i64(NEl - i, iExcitLevTest - i)
! write(stdout,*) ExcitLevBias(i)
end do
write(IUNIT, *) "Size of excitation level neglecting all symmetry: " &
, FullSpace
CALL neci_flush(IUNIT)
IF (iProcIndex == 0) THEN
OPEN (14, file="TruncSpaceMCStats", 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, int(CalcDetCycles)
!Create a random determinant up to excitation level iExcitLevTest from FDetiLut
!Returns det (iLut) and its excitation level, ExcitLevel, and the number of attempts
!needed to generate the symmetry allowed determinant.
CALL CreateRandomExcitLevDet(iExcitLevTest, FDet, FDetiLut, iLut, ExcitLev, Attempts)
TotalAttempts = TotalAttempts + Attempts
Accept = Accept + 1
!Add to correct bin for the excitation level
ExcitBin(ExcitLev) = ExcitBin(ExcitLev) + 1
IF (mod(i, int(CalcDetPrint)) == 0) THEN
!Write out statistics
call MPIReduce(Accept, MPI_SUM, AcceptAll)
call MPIReduce(TotalAttempts, MPI_SUM, TotalAttemptsAll)
call MPIReduce(ExcitBin(0:iExcitLevTest), MPI_SUM, ExcitBinAll(0:iExcitLevTest))
SymSpace = 0.0_dp
Frac = REAL(AcceptAll, dp) / REAL(TotalAttemptsAll, dp) !Fraction of the 'full' space which is symmetry allowed
do j = 0, iExcitLevTest
! write(stdout,*) REAL(ExcitBinAll(j),dp),REAL(AcceptAll,dp),Frac,FullSpace,ExcitLevBias(j)
SizeLevel(j) = ((REAL(ExcitBinAll(j), dp) / REAL(AcceptAll, dp)) * Frac * FullSpace) / ExcitLevBias(j)
SymSpace = SymSpace + SizeLevel(j)
end do
IF (iProcIndex == 0) THEN
write(14, "(2I16,2G35.15)", advance='no') i, AcceptAll, Frac, SymSpace
do j = 0, iExcitLevTest
write(14, "(F30.5)", advance='no') SizeLevel(j)
end do
write(14, "(A)") ""
end if
AcceptAll = 0
ExcitBinAll(0:iExcitLevTest) = 0
TotalAttemptsAll = 0
call ChangeVars(tDummy, tDummy2)
IF (tSoftExitFound) EXIT
end if
end do
call MPIReduce(Accept, MPI_SUM, AcceptAll)
call MPIReduce(TotalAttempts, MPI_SUM, TotalAttemptsAll)
call MPIReduce(ExcitBin(0:iExcitLevTest), MPI_SUM, ExcitBinAll(0:iExcitLevTest))
SymSpace = 0.0_dp
Frac = REAL(AcceptAll, dp) / REAL(TotalAttemptsAll, dp) !Fraction of the 'full' space which is symmetry allowed
do j = 0, iExcitLevTest
! write(stdout,*) REAL(ExcitBinAll(j),dp),REAL(AcceptAll,dp),Frac,FullSpace,ExcitLevBias(j)
SizeLevel(j) = ((REAL(ExcitBinAll(j), dp) / REAL(AcceptAll, dp)) * Frac * FullSpace) / ExcitLevBias(j)
SymSpace = SymSpace + SizeLevel(j)
end do
IF (iProcIndex == 0) THEN
write(14, "(2I16,2G35.15)", advance='no') i, AcceptAll, Frac, SymSpace
do j = 0, iExcitLevTest
write(14, "(F30.5)", advance='no') SizeLevel(j)
end do
write(14, "(A)") ""
close(14)
end if
write(IUNIT, *) "MC size of truncated space: ", SymSpace
write(IUNIT, *) "Individual excitation level contributions: "
do j = 0, iExcitLevTest
write(IUNIT, "(I5,F30.5)") j, SizeLevel(j)
end do
CALL neci_flush(IUNIT)
END SUBROUTINE FindSymMCSizeExcitLevel