check_start_rdm Subroutine

    subroutine check_start_rdm()

        ! This routine checks if we should start filling the RDMs -
        ! and does so if we should.

        use rdm_general, only: realloc_SpawnedParts
        use LoggingData, only: tReadRDMs, tTransitionRDMs
        logical :: tFullVaryshift

        tFullVaryShift = .false.

        if (all(.not. tSinglePartPhase)) tFullVaryShift = .true.

        ! If we're reading in the RDMs we've already started accumulating them in a previous calculation
        ! We don't want to put in an arbitrary break now!
        if (tReadRDMs) IterRDMonFly = 0

        if (tFullVaryShift .and. ((Iter - maxval(VaryShiftIter)) == (IterRDMonFly + 1))) then
            ! IterRDMonFly is the number of iterations after the shift has changed that we want
            ! to fill the RDMs.  If this many iterations have passed, start accumulating the RDMs!

            IterRDMStart = Iter + PreviousCycles
            IterRDM_HF = Iter + PreviousCycles

            if (tEN2) tEN2Started = .true.

            ! We have reached the iteration where we want to start filling the RDM.
            if (tExplicitAllRDM) then
                ! Explicitly calculating all connections - expensive...
                ! TODO: why are explicit RDMs not working with replicas?
                if (tPairedReplicas) call stop_all('check_start_rdm', "Cannot yet do replica RDM sampling with explicit RDMs. &
                    & e.g Hacky bit in Gen_Hist_ExcDjs to make it compile.")

                tFillingExplicRDMonFly = .true.
                if (tHistSpawn) NHistEquilSteps = Iter
            else

                ! If we are load balancing, this will disable the load balancer
                ! so we should do a last-gasp balance at this point.
                if (tLoadBalanceBlocks .and. .not. tSemiStochastic) &
                    call adjust_load_balance(iter_data_fciqmc)

                extract_bit_rep_avsign => extract_bit_rep_avsign_norm
                ! By default - we will do a stochastic calculation of the RDM.
                tFillingStochRDMonFly = .true.
                if (tTransitionRDMs) tTransitionRDMsStarted = .true.

                call realloc_SpawnedParts()
                ! The SpawnedParts array now needs to carry both the spawned parts Dj, and also it's
                ! parent Di (and it's sign, Ci). - We deallocate it and reallocate it with the larger size.
                ! Don't need any of this if we're just doing HF_Ref_Explicit calculation.
                ! This is all done in the add_rdm_hfconnections routine.
            end if

            if (RDMExcitLevel == 1) then
                write(stdout, '(A)') 'Calculating the 1 electron density matrix on the fly.'
            else
                write(stdout, '(A)') 'Calculating the 2 electron density matrix on the fly.'
            end if
            write(stdout, '(A,I10)') 'Beginning to fill the RDMs during iteration', Iter
        end if

    end subroutine check_start_rdm