read_rs_lmat_factors Subroutine

    subroutine read_rs_lmat_factors
        character(*), parameter :: filename_mo = "mos_in_r"
        character(*), parameter :: filename_ints = "x_w_ij_r"
        character(*), parameter :: filename_info = "n_grid_pts_mo_num"
        character(*), parameter :: this_routine = "read_rs_lmat_factors"

        integer :: iunit, ierr, i, j, k, l
        integer(int64) :: ii, num_mos
        real(dp) :: integral
        real(dp), allocatable :: array_mos(:,:)

        root_print "Reading in range-separated TC factors"

        iunit = get_free_unit()
        ! Read the number of grid points and MOs
        open(iunit, file=filename_info, status='old')
        read(iunit, *, iostat=ierr) ngrid
        read(iunit, *, iostat=ierr) num_mos
        close(iunit)

        root_print "with: ", ngrid, " grid points"
        root_print "and ", num_mos, " orbitals"

        ! read the MOs file
        allocate(array_mos(ngrid, num_mos), source=0.0_dp)
        iunit = get_free_unit()
        open(iunit, file=filename_mo, status='old')
        do
            read(iunit, *, iostat=ierr) i, j, array_mos(i, j)
            if (ierr < 0) then
                exit
            else if (ierr > 0) then
                call stop_all(this_routine, "error reading " // filename_mo)
            end if
        end do
        close(iunit)

        ! fill in the qwprod file:
        allocate(qwprod(ngrid, num_mos, num_mos), source=0.0_dp)
        do i = 1, int(num_mos)
            do j = 1, int(num_mos)
                qwprod(:, i, j) = array_mos(:, i) * array_mos(:, j)
            end do
        end do
        deallocate(array_mos)

        allocate(ycoulomb(3, ngrid, num_mos, num_mos), source=0.0_dp)
        iunit = get_free_unit()
        open(iunit, file=filename_ints, status='old')
        do
            read(iunit, *, iostat=ierr) i, j, k, l, integral
            if (ierr < 0) then
                exit
            else if (ierr > 0) then
                call stop_all(this_routine, "error reading " // filename_ints)
            end if
            ycoulomb(j, i, k, l) = integral
        end do
        close(iunit)

        if (t_hash_lmat_calc) then
            lMatIndMax = lMatIndSym(num_mos, num_mos, num_mos, num_mos, num_mos, num_mos)
            lMatCalcHSize = int( real(lMatCalcHFactor,dp) * real(lMatIndMax,dp), int64)

            root_print "Total Size of LMat: ", lMatIndMax
            root_print "Size of LMatCalc Hash Table: ", lMatCalcHSize

            allocate(lMatCalcHKeys(lMatCalcHSize))
            allocate(lMatCalcHVals(lMatCalcHSize))
            do ii = 1_int64, lMatCalcHSize
                lMatCalcHKeys(i) = -1
            end do
            lMatCalcHit = 0
            lMatCalcTot = 0
            lMatCalcHUsed = 0
        end if

        root_print "Done: reading in range-separated TC factors"

    end subroutine read_rs_lmat_factors