output_molcas_rdms Subroutine

    subroutine output_molcas_rdms(rdm_defs, rdm, rdm_trace)

        !! Print spin-free GUGA RDMs directly in Molcas format

        type(rdm_definitions_t), intent(in) :: rdm_defs
            !! Type contanining the number of RDMs sampled and which states
            !! contribute to each RDM.
        type(rdm_list_t), intent(in) :: rdm
            !! Stores RDMs as 1D lists whose elements can be accessed through
            !! a hash table.
        real(dp), intent(in) :: rdm_trace(rdm%sign_length)
            !! Trace of RDMs required for normalisation of sampled arrays.

        real(dp), allocatable :: psmat(:), pamat(:), dmat(:)
        real(dp), parameter :: thresh = 1e-12_dp
        integer :: iunit_psmat, iunit_pamat, iunit_dmat, i, irdm
        character(*), parameter :: this_routine = "output_molcas_rdms"

        if (rdm_defs%nrdms_transition > 0) then
            call stop_all(this_routine,"GUGA transition RDMs yet to be implemented")
        end if

        do irdm = 1, rdm_defs%nrdms_standard
            call fill_molcas_rdms(rdm, rdm_trace, irdm, psmat, pamat, dmat)
            if (iProcIndex == root) then
                open(newunit=iunit_psmat, file='PSMAT.'//str(irdm), &
                     status='replace')
                    do i = 1, size(psmat)
                        if (abs(psmat(i)) > thresh) then
                            write(iunit_psmat, '(I6, G25.17)') i, psmat(i)
                        end if
                    end do
                close(iunit_psmat)

                open(newunit=iunit_pamat, file='PAMAT.'//str(irdm), &
                     status='replace')
                do i = 1, size(pamat)
                    if (abs(pamat(i)) > thresh) then
                        write(iunit_pamat, '(I6, G25.17)') i, pamat(i)
                    end if
                end do
                close(iunit_pamat)

                open(newunit=iunit_dmat, file='DMAT.'//str(irdm), &
                     status='replace')
                do i = 1, size(dmat)
                    if (abs(dmat(i)) > thresh) then
                        write(iunit_dmat, '(I6, G25.17)') i, dmat(i)
                    end if
                end do
                close(iunit_dmat)
            end if
        end do

    end subroutine output_molcas_rdms