fcimc_initialisation.F90 Source File


Contents


Source Code

#include "macros.h"

module fcimc_initialisation

    use SystemData, only: tUseBrillouin, iRanLuxLev, tSpn, tHPHFInts, tHPHF, &
                          tRotateOrbs, tROHF, tFindCINatOrbs, nOccBeta, tHUB, &
                          nOccAlpha, tUHF, tBrillouinsDefault, ECore, tUEG, &
                          tNoSingExcits, tOddS_HPHF, tSpn, tNoBrillouin, G1, &
                          tAssumeSizeExcitgen, tMolproMimic, tMolpro, tFixLz, &
                          tRef_Not_HF, LzTot, LMS, tKPntSym, tReal, nBasisMax, &
                          tRotatedOrbs, MolproID, nBasis, arr, brr, nel, &
                          tPickVirtUniform, tGen_4ind_reverse, &
                          tGenHelWeighted, tGen_4ind_weighted, tLatticeGens, &
                          tGUGA, tUEGNewGenerator, &
                          tGen_4ind_2, tReltvy, t_new_real_space_hubbard, &
                          t_lattice_model, t_tJ_model, t_heisenberg_model, &
                          t_k_space_hubbard, t_3_body_excits, breathingCont, &
                          momIndexTable, t_trans_corr_2body, t_non_hermitian_2_body, &
                          tgen_guga_crude, t_impurity_excitgen, &
                          t_uniform_excits, t_mol_3_body,t_ueg_transcorr,t_ueg_3_body,tLatticeGens, &
                          irrepOrbOffset, nIrreps, t_trans_corr_hop, &
                          tTrcorrExgen, nClosedOrbs, irrepOrbOffset, nIrreps, &
                          nOccOrbs, tNoSinglesPossible, t_pcpp_excitgen, &
                          t_fci_pchb_excitgen, tGAS, t_guga_pchb, t_spin_dependent_transcorr, &
                          basisfn, t_mixed_excits, t_hole_focus_excits


    use tc_three_body_data, only: ptriples
    use SymExcitDataMod, only: tBuildOccVirtList, tBuildSpinSepLists
    use core_space_util, only: cs_replicas
    use dSFMT_interface, only: dSFMT_init

    use CalcData, only: G_VMC_Seed, MemoryFacPart, StepsSftImag, &
                        tCheckHighestPop, tSpatialOnlyHash, tStartCAS, &
                        InitialPart, tStartMP1, tReadPops, &
                        InitialPartVec, iReadWalkersRoot, SinglesBias, NMCYC, &
                        tTruncCAS, tTruncInitiator, DiagSft, tFCIMC, &
                        tTrialWavefunction, tSemiStochastic, OccCASOrbs, &
                        VirtCASOrbs, StepsSft, tStartSinglePart, InitWalkers, &
                        tShiftOnHFPop, tReadPopsRestart, tTruncNOpen, &
                        trunc_nopen_max, MemoryFacInit, MaxNoatHF, HFPopThresh, &
                        tAddToInitiator, InitiatorWalkNo, tRestartHighPop, &
                        tAllRealCoeff, tRealCoeffByExcitLevel, tTruncInitiator, &
                        RealCoeffExcitThresh, aliasStem, tPopsAlias, &
                        tDynamicCoreSpace, TargetGrowRate, &
                        TargetGrowRateWalk, InputTargetGrowRate, semistoch_shift_iter, &
                        InputTargetGrowRateWalk, tOrthogonaliseReplicas, &
                        use_spawn_hash_table, tReplicaSingleDetStart, RealSpawnCutoff, &
                        ss_space_in, trial_space_in, init_trial_in, trial_shift_iter, &
                        tContTimeFCIMC, tContTimeFull, tMultipleInitialRefs, &
                        initial_refs, trial_init_reorder, tStartTrialLater, tTrialInit, &
                        ntrial_ex_calc, tPairedReplicas, tMultiRefShift, tPreCond, &
                        tMultipleInitialStates, initial_states, &
                        t_back_spawn, &
                        t_trunc_nopen_diff, t_guga_back_spawn, &
                        t_back_spawn_option, t_back_spawn_flex_option, &
                        t_back_spawn_flex, back_spawn_delay, ScaleWalkers, tfixedN0, &
                        tReplicaEstimates, tDeathBeforeComms, pSinglesIn, pDoublesIn, pTriplesIn, pParallelIn, &
                        tSetInitFlagsBeforeDeath, tSetInitialRunRef, tEN2Init, &
                        tAutoAdaptiveShift, &
                        tInitializeCSF, S2Init, tWalkContgrow, tSkipRef, &
                        AAS_Cut, tTruncInitiator, &
                        tInitiatorSpace, i_space_in, tLinearAdaptiveShift, &
                        tAS_TrialOffset, ShiftOffset, &
                        tSpinProject

    use tau_main, only: tau_search_method, &
        tau_start_val, possible_tau_start, &
        max_death_cpt, min_tau, max_tau, tau, taufactor, &
        assign_value_to_tau, init_tau

    use adi_data, only: tReferenceChanged, tAdiActive, nExChecks, nExCheckFails, &
                        nRefUpdateInterval, SIUpdateInterval

    use Determinants, only: GetH0Element3, GetH0Element4, tDefineDet, &
                            get_helement, get_helement_det_only

    use hphf_integrals, only: hphf_diag_helement, hphf_spawn_sign, &
                              hphf_off_diag_helement_spawn

    use SymData, only: SymLabelList, SymLabelCounts, TwoCycleSymGens, &
                       SymClassSize, nSymLabels, sym_psi

    use dSFMT_interface, only: genrand_real2_dSFMT

    use DeterminantData, only: write_det, write_det_len, FDet

    use LoggingData, only: tTruncRODump, tCalcVariationalEnergy, tReadRDMs, &
                           tDiagAllSpaceEver, tFCIMCStats2, tCalcFCIMCPsi, &
                           tLogComplexPops, tHistExcitToFrom, tPopsFile, &
                           iWritePopsEvery, tRDMOnFly, tCoupleCycleOutput, StepsPrint, &
                           tDiagWalkerSubspace, tPrintOrbOcc, OrbOccs, &
                           tHistInitPops, OrbOccsTag, tHistEnergies, tMCOutput, &
                           HistInitPops, AllHistInitPops, OffDiagMax, &
                           OffDiagBinRange, iDiagSubspaceIter, t_store_ci_coeff, &
                           AllHistInitPopsTag, HistInitPopsTag, tHDF5PopsRead, &
                           tTransitionRDMs, tLogEXLEVELStats, t_no_append_stats, &
                           t_spin_measurements,  t_measure_local_spin, &
                           maxInitExLvlWrite, initsPerExLvl, AllInitsPerExLvl

    use DetCalcData, only: NMRKS, tagNMRKS, FCIDets, NKRY, NBLK, B2L, nCycle, &
                           ICILevel, det
    use IntegralsData, only: tPartFreezeCore, nHolesFrozen, tPartFreezeVirt, &
                             nVirtPartFrozen, nPartFrozen, nelVirtFrozen

    use bit_rep_data, only: NIfTot, NIfD, IlutBits, flag_initiator, &
                            flag_deterministic, extract_sign, nIfGUGA, &
                            test_flag_multi

    use bit_reps, only: encode_det, clear_all_flags, set_flag, encode_sign, &
                        decode_bit_det, nullify_ilut, encode_part_sign, &
                        extract_run_sign, &
                        get_initiator_flag, writebitdet, &
                        get_initiator_flag_by_run
    use hist_data, only: tHistSpawn, HistMinInd, HistMinInd2, Histogram, &
                         BeforeNormHist, InstHist, iNoBins, AllInstHist, &
                         HistogramEnergy, AllHistogramEnergy, AllHistogram, &
                         BinRange

    use hist, only: init_hist_excit_tofrom, clean_hist_excit_tofrom
    use PopsfileMod, only: FindPopsfileVersion, initfcimc_pops, &
                           ReadFromPopsfilePar, ReadPopsHeadv3, &
                           ReadPopsHeadv4, open_pops_head, checkpopsparams
    use HPHFRandExcitMod, only: gen_hphf_excit, FindDetSpinSym, exc_generator_for_HPHF
    use GenRandSymExcitNUMod, only: gen_rand_excit, init_excit_gen_store, &
                                    clean_excit_gen_store
    use replica_estimates, only: open_replica_est_file
    use procedure_pointers, only: generate_excitation, attempt_create, &
                                  get_spawn_helement, encode_child, &
                                  attempt_die, extract_bit_rep_avsign, &
                                  fill_rdm_diag_currdet, &
                                  get_conn_helement, scaleFunction, &
                                  generate_two_body_excitation, shiftFactorFunction, gen_all_excits
    use symrandexcit3, only: gen_rand_excit3
    use symrandexcit_Ex_Mag, only: gen_rand_excit_Ex_Mag
    use excit_gens_int_weighted, only: gen_excit_hel_weighted, &
                                       gen_excit_4ind_weighted, &
                                       gen_excit_4ind_reverse
    use hash, only: FindWalkerHash, add_hash_table_entry, init_hash_table, &
                    hash_table_lookup, RandomHash2
    use load_balance_calcnodes, only: DetermineDetNode, RandomOrbIndex
    use load_balance, only: tLoadBalanceBlocks, addNormContribution, &
                            AddNewHashDet, clean_load_balance, &
                            init_load_balance
    use matel_getter, only: get_diagonal_matel, get_off_diagonal_matel
    use SymExcit3, only: CountExcitations3, GenExcitations3
    use SymExcit4, only: CountExcitations4, GenExcitations4
    use HPHFRandExcitMod, only: ReturnAlphaOpenDet
    use FciMCLoggingMOD, only: InitHistInitPops
    use SymExcitDataMod, only: SymLabelList2, OrbClassCount, SymLabelCounts2, &
        tBuildSpinSepLists
    use rdm_general, only: init_rdms, dealloc_global_rdm_data, &
                           extract_bit_rep_avsign_no_rdm
    use rdm_filling, only: fill_rdm_diag_currdet_norm
    use DetBitOps, only: FindBitExcitLevel, CountBits, TestClosedShellDet, &
                         FindExcitBitDet, IsAllowedHPHF, DetBitEq, &
                         EncodeBitDet, DetBitLT
    use fcimc_pointed_fns, only: att_create_trunc_spawn_enc, &
                                 attempt_create_normal, &
                                 attempt_create_trunc_spawn, &
                                 new_child_stats_normal, &
                                 null_encode_child, attempt_die_normal, attempt_die_precond, &
                                 powerScaleFunction, expScaleFunction, negScaleFunction, &
                                 expCOScaleFunction, constShiftFactorFunction, &
                                 linearShiftFactorFunction, autoShiftFactorFunction

    use initial_trial_states, only: calc_trial_states_lanczos, &
                                    set_trial_populations, set_trial_states, calc_trial_states_direct
    use global_det_data, only: global_determinant_data, set_det_diagH, &
                               set_det_offdiagH, clean_global_det_data, &
                               init_global_det_data, set_spawn_rate, &
                               store_decoding, set_supergroup_idx
    use semi_stoch_gen, only: init_semi_stochastic, end_semistoch, &
                              enumerate_sing_doub_kpnt
    use semi_stoch_procs, only: return_mp1_amp_and_mp2_energy
    use initiator_space_procs, only: init_initiator_space
    use kp_fciqmc_data_mod, only: tExcitedStateKP
    use sym_general_mod, only: ClassCountInd
    use trial_wf_gen, only: init_trial_wf, end_trial_wf
    use ueg_excit_gens, only: gen_ueg_excit
    use gndts_mod, only: gndts
    use excit_gen_5, only: gen_excit_4ind_weighted2
    use tc_three_body_excitgen, only: gen_excit_mol_tc, setup_mol_tc_excitgen
    use pcpp_excitgen, only: gen_rand_excit_pcpp, init_pcpp_excitgen, finalize_pcpp_excitgen

    use fcimc_helper, only: CalcParentFlag, update_run_reference

    use cont_time_rates, only: spawn_rate_full, oversample_factors, &
                               secondary_gen_store, ostag

    use soft_exit, only: tSoftExitFound

    use get_excit, only: make_double

    use excitation_types, only: Excite_0_t
    use sltcnd_mod, only: sltcnd_excit
    use rdm_data, only: nrdms_transition_input, rdmCorrectionFactor, InstRDMCorrectionFactor, &
                        ThisRDMIter
    use rdm_data, only: nrdms_transition_input

    use Parallel_neci, only: MPI_2INTEGER, root, nProcessors, iProcIndex, &
        MPISumAll, MPIBarrier, MPIBCast, MPISum, MPI_MAXLOC, nNodes, &
        Sync_Time, MPIAllReduceDatatype

    use util_mod, only: stop_all, get_free_unit, neci_flush, &
        operator(.isclose.), operator(.div.), warning_neci

    use fortran_strings, only: str


    use sym_mod, only: getsym_wrapper, WRITESYM, GetLz, GetSym

    use constants, only: dp, int32, int64, n_int, stdout, stderr, &
        lenof_sign, inum_runs, EPS, bits_n_int, size_n_int, maxExcit, &
        MPIArg

    use MemoryManager, only: LogMemAlloc, LogMemDealloc

    use FciMCData, only: &
        Walker_Time, Annihil_Time, Sort_Time, &
        Comms_Time, ACF_Time, AnnSpawned_time, AnnMain_time, BinSearch_time, &
        SemiStoch_Comms_Time, SemiStoch_Multiply_Time, Trial_Search_Time, SemiStoch_Init_Time, SemiStoch_Hamil_Time, &
        SemiStoch_Davidson_Time, Trial_Init_Time, InitSpace_Init_Time, kp_generate_time, Stats_Comms_Time, &
        subspace_hamil_time, exact_subspace_h_time, subspace_spin_time, var_e_time, precond_e_time, &
        proj_e_time, rescale_time, death_time, hash_test_time, hii_test_time, &
        init_flag_time, iter_data_fciqmc, ll_node, GAS_PCHB_init_time

    use FciMCData, only: all_norms, tDetermSpawnedTo, core_run, &
        tSinglePartPhase, VaryShiftIter, con_send_buf, &
        popsfile_dets, alloc_popsfile_dets

    use FciMCData, only: &
        nWalkerHashes, nreplicas, MaxWalkersUncorrected, MaxSpawned, &
        HashLengthFrac, fcimc_excit_gen_store, tGenMatHEL, TempSpawnedPartsTag, &
        SpinInvBRRTag, pSingles, pDoubles, &
        pParallel, pSing_spindiff1, pDoub_spindiff1, pDoub_spindiff2, &
        tReplicaReferencesDiffer, old_norm_psi

    use FciMCData, only: &
        AvSignHFD, AvSign, CASMask, CoreMask, AllTruncatedWeight, &
        AllTotPartsOld, AllTotPartsLastOutput, allNValidExcits, &
        AllNoRemoved, AllNoNonInitWalk, AllNoNonInitDets, &
        tFillingExplicRDMonFly, tFillingStochRDMonFly, &
        tot_init_trial_denom, tot_init_trial_numerator, tot_trial_denom, &
        tot_trial_denom_inst, tot_trial_num_inst, tot_trial_numerator, &
        tPrintHighPop, trial_denom, trial_denom_inst, &
        pops_pert, &
        FreeSlot, InstAnnihil, AvAnnihil, AllAvAnnihil, AllInstAnnihil, &
        AttemptHist, SpawnHist, SinglesHist, DoublesHist, DoublesAttemptHist, &
        SinglesAttemptHist, SinglesHistOccOcc, SinglesHistVirtOcc, SinglesHistOccVirt, &
        SinglesHistVirtVirt, AllAttemptHist, AllSpawnHist, AllSinglesAttemptHist, &
         AllSinglesHist, AllDoublesAttemptHist, AllDoublesHist, AllSinglesHistOccOcc, &
        AllSinglesHistVirtOcc, AllSinglesHistOccVirt, AllSinglesHistVirtVirt, &
        TempSpawnedParts, HighestPopDet, refdetflip, ilutrefflip

    use FciMCData, only: &
        ProjectionE, SumENum, InitsENumCyc, SumNoatHF, Annihilated, Acceptances, AvDiagSft, SumDiagSft, &
        SumWalkersCyc, SumWalkersOut, NoAborted, NoRemoved, NoInitWalk, NoNonInitWalk, &
        AllSumENum, AllInitsENumCyc, AllNoatDoubs, AllEXLEVEL_WNorm, AllSumNoatHF, &
        AllGrowRate, AllGrowRateAbort, AllSumWalkersCyc, AllSumWalkersOut, AllNoBorn, &
        AllSpawnFromSing, AllNoDied, AllAnnihilated, AllENumCyc, AllENumOut, AllHFCyc, &
        AllHFOut, replica_overlaps_real, ValidSpawnedList
#ifdef CMPLX_
    use FciMCData, only: replica_overlaps_imag
#endif

    use FciMCData, only: HFDet_True, HFDet, tSpinCoupProjE, SpawnInfoVec2Tag, &
        SpawnInfoVecTag, iRefProc, SpawnInfoVec, SpawnInfoVec2, &
        SpawnVec2, SpawnVec2Tag, SpawnVec, SpawnVecTag, &
        HashIndex, AllNoAbortedOld, sfTag, OldAllHFCyc, OldAllAvWalkersCyc, &
        TotPartsOld, NoatHF, InitialSpawnedSlots, TotParts, &
        proje_ref_energy_offsets, AllTotParts, ProjEDet, VaryShiftCycles, &
        unitWalkerDiag, tZeroRef, truncatedWeight, TTruncSpace, &
        tTimeExit, tRestart, TotImagTime, Tot_Unique_Dets_Unit, &
        tfirst_cycle, tEScaleWalkers, TDebug, t_initialized_roi, &
        sum_proje_denominator, SpawnFromSing, sFBeta, proje_iter_tot, &
        proje_iter, PreviousCycles, nValidExcits, nInvalidExcits, &
        nSingles, nDoubles, norm_semistoch_squared, norm_semistoch, &
        norm_psi_squared, norm_psi, NoNonInitDets, NoInitDets, NoExtraInitDoubs, &
        NoDied, NoBorn, NoAtDoubs, NoAddedInitiators, n_core_non_init, &
        MaxTimeExit, maxdet, max_cyc_spawn, IterTime, iter, &
        iPopsTimers, iOffDiagNoBins, InstShift, InputDiagSft, &
        inits_proje_iter, InitRemoved, initiatorstats_unit, &
        ilutRef, iLutHF_True, iLutHF, iHighestPop, iBlockingIter, &
        Hii, HFSym, HFShift, HFOut, HFDetTag, HFCyc, HFConn, &
        hash_iter, Fii, fcimcstats_unit, fcimcstats_unit2, EXLEVELStats_unit, &
        DiagSftRe, DiagSftIm, ENumCyc, ENumCycAbs, ENumOut, exflag, &
        CASmin, CASmax, bloom_count, bloom_sizes, cyc_proje_denominator, &
        SpinInvBRR, ComplexStats_unit, all_cyc_proje_denominator, &
        all_n_core_non_init, all_norm_psi_squared, AllAvSign, AllAvSignHFD, &
        AllENumCycAbs, AllInitRemoved, allNInvalidExcits, AllNoAborted, &
        AllNoAddedInitiators, AllNoExtraInitDoubs, AllNoInitDets, &
        AllNoInitWalk, trial_num_inst, trial_numerator, tSinglePartPhase, &
        CurrentDets, AbsProjE, AccRat, &
        bloom_max, bloom_warn_string, max_calc_ex_level, con_space_size, &
        CurrentDets, init_trial_denom, init_trial_numerator, MaxInitPopNeg, &
        MaxInitPopPos, MaxWalkersPart, nhashes_spawn, &
        SpawnedParts, SpawnedParts2, SpawnInfo, SpawnInfo2, &
        WalkVecDets, WalkVecDetsTag, tPopsAlreadyRead, &
        OldAllNoatHF, TotWalkers, TotWalkersOld, AllNoatHF, AllTotWalkers, &
        AllTotWalkersOld, HFInd, InstNoatHF, OldAllNoatHF, &
        TotWalkers, TotWalkersOld, iEndFreeSlot, iStartFreeSlot, OldAllNoatHF, &
        SpawnInfoWidth, spawn_ht


    use sort_mod, only: sort

    ! use HElem


    use guga_bitRepOps, only: calcB_vector_nI, calcB_vector_ilut, convert_ilut_toNECI, &
                              convert_ilut_toGUGA, getDeltaB, write_det_guga, write_guga_list, &
                              calc_csf_i, CSF_Info_t, csf_ref, fill_csf_i

    use guga_main, only: generate_excitation_guga
    use guga_excitations, only: actHamiltonian
    use guga_matrixElements, only: calcDiagMatEleGUGA_ilut, calcDiagMatEleGUGA_nI

    use real_time_data, only: t_real_time_fciqmc

    use real_time_procs, only: attempt_create_realtime

    use adi_references, only: setup_reference_space, clean_adi

    use double_occ_mod, only: init_spin_measurements

    use back_spawn, only: init_back_spawn, setup_virtual_mask

    use real_space_hubbard, only: init_real_space_hubbard, init_get_helement_hubbard, &
                                  gen_excit_rs_hubbard, &
                                  gen_excit_rs_hubbard_transcorr_hole_focus, &
                                  gen_excit_rs_hubbard_transcorr_uniform, &
                                  gen_excit_rs_hubbard_transcorr, &
                                  gen_excit_rs_hubbard_spin_dependent_transcorr


    use sdt_amplitudes, only: init_ciCoeff

    use back_spawn_excit_gen, only: gen_excit_back_spawn, gen_excit_back_spawn_ueg, &
                                    gen_excit_back_spawn_hubbard, gen_excit_back_spawn_ueg_new

    use gasci_supergroup_index, only: lookup_supergroup_indexer

    use cepa_shifts, only: t_cepa_shift, init_cepa_shifts

    use tj_model, only: init_get_helement_tj, init_get_helement_heisenberg, &
                        init_get_helement_heisenberg_guga, init_get_helement_tj_guga, &
                        gen_excit_tJ_model, gen_excit_heisenberg_model

    use k_space_hubbard, only: init_get_helement_k_space_hub, init_k_space_hubbard, &
                               gen_excit_k_space_hub_transcorr, gen_excit_uniform_k_space_hub_transcorr, &
                               gen_excit_k_space_hub, &
                               gen_excit_uniform_k_space_hub, &
                               gen_excit_mixed_k_space_hub_transcorr

    use OneEInts, only: tmat2d

    use lattice_models_utils, only: gen_all_excits_k_space_hubbard, gen_all_excits_r_space_hubbard

    use impurity_models, only: setupImpurityExcitgen, clearImpurityExcitgen, gen_excit_impurity_model

    use guga_pchb_excitgen, only: init_guga_pchb_excitgen, finalize_pchb_excitgen_guga

    use symexcit3, only: gen_all_excits_default => gen_all_excits

    use local_spin, only: init_local_spin_measure

    use CAS_distribution_init, only: InitFCIMC_CAS

    use SD_spin_purification_mod, only: SD_spin_purification, possible_purification_methods, spin_pure_J

    use exc_gen_classes, only: init_exc_gen_class, finalize_exz_gen_class, class_managed

    use blas_interface_mod, only: dgeev
    implicit none

    external :: LargestBitSet

contains

    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

    ! This initialises the calculation, by allocating memory, setting up the
    ! initial walkers, and reading from a file if needed
    SUBROUTINE InitFCIMCCalcPar()
        INTEGER :: ierr, iunithead
        logical :: formpops, binpops, tStartedFromCoreGround
        integer(int64) :: MemoryAlloc
        INTEGER :: error, PopsVersion
        character(*), parameter :: t_r = 'InitFCIMCPar', this_routine = t_r
        integer :: PopBlockingIter
        real(dp) :: ExpectedMemWalk, read_tau, read_psingles, read_pparallel, read_ptriples
        integer(int64) :: read_walkers_on_nodes(0:nProcessors - 1)
        integer :: read_nnodes
        !Variables from popsfile header...
        logical :: tPop64Bit, tPopHPHF, tPopLz
        integer :: iPopLenof_sign, iPopNel, iPopIter, PopNIfD, PopNIfSgn, PopNIfFlag, PopNIfTot, Popinum_runs
        integer :: PopRandomHash(2056), PopBalanceBlocks
        integer(int64) :: iPopAllTotWalkers
        integer :: i, run
        real(dp) :: PopDiagSft(1:inum_runs)
        real(dp), dimension(lenof_sign) :: PopSumNoatHF
        HElement_t(dp) :: PopAllSumENum(1:inum_runs)
        integer :: perturb_ncreate, perturb_nannihilate
        integer :: nrdms_standard, nrdms_transition
        character(255) :: identifier
        ! set default pops version, this should not have any functional impact,
        ! but prevents using it uninitialized
        !default
        Popinum_runs = 1
        ! default version for popsfiles, this does not have any functional effect,
        ! but prevents it from using uninitialized
        PopsVersion = 4

        if (tPopsAlias) then
            identifier = aliasStem
        else
            identifier = 'POPSFILE'
        end if

        if (tReadPops .and. .not. (tPopsAlreadyRead .or. tHDF5PopsRead)) then
            call open_pops_head(iunithead, formpops, binpops, identifier)
            PopsVersion = FindPopsfileVersion(iunithead)
            if (iProcIndex == root) close(iunithead)
            write(stdout, *) "POPSFILE VERSION ", PopsVersion, " detected."
        end if

        if (tReadPops .and. (PopsVersion < 3) .and. &
            .not. (tPopsAlreadyRead .or. tHDF5PopsRead)) then
!Read in particles from multiple POPSFILES for each processor
            !Ugh - need to set up ValidSpawnedList here too...
            call SetupValidSpawned(int(InitWalkers, int64))
            write(stdout, *) "Reading in initial particle configuration from *OLD* POPSFILES..."
            CALL ReadFromPopsFilePar()
        ELSE
            !Scale walker number
            !This is needed to be done here rather than later,
            !because the arrays should be allocated with appropariate sizes
            if (tReadPops .and. .not. tPopsAlreadyRead) then
                InitWalkers = InitWalkers * ScaleWalkers
            end if

!initialise the particle positions - start at HF with positive sign
!Set the maximum number of walkers allowed
            if (tReadPops .and. .not. (tPopsAlreadyRead .or. tHDF5PopsRead)) then
                !Read header.
                call open_pops_head(iunithead, formpops, binpops, identifier)
                if (PopsVersion == 3) then
                    call ReadPopsHeadv3(iunithead, tPop64Bit, tPopHPHF, tPopLz, iPopLenof_Sign, iPopNel, &
                                        iPopAllTotWalkers, PopDiagSft, PopSumNoatHF, PopAllSumENum, iPopIter, &
                                        PopNIfD, PopNIfSgn, PopNIfFlag, PopNIfTot)

                    ! The following values were not read in...
                    read_tau = 0.0_dp
                    read_nnodes = 0
                    PopBalanceBlocks = -1
                else if (PopsVersion == 4) then
                    ! The only difference between 3 & 4 is just that 4 reads
                    ! in via a namelist, so that we can add more details
                    ! whenever we want.
                    call ReadPopsHeadv4(iunithead, tPop64Bit, tPopHPHF, tPopLz, iPopLenof_Sign, iPopNel, &
                                        iPopAllTotWalkers, PopDiagSft, PopSumNoatHF, PopAllSumENum, iPopIter, &
                                        PopNIfD, PopNIfSgn, Popinum_runs, PopNIfFlag, PopNIfTot, &
                                        read_tau, PopBlockingIter, PopRandomHash, read_psingles, &
                                        read_pparallel, read_ptriples, read_nnodes, read_walkers_on_nodes, &
                                        PopBalanceBlocks)
                else
                    call stop_all(this_routine, "Popsfile version invalid")
                end if

                ! Check the number of electrons created and annihilated by the
                ! perturbation operators, if any are being used.
                if (allocated(pops_pert)) then
                    perturb_ncreate = pops_pert(1)%ncreate
                    perturb_nannihilate = pops_pert(1)%nannihilate
                else
                    perturb_ncreate = 0
                    perturb_nannihilate = 0
                end if

                call CheckPopsParams(tPop64Bit, tPopHPHF, tPopLz, iPopLenof_Sign, iPopNel, &
                                     iPopAllTotWalkers, PopDiagSft, PopSumNoatHF, PopAllSumENum, iPopIter, &
                                     PopNIfD, PopNIfSgn, PopNIfTot, &
                                     MaxWalkersUncorrected, read_tau, PopBlockingIter, read_psingles, read_pparallel, &
                                     read_ptriples, perturb_ncreate, perturb_nannihilate)

                if (iProcIndex == root) close(iunithead)
            else
                MaxWalkersUncorrected = int(InitWalkers)
            end if

            MaxWalkersPart = NINT(MemoryFacPart * MaxWalkersUncorrected)
            ExpectedMemWalk = real( (NIfTot + 1) * size_n_int + 8, dp ) * &
                real(MaxWalkersPart, dp) / 2._dp**20
            if (ExpectedMemWalk < 20.0) then
                !Increase memory allowance for small runs to a min of 20mb
                MaxWalkersPart = ceiling( 20._dp * 2._dp**20 / &
                    real((NIfTot + 1) * size_n_int + 8, dp) )
                block
                    character(*), parameter :: mem_warning = "Low memory requested for walkers, so increasing memory to 20Mb to avoid memory errors"
                    if (iProcIndex == root) write(stderr, "(A)") mem_warning
                end block
            end if
            write(stdout, "(A,I14)") "Memory allocated for a maximum particle number per node of: ", MaxWalkersPart
            !Here is where MaxSpawned is set up - do we want to set up a minimum allocation here too?
            Call SetupValidSpawned(int(InitWalkers, int64))

!Put a barrier here so all processes synchronise
            CALL MPIBarrier(error)
!Allocate memory to hold walkers
            allocate(WalkVecDets(0:NIfTot, MaxWalkersPart), stat=ierr)
            CALL LogMemAlloc('WalkVecDets', MaxWalkersPart * (NIfTot + 1), size_n_int, this_routine, WalkVecDetsTag, ierr)
            WalkVecDets(0:NIfTot, 1:MaxWalkersPart) = 0
            MemoryAlloc = (NIfTot + 1) * MaxWalkersPart * size_n_int    !Memory Allocated in bytes

            if (alloc_popsfile_dets) allocate(popsfile_dets(0:NIfTot, MaxWalkersPart), stat=ierr)

            nrdms_standard = 0
            nrdms_transition = 0

            if (tRDMOnFly) then
                if (tPairedReplicas) then
                    nrdms_standard = lenof_sign.div.2
                else
                    nrdms_standard = lenof_sign
                end if
                if (tTransitionRDMs) then
                    ! nrdms_transition_input will have been read in from the user
                    ! input. But if we have two different replicas for each state
                    ! sampled, then there are two ways to form each transition RDMs.
                    nrdms_transition = nrdms_transition_input * nreplicas
                end if
            end if

            ! Allocate storage for persistent data to be stored alongside
            ! the current determinant list (particularly diagonal matrix
            ! elements, i.e. CurrentH; now global_determinant_data).
            call init_global_det_data(nrdms_standard, nrdms_transition)

            ! If we are doing cont time, then initialise it here
            call init_cont_time()

            ! set the dummies for trial wavefunction connected space
            ! load balancing before trial wf initialization
            if (tTrialWavefunction) then
                allocate(con_send_buf(0, 0))
                con_space_size = 0
            end if

            write(stdout, "(A,I12,A)") "Spawning vectors allowing for a total of ", MaxSpawned, &
                " particles to be spawned in any one iteration per core."
            write(stdout, *) "Memory requirement ", IlutBits%len_bcast * 8.0_dp * ( &
                MaxSpawned / 1048576.0_dp), "MB"

            allocate(SpawnVec(0:IlutBits%len_bcast, MaxSpawned), &
                      stat=ierr, source=0_n_int)
block
    use constants, only: int64
    use util_mod, only: operator(.div.)
    if (ierr /= 0) then
        call stop_all(this_routine, 'Error in allocation of SpawnVec.')
    end if
call LogMemAlloc("SpawnVec", size(SpawnVec, kind=int64), int(sizeof(SpawnVec), kind=int64) .div. size(SpawnVec, kind=int64),&
    & this_routine, SpawnVecTag, ierr)
end block
            allocate(SpawnVec2(0:IlutBits%len_bcast, MaxSpawned), &
                      stat=ierr, source=0_n_int)
block
    use constants, only: int64
    use util_mod, only: operator(.div.)
    if (ierr /= 0) then
        call stop_all(this_routine, 'Error in allocation of SpawnVec2.')
    end if
call LogMemAlloc("SpawnVec2", size(SpawnVec2, kind=int64), int(sizeof(SpawnVec2), kind=int64) .div. size(SpawnVec2, kind=int64),&
    & this_routine, SpawnVec2Tag, ierr)
end block

            if (use_spawn_hash_table) then
                nhashes_spawn = ceiling(0.8 * MaxSpawned)
                allocate(spawn_ht(nhashes_spawn), stat=ierr)
                call init_hash_table(spawn_ht)
            end if

!Point at correct spawning arrays
            SpawnedParts => SpawnVec
            SpawnedParts2 => SpawnVec2

            MemoryAlloc = MemoryAlloc + (NIfTot + 1) * MaxSpawned * 2 * size_n_int

            if (tAutoAdaptiveShift) then
                allocate(SpawnInfoVec(1:SpawnInfoWidth, MaxSpawned), stat=ierr)
block
    use constants, only: int64
    use util_mod, only: operator(.div.)
    if (ierr /= 0) then
        call stop_all(this_routine, 'Error in allocation of SpawnInfoVec.')
    end if
call LogMemAlloc("SpawnInfoVec", size(SpawnInfoVec, kind=int64), int(sizeof(SpawnInfoVec), kind=int64) .div. size(SpawnInfoVec,&
    & kind=int64), this_routine, SpawnInfoVecTag, ierr)
end block
                allocate(SpawnInfoVec2(1:SpawnInfoWidth, MaxSpawned), stat=ierr)
block
    use constants, only: int64
    use util_mod, only: operator(.div.)
    if (ierr /= 0) then
        call stop_all(this_routine, 'Error in allocation of SpawnInfoVec2.')
    end if
call LogMemAlloc("SpawnInfoVec2", size(SpawnInfoVec2, kind=int64), int(sizeof(SpawnInfoVec2), kind=int64) .div.&
    & size(SpawnInfoVec2, kind=int64), this_routine, SpawnInfoVec2Tag, ierr)
end block
                SpawnInfoVec(:, :) = 0
                SpawnInfoVec2(:, :) = 0
                SpawnInfo => SpawnInfoVec
                SpawnInfo2 => SpawnInfoVec2
                MemoryAlloc = MemoryAlloc + (SpawnInfoWidth) * MaxSpawned * 2 * size_n_int
            end if

            write(stdout, "(A)") "Storing walkers in hash-table. Algorithm is now formally linear scaling with walker number"
            write(stdout, "(A,I15)") "Length of hash-table: ", nWalkerHashes
            write(stdout, "(A,F20.5)") "Length of hash-table as a fraction of targetwalkers: ", HashLengthFrac
            ! TODO: Correct the memory usage.
            !MemTemp=2*(8*(nClashMax+1)*nWalkerHashes)+8*MaxWalkersPart
            !write(stdout,"(A,F10.3,A)") "This will use ",real(MemTemp,dp)/1048576.0_dp,&
            !    " Mb of memory per process, although this is likely to increase as it expands"
            allocate(HashIndex(nWalkerHashes), stat=ierr)
            if (ierr /= 0) call stop_all(this_routine, "Error in allocation")
            do i = 1, nWalkerHashes
                HashIndex(i)%ind = 0
            end do
            !Also need to allocate memory for the freeslot array
            allocate(FreeSlot(MaxWalkersPart), stat=ierr)
            if (ierr /= 0) call stop_all(this_routine, "Error in allocation")
            freeslot(:) = 0
            !MemoryAlloc=MemoryAlloc+MemTemp

            ! Allocate pointers to the correct walker arrays
            CurrentDets => WalkVecDets

            ! Get the (0-based) processor index for the HF det.
            do run = 1, inum_runs
                iRefProc(run) = DetermineDetNode(nel, ProjEDet(:, run), 0)
            end do
            write(stdout, "(A,I8)") "Reference processor is: ", iRefProc(1)
            write(stdout, "(A)", advance='no') "Initial reference is: "
            call write_det(stdout, ProjEDet(:, 1), .true.)

            TotParts(:) = 0.0
            TotPartsOld(:) = 0.0
            NoatHF(:) = 0.0_dp
            InstNoatHF(:) = 0.0_dp

            ! Has been moved to guarantee initialization before first load balancing
            ! Initialises RDM stuff for both explicit and stochastic calculations of RDM.

            tFillingStochRDMonFly = .false.
            tFillingExplicRDMonFly = .false.
            !One of these becomes true when we have reached the relevant iteration to begin filling the RDM.

            ! initialize excitation generator
            if (t_guga_pchb) call init_guga_pchb_excitgen()

            ! If we have a popsfile, read the walkers in now.
            if (tReadPops .and. .not. tPopsAlreadyRead) then
                call InitFCIMC_pops(iPopAllTotWalkers, PopNIfSgn, iPopNel, read_nnodes, &
                                    read_walkers_on_nodes, pops_pert, &
                                    PopBalanceBLocks, PopDiagSft)
            else
                if (tStartMP1) then
                    !Initialise walkers according to mp1 amplitude.
                    call InitFCIMC_MP1()

                else if (tStartCAS) then
                    !Initialise walkers according to a CAS diagonalisation.
                    call InitFCIMC_CAS()

                else if (tTrialInit .or. (tOrthogonaliseReplicas .and. &
                                          .not. tReplicaSingleDetStart)) then

                    call InitFCIMC_trial()

                else if (tInitializeCSF) then

                    call InitFCIMC_CSF()

                else !Set up walkers on HF det

                    if (tStartSinglePart) then
                        write(stdout, "(A,I10,A,F9.3,A,I15)") "Initial number of particles set to ", int(InitialPart), &
                            " and shift will be held at ", DiagSft(1), " until particle number gets to ", int(InitWalkers * nNodes)
                    else
                        write(stdout, "(A,I16)") "Initial number of walkers per processor chosen to be: ", nint(InitWalkers)
                    end if

                    call InitFCIMC_HF()

                end if   !tStartmp1
            end if
            if (t_ueg_3_body .and. tTrcorrExgen) tLatticeGens = .false.
            write(stdout, "(A,F14.6,A)") " Initial memory (without excitgens + temp arrays) consists of : ", &
                & REAL(MemoryAlloc, dp) / 1048576.0_dp, " Mb/Processor"
            write(stdout, *) "Only one array of memory to store main particle list allocated..."
            write(stdout, *) "Initial memory allocation sucessful..."
            write(stdout, *) "============================================="
            CALL neci_flush(stdout)

        end if   !End if initial walkers method

        ! There was no last output, use the same value as for the shift update
        AllTotPartsLastOutput = AllTotPartsOld
!Put a barrier here so all processes synchronise
        CALL MPIBarrier(error)

        call init_norm()

        IF (tPrintOrbOcc) THEN
            allocate(OrbOccs(nBasis), stat=ierr)
            CALL LogMemAlloc('OrbOccs', nBasis, 8, this_routine, OrbOccsTag, ierr)
            OrbOccs(:) = 0.0_dp
        end if

        IF (tHistInitPops) THEN
            CALL InitHistInitPops()
        end if
        tPrintHighPop = .false.
        MaxInitPopPos = 0.0
        MaxInitPopNeg = 0.0

        IF (abs(MaxNoatHF) < 1.0e-12_dp) THEN
            MaxNoatHF = InitWalkers * nNodes
            HFPopThresh = int(MaxNoatHF, int64)
        end if

        ! Initialise excitation generation storage
        call init_excit_gen_store(fcimc_excit_gen_store)


        if (t_pcpp_excitgen) call init_pcpp_excitgen()

        if (t_impurity_excitgen) call setupImpurityExcitgen()
        ! [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
        ! 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

        call init_tau()

        IF ((NMCyc /= 0) .and. (tRotateOrbs .and. (.not. tFindCINatOrbs))) then
            CALL Stop_All(this_routine, "Currently not set up to rotate and then go straight into a spawning &
            & calculation.  Ordering of orbitals is incorrect.  This may be fixed if needed.")
        end if

        if (tRDMonFly) then
            call init_rdms(nrdms_standard, nrdms_transition)
            !If the iteration specified to start filling the RDM has already been, want to
            !start filling as soon as possible.
            do run = 1, inum_runs
                if (.not. tSinglePartPhase(run)) VaryShiftIter(run) = 0
            end do
        end if

        ! Perform all semi-stochastic initialisation. This includes generating all the states in the
        ! deterministic space, finding their processors, ordering them, inserting them into
        ! CurrentDets, calculating and storing all Hamiltonian matrix elements and initalising all
        ! arrays required to store and distribute the vectors in the deterministic space later.

        ! in the real-time application, this is done after the initial state is set up
        if (tSemiStochastic .and. .not. t_real_time_fciqmc) then
            if (tDynamicCoreSpace .and. tRDMonFly) then
                tSemiStochastic = .false.
                semistoch_shift_iter = 1
            else
                call init_semi_stochastic(ss_space_in, tStartedFromCoreGround)
                if (tStartedFromCoreGround .and. tSetInitialRunRef) call set_initial_run_references()
            end if
        end if

        ! If the number of trial states to calculate hasn't been set by the
        ! user, then simply use the minimum number
        if ((tTrialWavefunction .or. tStartTrialLater) .and. (ntrial_ex_calc == 0)) then
            ntrial_ex_calc = inum_runs
        end if

        ! Initialise the trial wavefunction information which can be used for the energy estimator.
        ! This includes generating the trial space, generating the space connected to the trial space,
        ! diagonalising the trial space to find the trial wavefunction and calculating the vector
        ! in the connected space, required for the energy estimator.
        if (tRDMonFly .and. tDynamicCoreSpace .and. tTrialWavefunction) then
            tTrialWavefunction = .false.
            tStartTrialLater = .true.
            trial_shift_iter = 1
        end if
        if (tTrialWavefunction) then
            if (tPairedReplicas) then
                call init_trial_wf(trial_space_in, ntrial_ex_calc, inum_runs.div.2, .true.)
            else
                call init_trial_wf(trial_space_in, ntrial_ex_calc, inum_runs, .false.)
            end if
        else if (tStartTrialLater) then
            ! If we are going to turn on the use of a trial wave function
            ! later in the calculation, then zero the trial estimate arrays
            ! for now, to prevent junk being printed before then.
            trial_numerator = 0.0_dp
            tot_trial_numerator = 0.0_dp
            trial_denom = 0.0_dp
            tot_trial_denom = 0.0_dp

            init_trial_numerator = 0.0_dp
            tot_init_trial_numerator = 0.0_dp
            init_trial_denom = 0.0_dp
            tot_init_trial_denom = 0.0_dp

            trial_num_inst = 0.0_dp
            tot_trial_num_inst = 0.0_dp
            trial_denom_inst = 0.0_dp
            tot_trial_denom_inst = 0.0_dp
        end if

        replica_overlaps_real(:, :) = 0.0_dp
#ifdef CMPLX_
        replica_overlaps_imag(:, :) = 0.0_dp
#endif

        if (t_cepa_shift) call init_cepa_shifts()

        ! Set up the reference space for the adi-approach
        ! in real-time, we do this in the real-time init
        if (.not. t_real_time_fciqmc) call setup_reference_space(tReadPops)

        ! in fixed-n0, the variable shift mode and everything connected is
        ! controlled over the reference population
        if (tFixedN0) then
            if (tReadPops .and. .not. tWalkContGrow) then
                tSkipRef = .true.
                tSinglePartPhase = .false.
            else
                tSkipRef = .false.
                tSinglePartPhase = .true.
            end if
        end if

        if (tRDMonFly .and. tDynamicCoreSpace) call sync_rdm_sampling_iter()

        ! 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 (tInitiatorSpace) call init_initiator_space(i_space_in)

        if (tReplicaEstimates) then
            if (.not. tPairedReplicas) then
                call stop_all(this_routine, "The paired-replicas option must be used the logging &
                                            &block, in order to calculate replica estimates.)")
            end if

            if (tSemiStochastic) allocate (tDetermSpawnedTo( &
                                           cs_replicas(core_run)%determ_sizes(iProcIndex)))

            call open_replica_est_file()
        end if

        ! When should we perform death before communication?
        ! For integer walkers, do death before comms just so the tests don't fail (or need updating).
        if (.not. tAllRealCoeff) then
            tDeathBeforeComms = .true.
        end if
        if (t_back_spawn .or. t_back_spawn_flex) then
            tDeathBeforeComms = .true.
        end if

        ! For FCIQMC with preconditioning and a time step of 1, death will
        ! kill all walkers and remove them from the hash table. In this
        ! case, we must set the initiator flags for spawning to occupied
        ! determinants before this occurs.
        if (tPreCond .and. (tau.isclose.1.0_dp)) tSetInitFlagsBeforeDeath = .true.

        ! Make sure we are performing death *after* communication, in cases
        ! where this is essential.
        if (tPreCond .and. tDeathBeforeComms) then
            call stop_all(this_routine, "With preconditioning, death must &
                                        &be performed after communication.")
        end if
        if (tReplicaEstimates .and. tDeathBeforeComms) then
            call stop_all(this_routine, "In order to calculate replica estimates, &
                                        &death must be performed after communication.")
        end if
        if (tEN2Init .and. (.not. tTruncInitiator)) then
            call stop_all(this_routine, "Cannot calculate the EN2 correction to initiator &
                                        &error as the initiator method is not in use.")
        end if

    end subroutine InitFCIMCCalcPar

    subroutine init_fcimc_fn_pointers()
        character(*), parameter :: this_routine = "init_fcimc_fn_pointers"

        ! Almost all excitation generators in NECI are Full CI generators.
        gen_all_excits => gen_all_excits_default

        ! Select the excitation generator.
        if (tGAS) then
            call class_managed(generate_excitation, gen_all_excits)
        else if (t_3_body_excits .and. .not. (t_mol_3_body .or. t_ueg_3_body)) then
            if (t_uniform_excits) then
                generate_excitation => gen_excit_uniform_k_space_hub_transcorr
            else
                generate_excitation => gen_excit_k_space_hub_transcorr
            end if
        else if (t_ueg_3_body) then
            if (tTrcorrExgen) then
                generate_two_body_excitation => gen_ueg_excit
            else if (TLatticeGens) then
                generate_two_body_excitation => gen_rand_excit
            endif
            generate_excitation => gen_excit_mol_tc
        elseif(t_impurity_excitgen) then
            generate_excitation => gen_excit_impurity_model
        elseif ((t_back_spawn_option .or. t_back_spawn_flex_option)) then
            if (tHUB .and. tLatticeGens) then
                ! for now the hubbard + back-spawn still uses the old
                ! genrand excit gen
                generate_excitation => gen_excit_back_spawn_hubbard
            else if (tUEGNewGenerator .and. tLatticeGens) then
                generate_excitation => gen_excit_back_spawn_ueg_new
            else if (tUEG .and. tLatticeGens) then
                generate_excitation => gen_excit_back_spawn_ueg
            else
                generate_excitation => gen_excit_back_spawn
            end if
        else if (tUEGNewGenerator) then
            generate_excitation => gen_ueg_excit
        else if (tPickVirtUniform) then
            ! pick-uniform-random-mag is on
            if (tReltvy) then
                generate_excitation => gen_rand_excit_Ex_Mag
            else
                call stop_all(this_routine, "Excitation generator has not been set!")
            end if
        else if (tGenHelWeighted) then
            generate_excitation => gen_excit_hel_weighted
        else if (tGen_4ind_2) then
            generate_excitation => gen_excit_4ind_weighted2
        else if (tGen_4ind_weighted) then
            generate_excitation => gen_excit_4ind_weighted
        else if (tGen_4ind_reverse) then
            generate_excitation => gen_excit_4ind_reverse

        else if (tGUGA) then
            if (tgen_guga_crude) then
                if (t_k_space_hubbard) then
                    generate_excitation => gen_excit_k_space_hub
                else if (t_new_real_space_hubbard) then
                    generate_excitation => gen_excit_rs_hubbard
                else
                    generate_excitation => gen_excit_4ind_weighted2
                end if
            else
                generate_excitation => generate_excitation_guga
            end if

        else if (t_pcpp_excitgen) then
            generate_excitation => gen_rand_excit_pcpp
        else if (t_fci_pchb_excitgen) then
            call class_managed(generate_excitation, gen_all_excits)
        else if (t_k_space_hubbard) then
            if (t_3_body_excits) then
                if (t_uniform_excits) then
                    generate_excitation => gen_excit_uniform_k_space_hub_transcorr
                else if (t_mixed_excits) then
                    generate_excitation => gen_excit_mixed_k_space_hub_transcorr
                else
                    generate_excitation => gen_excit_k_space_hub_transcorr
                end if
            else
                if (t_uniform_excits) then
                    generate_excitation => gen_excit_uniform_k_space_hub
                else
                    generate_excitation => gen_excit_k_space_hub
                end if
            end if

        else if (t_new_real_space_hubbard) then
            if (t_trans_corr_hop) then
                if (t_hole_focus_excits) then
                    generate_excitation => gen_excit_rs_hubbard_transcorr_hole_focus
                else if (t_uniform_excits) then
                    generate_excitation => gen_excit_rs_hubbard_transcorr_uniform
                else
                    generate_excitation => gen_excit_rs_hubbard_transcorr
                end if
            else if (t_spin_dependent_transcorr) then
                generate_excitation => gen_excit_rs_hubbard_spin_dependent_transcorr
            else
                generate_excitation => gen_excit_rs_hubbard
            end if

        else if (t_tJ_model) then
            generate_excitation => gen_excit_tJ_model

        else if (t_heisenberg_model) then
            generate_excitation => gen_excit_heisenberg_model
        else
            generate_excitation => gen_rand_excit
        end if


        ! yes, fortran pointers work this way
            ! Pointer assignment with =>
                ! Fortran
                !> ptr1 => ptr2
                ! C
                !> T *ptr1, *ptr2;
                !> ptr1 = ptr2;
            ! Copy/value assignment with =
                ! Fortran
                !> ptr1 = ptr2
                ! C
                !> *ptr1 = *ptr2;

        ! if we are using the 3-body excitation generator, embed the chosen excitgen
        ! in the three-body one
        if (t_mol_3_body) then
            generate_two_body_excitation => generate_excitation
            generate_excitation => gen_excit_mol_tc
        end if
        ! Do the same for HPHF
        if (tHPHF) then
            exc_generator_for_HPHF => generate_excitation
            generate_excitation => gen_hphf_excit
        end if

        ! In the main loop, we only need to find out if a determinant is
        ! connected to the reference det or not (so no ex. level above 2 is
        ! required). Except in some cases where we need to know the maximum
        ! excitation level
        if (tTruncSpace .or. tHistSpawn .or. tCalcFCIMCPsi) then
            max_calc_ex_level = nel
        else
            if (t_3_body_excits) then
                max_calc_ex_level = 3
            else
                max_calc_ex_level = 2
            end if
        end if

        ! How many children should we spawn given an excitation?
        if (t_real_time_fciqmc) then
            attempt_create => attempt_create_realtime
        else if (tTruncCas .or. tTruncSpace .or. &
                 tPartFreezeCore .or. tPartFreezeVirt .or. tFixLz .or. &
                 (tUEG .and. .not. tLatticeGens) .or. tTruncNOpen .or. t_trunc_nopen_diff) then
            if (tHPHF .or. tSemiStochastic) then
                attempt_create => attempt_create_trunc_spawn
            else
                attempt_create => att_create_trunc_spawn_enc
            end if
        else
            attempt_create => attempt_create_normal
        end if

        ! In attempt create, how should we evaluate the off diagonal matrix
        ! elements between a parent and its (potentially) spawned offspring?
        if (tHPHF) then
            if (tGenMatHEL) then
                get_spawn_helement => hphf_spawn_sign
            else
                get_spawn_helement => hphf_off_diag_helement_spawn
            end if

            ! new guga addition: do not need to recalculate Helement
        else if (tGUGA) then
            ! use hphf_routine also, since it does exactly what needed
            get_spawn_helement => hphf_spawn_sign
        else
            get_spawn_helement => get_helement_det_only
        end if

        ! When calling routines to generate all possible connections, this
        ! routine is called to generate the corresponding Hamiltonian matrix
        ! elements.
        if (tHPHF) then
            get_conn_helement => hphf_off_diag_helement_spawn
        else
            get_conn_helement => get_helement_det_only
        end if

        ! Once we have generated the children, do we need to encode them?
        if (.not. (tHPHF .or. tGen_4ind_weighted .or. tGUGA)) then
            encode_child => FindExcitBitDet
        else
            encode_child => null_encode_child
        end if

        ! What message should we display for a particle bloom?
        if (tAddToInitiator) then
            bloom_warn_string = '("Bloom of more than n_add on ", a, " excit: &
                                &A max of ", f10.2, " particles created. ", &
                                &i8, " blooms occurred.")'
        else
            ! Use this variable to store the bloom cutoff level.
            InitiatorWalkNo = 3.0_dp
            bloom_warn_string = '("Bloom of more than 3 on ", a, " excit: &
                                &A max of ", f10.2, " particles created. ", &
                                &i8, " blooms occurred.")'
        end if
        bloom_max = 0

        if (tPreCond) then
            attempt_die => attempt_die_precond
        else
            attempt_die => attempt_die_normal
        end if

        extract_bit_rep_avsign => extract_bit_rep_avsign_no_rdm

        fill_rdm_diag_currdet => fill_rdm_diag_currdet_norm

        select case (sfTag)
        case (0)
            scaleFunction => powerScaleFunction
        case (1)
            scaleFunction => expScaleFunction
        case (2)
            scaleFunction => negScaleFunction
        case (3)
            scaleFunction => expCOScaleFunction
        case default
            call stop_all(this_routine, "Invalid scale function specified")
        end select

        ! if (tExpAdaptiveShift) then
        !     shiftFactorFunction => expShiftFactorFunction
        if (tLinearAdaptiveShift) then
            shiftFactorFunction => linearShiftFactorFunction
        else if (tAutoAdaptiveShift) then
            shiftFactorFunction => autoShiftFactorFunction
        else
            shiftFactorFunction => constShiftFactorFunction
        end if

        ! select the procedure that returns all connected determinants.
        if (t_k_space_hubbard) then
            gen_all_excits => gen_all_excits_k_space_hubbard
        else if (t_new_real_space_hubbard) then
            gen_all_excits => gen_all_excits_r_space_hubbard
        end if

    end subroutine init_fcimc_fn_pointers

    subroutine DeallocFCIMCMemPar()

        CHARACTER(len=*), PARAMETER :: this_routine = 'DeallocFciMCMemPar'
        type(ll_node), pointer :: Curr, Prev
        integer :: i, ierr

        deallocate(RandomHash2, stat=ierr)
        if (ierr /= 0) call stop_all(this_routine, "Err deallocating")

        ! Deallocate the linked list
        do i = 1, nWalkerHashes
            Curr => HashIndex(i)%Next
            Prev => HashIndex(i)
            nullify (Prev%Next)
            do while (associated(Curr))
                Prev => Curr
                Curr => Curr%Next
                deallocate(Prev)
                if (ierr /= 0) call stop_all(this_routine, "Err deallocating")
            end do
        end do
        deallocate(HashIndex, stat=ierr)
        if (ierr /= 0) call stop_all(this_routine, "Err deallocating")
        nullify (Curr)
        nullify (Prev)

        deallocate(FreeSlot, stat=ierr)
        if (ierr /= 0) call stop_all(this_routine, "Err deallocating")

        IF (tHistSpawn .or. tCalcFCIMCPsi) THEN
            DEallocate(Histogram)
            DEallocate(AllHistogram)
            IF (tHistSpawn) THEN
                DEallocate(InstHist)
                DEallocate(InstAnnihil)
                DEallocate(AvAnnihil)
            end if
            IF (iProcIndex == 0) THEN
                IF (tHistSpawn) THEN
                    DEallocate(AllInstHist)
                    DEallocate(AllAvAnnihil)
                    DEallocate(AllInstAnnihil)
                end if
            end if
        else if (tHistEnergies) THEN
            DEallocate(HistogramEnergy)
            DEallocate(AttemptHist)
            DEallocate(SpawnHist)
            DEallocate(SinglesHist)
            DEallocate(DoublesHist)
            DEallocate(DoublesAttemptHist)
            DEallocate(SinglesAttemptHist)
            DEallocate(SinglesHistOccOcc)
            DEallocate(SinglesHistVirtOcc)
            DEallocate(SinglesHistOccVirt)
            DEallocate(SinglesHistVirtVirt)
            IF (iProcIndex == Root) THEN
                DEallocate(AllHistogramEnergy)
                DEallocate(AllAttemptHist)
                DEallocate(AllSpawnHist)
                DEallocate(AllSinglesAttemptHist)
                DEallocate(AllSinglesHist)
                DEallocate(AllDoublesAttemptHist)
                DEallocate(AllDoublesHist)
                DEallocate(AllSinglesHistOccOcc)
                DEallocate(AllSinglesHistVirtOcc)
                DEallocate(AllSinglesHistOccVirt)
                DEallocate(AllSinglesHistVirtVirt)
            end if
        end if
        if (tHistExcitToFrom) call clean_hist_excit_tofrom()
        DEallocate(WalkVecDets)
        CALL LogMemDealloc(this_routine, WalkVecDetsTag)
        DEallocate(SpawnVec)
        CALL LogMemDealloc(this_routine, SpawnVecTag)
        DEallocate(SpawnVec2)
        CALL LogMemDealloc(this_routine, SpawnVec2Tag)
        if (tAutoAdaptiveShift) then
            DEallocate(SpawnInfoVec)
            CALL LogMemDealloc(this_routine, SpawnInfoVecTag)
            DEallocate(SpawnInfoVec2)
            CALL LogMemDealloc(this_routine, SpawnInfoVec2Tag)
        end if

        if (allocated(TempSpawnedParts)) then
            deallocate(TempSpawnedParts)
            log_dealloc(TempSpawnedPartsTag)
        end if
        DEallocate(HFDet)
        CALL LogMemDealloc(this_routine, HFDetTag)
        DEallocate(iLutHF)
        DEallocate(iLutRef)
        DEallocate(ProjEDet)
        DEallocate(iLutHF_True)
        DEallocate(HFDet_True)
        IF (ALLOCATED(HighestPopDet)) DEallocate(HighestPopDet)
        IF (ALLOCATED(RandomOrbIndex)) DEallocate(RandomOrbIndex)

        IF (ALLOCATED(SpinInvBrr)) THEN
            CALL LogMemDealloc(this_routine, SpinInvBRRTag)
            DEallocate(SpinInvBRR)
        end if
        IF (ALLOCATED(CoreMask)) THEN
            DEallocate(CoreMask)
            DEallocate(CASMask)
        end if
        IF (tPrintOrbOcc) THEN
            DEallocate(OrbOccs)
            CALL LogMemDeAlloc(this_routine, OrbOccsTag)
        end if

        IF (tHistInitPops) THEN
            if (allocated(HistInitPops)) then
                deallocate(HistInitPops)
                call LogMemDeAlloc(this_routine, HistInitPopsTag)
            end if
            IF (iProcIndex == 0) THEN
                if (allocated(AllHistInitPops)) then
                    deallocate(AllHistInitPops)
                    call LogMemDeAlloc(this_routine, AllHistInitPopsTag)
                end if
            end if
        end if

        if (tHub) then
            if (allocated(momIndexTable)) deallocate(momIndexTable)
            deallocate(breathingCont)
        end if

        if (tRDMonFly) call dealloc_global_rdm_data()

        if (allocated(refdetflip)) deallocate(refdetflip)
        if (allocated(ilutrefflip)) deallocate(ilutrefflip)
        if (allocated(ValidSpawnedList)) deallocate(ValidSpawnedList)
        if (allocated(InitialSpawnedSlots)) deallocate(InitialSpawnedSlots)

        ! Cleanup global storage
        call clean_global_det_data()

        ! Cleanup excitation generation storage
        call clean_excit_gen_store(fcimc_excit_gen_store)

        ! Cleanup cont time
        call clean_cont_time()

        ! Cleanup the load balancing
        call clean_load_balance()

        ! Cleanup adi caches
        call clean_adi()


        ! Cleanup excitation generator
        if (t_guga_pchb) then
            call finalize_pchb_excitgen_guga()
        end if

        if (t_pcpp_excitgen) call finalize_pcpp_excitgen()

        if(t_impurity_excitgen) call clearImpurityExcitgen()

        if (tSemiStochastic) call end_semistoch()

        if (tTrialWavefunction) call end_trial_wf()

        call finalize_exz_gen_class()


    end subroutine DeallocFCIMCMemPar

    subroutine InitFCIMC_CSF()
        implicit none
        integer(n_int), allocatable ::  initSpace(:, :)
        integer :: count, nUp, nOpen
        integer :: i, j, lwork, proc
        integer :: DetHash, pos, TotWalkersTmp, nI(nel), nJ(nel)
        integer(n_int) :: ilutJ(0:NIfTot)
        integer(n_int), allocatable :: openSubspace(:)
        real(dp), allocatable :: S2(:, :), eigsImag(:), eigs(:), evs(:, :), void(:, :), work(:)
        real(dp) :: normalization, rawWeight, HDiag, tmpSgn(lenof_sign)
        HElement_t(dp) :: HOffDiag
        integer :: err
        real(dp) :: HFWeight(inum_runs)
        logical :: tSuccess
        character(*), parameter :: t_r = "InitFCIMC_CSF"

        ! get the number of open orbitals
        nOpen = sum(nOccOrbs) - sum(nClosedOrbs)
        ! in a closed shell system, nothing to do
        if (nOpen == 0) then
            call InitFCIMC_HF()
            return
        end if
        ! first, set up the space considered for the CSF
        call generateInitSpace()
        if (allocated(openSubspace)) deallocate(openSubspace)
        ! we now have initSpace(:,:) with iluts belonging to all possible initial
        ! dets (i.e. all dets contributing to the target CSF) -> construct S2
        allocate(S2(count, count))
        do i = 1, count
            do j = 1, count
                S2(i, j) = S2Matel(initSpace(:, i), initSpace(:, j))
            end do
        end do

        ! prepare the diagonalization
        allocate(eigs(count))
        allocate(evs(count, count))
        allocate(eigsImag(count))
        allocate(work(1))
        allocate(void(0, 0))
        ! workspace query, get how much tmp memory we need
        call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, -1, err)
        ! allocate work array
        lwork = int(work(1))
        deallocate(work)
        allocate(work(lwork))
        ! diagonalize S2
        call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, lwork, err)
        deallocate(void)
        deallocate(work)
        deallocate(eigsImag)

        ! transfer the eigenvector
        do i = 1, count
            if (abs(S2Init * (S2Init + 1) - eigs(i)) < eps) exit
        end do
        if (i > count) then
            call stop_all(t_r, "Requested S2 eigenvalue does not exist")
        end if
        eigs = evs(:, i)
!      normalization = minval(eigs)
        rawWeight = sum(abs(eigs))
        normalization = InitialPart / rawWeight

!      refLoc = maxloc(eigs)
        eigs = eigs * normalization

        TotWalkers = 0
        iStartFreeSlot = 1
        iEndFreeSlot = 0
        do i = 1, count
            call decode_bit_det(nI, initSpace(:, i))
            proc = DetermineDetNode(nel, nI, 0)
            if (iProcIndex == proc) then
                HDiag = get_diagonal_matel(nI, initSpace(:, i))
                HOffDiag = get_off_diagonal_matel(nI, initSpace(:, i))
                DetHash = FindWalkerHash(nI, size(HashIndex))
                TotWalkersTmp = int(TotWalkers)
                tmpSgn = eigs(i)
                call encode_sign(initSpace(:, i), tmpSgn)
                if (tHPHF) then
                    call FindDetSpinSym(nI, nJ, nel)
                    call encodebitdet(nJ, ilutJ)
                    ! if initSpace(:,i) is not the det of the HPHF pair we are storing,
                    ! skip this - the correct contribution will be stored once
                    ! the spin-flipped version is stored
                    if (DetBitLT(initSpace(:, i), ilutJ, NIfD) == 1) cycle
                end if
                call AddNewHashDet(TotWalkersTmp, initSpace(:, i), DetHash, nI, HDiag, HOffDiag, pos, err)
                TotWalkers = TotWalkersTmp
            end if
            ! reset the reference?
        end do

        call hash_table_lookup(HFDet, ilutHF, nifd, HashIndex, CurrentDets, i, DetHash, tSuccess)
        if (tSuccess) then
            call extract_sign(CurrentDets(:, i), tmpSgn)
            do i = 1, inum_runs
                HFWeight(i) = mag_of_run(tmpSgn, i)
            end do
        else
            HFWeight = 0.0_dp
        end if

        AllTotParts = InitialPart
        AllTotPartsOld = InitialPart
        OldAllAvWalkersCyc = InitialPart
        OldAllHFCyc = HFWeight
        OldAllNoatHF = HFWeight

        ! cleanup
        deallocate(evs)
        deallocate(eigs)
        deallocate(S2)
        deallocate(initSpace)
    contains

        !------------------------------------------------------------------------------------------!

        subroutine generateOpenOrbIluts()
            use IntegralsData, only: nfrozen
            implicit none

            count = 0
            nUp = (nel + lms) / 2 - sum(nClosedOrbs) + nfrozen / 2
            do i = 1, 2**nOpen - 1
                if (popcnt(i) == nUp) then
                    count = count + 1
                    if (allocated(openSubspace)) openSubspace(count) = i
                end if
            end do
        end subroutine generateOpenOrbIluts

        !------------------------------------------------------------------------------------------!

        subroutine generateInitSpace()
            implicit none

            integer :: nI(nel), nIBase(nel)
            integer :: iEl, iElBase, iOpen
            integer :: openOrbList(nel)
            logical :: previousCont, nextCont, tClosed
            character(*), parameter :: t_r = "generateInitSpace"

            ! create a list of all open-shell determinants with the correct spin+orbs
            nUp = (nel + lms) / 2
            call generateOpenOrbIluts()
            allocate(openSubspace(count))
            call generateOpenOrbIluts()

            ! convert open-shell-only iluts into full iluts
            ! use the reference to determine which orbitals shall participate
            iElBase = 1
            nIBase = 0
            ! generate a list of open orbitals
            iOpen = 0
            openOrbList = 0
            do i = 1, nel
                if (i == 1) then
                    previousCont = .false.
                else
                    previousCont = FDet(i - 1) == FDet(i) - 1
                end if
                if (i == nel) then
                    nextCont = .false.
                else
                    nextCont = FDet(i + 1) == FDet(i) + 1
                end if
                tClosed = .true.
                ! identify open orbitals, using the ordering: does the next one belong
                ! to the same spatial orb?
                if (is_beta(FDet(i)) .and. .not. nextCont) then
                    iOpen = iOpen + 1
                    openOrbList(iOpen) = FDet(i)
                    tClosed = .false.
                end if
                ! or the previous one?
                if (is_alpha(FDet(i)) .and. .not. previousCont) then
                    iOpen = iOpen + 1
                    openOrbList(iOpen) = FDet(i)
                    tClosed = .false.
                end if
                ! if the orbital is not open, it is closed
                if (tClosed) then
                    nIBase(iElBase) = FDet(i)
                    iElBase = iElBase + 1
                end if
            end do
            if (iOpen /= nOpen) then
                write(stderr, *) "nOpen/iOpen conflict", FDet, openOrbList, iOpen, nOpen
                call stop_all(t_r, "Error in determining open shell orbitals")
            end if

            allocate(initSpace(0:NIfTot, count))
            ! now, add the open-shell contribution
            do i = 1, count
                ! start from the closed-shell base
                nI = nIBase
                iEl = iElBase
                ! for each open orb, add the electron
                do j = 1, nOpen
                    ! based on openSubspace(i), we are considering alpha/beta for this
                    ! orb
                    if (btest(openSubspace(i), j - 1)) then
                        nI(iEl) = get_alpha(openOrbList(j))
                    else
                        nI(iEl) = get_beta(openOrbList(j))
                    end if
                    iEl = iEl + 1
                end do
                ! encode the determinant to the initial space
                call sort(nI)
                call EncodeBitDet(nI, initSpace(:, i))
            end do

        end subroutine generateInitSpace

        function S2Matel(ilutA, ilutB) result(matel)
            implicit none
            integer(n_int), intent(in) :: ilutA(0:NIfTot), ilutB(0:NIfTot)
            integer(n_int) :: splus(0:NIfTot), sminus(0:NIfTot)
            real(dp) :: matel
            integer :: k, m, nI(nel), upOrb, downOrb

            matel = 0.0_dp
            if (DetBitEq(ilutA, ilutB, NIfD)) then
                matel = matel + real(lms * (lms + 2), dp) / 4.0_dp
            end if
            ! get the offdiag part of S2: S-S+
            call decode_bit_det(nI, ilutA)
            do k = 1, nel
                ! first, apply S- to all electrons (additively)
                if (is_beta(nI(k))) then
                    sminus = ilutA
                    downOrb = get_alpha(nI(k))
                    clr_orb(sminus, nI(k))
                    set_orb(sminus, downOrb)
                    ! now, apply S+ to all electrons (again, sum)
                    do m = 1, nel
                        ! check if it yields 0
                        if (is_alpha(nI(m)) .or. m == k) then
                            splus = sminus
                            ! if not, go on
                            if (m == k) then
                                clr_orb(splus, downOrb)
                            else
                                clr_orb(splus, nI(m))
                            end if
                            upOrb = get_beta(nI(m))
                            set_orb(splus, upOrb)
                            if (DetBitEq(splus, ilutB, NIFD)) matel = matel + 1.0_dp
                        end if
                    end do
                end if
            end do
        end function S2Matel
    end subroutine InitFCIMC_CSF

    !------------------------------------------------------------------------------------------!

    subroutine InitFCIMC_HF()

        integer :: run, DetHash
        real(dp), dimension(lenof_sign) :: InitialSign
        HElement_t(dp) :: h_temp

        if (tOrthogonaliseReplicas) then
            call InitFCIMC_HF_orthog()
            return
        end if

        InitialPartVec = 0.0_dp
        do run = 1, inum_runs
            InitialPartVec(min_part_type(run)) = InitialPart
#ifdef CMPLX_
            InitialPartVec(max_part_type(run)) = 0.0_dp
#endif
        end do

        !Setup initial walker local variables for HF walkers start
        IF (iProcIndex == iRefProc(1)) THEN

            ! Encode the reference determinant identification.
            call encode_det(CurrentDets(:, 1), iLutHF)

            !Point at the correct position for the first walker
            DetHash = FindWalkerHash(HFDet, nWalkerHashes)    !Find det hash position
            HashIndex(DetHash)%Ind = 1

            ! Clear the flags
            call clear_all_flags(CurrentDets(:, 1))

            ! Set reference determinant as an initiator if
            ! tTruncInitiator is set, for both imaginary and real flags
            ! in real-time calculations, the reference does not have any special role
            if (tTruncInitiator) then
                do run = 1, inum_runs
                    call set_flag(CurrentDets(:, 1), get_initiator_flag_by_run(run))
                end do
            end if

            ! If running a semi-stochastic simulation, set flag to specify the Hartree-Fock is in the
            ! deterministic space.
            if (tSemiStochastic) then
                do run = 1, inum_runs
                    call set_flag(CurrentDets(:, 1), flag_deterministic(run))
                end do
            end if

            ! if no reference energy is used, explicitly get the HF energy
            if (tZeroRef) then
                h_temp = get_diagonal_matel(HFDet, ilutHF)
            else
                ! HF energy is equal to 0 (when used as reference energy)
                h_temp = h_cast(0.0_dp)
            end if
            call set_det_diagH(1, real(h_temp, dp))
            call set_det_offdiagH(1, h_cast(0.0_dp))
            HFInd = 1

            call store_decoding(1, HFDet)

            if (associated(lookup_supergroup_indexer)) then
                call set_supergroup_idx(1, lookup_supergroup_indexer%idx_nI(HFDet))
            end if

            if (tContTimeFCIMC .and. tContTimeFull) &
                call set_spawn_rate(1, spawn_rate_full(HFDet, ilutHF))

            ! Obtain the initial sign
            InitialSign = 0.0_dp
            if (tStartSinglePart) then
                InitialSign(:) = InitialPartVec(:)
                TotParts(:) = InitialPartVec(:)
                TotPartsOld(:) = InitialPartVec(:)
            else
                do run = 1, inum_runs
                    InitialSign(min_part_type(run)) = InitWalkers
                    TotParts(min_part_type(run)) = real(InitWalkers, dp)
                    TotPartsOld(min_part_type(run)) = real(InitWalkers, dp)
#ifdef CMPLX_
                    TotParts(max_part_type(run)) = 0.0_dp
                    TotPartsOld(max_part_type(run)) = 0.0_dp
#endif
                end do
            end if

            ! set initial values for global control variables.

            TotWalkers = 1
            TotWalkersOld = 1
            NoatHF(:) = InitialSign(:)
            call encode_sign(CurrentDets(:, 1), InitialSign)
        ELSE
            NoatHF(:) = 0.0_dp
            TotWalkers = 0
            TotWalkersOld = 0
        end if

        OldAllNoatHF(:) = 0.0_dp
        AllNoatHF(:) = 0.0_dp
        IF (TStartSinglePart) THEN
            !Initialise global variables for calculation on the root node
            IF (iProcIndex == root) THEN
                OldAllNoatHF = InitialPartVec
                do run = 1, inum_runs
                    OldAllAvWalkersCyc(run) = sum(InitialPartVec( &
                                                  min_part_type(run):max_part_type(run)))
                end do
                AllNoatHF = InitialPartVec
                InstNoatHF = InitialPartVec
                AllTotParts = InitialPartVec
                AllTotPartsOld = InitialPartVec
                AllNoAbortedOld(:) = 0.0_dp
                iter_data_fciqmc%tot_parts_old = InitialPartVec
                AllTotWalkers = 1
                AllTotWalkersOld = 1
                do run = 1, inum_runs
                    OldAllHFCyc(run) = ARR_RE_OR_CPLX(InitialPartVec, run)
                end do
            end if
        ELSE
            !In this, only one processor has initial particles.
            IF (iProcIndex == Root) THEN
                AllTotWalkers = 1
                AllTotWalkersOld = 1
                do run = 1, inum_runs
                    iter_data_fciqmc%tot_parts_old(run) = real(InitWalkers, dp)
                    AllTotParts(run) = InitWalkers
                    AllTotPartsOld(run) = InitWalkers
                    AllNoAbortedOld(run) = 0.0_dp
                end do
            end if
        end if

    end subroutine InitFCIMC_HF

    subroutine InitFCIMC_HF_orthog()

        ! This is a reimplementation of InitFCIMC_HF to work with multiple
        ! different reference states.
        !
        ! In the end, we expect to be able to just substitute it back in
        ! for the original. It should give the same results
        ! TODO: Substitute back

        integer :: run, site, hash_val, i
        logical :: repeated
        HElement_t(dp) :: hdiag
        character(*), parameter :: this_routine = 'InitFCIMC_HF_orthog'

        ! Add some implementation guards

        ! Default values, unless overridder for individual procs
        NoatHF = 0.0_dp
        TotWalkers = 0
        TotWalkersOld = 0
        tRef_Not_HF = .true.
        tNoBrillouin = .true.

        !
        ! Initialise each of the runs separately.
        site = 0
        do run = 1, inum_runs

            ! If this run should have the reference on this site, then
            ! initialise it as appropriate.
            if (iProcIndex == iRefProc(run)) then

                ! Check if this reference is the same as any of the previous
                ! ones. If it is not, then at the end of the loop (i == site+1)
                repeated = .false.
                do i = 1, site
                    if (DetBitEQ(CurrentDets(:, i), ilutRef(:, run))) then
                        repeated = .true.
                        exit
                    end if
                end do
                site = i

                if (.not. repeated) then
                    ! Add the site to the main list (unless it is already there)
                    call encode_det(CurrentDets(:, site), ilutRef(:, run))
                    hash_val = FindWalkerHash(ProjEDet(:, run), nWalkerHashes)
                    call add_hash_table_entry(HashIndex, site, hash_val)

                    ! Clear all the flags and sign
                    call clear_all_flags(CurrentDets(:, site))
                    call nullify_ilut(CurrentDets(:, site))
                end if

                ! Set reference determinant as an initiator if tTruncInitiator
                if (tTruncInitiator) then
                    call set_flag(CurrentDets(:, site), get_initiator_flag_by_run(run))
                end if

                ! The global reference is the HF and is primary for printed
                ! energies.
                if (run == 1) HFInd = site
                hdiag = get_diagonal_matel(ProjEDet(:, run), ilutRef(:, run))
                call set_det_diagH(site, real(hdiag, dp) - Hii)
                call set_det_offdiagH(site, h_cast(0_dp))

                if (associated(lookup_supergroup_indexer)) then
                    call set_supergroup_idx(site, lookup_supergroup_indexer%idx_nI(ProjEDet(:, run)))
                end if

                ! store the determinant
                call store_decoding(site, ProjEDet(:, run))

                ! Obtain the initial sign
                if (.not. tStartSinglePart) &
                    call stop_all(this_routine, "Only startsinglepart supported")
                call encode_part_sign(CurrentDets(:, site), InitialPart, min_part_type(run))

                ! Initial control values
                TotWalkers = site
                TotWalkersOld = site
                NoatHF(min_part_type(run)) = InitialPart
                TotParts(min_part_type(run)) = real(InitialPart, dp)
                TotPartsOld(min_part_type(run)) = real(InitialPart, dp)
            end if
        end do

        ! Check to ensure that the following code is valid
        if (.not. tStartSinglePart) &
            call stop_all(this_routine, "Only startsinglepart supported")

        ! Initialise global variabes for calculation on the root node
        OldAllNoatHF = 0.0_dp
        AllNoatHF = 0.0_dp
        call MPISum(TotWalkers, AllTotWalkers)
        if (iProcIndex == root) then
            OldAllNoatHF(:) = InitialPart
            OldAllAvWalkersCyc(:) = InitialPart
            AllNoatHF(:) = InitialPart
            InstNoatHF(:) = InitialPart
            AllTotParts(:) = InitialPart
            AllTotPartsOld(:) = InitialPart
            AllNoAbortedOld(:) = InitialPart
            OldAllHFCyc(:) = InitialPart

            TotWalkersOld = TotWalkers
        end if

    end subroutine InitFCIMC_HF_orthog

    subroutine InitFCIMC_trial()

        ! Use the code generated for the KPFCIQMC excited state calculations
        ! to initialise the FCIQMC simulation.

!         integer :: nexcit, ndets_this_proc, i, det(nel)
        integer :: nexcit, ndets_this_proc, det(nel)
        integer(int64) :: i

        type(basisfn) :: sym
        real(dp) :: evals(inum_runs / nreplicas)
        HElement_t(dp), allocatable :: evecs_this_proc(:, :)
        integer(MPIArg) :: space_sizes(0:nProcessors - 1), space_displs(0:nProcessors - 1)
        character(*), parameter :: this_routine = 'InitFCIMC_trial'

        nexcit = inum_runs / nreplicas

        ! Create the trial excited states
        if (allocated(trial_init_reorder)) then
            call calc_trial_states_lanczos(init_trial_in, nexcit, ndets_this_proc, &
                                           SpawnedParts, evecs_this_proc, evals, &
                                           space_sizes, space_displs, trial_init_reorder)
        else
            call calc_trial_states_lanczos(init_trial_in, nexcit, ndets_this_proc, &
                                           SpawnedParts, evecs_this_proc, evals, &
                                           space_sizes, space_displs)
        end if
        ! Determine the walker populations associated with these states
        call set_trial_populations(nexcit, ndets_this_proc, evecs_this_proc)
        ! Set the trial excited states as the FCIQMC wave functions
        call set_trial_states(ndets_this_proc, evecs_this_proc, SpawnedParts, &
                              .false., tPairedReplicas)

        deallocate(evecs_this_proc)

        if (tSetInitialRunRef) call set_initial_run_references()

        ! Add an initialisation check on symmetries.
        if ((.not. tHub) .and. (.not. tUEG)) then
            do i = 1, TotWalkers
                call decode_bit_det(det, CurrentDets(:, i))
                call getsym_wrapper(det, sym)
                if (sym%sym%S /= HFSym%sym%S .or. sym%ml /= HFSym%Ml) &
                    call stop_all(this_routine, "Invalid det found")
            end do
        end if

    end subroutine InitFCIMC_trial

    subroutine set_initial_run_references()

        ! Analyse each of the runs, and set the reference determinant to the
        ! det with the largest coefficient, rather than the currently guessed
        ! one...

        HElement_t(dp) :: largest_coeff, sgn
        integer(n_int) :: largest_det(0:NIfTot)
!         integer :: run, j
        integer(int64) :: j
        integer :: run
        integer(int32) :: proc_highest
        integer(n_int) :: ilut(0:NIfTot)
        integer(int32) :: int_tmp(2)

        do run = 1, inum_runs

            if (tMultipleInitialRefs) then
                ! Use user specified reference states.
                call EncodeBitDet(initial_refs(:, run), ilut)
                call update_run_reference(ilut, run)
            else
                ! Find the largest det on this processor
                largest_coeff = h_cast(0.0_dp)
                do j = 1, TotWalkers
                    sgn = extract_run_sign(CurrentDets(:, j), run)
                    if (abs(sgn) > abs(largest_coeff)) then
                        largest_coeff = sgn
                        largest_det = CurrentDets(:, j)
                    end if
                end do

                ! Find the largest det on any processor (n.b. discard the
                ! non-integer part. This isn't all that important).
                ! [W.D. 15.5.2017:]
                ! for the test suite problems, maybe it is important..
                ! because there seems to be some compiler dependent
                ! differences..
                call MPIAllReduceDatatype( &
                    (/int(abs(largest_coeff), int32), int(iProcIndex, int32)/), 1, &
                    MPI_MAXLOC, MPI_2INTEGER, int_tmp)
                proc_highest = int_tmp(2)
                call MPIBCast(largest_det, NIfTot + 1, int(proc_highest))
!                 call MPIBCast(largest_det, NIfTot+1, proc_highest)

                write(stdout, *) 'Setting ref', run
                call writebitdet(stdout, largest_det, .true.)

                ! Set this det as the reference
                call update_run_reference(largest_det, run)

            end if
        end do

        if (tMultiRefShift) then
            do run = 1, inum_runs
                DiagSft(run) = proje_ref_energy_offsets(run)
            end do
        end if

    end subroutine set_initial_run_references

    !Routine to initialise the particle distribution according to the MP1 wavefunction.
    !This hopefully will help with close-lying excited states of the same sym.
    subroutine InitFCIMC_MP1()
        real(dp) :: TotMP1Weight, amp, MP2Energy, PartFac, rat, r, energy_contrib
        HElement_t(dp) :: HDiagtemp, HOffDiagtemp
        integer :: iExcits, exflag, Ex(2, maxExcit), nJ(NEl), DetIndex, iNode
        integer :: iInit, DetHash, ExcitLevel, run, part_type
        integer(n_int) :: iLutnJ(0:NIfTot)
        real(dp) :: NoWalkers, temp_sign(lenof_sign)
        logical :: tAllExcitsFound, tParity
        type(ll_node), pointer :: TempNode
        character(len=*), parameter :: this_routine = "InitFCIMC_MP1"
        integer(n_int) :: ilutG(0:nifguga)
        integer(n_int), allocatable :: excitations(:, :)
        integer :: i

#ifdef CMPLX_
        call stop_all(this_routine, "StartMP1 currently does not work with complex walkers")
#endif
        if (tReadPops) call stop_all(this_routine, "StartMP1 cannot work with with ReadPops")
        if (tStartSinglePart) call stop_all(this_routine, "StartMP1 cannot work with StartSinglePart")
        if (tRestartHighPop) call stop_all(this_routine, "StartMP1 cannot with with dynamically restarting calculations")

        write(stdout, *) "Initialising walkers proportional to the MP1 amplitudes..."

        if (tHPHF) then
            if (.not. TestClosedShellDet(iLutHF)) then
                call stop_all(this_routine, "Cannot use HPHF with StartMP1 if your reference is open-shell")
            end if
        end if

        if (tUEG) then
            !Parallel N^2M implementation of MP2 for UEG
            call CalcUEGMP2()
        end if

        !First, calculate the total weight - TotMP1Weight
        mp2energy = 0.0_dp
        TotMP1Weight = 1.0_dp
        iExcits = 0
        tAllExcitsFound = .false.
        if (tUEG) then
            exflag = 2
        else
            exflag = 3
        end if
        Ex(:, :) = 0
        ! for GUGA, this whole excitation generation is different!
        if (tGUGA) then
            ! should finally do some general routine, which does all this
            ! below...
            call convert_ilut_toGUGA(ilutHF, ilutG)
            call actHamiltonian(ilutG, CSF_Info_t(ilutG), excitations, iExcits)
            do i = 1, iExcits
                call convert_ilut_toNECI(excitations(:, i), ilutnJ)
                call decode_bit_det(nJ, iLutnJ)
                call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
                TotMP1Weight = TotMP1Weight + abs(amp)
                MP2Energy = MP2Energy + energy_contrib
            end do

        else
            do while (.true.)
                call GenExcitations3(HFDet, iLutHF, nJ, exflag, Ex, tParity, tAllExcitsFound, .false.)
                if (tAllExcitsFound) exit !All excits found

                call EncodeBitDet(nJ, iLutnJ)
                if (tHPHF) then
                    !Working in HPHF Space. Check whether determinant generated is an 'HPHF'
                    if (.not. IsAllowedHPHF(iLutnJ)) cycle
                end if
                iExcits = iExcits + 1

                call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
                TotMP1Weight = TotMP1Weight + abs(amp)
                MP2Energy = MP2Energy + energy_contrib
            end do
        end if

        if ((.not. tHPHF .and. .not. tGUGA) .and. (iExcits /= (nDoubles + nSingles))) then
            write(stderr, *) nDoubles, nSingles, iExcits
            call stop_all(this_routine, "Not all excitations accounted for in StartMP1")
        end if

        write(stdout, "(A,2G25.15)") "MP2 energy calculated: ", MP2Energy, MP2Energy + Hii

        if ((InitialPart.isclose.1._dp) .or. (InitialPart >= (InitWalkers * nNodes) - 50)) then
            !Here, all the walkers will be assigned to the MP1 wavefunction.
            !InitialPart = 1 by default
            write(stdout, "(A)") "All walkers specified in input will be distributed according to the MP1 wavefunction."
            write(stdout, "(A)") "Shift will be allowed to vary from the beginning"
            write(stdout, "(A)") "Setting initial shift to equal MP2 correlation energy"
            DiagSft = MP2Energy
            !PartFac is the number of walkers that should reside on the HF determinant
            !in an intermediate normalised MP1 wavefunction.
            PartFac = (real(InitWalkers, dp) * real(nNodes, dp)) / TotMP1Weight
        else
            !Here, not all walkers allowed will be initialised to the MP1 wavefunction.
            write(stdout, "(A,G15.5,A)") "Initialising ", InitialPart, " walkers according to the MP1 distribution."
            write(stdout, "(A,G15.5)") "Shift will remain fixed until the walker population reaches ", InitWalkers * nNodes
            !PartFac is the number of walkers that should reside on the HF determinant
            !in an intermediate normalised MP1 wavefunction.
            PartFac = real(InitialPart, dp) / TotMP1Weight
            tSinglePartPhase(:) = .true.
        end if

        !Now generate all excitations again, creating the required number of walkers on each one.
        ! puh... for GUGA this gets messy to change below.. damn

        DetIndex = 1
        TotParts = 0.0
        tAllExcitsFound = .false.
        if (tUEG) then
            exflag = 2
        else
            exflag = 3
        end if
        Ex(:, :) = 0
        ! figure out if the HF det gets store in the excitation list too?
        ! if yes I have to modify that all a bit, and maybe also in
        ! other parts of the NECI code ... todo
        if (tGUGA) call stop_all(this_routine, "deprecated option with GUGA!")

        do while (.true.)
            call GenExcitations3(HFDet, iLutHF, nJ, exflag, Ex, tParity, tAllExcitsFound, .false.)
            if (tAllExcitsFound) exit !All excits found
            call EncodeBitDet(nJ, iLutnJ)
            if (tHPHF) then
                !Working in HPHF Space. Check whether determinant generated is an 'HPHF'
                if (.not. IsAllowedHPHF(iLutnJ)) cycle
            end if

            iNode = DetermineDetNode(nel, nJ, 0)
            if (iProcIndex == iNode) then
                call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
                amp = amp * PartFac

                if (tRealCoeffByExcitLevel) ExcitLevel = FindBitExcitLevel(iLutnJ, iLutRef(:, 1), nEl)
                if (tAllRealCoeff .or. &
                    & (tRealCoeffByExcitLevel .and. (ExcitLevel <= RealCoeffExcitThresh))) then
                    NoWalkers = amp
                else
                    NoWalkers = int(amp)
                    rat = amp - real(NoWalkers, dp)
                    r = genrand_real2_dSFMT()
                    if (abs(rat) > r) then
                        if (amp < 0.0_dp) then
                            NoWalkers = NoWalkers - 1
                        else
                            NoWalkers = NoWalkers + 1
                        end if
                    end if
                end if

                if (abs(NoWalkers) > 1.0e-12_dp) then
                    call encode_det(CurrentDets(:, DetIndex), iLutnJ)
                    call clear_all_flags(CurrentDets(:, DetIndex))
                    do run = 1, inum_runs
                        temp_sign(run) = NoWalkers
                    end do
                    call encode_sign(CurrentDets(:, DetIndex), temp_sign)

                    ! Store the diagonal matrix elements
                    HDiagtemp = get_diagonal_matel(nJ, iLutnJ)
                    HOffDiagtemp = get_off_diagonal_matel(nJ, iLutnJ)
                    call set_det_diagH(DetIndex, real(HDiagtemp, dp) - Hii)
                    call set_det_offdiagH(DetIndex, HOffDiagtemp)

                    if (associated(lookup_supergroup_indexer)) then
                        call set_supergroup_idx(DetIndex, lookup_supergroup_indexer%idx_nI(nJ))
                    end if

                    ! store the determinant
                    call store_decoding(DetIndex, nJ)

                    if (tTruncInitiator) then
                        !Set initiator flag if needed (always for HF)
                        call CalcParentFlag(DetIndex, iInit)
                    end if

                    DetHash = FindWalkerHash(nJ, nWalkerHashes)
                    TempNode => HashIndex(DetHash)
                    ! If the first element in the list has not been used.
                    if (TempNode%Ind == 0) then
                        TempNode%Ind = DetIndex
                    else
                        do while (associated(TempNode%Next))
                            TempNode => TempNode%Next
                        end do
                        allocate(TempNode%Next)
                        nullify (TempNode%Next%Next)
                        TempNode%Next%Ind = DetIndex
                    end if
                    nullify (TempNode)

                    DetIndex = DetIndex + 1
                    do part_type = 1, lenof_sign
                        TotParts(part_type) = TotParts(part_type) + abs(NoWalkers)
                    end do
                end if
            end if   !End if desired node

        end do

        !Now for the walkers on the HF det
        if (iRefProc(1) == iProcIndex) then
            ! dont have to change this below, should also work for the GUGAc
            if (tAllRealCoeff .or. tRealCoeffByExcitLevel) then
                NoWalkers = PartFac
            else
                NoWalkers = int(PartFac)
                rat = PartFac - real(NoWalkers, dp)
                if (rat < 0.0_dp) &
                    call stop_all(this_routine, "Should not have negative weight on HF")
                r = genrand_real2_dSFMT()
                if (abs(rat) > r) NoWalkers = NoWalkers + 1
            end if
            if (abs(NoWalkers) > 1.0e-12_dp) then
                call encode_det(CurrentDets(:, DetIndex), iLutHF)
                call clear_all_flags(CurrentDets(:, DetIndex))
                do run = 1, inum_runs
                    temp_sign(run) = NoWalkers
                end do
                call encode_sign(CurrentDets(:, DetIndex), temp_sign)
                if (tTruncInitiator) then
                    !Set initiator flag (always for HF)
                    do run = 1, inum_runs
                        call set_flag(CurrentDets(:, DetIndex), get_initiator_flag(run))
                    end do
                end if
                call set_det_diagH(DetIndex, 0.0_dp)
                call set_det_offdiagH(DetIndex, h_cast(0_dp))

                ! store the determinant
                call store_decoding(DetIndex, HFDet)

                ! Now add the Hartree-Fock determinant (not with index 1).
                DetHash = FindWalkerHash(HFDet, nWalkerHashes)
                TempNode => HashIndex(DetHash)
                ! If the first element in the list has not been used.
                if (TempNode%Ind == 0) then
                    TempNode%Ind = DetIndex
                else
                    do while (associated(TempNode%Next))
                        TempNode => TempNode%Next
                    end do
                    allocate(TempNode%Next)
                    nullify (TempNode%Next%Next)
                    TempNode%Next%Ind = DetIndex
                end if
                nullify (TempNode)

                DetIndex = DetIndex + 1
                do run = 1, inum_runs
                    TotParts(run) = TotParts(run) + abs(NoWalkers)
                    NoatHF(run) = NoWalkers
                end do
            else
                call stop_all(this_routine, "No walkers initialised on the HF det with StartMP1")
            end if
        else
            NoatHF(:) = 0.0_dp
        end if

        TotWalkers = DetIndex - 1   !This is the number of occupied determinants on each node
        TotWalkersOld = TotWalkers

        !Set local&global variables
        TotPartsOld = TotParts
        call mpisumall(TotParts, AllTotParts)
        call mpisumall(NoatHF, AllNoatHF)
        call mpisumall(TotWalkers, AllTotWalkers)
        OldAllNoatHF = AllNoatHF
        do run = 1, inum_runs
            OldAllHFCyc(run) = AllNoatHF(run)
            OldAllAvWalkersCyc(run) = AllTotParts(run)
        end do
        AllTotWalkersOld = AllTotWalkers
        AllTotPartsOld = AllTotParts
        iter_data_fciqmc%tot_parts_old = AllTotPartsOld
        AllNoAbortedOld = 0.0_dp

    end subroutine InitFCIMC_MP1

    SUBROUTINE CheckforBrillouins()
        INTEGER :: i, j
        LOGICAL :: tSpinPair

!Standard cases.
        IF ((tHub .and. tReal) .or. (tRotatedOrbs) .or. ((LMS /= 0) .and. (.not. tUHF)) .or. tReltvy) THEN
!Open shell, restricted.
            tNoBrillouin = .true.
        ELSE
!Closed shell restricted, or open shell unrestricted are o.k.
            tNoBrillouin = .false.
            tUseBrillouin = .true.
        end if

!Special case of complex orbitals.
        IF (tFixLz .and. (.not. tNoBrillouin)) THEN
            write(stdout, *) "Turning Brillouins theorem off since we are using non-canonical complex orbitals"
            tNoBrillouin = .true.
        end if

        ! Special case of defining a det with LMS=0, but which is open shell.
        ! No Brillouins if it's a restricted HF calc.
        tSpinPair = .false.
        IF (tDefineDet .and. (LMS == 0) .and. (.not. tUHF)) THEN
            ! If we are defining our own reference determinant, we want to
            ! find out if it is open shell or closed to know whether or not
            ! brillouins theorem holds.
            !
            ! If LMS/=0, then it is easy and must be open shell, otherwise
            ! we need to consider the occupied orbitals.
            do i = 1, (NEl - 1), 2
                ! Assuming things will probably go alpha beta alpha beta,
                ! run through each alpha and see if there's a corresponding
                ! beta.
                tSpinPair = .false.
                IF (MOD(BRR(FDet(i)), 2) /= 0) THEN
!Odd energy, alpha orbital.
                    IF (BRR(FDet(i + 1)) /= (BRR(FDet(i)) + 1)) THEN
                        ! Check the next orbital to see if it's the beta (will
                        ! be alpha+1 when ordered by energy). If not, check
                        ! the other orbitals for the beta, as it's possible
                        ! the orbitals are ordered weird (?).
                        do j = 1, NEl
                            IF (BRR(FDet(j)) == (BRR(FDet(i)) + 1)) tSpinPair = .true.
                        end do
                    ELSE
                        tSpinPair = .true.
                    end if
                ELSE
!Even energy, beta orbital. The corresponding alpha will be beta-1.
                    IF (BRR(FDet(i + 1)) /= (BRR(FDet(i)) - 1)) THEN
                        do j = 1, NEl
                            IF (BRR(FDet(j)) == (BRR(FDet(i)) - 1)) tSpinPair = .true.
                        end do
                    ELSE
                        tSpinPair = .true.
                    end if
                end if
                IF (.not. tSpinPair) EXIT
            end do
            IF (.not. tSpinPair) THEN
!Open shell LMS=0 determinant.
!If restricted HF orbitals are being used, brillouins theorem does not hold.
                tNoBrillouin = .true.
                tUseBrillouin = .false.
                write(stdout, '(A)') " Using an open shell reference determinant in a basis of restricted HF orbitals; " &
                & //"Brillouins theorem is being turned off. "
            end if
        end if

    ENDSUBROUTINE CheckforBrillouins

    SUBROUTINE CalcApproxpDoubles()
        implicit none
        real(dp) :: denom
        INTEGER :: iTotal
        integer :: nSingles, nDoubles, nSing_spindiff1, nDoub_spindiff1, nDoub_spindiff2
        integer :: nTot
        character(*), parameter :: this_routine = "CalcApproxpDoubles"
        integer(n_int), allocatable :: dummy_list(:, :)

        nSingles = 0
        nDoubles = 0
        if (tReltvy) then
            nSing_spindiff1 = 0
            nDoub_spindiff1 = 0
            nDoub_spindiff2 = 0
            pDoub_spindiff1 = 0.0_dp
            pDoub_spindiff2 = 0.0_dp
        end if

!NSing=Number singles from HF, nDoub=No Doubles from HF

        write(stdout, "(A)") " Calculating approximate pDoubles for use with &
                       &excitation generator by looking a excitations from &
                       &reference."
        exflag = 3
        if (tReltvy) then
            write(stdout, *) "Counting magnetic excitations"
            ! subroutine CountExcitations4(nI, minRank, maxRank, minSpinDiff, maxSpinDiff, tot)
            call CountExcitations4(hfdet, 1, 1, 0, 0, nSingles)
            call CountExcitations4(hfdet, 1, 1, 1, 1, nSing_spindiff1)
            call CountExcitations4(hfdet, 2, 2, 0, 0, nDoubles)
            call CountExcitations4(hfdet, 2, 2, 1, 1, nDoub_spindiff1)
            call CountExcitations4(hfdet, 2, 2, 2, 2, nDoub_spindiff2)
            call CountExcitations4(hfdet, 1, 2, 0, 2, nTot)
            ASSERT(nTot == (nSingles + nSing_spindiff1 + nDoubles + nDoub_spindiff1 + nDoub_spindiff2))

            iTotal = nSingles + nDoubles + nSing_spindiff1 + nDoub_spindiff1 + nDoub_spindiff2

        else
            if (tKPntSym) THEN
                if (t_k_space_hubbard) then
                    ! change this to the new implementation
                    call gen_all_excits_k_space_hubbard(HFDet, nDoubles, dummy_list)
                else
                    call enumerate_sing_doub_kpnt(exFlag, .false., nSingles, nDoubles, .false.)
                end if
            else
                call CountExcitations3(hfdet, exflag, nSingles, nDoubles)
            end if
            iTotal = nSingles + nDoubles
        end if

        IF (tHub .or. tUEG) THEN
            IF (tReal) THEN
                write(stdout, *) "Since we are using a real-space hubbard model, only single excitations are connected &
                &   and will be generated."
                pDoubles = 0.0_dp
                if (tReltvy) then
                    pDoub_spindiff1 = 0.0_dp
                    pDoub_spindiff2 = 0.0_dp
                    pSingles = real(nSingles, dp) / real(nSingles + nSing_spindiff1, dp)
                    pSing_spindiff1 = 1.0_dp - pSingles
                else
                    pSingles = 1.0_dp
                end if
                return
            ELSE
                write(stdout, *) "Since we are using a momentum-space hubbard model/UEG, only double excitaitons &
     &                          are connected and will be generated."
                pSingles = 0.0_dp
                if (tReltvy) then
                    pSing_spindiff1 = 0.0_dp
                    pDoubles = real(nDoubles, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
                    pDoub_spindiff1 = real(nDoub_spindiff1, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
                    pDoub_spindiff2 = 1.0_dp - pDoubles - pDoub_spindiff1
                else
                    pDoubles = 1.0_dp
                end if
                return
            end if

        else if (tNoSingExcits) then
            pSingles = 0.0_dp
            if (tReltvy) then
                pSing_spindiff1 = 0.0_dp
                pDoubles = real(nDoubles, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
                pDoub_spindiff1 = real(nDoub_spindiff1, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
                pDoub_spindiff2 = 1.0_dp - pDoubles - pDoub_spindiff1
            else
                pDoubles = 1.0_dp
            end if
            write(stdout, *) "Only double excitations will be generated"
            return
        end if

        write(stdout, "(I7,A,I7,A)") nDoubles, " double excitations, and ", nSingles, &
            " single excitations found from reference. This will be used to calculate pDoubles."

        IF (abs(SinglesBias - 1.0_dp) > 1.0e-12_dp) THEN
            write(stdout, *) "Singles Bias detected. Multiplying single excitation connectivity of HF determinant by ", &
                SinglesBias, " to determine pDoubles."
        end if

        IF ((nSingles == 0) .or. (nDoubles == 0)) THEN
            write(stdout, *) "Number of singles or doubles found equals zero. pDoubles will be set to 0.95. Is this correct?"
            pDoubles = 0.95_dp
            pSingles = 0.05_dp
            return
        else if ((nSingles < 0) .or. (nDoubles < 0)) then
            call stop_all("CalcApproxpDoubles", &
                          "Number of singles, doubles or Yamanouchi symbols &
                          &found to be a negative number. Error here.")
        end if
        if (tReltvy) then
            denom = real(nSingles + nSing_spindiff1, dp) * SinglesBias &
                    + real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
            pSingles = real(nSingles, dp) * SinglesBias / denom
            pSing_spindiff1 = real(nSing_spindiff1, dp) * SinglesBias / denom
            pDoubles = real(nDoubles, dp) / denom
            pDoub_spindiff1 = real(nDoub_spindiff1, dp) / denom
            pDoub_spindiff2 = 1.0_dp - pSingles - pSing_spindiff1 &
                              - pDoubles - pDoub_spindiff1 - pDoub_spindiff2
        else
            denom = real(nSingles, dp) * SinglesBias + real(nDoubles, dp)
            pSingles = real(nSingles, dp) * SinglesBias / denom
            pDoubles = 1.0_dp - pSingles
        end if

        IF (abs(SinglesBias - 1.0_dp) > 1.0e-12_dp) THEN
            write (stdout, '("pDoubles set to ", f14.6, &
                       &" rather than (without bias): ", f14.6)') &
                       pDoubles, real(nDoubles, dp) / real(iTotal, dp)
            write (stdout, '("pSingles set to ", f14.6, &
                       &" rather than (without bias): ", f14.6)') &
                       pSingles, real(nSingles, dp) / real(iTotal, dp)

!            write(stdout,"(A,F14.6,A,F14.6)") "pDoubles set to: ",pDoubles, " rather than (without bias): ", &
!                & real(nDoub,dp)/real(iTotal,dp)
        ELSE
            if (tReltvy) then
                write(stdout, '(A)') " Where s and t are alpha or beta spin function labels: "
                write(stdout, '(A30,F14.6)') " pSingles(s->s) set to: ", pSingles
                write(stdout, '(A30,F14.6)') " pSingles(s->s') set to: ", pSing_spindiff1
                write(stdout, '(A30,F14.6)') " pDoubles(st->st) set to: ", pDoubles
                write(stdout, '(A30,F14.6)') " pDoubles(st->s't) set to: ", pDoub_spindiff1
                write(stdout, '(A30,F14.6)') " pDoubles(st->s't') set to: ", pDoub_spindiff2
            else
                write(stdout, '(A,F14.6)') " pDoubles set to: ", pDoubles
                write(stdout, '(A,F14.6)') " pSingles set to: ", pSingles
            end if
        end if

        if (allocated(pSinglesIn) .or. allocated(pDoublesIn) .or. allocated(pTriplesIn)) then
            if (allocated(pSinglesIn) .and. allocated(pDoublesIn)) then
                call stop_all(this_routine, 'It is not possible to define pSingles and pDoubles.')
            else if (.not. (allocated(pSinglesIn) .or. allocated(pDoublesIn))) then
                call stop_all(this_routine, 'One of pSingles or pDoubles is required.')
            end if
            if (t_mol_3_body) then
                ! We allow the users to input absolute values for the probabilities.
                ! Note that pSingles and pDoubles are internally conditional probabilities.
                ! Even for triple we still have that `pSingles + pDoubles .isclose. 1.0`.
                ! It first decides upon triple excitation or something else and then about singles or doubles.
                ! We have to convert the absolute probabilities into conditional ones.
                if (.not. allocated(pTriplesIn)) then
                    call stop_all(this_routine, "pTriples is required as input.")
                else
                    pTriples = pTriplesIn
                    if (allocated(pSinglesIn)) then
                        if (pTriples + pSinglesIn > 1.0_dp) call stop_all(this_routine, "pTriplesIn + pSinglesIn > 1.0_dp")
                        pSingles = pSinglesIn / (1.0_dp - pTriplesIn)
                        pDoubles = 1.0_dp - pSingles
                        write(stdout, '(" Using the input value of pSingles:",1x, f14.6)') pSinglesIn
                    else if (allocated(pDoublesIn)) then
                        if (pTriples + pDoublesIn > 1.0_dp) call stop_all(this_routine, "pTriplesIn + pDoublesIn > 1.0_dp")
                        pDoubles = pDoublesIn / (1.0_dp - pTriplesIn)
                        pSingles = 1.0_dp - pDoubles
                        write(stdout, '(" Using the input value of pDoubles:",1x, f14.6)') pDoublesIn
                    end if
                end if
            else
                if (allocated(pTriplesIn)) then
                    call stop_all(this_routine, "pTriples can only be given if triple excitations are performed.")
                else if (allocated(pSinglesIn)) then
                        pSingles = pSinglesIn
                        pDoubles = 1.0_dp - pSingles
                        write(stdout, '(" Using the input value of pSingles:",1x, f14.6)') pSingles
                else if (allocated(pDoublesIn)) then
                    pDoubles = pDoublesIn
                    pSingles = 1.0_dp - pDoubles
                    write(stdout, '(" Using the input value of pDoubles:",1x, f14.6)') pDoubles
                end if
            end if
        end if

        if (allocated(pParallelIn)) then
            write(stdout, '(" Using the input value of pParallel:",1x, f14.6)') pParallelIn
            pParallel = pParallelIn
        end if

    END SUBROUTINE CalcApproxpDoubles

    SUBROUTINE CreateSpinInvBRR()

        ! Create an SpinInvBRR containing spin orbitals,
        ! unlike 'createInvBRR' which only has spatial orbitals.
        ! This is used for the FixCASshift option in establishing whether or not
        ! a determinant is in the complete active space.
        ! In:
        !    BRR(i)=j: orbital i is the j-th lowest in energy.
        !    nBasis: size of basis
        ! SpinInvBRR is the inverse of BRR.  SpinInvBRR(j)=i: the j-th lowest energy
        ! orbital corresponds to the i-th orbital in the original basis.
        ! i.e the position in SpinInvBRR now corresponds to the orbital number and
        ! the value to the relative energy of this orbital.

        IMPLICIT NONE
        INTEGER :: I, t, ierr
        CHARACTER(len=*), PARAMETER :: this_routine = 'CreateSpinInvBrr'

        IF (ALLOCATED(SpinInvBRR)) return

        allocate(SpinInvBRR(NBASIS), STAT=ierr)
        CALL LogMemAlloc('SpinInvBRR', NBASIS, 4, this_routine, SpinInvBRRTag, ierr)

        SpinInvBRR(:) = 0

        t = 0
        do I = 1, NBASIS
            t = t + 1
            SpinInvBRR(BRR(I)) = t
        end do

        return

    END SUBROUTINE CreateSpinInvBRR

    subroutine SetupValidSpawned(WalkerListSize)

        use CalcData, only: MemoryFacSpawn

        implicit none
        integer(int64), intent(in) :: WalkerListSize
        integer ierr, j
        real(dp) Gap

        !When running normally, WalkerListSize will be equal to totalwalkers
        !However, when reading in (and not continuing to grow) it should be equal to the number of dets in the popsfile
        MaxSpawned = NINT(MemoryFacSpawn * WalkerListSize * inum_runs)
!            write(stdout,"(A,I14)") "Memory allocated for a maximum particle number per node for spawning of: ",MaxSpawned

!      write(stdout,"(A)") "*Direct Annihilation* in use...Explicit load-balancing disabled."
        allocate(ValidSpawnedList(0:nNodes - 1), stat=ierr)
        ! InitialSpawnedSlots is now filled later, once the number of particles
        ! wanted is known
        !(it can change according to the POPSFILE).
        allocate(InitialSpawnedSlots(0:nNodes - 1), stat=ierr)
        ! InitialSpawnedSlots now holds the first free position in the
        ! newly-spawned list for each processor, so it does not need to be
        ! reevaluated each iteration.
!      MaxSpawned=NINT(MemoryFacSpawn*InitWalkers)
        Gap = REAL(MaxSpawned, dp) / REAL(nNodes, dp)
        do j = 0, nNodes - 1
            InitialSpawnedSlots(j) = NINT(Gap * j) + 1
        end do
        ! ValidSpawndList now holds the next free position in the newly-spawned
        ! list, but for each processor.
        ValidSpawnedList(:) = InitialSpawnedSlots(:)

    end subroutine SetupValidSpawned

    subroutine sync_rdm_sampling_iter()
        use LoggingData, only: RDMEnergyIter, IterRDMOnfly
        use CalcData, only: coreSpaceUpdateCycle, semistoch_shift_iter
        implicit none
        integer :: frac
        ! first, adjust the offset to make the rdm sampling start right when a semistochastic
        ! update cycle ends
        IterRDMOnFly = IterRDMOnFly - &
                       mod(IterRDMOnFly - semistoch_shift_iter, coreSpaceUpdateCycle) - 1
        ! The -1 is just because the sampling starts one iteration after IterRDMOnFly
        ! If we subtracted too much, jump one cycle backwards
        if (IterRDMOnFly < semistoch_shift_iter) IterRDMOnFly = IterRDMOnFly + coreSpaceUpdateCycle
        write(stdout, *) "Adjusted starting iteration of RDM sampling to ", IterRDMOnFly

        ! Now sync the update cycles
        if (RDMEnergyIter > coreSpaceUpdateCycle) then
            RDMEnergyIter = coreSpaceUpdateCycle
            write(stdout, *) "The RDM sampling interval cannot be larger than the update " &
                //"interval of the semi-stochastic space. Reducing it to ", RDMEnergyIter
        end if
        if (mod(coreSpaceUpdateCycle, RDMEnergyIter) /= 0) then
            ! first, try to ramp up the RDMEnergyIter to meet the coreSpaceUpdateCycle
            frac = coreSpaceUpdateCycle / RDMEnergyIter
            RDMEnergyIter = coreSpaceUpdateCycle / frac
            write(stdout, *) "Update cycle of semi-stochastic space and RDM sampling interval" &
                //" out of sync. "
            write(stdout, *) "Readjusting RDM sampling interval to ", RDMEnergyIter

            ! now, if this did not succeed, adjust the coreSpaceUpdateCycle
            if (mod(coreSpaceUpdateCycle, RDMEnergyIter) /= 0) then
                coreSpaceUpdateCycle = coreSpaceUpdateCycle - &
                                       mod(coreSpaceUpdateCycle, RDMEnergyIter)
                write(stdout, *) "Adjusted update cycle of semi-stochastic space to ", &
                    coreSpaceUpdateCycle
            end if
        end if
    end subroutine sync_rdm_sampling_iter

    subroutine CalcUEGMP2()
        use SymExcitDataMod, only: kPointToBasisFn
        use SystemData, only: ElecPairs, NMAXX, NMAXY, NMAXZ, OrbECutOff, &
                              tMP2UEGRestrict, kiRestrict, kiMsRestrict, kjRestrict, kjMsRestrict, &
                              Madelung, tMadelung, tUEGFreeze, FreezeCutoff, kvec, tUEG2
        use Determinants, only: GetH0Element4, get_helement_excit
        integer :: Ki(3), Kj(3), Ka(3), LowLoop, HighLoop, X, i, Elec1Ind, Elec2Ind, K, Orbi, Orbj
        integer :: iSpn, FirstA, nJ(NEl), a_loc, Ex(2, maxExcit), kx, ky, kz, OrbB
        integer :: ki2, kj2
        logical :: tParity
        real(dp) :: Ranger, length
        HElement_t(dp) :: hel, H0tmp, mp2, mp2all

        !Divvy up the ij pairs
        Ranger = real(ElecPairs, dp) / real(nProcessors, dp)
        LowLoop = int(iProcIndex * Ranger) + 1
        Highloop = int((iProcIndex + 1) * Ranger)

        if ((iProcIndex + 1) == nProcessors) Highloop = ElecPairs
        if (iProcIndex == 0) then
            if (lowLoop /= 1) write(stderr, *) "Error here!"
        end if
        write(stdout, *) "Total ij pairs: ", ElecPairs
        write(stdout, *) "Considering ij pairs from: ", LowLoop, " to ", HighLoop
!        write(stdout,*) "HFDet: ",HFDet(:)

        do i = LowLoop, HighLoop   !Looping over electron pairs on this processor

            X = ElecPairs - i
            K = INT((SQRT(8.0_dp * REAL(X, dp) + 1.0_dp) - 1.0_dp) / 2.0_dp)
            Elec1Ind = NEl - 1 - K
            Elec2Ind = NEl - X + ((K * (K + 1)) / 2)
            Orbi = HFDet(Elec1Ind)
            Orbj = HFDet(Elec2Ind)
            Ki = G1(Orbi)%k
            Kj = G1(Orbj)%k
            !=======================================
            if (tUEG2) then
                Ki = kvec(Orbi, 1:3)
                Kj = kvec(Orbj, 1:3)
            end if
            !=======================================
            if (tUEGFreeze) then
                ki2 = ki(1)**2 + ki(2)**2