Gentlest ascent dynamics on manifolds defined by adaptively sampled point-clouds

The problem of finding saddle points of dynamical systems has one of its most notable applications in the search for transition states of chemical systems described at the atomistic level, since saddle points coincide with transition states ³ at the zero temperature limit. While finding local stable equilibria (sinks) is a relatively straightforward matter, finding saddle points is a more complicated endeavor for which a number of algorithms have been presented in the literature ⁴.

Saddle point search methods can be classified according to whether they require one single input state (usually the reactant, located at the minimum of a free energy well) or two states (reactant and product). The Gentlest Ascent Dynamics (GAD) ¹ belongs to the class of methods requiring a single reactant state as input and its applicability has been demonstrated with atomistic chemical systems ^{5, 6}. GAD can be regarded as a variant of the dimer method that is formulated as a continuous dynamical system whose integral curves with initial condition at the reactant state can lead to saddle points. Variants of GAD such as high-index saddle dynamics (HiSD) have been the subject of recent research efforts and applications ^{7, 8, 9, 2, 10, 11, 12, 13}.

While many search schemes attempt to find an optimal path between reactant and product (or between reactant state and transition state), it is interesting that there exist continuous curves joining the desired states in a variety of ways: following isoclines ^{14, 15}, gradient extremals ^{16, 17}, and the GAD studied in this paper. In most cases the study of these curves has been carried out in the Euclidean setting with some exceptions on the manifold of internal coordinates ^{18, 6} and on manifolds defined by the zeros of smooth maps ² (in these cases, the algorithms are formulated extrinsically on the ambient space). Our contribution is formulated intrinsically and is valid on arbitrary manifolds, not necessarily explicitly defined by an atlas or by the zeros of maps. Importantly, algorithms like the one presented here or our previous work ¹⁹ do not rely on a priori knowledge of good collective coordinates, but rather use manifold learning to find them on the fly. In our method, there is a feedback loop of data collection that drives progress towards a saddle point.

In this paper, we study an application of the GAD to manifolds defined by point-clouds. The manifold does not need to be characterized in advance either by the zeros of a smooth function or by an atlas, and it is only assumed that the user is capable of sampling the vicinity of arbitrary points on the manifold (e.g., umbrella sampling based on reduced local coordinates). The algorithm uses dimensionality reduction (namely, diffusion map coordinates ²⁰) to define a dynamical system intrinsically on the reduced coordinates that can lead to a saddle point. Since the saddle point in general is not expected to lie in the vicinity of the reactant, our algorithm works by iteratively sampling the manifold on the fly, resolving the path on the local chart, and repeatedly switching charts until convergence. Our approach shares algorithmic elements with our previous work ¹⁹ which, however, instead of GAD dynamics on manifolds, was following isoclines on manifolds.

Gentlest Ascent Dynamics (GAD) ¹ is an algorithm for finding saddle points of dynamical systems. We propose an extension of GAD to manifolds defined by point-clouds that finds saddle points combining nonlinear dimensionality reduction and adaptive sampling.

Let $U\colon\mathbb{R}^{n}\to\mathbb{R}$ be a smooth potential energy function and consider the associated gradient vector field $X$ in $\mathbb{R}^{n}$ given by $\dot{x}=X(x)$, where $X=-DU$ with $D$ denoting the gradient (or, equivalently, the Jacobian matrix). We restrict ourselves here, for the sake of simplicity, to the case of gradient systems and index-1 saddle points. The GAD algorithm consists of integrating the equations of motion of the related dynamical system $\hat{X}$ on an extended phase space $\mathbb{R}^{2n}$ given by

where $x,v\in\mathbb{R}^{n}$, $H(v)w=w-(2\,v\!\cdot\!w)v$ is the Householder reflection ²¹ of $w\in\mathbb{R}^{n}$ across the hyperplane $\langle v\rangle^{\perp}=\{z\in\mathbb{R}^{n}\mid v\cdot z=0\}$ and $r(x,v)=\|v\|^{-2}\,v\cdot D^{2}U(x)v$ is the Rayleigh quotient of the Hessian matrix $D^{2}U(x)$ corresponding to the vector $v\in\mathbb{R}^{n}$.

The intuition behind GAD is that the force $-DU$ can be written in a basis determined by the eigenvectors of the Hessian $D^{2}U$, which gives us the stable and unstable directions in the vicinity of an index-1 saddle point. The second differential equation in (1) acts as a (continuous) eigensolver yielding the unstable direction. Then, we can flip the sign of the component of the force corresponding to the unstable direction to make the resulting vector point towards the saddle point (see Figure 1).

Consider the manifold $M\subset\mathbb{R}^{n}$ and the gradient dynamical system $X\in\mathfrak{X}(M)$. We begin by drawing a total of $N$ samples from the manifold $M$ in the neighborhood of an initial point $p\in M$. This can be done in a variety of ways depending on the application. If $M$ is the inertial manifold of a dynamical system $X$, then a reasonable way to approach this problem is to generate $N$ distinct perturbations $\{p_{(i)}\in\mathbb{R}^{n}\mid i=1,\dotsc,N\}$ of $p$ and propagate them according to the flow of $X$ during a (short) time horizon $\tau>0$. Doing so, we obtain a data set $\mathscr{D}=\{q_{(i)}=\exp_{\tau}p_{(i)}\in M\mid i=1,\dotsc,N\}$ approximately on the manifold. Alternatively, one may numerically solve a stochastic differential equation such as the Brownian dynamics equation,

where the drift is the vector field $X$, $\sigma>0$ is a constant, and $B_{t}$ is a standard $n$-dimensional Brownian motion. Solving (3) (possibly with an added RMSD-based restraint around the initial conformation) up to a certain time $\tau>0$ and extracting an uncorrelated subset of the states at different time steps yields a data set $\mathscr{D}=\{q_{(i)}=q_{t_{i}}\mid i=1,\dotsc,N,0\leq t_{1}\leq\cdots\leq t_{N}\leq\tau\}$.

Next, we apply a dimensionality reduction algorithm on the data set $\mathscr{D}$ to obtain a set of reduced coordinates $\phi$. In our case, we use diffusion maps ^{20, 25} to obtain a set of vectors $\phi_{(i)}\in\mathbb{R}^{d}$ with $d\leq n$ but other methods, such as local tangent space alignment ³⁸, may be used as well. It is important to note that our dimensionality reduction method is applied to a local neighborhood of an initial point and, therefore, it is expected to yield a reasonable approximation to a chart on that neighborhood.

We fit a Gaussian process regressor $\phi$ to the pairs of points $(q_{(i)},\phi_{(i)})\in\mathbb{R}^{n}\times\mathbb{R}^{d}$ to obtain a smooth map $\phi\colon M\subset\mathbb{R}^{n}\to\mathbb{R}^{d}$ that will act as a system of coordinates (in particular, $\phi(q_{(i)})\approx\phi_{(i)}$). Proceeding in an analogous fashion, we compute the inverse mapping $\psi=\phi^{-1}$.

Note that one possible way of estimating the dimensionality is by computing the average of the approximate rank of the Jacobian matrix of $\phi$ at (a subset of) the data points and retaining the components that yield a local chart.

We consider the values $X_{q_{(i)}}\in TM$ of the vector field $X$ at the points in the data set $\mathscr{D}$ and map them via the system of coordinates in order to obtain the vector field $\phi_{\star}X_{q_{(i)}}=D\phi(q_{(i)})X_{q_{(i)}}$ in the new coordinates, where $D\phi$ is the Jacobian matrix of $\phi$ (note that the Jacobian-vector product can be computed either as a closed-form formula or efficiently using automatic differentiation —e.g., using jvp in JAX ³⁹). We fit another Gaussian process to the pushforward vector field $\phi_{\star}X$ in order to be able to evaluate it at arbitrary points in the new coordinates and we abuse notation in what follows by also referring to $\phi_{\star}X$ as $X$.

At this stage, we can readily compute the Riemannian metric $g$ as

where $g_{ij}=D\psi^{\top}D\psi$ for $i,j=1,\dotsc,d$. The metric induces an inner product, denoted by the bracket $\langle\cdot,\cdot\rangle$, between tangent vectors such that if $S=\sum_{i=1}^{d}S^{i}\frac{\partial}{\partial\phi^{i}}$ and $T=\sum_{i=1}^{d}T^{i}\frac{\partial}{\partial\phi^{i}}$ are the expressions in local coordinates of two tangent vectors $S,T\in T_{u}M$ at a point $u\in M$, then $\langle S,T\rangle=\sum_{i,j=1}^{d}g_{ij}S^{i}T^{i}$.

Using (4), we obtain the coefficients of the Levi-Civita connection ²³,

for $j,k,\ell\in\{1,\dotsc,d\}$, where $g^{ij}$ denotes the entries of the inverse of the matrix with components $g_{ij}$. This, in turn, allows us to take the covariant derivative and the Hessian, defined by

for arbitrary tangent vectors $S,T$. Observe that the Hessian is defined on the local chart, not on the ambient space. The eigenvector $V$, normalized with respect to $g$, with smallest eigenvalue of the Hessian then yields a vector field

such that the index-1 saddle points of $X$ become stable equilibria ¹ of $\hat{X}$. Consequently, an integral curve of $\hat{X}$ in the vicinity of a saddle point leads to the saddle point.

In general, in order to carry out the computation until convergence, we must frequently switch charts. This is due to the fact that at each step, we sample a point-cloud $\mathscr{D}$ in a small neighborhood of a given point and the Gaussian process regressor $\phi$ yields an approximation to $X$ on the chart that cannot extrapolate far away from the sampled points. Therefore, when we reach the confines of $\mathscr{D}$ (which can be determined by the density of points), we ought to map the latest point of our integral curve of $\hat{X}$ from the chart back to the ambient space using $\psi=\phi^{-1}$, as discussed earlier. After that, we start the whole procedure again: sample a new data set $\mathscr{D}$, compute $\phi$, trace an integral curve of $\hat{X}$, etc.

The factor that impacts the algorithm the most is the deviation of the sampled points in the data set from the manifold $M$. In other words, the farther away from $M$ our data points are, the more noisy is the estimate of $\hat{X}$, and, consequently, the harder it is to find the saddle points.

We have presented a formulation of GAD on manifolds defined by point-clouds that are meant to be sampled on-demand. Our formulation is intrinsic and does not require the specification of the manifolds either by a given atlas or by the zeros of a smooth map. The required charts are discovered through a data-driven, iterative process that only requires knowledge of an initial conformation of the reactant. We illustrated our approach with two simple examples and we expect the results to transfer to the high-dimensional dynamical systems of interest in computational statistical mechanics.

This work was supported by the US Air Force Office of Scientific Research (AFOSR) and the US Department of Energy DOE with IIT: SA22-0052-S001 and AFOSR-MURI: FA9550-21-1-0317.