real_time_init.F90 Source File


Source Code

Source Code

#include "macros.h"

! this is the initialization module for the real-time FCIQMC implementation

module real_time_init

    use real_time_data, only: t_real_time_fciqmc, gf_type, real_time_info, &
                              gf_overlap, temp_det_list, dyn_norm_red, &
                              temp_det_pointer, temp_det_hash, temp_freeslot, tOverpopulate, &
                              pert_norm, second_spawn_iter_data, DiagParts, stepsAlpha, &
                              DiagVec, normsize, valid_diag_spawns, tStabilizerShift, &
                              tDynamicAlpha, &
                              gf_count, &
                              allPopSnapshot, &
                              current_overlap, &
                              TotPartsStorage, t_rotated_time, TotPartsPeak, asymptoticShift, &
                              tau_imag, tau_real, elapsedRealTime, elapsedImagTime, tNewOverlap, &
                              TotWalkers_orig, dyn_norm_psi, gs_energy, shift_damping, &
                              tStaticShift, MaxSpawnedDiag, tDynamicCoreSpace, overlap_states, &
                              overlap_real, overlap_imag, allGfs, tRealTimePopsfile, &
                              tLimitShift, nspawnMax, shiftLimit, numCycShiftExcess, &
                              TotPartsLastAlpha, alphaDamping, tDynamicDamping, iterInit, &
                              etaDamping, tStartVariation, rotThresh, stabilizerThresh, &
                              tInfInit, popSnapshot, snapshotOrbs, phase_factors, tVerletSweep, &
                              numSnapshotOrbs, tLowerThreshold, t_kspace_operators, tVerletScheme, &
                              tLogTrajectory, tReadTrajectory, alphaCache, tauCache, trajFile, &
                              tGenerateCoreSpace, tGZero, wn_threshold, corespace_log_interval, &
                              alphaLog, alphaLogSize, alphaLogPos, tStaticShift, DiagVecTag, &
                              tOnlyPositiveShift, tHFOverlap
    use real_time_procs, only: create_perturbed_ground, setup_temp_det_list, &
                               calc_norm, clean_overlap_states, openTauContourFile
    use verlet_aux, only: backup_initial_state, setup_delta_psi
    use constants, only: dp, n_int, int64, lenof_sign, inum_runs, stdout
    use Parallel_neci
    use MPI_wrapper, only: iProcIndex, root, MPIbarrier, nNodes, MPI_SUM
    use util_mod, only: get_unique_filename, get_free_unit
    use Logging, only: tIncrementPops
    use MemoryManager, only: LogMemAlloc, LogMemDealloc
    use kp_fciqmc_data_mod, only: tMultiplePopStart, tScalePopulation, &
                                  tOverlapPert, overlap_pert, scaling_factor
    use CalcData, only: tChangeProjEDet, tReadPops, tRestartHighPop, tFCIMC, tJumpShift, &
                        tStartSinglePart, nmcyc, iPopsFileNoRead, tWritePopsNorm, &
                        tWalkContGrow, diagSft, pops_norm, InitWalkers, MemoryFacSpawn, &
                        StepsSft, tSemiStochastic, tTruncInitiator, tAddToInitiator
    use FciMCData, only: alloc_popsfile_dets, pops_pert, tPopsAlreadyRead, &
                         tSinglePartPhase, iter_data_fciqmc, iter, PreviousCycles, &
                         AllGrowRate, spawn_ht, pDoubles, pSingles, TotParts, &
                         MaxSpawned, TotWalkers, SumWalkersCyc, &
                         CurrentDets, popsfile_dets, MaxWalkersPart, WalkVecDets, &
                         SpawnedParts, core_run
    use core_space_util, only: cs_replicas
    use SystemData, only: lms, G1, nBasisMax, tHub, nel, tComplexWalkers_RealInts, &
         nBasis, tReal, t_k_space_hubbard, t_complex_ints
    use k_space_hubbard, only: init_k_space_hubbard
    use SymExcitDataMod, only: kTotal
    use sym_mod, only: MomPbcSym
    use perturbations, only: init_perturbation_annihilation, &
    use fcimc_initialisation, only: SetupParameters, InitFCIMCCalcPar, &
    use LoggingData, only: tZeroProjE, tFCIMCStats2
    use fcimc_output, only: write_fcimcstats2, WriteFciMCStatsHeader
    use replica_data, only: allocate_iter_data, set_initial_global_data
    use bit_rep_data, only: IlutBits, niftot, extract_sign, nifd
    use bit_reps, only: decode_bit_det
    use adi_references, only: setup_reference_space
    use tau_main, only: tau, assign_value_to_tau

    implicit none

    real(dp) :: benchmarkEnergy


    subroutine init_real_time_calc_single()
        use real_time_procs, only: get_tot_parts
        ! this routine takes care of the correct setup of the real-time
        ! calculation. like reading the popsfiles and preparing the start
        ! of the calculation and setting certain global variables
        implicit none
        character(*), parameter :: this_routine = "init_real_time_calc_single"

        integer :: ierr

        write(stdout, *) " Entering real-time FCIQMC initialisation "

        ! think about what variables have to be set for a succesful calc.

        ! also call the "normal" NECI setup routines to allow calculation
        call SetupParameters()

        ! for real-space hubbard, we conver to momentum operators if desired
        if (t_kspace_operators) call setup_momentum_operators()

        ! have to think about the the order about the above setup routines!
        ! within this Init a readpops is called
        ! this function already produces the correctly perturbed ground state
        call InitFCIMCCalcPar()

        ! also init pointer here, and think about what options and defaults
        ! should be set for a succsesfull init
        call init_fcimc_fn_pointers()

        ! setup the k-space hubbard if required (after pointers as some are
        ! overwritten
        if (t_k_space_hubbard) then
            call init_k_space_hubbard()
        end if
        ! then call the setup routine, which set all remaining needed quantities
        call setup_real_time_fciqmc()

        ! definetly read-in stored popsfile here.
        ! need to store both <y(0)| and also create a_j y(0)> during read-in!
!         call read_popsfile_real_time()
        ! actually the InitFCIMCCalcPar should do that now correctly already

        ! do an MPIbarrier here.. although don't quite know why
        call MPIBarrier(ierr)

        if (.not. tReadPops) call set_initial_global_data(TotWalkers, CurrentDets)

    end subroutine init_real_time_calc_single

    subroutine setup_real_time_fciqmc()
        ! this is the last setup routine, which depending on compilation,
        ! number of copies etc. sets up the final needed quantities to run
        ! a simulation
        implicit none
        character(*), parameter :: this_routine = "setup_real_time_fciqmc"
        integer :: ierr, run

        ! the new total momentum has to be constructed before the
        ! time-evolved state is read in, as the latter deletes the
        ! pops_pert, because perturbation and read-in are done in one
        ! function (dependencies...)
        if (tHub) then
            if (allocated(pops_pert)) then

                if (pops_pert(1)%nannihilate == 1) kTotal = kTotal &
                                                            - G1(pops_pert(1)%ann_orbs(1))%k
                if (pops_pert(1)%ncreate == 1) kTotal = kTotal &
                                                        + G1(pops_pert(1)%crtn_orbs(1))%k
                call MomPbcSym(kTotal, nBasisMax)
                write(stdout, *) "New total momentum", kTotal
            end if
        end if

        ! allocate the according quantities!
        ! n_time_steps have to be set here!
        write(stdout, *) " Allocating greensfunction and wavefunction norm arrays!"
        ! allocate an additional slot for initial values
        if (numSnapshotOrbs > 0) then
            allocate(popSnapshot(numSnapshotOrbs), stat=ierr)
            allocate(allPopSnapshot(numSnapshotOrbs), stat=ierr)
            popSnapshot = 0.0_dp
            allPopSnapshot = 0.0_dp
            allocate(popSnapshot(1), stat=ierr)
            allocate(allPopSnapshot(numSnapshotOrbs), stat=ierr)
            popSnapshot = 0.0_dp
            allPopSnapshot = 0.0_dp
        end if
        allocate(gs_energy(inum_runs), stat=ierr)
        allocate(temp_freeslot(MaxWalkersPart), stat=ierr)
        allocate(TotPartsPeak(inum_runs), stat=ierr)
        allocate(numCycShiftExcess(inum_runs), stat=ierr)
        ! allocate the buffer for storing previous values of alpha
        ! for now, take 50 values of alpha in the log
        alphaLogSize = 50
        alphaLogPos = 1
        allocate(alphaLog(alphaLogSize), stat=ierr)
        alphalog = 0.0_dp
        numCycShiftExcess = 0
        ! allocate spawn buffer for verlet scheme
        if (tVerletScheme) allocate(spawnBuf(0:niftot, 1:maxSpawned))

        TotPartsPeak = 0.0_dp
        gs_energy = benchmarkEnergy

        ! when projecting onto the perturbed reference, we obviously need to create
        ! a new state
        if (tHFOverlap) tNewOverlap = .true.

        call init_overlap_buffers()

        if (tRealTimePopsfile) call readTimeEvolvedState()

        ! check for set lms.. i think that does not quite work yet
        write(stdout, *) "mz spin projection: ", lms

        write(stdout, *) "tSinglePartPhase?:", tSinglePartPhase
        write(stdout, *) "tWalkContGrow?", tWalkContGrow
        write(stdout, *) "diagSft:", diagSft

        ! intialize the 2nd temporary determinant list needed in the
        ! real-time fciqmc

        ! also maybe use the spawn_ht hash table, so allocated it here!
        call setup_temp_det_list()

        write(stdout, *) "allocated(temp_det_list)?", allocated(temp_det_list)
        write(stdout, *) "associated(temp_det_pointer)?", associated(temp_det_pointer)
        write(stdout, *) "associated(temp_det_hash)?", associated(temp_det_hash)

        write(stdout, *) "associated(spawn_ht)?", associated(spawn_ht)

        write(stdout, *) "Allgrowrate: ", AllGrowRate
        ! print out the first infos on the calculation..
        ! although that definetly has to be changed for the real-time fciqm

        ! use new output format!
        tFCIMCStats2 = .true.

        if (tFCIMCStats2) then
            call write_fcimcstats2(iter_data_fciqmc, initial=.true.)
            call WriteFciMCStatsHeader()
        end if

        ! set the iter variable to 0 probably
        iter = 0

        ! and also the PreviousCycles var. since its essentially regarded as
        ! a new calulcation
        PreviousCycles = 0

        ! for intermediate test_purposes turn off spawning to check if the
        ! diagonal step works as intented
!         pSingles = 0.0_dp
!         pDoubles = 0.0_dp

        ! also initialize the second_spawn_iter_data type
        call allocate_iter_data(second_spawn_iter_data)

        ! and also initialize the values:
        second_spawn_iter_data%ndied = 0.0_dp
        second_spawn_iter_data%nborn = 0.0_dp
        second_spawn_iter_data%nannihil = 0.0_dp
        second_spawn_iter_data%naborted = 0.0_dp
        second_spawn_iter_data%nremoved = 0.0_dp
        second_spawn_iter_data%update_growth = 0.0_dp
        second_spawn_iter_data%update_growth_tot = 0.0_dp
        second_spawn_iter_data%tot_parts_old = TotParts
        second_spawn_iter_data%update_iters = 0

        TotPartsStorage = TotParts
        TotPartsLastAlpha = TotParts

        ! also intitialize the 2nd spawning array to deal with the
        ! diagonal death step in the 2nd rt-fciqmc loop
        allocate(DiagVec(0:IlutBits%len_bcast, MaxWalkersPart), stat=ierr)
        call LogMemAlloc('DiagVec', MaxWalkersPart * (1 + IlutBits%len_bcast), size_n_int, &
                         this_routine, DiagVecTag, ierr)

        DiagVec = 0

        DiagParts => DiagVec

        ! and the initial_spawn_slots equivalent
        ! although i think i can reuse the initialSpawnedSlots..
!         allocate(initial_diag_spawn_list(0:nNodes-1), stat = ierr)

        valid_diag_spawns = 1

        do run = 1, inum_runs
            SumWalkersCyc(run) = SumWalkersCyc(run) + &
        end do

        tVerletSweep = .false.
        if (tVerletScheme) then
            call setup_delta_psi()
            call backup_initial_state()
            call assign_value_to_tau(tau / iterInit, this_routine)
        end if

        if (tStaticShift) DiagSft = asymptoticShift

        if (tGenerateCoreSpace) call initialize_corespace_construction()

        if (tReadTrajectory) call read_in_trajectory()

        ! Set up the reference space for the adi-approach
        call setup_reference_space(tReadPops)

        call rotate_time()

    end subroutine setup_real_time_fciqmc

    subroutine rotate_time()
        implicit none
        ! to avoid code multiplication
        if (t_rotated_time) then
            tau_imag = -sin(real_time_info%time_angle) * tau
            tau_real = cos(real_time_info%time_angle) * tau
            tau_imag = 0.0_dp
            tau_real = tau
        end if
    end subroutine rotate_time

    subroutine init_overlap_buffers
        use CalcData, only: tSemiStochastic, ss_space_in
        use semi_stoch_gen, only: init_semi_stochastic
        use real_time_procs, only: reset_tot_parts
        implicit none
        ! this subroutine sets up everything required to compute green's functions
        integer :: ierr, j, i
        complex(dp), allocatable :: norm_buf(:)
        logical :: tStartedFromCoreGround

        normsize = inum_runs**2
        allocate(overlap_real(gf_count), overlap_imag(gf_count))
        allocate(gf_overlap(normsize, gf_count), stat=ierr)
        allocate(pert_norm(normsize, gf_count), stat=ierr)
        allocate(dyn_norm_psi(normsize), stat=ierr)
        dyn_norm_psi = 1.0_dp
        allocate(dyn_norm_red(normsize, gf_count), stat=ierr)
        allocate(current_overlap(normsize, gf_count), stat=ierr)
        dyn_norm_red = 1.0_dp
        gf_overlap = 0.0_dp

        ! also need to create the perturbed ground state to calculate the
        ! overlaps to |y(t)>
        call create_perturbed_ground()

        if (tSemiStochastic) call init_semi_stochastic(ss_space_in, tStartedFromCoreGround)

        ! If only the corespace time-evolution is to be taken, truncate the
        ! initial wavefunction to the corespace
        ! We currently do not truncate the overlap state too, but it might come
        ! later
        if (tGZero) then
            call truncate_initial_state()
            call truncate_overlap_states()
            call reset_tot_parts()
            TotWalkers = cs_replicas(core_run)%determ_sizes(iProcIndex)
        end if

        ! if a ground-state POPSFILE is used, we need to ensure coherence between the replicas
        if (.not. tRealTimePopsfile) call equalize_initial_phase()

        allocate(norm_buf(normsize), stat=ierr)
        ! to avoid dividing by 0 if not all entries get filled
        pert_norm = 1.0_dp
        norm_buf = calc_norm(CurrentDets, int(TotWalkers))
        call MPIReduce(norm_buf, MPI_SUM, dyn_norm_psi)
        do j = 1, gf_count
            ! calc. the norm of the perturbed ground states
            norm_buf = calc_norm(overlap_states(j)%dets, overlap_states(j)%nDets)
            write(stdout, *) "Number of walkers / in overlap state", TotWalkers, overlap_states(j)%nDets
            ! the norm (squared) can be obtained by reduction over all processes
            call MPIReduce(norm_buf, MPI_SUM, pert_norm(:, j))
            ! for diagonal green's functions, this is the same as pert_norm, but
            ! in general, this general normalization is required.
            do i = 1, normsize
                dyn_norm_red(i, j) = sqrt(pert_norm(i, j) * dyn_norm_psi(i))
            end do
        end do


        if (.not. allocated(shift_damping)) then
            allocate(shift_damping(inum_runs), stat=ierr)
            shift_damping = 0.0_dp
        end if

        if (tLogTrajectory) call openTauContourFile()

    end subroutine init_overlap_buffers

    ! write a routine which sets the defauls values, and maybe even turns off
    ! certain otherwise non-compatible variable to the appropriate values
    ! problem: this is called after most of the input is read
    subroutine set_real_time_defaults()
        implicit none

        ! by default, use runge-kutta instead of verlet
        tVerletScheme = .false.
        iterInit = 1

        ! todo: figure out quantities
        ! maximum number of spawn attempts per determinant if the number of
        ! total spawns exceeds some threshold
        nspawnMax = 1000

        ! energy offset
        benchmarkEnergy = 0.0_dp

        ! also set readpops to get the <y(0)| reference from the "normal"
        ! neci init routines
        tReadPops = .true.
        ! By default, apply a seperate perturbation to the overlap_states
        tNewOverlap = .true.

        ! startsinglepart does not work with popsfile and is not wanted too
        tStartSinglePart = .false.
        ! also, the shift is something infamous in real-time, so definitely no
        ! jumpshift
        tJumpShift = .false.

        ! but to ensure that the shift does not vary anymore, since there is
        ! no such concept as the varying shift in the real-time fciqmc
        ! exception: when using rotated times, the shift still has to be considered
        tWalkContGrow = .true.

        ! probably not change reference anymore.. but check
        tChangeProjEDet = .true.

        ! and dont only restart on the highly populated dets
        tRestartHighPop = .false.

        ! nick also has some scale_population routine, think about
        ! implementing this also for the real-time restart
        tScalePopulation = .false.
        scaling_factor = 1.0_dp

        ! nick also has some multiple popsfile start routines..
        ! set the multiple popsstart with the number of replicas of mneci
        ! provided
        tMultiplePopStart = .false.

        ! by default, do not output any population of orbitals
        numSnapshotOrbs = 0

        ! from my way of outputting popsfiles i always do it in popsfile.n
        ! format -> so i probably have to set tIncrementPops and the count
        tIncrementPops = .true.
        iPopsFileNoRead = -1
        ! this will always result in *.0 being chosen as name, was rethought and
        ! decided to be good - this way, no files will be overwritten and
        ! the read-in file is always the same

        ! have to set, that popsfile is not yet read in:
        tPopsAlreadyRead = .false.

        ! overwrite tfcimc flag to not enter the regular fcimc routine
        tFCIMC = .false.

        ! Default to complex ints, this is then turned off when running real-time
        !t_complex_ints = .true.

        ! probably should zero the projected energy, since its a total
        ! different system
        tZeroProjE = .false.

        ! setup_rotated_time: by default, pure real time is used
        t_rotated_time = .false.
        ! usually, the walker number shall be controlled
        tStaticShift = .false.
        asymptoticShift = 0.0_dp

        ! and in case of semistochastic approach, the core space shall be static
        tDynamicCoreSpace = .false.

        ! usually, systems with real integrals will be considered, but the walkers will
        ! always be complex
        tComplexWalkers_RealInts = .true.

        ! by default, the perturbation operators are in the same basis as the wavefunction
        t_kspace_operators = .false.

        ! if no gf kind is specified, only the overlap with the initial state will
        ! be considered -> only one overlap is obtained
        gf_count = 1
        allGfs = 0
        ! normally, take the multi-replica overlap
        tHFOverlap = .false.

        ! if starting a new calculation, we start at time 0 (arbitrary)
        elapsedRealTime = 0.0_dp
        elapsedImagTime = 0.0_dp

        ! by default, the initial state is taken from an ordinary popsfile
        ! if a time evolved state is desired, a second popsfile has to be supplied
        tRealTimePopsfile = .false.
        tStabilizerShift = .false.
        ! threshold population (in relation to the peak walker number) for activation of
        ! stabilizer shift
        stabilizerThresh = 0.8
        ! the merging of spawning events is done entirely automatically and therfore can not
        ! be switched on manually
        tLimitShift = .false.
        shiftLimit = 0.7_dp
        ! by default, negative shifts are not allowed as they can easily lead to instability
        tOnlyPositiveShift = .true.

        ! default values for dynamic rotation angle updating (it is not enabled by default)
        tDynamicAlpha = .false.
        stepsAlpha = 10
        alphaDamping = 0.05
        rotThresh = 0
        tLowerThreshold = .false.
        tOverpopulate = .false.
        tStartVariation = .false.

        ! the damping is constant by default
        tDynamicDamping = .false.
        etaDamping = 0.01

        ! and we allow normal usage of inititators
        tInfInit = .false.

        ! and we do not print out the contour separately
        tLogTrajectory = .false.
        tReadTrajectory = .false.

        tGenerateCoreSpace = .false.
        wn_threshold = 0.01
        corespace_log_interval = 300
        ! Get the full Green's function, not only in the corespace
        tGZero = .false.

        ! Do not use additional damping quadratic in energy
        real_time_info%quad_damp_fac = 0.0_dp
    end subroutine set_real_time_defaults

    ! need a specific popsfile read function for the real-time calculation
    ! based upon georges previous implementation, but changed to handle
    ! new code base
    ! need a routine which prepares the converged groundstates from an
    ! imaginary-time FCIQMC calculation

    ! use code provided in the NECI GreensFuncs branch by George Booth
    ! in general reuse as much of the already provided functionality!

    ! need a setup routine in the regular imag-time neci routine, which prints
    ! out the specified amount of GS wavefunctions
    ! and which sets the necessary flags and stuff

    ! i want through the CHANGEVARS facility start the print out of
    ! the input specified print out of POPSFILES between certain intervals

    subroutine read_popsfile_real_time()
        use PopsfileMod, only: open_pops_head, FindPopsfileVersion, ReadPopsHeadv4
        implicit none

        integer :: iunit, popsversion, iPopLenof_Sign, iPopNel, iPopIter, &
                   PopNIfD, PopNIfSgn, PopNIfFlag, PopNIfTot, &
                   PopBlockingIter, Popinum_runs, PopRandomHash(2056), &
                   read_nnodes, PopBalanceBlocks
        logical :: formpops, binpops, tPop64Bit, tPopHPHF, tPopLz
        integer(int64) :: iPopAllTotWalkers, read_walkers_on_nodes(0:nProcessors - 1)
        real(dp) :: PopDiagSft(inum_runs), read_tau, PopSumNoatHF(lenof_sign), &
                    read_psingles, read_pparallel
        HElement_t(dp) :: PopAllSumENum(inum_runs)
        character(255) :: popsfile
        real(dp) :: unused_triples

        character(*), parameter :: this_routine = "read_popsfile_real_time"

        call open_pops_head(iunit, formpops, binpops)

        if (formpops) then
            ! this is the NON-binary read
            popsversion = FindPopsfileVersion(iunit)
            if (popsversion /= 4) then
                call stop_all(this_routine, "wrong POPSFILE version /= 4!")
            end if
            call ReadPopsHeadv4(iunit, 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, unused_triples, read_nnodes, read_walkers_on_nodes, PopBalanceBlocks)

            ! if popsfiles are stored in binary! there are seperate files for
            ! the header and the actual population stats
            if (iProcIndex == root) then
                call get_unique_filename('POPSFILEBIN', tIncrementPops, .false., &
                                         iPopsFileNoRead, popsfile)
                open(iunit, file=popsfile, status='old', form='unformatted')
            end if
        end if

    end subroutine read_popsfile_real_time

    subroutine readTimeEvolvedState()
        use PopsfileMod, only: FindPopsfileVersion, ReadPopsHeadv4, InitFCIMC_pops, &
        use semi_stoch_gen, only: init_semi_stochastic
        use semi_stoch_procs, only: end_semistoch
        use real_time_procs, only: reset_core_space
        use CalcData, only: ss_space_in
        use FciMCData, only: TotImagTime
        implicit none

        integer :: iunit, popsversion, iPopLenof_Sign, iPopNel, iPopIter, &
                   PopNIfD, PopNIfSgn, PopNIfFlag, PopNIfTot, &
                   PopBlockingIter, Popinum_runs, PopRandomHash(2056), &
                   read_nnodes, PopBalanceBlocks
        logical :: formpops, binpops, tPop64Bit, tPopHPHF, tPopLz
        integer(int64) :: iPopAllTotWalkers, read_walkers_on_nodes(0:nProcessors - 1)
        real(dp) :: PopDiagSft(inum_runs), read_tau, PopSumNoatHF(lenof_sign), &
                    read_psingles, read_pparallel
        HElement_t(dp) :: PopAllSumENum(inum_runs)
        integer :: ierr
        logical :: tStartedFromCoreGround
        character(255) :: rtPOPSFILE_name
        character(*), parameter :: this_routine = "readTimeEvolvedState"
        real(dp) :: unused_triples

        if (tSemiStochastic) then
            ! if semi-stochastic mode is enabled, it has to be disabled for read-in again
            ! as load balancing has to be performed
            call end_semistoch()
            call reset_core_space()
        end if

        binpops = .false.

        rtPOPSFILE_name = 'TIME_EVOLVED_POP'

        ! get the file containing the time evolved state
        call open_pops_head(iunit, formpops, binpops, rtPOPSFILE_name)

        popsversion = FindPopsfileVersion(iunit)
        if (popsversion /= 4) call stop_all(this_routine, "wrong popsfile version of TIME_EVOLVED_POP")

        call ReadPopsHeadv4(iunit, 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, unused_triples, read_nnodes, read_walkers_on_nodes, PopBalanceBlocks)

        ! at this point, we do not want to perturb the state and have no use for the
        ! pops_pert variable anymore -> deallocate it

        ! read in the hacked shift_damping
        if (.not. allocated(shift_damping)) allocate(shift_damping(inum_runs), stat=ierr)
        shift_damping = PopSumNoatHF(1:inum_runs)

        call clear_pops_pert(pops_pert)

        ! read in the time evolved state and use it as initial state
        call InitFCIMC_pops(iPopAllTotWalkers, PopNIfSgn, iPopNel, read_nnodes, &
                            read_walkers_on_nodes, pops_pert, &
                            PopBalanceBLocks, PopDiagSft, rtPOPSFILE_name)

        call set_initial_times(read_tau, TotImagTime, PopDiagSft(1))

        ! if we disabled semi-stochastic mode temporarily, reenable it now
        if (tSemiStochastic) call init_semi_stochastic(ss_space_in, tStartedFromCoreGround)

    end subroutine readTimeEvolvedState

    subroutine set_initial_times(real_time, imag_time, alpha)
        implicit none
        real(dp), intent(in) :: real_time, imag_time, alpha

        ! with the inclusion of dynamic alpha, both, the real and the imaginary part
        ! have to be stored in the time-evolved popsfile
        ! the
        elapsedImagTime = imag_time
        elapsedRealTime = real_time
        ! also, the previous value of alpha is to be loaded
        real_time_info%time_angle = alpha
    end subroutine set_initial_times

    subroutine equalize_initial_phase()
        use FciMCData, only: AllSumNoatHF
        use bit_rep_data, only: extract_sign
        use bit_reps, only: encode_sign
        use util_mod, only: near_zero
        implicit none

        integer :: signs(lenof_sign), iGf
        integer(int64) :: i
        real(dp) :: tmp_sgn(lenof_sign)

        signs = 1

        if(.not. near_zero(AllSumNoatHF(1))) then
            do i = 1, lenof_sign
                if (AllSumNoatHF(i) / AllSumNoatHF(1) < 0) then
                    signs(i) = -1
                    signs(i) = 1
                end if
            end do
        end if

        if (any(signs < 0)) then
            do i = 1, TotWalkers
                call extract_sign(CurrentDets(:, i), tmp_sgn)
                tmp_sgn = tmp_sgn * signs
                call encode_sign(CurrentDets(:, i), tmp_sgn)
            end do
            do iGf = 1, gf_count
                do i = 1, overlap_states(iGf)%nDets
                    call extract_sign(overlap_states(iGf)%dets(:, i), tmp_sgn)
                    tmp_sgn = tmp_sgn * signs
                    call encode_sign(overlap_states(iGf)%dets(:, i), tmp_sgn)
                end do
            end do
        end if

    end subroutine equalize_initial_phase

    ! here, we convert the real-space annihilation/creation operators to momentum space
    subroutine setup_momentum_operators()
        use SystemData, only: G1
        implicit none

        integer :: k(3), state, momentum_state, spin, i, j

        if (gf_type == -1) then
            momentum_state = pops_pert(1)%ann_orbs(1)
        else if (gf_type == 1) then
            momentum_state = pops_pert(1)%crtn_orbs(1)
            call stop_all("Equalize initial phase", "Invalid momentum space green's function")
        end if
        allocate(phase_factors(nBasis / 2))
        k = G1(momentum_state)%k
        spin = G1(momentum_state)%ms
        call clear_pops_pert(pops_pert)
        call clear_pops_pert(overlap_pert)
        allocate(pops_pert(nBasis / 2))
        allocate(overlap_pert(nBasis / 2))
        i = 0
        do state = 1, nBasis
            if (G1(state)%ms /= spin) cycle
            i = i + 1
            phase_factors(i) = 0
            do j = 1, 3
                phase_factors(i) = phase_factors(i) + G1(state)%k(j) * k(j)
            end do
            call setup_single_perturbation(i, state, gf_type)
        end do
    end subroutine setup_momentum_operators

    subroutine setup_single_perturbation(pos, orbital, type)
        implicit none

        integer, intent(in) :: pos, orbital, type

        if (type == 1) then
            overlap_pert(pos)%ncreate = 1
            pops_pert(pos)%ncreate = 1
            pops_pert(pos)%crtn_orbs(1) = orbital
            overlap_pert(pos)%crtn_orbs(1) = orbital
            call init_perturbation_creation(overlap_pert(pos))
            call init_perturbation_creation(pops_pert(pos))
            overlap_pert(pos)%nannihilate = 1
            pops_pert(pos)%nannihilate = 1
            pops_pert(pos)%ann_orbs(1) = orbital
            overlap_pert(pos)%ann_orbs(1) = orbital
            call init_perturbation_annihilation(overlap_pert(pos))
            call init_perturbation_annihilation(pops_pert(pos))
        end if
    end subroutine setup_single_perturbation


    subroutine dealloc_real_time_memory
        use replica_data, only: clean_iter_data
        implicit none

        integer :: ierr
        character(*), parameter :: this_routine = "dealloc_real_time_memory"

        if (numSnapshotOrbs > 0) deallocate(snapShotOrbs, stat=ierr)
        if (allocated(numCycShiftExcess)) deallocate(numCycShiftExcess, stat=ierr)
        deallocate(DiagVec, stat=ierr)
        call LogMemDealloc(this_routine, DiagVecTag)
        call clean_iter_data(second_spawn_iter_data)
        deallocate(shift_damping, stat=ierr)
        deallocate(temp_freeslot, stat=ierr)
        deallocate(current_overlap, stat=ierr)
        deallocate(gs_energy, stat=ierr)
        deallocate(dyn_norm_red, stat=ierr)
        deallocate(dyn_norm_psi, stat=ierr)
        deallocate(pert_norm, stat=ierr)
        deallocate(gf_overlap, stat=ierr)
        deallocate(TotPartsPeak, stat=ierr)
        deallocate(temp_det_list, stat=ierr)
        call clean_overlap_states()
        call clear_pops_pert(pops_pert)
        call clear_pops_pert(overlap_pert)

        if (allocated(tauCache)) deallocate(tauCache)

    end subroutine dealloc_real_time_memory


    subroutine clear_pops_pert(perturbs)
        use FciMCData, only: perturbation
        implicit none

        type(perturbation), intent(inout), allocatable :: perturbs(:)
        integer :: j

        if (allocated(perturbs)) then
            do j = 1, gf_count
                if (allocated(perturbs(j)%crtn_orbs)) deallocate(perturbs(j)%crtn_orbs)
                if (allocated(perturbs(j)%ann_orbs)) deallocate(perturbs(j)%ann_orbs)
            end do
        end if
    end subroutine clear_pops_pert


    subroutine read_from_contour_file(iunit)
        use CalcData, only: nmcyc
        use real_time_data, only: tauCache, alphaCache
        implicit none
        integer, intent(in) :: iunit
        integer :: i

        ! Then, the cache for the values of alpha and tau is allocated
        if (allocated(tauCache)) deallocate(tauCache)
        if (allocated(alphaCache)) deallocate(alphaCache)

        ! And then the values for alpha/tau
        open(iunit, file=trajFile, status='old')
        ! Here, we read in the timesteps and alphas that shall be used
        do i = 1, nmcyc
            read(iunit, *) tauCache(i), alphaCache(i)
        end do
    end subroutine read_from_contour_file


    subroutine read_in_trajectory
        use CalcData, only: nmcyc
        use real_time_data, only: tauCache, alphaCache
        implicit none
        integer :: i, iunit, eof
        real(dp) :: x, y
        logical :: checkTraj
        character(*), parameter :: this_routine = "read_in_trajectory"

        call get_unique_filename('tauContour', .false., .false., 0, trajFile)
        iunit = get_free_unit()
        inquire (file=trajFile, exist=checkTraj)
        if (.not. checkTraj) call stop_all(this_routine, "No tauContour file detected.")

        ! We first need to read in the number of cycles
        open(iunit, file=trajFile, status='old')
        nmcyc = 0
        ! check the number of lines
            read(iunit, *, iostat=eof) x, y
            if (eof /= 0) exit
            nmcyc = nmcyc + 1
        end do

        call read_from_contour_file(iunit)

    end subroutine read_in_trajectory


    subroutine truncate_initial_state
        use semi_stoch_procs, only: check_determ_flag
        use hash, only: hash_table_lookup, remove_hash_table_entry
        use FciMCData, only: HashIndex
        use bit_reps, only: nullify_ilut, decode_bit_det
        implicit none
        integer :: nI(nel), DetHash, PartInd
        integer(int64) :: i
        logical :: tSuccess

        do i = 1, TotWalkers
            if (.not. check_determ_flag(CurrentDets(:, i))) then
                call nullify_ilut(CurrentDets(:, i))
                call decode_bit_det(nI, CurrentDets(:, i))
                call hash_table_lookup(nI, CurrentDets(:, i), nifd, HashIndex, CurrentDets, &
                                       PartInd, DetHash, tSuccess)
                if (tSuccess) call remove_hash_table_entry(HashIndex, nI, PartInd)
            end if
        end do
    end subroutine truncate_initial_state


    subroutine truncate_overlap_states
        use semi_stoch_procs, only: check_determ_flag
        use FciMCData, only: HashIndex
        use bit_reps, only: nullify_ilut
        use hash, only: hash_table_lookup
        implicit none
        integer :: i, iGf, nI(nel), DetHash, PartInd
        logical :: tSuccess

        do iGf = 1, gf_count
            ! For each gf, truncate the corresponding overlap state
            do i = 1, overlap_states(iGf)%ndets
                call decode_bit_det(nI, overlap_states(iGf)%dets(:, i))
                call hash_table_lookup(nI, overlap_states(iGf)%dets(:, i), nifd, HashIndex, &
                                       CurrentDets, PartInd, DetHash, tSuccess)
                if (tSuccess) then
                    ! In principle, there are no non-core determinants left when this
                    ! is called, but we do not want to depend on the order in which the
                    ! truncation is carried out
                    if (.not. check_determ_flag(CurrentDets(:, PartInd))) call nullify_ilut( &
                        overlap_states(iGf)%dets(:, i))
                    call nullify_ilut(overlap_states(iGf)%dets(:, i))
                end if
            end do
        end do
    end subroutine truncate_overlap_states


    subroutine set_deterministic_flag(i)
        use FciMCData, only: HashIndex, core_run
        use bit_rep_data, only: flag_deterministic
        use hash, only: hash_table_lookup
        use bit_reps, only: set_flag
        use core_space_util, only: cs_replicas
        implicit none
        integer, intent(in) :: i
        integer ::  nI(nel), DetHash, PartInd
        logical :: tSuccess
        character(*), parameter :: this_routine = "set_deterministic_flag"

        ! For each core-state, we check if it is in the CurrentDets (which should
        ! always be the case in the initialization
        associate(rep => cs_replicas(core_run))
            call decode_bit_det(nI, rep%core_space(:, i))
            call hash_table_lookup(nI, rep%core_space(:, i), nifd, HashIndex, CurrentDets, PartInd, &
                                   DetHash, tSuccess)
            if (tSuccess) then
                ! And then we set the deterministic flag
                call set_flag(CurrentDets(:, PartInd), flag_deterministic(core_run))
                call stop_all(this_routine, "Deterministic state not present in CurrentDets")
            end if
        end associate
    end subroutine set_deterministic_flag


    subroutine initialize_corespace_construction
        use hash, only: init_hash_table
        use FciMCData, only: nWalkerHashes, MaxWalkersPart
        use real_time_data, only: ssht, core_space_buf, csbuf_size
        implicit none
        integer :: ierr

        ! allocate the buffer for the corespace
        allocate(core_space_buf(0:niftot, MaxWalkersPart), stat=ierr)
        ! set the position to 0
        csbuf_size = 0
        ! setup the buffer hashtable
        allocate(ssht(nWalkerHashes), stat=ierr)
        call init_hash_table(ssht)

    end subroutine initialize_corespace_construction

end module real_time_init