subroutine perform_kp_fciqmc(kp)
use orthogonalise, only: calc_replica_overlaps
type(kp_fciqmc_data), intent(inout) :: kp
integer :: iiter, idet, ireplica, ispawn, ierr, err
integer :: iconfig, irepeat, ivec, nlowdin
integer :: nspawn, parent_flags, unused_flags
integer :: ex_level_to_ref, ex_level_to_hf
integer :: TotWalkersNew, determ_ind, ic, ex(2, maxExcit), ms_parent
integer :: nI_parent(nel), nI_child(nel), MaxIndex
integer(n_int) :: ilut_child(0:NIfTot)
integer(n_int), pointer :: ilut_parent(:)
real(dp) :: prob, unused_rdm_real, parent_hdiag
real(dp) :: child_sign(lenof_sign), parent_sign(lenof_sign)
real(dp) :: unused_sign(lenof_sign), precond_fac
real(dp), allocatable :: lowdin_evals(:, :)
logical :: tChildIsDeterm, tParentIsDeterm, tParentUnoccupied
logical :: tParity, tSingBiasChange, tWritePopsFound
HElement_t(dp) :: HElGen, parent_hoffdiag
! Stores of the overlap and projected Hamiltonian matrices.
real(dp), pointer :: overlap_matrices(:, :, :)
real(dp), pointer :: hamil_matrices(:, :, :)
! Pointers to the matrices for a given repeat only.
real(dp), pointer :: overlap_matrix(:, :)
real(dp), pointer :: hamil_matrix(:, :)
! Variables to hold information output for the test suite.
real(dp) :: s_sum, h_sum
integer(n_int) :: int_sign(lenof_all_signs)
real(dp) :: test_sign(lenof_all_signs)
type(ll_node), pointer :: temp_node
! Unused factor
precond_fac = 1.0_dp
call init_kp_fciqmc(kp)
if (.not. tAllSymSectors) ms_parent = lms
allocate(overlap_matrices(kp%nvecs, kp%nvecs, kp%nrepeats), stat=ierr)
allocate(hamil_matrices(kp%nvecs, kp%nvecs, kp%nrepeats), stat=ierr)
allocate(lowdin_evals(kp%nvecs, kp%nvecs), stat=ierr)
overlap_matrices = 0.0_dp
hamil_matrices = 0.0_dp
outer_loop: do iconfig = 1, kp%nconfigs
do irepeat = 1, kp%nrepeats
! Point to the regions of memory where the projected Hamiltonian
! and overlap matrices for this repeat will be accumulated and stored.
overlap_matrix => overlap_matrices(:, :, irepeat)
hamil_matrix => hamil_matrices(:, :, irepeat)
overlap_matrix(:, :) = 0.0_dp
hamil_matrix(:, :) = 0.0_dp
call init_kp_fciqmc_repeat(iconfig, irepeat, kp%nrepeats, kp%nvecs, iter_data_fciqmc)
do ivec = 1, kp%nvecs
! Copy the current state of CurrentDets to krylov_vecs.
call store_krylov_vec(ivec)
! Calculate the overlap of the perturbed ground state vector
! with the new Krylov vector, if requested.
if (tOverlapPert) call calc_perturbation_overlap(ivec)
do iiter = 1, kp%niters(ivec)
call set_timer(walker_time)
iter = iter + 1
call init_kp_fciqmc_iter(iter_data_fciqmc, determ_ind)
!if (iter < 10) then
! write(stdout,*) "CurrentDets before:"
! do idet = 1, int(TotWalkers)
! call extract_bit_rep(CurrentDets(:, idet), nI_parent, parent_sign, unused_flags, &
! fcimc_excit_gen_store)
! if (tUseFlags) then
! write(stdout,'(i7, i12, 4x, f18.7, 4x, f18.7, 4x, l1)') idet, CurrentDets(0,idet), parent_sign, &
! test_flag(CurrentDets(:,idet), flag_deterministic)
! else
! write(stdout,'(i7, i12, 4x, f18.7, 4x, f18.7)') idet, CurrentDets(0,idet), parent_sign
! end if
! end do
! write(stdout,"(A)") "Hash Table: "
! do idet = 1, nWalkerHashes
! temp_node => HashIndex(idet)
! if (temp_node%ind /= 0) then
! write(stdout,'(i9)',advance='no') idet
! do while (associated(temp_node))
! write(stdout,'(i9)',advance='no') temp_node%ind
! temp_node => temp_node%next
! end do
! write(stdout,'()',advance='yes')
! end if
! end do
!end if
do idet = 1, int(TotWalkers)
! The 'parent' determinant from which spawning is to be attempted.
ilut_parent => CurrentDets(:, idet)
parent_flags = 0_n_int
! Indicate that the scratch storage used for excitation generation from the
! same walker has not been filled (it is filled when we excite from the first
! particle on a determinant).
fcimc_excit_gen_store%tFilled = .false.
call extract_bit_rep(ilut_parent, nI_parent, parent_sign, unused_flags, &
idet, fcimc_excit_gen_store)
ex_level_to_ref = FindBitExcitLevel(iLutRef(:, 1), ilut_parent, &
max_calc_ex_level)
if (tRef_Not_HF) then
ex_level_to_hf = FindBitExcitLevel(iLutHF_true, ilut_parent, max_calc_ex_level)
else
ex_level_to_hf = ex_level_to_ref
end if
tParentIsDeterm = check_determ_flag(ilut_parent, core_run)
tParentUnoccupied = IsUnoccDet(parent_sign)
! If this determinant is in the deterministic space then store the relevant
! data in arrays for later use.
if (tParentIsDeterm) then
! Store the index of this state, for use in annihilation later.
associate(rep => cs_replicas(core_run))
rep%indices_of_determ_states(determ_ind) = idet
! Add the amplitude to the deterministic vector.
rep%partial_determ_vecs(:, determ_ind) = parent_sign
end associate
determ_ind = determ_ind + 1
! The deterministic states are always kept in CurrentDets, even when
! the amplitude is zero. Hence we must check if the amplitude is zero
! and, if so, skip the state.
if (tParentUnoccupied) cycle
end if
! If this slot is unoccupied (and also not a core determinant) then add it to
! the list of free slots and cycle.
if (tParentUnoccupied) then
iEndFreeSlot = iEndFreeSlot + 1
FreeSlot(iEndFreeSlot) = idet
cycle
end if
! The current diagonal matrix element is stored persistently.
parent_hdiag = det_diagH(idet)
parent_hoffdiag = det_offdiagH(idet)
if (tTruncInitiator) call CalcParentFlag(idet, parent_flags)
call SumEContrib(nI_parent, ex_level_to_ref, parent_sign, ilut_parent, &
parent_hdiag, parent_hoffdiag, 1.0_dp, tPairedReplicas, idet)
! If we're on the Hartree-Fock, and all singles and
! doubles are in the core space, then there will be
! no stochastic spawning from this determinant, so
! we can the rest of this loop.
if (ss_space_in%tDoubles .and. ex_level_to_hf == 0 .and. tDetermHFSpawning) cycle
if (tAllSymSectors) then
ms_parent = return_ms(ilut_parent)
nOccAlpha = (nel + ms_parent) / 2
nOccBeta = (nel - ms_parent) / 2
end if
! If this condition is not met (if all electrons have spin up or all have spin down)
! then there will be no determinants to spawn to, so don't attempt spawning.
if (abs(ms_parent) /= nbasis / 2) then
do ireplica = 1, lenof_sign
call decide_num_to_spawn(parent_sign(ireplica), AvMCExcits, nspawn)
do ispawn = 1, nspawn
! Zero the bit representation, to ensure no extraneous data gets through.
ilut_child = 0_n_int
call generate_excitation(nI_parent, ilut_parent, nI_child, &
ilut_child, exFlag, ic, ex, tParity, prob, &
HElGen, fcimc_excit_gen_store)
! If a valid excitation.
if (.not. IsNullDet(nI_child)) then
call encode_child(ilut_parent, ilut_child, ic, ex)
ilut_child(IlutBits%ind_flag) = 0_n_int
if (tSemiStochastic) then
tChildIsDeterm = is_core_state(ilut_child, nI_child)
! Is the parent state in the core space?
if (tParentIsDeterm) then
! If spawning is both from and to the core space, cancel it.
if (tChildIsDeterm) cycle
call set_flag(ilut_child, flag_determ_parent)
else
if (tChildIsDeterm) call set_flag(ilut_child, flag_deterministic(core_run))
end if
end if
child_sign = attempt_create(nI_parent, ilut_parent, parent_sign, &
nI_child, ilut_child, prob, HElGen, ic, ex, tParity, &
ex_level_to_ref, ireplica, unused_sign, &
AvMCExcits, unused_rdm_real, precond_fac)
else
child_sign = 0.0_dp
end if
! If any (valid) children have been spawned.
if (.not. (all(near_zero(child_sign) .or. ic == 0 .or. ic > 2))) then
call new_child_stats_normal(iter_data_fciqmc, ilut_parent, &
nI_child, ilut_child, ic, ex_level_to_ref, &
child_sign, parent_flags, ireplica)
call create_particle_with_hash_table(nI_child, ilut_child, child_sign, &
ireplica, ilut_parent, iter_data_fciqmc, ierr)
end if ! If a child was spawned.
end do ! Over mulitple particles on same determinant.
end do ! Over the replicas on the same determinant.
end if ! If connected determinants exist to spawn to.
! If this is a core-space determinant then the death step is done in
! determ_projection.
if (.not. tParentIsDeterm) then
call walker_death(iter_data_fciqmc, nI_parent, ilut_parent, parent_hdiag, &
parent_sign, idet, ex_level_to_ref)
end if
end do ! Over all determinants.
if (tSemiStochastic) call determ_projection()
TotWalkersNew = int(TotWalkers)
call end_iter_stats(TotWalkersNew)
call end_iteration_print_warn(TotWalkersNew)
call halt_timer(walker_time)
call set_timer(annihil_time)
call communicate_and_merge_spawns(MaxIndex, iter_data_fciqmc, .false.)
call DirectAnnihilation(TotWalkersNew, MaxIndex, iter_data_fciqmc, ierr)
TotWalkers = int(TotWalkersNew, int64)
call halt_timer(annihil_time)
if (tPrintReplicaOverlaps) call calc_replica_overlaps()
call update_iter_data(iter_data_fciqmc)
if (mod(iter, StepsSft) == 0) then
call set_timer(Stats_Comms_Time)
call calculate_new_shift_wrapper(iter_data_fciqmc, TotParts, tPairedReplicas)
call halt_timer(Stats_Comms_Time)
call ChangeVars(tSingBiasChange, tWritePopsFound)
if (tWritePopsFound) call WriteToPopsfileParOneArr(CurrentDets, TotWalkers)
if (tSingBiasChange) call CalcApproxpDoubles()
if (tSoftExitFound) exit outer_loop
end if
end do ! Over all iterations between Krylov vectors.
end do ! Over all Krylov vectors.
call calc_overlap_matrix(kp%nvecs, krylov_vecs, TotWalkersKP, overlap_matrix)
if (tExactHamil) then
call calc_hamil_exact(kp%nvecs, krylov_vecs, TotWalkersKP, hamil_matrix, krylov_helems)
else
if (tSemiStochastic) then
call calc_projected_hamil(kp%nvecs, krylov_vecs, krylov_vecs_ht, TotWalkersKP, hamil_matrix, &
partial_determ_vecs_kp, full_determ_vecs_kp, krylov_helems)
else
call calc_projected_hamil(kp%nvecs, krylov_vecs, krylov_vecs_ht, TotWalkersKP, hamil_matrix, &
h_diag=krylov_helems)
end if
end if
! Sum the overlap and projected Hamiltonian matrices from the various processors.
call communicate_kp_matrices(overlap_matrix, hamil_matrix)
if (iProcIndex == root) call output_kp_matrices_wrapper(iconfig, overlap_matrices, hamil_matrices)
end do ! Over all repeats for a fixed initial walker configuration.
if (tOverlapPert) call average_and_comm_pert_overlaps(kp%nrepeats)
if (iProcIndex == root) then
call average_kp_matrices_wrapper(iconfig, kp%nrepeats, overlap_matrices, hamil_matrices, &
kp_overlap_mean, kp_hamil_mean, kp_overlap_se, kp_hamil_se)
call find_and_output_lowdin_eigv(iconfig, kp%nvecs, kp_overlap_mean, kp_hamil_mean, nlowdin, lowdin_evals, .true.)
! Calculate data for the testsuite.
s_sum = sum(kp_overlap_mean)
h_sum = sum(kp_hamil_mean)
end if
end do outer_loop ! Over all initial walker configurations.
deallocate(overlap_matrices)
deallocate(hamil_matrices)
deallocate(lowdin_evals)
if (tPopsFile) call WriteToPopsfileParOneArr(CurrentDets, TotWalkers)
if (iProcIndex == root) call write_kpfciqmc_testsuite_data(s_sum, h_sum)
end subroutine perform_kp_fciqmc