ReadFromPopsfilePar Subroutine

    SUBROUTINE ReadFromPopsfilePar()

        use rdm_data, only: rdm_definitions

        LOGICAL :: exists, tBinRead
        Real(dp) :: NodeSumNoatHF(nProcessors)
        INTEGER :: WalkerstoReceive(nProcessors)
        integer(int64) :: AvWalkers
        integer(int64) :: TempCurrWalkers
        INTEGER :: TempInitWalkers, error, i, j, l, total, ierr, MemoryAlloc, Tag, Proc, CurrWalkers, ii
        INTEGER, DIMENSION(lenof_sign) :: TempSign
        integer(int64) :: TempSign64
        real(dp) :: RealTempSign(lenof_sign)
        integer(int64) :: iLutTemp64(0:nBasis / 64 + 1)
        INTEGER :: iLutTemp32(0:nBasis / 32 + 1)
        INTEGER(KIND=n_int) :: iLutTemp(0:NIfTot)
        INTEGER :: AvSumNoatHF, IntegerPart, TempnI(NEl), ExcitLevel
        integer :: NIfWriteOut, pos, orb, PopsVersion, iunit, run
        real(dp) :: r, FracPart, Gap, DiagSftTemp, tmp_dp
        HElement_t(dp) :: HElemTemp, HElemTemp2
        character(255) :: popsfile, FirstLine, stem
        character(len=24) :: junk, junk2, junk3, junk4
        LOGICAL :: tPop64BitDets, tPopHPHF, tPopLz, tPopInitiator
        integer(n_int) :: ilut_largest(0:NIfTot, inum_runs)
        real(dp) :: sign_largest(inum_runs)
        character(*), parameter :: this_routine = 'ReadFromPopsfilePar'

        write (stdout, *) "THIS IS THE POPSFILE ROUTINE WE'RE USING"

        if (lenof_sign /= 1) &
            call Stop_All(this_routine, "Popsfile V.2 does not work with " &
                                        // "complex walkers")

        PreviousCycles = 0    !Zero previous cycles
        SumENum = 0.0_dp
        TotParts = 0
        SumNoatHF = 0
        DiagSft = 0.0_dp
        Tag = 124             !Set Tag

        if (tPopsAlias) then
            stem = aliasStem
        else
            stem = 'POPSFILE'
        endif

        call get_unique_filename(stem, tIncrementPops, .false., iPopsFileNoRead, popsfile)
        iunit = get_free_unit()
        INQUIRE (FILE=popsfile, EXIST=exists)
        IF (exists) THEN
            open (iunit, FILE=popsfile, Status='old')
            tBinRead = .false.
        ELSE
            tBinRead = .true.
            call stop_all("READ PAR", "BINARY")
            call get_unique_filename('POPSFILEHEAD', tIncrementPops, .false., iPopsFileNoRead, popsfile)
            INQUIRE (FILE=popsfile, EXIST=exists)
            IF (.not. exists) THEN
                call get_unique_filename('POPSFILEBIN', tIncrementPops, .false., iPopsFileNoRead, popsfile)
                INQUIRE (FILE=popsfile, EXIST=exists)
                IF (.not. exists) THEN
                    CALL Stop_All(this_routine, "No POPSFILEs of any kind found.")
                ELSE
                    CALL Stop_All(this_routine, "POPSFILEBIN(.x) found, but POPSFILEHEAD(.x) also needed for header information")
                ENDIF
            ELSE
                call get_unique_filename('POPSFILEBIN', tIncrementPops, .false., iPopsFileNoRead, popsfile)
                INQUIRE (FILE=popsfile, EXIST=exists)
                IF (.not. exists) THEN
                    CALL Stop_All(this_routine, "POPSFILEHEAD(.x) found, but no POPSFILEBIN(.x) " &
                     //"for particle information - this is also needed")
                ELSE
                    call get_unique_filename('POPSFILEHEAD', tIncrementPops, .false., iPopsFileNoRead, popsfile)
                    open (iunit, FILE=popsfile, Status='old')
                ENDIF
            ENDIF
        ENDIF

        READ (iunit, '(a255)') FirstLine

        IF (INDEX(FirstLine, 'VERSION') == 0) THEN
!No version number to be found
            PopsVersion = 1
            REWIND (iunit)
        ELSE
            !Found version - which number is it?
            REWIND (iunit)
            READ (iunit, *) FirstLine, FirstLine, FirstLine, PopsVersion
        ENDIF
        write (stdout, "(A,I5,A)") "Version", PopsVersion, " POPSFILE detected"

!Read in initial data on processors which have a popsfile
        IF (PopsVersion == 2) THEN
            READ (iunit, '(A12,L5,A8,L5,A8,L5,A12,L5)') junk, tPop64BitDets, junk2, tPopHPHF, junk3, tPopLz, junk4, tPopInitiator
        ELSE
            write (stdout, '(A)') "Reading in from deprecated POPSFILE - assuming that parameters " &
             //"are the same as when POPSFILE was written"
        ENDIF
        READ (iunit, *) tmp_dp
        AllTotWalkers = int(tmp_dp, int64)
        READ (iunit, *) DiagSftTemp
        READ (iunit, *) AllSumNoatHF
        READ (iunit, *) AllSumENum
        READ (iunit, *) PreviousCycles

        IF (iProcIndex == Root) THEN
            IF (abs(iWeightPopRead) > 1.0e-12_dp) THEN
                write (stdout, "(A,I15,A,I4,A)") "Although ", AllTotWalkers, &
                    " configurations will be read in, only determinants with a weight of over ", iWeightPopRead, " will be stored."
            ENDIF
        ENDIF

        IF (.not. tWalkContGrow) THEN
!If we want the walker number to continue growing, then take the diagonal
! shift from the input, rather than the POPSFILE.
            DiagSft = DiagSftTemp
        ENDIF

        IF (abs(DiagSftTemp) < 1.0e-12_dp) THEN
            tWalkContGrow = .true.
            DiagSft = DiagSftTemp
        ENDIF

        IF (tBinRead) THEN
!Test for the end of the file.
!If this is not the end of the file, there is one more keyword that tells us
! the calculation had not entered variable shift mode yet.
!Want to put this test at the end of the non-binary file too.
            close (iunit)
            call get_unique_filename('POPSFILEBIN', tIncrementPops, .false., iPopsFileNoRead, popsfile)
            open (iunit, FILE=popsfile, Status='old', form='unformatted')
        ENDIF

        IF (iProcIndex == Root) THEN

            write (stdout, *) "Number of cycles in previous simulation: ", PreviousCycles
            IF (NEquilSteps > 0) THEN
                write (stdout, *) "Removing equilibration steps since reading in from POPSFILE."
                NEquilSteps = 0
            ENDIF
            IF (TZeroProjE) THEN
!Reset energy estimator
                write (stdout, *) "Resetting projected energy counters to zero..."
                AllSumENum = 0.0_dp
                AllSumNoatHF = 0
            ENDIF

!Need to calculate the number of walkers each node will receive...
            AvWalkers = NINT(AllTotWalkers / real(nProcessors, dp), int64)

!Divide up the walkers to receive for each node
            do i = 1, nProcessors - 1
                WalkerstoReceive(i) = int(AvWalkers)
            enddo
!The last processor takes the 'remainder'
            WalkerstoReceive(nProcessors) = int(AllTotWalkers - (AvWalkers * (nProcessors - 1)))

!Quick check to ensure we have all walkers accounted for
            total = 0
            do i = 1, nProcessors
                total = total + WalkerstoReceive(i)
            enddo
            if (total /= AllTotWalkers) then
                CALL Stop_All("ReadFromPopsfilePar", "All Walkers not accounted for when reading in from POPSFILE")
            endif

!InitWalkers needs to be reset for the culling criteria
            IF (.not. tWalkContGrow) THEN
!Now, let the total space allocated for storing walkers which have been read in to be equal to the initwalkers from the input file.
!                InitWalkers=AvWalkers
            ELSE
                TSinglePartPhase = .true.
            ENDIF
            SumENum = AllSumENum / REAL(nProcessors, dp)     !Divide up the SumENum over all processors
            AvSumNoatHF = int(AllSumNoatHF(1) / nProcessors) !This is the average Sumnoathf
            do i = 1, nProcessors - 1
                NodeSumNoatHF(i) = real(INT(AvSumNoatHF, int64), dp)
            enddo
            NodeSumNoatHF(nProcessors) = real(AllSumNoatHF(1) - INT((AvSumNoatHF * (nProcessors - 1)), int64), dp)

            ProjectionE = AllSumENum / real(AllSumNoatHF(1), dp)

!Reset the global variables
            AllSumENum = 0.0_dp
            AllSumNoatHF = 0

        ENDIF

        CALL MPIBarrier(error)  !Sync

!Now we need to scatter the WalkerstoReceive to each node, and allocate the desired memory to each node...
!Broadcast info which needs to go to all processors
        CALL MPIBCast(DiagSft)
        CALL MPIBCast(SumENum)
        CALL MPIBCast(InitWalkers)
        CALL MPIBCast(NEquilSteps)
        CALL MPIBCast(NShiftEquilSteps)
        CALL MPIBCast(TSinglePartPhase)
!        CALL MPI_BCast(tChangenProcessors,1,MPI_LOGICAL,root,MPI_COMM_WORLD,error)
!Scatter the number of walkers each node will receive to TempInitWalkers, and
! the SumNoatHF for each node which is distributed approximatly equally
        CALL MPIScatter(WalkerstoReceive, TempInitWalkers, error)
        CALL MPIScatter(NodeSumNoatHF, SumNoatHF(1), error)

        IF (MemoryFacPart <= 1.0_dp) THEN
            write (stdout, *) 'MemoryFacPart must be larger than 1.0 when reading in a POPSFILE - increasing it to 1.50.'
            MemoryFacPart = 1.50
        ENDIF

!Now we want to allocate memory on all nodes.
!InitWalkers here is simply the average number of walkers per node, not actual
        MaxWalkersPart = NINT(MemoryFacPart * (NINT(InitWalkers * ScaleWalkers)))
        MaxSpawned = NINT(MemoryFacSpawn * (NINT(InitWalkers * ScaleWalkers * inum_runs)))

        Gap = REAL(MaxSpawned, dp) / REAL(nProcessors, dp)
        do i = 0, nProcessors - 1
            InitialSpawnedSlots(i) = NINT(Gap * i) + 1
        enddo
!ValidSpawndList now holds the next free position in the newly-spawned list,
! but for each processor.
        ValidSpawnedList(:) = InitialSpawnedSlots(:)

        CALL MPIBarrier(error)  !Sync

!Allocate memory to hold walkers at least temporarily
        ALLOCATE (WalkVecDets(0:NIfTot, MaxWalkersPart), stat=ierr, source=0_n_int)
        CALL LogMemAlloc('WalkVecDets', MaxWalkersPart * (NIfTot + 1), size_n_int, this_routine, WalkVecDetsTag, ierr)
        MemoryAlloc = (NIfTot + 1) * MaxWalkersPart * size_n_int    !Memory Allocated in bytes

!Just allocating this here, so that the SpawnParts arrays can be used for sorting the determinants when using direct annihilation.
        write (stdout, "(A,I12,A)") " Spawning vectors allowing for a total of ", MaxSpawned, &
            " particles to be spawned in any one iteration."
        ALLOCATE (SpawnVec(0:NIfTot, MaxSpawned), stat=ierr, source=0_n_int)
        CALL LogMemAlloc('SpawnVec', MaxSpawned * (NIfTot + 1), size_n_int, this_routine, SpawnVecTag, ierr)
        ALLOCATE (SpawnVec2(0:NIfTot, MaxSpawned), stat=ierr, source=0_n_int)
        CALL LogMemAlloc('SpawnVec2', MaxSpawned * (NIfTot + 1), size_n_int, this_routine, SpawnVec2Tag, ierr)
        MemoryAlloc = MemoryAlloc + MaxSpawned * (1 + NIfTot) * 2 * size_n_int

!Point at correct spawning arrays
        SpawnedParts => SpawnVec
        SpawnedParts2 => SpawnVec2
!Allocate pointer to the correct walker array...
        CurrentDets => WalkVecDets

        ! 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(rdm_definitions%nrdms_standard, rdm_definitions%nrdms_transition)

! The hashing will be different in the new calculation from the one where the
!  POPSFILE was produced, this means we must recalculate the processor each
! determinant wants to go to.
! This is done by reading in all walkers to the root and then distributing
! them in the same way as the spawning steps are done - by finding the
!  determinant and sending it there.
        IF ((PopsVersion /= 1) .and. tHPHF .and. (.not. tPopHPHF)) THEN
            CALL Stop_All(this_routine, "HPHF on, but HPHF was not used in creation of the POPSFILE")
        ENDIF
        IF ((PopsVersion /= 1) .and. tFixLz .and. (.not. tPopLz)) THEN
            CALL Stop_All(this_routine, "Lz on, but Lz was not used in creation of the POPSFILE")
        ENDIF
        IF ((PopsVersion /= 1) .and. (.not. tHPHF) .and. tPopHPHF) THEN
            CALL Stop_All(this_routine, "HPHF off, but HPHF was used for creation of the POPSFILE")
        ENDIF
        IF ((PopsVersion /= 1) .and. (.not. tFixLz) .and. tPopLz) THEN
            CALL Stop_All(this_routine, "Lz off, but Lz was used for creation of the POPSFILE")
        ENDIF
        IF (PopsVersion == 1) THEN
            tPop64BitDets = .false.
            NIfWriteOut = nBasis / 32
        ELSE
            IF (.not. tPop64BitDets) THEN
                NIfWriteOut = nBasis / 32
            ELSE
                NIfWriteOut = nBasis / 64
            ENDIF
        ENDIF

        CurrWalkers = 0
        sign_largest = 0
        ilut_largest = 0
        do i = 1, int(AllTotWalkers)
            iLutTemp(:) = 0
            IF (PopsVersion /= 1) THEN
                IF (tBinRead) THEN
                    IF (tPop64BitDets) THEN
                        READ (iunit) iLutTemp64(0:NIfWriteOut), TempSign
                    ELSE
                        READ (iunit) iLutTemp32(0:NIfWriteOut), TempSign
                    ENDIF
                ELSE
                    IF (tPop64BitDets) THEN
                        READ (iunit, *) iLutTemp64(0:NIfWriteOut), TempSign
                    ELSE
                        READ (iunit, *) iLutTemp32(0:NIfWriteOut), TempSign
                    ENDIF
                ENDIF
            ELSE
                !POPSFILE v. 1 only printed out 32 bit determinant strings.
                IF (tBinRead) THEN
                    READ (iunit) iLutTemp32(0:NIfWriteOut), TempSign
                ELSE
                    READ (iunit, *) iLutTemp32(0:NIfWriteOut), TempSign
                ENDIF
            ENDIF
            do j = 1, lenof_sign
                TempSign64 = TempSign(j)
                RealTempSign(j) = transfer(TempSign64, RealTempSign(j))
            enddo

#ifdef INT64_
            if (.not. tPop64BitDets) then
                ! If we are using 64 bit integers, but have read in 32 bit
                ! integers, then we need to convert them.
                call stop_all(this_routine, &
                              "WD: I removed nifdbo as i assumes its always = nifd.. this is a problem here")
                ! do ii=0,nBasis/32
                !     do j=0,31
                !         if(btest(iLutTemp32(ii),j)) then
                !             orb=(ii*32)+j+1
                !             pos=(orb-1)/bits_n_int
                !             iLutTemp(pos)=ibset(iLutTemp(pos),mod(orb-1,bits_n_int))
                !         endif
                !     enddo
                ! enddo
                ! iLutTemp(NIfD+1:NIfDBO) = iLutTemp64(ii:NIfWriteOut)
            else
                iLutTemp(0:nifd) = iLutTemp64(0:nifd)
            endif

#else
            ! If we are using 32 bit integers, but have read in 64 bit
            ! integers, then we need to convert them.
            if (tPop64BitDets) then
                call stop_all(this_routine, &
                              "WD: I removed nifdbo as i assumes its always = nifd.. this is a problem here")
                ! do ii=0,nBasis/64
                !     do j=0,63
                !         if(btest(iLutTemp64(ii),j)) then
                !             orb=(ii*64)+j+1
                !             pos=(orb-1)/bits_n_int
                !             iLutTemp(pos)=ibset(iLutTemp(pos),mod(orb-1,bits_n_int))
                !         endif
                !     enddo
                ! enddo
                ! iLutTemp(NIfD+1:NIfDBO) = iLutTemp32(ii:NIfWriteOut)
            else
                iLutTemp(0:nifd) = iLutTemp32(0:nifd)
            endif

#endif
            call decode_bit_det(TempnI, iLutTemp)
            Proc = DetermineDetNode(nel, TempnI, 0)
            IF ((Proc == iNodeIndex) .and. (abs(RealTempSign(1)) >= iWeightPopRead)) THEN
                CurrWalkers = CurrWalkers + 1
                !Do not need to send a flag here...
                call encode_bit_rep(CurrentDets(:, CurrWalkers), iLutTemp(0:nifd), RealTempSign, 0)
                !TODO: Add flag for complex walkers to read in both

!>>>"                CurrentH(1:1+2*lenof_sign,CurrWalkers)=CurrentHEntry(1:1+2*lenof_sign)
            ENDIF

            ! Keep track of what the most highly weighted determinant is
            do run = 1, inum_runs
                if (abs(RealTempSign(run)) > sign_largest(run)) then
                    sign_largest(run) = abs(RealTempSign(run))
                    ilut_largest(:, run) = iLutTemp
                endif
            end do
        enddo
        close (iunit)
        TempCurrWalkers = int(CurrWalkers, int64)

        ! Sort the lists so that they are in order if we change the number
        ! of processors.
        call sort(currentdets(:, 1:CurrWalkers), &
                  global_determinant_data(:, 1:CurrWalkers))

        ! Check that the bit-det comparisons agree that it is in order.
        do i = 2, currwalkers
            if (DetBitLT(CurrentDets(:, i), CurrentDets(:, i - 1), nifd) == 1) then
                print *, 'Walkers: ', i - 1, i
                print *, 'bit reps: '
                print *, currentdets(:, i - 1)
                print *, currentdets(:, i)
                call stop_all(this_routine, 'Main list out of order')
            endif
        enddo

        CALL MPIBarrier(error)  !Sync
        CALL MPIAllReduce(TempCurrWalkers, MPI_SUM, AllTotWalkers)

       IF (iProcIndex == root) write (stdout, '(I10,A)') INT(AllTotWalkers, int64), " configurations read in from POPSFILE and distributed."

        IF (abs(ScaleWalkers - 1) > 1.0e-12_dp) THEN

            write (stdout, *) "Rescaling walkers by a factor of: ", ScaleWalkers

! CurrWalkers is the number of determinants on a particular node, AllTotWalkers is the total over all nodes.
            IntegerPart = INT(ScaleWalkers)
            FracPart = ScaleWalkers - REAL(IntegerPart, dp)

            do l = 1, CurrWalkers
                call extract_sign(CurrentDets(:, l), RealTempSign)
                RealTempSign = RealTempSign * IntegerPart
                r = genrand_real2_dSFMT()
                IF (r < FracPart) THEN
!Stochastically create another particle
                    IF (RealTempSign(1) < 0) THEN
                        RealTempSign(1) = RealTempSign(1) - 1
                    ELSE
                        RealTempSign(1) = RealTempSign(1) + 1
                    ENDIF
                ENDIF
                call encode_sign(CurrentDets(:, l), RealTempSign)
            enddo

            InitWalkers = NINT(InitWalkers * ScaleWalkers)  !New (average) number of initial particles for culling criteria
!Other parameters don't change (I think) because the number of determinants isn't changing.
            TotWalkers = CurrWalkers
            TotWalkersOld = CurrWalkers
            IF (iProcIndex == root) THEN
!                AllTotWalkers=TotWalkers
                AllTotWalkersOld = AllTotWalkers
                write (stdout, '(A,I10)') " Number of initial walkers on this processor is now: ", INT(TotWalkers, int64)
            ENDIF

        ELSE
!We are not scaling the number of walkers...

            TotWalkers = CurrWalkers
            TotWalkersOld = CurrWalkers
            IF (iProcIndex == root) THEN
!                AllTotWalkers=TotWalkers
                AllTotWalkersOld = AllTotWalkers
                write (stdout, '(A,I10)') " Number of initial walkers on this processor is now: ", INT(TotWalkers, int64)
            ENDIF

        ENDIF

        write (stdout, *) "Initial Diagonal Shift (ECorr guess) is now: ", DiagSft(1)
        write (stdout, "(A,F14.6,A)") " Initial memory (without excitgens) consists of : ", REAL(MemoryAlloc, dp) / 1048576.0_dp, " Mb"
        write (stdout, *) "Initial memory allocation successful..."
        write (stdout, *) "Excitgens will be regenerated when they are needed..."
        CALL neci_flush(stdout)

        ! If we are changing the reference determinant to the largest
        ! weighted one in the file, do it here
        if (tReadPopsChangeRef .or. tReadPopsRestart) then
            do run = 1, inum_runs

                ! If using the same reference for all runs, then don't consider
                ! the populations seperately.
                if (run /= 1 .and. .not. tReplicaReferencesDiffer) &
                    exit

                if (.not. DetBitEq(ilut_largest(:, run), iLutRef(:, run))) &
                    call update_run_reference(ilut_largest(:, run), run)

            end do
        endif

        ! Calculate the stored data for the particlse that have been read in
        TotParts = 0
        do j = 1, int(TotWalkers)
            call decode_bit_det(TempnI, currentDets(:, j))
            ! note on GUGA: here it is fine since working out excit level = 0
            ! works and thats all what is necessary here!
            Excitlevel = FindBitExcitLevel(iLutHF, CurrentDets(:, j), 2)
            IF (Excitlevel == 0) THEN
                call set_det_diagH(j, 0.0_dp)
            ELSE
                HElemTemp = get_diagonal_matel(TempnI, CurrentDets(:,j))
                call set_det_diagH(j, real(HElemTemp, dp) - Hii)
            ENDIF
            HElemTemp2 = get_off_diagonal_matel (TempnI, CurrentDets(:,j))
            call set_det_offdiagH(j, HElemTemp2)
            call store_decoding(j, TempnI)
            call extract_sign(CurrentDets(:, j), RealTempSign)
            TotParts(1) = TotParts(1) + abs(RealTempSign(1))
            TotParts(inum_runs) = TotParts(1)
        enddo

        CALL MPIBarrier(error)  !Sync
        CALL MPIReduce(TotParts, MPI_SUM, AllTotParts)

        IF (iProcIndex == root) then
            AllTotPartsOld = AllTotParts
            iter_data_fciqmc%tot_parts_old = AllTotParts
        endif
        write (stdout, '(A,i20)') ' The total number of particles read from the POPSFILE is: ', AllTotParts(1)

        if (tReadPopsRestart) then
            tPopsAlreadyRead = .true.
            call ChangeRefDet(ProjEDet(:, 1))
            tPopsAlreadyRead = .false.
        endif

    END SUBROUTINE ReadFromPopsfilePar