” CurrentH(1:1+2lenof_sign,CurrWalkers)=CurrentHEntry(1:1+2lenof_sign)
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)
IF (ierr /= 0) CALL Stop_All(this_routine, 'Problem allocating WalkVecDets array.')
CALL LogMemAlloc('WalkVecDets', MaxWalkersPart * (NIfTot + 1), size_n_int, this_routine, WalkVecDetsTag, ierr)
WalkVecDets(:, :) = 0
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)
CALL LogMemAlloc('SpawnVec', MaxSpawned * (NIfTot + 1), size_n_int, this_routine, SpawnVecTag, ierr)
SpawnVec(:, :) = 0
ALLOCATE (SpawnVec2(0:NIfTot, MaxSpawned), stat=ierr)
CALL LogMemAlloc('SpawnVec2', MaxSpawned * (NIfTot + 1), size_n_int, this_routine, SpawnVec2Tag, ierr)
SpawnVec2(:, :) = 0
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