subroutine calc_ests_and_set_init_flags(ValidSpawned, proj_energy)
! This routine calculates the various energy estimates to be printed
! to the file for the preconditioned approach.
! *IMPORTANT* It also sets the initiator flag for all spawnings to
! already-occupied determinants. This is because if tau=1, then
! all walkers will die in the death step to follow, which causes
! the initiator criteria to not be properly applied (as the
! AnnihilateSpawnedParts routine thinks these are spawnings to
! unoccupied determinants). So setting the initiator flag here
! allows everything to proceed as normal.
! This second point means that this routine *does* affect the
! walker dynamics - it doesn't just calculate some energies
! without affecting walkers. This step is done here for efficiency.
use CalcData, only: tEN2Init, tEN2Rigorous, tTruncInitiator
use global_det_data, only: det_diagH, replica_est_len
use semi_stoch_procs, only: core_space_pos, check_determ_flag
integer, intent(in) :: ValidSpawned
real(dp), intent(in) :: proj_energy(lenof_sign)
character(*), parameter :: this_routine = 'calc_ests_and_set_init_flags'
#if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_)
integer :: i, j, PartInd, DetHash, determ_pos, ierr
integer :: nI_spawn(nel)
real(dp) :: spwnsign(lenof_sign), cursign(lenof_sign)
real(dp) :: hdiag, mean_energy(replica_est_len)
logical :: tCoreDet, tSuccess, abort(lenof_sign)
! Allow room to send up to 1000 elements.
real(dp) :: send_arr(6 * replica_est_len)
! Allow room to receive up to 1000 elements.
real(dp) :: recv_arr(6 * replica_est_len)
var_e_num = 0.0_dp
rep_est_overlap = 0.0_dp
var_e_num_all = 0.0_dp
rep_est_overlap_all(1) = 0.0_dp
e_squared_num = 0.0_dp
e_squared_num_all = 0.0_dp
en2_pert = 0.0_dp
en2_pert_all = 0.0_dp
precond_e_num = 0.0_dp
precond_denom = 0.0_dp
precond_e_num_all = 0.0_dp
precond_denom_all = 0.0_dp
do i = 1, replica_est_len
mean_energy(i) = Hii + (proj_energy(2 * i - 1) &
+ proje_ref_energy_offsets(2 * i - 1) + proj_energy(2 * i) &
+ proje_ref_energy_offsets(2 * i)) / 2.0_dp
end do
tCoreDet = .false.
tSuccess = .false.
abort = .false.
tDetermSpawnedTo = .false.
associate(rep => cs_replicas(core_run))
! Contributions from diagonal where necessary
do i = 1, int(TotWalkers)
hdiag = det_diagH(i) + Hii
call extract_sign(CurrentDets(:, i), cursign)
do j = 1, replica_est_len
rep_est_overlap(j) = rep_est_overlap(j) &
+ cursign(2 * j - 1) * cursign(2 * j)
var_e_num(j) = var_e_num(j) &
+ hdiag * cursign(2 * j - 1) * cursign(2 * j)
e_squared_num(j) = e_squared_num(j) &
+ (hdiag**2) * cursign(2 * j - 1) * cursign(2 * j)
end do
end do
do i = 1, ValidSpawned
call decode_bit_det(nI_spawn, SpawnedParts(:, i))
call extract_sign(SpawnedParts(:, i), spwnsign)
! Now add in the diagonal elements
call hash_table_lookup(nI_spawn, SpawnedParts(:, i), nifd, HashIndex, &
CurrentDets, PartInd, DetHash, tSuccess)
if (tSuccess) then
call extract_sign(CurrentDets(:, PartInd), cursign)
! Set initiator flags for the spawning, before the currently
! occupied determinant is potentially killed in the death step.
if (tTruncInitiator) then
do j = 1, lenof_sign
if (abs(cursign(j)) > 1.e-12_dp) then
call set_flag(SpawnedParts(:, i), get_initiator_flag(j))
end if
end do
end if
hdiag = det_diagH(PartInd) + Hii
if (tSemiStochastic) then
tCoreDet = check_determ_flag(CurrentDets(:, PartInd), core_run)
if (tCoreDet) then
determ_pos = core_space_pos(SpawnedParts(:, i), nI_spawn) - rep%determ_displs(iProcIndex)
tDetermSpawnedTo(determ_pos) = .true.
spwnsign = spwnsign + rep%partial_determ_vecs(:, determ_pos)
end if
end if
do j = 1, replica_est_len
var_e_num(j) = var_e_num(j) &
- (spwnsign(2 * j - 1) * cursign(2 * j) &
+ spwnsign(2 * j) * cursign(2 * j - 1)) / (2.0_dp * tau)
precond_e_num(j) = precond_e_num(j) &
- (spwnsign(2 * j - 1) * hdiag * cursign(2 * j) &
+ spwnsign(2 * j) * hdiag * cursign(2 * j - 1)) / &
(2.0_dp * (mean_energy(j) - hdiag))
precond_denom(j) = precond_denom(j) &
- (spwnsign(2 * j - 1) * cursign(2 * j) &
+ spwnsign(2 * j) * cursign(2 * j - 1)) &
/ (2.0_dp * (mean_energy(j) - hdiag))
e_squared_num(j) = e_squared_num(j) &
- (spwnsign(2 * j - 1) * hdiag * cursign(2 * j) &
+ spwnsign(2 * j) * hdiag * cursign(2 * j - 1)) / tau
end do
end if
! Add in the contributions corresponding to off-diagonal
! elements of the Hamiltonian
do j = 1, replica_est_len
if (abs(spwnsign(2 * j - 1)) > 1.e-12_dp .and. abs(spwnsign(2 * j)) > 1.e-12_dp) then
hdiag = real(extract_spawn_hdiag(SpawnedParts(:, i)), dp)
precond_e_num(j) = precond_e_num(j) &
+ spwnsign(2 * j - 1) * spwnsign(2 * j) &
/ (tau * (mean_energy(j) - hdiag))
e_squared_num(j) = e_squared_num(j) &
+ spwnsign(2 * j - 1) * spwnsign(2 * j) / (tau**2)
! Only get EN2 contribution if we're due to cancel this
! spawning on both replicas
if (tEN2Init .and. (.not. tSuccess)) then
abort(2 * j - 1) = test_abort_spawn(SpawnedParts(:, i), 2 * j - 1)
abort(2 * j) = test_abort_spawn(SpawnedParts(:, i), 2 * j)
if (abort(2 * j - 1) .and. abort(2 * j)) then
en2_pert(j) = en2_pert(j) &
+ spwnsign(2 * j - 1) * spwnsign(2 * j) &
/ ((tau**2) * (mean_energy(j) - hdiag))
end if
end if
end if
end do
end do
! Contribution from deterministic states that were not spawned to
if (tSemiStochastic) then
do i = 1, rep%determ_sizes(iProcIndex)
if (.not. tDetermSpawnedTo(i)) then
call extract_sign(CurrentDets(:, rep%indices_of_determ_states(i)), cursign)
do j = 1, replica_est_len
! Variational energy terms
var_e_num(j) = var_e_num(j) &
- (rep%partial_determ_vecs(2 * j - 1, i) * cursign(2 * j) &
+ rep%partial_determ_vecs(2 * j, i) &
* cursign(2 * j - 1)) / (2.0_dp * tau)
! Preconditioned estimator
precond_e_num(j) = precond_e_num(j) &
+ rep%partial_determ_vecs(2 * j - 1, i) &
* rep%partial_determ_vecs(2 * j, i) &
/ (tau * (mean_energy(j) - rep%core_ham_diag(i) - Hii))
precond_e_num(j) = precond_e_num(j) &
- (rep%partial_determ_vecs(2 * j - 1, i) &
* (rep%core_ham_diag(i) + Hii) * cursign(2 * j) &
+ rep%partial_determ_vecs(2 * j, i) &
* (rep%core_ham_diag(i) + Hii) * cursign(2 * j - 1)) &
/ (2.0_dp * (mean_energy(j) - rep%core_ham_diag(i) - Hii))
precond_denom(j) = precond_denom(j) &
- (rep%partial_determ_vecs(2 * j - 1, i) &
* cursign(2 * j) + rep%partial_determ_vecs(2 * j, i) &
* cursign(2 * j - 1)) / (2.0_dp * (mean_energy(j) &
- rep%core_ham_diag(i) - Hii))
! E squared terms
e_squared_num(j) = e_squared_num(j) &
+ rep%partial_determ_vecs(2 * j - 1, i) &
* rep%partial_determ_vecs(2 * j, i) / (tau**2)
e_squared_num(j) = e_squared_num(j) &
- (rep%partial_determ_vecs(2 * j - 1, i) &
* (rep%core_ham_diag(i) + Hii) * cursign(2 * j) &
+ rep%partial_determ_vecs(2 * j, i) &
* (rep%core_ham_diag(i) + Hii) * cursign(2 * j - 1)) / tau
end do
end if
end do
end if
! ---- MPI communication --------------------------------
send_arr(0 * replica_est_len + 1:1 * replica_est_len) = var_e_num
send_arr(1 * replica_est_len + 1:2 * replica_est_len) = rep_est_overlap
send_arr(2 * replica_est_len + 1:3 * replica_est_len) = en2_pert
send_arr(3 * replica_est_len + 1:4 * replica_est_len) = precond_e_num
send_arr(4 * replica_est_len + 1:5 * replica_est_len) = precond_denom
send_arr(5 * replica_est_len + 1:6 * replica_est_len) = e_squared_num
call MPIBarrier(ierr)
call MPISumAll(send_arr, recv_arr)
var_e_num_all = recv_arr(0 * replica_est_len + 1:1 * replica_est_len)
rep_est_overlap_all = recv_arr(1 * replica_est_len + 1:2 * replica_est_len)
en2_pert_all = recv_arr(2 * replica_est_len + 1:3 * replica_est_len)
precond_e_num_all = recv_arr(3 * replica_est_len + 1:4 * replica_est_len)
precond_denom_all = recv_arr(4 * replica_est_len + 1:5 * replica_est_len)
e_squared_num_all = recv_arr(5 * replica_est_len + 1:6 * replica_est_len)
! -------------------------------------------------------
! Use a slightly more rigorous expression for the variational plus
! EN2 energy.
if (tEN2Rigorous) then
en2_new = 0.0_dp
en2_new_all = 0.0_dp
mean_energy = var_e_num_all / rep_est_overlap_all
do i = 1, int(TotWalkers)
hdiag = det_diagH(i) + Hii
call extract_sign(CurrentDets(:, i), cursign)
do j = 1, replica_est_len
en2_new(j) = en2_new(j) &
+ hdiag * cursign(2 * j - 1) * cursign(2 * j)
end do
end do
do i = 1, ValidSpawned
call decode_bit_det(nI_spawn, SpawnedParts(:, i))
call extract_sign(SpawnedParts(:, i), spwnsign)
! Now add in the diagonal elements
call hash_table_lookup(nI_spawn, SpawnedParts(:, i), nifd, HashIndex, &
CurrentDets, PartInd, DetHash, tSuccess)
if (tSuccess .and. tSemiStochastic) then
tCoreDet = check_determ_flag(CurrentDets(:, PartInd), core_run)
if (tCoreDet) then
determ_pos = core_space_pos(SpawnedParts(:, i), nI_spawn) - rep%determ_displs(iProcIndex)
tDetermSpawnedTo(determ_pos) = .true.
spwnsign = spwnsign + rep%partial_determ_vecs(:, determ_pos)
end if
end if
do j = 1, replica_est_len
if (abs(spwnsign(2 * j - 1)) > 1.e-12_dp .and. abs(spwnsign(2 * j)) > 1.e-12_dp) then
hdiag = real(extract_spawn_hdiag(SpawnedParts(:, i)), dp)
en2_new(j) = en2_new(j) &
+ spwnsign(2 * j - 1) * spwnsign(2 * j) &
/ ((tau**2) * (mean_energy(j) - hdiag))
end if
end do
end do
if (tSemiStochastic) then
do i = 1, rep%determ_sizes(iProcIndex)
if (.not. tDetermSpawnedTo(i)) then
do j = 1, replica_est_len
en2_new(j) = en2_new(j) &
+ rep%partial_determ_vecs(2 * j - 1, i) &
* rep%partial_determ_vecs(2 * j, i) &
/ ((tau)**2 * (mean_energy(j) &
- rep%core_ham_diag(i) - Hii))
end do
end if
end do
end if
call MPISumAll(en2_new, en2_new_all)
end if
if (iProcIndex == 0) then
write(replica_est_unit, '(1x,i13)', advance='no') Iter + PreviousCycles
do j = 1, replica_est_len
write(replica_est_unit, '(2(3x,es20.13))', advance='no') var_e_num_all(j), e_squared_num_all(j)
if (tEN2Init) then
write(replica_est_unit, '(2(3x,es20.13))', advance='no') en2_pert_all(j), var_e_num_all(j) + en2_pert_all(j)
end if
if (tEN2Rigorous) then
write(replica_est_unit, '(1(3x,es20.13))', advance='no') en2_new_all(j)
end if
write(replica_est_unit, '(1(3x,es20.13))', advance='no') rep_est_overlap_all(j)
write(replica_est_unit, '(2(3x,es20.13))', advance='no') precond_e_num_all(j), precond_denom_all(j)
end do
write(replica_est_unit, '()')
end if
end associate
#else
unused_var(ValidSpawned); unused_var(proj_energy)
call stop_all(this_routine, "Should not be here")
#endif
end subroutine calc_ests_and_set_init_flags