subroutine SumEContrib(nI, ExcitLevel, RealWSign, ilut, HDiagCurr, &
HOffDiagCurr, dProbFin, tPairedReplicas, ind)
use CalcData, only: qmc_trial_wf
use searching, only: get_con_amp_trial_space
use guga_excitations, only: calc_off_diag_guga_ref_direct
integer, intent(in) :: nI(nel), ExcitLevel
real(dp), intent(in) :: RealwSign(lenof_sign)
integer(n_int), intent(in) :: ilut(0:NIfTot)
real(dp), intent(in) :: HDiagCurr, dProbFin
HElement_t(dp), intent(in) :: HOffDiagCurr
logical, intent(in) :: tPairedReplicas
integer, intent(in), optional :: ind
integer :: i, ExcitLevel_local, ExcitLevelSpinCoup
integer :: run, tmp_exlevel, step
HElement_t(dp) :: HOffDiag(inum_runs), tmp_off_diag(inum_runs), tmp_diff(inum_runs)
character(*), parameter :: this_routine = 'SumEContrib'
#ifdef CMPLX_
complex(dp) :: CmplxwSign
#endif
real(dp), allocatable :: amps(:)
real(dp) :: w(0:2)
if (allocated(current_trial_amps)) allocate (amps(size(current_trial_amps, 1)))
if (tReplicaReferencesDiffer) then
call SumEContrib_different_refs(nI, realWSign, ilut, dProbFin, tPairedReplicas, ind)
return
end if
HOffDiag = 0.0_dp
! Add in the contributions to the numerator and denominator of the trial
! estimator, if it is being used.
#ifdef CMPLX_
CmplxwSign = ARR_RE_OR_CPLX(realwsign, 1)
if (tTrialWavefunction .and. present(ind)) then
if (test_flag(ilut, flag_trial)) then
if (ntrial_excits == 1) then
trial_denom = trial_denom + conjg(current_trial_amps(1, ind)) * CmplxwSign
! this does somehow not support kmneci
if (test_flag(ilut, get_initiator_flag_by_run(1))) &
init_trial_denom = init_trial_denom + conjg( &
current_trial_amps(1, ind)) * CmplxwSign
else if (ntrial_excits == lenof_sign) then
call stop_all(this_routine, 'ntrial_excits has to be 1 currently for complex')
end if
if (qmc_trial_wf) then
call stop_all(this_routine, 'qmc_trial_wf currently not implemented for complex')
end if
else if (test_flag(ilut, flag_connected)) then
if (ntrial_excits == 1) then
trial_numerator = trial_numerator + conjg(current_trial_amps(1, ind)) * cmplxwsign
if (test_flag(ilut, get_initiator_flag_by_run(1))) &
init_trial_numerator = init_trial_numerator + conjg( &
current_trial_amps(1, ind)) * CmplxwSign
else if (ntrial_excits == lenof_sign) then
call stop_all(this_routine, 'ntrial_excits has to be 1 currently for complex')
end if
end if
end if
#else
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) * RealwSign
trial_denom_inst = trial_denom_inst + current_trial_amps(1, ind) * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
init_trial_denom(run) = init_trial_denom(run) + &
current_trial_amps(1, ind) * RealwSign(run)
end do
else if (ntrial_excits == lenof_sign) then
trial_denom = trial_denom + current_trial_amps(:, ind) * RealwSign
trial_denom_inst = trial_denom_inst + current_trial_amps(:, ind) * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
init_trial_denom(run) = init_trial_denom(run) + &
current_trial_amps(run, ind) * RealwSign(run)
end do
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) * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
init_trial_numerator(run) = init_trial_numerator(run) + &
amps(1) * RealwSign(run)
end do
else if (ntrial_excits == lenof_sign) then
trial_numerator = trial_numerator + amps * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
init_trial_numerator(run) = init_trial_numerator(run) + &
amps(run) * RealwSign(run)
end do
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) * RealwSign
trial_num_inst = trial_num_inst + current_trial_amps(1, ind) * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
! this is the real case, so inum_runs == lenof_sign
init_trial_numerator(run) = init_trial_numerator(run) + &
current_trial_amps(1, ind) * RealwSign(run)
end do
else if (ntrial_excits == lenof_sign) then
trial_numerator = trial_numerator + current_trial_amps(:, ind) * RealwSign
trial_num_inst = trial_num_inst + current_trial_amps(:, ind) * RealwSign
do run = 1, inum_runs
if (test_flag(ilut, get_initiator_flag_by_run(run))) &
init_trial_numerator(run) = init_trial_numerator(run) + &
current_trial_amps(run, ind) * RealwSign(run)
end do
end if
end if
end if
#endif
! ExcitLevel indicates the excitation level between the det and
! *one* of the determinants in an HPHF/MomInv function. If needed,
! calculate the connection between it and the other one. If either
! is connected, then it has to be counted. Since the excitation
! level is the same to either det, we don't need to consider the
! spin-coupled det of both reference and current HPHFs.
!
! For determinants, set ExcitLevel_local == ExcitLevel.
ExcitLevel_local = ExcitLevel
if (tSpinCoupProjE(1) .and. (ExcitLevel /= 0)) then
ExcitLevelSpinCoup = FindBitExcitLevel(iLutRefFlip(:, 1), &
ilut, 2)
if (ExcitLevelSpinCoup <= 2 .or. ExcitLevel <= 2) &
ExcitLevel_local = 2
end if
! this is the first important change to make the triples run!!
! consider the matrix elements of triples!
! for the real-time i have to distinguish between the first and
! second RK step, if i want to keep track of the statistics
! seperately: in the first loop i analyze the the wavefunction
! from on step behind.. so store it in the "normal" noathf var
if (ExcitLevel_local == 0) then
if (.not. t_real_time_fciqmc .or. runge_kutta_step == 2) then
HFCyc(1:lenof_sign) = HFCyc(1:lenof_sign) + RealwSign
HFOut(1:lenof_sign) = HFOut(1:lenof_sign) + RealwSign
NoatHF(1:lenof_sign) = NoatHF(1:lenof_sign) + RealwSign
if (iter > NEquilSteps) then
SumNoatHF(1:lenof_sign) = SumNoatHF(1:lenof_sign) + RealwSign
end if
end if
end if
! Perform normal projection onto reference determinant
if (t_adjoint_replicas) then
if ( ExcitLevel_local == 2 &
.or. (ExcitLevel_local == 1 .and. tNoBrillouin) &
.or. (ExcitLevel_local == 3 &
.and. (t_3_body_excits .or. t_ueg_3_body .or. t_mol_3_body))) then
! Obtain off-diagonal element
do run = 1, inum_runs
if(t_evolve_adjoint(part_type_to_run(run))) then
HOffDiag(run) = &
get_helement(ProjEDet(:, 1), nI, ExcitLevel, ilutRef(:, 1), ilut)
else
HOffDiag(run) = &
get_helement(nI, ProjEDet(:,1), ExcitLevel, ilut, ilutRef(:, 1))
end if
end do
end if
else
HOffDiag(1:inum_runs) = HOffDiagCurr
endif
! For the real-space Hubbard model, determinants are only
! connected to excitations one level away, and Brillouins
! theorem cannot hold.
!
! For Rotated orbitals, Brillouins theorem also cannot hold,
! and energy contributions from walkers on singly excited
! determinants must also be included in the energy values
! along with the doubles
! RT_M_Merge: Adjusted to kmneci
! rmneci_setup: Added multirun functionality for real-time
if (ExcitLevel_local == 2) then
do run = 1, inum_runs
if (.not. t_real_time_fciqmc .or. runge_kutta_step == 2) then
#if defined(CMPLX_)
NoatDoubs(run) = NoatDoubs(run) + sum(abs(RealwSign &
(min_part_type(run):max_part_type(run))))
#else
NoatDoubs(run) = NoatDoubs(run) + abs(RealwSign(run))
#endif
end if
end do
end if
! L_{0,1,2} norms of walker weights by excitation level.
if (tLogEXLEVELStats) then
do run = 1, inum_runs
#ifdef CMPLX_
w(0) = real(1 + max_part_type(run) - min_part_type(run), dp)
w(1) = sum(abs(RealwSign(min_part_type(run): &
max_part_type(run))))
w(2) = sum(RealwSign(min_part_type(run): &
max_part_type(run))**2)
#else
w(0) = 1_dp
w(1) = abs(RealwSign(run))
w(2) = RealwSign(run)**2
#endif
EXLEVEL_WNorm(0:2, ExcitLevel_local, run) = &
EXLEVEL_WNorm(0:2, ExcitLevel_local, run) + w(0:2)
end do ! run
end if ! tLogEXLEVELStats
! if in the real-time fciqmc: when we are in the 2nd loop
! return here since, the energy got already calculated in the
! first RK step, and doing it on the intermediate step would
! be meaningless
! Sum in energy contribution
do run = 1, inum_runs
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())
if (test_flag(ilut, get_initiator_flag_by_run(run))) then
InitsENumCyc(run) = InitsENumCyc(run) + enum_contrib()
end if
end do
! -----------------------------------
! HISTOGRAMMING
! -----------------------------------
if ((tHistSpawn .or. (tCalcFCIMCPsi .and. tFCIMC)) .and. &
(iter >= NHistEquilSteps)) then
! Histogram particles by determinant
call add_hist_spawn(ilut, RealwSign, ExcitLevel_local, dProbFin)
else if (tHistEnergies) then
! Histogram particles by energy
call add_hist_energies(ilut, RealwSign, HDiagCurr)
end if
! Maintain a list of the degree of occupation of each orbital
if (tPrintOrbOcc .and. (iter >= StartPrintOrbOcc)) then
if (iter == StartPrintOrbOcc .and. &
DetBitEq(ilut, ilutHF, nifd)) then
write(stdout, *) 'Beginning to fill the HF orbital occupation list &
&during iteration', iter
if (tPrintOrbOccInit) &
write(stdout, *) 'Only doing so for initiator determinants'
end if
if ((tPrintOrbOccInit .and. test_flag(ilut, get_initiator_flag(1))) &
.or. .not. tPrintOrbOccInit) then
forall (i=1:nel) OrbOccs(nI(i)) &
= OrbOccs(nI(i)) + (RealwSign(1) * RealwSign(1))
end if
end if
contains
function enum_contrib() result(dE)
HElement_t(dp) :: dE
dE = (HOffDiag(run) * ARR_RE_OR_CPLX(RealwSign, run)) / dProbFin
end function enum_contrib
end subroutine SumEContrib