print_1rdms_from_2rdms_wrapper Subroutine

    subroutine print_1rdms_from_2rdms_wrapper(rdm_defs, one_rdms, two_rdms, open_shell)

        ! Wrapper function to calculate and print 1-RDMs from 2-RDMs, as might
        ! be useful if the user forgot to print 1-RDMs in a calculation.

        use Parallel_neci, only: MPISumAll, iProcIndex
        use rdm_data, only: rdm_definitions_t, one_rdm_t, rdm_list_t
        use rdm_estimators, only: calc_rdm_trace
        use rdm_finalising, only: calc_1rdms_from_2rdms, write_1rdm, finalise_1e_rdm
        use SystemData, only: nel

        type(rdm_definitions_t), intent(in) :: rdm_defs
        type(one_rdm_t), intent(inout) :: one_rdms(:)
        type(rdm_list_t), intent(in) :: two_rdms
        logical, intent(in) :: open_shell

        integer :: irdm
        real(dp) :: rdm_trace(two_rdms%sign_length), rdm_trace_all(two_rdms%sign_length)
        real(dp) :: rdm_norm_all(two_rdms%sign_length), norm_1rdm(size(one_rdms))
        character(len=*), parameter :: t_r = 'print_1rdms_from_2rdms_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

        !call calc_rdm_trace(two_rdms, rdm_trace)
        !call MPISumAll(rdm_trace, rdm_trace_all)
        ! RDMs are normalised so that their trace is nel*(nel-1)/2.
        !rdm_norm_all = rdm_trace_all*2.0_dp/(nel*(nel-1))

        ! The read-in spinfree 2-RDMs should already be normalised.
        rdm_norm_all = 1.0_dp
        norm_1rdm = 1.0_dp

        call calc_1rdms_from_2rdms(rdm_defs, one_rdms, two_rdms, rdm_norm_all, open_shell)

        ! we need to finalise the 1rdms
        call finalise_1e_rdm(rdm_defs, one_rdms, norm_1rdm)

        if (iProcIndex == 0) then
            do irdm = 1, size(one_rdms)
                call write_1rdm(rdm_defs, one_rdms(irdm)%matrix, irdm, norm_1rdm(irdm), .true.)
            end do
        end if

    end subroutine print_1rdms_from_2rdms_wrapper