Transform2ElIntsMemSave_RDM Subroutine

    subroutine Transform2ElIntsMemSave_RDM(one_rdm, sym_list)

        ! This is an M^5 transform, which transforms all the two-electron
        ! integrals into the new basis described by the Coeff matrix.
        ! This is very memory inefficient and currently does not use any
        ! spatial symmetry information.

        use IntegralsData, only: umat
        use MemoryManager, only: LogMemAlloc, LogMemDealloc
        use RotateOrbsMod, only: FourIndInts
        use UMatCache, only: UMatInd

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

        integer :: i, j, k, l, a, b, g, d, ierr, Temp4indintsTag, a2, d2, b2, g2
        real(dp), allocatable :: Temp4indints(:, :)
        character(len=*), parameter :: t_r = 'Transform2ElIntsMemSave_RDM'
#ifdef CMPLX_
        call stop_all('Transform2ElIntsMemSave_RDM', &
                      'Rotating orbitals not implemented for complex orbitals.')
#endif

        ! Zero arrays from previous transform.

        allocate(Temp4indints(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('Temp4indints', NoOrbs**2, 8, &
                         'Transform2ElIntsMemSave_RDM', Temp4indintsTag, ierr)
        if (ierr /= 0) call Stop_All('Transform2ElIntsMemSave_RDM', &
                                     'Problem allocating memory to Temp4indints.')

        FourIndInts(:, :, :, :) = 0.0_dp

        ! Calculating the two-transformed, four index integrals.

        ! The untransformed <alpha beta | gamma delta> integrals are found from
        ! UMAT(UMatInd(i,j,k,l)
        ! All our arrays are in spin orbitals - if tStoreSpinOrbs is false,
        ! UMAT will be in spatial orbitals - need to account for this.

        ! Running through 1,NoOrbs - the actual orbitals corresponding to that
        ! index are given by sym_list.

        do b = 1, NoOrbs
            b2 = sym_list(b)
            do d = 1, NoOrbs
                d2 = sym_list(d)
                do a = 1, NoOrbs
                    a2 = sym_list(a)
                    do g = 1, NoOrbs
                        g2 = sym_list(g)
                        ! if tStoreSpinOrbs is false but NoOrbs == nBasis, the UMAT
                        ! indices have to be converted
                        if (tOpenSpatialOrbs) then
                            a2 = gtID(a2)
                            b2 = gtID(b2)
                            g2 = gtID(g2)
                            d2 = gtID(d2)
                        end if
                        ! UMatInd in physical notation, but FourIndInts in
                        ! chemical (just to make it more clear in these
                        ! transformations). This means that here, a and g are
                        ! interchangable, and so are b and d.
                        FourIndInts(a, g, b, d) = real(UMAT(UMatInd(a2, b2, g2, d2)), dp)
                    end do
                end do

                Temp4indints(:, :) = 0.0_dp
                call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, one_rdm, NoOrbs, &
                           FourIndInts(1:NoOrbs, 1:NoOrbs, b, d), NoOrbs, 0.0_dp, &
                           Temp4indints(1:NoOrbs, 1:NoOrbs), NoOrbs)

                ! Temp4indints(i,g) comes out of here, so to transform g to k,
                ! we need the transpose of this.

                call dgemm('T', 'T', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, one_rdm, NoOrbs, &
                           Temp4indints(1:NoOrbs, 1:NoOrbs), NoOrbs, 0.0_dp, &
                           FourIndInts(1:NoOrbs, 1:NoOrbs, b, d), NoOrbs)
                ! Get Temp4indits02(i,k).
            end do
        end do

        ! Calculating the 3 transformed, 4 index integrals.
        ! 01=a untransformed,02=b,03=g,04=d.
        do i = 1, NoOrbs
            do k = 1, NoOrbs

                Temp4indints(:, :) = 0.0_dp
                call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, one_rdm, NoOrbs, &
                           FourIndInts(i, k, 1:NoOrbs, 1:NoOrbs), NoOrbs, 0.0_dp, &
                           Temp4indints(1:NoOrbs, 1:NoOrbs), NoOrbs)

                call dgemm('T', 'T', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, one_rdm, &
                           NoOrbs, Temp4indints(1:NoOrbs, 1:NoOrbs), NoOrbs, 0.0_dp, &
                           FourIndInts(i, k, 1:NoOrbs, 1:NoOrbs), NoOrbs)
            end do
        end do

        deallocate(Temp4indints)
        call LogMemDeAlloc('Transform2ElIntsMemSave_RDM', Temp4indintsTag)

    end subroutine Transform2ElIntsMemSave_RDM