subroutine Standalone_Errors()
use sym_mod, only: getsym
USE MolproPlugin, only: MolproPluginResult
real(dp) :: mean_ProjE_re, mean_ProjE_im, mean_Shift
real(dp) :: ProjE_Err_re, ProjE_Err_im, Shift_Err
logical :: tNoProjEValue, tNoShiftValue
integer :: iroot, isymh
TYPE(BasisFn) RefSym
HElement_t(dp) :: h_tmp
real(dp) :: Hii, BestEnergy, EnergyDiff
!Automatic error analysis
call error_analysis(tSinglePartPhase(1), iBlockingIter(1), mean_ProjE_re, ProjE_Err_re, &
mean_ProjE_im, ProjE_Err_im, mean_Shift, Shift_Err, tNoProjEValue, tNoShiftValue, &
equilshift=iBlockEquilShift, equilproje=iBlockEquilProjE)
call MPIBCast(ProjectionE)
call MPIBCast(mean_ProjE_re)
call MPIBCast(ProjE_Err_re)
call MPIBCast(mean_ProjE_im)
call MPIBCast(ProjE_Err_im)
call MPIBCast(mean_Shift)
call MPIBCast(Shift_Err)
call MPIBCast(tNoProjEValue)
call MPIBCast(tNoShiftValue)
h_tmp = get_helement(FDet, FDet, 0)
Hii = real(h_tmp, dp)
iroot = 1
CALL GetSym(FDet, NEl, G1, NBasisMax, RefSym)
isymh = int(RefSym%Sym%S) + 1
write(stdout, 10101) iroot, isymh
10101 format(//'RESULTS FOR STATE', i2, '.', i1/'====================='/)
write(stdout, '('' Current reference energy'',T52,F19.12)') Hii
if (tNoProjEValue) then
write(stdout, '('' No projected energy value could be obtained'',T52)')
else
write(stdout, '('' Projected correlation energy'',T52,F19.12)') mean_ProjE_re
write(stdout, '('' Estimated error in Projected correlation energy'',T52,F19.12)') ProjE_Err_re
end if
if (lenof_sign == 2 .and. inum_runs == 1) then
write(stdout, '('' Projected imaginary energy'',T52,F19.12)') mean_ProjE_im
write(stdout, '('' Estimated error in Projected imaginary energy'',T52,F19.12)') ProjE_Err_im
end if
if (tNoShiftValue) then
write(stdout, '('' No shift energy value could be obtained'',T52)')
else
write(stdout, '('' Shift correlation energy'',T52,F19.12)') mean_Shift
write(stdout, '('' Estimated error in shift correlation energy'',T52,F19.12)') shift_err
end if
!Do shift and projected energy agree?
write(stdout, "(A)")
if (tNoProjEValue .and. tNoShiftValue) return
EnergyDiff = abs(mean_Shift - mean_ProjE_re)
if (EnergyDiff <= sqrt(shift_err**2 + ProjE_Err_re**2)) then
write(stdout, "(A,F15.8)") " Projected and shift energy estimates agree " &
& //"within errorbars: EDiff = ", EnergyDiff
else if (EnergyDiff <= sqrt((max(shift_err, ProjE_Err_re) * 2)**2 + min(shift_err, ProjE_Err_re)**2)) then
write(stdout, "(A,F15.8)") " Projected and shift energy estimates agree to within two sigma " &
& //"of largest error: EDiff = ", EnergyDiff
else
write(stdout, "(A,F15.8)") " Projected and shift energy estimates do not agree to " &
& //"within approximate errorbars: EDiff = ", EnergyDiff
end if
if (ProjE_Err_re < shift_err) then
BestEnergy = mean_ProjE_re + Hii
else
BestEnergy = mean_shift + Hii
end if
write(stdout, "(A)")
write(stdout, "(A,F20.8,A,G15.6)") " Total projected energy ", &
mean_ProjE_re + Hii, " +/- ", ProjE_Err_re
write(stdout, "(A,F20.8,A,G15.6)") " Total shift energy ", &
mean_shift + Hii, " +/- ", shift_err
CALL MolproPluginResult('ENERGY', [BestEnergy])
CALL MolproPluginResult('FCIQMC_ERR', [min(ProjE_Err_re, shift_err)])
write(stdout, "(/)")
end subroutine Standalone_Errors