SUBROUTINE SetupParameters()
INTEGER :: ierr, i, j, HFDetTest(NEl), Seed, alpha, beta, symalpha, symbeta, endsymstate
INTEGER :: LargestOrb, nBits, HighEDet(NEl), orb
INTEGER(KIND=n_int) :: iLutTemp(0:NIfTot)
HElement_t(dp) :: TempHii
real(dp) :: UpperTau, r
CHARACTER(len=*), PARAMETER :: t_r = 'SetupParameters'
CHARACTER(*), parameter :: this_routine = t_r
CHARACTER(len=12) :: abstr
character(len=40) :: filename
#ifndef PROG_NUMRUNS_
character(len=40) :: filename2
#endif
LOGICAL :: tSuccess, tFoundOrbs(nBasis), FoundPair, tSwapped, tAlreadyOcc
INTEGER :: HFLz, ChosenOrb, step, SymFinal, run
integer(int64) :: SymHF
integer(n_int), allocatable :: dummy_list(:, :)
! CALL MPIInit(.false.) !Initialises MPI - now have variables iProcIndex and nProcessors
write(stdout, *)
if (nProcessors > 1) then
write(stdout, *) "Performing Parallel FCIQMC...."
else
write(stdout, *) "Performing single-core FCIQMC...."
end if
write(stdout, *)
!Set timed routine names
Walker_Time%timer_name = 'WalkerTime'
Annihil_Time%timer_name = 'AnnihilTime'
Sort_Time%timer_name = 'SortTime'
Comms_Time%timer_name = 'CommsTime'
ACF_Time%timer_name = 'ACFTime'
AnnSpawned_time%timer_name = 'AnnSpawnedTime'
AnnMain_time%timer_name = 'AnnMainTime'
BinSearch_time%timer_name = 'BinSearchTime'
Sync_Time%timer_name = 'SyncTime'
SemiStoch_Comms_Time%timer_name = 'SemiStochCommsTime'
SemiStoch_Multiply_Time%timer_name = 'SemiStochMultiplyTime'
Trial_Search_Time%timer_name = 'TrialSearchTime'
SemiStoch_Init_Time%timer_name = 'SemiStochInitTime'
SemiStoch_Hamil_Time%timer_name = 'SemiStochHamilTime'
SemiStoch_Davidson_Time%timer_name = 'SemiStochDavidsonTime'
Trial_Init_Time%timer_name = 'TrialInitTime'
InitSpace_Init_Time%timer_name = 'InitSpaceTime'
kp_generate_time%timer_name = 'KPGenerateTime'
Stats_Comms_Time%timer_name = 'StatsCommsTime'
subspace_hamil_time%timer_name = 'SubspaceHamilTime'
exact_subspace_h_time%timer_name = 'ExactSubspace_H_Time'
subspace_spin_time%timer_name = 'SubspaceSpinTime'
var_e_time%timer_name = 'VarEnergyTime'
precond_e_time%timer_name = 'PreCondEnergyTime'
proj_e_time%timer_name = 'ProjEnergyTime'
rescale_time%timer_name = 'RescaleTime'
death_time%timer_name = 'DeathTime'
hash_test_time%timer_name = 'HashTestTime'
hii_test_time%timer_name = 'HiiTestTime'
init_flag_time%timer_name = 'InitFlagTime'
GAS_PCHB_init_time%timer_name = 'GAS_PCHB_init_time'
! Initialise allocated arrays with input data
TargetGrowRate(:) = InputTargetGrowRate
TargetGrowRateWalk(:) = InputTargetGrowRateWalk
! Initialize
AllTotParts = 0.0_dp
AllTotPartsOld = 0.0_dp
AllTotPartsLastOutput = 0.0_dp
IF (TDebug) THEN
!This will open a file called LOCALPOPS-"iprocindex" on unit number 11 on every node.
abstr = 'LOCALPOPS-'//str(iProcIndex)
open(11, FILE=abstr, STATUS='UNKNOWN')
end if
IF (iProcIndex == Root) THEN
if (.not. tFCIMCStats2) then
fcimcstats_unit = get_free_unit()
if (tReadPops .and. .not. t_no_append_stats) then
! Restart calculation. Append to stats file (if it exists).
if (tMolpro .and. .not. tMolproMimic) then
filename = 'FCIQMCStats_'//adjustl(MolproID)
open(fcimcstats_unit, file=filename, status='unknown', position='append')
else
open(fcimcstats_unit, file='FCIMCStats', status='unknown', position='append')
end if
else
call MoveFCIMCStatsFiles() !This ensures that FCIMCStats files are not overwritten
if (tMolpro .and. .not. tMolproMimic) then
filename = 'FCIQMCStats_'//adjustl(MolproID)
open(fcimcstats_unit, file=filename, status='unknown')
else
open(fcimcstats_unit, file='FCIMCStats', status='unknown')
end if
end if
end if
#ifndef PROG_NUMRUNS_
if (inum_runs == 2) then
fcimcstats_unit2 = get_free_unit()
if (tReadPops .and. .not. t_no_append_stats) then
! Restart calculation. Append to stats file (if it exists).
if (tMolpro .and. .not. tMolproMimic) then
filename2 = 'FCIQMCStats2_'//adjustl(MolproID)
open(fcimcstats_unit2, file=filename2, status='unknown', position='append')
else
open(fcimcstats_unit2, file='FCIMCStats2', status='unknown', position='append')
end if
else
if (tMolpro .and. .not. tMolproMimic) then
filename2 = 'FCIQMCStats2_'//adjustl(MolproID)
open(fcimcstats_unit2, file=filename2, status='unknown')
else
open(fcimcstats_unit2, file='FCIMCStats2', status='unknown')
end if
end if
end if
#endif
IF (tTruncInitiator .and. (.not. tFCIMCStats2)) THEN
initiatorstats_unit = get_free_unit()
if (tReadPops .and. .not. t_no_append_stats) then
! Restart calculation. Append to stats file (if it exists)
open(initiatorstats_unit, file='INITIATORStats', status='unknown', form='formatted', position='append')
else
open(initiatorstats_unit, file='INITIATORStats', status='unknown', form='formatted')
end if
end if
IF (tLogComplexPops) THEN
ComplexStats_unit = get_free_unit()
open(ComplexStats_unit, file='COMPLEXStats', status='unknown')
end if
if (tLogEXLEVELStats) then
EXLEVELStats_unit = get_free_unit()
open (EXLEVELStats_unit, file='EXLEVELStats', &
status='unknown', position='append')
end if
end if
!Store information specifically for the HF determinant
allocate(HFDet(NEl), stat=ierr)
CALL LogMemAlloc('HFDet', NEl, 4, t_r, HFDetTag)
allocate(iLutHF(0:NIfTot), stat=ierr)
IF (ierr /= 0) CALL Stop_All(t_r, "Cannot allocate memory for iLutHF")
do i = 1, NEl
HFDet(i) = FDet(i)
end do
CALL EncodeBitDet(HFDet, iLutHF)
if (tHPHF .and. TestClosedShellDet(iLutHF) .and. tOddS_HPHF .and. TwoCycleSymGens) then
!This is not a compatible reference function.
!Create single excitation of the correct symmetry
!Use this as the reference.
write(stdout, "(A)") "Converging to ODD S eigenstate"
SymFinal = int((G1(HFDet(nel))%Sym%S)) + 1
tAlreadyOcc = .false.
do i = SymLabelCounts(1, SymFinal), SymLabelCounts(1, SymFinal) + SymLabelCounts(2, SymFinal) - 1
if (G1(HFDet(nel))%Ms == -1) then
!Choose beta ones
orb = (2 * SymLabelList(i)) - 1
else
orb = (2 * SymLabelList(i))
end if
tAlreadyOcc = .false.
do j = 1, nel
if (orb == HFDet(j)) then
tAlreadyOcc = .true.
exit
end if
end do
if (.not. tAlreadyOcc) then
HFDet(nel) = orb
call sort(HFDet)
exit
end if
end do
if (tAlreadyOcc) &
call stop_all(t_r, "Cannot automatically detect open-shell determinant for reference to use with odd S")
call EncodeBitDet(HFDet, iLutHF)
if (TestClosedShellDet(iLutHF)) &
call stop_all(t_r, "Fail to create open-shell determinant for reference to use with odd S")
write(stdout, *) "Reference determinant changed to the open-shell:"
call write_det(stdout, HFDet, .true.)
end if
!iLutRef is the reference determinant for the projected energy.
!Initially, it is chosen to be the same as the inputted reference determinant
call setup_adi()
allocate(iLutRef(0:NIfTot, inum_runs), stat=ierr)
ilutRef = 0
allocate(ProjEDet(NEl, inum_runs), stat=ierr)
IF (ierr /= 0) CALL Stop_All(t_r, "Cannot allocate memory for iLutRef")
! The reference / projected energy determinants are the same as the
! HF determinant.
call assign_reference_dets()
allocate(iLutHF_True(0:NIfTot), stat=ierr)
IF (ierr /= 0) CALL Stop_All(t_r, "Cannot allocate memory for iLutHF_True")
allocate(HFDet_True(NEl), stat=ierr)
IF (ierr /= 0) CALL Stop_All(t_r, "Cannot allocate memory for HFDet_True")
!RDM uses HFDet_True in some parts but ilutRef in others
!Sorry here we make them the same to avoid errors there.
!Let's hope that unkonwn parts of the code do not depend on HFDet_True
if (tRef_Not_HF) then
do i = 1, NEl
HFDet_True(i) = BRR(i)
end do
call sort(HFDet_True(1:NEl))
CALL EncodeBitDet(HFDet_True, iLutHF_True)
else
iLutHF_True = iLutHF
HFDet_True = HFDet
end if
if (tGUGA) then
do run = 1, inum_runs
call fill_csf_i(ilutRef(:, run), csf_ref(run))
end do
end if
if (tHPHF) then
allocate(RefDetFlip(NEl, inum_runs), &
ilutRefFlip(0:NifTot, inum_runs))
do run = 1, inum_runs
if (.not. TestClosedShellDet(ilutRef(:, run))) then
! If the reference determinant corresponds to an open shell
! HPHF, then we need to specify the paired determinant and
! mark that this needs to be considered in calculating
! the projected energy.
tSpinCoupProjE(run) = .true.
call ReturnAlphaOpenDet(ProjEDet(:, run), &
RefDetFlip(:, run), &
ilutRef(:, run), &
ilutRefFlip(:, run), &
.true., .true., tSwapped)
if (tSwapped) &
write(stdout, *) 'HPHF used, and open shell determinant &
&for run ', run, ' spin-flippd for &
&consistency.'
write(stdout, *) "Two *different* determinants contained in &
&initial HPHF for run ", run
write(stdout, *) "Projected energy will be calculated as &
&projection onto both of these."
else
tSpinCoupProjE(run) = .false.
end if
end do
else
tSpinCoupProjE(:) = .false.
end if
!Init hash shifting data
hash_iter = 0
IF (tFixLz) THEN
CALL GetLz(HFDet, NEl, HFLz)
write(stdout, "(A,I5)") "Ml value of reference determinant is: ", HFLz
IF (HFLz /= LzTot) THEN
CALL Stop_All("SetupParameters", "Chosen reference determinant does not have the " &
& //"same Lz value as indicated in the input.")
end if
end if
!Do a whole lot of tests to see if we can use Brillouins theorem or not.
IF (tBrillouinsDefault) CALL CheckforBrillouins()
!test the encoding of the HFdet to bit representation.
! Test that the bit operations are working correctly...
! TODO: Move this to using the extract_bit_det routines to test those
! too...
call decode_bit_det(HFDetTest, iLutHF)
do i = 1, NEl
IF (HFDetTest(i) /= HFDet(i)) THEN
write(stdout, *) "HFDet: ", HFDet(:)
write(stdout, *) "HFDetTest: ", HFDetTest(:)
CALL Stop_All(t_r, "HF Determinant incorrectly decoded.")
end if
end do
CALL LargestBitSet(iLutHF, NIfD, LargestOrb)
IF (LargestOrb /= hfdet(nel)) then
write(stdout, *) 'ilut HF', ilutHF
write(stdout, *) 'largest orb', largestorb
write(stdout, *) 'HFDet', hfdet
CALL Stop_All(t_r, "LargestBitSet FAIL")
end if
nBits = CountBits(iLutHF, NIfD, NEl)
IF (nBits /= NEl) THEN
CALL Stop_All(t_r, "CountBits FAIL")
end if
allocate(HighestPopDet(0:NIfTot, inum_runs), stat=ierr)
IF (ierr /= 0) CALL Stop_All(t_r, "Cannot allocate memory for HighestPopDet")
HighestPopDet(:, :) = 0
iHighestPop = 0
!Check that the symmetry routines have set the symmetry up correctly...
tSuccess = .true.
tFoundOrbs(:) = .false.
IF ((.not. tHub) .and. (.not. tUEG) .and. TwoCycleSymGens) THEN
do i = 1, nSymLabels
EndSymState = SymLabelCounts(1, i) + SymLabelCounts(2, i) - 1
do j = SymLabelCounts(1, i), EndSymState
Beta = (2 * SymLabelList(j)) - 1
Alpha = (2 * SymLabelList(j))
SymAlpha = int((G1(Alpha)%Sym%S))
SymBeta = int((G1(Beta)%Sym%S))
IF (.not. tFoundOrbs(Beta)) THEN
tFoundOrbs(Beta) = .true.
ELSE
CALL Stop_All("SetupParameters", "Orbital specified twice")
end if
IF (.not. tFoundOrbs(Alpha)) THEN
tFoundOrbs(Alpha) = .true.
ELSE
CALL Stop_All("SetupParameters", "Orbital specified twice")
end if
IF (G1(Beta)%Ms /= -1) THEN
tSuccess = .false.
else if (G1(Alpha)%Ms /= 1) THEN
tSuccess = .false.
else if ((SymAlpha /= (i - 1)) .or. (SymBeta /= (i - 1))) THEN
tSuccess = .false.
end if
end do
end do
do i = 1, nBasis
IF (.not. tFoundOrbs(i)) THEN
write(stderr, *) "Orbital: ", i, " not found."
CALL Stop_All("SetupParameters", "Orbital not found")
end if
end do
end if
IF (.not. tSuccess) THEN
write(stderr, *) "************************************************"
write(stderr, *) "** WARNING!!! **"
write(stderr, *) "************************************************"
write(stderr, *) "Symmetry information of orbitals not the same in alpha and beta pairs."
write(stderr, *) "Symmetry now set up in terms of spin orbitals"
write(stderr, *) "I strongly suggest you check that the reference energy is correct."
end if
! From now on, the orbitals are contained in symlabellist2 and
! symlabelcounts2 rather than the original arrays. These are stored
! using spin orbitals. Assume that if we want to use the non-uniform
! random excitation generator, we also want to use the NoSpinSym full
! excitation generators if they are needed.
CALL GetSym(HFDet, NEl, G1, NBasisMax, HFSym)
Sym_Psi = int(HFSym%Sym%S) !Store the symmetry of the wavefunction for later
write(stdout, "(A,I10)") "Symmetry of reference determinant is: ", int(HFSym%Sym%S)
if (TwoCycleSymGens) then
SymHF = 0
do i = 1, NEl
SymHF = IEOR(SymHF, G1(HFDet(i))%Sym%S)
end do
write(stdout, "(A,I10)") "Symmetry of reference determinant from spin orbital symmetry info is: ", SymHF
if (SymHF /= HFSym%Sym%S) then
call stop_all(t_r, "Inconsistency in the symmetry arrays.")
end if
end if
!If using a CAS space truncation, write out this CAS space
IF (tTruncCAS) THEN
write(stdout, *) "Truncated CAS space detected. Writing out CAS space..."
write(stdout, '(A,I2,A,I2,A)') " In CAS notation, (spatial orbitals, electrons), this has been chosen as: (", &
(OccCASOrbs + VirtCASOrbs) / 2, ",", OccCASOrbs, ")"
do I = NEl - OccCASorbs + 1, NEl
write(stdout, '(6I7)', advance='no') I, BRR(I), G1(BRR(I))%K(1), G1(BRR(I))%K(2), G1(BRR(I))%K(3), G1(BRR(I))%MS
CALL WRITESYM(stdout, G1(BRR(I))%SYM, .FALSE.)
write(stdout, '(I4)', advance='no') G1(BRR(I))%Ml
write(stdout, '(2F19.9)') ARR(I, 1), ARR(BRR(I), 2)
end do
write(stdout, '(A)') " ================================================================================================="
do I = NEl + 1, NEl + VirtCASOrbs
write(stdout, '(6I7)', advance='no') I, BRR(I), G1(BRR(I))%K(1), G1(BRR(I))%K(2), G1(BRR(I))%K(3), G1(BRR(I))%MS
CALL WRITESYM(stdout, G1(BRR(I))%SYM, .FALSE.)
write(stdout, '(I4)', advance='no') G1(BRR(I))%Ml
write(stdout, '(2F19.9)') ARR(I, 1), ARR(BRR(I), 2)
end do
else if (tTruncInitiator) THEN
write(stdout, '(A)') "*********** INITIATOR METHOD IN USE ***********"
write(stdout, '(A)') "Starting with only the reference determinant in the fixed initiator space."
end if
! Setup excitation generator for the HF determinant. If we are using
! assumed sized excitgens, this will also be assumed size.
IF (tUEG .or. tHub .or. tNoSingExcits) THEN
exflag = 2
ELSE
exflag = 3
end if
IF (.not. tKPntSym) THEN
!Count all possible excitations - put into HFConn
!TODO: Get CountExcitations3 working with tKPntSym
CALL CountExcitations3(HFDet, exflag, nSingles, nDoubles)
ELSE
if (t_k_space_hubbard) then
! use my gen_all_excits_k_space_hubbard routine from the
! unit tests.. but i might have to set up some other stuff
! for this to work and also make sure this works with my
! new symmetry implementation
if (.not. t_trans_corr_2body) then
call gen_all_excits_k_space_hubbard(HFDet, nDoubles, dummy_list)
end if
nSingles = 0
else
! Use Alex's old excitation generators to enumerate all excitations.
call enumerate_sing_doub_kpnt(exflag, .false., nSingles, nDoubles, .false.)
end if
end if
HFConn = nSingles + nDoubles
if (tAutoAdaptiveShift .and. AAS_Cut < 0.0) then
!The user did not specify the value, use this as a default
if (HFConn > 0) then
AAS_Cut = 1.0_dp / HFConn
else
! if the HF is disconnected (can happen in rare corner cases), set it to 0
AAS_Cut = 0.0_dp
end if
end if
! Initialise random number seed - since the seeds need to be different
! on different processors, subract processor rank from random number
if (.not. tRestart) then
Seed = abs(G_VMC_Seed - iProcIndex)
write(stdout, "(A,I12)") "Value for seed is: ", Seed
!Initialise...
CALL dSFMT_init(Seed)
if (tMolpro) then
if ((NMCyc == -1) .and. (.not. tTimeExit)) then
!No iteration number, or TIME option has been specified.
call warning_neci(t_r, &
"No iteration number specified. Only running for 100 iterations initially. Change with ITERATIONS option.")
NMCyc = 100 !Only run for 100 iterations.
else if (tTimeExit .and. (NMCyc == -1)) then
write(stdout, "(A,F10.3,A)") "Running FCIQMC for ", MaxTimeExit / 60.0_dp, " minutes."
else if (tTimeExit .and. (NMCyc /= -1)) then
write(stdout, "(A,F10.3,A,I15,A)") "Running FCIQMC for ", MaxTimeExit / 60.0_dp, " minutes OR ", NMCyc, " iterations."
else if ((.not. tTimeExit) .and. (NMCyc > 0)) then
write(stdout, "(A,I15,A)") "Running FCIQMC for ", NMCyc, " iterations."
else
call stop_all(t_r, "Iteration number/Time unknown for simulation - contact ghb")
end if
end if
end if
! Option tRandomiseHashOrbs has now been removed.
! Its behaviour is now considered default
! --> Create a random mapping for the orbitals
allocate(RandomOrbIndex(nBasis), stat=ierr)
IF (ierr /= 0) THEN
CALL Stop_All(t_r, "Error in allocating RandomOrbIndex")
end if
RandomOrbIndex(:) = 0
! We want another independent randomizing array for the hash table, so
! we do not introduce correlations between the two
allocate(RandomHash2(nBasis), stat=ierr)
if (ierr /= 0) &
call stop_all(t_r, "Error in allocating RandomHash2")
RandomHash2(:) = 0
IF (iProcIndex == root) THEN
do i = 1, nBasis
! If we want to hash only by spatial orbitals, then the
! spin paired orbitals must be set equal
if (tSpatialOnlyHash) then
if (.not. btest(i, 0)) then
RandomOrbIndex(i) = RandomOrbIndex(i - 1)
cycle
end if
end if
! Ensure that we don't set two values to be equal accidentally
FoundPair = .false.
do while (.not. FoundPair)
r = genrand_real2_dSFMT()
ChosenOrb = INT(nBasis * r * 1000) + 1
! Check all values which have already been set.
do j = 1, nBasis
IF (RandomOrbIndex(j) == ChosenOrb) EXIT
end do
! If not already used, then we can move on
if (j == nBasis + 1) FoundPair = .true.
RandomOrbIndex(i) = ChosenOrb
end do
end do
!Do again for RandomHash2
do i = 1, nBasis
! If we want to hash only by spatial orbitals, then the
! spin paired orbitals must be set equal
if (tSpatialOnlyHash) then
if (.not. btest(i, 0)) then
RandomHash2(i) = RandomHash2(i - 1)
cycle
end if
end if
! Ensure that we don't set two values to be equal accidentally
FoundPair = .false.
do while (.not. FoundPair)
r = genrand_real2_dSFMT()
ChosenOrb = INT(nBasis * r * 1000) + 1
! Check all values which have already been set.
do j = 1, nBasis
IF (RandomHash2(j) == ChosenOrb) EXIT
end do
! If not already used, then we can move on
if (j == nBasis + 1) FoundPair = .true.
RandomHash2(i) = ChosenOrb
end do
end do
if (tSpatialOnlyHash) then
step = 2
else
step = 1
end if
do i = 1, nBasis
IF ((RandomOrbIndex(i) == 0) .or. (RandomOrbIndex(i) > nBasis * 1000)) THEN
CALL Stop_All(t_r, "Random Hash incorrectly calculated")
end if
IF ((RandomHash2(i) == 0) .or. (RandomHash2(i) > nBasis * 1000)) THEN
CALL Stop_All(t_r, "Random Hash 2 incorrectly calculated")
end if
do j = i + step, nBasis, step
IF (RandomOrbIndex(i) == RandomOrbIndex(j)) THEN
CALL Stop_All(t_r, "Random Hash incorrectly calculated")
end if
IF (RandomHash2(i) == RandomHash2(j)) THEN
CALL Stop_All(t_r, "Random Hash 2 incorrectly calculated")
end if
end do
end do
end if
!Now broadcast to all processors
CALL MPIBCast(RandomOrbIndex, nBasis)
call MPIBCast(RandomHash2, nBasis)
t_initialized_roi = .true.
call init_load_balance()
IF (tHPHF) THEN
!IF(tLatticeGens) CALL Stop_All("SetupParameters","Cannot use HPHF with model systems currently.")
IF (tROHF .or. (LMS /= 0)) CALL Stop_All("SetupParameters", "Cannot use HPHF with high-spin systems.")
tHPHFInts = .true.
end if
if (t_lattice_model) then
if (t_tJ_model) then
if (tGUGA) then
call init_get_helement_tj_guga()
else
call init_get_helement_tj()
end if
else if (t_heisenberg_model) then
if (tGUGA) then
call init_get_helement_heisenberg_guga
else
call init_get_helement_heisenberg()
end if
else if (t_new_real_space_hubbard) then
call init_get_helement_hubbard()
else if (t_k_space_hubbard) then
call init_get_helement_k_space_hub()
end if
end if
!Calculate Hii (unless suppressed)
if (tZeroRef) then
TempHii = 0.0_dp
else IF (tHPHF) THEN
TempHii = hphf_diag_helement(HFDet, iLutHF)
ELSE
TempHii = get_helement(HFDet, HFDet, 0)
end if
Hii = REAL(TempHii, dp)
write(stdout, "(A,F20.10)") "Reference Energy set to: ", Hii
if (tUEG) then
!We require calculation of the sum of fock eigenvalues,
!without knowing them - calculate from the full 1e matrix elements
!of full hamiltonian removing two electron terms.
TempHii = GetH0Element4(HFDet, HFDet)
else
!Know fock eigenvalues, so just use these.
TempHii = GetH0Element3(hfdet)
end if
Fii = REAL(TempHii, dp)
!Find the highest energy determinant...
IF (LMS == 0) THEN
do i = 1, NEl
HighEDet(i) = Brr(nBasis - (i - 1))
end do
call sort(HighEDet)
IF (tHPHF) THEN
call EncodeBitDet(HighEDet, iLutTemp)
TempHii = hphf_diag_helement(HighEDet, iLutTemp)
ELSE
TempHii = get_helement(HighEDet, HighEDet, 0)
end if
if (abs(TempHii - Hii) > EPS) then
UpperTau = 1.0_dp / REAL(TempHii - Hii, dp)
else
UpperTau = 0.0_dp
end if
write(stdout, "(A,G25.15)") "Highest energy determinant is (approximately): ", REAL(TempHii, dp)
write(stdout, "(a,g25.15)") "Corresponding to a correlation energy of: ", real(temphii - hii, dp)
write(stdout, "(A,F25.15)") "This means tau should be no more than about ", UpperTau
write(stdout, *) "Highest energy determinant is: ", HighEDet(:)
associate(deterministic_max_tau => UpperTau * 1.1_dp)
if (deterministic_max_tau < max_tau) then
write(stdout, "(A)") "The deterministic tau is smaller than max_tau."
write(stdout, "(A, F25.15)") "We will limit max_tau to:", deterministic_max_tau
max_tau = deterministic_max_tau
if (tau > max_tau) then
call assign_value_to_tau(max_tau, 'Initial tau was higher than deterministic tau limit.')
end if
end if
end associate
else
UpperTau = 0.0_dp
end if
if (tau_start_val == possible_tau_start%deterministic) then
call assign_value_to_tau(UpperTau, 'Deterministically approximated value 1 / (E_max - E_0)')
end if
! Initialise DiagSft according to the input parameters. If we have
! multiple projected-energy references, then the shift on each of the
! runs should be adjusted so that it is still relative to the first
! replica, but is offset by the replica's reference's diagonal energy.
DiagSft = InputDiagSft
proje_ref_energy_offsets = 0.0_dp
if (tOrthogonaliseReplicas) then
do run = 1, inum_runs
if (tHPHF) then
TempHii = hphf_diag_helement(ProjEDet(:, run), ilutRef(:, run))
else
TempHii = get_helement(ProjEDet(:, run), ProjEDet(:, run), 0)
end if
proje_ref_energy_offsets(run) = real(TempHii, dp) - Hii
if (tMultiRefShift) DiagSft(run) = proje_ref_energy_offsets(run)
end do
end if
IF (tHub) THEN
IF (tReal) THEN
!We also know that in real-space hubbard calculations, there are only single excitations.
exFlag = 1
ELSE
!We are doing a momentum space hubbard calculation - set exFlag to 2 since only doubles are connected for momentum conservation.
exFlag = 2
end if
end if
IF (LMS /= 0) THEN
IF (tNoBrillouin .or. (tHub .and. tReal) .or. tRotatedOrbs) THEN
write(stdout, *) "No brillouin theorem assumed. Single excitations also used to calculate projected energy."
else if (tUHF) THEN
write(stdout, *) "High spin calculation - but single excitations will *NOT* be used to calculate energy as "&
& //"this is an unrestricted calculation."
ELSE
CALL Stop_All("SetupParameters", "High-spin, restricted calculation detected, but single excitations are "&
& //"not being used to calculate the energy. &
& Either use the UHF keyword, or turn off brillouins theorem using NOBRILLOUINS, ROHF or ROTATEDORBS.")
end if
! tRotatedOrbs=.true.
! else if(LMS.ne.0) THEN
! CALL Stop_All(t_r,"Ms not equal to zero, but tSpn is false. Error here")
end if
!Initialise variables for calculation on each node
iter = 0 !This is set so that calls to CalcParentFlag in the initialisation are ok with the logging.
iPopsTimers = 1 !Number of timed popsfiles written out
iBlockingIter = 0
IterTime = 0.0
ProjectionE(:) = 0.0_dp
AvSign = 0.0_dp
AvSignHFD = 0.0_dp
SumENum(:) = 0.0_dp
InitsENumCyc(:) = 0.0_dp
SumNoatHF(:) = 0.0_dp
NoatHF(:) = 0.0_dp
InstNoatHF(:) = 0.0_dp
Annihilated(:) = 0.0_dp
Acceptances(:) = 0.0_dp
PreviousCycles = 0
NoBorn = 0.0_dp
SpawnFromSing = 0
max_cyc_spawn = 0.0_dp
! in case tau-search is off
max_death_cpt = 0.0_dp
NoDied = 0
HFCyc = 0.0_dp
HFOut = 0.0_dp
ENumCyc = 0.0_dp
ENumOut = 0.0_dp
ENUmCycAbs = 0.0_dp
VaryShiftCycles = 0
AvDiagSft(:) = 0.0_dp
SumDiagSft(:) = 0.0_dp
SumWalkersCyc(:) = 0.0_dp
SumWalkersOut(:) = 0.0_dp
! SumDiagSftAbort=0.0_dp
! AvDiagSftAbort=0.0_dp
NoAborted(:) = 0.0_dp
NoRemoved(:) = 0.0_dp
NoAddedInitiators = 0
NoInitDets = 0
NoNonInitDets = 0
NoInitWalk(:) = 0.0_dp
NoNonInitWalk(:) = 0.0_dp
NoExtraInitDoubs = 0
InitRemoved = 0
TotImagTime = 0.0_dp
DiagSftRe = 0.0_dp
DiagSftIm = 0.0_dp
sum_proje_denominator = 0
cyc_proje_denominator = 0
!Also reinitialise the global variables - should not necessarily need to do this...
AllSumENum(:) = 0.0_dp
AllInitsENumCyc(:) = 0.0_dp
AllNoatHF(:) = 0.0_dp
AllNoatDoubs(:) = 0.0_dp
NoAtDoubs = 0.0_dp
if (tLogEXLEVELStats) AllEXLEVEL_WNorm(:, :, :) = 0.0_dp
AllSumNoatHF(:) = 0.0_dp
AllGrowRate(:) = 0.0_dp
AllGrowRateAbort(:) = 0
! AllMeanExcitLevel=0.0_dp
AllSumWalkersCyc(:) = 0
AllSumWalkersOut(:) = 0.0_dp
AllAvSign = 0.0_dp
AllAvSignHFD = 0.0_dp
AllNoBorn(:) = 0
AllSpawnFromSing(:) = 0
AllNoDied(:) = 0
AllAnnihilated(:) = 0
AllENumCyc(:) = 0.0_dp
AllENumOut(:) = 0.0_dp
AllENumCycAbs = 0.0_dp
AllHFCyc(:) = 0.0_dp
all_cyc_proje_denominator = 1.0_dp
AllHFOut(:) = 0.0_dp
! AllDetsNorm=0.0_dp
AllNoAborted = 0
AllNoRemoved = 0
AllNoAddedInitiators = 0
AllNoInitDets = 0
AllNoNonInitDets = 0
AllNoInitWalk = 0.0_dp
AllNoNonInitWalk = 0.0_dp
AllNoExtraInitDoubs = 0
AllInitRemoved = 0
all_n_core_non_init = 0
n_core_non_init = 0
proje_iter = 0
inits_proje_iter = 0.0_dp
AccRat = 0
HFShift = 0
bloom_count = 0
InstShift = 0
AbsProjE = 0
norm_semistoch = 0
norm_psi = 0
bloom_sizes = 0
proje_iter_tot = 0.0_dp
! initialize as one (kind of makes sense for a norm)
norm_psi_squared = 1.0_dp
all_norm_psi_squared = 1.0_dp
norm_semistoch_squared = 1.0_dp
tSoftExitFound = .false.
tReferenceChanged = .false.
! Set the flag to indicate that no shift adjustment has been made
tfirst_cycle = .true.
! Initialise the fciqmc counters
iter_data_fciqmc%update_growth = 0.0_dp
iter_data_fciqmc%update_iters = 0
nExChecks = 0
nExCheckFails = 0
! 0-initialize truncated weight
truncatedWeight = 0.0_dp
AllTruncatedWeight = 0.0_dp
! RDMs are taken as they are until we have some data on the f-function
! of the adaptive shift
rdmCorrectionFactor = 0.0_dp
InstRDMCorrectionFactor = 0.0_dp
ThisRDMIter = 0.0_dp
! initialize excitation number trackers
nInvalidExcits = 0
nValidExcits = 0
allNInvalidExcits = 0
allNValidExcits = 0
if (tEScaleWalkers) then
if (abs(RealSpawnCutoff - sFBeta) > eps) then
write(stdout, *) "Warning: Overriding RealSpawnCutoff with scale function parameter"
RealSpawnCutoff = sFBeta
end if
end if
tNoSinglesPossible = t_k_space_hubbard .or. tUEG .or. tNoSinglesPossible
if (.not. allocated(allInitsPerExLvl)) allocate(allInitsPerExLvl(maxInitExLvlWrite))
if (.not. allocated(initsPerExLvl)) allocate(initsPerExLvl(maxInitExLvlWrite))
initsPerExlvl = 0
allInitsPerExLvl = 0
IF (tHistSpawn .or. (tCalcFCIMCPsi .and. tFCIMC)) THEN
allocate(HistMinInd(NEl))
allocate(HistMinInd2(NEl))
maxdet = 0
do i = 1, nel
maxdet = maxdet + 2**(nbasis - i)
end do
IF (.not. allocated(FCIDets)) THEN
CALL Stop_All(t_r, "A Full Diagonalization is required before histogramming can occur.")
end if
write(stdout, *) "Histogramming spawning wavevector, with Dets=", Det
allocate(Histogram(1:lenof_sign, 1:det), stat=ierr)
IF (ierr /= 0) THEN
CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays ")
end if
Histogram(:, :) = 0.0_dp
allocate(AllHistogram(1:lenof_sign, 1:det), stat=ierr)
allocate(BeforeNormHist(1:det), stat=ierr)
IF (ierr /= 0) CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays")
IF (tHistSpawn) THEN
allocate(InstHist(1:lenof_sign, 1:det), stat=ierr)
IF (ierr /= 0) THEN
CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays")
end if
InstHist(:, :) = 0.0_dp
allocate(AvAnnihil(1:lenof_sign, 1:det), stat=ierr)
IF (ierr /= 0) THEN
CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays")
end if
AvAnnihil(:, :) = 0.0_dp
allocate(InstAnnihil(1:lenof_sign, 1:det), stat=ierr)
IF (ierr /= 0) THEN
CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays")
end if
InstAnnihil(:, :) = 0.0_dp
end if
IF (iProcIndex == 0) THEN
IF (tHistSpawn) THEN
Tot_Unique_Dets_Unit = get_free_unit()
open(Tot_Unique_Dets_Unit, FILE='TOTUNIQUEDETS', STATUS='UNKNOWN')
if (tCalcVariationalEnergy .and. tDiagAllSpaceEver) then
write(Tot_Unique_Dets_Unit, "(A)") "# Iter UniqueDetsEver AvVarEnergy InstVarEnergy GroundE_Ever"
else if (tCalcVariationalEnergy) then
write(Tot_Unique_Dets_Unit, "(A)") "# Iter UniqueDetsEver AvVarEnergy InstVarEnergy"
else if (tDiagAllSpaceEver) then
write(Tot_Unique_Dets_Unit, "(A)") "# Iter UniqueDetsEver GroundE_Ever"
else
write(Tot_Unique_Dets_Unit, "(A)") "# Iter UniqueDetsEver"
end if
allocate(AllInstHist(1:lenof_sign, 1:det), stat=ierr)
allocate(AllInstAnnihil(1:lenof_sign, 1:det), stat=ierr)
allocate(AllAvAnnihil(1:lenof_sign, 1:det), stat=ierr)
end if
IF (ierr /= 0) THEN
CALL Stop_All("SetupParameters", "Error assigning memory for histogramming arrays")
end if
end if
else if (tHistEnergies) THEN
write(stdout, *) "Histogramming the energies of the particles, with iNoBins=", iNoBins, " and BinRange=", BinRange
write(stdout, *) "Histogramming spawning events from ", -OffDiagMax, " with BinRange = ", OffDiagBinRange
iOffDiagNoBins = INT((2.0_dp * OffDiagMax) / OffDiagBinRange) + 1
write(stdout, *) "This gives ", iOffDiagNoBins, " bins to histogram the off-diagonal matrix elements."
allocate(HistogramEnergy(1:iNoBins))
allocate(AttemptHist(1:iNoBins))
allocate(SpawnHist(1:iNoBins))
allocate(SinglesHist(1:iOffDiagNoBins))
allocate(SinglesAttemptHist(1:iOffDiagNoBins))
allocate(SinglesHistOccOcc(1:iOffDiagNoBins))
allocate(SinglesHistOccVirt(1:iOffDiagNoBins))
allocate(SinglesHistVirtOcc(1:iOffDiagNoBins))
allocate(SinglesHistVirtVirt(1:iOffDiagNoBins))
allocate(DoublesHist(1:iOffDiagNoBins))
allocate(DoublesAttemptHist(1:iOffDiagNoBins))
HistogramEnergy(:) = 0.0_dp
AttemptHist(:) = 0.0_dp
SpawnHist(:) = 0.0_dp
SinglesHist(:) = 0.0_dp
SinglesAttemptHist(:) = 0.0_dp
SinglesHistOccOcc(:) = 0.0_dp
SinglesHistOccVirt(:) = 0.0_dp
SinglesHistVirtOcc(:) = 0.0_dp
SinglesHistVirtVirt(:) = 0.0_dp
DoublesHist(:) = 0.0_dp
DoublesAttemptHist(:) = 0.0_dp
IF (iProcIndex == Root) THEN
allocate(AllHistogramEnergy(1:iNoBins))
allocate(AllAttemptHist(1:iNoBins))
allocate(AllSpawnHist(1:iNoBins))
allocate(AllSinglesHist(1:iOffDiagNoBins))
allocate(AllDoublesHist(1:iOffDiagNoBins))
allocate(AllSinglesAttemptHist(1:iOffDiagNoBins))
allocate(AllDoublesAttemptHist(1:iOffDiagNoBins))
allocate(AllSinglesHistOccOcc(1:iOffDiagNoBins))
allocate(AllSinglesHistOccVirt(1:iOffDiagNoBins))
allocate(AllSinglesHistVirtOcc(1:iOffDiagNoBins))
allocate(AllSinglesHistVirtVirt(1:iOffDiagNoBins))
end if
end if
if ((iProcIndex == Root) .and. tDiagWalkerSubspace) then
write(stdout, '(A,I9,A)') "Diagonalising walker subspace every ", iDiagSubspaceIter, " iterations"
unitWalkerDiag = get_free_unit()
open(unitWalkerDiag, file='WalkerSubspaceDiag', status='unknown')
if (tTruncInitiator) then
write(unitWalkerDiag, '(A)') "# Iter NoInitDets NoOccDets InitiatorSubspaceEnergy &
& FullSubspaceEnergy ProjInitEnergy ProjFullEnergy"
else
write(unitWalkerDiag, '(A)') "# Iter NoOccDets InitiatorSubspaceEnergy &
& FullSubspaceEnergy ProjFullEnergy"
end if
end if
if (tHistExcitToFrom) &
call init_hist_excit_tofrom()
IF (tUseBrillouin) THEN
write(stdout, "(A)") "Brillouin theorem in use for calculation of projected energy."
end if
if (.not. (t_k_space_hubbard .and. t_trans_corr_2body)) then
! for too big lattices my implementation breaks, due to
! memory limitations.. but i think we do not actually need it.
CALL CalcApproxpDoubles()
end if
IF (tau_start_val == possible_tau_start%tau_factor) THEN
write(stdout, *) "TauFactor detected. Resetting Tau based on connectivity of: ", HFConn
call assign_value_to_tau(TauFactor / REAL(HFConn, dp), 'Initialization from Tau-Factor.')
end if
if (t_store_ci_coeff) then
call init_ciCoeff()
end if
! [W.D.] I guess I want to initialize that before the tau-search,
! or otherwise some pgens get calculated incorrectly
if (t_back_spawn .or. t_back_spawn_flex) then
call init_back_spawn()
end if
if (t_guga_back_spawn) then
call setup_virtual_mask()
end if
! also i should warn the user if this is a restarted run with a
! set delay in the back-spawning method:
! is there actually a use-case where someone really wants to delay
! a back-spawn in a restarted run?
if (tReadPops .and. back_spawn_delay /= 0) then
call Warning_neci(t_r, &
"Do you really want a delayed back-spawn in a restarted run?")
end if
! can i initialize the k-space hubbard here already?
! because we need information for the tau-search already..
if (t_k_space_hubbard) then
call init_k_space_hubbard()
end if
if (t_new_real_space_hubbard) then
call init_real_space_hubbard
end if
if (t_spin_measurements) then
call init_spin_measurements()
end if
if (t_measure_local_spin) then
call init_local_spin_measure()
end if
IF (abs(StepsSftImag) > 1.0e-12_dp) THEN
write(stdout, *) "StepsShiftImag detected. Resetting StepsShift."
StepsSft = NINT(StepsSftImag / Tau)
IF (StepsSft == 0) StepsSft = 1
write(stdout, *) "StepsShift set to: ", StepsSft
end if
! StepsPrint < 1 while not coupling update and output cycle means no std output
if (.not. tCoupleCycleOutput .and. StepsPrint < 1) then
tMCOutput = .false.
! But there shall be output in the FCIMCStats file
! There is no specified output cycle, so we default to the shift cycle -> couple them
tCoupleCycleOutput = .true.
end if
! Coupling output and shift update means these two are the same
if (tCoupleCycleOutput) StepsPrint = StepsSft
if (StepsPrint == StepsSft) tCoupleCycleOutput = .true.
IF (TPopsFile) THEN
IF (mod(iWritePopsEvery, StepsSft) /= 0) then
CALL Warning_neci(t_r, "POPSFILE writeout should be a multiple of the update cycle length.")
end if
end if
if (TReadPops) then
if (tStartSinglePart .and. .not. tReadPopsRestart) then
if (iProcIndex == root) call warning_neci(t_r, &
"ReadPOPS cannot work with StartSinglePart: ignoring StartSinglePart")
tStartSinglePart = .false.
end if
end if
IF (.not. TReadPops) THEN
write(stdout, "(A,F20.10)") "Initial Diagonal Shift: ", DiagSft(1)
end if
if (tShiftonHFPop) then
write(stdout, *) "Shift will be varied in order to keep the population on the reference determinant fixed"
end if
write(stdout, *) "Connectivity of HF determinant is: ", HFConn
IF (TStartSinglePart) THEN
TSinglePartPhase(:) = .true.
ELSE
TSinglePartPhase(:) = .false.
end if
IF (ICILevel /= 0) THEN
!We are truncating the excitations at a certain value
TTruncSpace = .true.
write(stdout, '(A,I4)') "Truncating the S.D. space at determinants will an excitation level w.r.t. HF of: ", ICILevel
end if
IF (tTruncCAS .or. tStartCAS) THEN
! We are truncating the FCI space by only allowing excitations
! in a predetermined CAS space.
! The following line has already been written out if we are doing
! a CAS calculation.
! The SpinInvBRR array is required for the tTruncCAS option. Its
! properties are explained more fully in the subroutine.
CALL CreateSpinInvBRR()
! CASmax is the max spin orbital number (when ordered
! energetically) within the chosen active space.
! Spin orbitals with energies larger than this maximum value must
! be unoccupied for the determinant to be in the active space.
CASmax = NEl + VirtCASorbs
! CASmin is the max spin orbital number below the active space.
! As well as the above criteria, spin orbitals with energies
! equal to, or below that of the CASmin orbital must be completely
! occupied for the determinant to be in the active space.
CASmin = NEl - OccCASorbs
IF (OccCASOrbs > NEl) then
CALL Stop_All("SetupParameters", "Occupied orbitals in CAS space specified is greater than number of electrons")
end if
IF (VirtCASOrbs > (nBasis - NEl)) then
CALL Stop_All("SetupParameters", "Virtuals in CAS space specified greater than number of unoccupied orbitals")
end if
!Create the bit masks for the bit calculation of these properties.
allocate(CASMask(0:NIfD))
allocate(CoreMask(0:NIfD))
CASMask(:) = 0
CoreMask(:) = 0
do i = 1, nBasis
IF (SpinInvBRR(i) > CASmax) THEN
!Orbital is in external space
CASMask((SpinInvBRR(i) - 1) / bits_n_int) = ibset(CASMask((i - 1) / bits_n_int), MOD((i - 1), bits_n_int))
else if (SpinInvBRR(i) <= CASmin) THEN
!Orbital is in core space
CoreMask((SpinInvBRR(i) - 1) / bits_n_int) = ibset(CoreMask((i - 1) / bits_n_int), MOD((i - 1), bits_n_int))
CASMask((SpinInvBRR(i) - 1) / bits_n_int) = ibset(CASMask((i - 1) / bits_n_int), MOD((i - 1), bits_n_int))
end if
end do
end if
IF (tPartFreezeCore) THEN
write(stdout, '(A,I4,A,I5)') 'Partially freezing the lowest ', NPartFrozen, ' spin orbitals so that no more than ', &
NHolesFrozen, ' holes exist within this core.'
CALL CreateSpinInvBRR()
end if
IF (tPartFreezeVirt) THEN
write(stdout, '(A,I4,A,I5)') 'Partially freezing the highest ', NVirtPartFrozen, &
' virtual spin orbitals so that no more than ', NElVirtFrozen, ' electrons occupy these orbitals.'
CALL CreateSpinInvBRR()
end if
if (tTruncNOpen) then
write (stdout, '("Truncating determinant space at a maximum of ",i3," &
&unpaired electrons.")') trunc_nopen_max
end if
! for the (uniform) 3-body excitgen, the generation probabilities are uniquely given
! by the number of alpha and beta electrons and the number of orbitals
! and can hence be precomputed
if (t_mol_3_body .or. t_ueg_3_body) call setup_mol_tc_excitgen()
if (allocated(SD_spin_purification)) then
if (allocated(spin_pure_J)) then
if (SD_spin_purification == possible_purification_methods%TRUNCATED_LADDER) then
write(stdout, *)
write(stdout, '(A)') 'Spin purification of Slater Determinants &
&with off-diagonal $ J * S_{+} S_{-} $ correction.'
write(stdout, '(A, 1x, E10.5)') 'J =', spin_pure_J
write(stdout, *)
else if (SD_spin_purification == possible_purification_methods%ONLY_LADDER) then
write(stdout, *)
write(stdout, '(A)') 'Spin purification of Slater Determinants &
&with $ J * S_{+} S_{-} $ correction.'
write(stdout, '(A, 1x, E10.5)') 'J =', spin_pure_J
write(stdout, *)
else
write(stdout, *)
write(stdout, '(A)') 'Spin purification of Slater Determinants &
&with full $ J *S^2 $ spin penalty.'
write(stdout, '(A, 1x, E10.5)') 'J =', spin_pure_J
write(stdout, *)
end if
else
call stop_all(this_routine, "spin_pure_J not allocated")
end if
end if
call init_exc_gen_class()
END SUBROUTINE SetupParameters