  Constructing accurate, high dimensional molecular potential energy surfaces
(PESs) for polyatomic molecules is challenging. Reproducing Kernel Hilbert
space (RKHS) interpolation is an efficient way to construct such PESs. However,
the scheme is most effective when the input energies are available on a regular
grid. Thus the number of reference energies required can become very large even
for penta-atomic systems making such an approach computationally prohibitive
when using high-level electronic structure calculations. Here an efficient and
robust scheme is presented to overcome these limitations and is applied to
constructing high dimensional PESs for systems with up to 10 atoms. Using
energies as well as gradients reduces the number of input data required and
thus keeps the number of coefficients at a manageable size. Correct
implementation of permutational symmetry in the kernel products is tested and
explicitly demonstrated for the highly symmetric CH$_4$ molecule.
