RefillUMATandTMAT2D_RDM Subroutine

    subroutine RefillUMATandTMAT2D_RDM(one_rdm, sym_list)

        ! This routine refills these to more easily write out the ROFCIDUMP,
        ! and originally to be able to continue a calculation (although I doubt
        ! this works at the moment).

        ! UMat is in spin or spatial orbitals, TMAT2D only spin.

        use IntegralsData, only: umat
        use MemoryManager, only: LogMemAlloc, LogMemDealloc
        use OneEInts, only: TMAT2D
        use rdm_data, only: tRotatedNOs
        use RotateOrbsMod, only: FourIndInts
        use SystemData, only: nbasis
        use UMatCache, only: UMatInd

        real(dp), intent(in) :: one_rdm(:, :)
        integer, intent(in) :: sym_list(:)

        integer :: l, k, j, i, a, b, g, d, c
        integer :: nBasis2, TMAT2DPartTag, ierr
        real(dp) :: NewTMAT
        real(dp), allocatable :: TMAT2DPart(:, :)
        character(len=*), parameter :: t_r = 'RefillUMATandTMAT2D_RDM'

#ifdef CMPLX_
        call stop_all('RefillUMATandTMAT2D_RDM', &
                      'Rotating orbitals not implemented for complex orbitals.')
#endif

        ! TMAT2D is always in spin orbitals.
        allocate(TMAT2DPart(nBasis, nBasis), stat=ierr)
        if (ierr /= 0) call Stop_All(t_r, 'Problem allocating TMAT2DPart array,')
        call LogMemAlloc('TMAT2DPart', nBasis * nBasis, 8, &
                         'RefillUMAT_RDM', TMAT2DPartTag, ierr)
        TMAT2DPart = 0.0_dp

        ! Make the UMAT elements the four index integrals.
        ! These are calculated by transforming the HF orbitals using the
        ! coefficients that have been found.
        do l = 1, NoOrbs
            d = sym_list(l)
            do k = 1, NoOrbs
                b = sym_list(k)
                do j = 1, NoOrbs
                    g = sym_list(j)
                    do i = 1, NoOrbs
                        a = sym_list(i)
                        if (tOpenSpatialOrbs) then
                            a = gtID(a)
                            b = gtID(b)
                            g = gtID(g)
                            d = gtID(d)
                        end if
                        ! The FourIndInts are in chemical notation, the UMatInd
                        ! in physical.
                        UMAT(UMatInd(a, b, g, d)) = FourIndInts(i, j, k, l)
                    end do
                end do
            end do
        end do

        ! Also calculate the 2 index integrals, and make these the elements
        ! of the TMAT2D matrix. TMAT2D is in spin orbitals.

        do a = 1, nBasis
            do k = 1, NoOrbs
                i = sym_list(k)
                NewTMAT = 0.0_dp
                do b = 1, NoOrbs
                    d = sym_list(b)
                    if (tOpenShell) then
                        NewTMAT = NewTMAT + (one_rdm(b, k) * real(TMAT2D(d, a), dp))
                    else
                        NewTMAT = NewTMAT + (one_rdm(b, k) * real(TMAT2D(2 * d, a), dp))
                    end if
                end do
                if (tOpenShell) then
                    TMAT2DPart(i, a) = NewTMAT
                else
                    if (mod(a, 2) == 0) then
                        TMAT2DPart(2 * i, a) = NewTMAT
                    else
                        TMAT2DPart((2 * i) - 1, a) = NewTMAT
                    end if
                end if
            end do
        end do

        do k = 1, nBasis
            do l = 1, NoOrbs
                j = sym_list(l)
                NewTMAT = 0.0_dp
                do a = 1, NoOrbs
                    c = sym_list(a)
                    if (tOpenShell) then
                        NewTMAT = NewTMAT + (one_rdm(a, l) * TMAT2DPart(k, c))
                    else
                        NewTMAT = NewTMAT + (one_rdm(a, l) * TMAT2DPart(k, 2 * c))
                    end if
                end do
                if (tOpenShell) then
                    TMAT2D(k, j) = NewTMAT
                else
                    if (mod(k, 2) == 0) then
                        TMAT2D(k, 2 * j) = NewTMAT
                    else
                        TMAT2D(k, (2 * j) - 1) = NewTMAT
                    end if
                end if
            end do
        end do

        deallocate(TMAT2DPart)
        call LogMemDeAlloc('RefillUMAT_RDM', TMAT2DPartTag)

        if (.not. tRotatedNOs) then
            call PrintROFCIDUMP_RDM("ROFCIDUMP")
        end if

    end subroutine RefillUMATandTMAT2D_RDM