Efficient discovery of multiple minimum action pathways using Gaussian process

Finding multiple transition pathways of phase transitions, chemical reactions, and conformational transitions remains challenging in physics and chemistry. With conventional molecular dynamics (MD) or Monte-Carlo (MC) approaches, i.e., initial-value formulation, it takes a considerably long time to simulate an evolution over energy barriers starting from an initial state. To overcome such limitations, numerous single-end methods have been suggested [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 12, 13, 14, 15, 16]. Methods include using Hessian to climb up the energy barrier [1, 2], using the eigenvalues of a Hessian matrix [3, 4], repeating the process of climbing up and relaxing to locate multiple local minima [5], and using the knowledge, which atoms are more likely to form or to break a bond, to guide the transition [12, 13, 14, 16], were suggested. The referred methods find the transition state efficiently. However, these methods do not guarantee finding a transition pathway reaching a targeted product state.

Another class of approaches, biased sampling methods, were introduced to sample a broad range of configurational space [8, 9, 10, 11]. By applying additional potential on the region already sampled, it forces a system to explore over energy barriers [8, 9]. They efficiently estimate the probabilities of sampled local minima or observable properties of an ensemble. However, they cannot provide the probability estimates of multiple transition pathways and associated kinetic information.

To overcome the limitations of single-end methods, various double-end, or fixed-boundary-value, sampling approaches have been developed [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. The fixed-boundary-value formulation guarantees to find virtual trajectories connecting two end states if a calculation converges successfully, while the initial-value formulation does not. Berkowitz et al. suggested variational formula for the system with frictional resistance [17], Elber and Karplus suggested the Gaussian chain approach that minimizes the line integration of a potential along a trajectory [19]. Czerminski and Elber added a repulsion constraint to the Gaussian chain approach to prevent replicas from aggregation [20]. To accelerate the convergence, the nudged elastic band (NEB) method was suggested [24]. NEB applies an artificial force to the tangential component of a trajectory chain to remove the corner-cutting problem. These methods find the minimum energy pathway (MEP) with the lowest energy barrier along the transition trajectory. However, the concept of a minimum energy pathway is not well-defined when multiple energy barriers exist between end states on a highly rugged energy landscape.

Contrary to MEP-based approaches, pathway sampling methods based on the least action principles were proposed [17, 21, 22, 23, 25, 37, 29, 30, 31, 38, 34, 35, 36]. The variational Verlet algorithm [21] and the Fourier component relaxation method that minimize the classical action of the trajectory expressed as Fourier (sine) components [22] were suggested. Olender and Elber suggested using Onsager Machlup action to investigate the long-time molecular trajectory [23].

Passerone and Parrinello [25] developed an action-derived molecular dynamics (ADMD) method, which finds low-action pathways via local optimization of a modified classical action with an energy conservation restraint term to keep the system’s total energy constant throughout a reaction. Lee and coworker [29] developed a kinetic controlled ADMD method introducing a restraining term to the modified action that satisfies the equipartition theorem. However, these classical-action-based methods have an inherent limitation: the classical principle of least action is an extremal, not a minimum, principle. In other words, the classical action can be either minimized or maximized, which makes its computational outcome ambiguous.

As a way to overcome this limitation of classical-action-based methods, pathway sampling approaches based on the Onsager-Machlup (OM) action [39, 40, 41] have been suggested. the probability $p(\mathbf{x}_{B}\|\mathbf{x}_{A};t)$ for a state $\mathbf{x}_{A}$ to become another state $\mathbf{x}_{B}$ within a given time interval $t$, can be described as follows [30, 39, 40, 41]:

where $S_{\mathrm{OM}}$ is the OM action and $\mathcal{D}$ is the notation for integrating all the possible pathways connecting two end states. The probability of a transition from the state $\mathbf{x}_{A}$ to $\mathbf{x}_{B}$ is calculated by summing up the probability of all pathways $\sum_{s}e^{-S_{\mathrm{OM}}\left(\mathbf{x}\left(\tau\right)\right)/k_{\mathrm{B}}T}$ [39, 40, 42, 43]. For most cases, a state transition usually incorporates multiple pathways, such as multiple protein folding pathways. The most dominant pathway is determined by the path with the lowest OM action, corresponding to the maximum probability of a transition [42, 37, 30, 41, 44, 31, 36, 35]. Lee and coworkers utilized an efficient global optimization method to identify the most dominant transition pathway corresponding to the global minimum of OM action and sample multiple less prevalent pathways [36].

A major bottleneck of the existing pathway search methods is their computational cost. Calculating the action of a pathway requires the energy and gradient of every replica of a pathway. A small step size between replicas is necessary to generate smooth pathways. Thus, there is a trade-off between trajectories’ accuracy and computational efficiency.

To alleviate this computational cost issue, Jónsson and coworkers [45, 46] suggested the one-image evaluation (OIE) method by combining the Gaussian process (GP)  [47, 48] with NEB [24]. The OIE method finds the MEP by building a Gaussian potential energy surface (PES) approximating the actual PES iteratively. They reduced the number of function calls reduced by two orders of magnitude on an artificial 2D PES [46].

Another major limitation of existing pathway search methods, such as NEB, is that they are local optimization methods whose final results heavily depend on initial guesses on pathways. Generally, conventional pathway search methods start from the linear interpolation between end states and perform local optimization of the initial pathway. As a result, these approaches can find the correct minimum energy or action pathways only when a given PES is simple and smooth. However, many complex systems have highly rugged PESs, which do not guarantee to locate the most dominant pathway near the linear interpolation between two states. An extensive search on a pathway space is essential to avoid this problem.

To overcome these limitations of the pathway search methods, we introduce a novel transition-pathway-sampling method, Gaussian process action optimization (GPAO), by combining the GP algorithm with a global action optimization approach. Inspired by the OIE method, we only calculate the most uncertain images measured by the variance obtained with GP. Then the energy and forces of the images are used to construct the Gaussian PES. Afterward, the intermediate images’ potential energies, forces, and Hessian are estimated using an updated Gaussian PES. Pathways are optimized to have the minimum action using the predicted potential, forces, and Hessian. Our method conducts direct OM action optimization with a total-energy-conservation restraint without using the original system’s Hessian.

The most significant advantage of GPAO is a dramatic reduction of the number of force calls by four-to-five orders of magnitude than the conventional action optimization method while preserving their accuracy. The significant reduction of the computational cost by GPAO allows us to find multiple transition pathways of a molecule—here, we choose alanine dipeptide—with a quantum-mechanical (QM) potential. Additionally, we investigated the isomerization pathways of small organic molecules and the rearrangement pathways of atoms on a metallic surface. Throughout these examples, it is clearly identified that GPAO successfully identifies multiple physically accessible trajectories between the two end-states using QM potential energy. To the best of our knowledge, this work is the first report, which identifies multiple conformational transition pathways of alanine dipeptide with an ab initio potential.

The Müller-Brown (MB) potential is a 2-dimensional model potential used to benchmark double-end methods[49]. To benchmark the efficiency of GPAO, we performed three calculations using GPAO, direct action optimization (DAO) without using a Gaussian potential, and the NEB method on the MB potential and compared their efficiency and accuracy. The number of intermediate images, range of hyperparameters, and Gaussian kernel are set the same as previous studies [25, 46]. The MB potential is defined as follows:

where $x_{\mu}^{0}=\left(1,\ 0,\ -0.5,\ -1\right)$, $y_{\mu}^{0}=\left(0,\ 0.5,\ 1.5,\ 1\right)$, $a_{\mu}=\left(-1,\ -1,\ -6.5,\ 0.7\right)$, $b_{\mu}=\left(0,\ 0,\ 11,\ 0.6\right)$, $c_{\mu}=\left(-10,\ -10,\ -6.5,\ 0.7\right)$, and $A_{\mu}=\left(-2,-1,-1.7,\ 0.15\right)$.

The process of how GPAO minimizes the MEP’s uncertainty and approximates the MB potential is illustrated in Figure 1. At iteration 1, the pathway is optimized with the three data points: the initial, the final, and a random intermediate state between the two. Due to the lack of data points, potential energy values estimated along the trajectory are highly uncertain. The uncertainty range (light blue intervals) of the estimation at iteration 1 is much wider than iteration 9 (Figure 1b). The main idea is to calculate the energy and forces of the point that minimizes the uncertainty of a trajectory most. Thus, at iteration 2, GPAO collects data at the most uncertain region of the previous iteration (marked with a red cross in Figure 1b). After the update, the OM action of a pathway is minimized on the updated Gaussian PES, whose hyperparameters are optimized with new data points. After hyperparameter optimization, at iteration 3, the energy and forces of the most uncertain points in the pathway are added to the database (red crosses in Figure 1b), and the process is iterated. At iteration 9, the pathway is considered as converged since the uncertainty criteria, $\mathbf{\Sigma}^{\mathrm{(max)}}<0.05$, are satisfied. In total, GPAO requires 11 force calls.

The targeted total energy of the system, $E_{\mathrm{t}}$, one of the hyperparameters of the method, is adjusted during iterations based on the following equation:

where $E_{\mathrm{t}}^{(i)}$ is the target energy of $i$th iteration, $E_{\mathrm{t}}^{(i-1)}$ is the target energy of iteration $i-1$, and $\boldsymbol{\mu}^{\mathrm{(max)}}$ is the estimated maximum potential energy of iteration $i-1$. Figure 1(c) shows how GPAO approximates $E_{\mathrm{t}}$, iteratively. At iteration 1, $E_{\mathrm{t}}$ is set to $-1.47$, the minimum energy between the initial and the final image. GPAO sets $\boldsymbol{\mu}^{\mathrm{(max)}}=-0.10$ from the trajectory optimized with three data points. Using this information, at iteration 2, $E_{\mathrm{t}}$ is raised to $-0.81$, which is the half of $E_{\mathrm{t}}+\boldsymbol{\mu}^{\mathrm{(max)}}$ of iteration 1. Thus, GPAO estimates $\boldsymbol{\mu}^{\mathrm{(max)}}$ to be $-0.20$ with the four data points. Similarly, at iteration 3, $E_{\mathrm{t}}$ is set to the half of $E_{\mathrm{t}}+\boldsymbol{\mu}^{\mathrm{(max)}}$ of iteration 2, $-0.368$.

Using the target energy obtained from GPAO, we also obtain pathways using DAO with three actions: classical action with total energy conservation restraint $S^{\mathrm{cls}}$ (blue dashed line in Figure 2), OM action $S^{\mathrm{OM}}$ (orange dash-dotted line in Figure 2), and OM action with total energy conservation restraint ${\Theta}^{\mathrm{OM}}$ (solid green line in Figure 2). The trajectories obtained with direct optimization of three actions are denoted as $\mathcal{C}_{\Theta}^{\mathrm{cls}}$, $\mathcal{C}_{S}^{\mathrm{OM}}$, and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$, and the trajectories obtained with GPAO are denoted as GP-$\mathcal{C}_{\Theta}^{\mathrm{cls}}$, GP-$\mathcal{C}_{S}^{\mathrm{OM}}$, and GP-$\mathcal{C}_{\Theta}^{\mathrm{OM}}$.

The most remarkable result is that GPAO requires significantly fewer force calls than DAO (Table 1). GP-$\mathcal{C}_{\Theta}^{\mathrm{cls}}$ requires 12 force calls, while $\mathcal{C}_{\Theta}^{\mathrm{cls}}$ requires around 1,200,000 force calls. Similarly, GP-$\mathcal{C}_{S}^{\mathrm{OM}}$ and GP-$\mathcal{C}_{\Theta}^{\mathrm{OM}}$ require 14 and 12 force calls respectively, while $\mathcal{C}_{\Theta}^{\mathrm{cls}}$ and $\mathcal{C}_{\Theta}^{\mathrm{cls}}$ require 600,000 and 900,000 force calls respectively. Overall, GPAO reduced the number of force calls in four to five orders of magnitude.

The GPAO results are highly similar to those obtained with DAO with the significantly reduced numbers of force calls. Overall, the Frétchet distances between the pathways obtained with GPAO and DAO are less than 0.08% of the pathways’ total distances. The distance between $\mathcal{C}_{\Theta}^{\mathrm{cls}}$ and GP-$\mathcal{C}_{\Theta}^{\mathrm{cls}}$ is 0.006. Similarly, the distances between $\mathcal{C}_{S}^{\mathrm{OM}}$ and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ and their counterpart obtained with GPAO are 0.082 and 0.042, respectively. GP-$\mathcal{C}_{\Theta}^{\mathrm{cls}}$ deviates from $\mathcal{C}_{\Theta}^{\mathrm{cls}}$ less than a Fréchet distance of 0.01. All four GP-assisted pathways show a deviation of less than a Fréchet distance of 0.1 than the DAO results. These results clearly demonstrate that approximating PES with GP is accurate enough to locate the transition pathways with low action on a given PES.

The convergence of GPAO depends on the number of data points near the MEP, which should be large enough to approximate the PES accurately. Overall, GPAO requires similar numbers of force calls regardless of the selection of the objective function (Table 1). For example, the results of the GP-assisted NEB method, GP-$\mathcal{C}^{\mathrm{NEB}}$, showed similar numbers of force calls to reproduce the result of the conventional action optimization method. Indeed, the number of data points required to find an accurate MEP on the MB potential is also between 10 and 20. These results imply that the type of local relaxation method does not significantly affect the efficiency of GPAO. In contrast to GPAO, the calculation cost of DAO is strictly proportional to the number of images in a pathway.

After confirming that GPAO successfully optimizes a given action of a pathway, we compare the characteristic of $\mathcal{C}_{\Theta}^{\mathrm{cls}}$, $\mathcal{C}_{S}^{\mathrm{OM}}$, and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ (Figure 2a). The results demonstrate that $\Theta_{\mathrm{OM}}$ successfully conserves the system’s total energy throughout the pathway (Figure 2b). The magnitude of total energy fluctuations of $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ is reduced by 98.7% than that of $\mathcal{C}_{S}^{\mathrm{OM}}$. The action value of $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ increases only about 19%, and the Fréchet distance to $\mathcal{C}_{S}^{\mathrm{OM}}$ is less than 0.067 (Table 1). Also, 95% of the intermediate images of $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ are located within a distance of 0.013 from $\mathcal{C}_{S}^{\mathrm{OM}}$. This result demonstrates that the total energy conservation term affects little to the shape of the pathway while conserving the total energy of the system. The results indicate that enforcing the energy conservation restraint on $\mathcal{C}_{S}^{\mathrm{OM}}$ gives only little effect on finding low $V^{\mathrm{(max)}}$ and $S_{\mathrm{OM}}$.

Figure 2b shows the total energy profiles of three trajectories and their energy gaps, the energy differences between the maximum and the minimum total energies of pathways. The energy gaps of $\mathcal{C}_{\Theta}^{\mathrm{cls}}$, $\mathcal{C}_{S}^{\mathrm{OM}}$ and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ are 0.09, 3.85 and 0.05, respectively. The large spikes of the total energy are observed between the time steps 0.8 $\sim$ 1.1, where the trajectory climbs up the energy barriers of the PES. The main reason for this huge spike of the total energy is the second term of the $S_{\mathrm{OM}}$ equation (Eq. (7.2)). The first term in Eq. (7.2) prefers to minimize atoms’ forces, and the third term keeps the velocities between images low. The second term in Eq. (7.2) favors following a concave potential energy surface along its ongoing direction. More specifically, if a pathway has to climb up a PES during a transition, it prefers to climb up the surface fast to minimize action. This tendency is clearly manifested in Figure 2a. Right before $\mathcal{C}_{S}^{\mathrm{OM}}$ reaches the saddle point with the maximum potential energy of the PES, the images are largely separated because of the large velocities of corresponding images. In contrast, the images of $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ near the saddle point are almost continuous and uniformly distributed. Uniform distances between images prevent large fluctuations of total energy. In other words, GPAO avoids the occurrence of unphysical energy spikes observed with $\mathcal{C}_{S}^{\mathrm{OM}}$, but still maintains $S^{\mathrm{OM}}$ low. The action values of $\mathcal{C}_{\Theta}^{\mathrm{cls}}$, $\mathcal{C}_{S}^{\mathrm{OM}}$ and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ are 14.49, 1.387, and 1.652, respectively. The potential energies at the saddle point, $V^{\mathrm{(max)}}$, of $\mathcal{C}_{\Theta}^{\mathrm{cls}}$, $\mathcal{C}_{S}^{\mathrm{OM}}$, and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$ are $-0.400$, $-0.406$, and $-0.406$, respectively (Table 1). Our results demonstrate that unphysical increases of atomic velocities to lower the second term can be suppressed effectively with the total energy conservation restraint without changing trajectories. Furthermore, two pathways, $\mathcal{C}_{S}^{\mathrm{OM}}$ and $\mathcal{C}_{\Theta}^{\mathrm{OM}}$, overlap significantly each other.

To further validate the advantages of GPAO in application to real molecular systems, we benchmark the computational cost and characteristics of pathways obtained with DAO and GPAO using density functional theory (DFT) on the isomerization of small molecules: formaldehyde, formic acid, and propyne (Figure 3). We perform GPAO first, and the resulting target total energy is used for DAO calculations. Additionally, the final trajectory obtained from a GPAO calculation is used as the initial trajectory for DAO to reduce computational cost. In this way, the difference between trajectories approximated with GPAO and true ground state obtained with DAO is compared.

Atomic configuration of initial states (IS), transition state (TS), and final state (FS) of each molecule are drawn in Figure 3 with the arrows pointing to corresponding reaction coordinates. GPAO begins with the five data, two initial and final states and three points between them. During the process, the potential profile of the trajectory quickly converges as the data accumulate throughout iterations. In Figure. 3, the potential energy profiles of GPAO (solid blue line) and their variance (light blue shade) slowly converge to those of DAOs (solid orange line). In the rightmost graph in Figure 3, Gaussian PES constructed with this process are compared with true DFT data sampled along the DAOs trajectory (orange dot). In all three cases, using only, comparably a few data points sampled with the GPAOs process (black cross), root mean square error (RMSE) is less than $0.06$ eV, or $0.01$ eV/atom.

A comparison of the number of force calls, energy barrier, and actions obtained with GPAO and DAO is summarized in Table 2. The number of force calls for DAO of formaldehyde, formic acid, propyne are $3.52\times 10^{5}$, $3.70\times 10^{5}$, and $4.73\times 10^{5}$, respectively. Meanwhile, the numbers of force calls conducting GPAO counterpart are $18$, $30$, and $49$, respectively. Thus, the overall computational cost of GPAO is less than that of DAO by two to three orders of magnitude, while the error of OM action and energy barrier is suppressed under 5%. All the example shows consistent improvement in computational efficiency with comparable accuracy of GPAO.

One of the major challenges in performing DAO with DFT is an enormous computational cost for Hessian calculations. The Hessian matrix can be obtained analytically for simple potentials such as the MB potential. However, obtaining a Hessian matrix is not straightforward for a DFT potential in general. Thus, a finite difference method is used in this study. To obtain a Hessian matrix using the central difference method, $D$ additional force calls are needed, where $D$ is the degree of freedom. For example, for the Hessian calculations of formic acid and propyne, 12 and 18 additional force calls are required.

The Cartesian coordinates of atoms are used as the input for the Gaussian PES kernel. To satisfy the rotational and translational invariance of the GP process, the position of an arbitrary atom of a system is fixed. When an atom is not fixed, GPAO keeps searching the same geometry but slightly translated position, making the convergence of GPAO difficult.

In the third class of examples, the diffusion and rearrangement of atoms on a metallic surface are studied. The potential energies of metallic systems are calculated using effective medium theory (EMT) in the Atomic Simulation Environment (ASE) package [50]. Using the EMT potential, we compare the characteristics of DAO and GPAO results on two examples: ‘Au hopping’ and ‘Pt rearrangement’ [46].

In the ‘Au hopping’ example, a gold atom jumps to its neighboring site on the aluminum fcc (100) surface. To visualize the diffusive pathway, the atomic configurations at the key states are plotted in Figure 4. From left to right, the atomic configurations of the initial state (IS), transition state (TS), and the final state (FS) watching from the top and the side are illustrated. For both IS and FS, the gold atom is propped in the middle of four aluminum atoms. In TS, the gold atom is on the saddle point supported by two aluminum atoms. Two aluminum atoms right under the gold atom sink slightly, while the other two aluminum atoms, far from the gold atom, slightly protrude along the z-axis. This shift lowers the energy barrier to 0.4 eV. The GPAO trajectory successfully reconstructs the shift of potential energy barrier, leading to an almost identical trajectory to the DAO trajectory.

Below the atomic configuration, energy profiles at iterations 1, 3, and 10 (Iter: 1, 3, and 10 in Figure 4) are plotted. While optimizing action, only the gold atom and its neighboring four aluminum atoms are allowed to move and the other atoms are fixed reducing the degree of freedom $D=15$. At iteration 1, GPAO optimizes a trajectory on a Gaussian PES, which is constructed with five initial data points. The initial data are IS, FS, and additional three random points between the two end states. The energy profile at this stage highly deviates from that of DAO, and it has high uncertainty. As the calculation proceeds, the GPAO trajectory and energy profile converge to those of DAO. At iteration 10, the trajectory converges to the DAO trajectory with a covariance less than 0.05 eV. As a result, using only 13 data points, GPAO accurately predicted the saddle points. The EMT data and predicted potential with GP of the images of the DAO trajectory are compared. For 15 images, the RMSE of potential prediction is 70 meV, or 5 (meV/atom). Compare to the common neural network potential, the number of data required to find the energy barrier is relatively small. This implies that our iterative searching scheme is highly effective for investigating surface hopping problems.

In the ‘Pt rearrangement’ example, seven platinum atoms packed as a hexagonal shape on the fcc (111) surface are rearranged to a parallelogram shape. Among 71 atoms, seven atoms on the surface are set to move freely, setting the degree of the freedom as $D=21$. To lower the calculation cost for the Gaussian kernel, surface platinum atoms under the seven hexagon atoms are fixed. In Figure 5, from left to right, IS, TS, and FS of ‘Pt rearrangements’ are drawn emphasizing moving atoms with a blue edge. In IS and FS, seven platinum atoms on the surface are located at symmetry sites. In TS, the center atom is misplaced from its symmetric sites, while maintaining a two-fold symmetry. This lowers the energy barrier to 1.0 eV.

Below the atomic configurations in Figure 5, the energy profiles of GPAO at iterations 1, 9, and 186 (Iter:1, 9, and 186) are illustrated. At iteration 1, GPAO optimizes the trajectory from five data points: IS, FS, and three linearly interpolated points between the two. The energy barrier at this stage is nearly 2 eV, and the confidence for this prediction is low as expected. As low energy configurations accumulate, the trajectory converges to the DAO trajectory. The black cross mark stacked at the RMSE graph is these attempts to correctly locate the configuration near-optimal pathways. At iteration 186, GPAO converges and successfully locates the MEP. A comparison between the predicted energies on the Gaussian PES and EMT data sampled from the DAO trajectory is plotted on the right side of Figure 5. The RMSE between two calculations is 95 meV, or 5 meV/atom. These results show that GPAO predictions are in great agreement with the EMT results.

To validate the advantage of the GPAO method, the results of both DAO and GPAO are summarized in Table 3. The degree of freedom, a number of force calls, energy barriers, and the OM actions of final trajectories obtained with GPAO and DAO are tabulated. The degrees of freedom of ‘Au hopping’ and ‘Pt rearranges’ are 15 and 21, respectively. Because DAO requires Hessian calculations, the numbers of force calls for DAO calculations of ‘Au hopping’ and ‘Pt rearrange’ are 420,750 and 768,600 respectively, while those of GPAO are 15 and 115. GPAO predicts the energy barriers of three examples as 0.37 eV and 1.06 eV while DAO, which are the references for GPAO, is 0.37 eV and 1.06 eV respectively. These result show that GPAO and DAO calculations match well. The surface reaction examples demonstrate that GPAO reduces computational cost more than DAO by three to four orders of magnitude while giving suppressing the errors of energy barriers under 0.01 eV/atom.

Encouraged by the significant gain of computational efficiency by GPAO in five orders of magnitude, we investigate multiple conformational transition paths of alanine dipeptide through ab initio calculations. Alanine dipeptide has been widely used as a test system for computational methods in biophysics [36]. It has two rotatable dihedral angles, the $\phi$ and $\psi$ angles and two stable conformations, $C7_{\mathrm{eq}}$ and $C7_{\mathrm{ax}}$ (Figure 6). The PES’s shape shows that multiple conformational transition pathways from one conformation to the other are possible (Figure 6). We combine GPAO with the Action-CSA method that finds multiple transition pathways through the global optimization of action on a trajectory space and call it GP-CSA. The multiple transition pathways of alanine dipeptide are sampled using GP-CSA calculations with three different total transition times, $t_{\mathrm{total}}=Ndt=$ 3, 6, and 9 ps. Each pathway converged locally using GPAO. Further details on the converging process of alanine dipeptide are provided in Figure LABEL:subfig:ADGPAO.

The potential energies and forces of alanine dipeptide are calculated using the ab initio based density-functional-theory package VASP [51, 52]. The constrained MD method of VASP is used to optimize the conformations of $C7_{\mathrm{ax}}$ and $C7_{\mathrm{eq}}$ with fixed $\phi$ and $\psi$ angles. Since the structures of alanine dipeptide are almost exclusively defined by $\phi$ and $\psi$ angles, its potential can be expressed as follows: $V\left(\left\{\mathbf{x}^{\left(n\right)}\right\}\right)=V\left(\phi^{\left(n\right)},\psi^{\left(n\right)}\right)$. This reduction on coordinates $\mathbf{x}^{\left(n\right)}=\left[\phi^{\left(n\right)},\psi^{\left(n\right)}\right],\mathbf{v}^{\left(n\right)}=\left[\phi^{\left(n+1\right)}-\phi^{\left(n\right)},\psi^{\left(n+1\right)}-\psi^{\left(n\right)}\right]/{dt}$ is extended on the formulation of the modified OM action for alanine dipeptide. Our $\Theta_{\mathrm{OM}}$ for alanine dipeptide is

where $\nabla V\left(\left\{\mathbf{x}^{\left(n\right)}\right\}\right)=\left[\partial_{\phi}V\left(\left\{\mathbf{x}^{\left(n\right)}\right\}\right),\ \ \partial_{\psi}V\left(\left\{\mathbf{x}^{\left(n\right)}\right\}\right)\right]$ indicates the potential difference with respect to $\phi$ and $\psi$, $I_{\phi}^{\left(n\right)}\left(\phi,\psi\right)$ and $I_{\psi}^{\left(n\right)}\left(\phi,\psi\right)$ indicate the moments of inertia associated with $\phi$ and $\psi$ angles of $n$’th image, respectively. We use the following periodic kernel to describe the PES of alanine dipeptide with GP:

where $\mathbf{x}=(\phi,\psi)$ is a two-dimensional reaction coordinate ($D=2$), $d$ is a dimensional index, and $l,\sigma_{f}$ are isotropic hyperparameters.

Among the pathways obtained with GP-CSA, distinctive pathways with low action values are plotted in Figure 7. Path A in Figure 7 has the lowest $S_{\mathrm{OM}}$ value for all transition times, suggesting that Path A is the most dominant pathway between $C7_{\mathrm{eq}}$ and $C7_{\mathrm{ax}}$ conformations. This result is consistent with previous studies conducted with molecular mechanics potentials [36, 53]. The maximum potential energy of a pathway estimated with GP, $\mu^{\mathrm{(max)}}$ is near -128.48 eV (Table. 4). For all three different $t_{\mathrm{total}}$, GP-CSA consistently locates the correct saddle points.

The $S_{\mathrm{OM}}$ value of Path D with a $t_{\mathrm{total}}$ of 3 ps is 11.282 eV, which is high enough to make the path improbable. However, the action value decreases drastically to 5.551 eV and 4.887 eV when $t_{\mathrm{total}}$ elongates to 6 ps and 9 ps. As $t_{\mathrm{total}}$ increases, Path D becomes longer, allowing the pathway to detour the high potential region and to seek a lower valley of the PES (see Path D in Figure 7b and c). This tendency is also observed in the $\mu^{\mathrm{(max)}}$ values in Table 4. $\mu^{\mathrm{(max)}}$ of Path D at 3 ps is 0.2 eV higher than the other transition times.

It is noticeable that the relatively longer pathways, i.e, Path E, F, and G, are not observed in simulations with $t_{\mathrm{total}}=3$ ps. As $t_{\mathrm{total}}$ becomes longer, GPAO samples additional pathways, which require longer $t_{\mathrm{total}}$ to occur. $S_{\mathrm{OM}}$ values of Path E obtained with $t_{\mathrm{total}}=$ 6 ps and 9 ps are 8.547 eV and 6.745 eV, and $\mu^{\mathrm{(max)}}$ values are $-129.237$ eV and $-129.316$ eV, respectively. The maximum potential energy of Path E of $t_{\mathrm{total}}=6$ ps is 0.1 eV higher than that of Path E with $t_{\mathrm{total}}=9$ ps. This difference arises because Path E with $t_{\mathrm{total}}=6$ ps did not pass the true saddle point due to a short transition time. However, Path E with $t_{\mathrm{total}}=9$ ps correctly passed the saddle point (see E in the Figure 7c).

With $t_{\mathrm{total}}=6$ ps, two additional pathways, Path F and Path G are found, which are similar to the most dominant pathway, Path A. These two pathways are longer pathways with higher $S_{\mathrm{OM}}$ values, 4.56 eV and 5.293 eV, than Path A. These results clearly indicate that our GPAO method finds multiple transition pathways of the conformational change of alanine dipeptide accurately with a QM potential, which has not been reported previously.

GPAO also improves the Hessian’s accuracy by accumulating data points, which allows one to utilize potential functions whose Hessians are hard to calculate for path sampling. The conventional way of obtaining the Hessian of a DFT potential is to utilize a finite difference calculation. That is, to acquire a Hessian at an arbitrary point $\phi^{(n)},\psi^{(n)}$, one needs to compute $\nabla V(\phi^{(n)},\psi^{(n)}),\nabla V(\phi^{(n)}+\varepsilon,\psi^{(n)})$ and $\nabla V(\phi^{(n)},\psi^{(n)}+\varepsilon)$. These two additional calculations triple the computation cost to evaluate $\Theta_{\mathrm{OM}}$ preventing the use of finite difference calculation from being applied to relatively large systems. Thus, a fast evaluation of a potential using GP makes using Hessian for practical problems feasible.

In this study, we demonstrate that the GP algorithm dramatically enhances the efficiency of searching minimum action pathways by approximating a PES accurately with much-reduced numbers of energy evaluations with diverse examples: the MB potential, isomerization of small molecules, surface reactions, and conformational transition of alanine dipeptide. For most cases, the transition pathways obtained with GP and DAO show little discrepancies. The errors of energy barriers are also marginal. These results indicate that transition pathway search using the modified OM action with a total energy conservation restraint on an approximate PES is an efficient and accurate approach. In addition, compared to the direct optimization of the OM action itself, large fluctuations of kinetic energies are effectively removed without affecting the OM action of a pathway.

The gradients of OM action require the Hessian calculation of a potential function, which is computationally expensive to obtain in general. Thus, combining GP with action optimization gives not only superior computational efficiency also accurate approximations to the gradients of the OM action without actual Hessian calculations. This advantage is especially significant for the case where the evaluation of potential energies and their Hessians are challenging to obtain, i.e., the Hessian of DFT-based potentials is often incorrect. On the contrary, GPAO accumulates gradient information at multiple data points and leads to the accurate estimation of the Hessian from a Gaussian PES. Since the Hessian of a Gaussian PES can be estimated accurately, the OM action of a pathway can be readily obtained with DFT-based potentials with enhanced computational efficiency by five orders of magnitude. Additionally, we also showed that our GPAO method automatically finds suitable hyperparameters, such as target energy.

This study also demonstrates that combining GP with Action-CSA [36] enables the sampling of multiple pathways. Due to GPAO’s efficiency, multiple conformational transition pathways between the two stable conformations of alanine dipeptide are successfully obtained using a DFT-based potential. This is the first work that sampled multiple low action pathways of alanine dipeptide with a DFT potential to the best of our knowledge.

Tackling higher-dimensional problems in an efficient way is still a daunting challenge in machine learning fields and other related fields. This paper proves that our GPAO method properly works for relatively small-dimensional problems, whose degree of freedom is less than 30. However, the method does not guarantee its efficiency on problems with a large degree of freedom yet. Fortunately, many methods have been suggested to solve problems in a high dimensional space using various types of descriptors [54, 55, 56, 57]. Furthermore, finding reaction pathways of the systems with many degrees of freedom can be applied to various problems such as finding drug binding pathways, battery, and protein folding. Thus, we believe that GPAO will open up new possibilities for such problems.

The source code of GPAO can be found in the GitHub repository
(https://github.com/schinavro/taps). GPAO consists of multiple tools in the open-source package TAPS. TAPS incorporates tools needed for data-driven pathway search methods such as a database or atomic calculator. For atomic calculations, we use the atomic simulation environment (ASE)  [50] package, which can easily link with atomic calculators such as VASP [51, 52], Quantumespresso [59], OpenMX [60], and SIESTA [61].