subroutine calc_ests_simple_initiator(ValidSpawned, proj_energy)
! This routine calculates the various energy estimates to be printed
! to the file for the preconditioned approach.
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(len=*), parameter :: t_r = 'calc_ests_simple_initiator'
#if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_)
integer :: i, j, PartInd, DetHash, determ_pos, nrepeats, ierr
integer :: nI_spawn(nel)
real(dp) :: spwnsign_init(lenof_sign), spwnsign_non(lenof_sign)
real(dp) :: spwnsign(lenof_sign), cursign(lenof_sign)
real(dp) :: hdiag, mean_energy(replica_est_len)
logical :: tCoreDet, tSuccess, abort(lenof_sign)
real(dp) :: b_term(replica_est_len), b_term_all(replica_est_len)
real(dp) :: c_term(replica_est_len), c_term_all(replica_est_len)
! 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
en2_pert = 0.0_dp
en2_pert_all = 0.0_dp
e_squared_num = 0.0_dp
e_squared_num_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
i = 1
do
call decode_bit_det(nI_spawn, SpawnedParts(:, i))
! How many times is this particular determinant repeated in the
! spawning array (either 0 or 1)?
nrepeats = 0
spwnsign_init = 0.0_dp
spwnsign_non = 0.0_dp
! Only need to check initiator status for one replica - if
! one is an initiator then they all are, when using the simple
! initiator approximation.
if (test_flag(SpawnedParts(:, i), get_initiator_flag(1))) then
! This is an initiator
call extract_sign(SpawnedParts(:, i), spwnsign_init)
if (i + 1 <= ValidSpawned) then
if (DetBitEq(SpawnedParts(:, i), SpawnedParts(:, i + 1), nifd)) then
! This next state is a different state, and so will
! be a non-initiator
call extract_sign(SpawnedParts(:, i + 1), spwnsign_non)
nrepeats = 1
end if
end if
else
call extract_sign(SpawnedParts(:, i), spwnsign_non)
end if
spwnsign = spwnsign_init + spwnsign_non
! 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)
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
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) then
en2_pert(j) = en2_pert(j) &
+ spwnsign_non(2 * j - 1) * spwnsign_non(2 * j) &
/ ((tau**2) * (mean_energy(j) - hdiag))
end if
end do
i = i + 1 + nrepeats
if (i > ValidSpawned + 1) call stop_all(t_r, "The spawning array index is larger than it should be.")
if (i == ValidSpawned + 1) exit
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
b_term = 0.0_dp
b_term_all = 0.0_dp
c_term = 0.0_dp
c_term_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)
b_term(j) = b_term(j) + cursign(2 * j - 1) * cursign(2 * j)
c_term(j) = c_term(j) - cursign(2 * j - 1) * cursign(2 * j) / (mean_energy(j) - hdiag)
end do
end do
i = 1
do
call decode_bit_det(nI_spawn, SpawnedParts(:, i))
! How many times is this particular determinant repeated in the
! spawning array (either 0 or 1)?
nrepeats = 0
spwnsign_init = 0.0_dp
spwnsign_non = 0.0_dp
! Only need to check initiator status for one replica - if
! one is an initiator then they all are, when using the simple
! initiator approximation.
if (test_flag(SpawnedParts(:, i), get_initiator_flag(1))) then
! This is an initiator
call extract_sign(SpawnedParts(:, i), spwnsign_init)
if (i + 1 <= ValidSpawned) then
if (DetBitEq(SpawnedParts(:, i), SpawnedParts(:, i + 1), nifd)) then
! This next state is a different state, and so will
! be a non-initiator
call extract_sign(SpawnedParts(:, i + 1), spwnsign_non)
nrepeats = 1
end if
end if
else
call extract_sign(SpawnedParts(:, i), spwnsign_non)
end if
spwnsign = spwnsign_init + spwnsign_non
! 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)
hdiag = det_diagH(PartInd) + Hii
if (tSemiStochastic) then
tCoreDet = check_determ_flag(CurrentDets(:, PartInd))
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
b_term(j) = b_term(j) &
+ (spwnsign(2 * j - 1) * cursign(2 * j) &
+ spwnsign(2 * j) * cursign(2 * j - 1)) &
/ (2.0_dp * tau * (mean_energy(j) - hdiag))
end do
end if
do j = 1, replica_est_len
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 do
i = i + 1 + nrepeats
if (i > ValidSpawned + 1) call stop_all(t_r, "The spawning array index is larger than it should be.")
if (i == ValidSpawned + 1) exit
end do
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
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))
b_term(j) = b_term(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 * (mean_energy(j) - rep%core_ham_diag(i) - Hii))
end do
end if
end do
end if
call MPISumAll(en2_new, en2_new_all)
call MPISumAll(b_term, b_term_all)
call MPISumAll(c_term, c_term_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)
if (tEN2Rigorous) then
write(replica_est_unit, '(2(3x,es20.13))', advance='no') b_term_all(j), c_term_all(j)
end if
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(t_r, "Should not be here")
#endif
end subroutine calc_ests_simple_initiator