InitFCIMC_MP1 Subroutine

    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 (allocated(lookup_property_indexer)) then
                        call set_prop_vec_idx(DetIndex, lookup_property_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