Gaussian process regression for transition state search

The idea of GPR is described at length in the literature⁸ and we will only briefly review the basic idea. One can build a surrogate model for the PES using $N$ energies $E_{1},E_{2},...,E_{N}$ at certain configurations of the molecule $\vec{x}_{1},\vec{x}_{2},...,\vec{x}_{N}\in\mathbb{R}^{d}$. These configurations are called training points. In Cartesian coordinates the dimension $d$ of the system is $d=3n$, while $n$ is the number of atoms in the system. The key element of the GPR scheme is the covariance function $k(\vec{x}_{i},\vec{x}_{j})$. The covariance function describes the covariance between two random variables at the points $\vec{x}_{i}$ and $\vec{x}_{j}$ that take on the value of the energies. In simplified terms, the covariance function is a similarity measure between these two energies. In our case we use a form of the Matérn covariance function³⁰

in which we abbreviated $r=|\vec{x}_{i}-\vec{x}_{j}|$. The parameter $l$ describes a characteristic length scale of the Gaussian process (GP). It determines how strongly the covariance between two random variables (describing energies) decreases with distance.

Given a prior estimate $E_{\text{prior}}(\vec{x})$ of the PES, before we have included the training points in the interpolation, one can build the GPR-PES as follows.

The vector $\vec{w}=(w_{1}\ w_{2}\ ...\ w_{N})^{T}$ is the solution of the linear system

for all $m=1,2,...,N$. The elements

are the entries of the so called covariance matrix $K$ in which $\delta_{mn}$ is the Kronecker delta. The parameter $\sigma_{\text{e}}$ describes a noise that one assumes on the energy data. If one includes gradients in the scheme, one can introduce an additional parameter $\sigma_{\text{g}}$ for the noise on the gradient data.

Since the kernel function is the only dependency of $x$, we can obtain gradients and Hessians of the GPR-PES. For the first derivative in dimension $k$, i.e. in the direction of the $k$-th unit vector we get

and similarly

for second derivatives, in dimensions $k$ and $l$. In a previous paper we explicitly describe how one can include gradient information into this scheme.²⁰ In this paper we always use energy and gradient information at the training points. In order to build the GPR-PES we then have to solve a linear system with a covariance matrix of size $N(d+1)$. We solve this linear system via Cholesky decomposition. This yields a scaling of $\mathcal{O}\left(N^{3}d^{3}\right)$. In our case we can decrease the formal scaling to $\mathcal{O}\left(d^{3}\right)$ with a multi-level approach described below.

In our optimization algorithm  we take a similar approach as the P-RFO optimizer and built up a surrogate model for the PES. In contrast to P-RFO we only use energy and gradient information. We do not need any second derivatives. We use GPR to build the surrogate model by interpolating the energy and gradient information we obtain along the optimization procedure. On the resulting GPR-PES we can cheaply obtain the Hessian and therefore, also perform a P-RFO optimization. Our algorithm can be seen as a GPR-based extension to P-RFO to dispense the need for Hessian evaluations.  The result is used to predict a SP on the real PES. We perform all optimizations in Cartesian coordinates.

We call the vector that points from the last estimate of the TS to the current estimate of the TS the step vector $\vec{s}$, while the gradient at the current estimate of the TS is referred to as $\vec{g}$. DL-FIND uses a combination of convergence criteria, given a single tolerance value, $\delta$. The highest absolute value of all entries and the Euclidean norm of both, $\vec{g}$ and $\vec{s}$, have to drop below certain thresholds:

In these equations $d$ stands for the number of dimensions in the system.  If all these criteria are fulfilled, the TS is considered to be converged. These are the same criteria that are used for other TS optimizers in DL-FIND.

We chose a length scale of $l=20$ during all optimization runs. The noise parameters were chosen as $\sigma_{\text{e}}=\sigma_{\text{g}}=10^{-7}$. Despite accounting for numerical errors in the electronic structure calculations they also function as regularization parameters to guarantee convergence of the solution in Equation (3). We chose the noise parameters as small as possible without compromising the stability of the system. We generally found small values for the noise parameters to work better for almost all systems.  As the prior mean, $E_{\text{prior}}(\vec{x})$, of Equation (2), we set the mean value of all training points.

The parameter $s_{\text{max}}$, the maximum step size, must be specified by the user.

We start the optimization at an initial guess $\vec{x}^{\text{trans}}_{0}$ for the SP. At this point we obtain an approximate Hessian from a GPR-PES constructed from gradient calculations. This is done in such a way that we try to converge the eigenvector to the smallest eigenvalue of the Hessian, the transition mode $\vec{v}$. An estimate of that transition mode at the point $\vec{x}^{\text{trans}}_{0}$ is found according to the following procedure, which is the equivalent of dimer rotations in the dimer method.

1. 1.

   Calculate the energy and the gradient at the point $\vec{x}^{\text{trans}}_{0}$ and also at the point

   $$\vec{x}^{\text{rot}}_{1}=\vec{x}^{\text{trans}}_{0}+\frac{\Delta}{|\vec{v}_{0}|}\vec{v}_{0}$$

   (12)

   with $\vec{v}_{0}$ arbitrarily chosen. We generally chose $\vec{v}_{0}=(1\ 1\ ...\ 1)^{T}$. The results are included as training points in a new GPR-PES. We set $\Delta=0.1$ for our optimizer. Let $i=1$.

2. 2.

   Evaluate the Hessian $H_{i}(x^{\text{trans}}_{0})$ of the resulting GPR-PES at the point $\vec{x}^{\text{trans}}_{0}$.

3. 3.

   Compute the smallest eigenvalue of $H_{i}$ and the corresponding eigenvector $\vec{v}_{i}$.

4. 4.

   As soon as

   $$\left\lvert\frac{(\vec{v}_{i}\cdot\vec{v}_{i-1})}{|\vec{v}_{i}|\cdot|\vec{v}_{i-1}|}\right\rvert>1-\delta_{\text{rtol}}$$

   (13)

   we assume that the transition mode is converged, this procedure is terminated, and we move $\vec{x}^{\text{trans}}_{0}$ on the PES as we describe in the next section.

5. 5.

   If the transition mode is not converged, calculate the energy and gradient at the point

   $$\vec{x}^{\text{rot}}_{i+1}=\vec{x}^{\text{trans}}_{0}+\frac{\Delta}{|\vec{v}_{i}|}\vec{v}_{i}$$

   (14)

   and include the results to build a new GPR-PES. Increment $i$ by one and go back to step $2$.

The procedure to converge the transition mode is not repeated after each movement of $\vec{x}^{\text{trans}}$ on the PES but only after $50$ such moves. In these cases we start at the second step of the described procedure with the evaluation of the Hessian at the respective point. The initial guess for the transition mode after $50$ steps is the vector $\vec{v}_{i}$ of the last optimization of the transition mode. We use $\delta_{\text{rtol}}=10^{-4}$ (an angle smaller than $0.81\mathrm{°}$) at the very first point $x_{0}$ and $\delta_{\text{rtol}}=10^{-3}$ (an angle smaller than $2.56\mathrm{°}$) at every subsequent point at which we want to converge the transition mode.

In agreement with the minimum-mode following methods we assume that we have enough Hessian-information available to move on the PES to the SP as soon as the transition mode is found. We call the points that are a result from these movements on the PES $\vec{x}^{\text{trans}}_{i}$, with $i=0,1,...$. The points $\vec{x}^{\text{trans}}_{i}$ correspond to the midpoints in the dimer method. We use a user-defined parameter $s_{\text{max}}$ to limit the step size along the optimization. The step size can never be larger than $s_{\text{max}}$. The steps on the PES, starting at the point $\vec{x}^{\text{trans}}_{j}$, are performed according to the following procedure.

1. 1.

   Find the SP on the GPR-PES using a P-RFO optimizer. This optimization on the GPR-PES is stopped if one of the following criteria is fulfilled.

   • •

     The step size of this optimization is below $\delta_{\text{max}(s)}/50$.

   • •

     We found a negative eigenvalue (smaller than $-10^{-10}$) and the highest absolute value of all entries of the gradient on the GPR-PES drops below $\delta_{\text{max}(g)}/100$.

   • •

     The Euclidean distance between the currently estimated TS and the starting point of the P-RFO optimization is larger than $2s_{\text{max}}$.

   If none of these are fulfilled after $100$ P-RFO steps, we use a simple dimer translation. The dimer translation is done until one of the criteria above is fulfilled or the Euclidean distance between the currently estimated TS and the starting point is larger than $s_{\text{max}}$. The P-RFO method converged in less than $100$ steps in all of the presented test cases of this paper.

2. 2.

   Overshooting the estimated step, according to the overshooting procedures described in the next section, resulting in $\vec{x}^{\text{trans}}_{j+1}$.

3. 3.

   Calculate the energy and gradient at $\vec{x}^{\text{trans}}_{j+1}$ and build a new GPR-PES.

The used P-RFO implementation for this procedure is the same as used in DL-FIND (with Hessians of the GPR-PES re-calculated in each step). P-RFO tries to find the modes corresponding to overall translation and rotation of the system and ignores them. Numerically, it is not always clear which modes these are. Therefore, we found it to be very beneficial to project the translational modes out of the inferred Hessian. This yields translational eigenvalues that are numerically zero. Otherwise, the TS search via P-RFO tends to translate the whole molecule, leading to an unnecessary large step size.

As we have done in the optimization algorithm for minimum search,²⁰ we try to shift the area in which we predict the SP to an interpolation regime rather than an extrapolation regime. This is because machine learning methods perform poorly in extrapolation. The overshooting is done, dependent on the angle between the vectors along the optimization: let $\vec{s}_{N}^{\,\prime}$ be the vector from the point $\vec{x}^{\text{trans}}_{N-1}$ to the next estimate for the TS according to the first step in the procedure described in the previous section. Let $\vec{s}_{N-1}$ be the vector pointing from $\vec{x}^{\text{trans}}_{N-2}$ to $\vec{x}^{\text{trans}}_{N-1}$. If $\vec{x}^{\text{trans}}_{N-1}=\vec{x}_{0}$, the point $\vec{x}^{\text{trans}}_{N}$ is simply calculated according to the procedure described in the previous section with no overshooting. We calculate an angle

and introduce a scaling factor

that scales

as soon as $\alpha_{N}>0.9$. The scaling limit $\lambda_{\text{max}}$ is chosen to be $5$. But it is reduced if the algorithm is close to convergence and it is increased if $\alpha_{N}>0.9$ for two consecutive steps, i.e. we overshoot more than once in a row. We refer to our previous work for a more detailed description of the calculation of $\lambda_{\text{max}}$.²⁰ We also use the separate dimension overshooting procedure from that paper: if the value of a coordinate along the steps $\vec{x}^{\text{trans}}_{i}$ changed monotonically for the last $20$ steps, a one dimensional GP is used to represent the development of this coordinate along the optimization. It is optimized, independent from the other coordinates. To account for coupling between the coordinates in succeeding optimization steps this procedure is suspended for $20$ steps after it was performed. After these two overshooting procedures the step is scaled down to $s_{\text{max}}$. The chosen overshooting procedures might seem too aggressive at first glance. But the accuracy of the GPR-PES is largely improved if we overshoot the real TS. This is a crucial difference to conventional optimizers: a step that is too large/in the wrong direction can still be used to improve the next estimate of the TS.

Just like in our previous work²⁰ we include a multi-level scheme to reduce the computational effort of the GPR. A more detailed explanation of the multi-level approach can also be found there. The most demanding computational step in the GPR scheme is solving Equation (3). The computational effort to solve this system can be reduced by limiting the size of the covariance matrix. Our multi-level scheme achieves this by solving several smaller systems instead of one single large system. This eliminates  the formal scaling with the number of steps in the optimization history. To achieve this, we take the oldest $m$ training points in the GP to build a separate GP called $GP_{1}$. This is done as soon as the number of training points reaches $N_{\text{max}}$. The predicted GPR-PES from $GP_{1}$ is used as $E_{\text{prior}}(\vec{x})$ in Equation (2) for the GP built with  the remaining $N_{\text{max}}-m$ training points that we call $GP_{0}$. The new training points are added to $GP_{0}$. When $GP_{0}$ eventually has more than $N_{\text{max}}$ training points we rename $GP_{1}$ to $GP_{2}$ and use the $m$ oldest training points in $GP_{0}$ to build a new $GP_{1}$. We always take $GP_{i+1}$ as the prior to $GP_{i}$. The $GP_{q}$ with the highest index $q$ just uses the mean of all energy values at the contained training points as its prior. The number of levels increases along the optimization but since the number of points in all $GP_{i}$ is kept below a certain constant the formal scaling of our GPR scheme is $\mathcal{O}\left(d^{3}\right)$. Splitting the points which were used to converge the smallest eigenvalue of the Hessian (according to the procedure for the converging of the transition mode described above) leads to a loss in accuracy of the second order information of the GPR-PES and consequently, to inaccurate predictions of the step direction by the underlying P-RFO.  Therefore, we increase $N_{\text{max}}$ and $m$ by one if those points would be split into different levels and try the splitting again after including the next training point. When we successfully performed the splitting, the original values of $N_{\text{max}}$ and $m$ are restored. We set the values $N_{\text{max}}=60$ and $m=10$ for all tests performed in this paper. We also want to guarantee that we always have information about the second derivative in $GP_{0}$. Therefore, we start our Hessian approximation via Procedure 1 after $N_{\text{max}}-m=50$ steps on the PES. As a result, the points that belong to the last Hessian approximation are always included in $GP_{0}$. The method described so far is referred to as GPRTS in the following.

If one intends to find a TS that connects two known minima, the algorithm tries to find a good guess for the starting point of the TS search automatically. To that end we optimize a nudged elastic band (NEB) ^{21, 22} in the image dependent pair potential (IDPP).²⁵ The construction of this potential does not need any electronic structure calculations, but is a good first guess for a NEB path based on geometrical considerations. Given $M$ images $\vec{x}^{\text{IDPP}}_{i},\ i=1,...,M$ of the optimized NEB in the IDPP, our algorithm calculates the energy and gradient of the real PES at image $\vec{x}^{\text{IDPP}}_{j}$ with $j={M/2}$ for even and $j=(M-1)/2$ for odd $M$. Only from this point on we perform energy evaluations of the real PES.  Then we calculate the scalar product between the gradient $\vec{\gamma}_{j}$ at point $\vec{x}^{\text{IDPP}}_{j}$ and the vector to the next image on the NEB path.

If $\beta>0$, we calculate the energy and the gradient at $\vec{x}^{\text{IDPP}}_{j+1}$, otherwise at $\vec{x}^{\text{IDPP}}_{j-1}$ in order to aim at a maximum in the energy. This is repeated until we change the direction along the NEB and would go back to a point $\vec{x}^{\text{IDPP}}_{\text{best}}$ at which we already know the energy and the gradient. The point $\vec{x}^{\text{IDPP}}_{\text{best}}$ we found in this procedure is assumed to be the best guess for the TS on the NEB and is the starting point for our optimizer. All the energy and gradient information acquired in this procedure is used to build the GPR-PES. The optimization technique including the usage of the NEB in the IDPP potential will be abbreviated as GPRPP.

We optimized the TSs of our test set using the new GPRTS optimizer and also the established dimer method and the P-RFO method in DL-FIND. We compare the number of energy evaluations that all three methods need until convergence in table 1. Note that analytic Hessians are available for systems $26$ and $27$. P-RFO uses $8$ ($4$) analytic Hessians for system $26$ ($27$). Those are not counted as energy evaluations. In the other systems the Hessians are calculated using central difference approximation via gradients.  The required $2d$ gradients are counted as energy evaluations in the table. Note that every energy evaluation in our tables always refers to evaluating the energy and the gradient of the PES.  We chose a maximum step size of $s_{\text{max}}=$0.3$$ and a tolerance $\delta=$3\text{\times}{10}^{-4}$$ for all optimizations.

The P-RFO method calculates the Hessian at the first point and then only every $50$ following steps. All other Hessians are inferred via the update mechanism by Bofill. ³⁸

We see that the GPRTS method generally requires fewer energy evaluations than the other methods and has the least convergence problems. In table 1 we also show the energy differences (the energy from the respective method minus the energy from GPRTS) and root mean squared deviations (RMSD) of the found TSs, compared to the TSs found by GPRTS. Looking at the specific systems with relevant deviations for the found TSs we observe the following differences, compare Fig. 1:

• •

  System 3: P-RFO finds a different TS in which $\ce{H2}$ is completely abstracted. Its energy is $9.34\times 10^{-2}$ Hartree lower than that of the TS found by GPRTS.

• •

  System 4: GPRTS and P-RFO yield a smaller distance of the $\ce{H}$-Atom to the remaining $\ce{O}$ atom. which corresponds to a different transition than intended: a rotation of the $\ce{OH}$ group.

• •

  System 6: GPRTS finds a different angle H–C–C at the carbene end.

• •

  System 7: P-RFO finds an opening of the ring structure which remains closed in the other cases.

• •

  System 10: GPRTS finds a TS in that $\ce{N2}$ is abstracted, the dimer method finds the depicted opening of the ring. The opening of the ring is also found via GPRTS if one uses a smaller step size. P-RFO finds a closed ring structure as the TS.

• •

  System 11: P-RFO finds a more planar structure.

• •

  System 14: The dimer method finds a mirror image of the TS found by the other methods, thus the high RMSD.

• •

  System 20: P-RFO finds a different orientation of $\ce{NH3}$, rotated relative to $\ce{HCO}$, and also with a slightly different distance.

• •

  System 22: P-RFO finds a different angle of the attached $\ce{OH}$ group. The resulting molecule is not planar. The resulting structures of GPRTS and the dimer method are planar.

• •

  System 26: P-RFO does not find the $\ce{H}$-transfer but an opening of the ring.

Furthermore, we also compare the evolution of the residual error of GPRTS, P-RFO, and the dimer method during the TS search for the four largest test systems in the Baker test set in Fig. 3 and for systems $26$ and $27$ in Fig. 4. In both plots we only count the number of steps on the PES that are taken along the optimization runs. The energy evaluations performed to converge the transition mode, to rotate the dimer, and to calculate the Hessian in P-RFO are excluded. We see a significant speedup with the new GPRTS optimizer. In addition, the new method needs also fewer gradients than the dimer method to approximate the Hessian appropriately. In our optimizer we often see jumps of the gradient norm at around $50$ steps. This is where the multi-level approach creates the first separate GP. Therefore, we include a possibility to increase the number of training points in the last level $GP_{0}$ in the code to get better performance in high-precision calculations. We do not make use of this possibility in all test systems presented in the paper.

As described in the previous section about finding a starting point via IDPP we use an approximated NEB path between two known minima to find a good starting point for the TS search. To prepare the input for our test systems we proceeded as follows:

• •

  We start from the TSs shown in Fig. 1 and Fig. 2.

• •

  We performed IRC path searches to find the two minima connected by the TS.

• •

  The geometries of these minima are used as starting points for GPRPP.

• •

  Additionally, GPRTS, dimer and P-RFO searches are started from $\vec{x}^{\text{IDPP}}_{\text{best}}$. In GPRPP, the points on the NEB path are included in the GPR, while in GPRTS they are excluded.

In the benchmarks reported in table 2 the additional points to find $\vec{x}^{\text{IDPP}}_{\text{best}}$ are included in the number of steps that the GPRPP optimizer needs, but not in the numbers for the other methods. Our approximated NEB path consists of $10$ images for all test cases. The number of additional energy evaluations that GPRPP needs to find $\vec{x}^{\text{IDPP}}_{\text{best}}$ is $2$ or $3$ for all test systems. We set the same maximum step size and tolerance value as in the previous section. We compare the number of energy evaluations all methods need until convergence in table 2. We also show the energy differences (the energy from the respective method minus the energy from GPRPP) and RMSD values compared to GPRPP. Note that analytic Hessians are available for systems $26$ and $27$. P-RFO uses $2$  analytic Hessians for system $26$ and $5$ for system $27$.  Those are not counted as energy evaluations. In the other systems the Hessians are calculated via finite differences and the required $2d$ gradients are counted as energy evaluations in the table. We observe that GPRPP and GPRTS yield consistently good results. Comparing GPRPP and GPRTS we see no clear tendency whether it is useful to include the points that lead to $\vec{v}_{\text{best}}$ in order to improve the speed of the convergence. On the other hand, we have one system that only converges when these points are included. We conclude that the stability might be improved using the GPRPP method with the additional training points in comparison to pure GPRTS. This will strongly depend on the choice of the starting point for GPRTS.