subroutine init_kp_fciqmc(kp)
use CalcData, only: tSemiStochastic, tUseRealCoeffs, AvMCExcits, tPairedReplicas
use fcimc_initialisation, only: SetupParameters, InitFCIMCCalcPar, init_fcimc_fn_pointers
use FciMCData, only: tPopsAlreadyRead, nWalkerHashes, SpawnVecKP
use FciMCData, only: SpawnVecKP2, MaxSpawned
use FciMCData, only: SpawnedPartsKP, SpawnedPartsKP2, MaxWalkersUncorrected
use FciMCData, only: iter_data_fciqmc, spawn_ht, nhashes_spawn, tReplicaReferencesDiffer
use hash, only: init_hash_table
use LoggingData, only: tFCIMCStats2
use Parallel_neci, only: MPIBarrier
use SystemData, only: G1, tHeisenberg, tAllSymSectors
use util_mod, only: int_fmt
type(kp_fciqmc_data), intent(inout) :: kp
integer :: ierr, krylov_vecs_mem, krylov_ht_mem, matrix_mem, spawn_ht_mem, i
character(len=*), parameter :: t_r = "init_kp_fciqmc"
! Checks.
if (.not. tUseRealCoeffs) call stop_all(t_r, 'kp-fciqmc can only be run using &
&real coefficients).')
if (tExactHamil .and. nProcessors /= 1) call stop_all(t_r, 'The exact-hamil &
&option can only be used when running with one processor.')
if (theisenberg .and. tAllSymSectors) call stop_all(t_r, 'The option to use all &
&symmetry sectors at once has not been implemented with the Heisenberg model.')
if (n_int == 4) call stop_all(t_r, 'Use of RealCoefficients does not work with 32 bit &
&integers due to the use of the transfer operation from dp reals to 64 bit integers.')
if (tExcitedStateKP) then
if (tPairedReplicas .and. inum_runs /= 2 * kp%nvecs) then
call stop_all(t_r, 'When using the ExcitedStateKP option with paired replicas, the &
&number of replicas must be twice the number of states to be calculated.')
else if ((.not. tPairedReplicas) .and. inum_runs /= kp%nvecs) then
call stop_all(t_r, 'When using the ExcitedStateKP option without paired replicas, the &
&number of replicas must be equal to the number of states to be calculated.')
end if
end if
! Call external NECI initialisation routines.
tPopsAlreadyRead = .false.
call SetupParameters()
call InitFCIMCCalcPar()
call init_fcimc_fn_pointers()
write(stdout, '(/,12("="),1x,a9,1x,12("="))') "KP-FCIQMC"
if (tSemiStochastic .and. (.not. tFullyStochasticHamil)) then
tSemiStochasticKPHamil = .true.
else
tSemiStochasticKPHamil = .false.
end if
#if defined(PROG_NUMRUNS_)
! The number of *repeated* replicas for a particular state.
! i.e. if tExcitedState = .true. then this is the number of repeated
! replicas for each excited state. If tExcitedState = .false. (doing
! the old KP-FCIQMC algorithm), this is the number of repeats for the
! KP-FCIQMC wave function.
#ifndef CMPLX_
if (tPairedReplicas) then
lenof_sign_kp = 2
else
lenof_sign_kp = 1
end if
#else
if (tPairedReplicas) then
lenof_sign_kp = 4
else
lenof_sign_kp = 2
end if
#endif
#endif
! The number of elements required to store all replicas of all Krylov vectors.
lenof_all_signs = lenof_sign_kp * kp%nvecs
! The total length of a bitstring containing all Krylov vectors.
! +1 one is for the flag integer
NIfTotKP = nifd + lenof_all_signs + 1
! Set up kp_ind_* arrays.
allocate(kp_ind_1(kp%nvecs))
allocate(kp_ind_2(kp%nvecs))
if (tPairedReplicas) then
do i = 1, kp%nvecs
kp_ind_1(i) = 2 * i - 1
kp_ind_2(i) = 2 * i
end do
else
do i = 1, kp%nvecs
kp_ind_1(i) = i
kp_ind_2(i) = i
end do
end if
if (.not. tExcitedStateKP) then
! Allocate all of the KP arrays.
nhashes_kp = nWalkerHashes
TotWalkersKP = 0
krylov_vecs_length = nint(MaxWalkersUncorrected * memory_factor_kp * kp%nvecs)
nkrylov_amp_elems_tot = lenof_sign * kp%nvecs * krylov_vecs_length
! Allocate the krylov_vecs array.
! The number of MB of memory required to allocate krylov_vecs.
krylov_vecs_mem = krylov_vecs_length * (NIfTotKP + 1) * size_n_int / 1000000
write(stdout, '(a73,'//int_fmt(krylov_vecs_mem, 1)//')') "About to allocate array to hold all Krylov vectors. &
&Memory required (MB):", krylov_vecs_mem
write(stdout, '(a13)', advance='no') "Allocating..."; call neci_flush(stdout)
allocate(krylov_vecs(0:NIfTotKP, krylov_vecs_length), stat=ierr)
if (ierr /= 0) then
write(stdout, '(1x,a11,1x,i5)') "Error code:", ierr
call stop_all(t_r, "Error allocating krylov_vecs array.")
else
write(stdout, '(1x,a5)') "Done."
end if
call neci_flush(stdout)
krylov_vecs = 0_n_int
! Allocate the krylov_helems array.
! The number of MB of memory required to allocate krylov_helems.
krylov_vecs_mem = krylov_vecs_length * size_n_int / 1000000
write(stdout, '(a103,'//int_fmt(krylov_vecs_mem, 1)//')') "About to allocate array to hold diagonal Hamiltonian &
&elements for Krylov vectors. Memory required (MB):", krylov_vecs_mem
write(stdout, '(a13)', advance='no') "Allocating..."; call neci_flush(stdout)
allocate(krylov_helems(krylov_vecs_length), stat=ierr)
if (ierr /= 0) then
write(stdout, '(1x,a11,1x,i5)') "Error code:", ierr
call stop_all(t_r, "Error allocating krylov_helems array.")
else
write(stdout, '(1x,a5)') "Done."
end if
call neci_flush(stdout)
krylov_helems = 0.0_dp
! Allocate the hash table to krylov_vecs.
! The number of MB of memory required to allocate krylov_vecs_ht.
! Each node requires 16 bytes.
krylov_ht_mem = nhashes_kp * 16 / 1000000
write(stdout, '(a78,'//int_fmt(krylov_ht_mem, 1)//')') "About to allocate hash table to the Krylov vector array. &
&Memory required (MB):", krylov_ht_mem
write(stdout, '(a13)', advance='no') "Allocating..."; call neci_flush(stdout)
allocate(krylov_vecs_ht(nhashes_kp), stat=ierr)
if (ierr /= 0) then
write(stdout, '(1x,a11,1x,i5)') "Error code:", ierr
call stop_all(t_r, "Error allocating krylov_vecs_ht array.")
else
write(stdout, '(1x,a5)') "Done."
write(stdout, '(a106)') "Note that the hash table uses linked lists, and the memory usage will &
&increase as further nodes are added."
end if
call init_hash_table(krylov_vecs_ht)
allocate(SpawnVecKP(0:IlutBits%ind_pop + lenof_all_signs - 1, MaxSpawned), stat=ierr)
allocate(SpawnVecKP2(0:IlutBits%ind_pop + lenof_all_signs - 1, MaxSpawned), stat=ierr)
SpawnVecKP(:, :) = 0_n_int
SpawnVecKP2(:, :) = 0_n_int
SpawnedPartsKP => SpawnVecKP
SpawnedPartsKP2 => SpawnVecKP2
if (tSemiStochastic) then
associate(rep => cs_replicas(core_run))
allocate(partial_determ_vecs_kp(lenof_all_signs, rep%determ_sizes(iProcIndex)), stat=ierr)
allocate(full_determ_vecs_kp(lenof_all_signs, rep%determ_space_size), stat=ierr)
end associate
partial_determ_vecs_kp = 0.0_dp
full_determ_vecs_kp = 0.0_dp
end if
end if
! Allocate the hash table to the spawning array.
! The number of MB of memory required to allocate spawn_ht.
! Each node requires 16 bytes.
nhashes_spawn = int(0.8 * MaxSpawned)
spawn_ht_mem = nhashes_spawn * 16 / 1000000
write(stdout, '(a78,'//int_fmt(spawn_ht_mem, 1)//')') "About to allocate hash table to the spawning array. &
&Memory required (MB):", spawn_ht_mem
write(stdout, '(a13)', advance='no') "Allocating..."; call neci_flush(stdout)
allocate(spawn_ht(nhashes_spawn), stat=ierr)
if (ierr /= 0) then
write(stdout, '(1x,a11,1x,i5)') "Error code:", ierr
call stop_all(t_r, "Error allocating spawn_ht array.")
else
write(stdout, '(1x,a5)') "Done."
write(stdout, '(a106)') "Note that the hash table uses linked lists, and the memory usage will &
&increase as further nodes are added."
end if
call init_hash_table(spawn_ht)
! (2*kp%nrepeats+16) arrays with (kp%nvecs**2) 8-byte elements each.
matrix_mem = (2 * kp%nrepeats + 16) * (kp%nvecs**2) * 8 / 1000000
write(stdout, '(a66,'//int_fmt(matrix_mem, 1)//')') "About to allocate various subspace matrices. &
&Memory required (MB):", matrix_mem
write(stdout, '(a13)', advance='no') "Allocating..."
call neci_flush(stdout)
if (tOverlapPert) then
allocate(pert_overlaps(kp%nvecs))
allocate(kp_all_pert_overlaps(kp%nvecs))
end if
allocate(kp_hamil_mean(kp%nvecs, kp%nvecs))
allocate(kp_overlap_mean(kp%nvecs, kp%nvecs))
allocate(kp_hamil_se(kp%nvecs, kp%nvecs))
allocate(kp_overlap_se(kp%nvecs, kp%nvecs))
allocate(kp_overlap_eigv(kp%nvecs))
allocate(kp_init_overlaps(kp%nvecs))
allocate(kp_overlap_eigenvecs(kp%nvecs, kp%nvecs))
allocate(kp_transform_matrix(kp%nvecs, kp%nvecs))
allocate(kp_inter_matrix(kp%nvecs, kp%nvecs))
allocate(kp_eigenvecs_krylov(kp%nvecs, kp%nvecs))
write(stdout, '(1x,a5)') "Done."
! If av_mc_excits_kp hasn't been set by the user, just use AvMCExcits.
if (near_zero(av_mc_excits_kp)) av_mc_excits_kp = AvMCExcits
! Initialize
kp_overlap_mean = 0.0_dp
kp_hamil_mean = 0.0_dp
MaxSpawnedEachProc = int(0.88_dp * real(MaxSpawned, dp) / nProcessors)
if (tFCIMCStats2) then
call write_fcimcstats2(iter_data_fciqmc, initial=.true.)
else
call WriteFciMCStatsHeader()
end if
call MPIBarrier(ierr)
! Store the initial state of tSinglePartPhase so that we can stop the
! shift from varying on subsequent repeats.
tSinglePartPhaseKPInit = tSinglePartPhase
! Allow different replicas to have different references in this case.
if (tExcitedStateKP) tReplicaReferencesDiffer = .true.
end subroutine init_kp_fciqmc