CalcFOCKMatrix_RDM Subroutine

    subroutine CalcFOCKMatrix_RDM(rdm)

        ! Calculate the fock matrix in the natural orbital basis.

        use MemoryManager, only: LogMemAlloc, LogMemDealloc
        use rdm_data, only: one_rdm_t
        use SystemData, only: nbasis, ARR, BRR

        type(one_rdm_t), intent(in) :: rdm

        integer :: i, j, k, l, a, b, ierr, ArrDiagNewTag
        real(dp) :: FOCKDiagSumHF, FOCKDiagSumNew
        character(len=*), parameter :: t_r = 'CalcFOCKMatrix_RDM'
        real(dp), allocatable :: ArrDiagNew(:)

        ! This subroutine calculates and writes out the fock matrix for the
        ! transformed orbitals.
        ! ARR is originally the fock matrix in the HF basis.
        ! ARR(:,1) - ordered by energy, ARR(:,2) - ordered by spin-orbital index.

        ! When transforming the orbitals into approximate natural orbitals, we
        ! want to save memory, so don't bother calculating the whole matrix,
        ! just the diagonal elements that we actually need.

        allocate(ArrDiagNew(nBasis), stat=ierr)
        if (ierr /= 0) call Stop_All(t_r, 'Problem allocating ArrDiagNew array,')
        call LogMemAlloc('ArrDiagNew', nBasis, 8, t_r, ArrDiagNewTag, ierr)
        ArrDiagNew(:) = 0.0_dp

        ! First calculate the sum of the diagonal elements, ARR.
        ! Check if this is already being done.
        FOCKDiagSumHF = 0.0_dp
        do a = 1, nBasis
            FOCKDiagSumHF = FOCKDiagSumHF + Arr(a, 2)
        end do

        ! Then calculate the fock matrix in the transformed basis, and the sum
        ! of the new diagonal elements.

        FOCKDiagSumNew = 0.0_dp
        do j = 1, NoOrbs
            l = rdm%sym_list_no(j)
            if (tOpenShell) then
                ArrDiagNew(l) = 0.0_dp
            else
                ArrDiagNew(2 * l) = 0.0_dp
                ArrDiagNew((2 * l) - 1) = 0.0_dp
            end if
            do a = 1, NoOrbs
                b = rdm%sym_list_no(a)
                if (tOpenShell) then
                    ArrDiagNew(l) = ArrDiagNew(l) + (rdm%matrix(a, j) * ARR(b, 2) * rdm%matrix(a, j))
                else
                    ArrDiagNew(2 * l) = ArrDiagNew(2 * l) + (rdm%matrix(a, j) * ARR(2 * b, 2) * rdm%matrix(a, j))
                    ArrDiagNew((2 * l) - 1) = ArrDiagNew((2 * l) - 1) + (rdm%matrix(a, j) * ARR((2 * b) - 1, 2) * rdm%matrix(a, j))
                end if
            end do
            if (tOpenShell) then
                FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(l))
            else
                FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(2 * l))
                FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew((2 * l) - 1))
            end if
        end do
        ! If we are truncation the virtual space, only the unfrozen entries will
        ! be transformed.

        ! Refill ARR(:,1) (ordered in terms of energies), and ARR(:,2)
        ! (ordered in terms of orbital number).
        ! ARR(:,2) needs to be ordered in terms of symmetry and then energy
        ! (like SymLabelList), so currently this ordering will not be correct
        ! when reading in qchem INTDUMPs as the orbital number ordering is by energy.

        do j = 1, nBasis
            ARR(j, 2) = ArrDiagNew(j)
            ARR(j, 1) = ArrDiagNew(BRR(j))
        end do

        deallocate(ArrDiagNew)
        call LogMemDealloc(t_r, ArrDiagNewTag)

    end subroutine CalcFOCKMatrix_RDM