subroutine SumEContrib_different_refs(nI, sgn, ilut, dProbFin, tPairedReplicas, ind)
use CalcData, only: qmc_trial_wf
use searching, only: get_con_amp_trial_space
! This is a modified version of SumEContrib for use where the
! projected energies need to be calculated relative to differing
! references.
!
! Some of the arguments to SumeEContrib have been dropped, as they
! only refer to the first of the replicas, and thus cannot be relied on
! n.b. We don't want to just modify SumEContrib for this, as performing
! the calculations _inside_ a sum over runs would radically slow
! down any simulations that were not using differing references
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilut(0:NifTot)
real(dp), intent(in) :: sgn(lenof_sign), dProbFin
logical, intent(in) :: tPairedReplicas
integer, intent(in), optional :: ind
integer :: run, exlevel
HElement_t(dp) :: sgn_run
HElement_t(dp) :: hoffdiag
character(*), parameter :: this_routine = 'SumEContrib_different_refs'
HElement_t(dp) :: amps(size(current_trial_amps, 1))
if (tHistSpawn .or. &
(tCalcFCIMCPsi .and. tFCIMC) .or. tHistEnergies .or. tPrintOrbOcc) then
call stop_all(this_routine, "Not yet supported: Turn off HISTSPAWN,&
& PRINTFCIMCPSI, PRINTORBOCCS, HISTPARTENERGIES, HIST-SPIN-DIST")
end if
! Add in the contributions to the numerator and denominator of the trial
! estimator, if it is being used.
if (tTrialWavefunction .and. present(ind)) then
if (test_flag(ilut, flag_trial)) then
if (ntrial_excits == 1) then
trial_denom = trial_denom + current_trial_amps(1, ind) * sgn
else
if (tPairedReplicas) then
#if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_)
do run = 2, inum_runs, 2
#ifdef CMPLX_
trial_denom(run - 1) = trial_denom(run - 1) + current_trial_amps(run / 2, ind) * &
cmplx(sgn(min_part_type(run - 1)), sgn(max_part_type(run - 1)), dp)
trial_denom(run) = trial_denom(run) + current_trial_amps(run / 2, ind) * &
cmplx(sgn(min_part_type(run)), sgn(max_part_type(run)), dp)
#else
trial_denom(run - 1:run) = trial_denom(run - 1:run) + current_trial_amps(run / 2, ind) * sgn(run - 1:run)
#endif
end do
#else
call stop_all(this_routine, "INVALID")
#endif
else
#ifdef CMPLX_
do run = 1, inum_runs
trial_denom(run) = trial_denom(run) + current_trial_amps(run, ind) * &
cmplx(sgn(min_part_type(run)), sgn(max_part_type(run)), dp)
end do
#else
trial_denom = trial_denom + current_trial_amps(:, ind) * sgn(:)
#endif
end if
end if
if (qmc_trial_wf) then
call get_con_amp_trial_space(ilut, amps)
if (ntrial_excits == 1) then
trial_numerator = trial_numerator + amps(1) * sgn
else
if (tPairedReplicas) then
#if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_)
do run = 2, inum_runs, 2
trial_numerator(run - 1:run) = trial_numerator(run - 1:run) + amps(run / 2) * sgn(run - 1:run)
end do
#else
call stop_all(this_routine, "INVALID")
#endif
else
trial_numerator = trial_numerator + amps * sgn
end if
end if
end if
else if (test_flag(ilut, flag_connected)) then
! Note, only attempt to add in a contribution from the
! connected space if we're not also in the trial space.
if (ntrial_excits == 1) then
trial_numerator = trial_numerator + current_trial_amps(1, ind) * sgn
else
if (tPairedReplicas) then
#if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_)
do run = 2, inum_runs, 2
#ifdef CMPLX_
trial_numerator(run - 1) = trial_numerator(run - 1) + current_trial_amps(run / 2, ind) * &
cmplx(sgn(min_part_type(run - 1)), sgn(max_part_type(run - 1)), dp)
trial_numerator(run) = trial_numerator(run) + current_trial_amps(run / 2, ind) * &
cmplx(sgn(min_part_type(run)), sgn(max_part_type(run)), dp)
#else
trial_numerator(run - 1:run) = trial_numerator(run - 1:run) + current_trial_amps(run / 2, ind) * sgn(run - 1:run)
#endif
end do
#else
call stop_all(this_routine, "INVALID")
#endif
else
#ifdef CMPLX_
do run = 1, inum_runs
trial_numerator(run) = trial_numerator(run) + current_trial_amps(run, ind) * &
cmplx(sgn(min_part_type(run)), sgn(max_part_type(run)), dp)
end do
#else
trial_numerator = trial_numerator + current_trial_amps(:, ind) * sgn
#endif
end if
end if
end if
end if
! This is the normal projected energy calculation, but split over
! multiple runs, rather than done in one go.
do run = 1, inum_runs
! We need to use the excitation level relevant for this run
exlevel = FindBitExcitLevel(ilut, ilutRef(:, run), t_hphf_ic=.true.)
if (tSpinCoupProjE(run) .and. exlevel /= 0) then
if (exlevel <= 2) then
exlevel = 2
else if (FindBitExcitLevel(ilut, ilutRefFlip(:, run)) <= 2) then
exlevel = 2
end if
end if
#ifdef CMPLX_
sgn_run = cmplx(sgn(min_part_type(run)), sgn(max_part_type(run)), dp)
#else
sgn_run = sgn(run)
#endif
hoffdiag = 0.0_dp
if (exlevel == 0 .and. (iter > nEquilSteps)) then
call add_sign_on_run(SumNoatHF)
call add_sign_on_run(NoatHF)
call add_sign_on_run(HFCyc)
call add_sign_on_run(HFOut)
end if
if (tGUGA) then
if (exLevel /= 0) then
! TODO(@Oskar): Perhaps keep csf_i calculated?
hoffdiag = calc_off_diag_guga_ref(ilut, CSF_Info_t(ilut), run, exlevel)
end if
else
if (exlevel == 2 .or. (exlevel == 1 .and. tNoBrillouin)) then
! Obtain the off-diagonal elements
if (tHPHF) then
hoffdiag = hphf_off_diag_helement(ProjEDet(:, run), nI, &
iLutRef(:, run), ilut)
else
if(t_evolve_adjoint(part_type_to_run(run))) then
hoffdiag = get_helement(ProjEDet(:, run), nI, exlevel, &
ilutRef(:, run), ilut)
else
hoffdiag = get_helement(nI, ProjEDet(:, run), exlevel, &
ilut, ilutRef(:, run))
end if
end if
else if (exlevel == 3 .and. &
(t_3_body_excits .or. t_ueg_3_body .or. t_mol_3_body)) then
if(t_evolve_adjoint(part_type_to_run(run))) then
hoffdiag = get_helement(ProjEDet(:, run), nI, exlevel, &
iLutRef(:, run), ilut)
else
hoffdiag = get_helement(nI, ProjEDet(:, run), exlevel, &
ilut, iLutRef(:, run))
end if
end if
end if
if (exlevel == 2) then
NoatDoubs(run) = NoatDoubs(run) + mag_of_run(sgn, run)
end if
! Sum in energy contributions
if (iter > nEquilSteps) then
SumENum(run) = SumENum(run) + enum_contrib()
end if
ENumCyc(run) = ENumCyc(run) + enum_contrib()
ENumOut(run) = ENumOut(run) + enum_contrib()
ENumCycAbs(run) = ENumCycAbs(run) + abs(enum_contrib())
end do
contains
subroutine add_sign_on_run(var)
real(dp), intent(inout) :: var(lenof_sign)
#ifdef CMPLX_
var(min_part_type(run)) = var(min_part_type(run)) + real(sgn_run)
var(max_part_type(run)) = var(max_part_type(run)) + aimag(sgn_run)
#else
var(run) = var(run) + sgn_run
#endif
end subroutine add_sign_on_run
function enum_contrib() result(dE)
HElement_t(dp) :: dE
dE = (hoffdiag * sgn_run) / dProbFin
end function enum_contrib
end subroutine SumEContrib_different_refs