readLMatFactors Subroutine

    subroutine readLMatFactors()

        character(*), parameter :: filename = "tcfactors.h5"
        character(*), parameter :: nm_grp = "tcfactors", nm_nBasis = "nBasis", nm_nGrid = "nGrid", &
                                   nm_weights = "weights", nm_mo_vals = "mo_vals", nm_ycoulomb = "ycoulomb"
        real(dp), allocatable :: mo_vals(:, :), weights(:)
#ifdef USE_HDF5_
        integer(hid_t) :: err, file_id, grp_id, dataset, type_id
        integer(hsize_t) :: weights_dims(1), mo_vals_dims(2), ycoulomb_dims(4)
#endif
        integer i, a, b
        character(*), parameter :: this_routine = "readLMatFactors"
#ifdef USE_HDF5_
        write(stdout, *) "**************** Reading TcFactors File ****************"
        call h5open_f(err)

        ! open the file
        call h5fopen_f(filename, H5F_ACC_RDONLY_F, file_id, err)

        call h5gopen_f(file_id, nm_grp, grp_id, err)

        call read_int64_attribute(grp_id, nm_nBasis, nBasis, required=.true.)
        write(stdout, *) "Number of basis function: ", nBasis

        call read_int64_attribute(grp_id, nm_nGrid, nGrid, required=.true.)
        write(stdout, *) "Number of grid points: ", nGrid

        ! load weights
        allocate(weights(nGrid))
        weights_dims = size(weights)
        call h5dopen_f(grp_id, nm_weights, dataset, err)
        call h5dget_type_f(dataset, type_id, err)
        call h5dread_f(dataset, type_id, weights, weights_dims, err)
        call h5tclose_f(type_id, err)
        call h5dclose_f(dataset, err)

        ! load mo_vals
        allocate(mo_vals(nGrid, nBasis))
        mo_vals_dims = size(mo_vals)
        call h5dopen_f(grp_id, nm_mo_vals, dataset, err)
        call h5dget_type_f(dataset, type_id, err)
        call h5dread_f(dataset, type_id, mo_vals, mo_vals_dims, err)
        call h5tclose_f(type_id, err)
        call h5dclose_f(dataset, err)

        ! load ycoulomb
        allocate(ycoulomb(3, nGrid, nBasis, nBasis))
        ycoulomb_dims = size(ycoulomb)
        call h5dopen_f(grp_id, nm_ycoulomb, dataset, err)
        call h5dget_type_f(dataset, type_id, err)
        call h5dread_f(dataset, type_id, ycoulomb, ycoulomb_dims, err)
        call h5tclose_f(type_id, err)
        call h5dclose_f(dataset, err)

        ! close the file, finalize hdf5
        call h5gclose_f(grp_id, err)
        call h5fclose_f(file_id, err)
        call h5close_f(err)

        !Combine weights and molecular oribtals
        allocate(qwprod(nGrid, nBasis, nBasis))
        do a = 1, nBasis
            do b = 1, nBasis
                qwprod(:, a, b) = mo_vals(:, a) * mo_vals(:, b) * weights
            end do
        end do
        deallocate(mo_vals)
        deallocate(weights)

        !ycoulom and ycoulom2 are symmetric in MO. However, only the upper part is filled
        !in the script producing the HDF5 file. Therefore, we fill the lower part here
        do a = 1, nBasis
            do b = a + 1, nBasis
                ycoulomb(:, :, a, b) = ycoulomb(:, :, b, a)
            end do
        end do

        lMatIndMax = lMatIndSym(nBasis, nBasis, nBasis, nBasis, nBasis, nBasis)
        lMatCalcHSize = lMatCalcHFactor * lMatIndMax

        write(stdout, *) "Total Size of LMat: ", lMatIndMax
        write(stdout, *) "Size of LMatCalc Hash Table: ", lMatCalcHSize

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

        write(stdout, *) "********************************************************"
#else
        call stop_all(this_routine, 'HDF5 support not enabled at compile time')
        unused_var(filename)
#endif
    end subroutine readLMatFactors