Modelling potential energy surfaces for small clusters using Shepard interpolation with Gaussian-form nodal functions 利用谢泼德插值法和高斯形式节点函数为小集群的势能面建模
Haina Wang (DD)a and Ryan P. A. Bettens (DD) 王海娜（副司长）a 和 Ryan P. A. Bettens（副司长）
Abstract 摘要
The potential energy surface (PES) of a chemical system is an analytical function that outputs the potential energy of the system when a nuclear configuration is given as input. The PESs of small atmospheric clusters have theoretical as well as environmental significance. A common method used to generate analytical PESs is the Shepard interpolation, where the PES is a weighed sum of Taylor series expansions (nodal functions) at ab initio sample points. Based on this, in this study we present a new method based on the Shepard interpolation, where the nodal functions are composed of a symmetric Gaussian term and an asymmetric exponential term in each dimension. Corresponding sampling methods were also developed. We tested the method on several atmospheric bimolecular clusters and achieved root mean square errors (RMSE) below 0.13kJmol^(-1)0.13 \mathrm{~kJ} \mathrm{~mol}^{-1} in 150 samples for Ar-rigid H_(2)O\mathrm{H}_{2} \mathrm{O} and Ne -rigid CO_(2)\mathrm{CO}_{2}, and below 0.39kJmol^(-1)0.39 \mathrm{~kJ} \mathrm{~mol}^{-1} in 1800 samples for rigid N_(2)\mathrm{N}_{2}-rigid CO_(2)\mathrm{CO}_{2}. 化学系统的势能面（PES）是一个分析函数，当输入核构型时，它输出系统的势能。小型大气团块的势能面具有理论和环境意义。生成分析 PES 的常用方法是 Shepard 插值法，其中 PES 是在 ab initio 样本点上的泰勒级数展开（节点函数）的加权和。在此基础上，我们在本研究中提出了一种基于 Shepard 插值的新方法，其中节点函数由每个维度上的对称高斯项和非对称指数项组成。同时还开发了相应的采样方法。我们在几个大气双分子簇上对该方法进行了测试，结果表明氩-刚性 H_(2)O\mathrm{H}_{2} \mathrm{O} 和氖-刚性 CO_(2)\mathrm{CO}_{2} 的均方根误差（RMSE）在 150 个样本内低于 0.13kJmol^(-1)0.13 \mathrm{~kJ} \mathrm{~mol}^{-1} ，刚性 N_(2)\mathrm{N}_{2} 和刚性 CO_(2)\mathrm{CO}_{2} 的均方根误差在 1800 个样本内低于 0.39kJmol^(-1)0.39 \mathrm{~kJ} \mathrm{~mol}^{-1} 。
Introduction 简介
The potential energy surface (PES) of a chemical system is an important concept in theoretical chemistry that has wide applications including spectral analysis, ^(1){ }^{1} protein folding, ^(2-4){ }^{2-4} and reaction dynamics. ^(5){ }^{5} Based on the Born-Oppenheimer approximation, the PES refers to a function that outputs the electronic potential energy of the system when a nuclear configuration is given as input. For a nonlinear molecule with NN atoms, the PES is a function of 3N-63 N-6 independent degrees of freedom. ^(6,7){ }^{6,7} However, “redundant dimensions” can be present depending on the choice of the coordinate system that describes the nuclear configuration. In the 1990s and early 2000s, there has been controversy or reservation about using all internuclear distances as coordinates because of this “redundancy”, and there were suggestions of using different sets of 3N-63 N-6 internal coordinates in different regions of the configuration space, e.g. in ref. 8. However, starting with Bowman’s group, recent decades have seen successful uses of all internuclear distances, or “Morse variables” derived therefrom. ^(9,10){ }^{9,10} 化学体系的势能面（PES）是理论化学中的一个重要概念，应用广泛，包括光谱分析、 ^(1){ }^{1} 蛋白质折叠、 ^(2-4){ }^{2-4} 和反应动力学。 ^(5){ }^{5} 基于玻恩-奥本海默近似，势能面是指在输入核构型时输出系统电子势能的函数。对于 NN 原子的非线性分子，PES 是 3N-63 N-6 独立自由度的函数。 ^(6,7){ }^{6,7} 然而，"冗余维度 "可能存在，这取决于描述核构型的坐标系的选择。在 20 世纪 90 年代和 21 世纪初，由于这种 "冗余性"，对使用所有核间距作为坐标一直存在争议或保留意见，有人建议在构型空间的不同区域使用不同的 3N-63 N-6 内部坐标集，例如参考文献 8。8.然而，从鲍曼研究小组开始，近几十年来所有核间距或由此衍生的 "莫尔斯变量 "都得到了成功应用。 ^(9,10){ }^{9,10}
Recently, due to the increasing awareness of human impact on the atmosphere, there has been extensive theoretical interest in the PESs of atmospheric clusters. ^(11-13){ }^{11-13} These are small systems of common atmospheric molecules, e.g. N_(2),O_(2),O_(3)\mathrm{N}_{2}, \mathrm{O}_{2}, \mathrm{O}_{3}, 最近，由于人们越来越意识到人类对大气的影响，理论界对大气团块的 PES 产生了广泛的兴趣。 ^(11-13){ }^{11-13} 这些是由常见大气分子（如 N_(2),O_(2),O_(3)\mathrm{N}_{2}, \mathrm{O}_{2}, \mathrm{O}_{3} ）组成的小系统、
H_(2)O\mathrm{H}_{2} \mathrm{O}, Ar and CO_(2)\mathrm{CO}_{2}, interacting through van der Waals forces or hydrogen bonds. Careful investigations of these clusters can shed new light on current environmental issues concerning atmospheric chemistry. For example, it has been predicted that as CO_(2)\mathrm{CO}_{2} forms clusters, its usually symmetric stretch mode nu_(1)\nu_{1} can be significantly asymmetric and absorb infrared radiation, contributing to the greenhouse effect together with its usually asymmetric vibrational modes, nu_(2)\nu_{2} and nu_(3).^(11)\nu_{3} .{ }^{11} H_(2)O\mathrm{H}_{2} \mathrm{O} 、Ar 和 CO_(2)\mathrm{CO}_{2} ，通过范德华力或氢键相互作用。对这些团簇的仔细研究可以为当前有关大气化学的环境问题提供新的启示。例如，据预测，当 CO_(2)\mathrm{CO}_{2} 形成簇时，其通常对称的伸展模式 nu_(1)\nu_{1} 可能会明显不对称，并吸收红外辐射，与通常不对称的振动模式 nu_(2)\nu_{2} 和 一起造成温室效应。 nu_(3).^(11)\nu_{3} .{ }^{11}
Due to the small sizes of atmospheric clusters, ab initio single point energy (SPE) evaluations with very high accuracy (1.0(kJ)mol^(-1))\left(1.0 \mathrm{~kJ} \mathrm{~mol}^{-1}\right) can be achieved. ^(14){ }^{14} However, a thorough understanding of their possible interactions demands the generation of accurate analytical PESs over a global range of chemically and physically interesting configurations using a finite set of ab initio data, which remains a challenging task. The advantages of an analytical PES over pointwise evaluations include faster SPE retrieval, easy mathematical treatment and convenience of visualization. Typically, generating an analytical PES involves three aspects: the basic form of the function, the sampling of the configuration space (for aba b initio calculations), and finally, the fitting or interpolation method that gives the global PES. ^(6){ }^{6} 由于大气团簇的尺寸较小，因此可以实现精度极高的非初始单点能量（SPE）评估 (1.0(kJ)mol^(-1))\left(1.0 \mathrm{~kJ} \mathrm{~mol}^{-1}\right) 。 ^(14){ }^{14} 然而，要彻底了解它们之间可能存在的相互作用，就必须利用有限的 ab initio 数据集，在全局范围内对化学和物理上有趣的构型生成精确的分析 PES，这仍然是一项具有挑战性的任务。与点式评估相比，解析 PES 的优势在于 SPE 检索速度更快、数学处理简单、可视化方便。通常情况下，生成分析型 PES 涉及三个方面：函数的基本形式、构型空间采样（用于 aba b initio 计算）以及最后给出全局 PES 的拟合或插值方法。 ^(6){ }^{6}
There are many milestones for all these aspects. The requirement for the PES to be symmetric with respect to permutations of like atoms has been explicitly stated by Murrell and colleagues in their classic monograph. ^(15){ }^{15} To incorporate this permutational symmetry, theories of symmetrized coordinates ^(16){ }^{16} and invariant fitting bases ^(17){ }^{17} have been developed for triatomic and homonuclear tetra-atomic systems. More recently, the least squares fitting of permutationally invariant polynomials (PIP) was developed for high dimensionality 在所有这些方面都有许多里程碑式的成果。穆雷尔及其同事在他们的经典专著中明确提出了 PES 对于同类原子排列对称的要求。 ^(15){ }^{15} 为了纳入这种排列对称性，针对三原子和同核四原子系统开发了对称坐标理论 ^(16){ }^{16} 和不变拟合基 ^(17){ }^{17} 。最近，针对高维度的包络不变多项式 (PIP) 的最小二乘法拟合也得到了发展。
by Braams and Bowman. ^(18){ }^{18} The many-body expansion (MBE), first applied by Varandas and Murrell to PESs of H_(n)\mathrm{H}_{\mathrm{n}} systems, ^(16){ }^{16} has become a widely used approach to model asymptotic regions of configuration spaces. ^(6,15,19-21){ }^{6,15,19-21} Common sampling methods include random sampling from a grid of bond lengths and angles. There are, however, more sophisticated sampling methods such as Mutually Orthogonal Latin Squares that aim to cover a wide range of a high dimensional space. ^(4){ }^{4} The Bowman group has used scattered sampling methods, where samples are generated from molecular dynamics simulation data with widely different trajectories, effectively sampling a large range of configurations and energies. ^(22){ }^{22} 由 Braams 和 Bowman 合著。 ^(18){ }^{18} 多体展开（MBE）最早由 Varandas 和 Murrell 应用于 H_(n)\mathrm{H}_{\mathrm{n}} 系统的 PES， ^(16){ }^{16} ，现已成为一种广泛应用于构型空间渐近区域建模的方法。 ^(6,15,19-21){ }^{6,15,19-21} 常见的取样方法包括从键长和角度网格中随机取样。不过，还有更复杂的取样方法，例如互正交拉丁方阵，旨在覆盖高维空间的广泛范围。 ^(4){ }^{4} 鲍曼研究小组使用了分散采样法，即从分子动力学模拟数据中生成轨迹差异很大的样本，从而有效地对大量构型和能量进行采样。 ^(22){ }^{22}
One of the most successful interpolation methods is the modified Shepard interpolation proposed by Collins and colleagues in 1994. ^(8,23){ }^{8,23} It represents the PES as a weighed sum of Taylor expansions in inverse internuclear distances up to the second order on each sample. The weight of a sample at a configuration is inversely proportional to the distance between the configuration and the sample. In addition, the group proposes an iterative sampling scheme that outperforms random sampling, adding in each round new samples that are distant from existing samples. Problems of the modified Shepard method include bump- or step-shaped artifacts where distinct samples compete for dominance. Various studies have been focusing on reducing such effect. ^(20,24){ }^{20,24} Another problem with the modified Shepard interpolation is that the number of ab initio points needed to calculate the second derivatives grows as the square of the number of dimensions, making the method computationally expensive for large systems. 最成功的内插法之一是柯林斯及其同事于 1994 年提出的修正谢泼德内插法。 ^(8,23){ }^{8,23} 它将 PES 表示为每个样本上核间距逆二阶泰勒展开的加权和。某个构型的样本权重与构型和样本之间的距离成反比。此外，研究小组还提出了一种优于随机抽样的迭代抽样方案，即在每一轮中增加与现有样本距离较远的新样本。修改后的谢泼德方法存在的问题包括凹凸或阶梯形伪影，在这些伪影中，不同的样本会争夺主导地位。各种研究都在关注如何减少这种影响。 ^(20,24){ }^{20,24} 改良谢泼德插值法的另一个问题是，计算二阶导数所需的原子序数点数随着维数的平方而增长，这使得该方法在大型系统中的计算成本非常昂贵。
More recently, some groups have studied the application of Gaussian Process (GP) in PES generation. ^(25-30){ }^{25-30} As another method of interpolation, it uses Gaussians instead of Taylor expansion as nodal functions. Being a less dramatically varying function than the polynomial, the Gaussian has promises to reduce artifacts. It also provides a natural estimation of error at not-yet sampled configurations, thereby providing an iterative sampling scheme where points with large estimated errors are added to the sample set. GP has been used successfully in modeling three-dimensional PESs of proton transfer in crystals. ^(25){ }^{25} Very recently, PIP-GP approaches have been developed to incorporate permutation invariance into GP. ^(30){ }^{30} The application of GP to higher dimensional gaseous systems, however, has been largely limited by the inflexible width of the Gaussian nodal functions and the fact that Gaussians are inevitably symmetric at samples, making it hard to account for first derivatives. 最近，一些研究小组研究了高斯过程（GP）在 PES 生成中的应用。 ^(25-30){ }^{25-30} 作为另一种插值方法，它使用高斯而不是泰勒展开作为节点函数。与多项式相比，高斯函数的变化幅度较小，因此有望减少伪影。它还能自然地估计尚未采样配置的误差，从而提供一种迭代采样方案，将估计误差较大的点添加到采样集中。GP 已成功用于晶体中质子传递的三维 PES 建模。 ^(25){ }^{25} 最近，人们开发了 PIP-GP 方法，将包换不变性纳入 GP。 ^(30){ }^{30} 然而，由于高斯节点函数的宽度不够灵活，而且高斯在样本处不可避免地具有对称性，因此很难考虑一阶导数，这在很大程度上限制了 GP 在更高维度气态系统中的应用。
In this study, we developed a new form of nodal functions that aims to tackle the limitations of GP. It composes of a symmetric Gaussian part and an asymmetric inverse-exponential part for each dimension, both having flexible widths and amplitudes. We also proposed an iterative sampling scheme compatible with these nodal functions. We tested our method on three atmospheric clusters: Ar-rigid H_(2)O,^(31)\mathrm{H}_{2} \mathrm{O},{ }^{31} Ne-rigid CO_(2),^(28)\mathrm{CO}_{2},{ }^{28} and rigid CO_(2)\mathrm{CO}_{2}-rigid N_(2).^(11)\mathrm{N}_{2} .{ }^{11} The focus for this present paper is on non-covalent interactions, with expectations of possible application to atmospheric science. For example, the PES for the clusters can serve as a perturbation term for the rovibration spectrum of interacting atmospheric molecules, with the monomers as the “base case”, in a similar spirit to ref. 32. Nevertheless, we believe that the method itself is general enough to also deal with reacting molecules. 在这项研究中，我们开发了一种新的节点函数形式，旨在解决 GP 的局限性。它由每个维度的对称高斯部分和非对称反指数部分组成，两者都具有灵活的宽度和振幅。我们还提出了一种与这些节点函数兼容的迭代采样方案。我们在三个大气集群上测试了我们的方法： H_(2)O,^(31)\mathrm{H}_{2} \mathrm{O},{ }^{31}CO_(2),^(28)\mathrm{CO}_{2},{ }^{28}CO_(2)\mathrm{CO}_{2}N_(2).^(11)\mathrm{N}_{2} .{ }^{11} 本文的重点是非共价相互作用，希望能应用于大气科学。例如，以单体为 "基例"，簇的 PES 可以作为大气分子相互作用振动谱的扰动项，这与参考文献 32 的精神类似。32.尽管如此，我们相信该方法本身的通用性足以处理反应分子。
We focus only on interpolation of ab initio data in this current study, while noting that the dissociation regime can be well described by MBE methods, and that potential functions accurate for different regimes can be combined to give a global PES, e.g. with the energy switching (ES) approach developed by Varandas. ^(33-35){ }^{33-35} 在目前的研究中，我们只关注 ab initio 数据的插值，同时注意到解离机制可以用 MBE 方法很好地描述，而且不同机制的精确势函数可以结合起来给出全局 PES，例如用 Varandas 开发的能量转换 (ES) 方法。 ^(33-35){ }^{33-35}
We will first explain the general methodology and then discuss the specific implementations and results for the different clusters. 我们将首先解释一般方法，然后讨论不同集群的具体实施和结果。
Methodology 方法
Nodal functions 节点功能
We use internuclear distances to describe configurations. As all clusters we study are rigid, only intermolecular internuclear distances are included. For example, in the cluster Ar-rigid H_(2)O\mathrm{H}_{2} \mathrm{O}, configurations are defined by the internuclear distances ( Ar-O,Ar-H^(1),Ar-H^(2)\mathrm{Ar}-\mathrm{O}, \mathrm{Ar}-\mathrm{H}^{1}, \mathrm{Ar}-\mathrm{H}^{2} ). Given the configuration vector r\mathbf{r} of a sample with ab initio energy E_(r)E_{\mathbf{r}}, the nodal function V_(r)V_{\mathbf{r}} is expressed by a Gaussian symmetric term (resembling the s orbital) and an asymmetric term (resembling the p orbital) in each dimension, plus the constant E_(r)E_{\mathbf{r}}, ensuring that the nodal function passes through the sample, i.e. V_(r)(r)=E_(r)V_{\mathbf{r}}(\mathbf{r})=E_{\mathbf{r}}, as shown in (1). 我们使用核间距来描述构型。由于我们研究的所有原子团都是刚性的，因此只包括分子间核间距。例如，在氩刚性原子团簇 H_(2)O\mathrm{H}_{2} \mathrm{O} 中，构型由核间距 ( Ar-O,Ar-H^(1),Ar-H^(2)\mathrm{Ar}-\mathrm{O}, \mathrm{Ar}-\mathrm{H}^{1}, \mathrm{Ar}-\mathrm{H}^{2} ) 定义。给定样本的构型矢量 r\mathbf{r} ，初始能量为 E_(r)E_{\mathbf{r}} ，结点函数 V_(r)V_{\mathbf{r}} 由每个维度上的一个高斯对称项（类似于 s 轨道）和一个不对称项（类似于 p 轨道）加上常数 E_(r)E_{\mathbf{r}} 来表示，以确保结点函数穿过样本，即 V_(r)(r)=E_(r)V_{\mathbf{r}}(\mathbf{r})=E_{\mathbf{r}} ，如 (1) 所示。
where tilde(r)\tilde{\mathbf{r}} is any arbitrary configuration in the vicinity of r\mathbf{r}. For each dimension, the parameters A_(i)A_{i} and alpha_(i)\alpha_{i} are amplitude and width of the asymmetric term. B_(i)B_{i} and beta_(i)\beta_{i} are amplitude and width of the symmetric term. alpha_(i)\alpha_{i} and beta_(i)\beta_{i} must be positive. 其中 tilde(r)\tilde{\mathbf{r}} 是 r\mathbf{r} 附近的任意配置。对于每个维度，参数 A_(i)A_{i} 和 alpha_(i)\alpha_{i} 是不对称项的振幅和宽度。 B_(i)B_{i} 和 beta_(i)\beta_{i} 是对称项的振幅和宽度。 alpha_(i)\alpha_{i} 和 beta_(i)\beta_{i} 必须为正值。
All parameters A_(i),B_(i),alpha_(i),beta_(i)A_{i}, B_{i}, \alpha_{i}, \beta_{i} need to be determined. The most obvious way is to fit the parameters based on ab initio calculations on a set of configurations surrounding each sample, which we terms as co-samples. We then optimize the width parameters with a rather “brute force” line search, first letting alpha_(1)\alpha_{1} increase from 0.01 to 50 with increment 0.01 , while all other widths are initially fixed at 0.01 . For each value of alpha_(1)\alpha_{1}, the amplitude parameters A_(i)A_{i} and B_(i)B_{i} are obtained by linear regression; a least-squares regression error is thus generated. The value of alpha_(1)\alpha_{1} that gives the smallest regression error is taken to be optimal and fixed, while the next width parameter, e.g. alpha_(2)\alpha_{2}, undergoes optimization, and so on for all alpha_(i)\alpha_{i}, beta_(i)\beta_{i}. After all widths undergo optimization, the whole process starts over from alpha_(1)\alpha_{1} again and repeats for totally three rounds. This regression scheme apparently does not search for the global optimum of alpha_(i)\alpha_{i} and beta_(i)\beta_{i}, but as we will see in the next section, it turns out to be a rather practical scheme. 所有参数 A_(i),B_(i),alpha_(i),beta_(i)A_{i}, B_{i}, \alpha_{i}, \beta_{i} 都需要确定。最明显的方法是根据围绕每个样本的一组构型（我们称之为共样本）的 ab initio 计算来拟合参数。然后，我们用一种相当 "粗暴 "的直线搜索来优化宽度参数，首先让 alpha_(1)\alpha_{1} 从 0.01 增加到 50，增量为 0.01，而所有其他宽度最初都固定为 0.01。对于 alpha_(1)\alpha_{1} 的每个值，振幅参数 A_(i)A_{i} 和 B_(i)B_{i} 都是通过线性回归得到的；因此会产生最小二乘回归误差。得出最小回归误差的 alpha_(1)\alpha_{1} 值被视为最优值并固定下来，同时对下一个宽度参数（如 alpha_(2)\alpha_{2} ）进行优化，依此类推，对所有 alpha_(i)\alpha_{i} , beta_(i)\beta_{i} 进行优化。所有宽度参数优化后，整个过程从 alpha_(1)\alpha_{1} 重新开始，共重复三轮。这种回归方案显然没有搜索 alpha_(i)\alpha_{i} 和 beta_(i)\beta_{i} 的全局最优值，但正如我们将在下一节看到的，它被证明是一种相当实用的方案。
The permutational invariance, i.e. the requirement that when two atoms of the identical element switch position, the potential energy should remain the same, can be incorporated easily. Whenever a sample is chosen, we also add its “permutational twins” into the sample set. The nodal functions of the twin samples are obtained by permuting the terms of the nodal function of 我们可以很容易地加入排列不变性，即当两个相同元素的原子交换位置时，其势能应保持不变。每当选择一个样本时，我们也会将其 "孪生包络 "加入样本集中。孪生样本的节点函数是通过对
the original sample. For example, if a sample of Ar-rigid H_(2)O\mathrm{H}_{2} \mathrm{O} is given by the internuclear distances (Ar-O,Ar-H^(1),Ar-H^(2))=\left(\mathrm{Ar}-\mathrm{O}, \mathrm{Ar}-\mathrm{H}^{1}, \mathrm{Ar}-\mathrm{H}^{2}\right)=(r_(1),r_(2),r_(3))=r\left(r_{1}, r_{2}, r_{3}\right)=\mathbf{r} with nodal function (1), then it has a twin sample r^(')=(r_(1)^('),r_(2)^('),r_(3)^('))=(r_(1),r_(2),r_(3))\mathbf{r}^{\prime}=\left(r_{1}{ }^{\prime}, r_{2}{ }^{\prime}, r_{3}{ }^{\prime}\right)=\left(r_{1}, r_{2}, r_{3}\right) with the nodal function (2). 原始样本。例如，如果一个 Ar-rigid H_(2)O\mathrm{H}_{2} \mathrm{O} 样本由核间距离 (Ar-O,Ar-H^(1),Ar-H^(2))=\left(\mathrm{Ar}-\mathrm{O}, \mathrm{Ar}-\mathrm{H}^{1}, \mathrm{Ar}-\mathrm{H}^{2}\right)=(r_(1),r_(2),r_(3))=r\left(r_{1}, r_{2}, r_{3}\right)=\mathbf{r} 与节点函数 (1) 给出，那么它就有一个节点函数 (2) 的孪生样本 r^(')=(r_(1)^('),r_(2)^('),r_(3)^('))=(r_(1),r_(2),r_(3))\mathbf{r}^{\prime}=\left(r_{1}{ }^{\prime}, r_{2}{ }^{\prime}, r_{3}{ }^{\prime}\right)=\left(r_{1}, r_{2}, r_{3}\right) 。