subroutine print_1rdms_from_sf2rdms_wrapper(rdm_defs, one_rdms, two_rdms)
! Wrapper routine to calculate and print spinfree 1-RDMs from
! spinfree 2-RDMs, as read in directly from a spinfree_TwoRDM file
! which was output by NECI. This routine will *not* work on standard
! format NECI RDM objects, nor is it capable of outputting spinned
! 1-RDMs in open shell systems, only spinfree 1-RDMs.
use Parallel_neci, only: MPISumAll
use rdm_data, only: rdm_definitions_t, one_rdm_t, rdm_list_t
use rdm_finalising, only: calc_1rdms_from_spinfree_2rdms, write_1rdm
type(rdm_definitions_t), intent(in) :: rdm_defs
type(one_rdm_t), intent(inout) :: one_rdms(:)
type(rdm_list_t), intent(in) :: two_rdms
integer :: irdm
real(dp) :: rdm_norm_all(two_rdms%sign_length), norm_1rdm
character(len=*), parameter :: t_r = 'print_1rdms_from_sf2rdms_wrapper'
if (size(one_rdms) == 0) then
call stop_all(t_r, "You have asked to print 1-RDMs but they are not allocated. Make sure &
&that you have asked for both 1-RDMs and 2-RDMs to be calculated by &
&seting the RDMExcitLevel to 3, i.e. 'CALCRDMONFLY 3 ...' in input options.")
end if
! The read-in spinfree 2-RDMs should already be normalised.
rdm_norm_all = 1.0_dp
call calc_1rdms_from_spinfree_2rdms(one_rdms, two_rdms, rdm_norm_all)
do irdm = 1, size(one_rdms)
call write_1rdm(rdm_defs, one_rdms(irdm)%matrix, irdm, 1.0_dp, .true.)
end do
end subroutine print_1rdms_from_sf2rdms_wrapper