subroutine write_rdm_estimates(rdm_defs, est, final_output, write_to_separate_file, &
tInitsRDM)
! Write RDM estimates to the RDMEstimates file. Specifically, the
! numerator of the energy and spin^2 estimators are output, as is the
! trace (the denominator of these estimator).
! Also, if final_output is true, then output the final total estimates
! to standard output, and close the RDMEstimates unit.
use CalcData, only: tEN2
use FciMCData, only: Iter, PreviousCycles
use LoggingData, only: tRDMInstEnergy, tCalcPropEst, iNumPropToEst
use rdm_data, only: rdm_estimates_t, rdm_definitions_t
use util_mod, only: int_fmt
type(rdm_definitions_t), intent(in) :: rdm_defs
type(rdm_estimates_t), intent(in) :: est
logical, intent(in) :: final_output, write_to_separate_file, tInitsRDM
integer :: irdm, iprop
write(stdout, *) "Writing RDMs to file at iteration ", iter
if (write_to_separate_file) then
if (tRDMInstEnergy) then
write(est%write_unit, '(1x,i13)', advance='no') Iter + PreviousCycles
do irdm = 1, est%nrdms_standard
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%energy_num_inst(irdm)
if (.not. tGUGA) then
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%spin_num_inst(irdm)
end if
if (tEN2) then
write(est%write_unit, '(2(3x,es20.13))', advance='no') &
est%energy_pert_inst(irdm), est%energy_pert_inst(irdm) + est%energy_num_inst(irdm)
end if
if (tCalcPropEst) then
do iprop = 1, iNumPropToEst
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%property_inst(iprop, irdm)
end do
end if
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%norm_inst(irdm)
end do
do irdm = est%nrdms_standard + 1, est%nrdms_standard + est%nrdms_transition
if (tCalcPropEst) then
do iprop = 1, iNumPropToEst
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%property_inst(iprop, irdm)
end do
end if
end do
write(est%write_unit, '()')
else
write(est%write_unit, '(1x,i13)', advance='no') Iter + PreviousCycles
do irdm = 1, est%nrdms_standard
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%energy_num(irdm)
if (.not. tGUGA) then
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%spin_num(irdm)
end if
if (tEN2) then
write(est%write_unit, '(2(3x,es20.13))', advance='no') &
est%energy_pert(irdm), est%energy_pert(irdm) + est%energy_num(irdm)
end if
if (tCalcPropEst) then
do iprop = 1, iNumPropToEst
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%property(iprop, irdm)
end do
end if
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%norm(irdm)
end do
do irdm = est%nrdms_standard + 1, est%nrdms_standard + est%nrdms_transition
if (tCalcPropEst) then
do iprop = 1, iNumPropToEst
write(est%write_unit, '(3x,es20.13)', advance='no') &
est%property(iprop, irdm)
end do
end if
end do
write(est%write_unit, '()')
end if
call neci_flush(est%write_unit)
end if
if (final_output) then
! Banner for the start of the 2-RDM section in output.
if (tInitsRDM) then
write(stdout, '(1x,2("="),1x,"INFORMATION FOR FINAL 2-RDMS (Initiators)",1x,57("="),/)')
else if (tAdaptiveShift) then
write(stdout, '(1x,2("="),1x,"INFORMATION FOR FINAL 2-RDMS (Lagrangian)",1x,57("="),/)')
else
write(stdout, '(1x,2("="),1x,"INFORMATION FOR FINAL 2-RDMS",1x,57("="),/)')
end if
do irdm = 1, est%nrdms_standard
write(stdout, '(1x,"2-RDM ESTIMATES FOR STATE",1x,'//int_fmt(irdm)//',":",)') irdm
write(stdout, '(1x,"Trace of 2-el-RDM before normalisation:",1x,es17.10)') est%trace(irdm)
if (tGUGA) then
write(stdout, '(1x,"Trace of 2-el-RDM after normalisation:",1x,es17.10)') est%trace(irdm) / est%norm(irdm) / 2.0_dp
else
write(stdout, '(1x,"Trace of 2-el-RDM after normalisation:",1x,es17.10)') est%trace(irdm) / est%norm(irdm)
end if
write(stdout, '(1x,"Energy contribution from the 1-RDM:",1x,es17.10)') est%energy_1_num(irdm) / est%norm(irdm)
write(stdout, '(1x,"Energy contribution from the 2-RDM:",1x,es17.10)') est%energy_2_num(irdm) / est%norm(irdm)
write(stdout, '(1x,"*TOTAL ENERGY* CALCULATED USING THE *REDUCED DENSITY MATRICES*:",1x,es20.13,/)') &
est%energy_num(irdm) / est%norm(irdm)
if (tEN2) then
write (stdout, '(1x,"EN2 corrections are below. Note that these may have a much &
&larger error bar than the",/," variational energy above. Please do a &
&blocking analysis rather than just using the energies below.")')
write(stdout, '(1x,"EN2 energy correction:",1x,es17.10)') est%energy_pert(irdm) / est%norm(irdm)
write(stdout, '(1x,"*TOTAL ENERGY* including the EN2 correction:",1x,es17.10,/)') &
(est%energy_num(irdm) + est%energy_pert(irdm)) / est%norm(irdm)
end if
! Hermiticity error measures.
write(stdout, '(1x,"Hermiticty error estimates:")')
write(stdout, '(1x,i15,f30.20,5x,a41)') Iter + PreviousCycles, est%max_error_herm(irdm), &
'(Iteration, MAX ABS ERROR IN HERMITICITY)'
write(stdout, '(1x,i15,f30.20,5x,a41,/)') Iter + PreviousCycles, est%sum_error_herm(irdm), &
'(Iteration, SUM ABS ERROR IN HERMITICITY)'
end do
do irdm = est%nrdms_standard + 1, est%nrdms
associate(state_labels => rdm_defs%state_labels, repeat_label => rdm_defs%repeat_label)
write(stdout, '(1x,"2-RDM ESTIMATES FOR TRANSITION",1x,'//int_fmt(state_labels(2, irdm))//'," -> ",&
&'//int_fmt(state_labels(1, irdm))//',1x,"(",i1,")",)') &
state_labels(2, irdm), state_labels(1, irdm), repeat_label(irdm)
end associate
write(stdout, '(1x,"Trace of 2-el-RDM before normalisation:",1x,es17.10)') est%trace(irdm)
write(stdout, '(1x,"Trace of 2-el-RDM after normalisation:",1x,es17.10,/)') est%trace(irdm) / est%norm(irdm)
! Hermiticity difference measures - these shouldn't be zero for
! transition RDMs, but it useful to give them to the test
! suite, to make sure somebody doesn't change something to
! start adding an RDM element on the wrong side of the diagonal.
write(stdout, '(1x,"Hermiticty difference estimates, for test suite:")')
write(stdout, '(1x,i15,f30.20,5x,a41)') Iter + PreviousCycles, est%max_error_herm(irdm), &
'(Iteration, MAX ABS DIFF IN HERMITICITY)'
write(stdout, '(1x,i15,f30.20,5x,a41,/)') Iter + PreviousCycles, est%sum_error_herm(irdm), &
'(Iteration, SUM ABS DIFF IN HERMITICITY)'
end do
! Banner for the end of the 2-RDM section in output.
write(stdout, '(1x,89("="))')
end if
end subroutine write_rdm_estimates