Non-empirical weighted Langevin mechanics for the potential escape problem: parallel algorithm and application to the Argon clusters

TL;DR

This paper introduces a parallelized, non-empirical stochastic algorithm for accurately finding minimum-energy escape paths in potential surfaces, demonstrated on Lennard-Jones argon clusters, enabling detailed rare reaction simulations.

Contribution

It develops an MPI-parallelized implementation of a non-empirical Langevin-based algorithm for MEP search, with improved efficiency and automatic initialization, applied to argon clusters.

Findings

Successful tracking of reaction paths in argon clusters with 7-38 atoms

Reduced number of walkers needed for accurate MEP detection

Demonstrated potential for non-empirical simulation of rare reactions

Abstract

Recently a non-empirical stochastic walker algorithm has been developed to search for the minimum-energy escape paths (MEP) from the minima of the potential surface [J. Phys. Soc. Jpn. 87, 063801 (2018)]. This algorithm is novel in that it tracks the MEP monotonically and does not use the whole Hessian matrix but only gradient and Laplacian of the potential. In this work, we implement an MPI-parallelized version of this algorithm in a simple way. We also explore efficient ways to reduce the number of walkers required for the accurate tracking of the MEP and generate initial positions automatically. We apply the whole scheme to the Lennard-Jones argon cluster with 7-38 atoms to demonstrate the successful tracking of the reaction paths. This achievement paves the path to non-empirical simulation of rare reactions without coarse-graining or artificial potential.

Figures

Equations

$$d{x}_i=-\frac{{\partial }_{i}U(\mathbf{x})}{\Gamma }dt+\sqrt{\frac{2k_{\mathrm{B}}Tdt}{\Gamma }}{W}_i,$$

$${\partial }_{t}p(\mathbf{x},t)$$

$$\equiv$$

$${\partial }_{t}q(\mathbf{x},t)=\{{\hat{L}}_{\mathrm{Sm}}^{\prime }+{\hat{L}}_{\mathrm{rate}}\}q(\mathbf{x},t)$$

$${\hat{L}}_{\mathrm{Sm}}^{\prime }$$

$${\hat{L}}_{\mathrm{rate}}\equiv L_{\mathrm{rate}}(\mathbf{x},t)$$

$${\partial }_t\ln C(t)$$

$$F(\mathbf{x})$$

$$q(\mathbf{x},t+\tau )=\exp \{[{\hat{L}}_{\mathrm{Sm}}^{\prime }+{\hat{L}}_{\mathrm{rate}}]\tau \}q(\mathbf{x},t).$$

$$\exp \{[{\hat{L}}_{\mathrm{Sm}}^{\prime }+{\hat{L}}_{\mathrm{rate}}]\tau \}$$

$$\simeq \exp \{{\hat{L}}_{\mathrm{Sm}}^{\prime }\tau /2\}\exp \{{\hat{L}}_{\mathrm{rate}}\tau \}\exp \{{\hat{L}}_{\mathrm{Sm}}^{\prime }\tau /2\}.$$

$$\mathbf{x}(t+\tau )-\mathbf{x}(t)$$

$$=-\frac{1}{\Gamma }\nabla (U(\mathbf{x}(t))-2V(\mathbf{x}(t)))\tau +\sqrt{\frac{2k_{\mathrm{B}}T\tau }{\Gamma }}\mathbf{W}$$

$$⟨f(\mathbf{x}){⟩}_{q(t)}\equiv \int d\mathbf{x}f(\mathbf{x})q(\mathbf{x},t)\simeq \frac{1}{N_{\mathrm{w}}}\sum _{i=1}^{N_{\mathrm{w}}}f({\mathbf{x}}_i(t)).$$

$$V(\mathbf{x})=(1-\delta )U(\mathbf{x}).$$

$${\tilde{p}}_{iw}=\frac{p_{iw}}{{\sum }_{jw}{p}_{jw}};p_{iw}=\exp \{L_{\mathrm{rate}}({\mathbf{x}}^{(iw)},t)\tau \}.$$

Full text

††journal: Computer Physics Communications

1 Introduction

Efficient simulation of structural transitions and chemical reactions of matter in atomic scale–which concerns folding of proteins, molecular reactions, diffusion of impurities in solids, nucleation, etc.–is a common issue across a wide range of research fields - chemistry, physics, biology and their intersections. The difficulty of this task originates from the rarity of these phenomena in the simulation time scale; the time step $\Delta t$ for the simulation is required to be as small as picoseconds or even sub-femtoseconds, whereas the typical timescale of the phenomena is more than microseconds. This difference in the time scales renders their straightforward simulations almost impossible and alternative ways are usually sought. This paper concerns a development of the method to manage the simulation of a class of the rare events described by infrequent transitions between the minima of the potential surface (or free-energy surface) defined in the configuration space, induced by thermal fluctuation.

The frequency of the transition across the potential minima under the thermal fluctuation is dominated by the minimum potential energy barrier between them, and therefore the seeking of the paths connecting the potential minima with minimum energy barrier (minimum energy paths; MEPs) is of particular importance. Many existing methods can be applied when either the product (destination of the transition in the configuration space) or the reaction coordinate (collective coordinates which well characterize the transition) is given a priori, for example, the blue-moon [1, 2], targeted dynamics [3], steered dynamics [4] hyperdynamics [5, 6], metadynamics [7], adaptive biasing force [8] and many others. In the former case, various path optimization or sampling methods are available, whereas in the latter, one can add artificial force with respect to the reaction coordinates to the potential force  [9] or execute a weighted sampling dynamics referring to the coordinates [10, 11, 12, 13].

Recently, the authors has proposed an atomistic simulation method that realize the transition between the potential minima by tracking the path from a minimum to the saddle point through the valley of the potential surface [14]. This method does rely neither on the a priori knowledge of the product nor the reaction coordinates. The goal of this development is to establish a non-empirical simulation method that is applicable to the situation when the possible products are unknown. Although there are also methods aiming at the same goal such as the eigenvector following [15, 16], gentlest ascent [17] and dimer methods [18], the present method has advantages and disadvantages: It does not require us to calculate the Hessian matrix of the potential and its eigenpairs and yields monotonic tracking of the MEPs. Instead, we need to run a number of walkers interacting through a stochastic process. Further development of the method could therefore give us an interesting alternative for the abovementioned goal.

The method has been demonstrated only with a very simple two-dimensional potential surface. In this study, we seek for the possible application to larger systems. In the path search in the higher dimensional systems, the entire volume of the configuration space diverges, which is why the structure and path search is generally difficult. Nevertheless, the dimension of the MEPs is in principle always one and the number of the relevant paths is, in the low temperature, a few. With these contrasting facts we conceive that the path search that tracks the valley lines to the saddle points could be possible with computational cost and time that slowly increase with respect to the dimension of the system.

2 Non-empirical weighted Langevin mechanics for tracking minimum energy escape paths

Here in a brief we show concept, mathematical model and then main algorithmic steps of our method published in Ref. [14]. We address the molecular mechanics described by the following Langevin equation

[TABLE]

where $U({\bf x})$, $k_{B}$ and $T$ are the potential, Boltzmann constant and temperature, respectively. $i$ is the index for the degrees of freedom. ${\bf W}$ is the vector whose components are randomly generated from the standard normal distribution at each step. $\Gamma$ is the friction constant. This equation represents the motion of the degrees of freedom ${\bf x}=\{x_{1},x_{2},\dots\}$ under the potential $U({\bf x})$ and thermal fluctuation. When the temperature $T$ is low, the trajectories escaping from one minimum is observed with extremely low probability. In order to see how such trajectories can be rendered more frequent, we first move to the Smoluchowski equation describing the distribution function of ${\bf x}$, $p({\bf x},t)$, given by [19]

[TABLE]

This equation has a general intriguing property that, regardless of the initial distribution $p({\bm{x}},t=0)$, it converges to the Boltzmann distribution $p_{\rm B}({\bm{x}})$ in the infinite $t$ limit. This means that, starting from any localized initial distribution $p({\bm{x}},t=0)~{}\delta({\bm{x}}-{\bm{x}}_{0})$, there must be a component in $p({\bm{x}},t)$ that climbs up the potential slope. An intuitive way to enhance such component is to divide it by $p_{\rm B}({\bm{x}})$.

Along this line of consideration, the property of the modified distribution $q({\bm{x}},t)\propto p({\bm{x}},t)/[{\rm exp}(-(1-\delta)U({\bf x})/(k_{\rm B}T))]$ has been examined with one-dimensional harmonic potential $U(x)=\alpha x^{2}/2$ [14]. The analytical solution of $q({\bm{x}},t)\propto p({\bm{x}},t)$ with the initial condition $p({\bm{x}},0)=\delta({\bm{x}}-{\bm{x}}_{0})$ has been found to have a property as follows: (i) It has the Gaussian form and its width gradually spreads with $t$, (ii) if $\delta<1/2$, its center first increases by $\propto x_{0}t$ and, as expected from the stationary solution $p_{\rm B}({\bm{x}})$, converges to the center of the parabola in a long time. A scenario is conceived from these features: In many-dimensional space, If ${\bm{x}}_{0}$ is placed along the “valley line” of the potential, where in one direction the potential has nonzero gradient and in others zero, the center of $q({\bm{x}},t)$ goes farther along the line upward for a while. By tracking the center of $q({\bm{x}},t)$, we can therefore draw the MEPs as the valley lines extending from a potential minimum. This simplistic argument has been actually validated in the two-dimensional case [14] with the algorithm appended later.

We then construct a microscopic stochastic algorithm to reproduce $q({\bm{x}},t)$. From the Smoluchowski equation [Eq. (3)], the corresponding Master equation for $q({\bm{x}},t)$ with general transformation $p({\bm{x}},t)=C(t){\rm exp}[-V({\bm{x}})/k_{\rm B}T]q({\bm{x}},t)$ is given by [20, 14]

[TABLE]

with

[TABLE]

Coefficient $C(t)$ is introduced so that the norm of $q$ is kept to unity. We define the average of the function $f({\bm{x}})$ by $\langle f\rangle_{q(t)}=\int d{\bm{x}}q({\bm{x}},t)f({\bm{x}})$. The finite time evolution of $q$ is formally represented as

[TABLE]

We then construct the stochastic process that operates on the walkers having individual values of ${\bm{x}}$ and reproduces $q({\bm{x}},t)$ as their aseembly. The Suzuki-Trotter decomposition [21, 22, 23] of Eq. (4) yields:

[TABLE]

In terms of the walker dynamics, operator ${\rm exp}\{\hat{L}^{\prime}_{\rm Sm}\tau/2\}$ is the time evolution

[TABLE]

applied to all the walkers. Operator ${\rm exp}\{\hat{L}_{\rm rate}\tau\}$, being simply a multiplication of a function of ${\bm{x}}$ on $q({\bm{x}},t)$ to $q({\bm{x}},t)$, is implemented as the importance sampling [24] as described below. The statistical average of physical quantity $f({\bm{x}})$ with respect to $q({\bm{x}},t)$, which is used in Eqs. (6) and (7), is then recast to the average with respect to the walkers

[TABLE]

We here append the specific procedure for the time evolution Eq. (10) for each time step:

${\rm exp}\{\hat{L}^{\prime}_{\rm Sm}\tau/2\}$:

Evolve the walkers independently by the Langevin equation Eq. (11) 2. 2.

${\rm exp}\{\hat{L}_{\rm rate}\tau\}$:

Do the following for all the walkers $i=1$–$N_{\rm w}$ independently. If $r\equiv{\rm exp}\{L_{\rm rate}({\bm{x}}_{i}(t),t)\tau\}<1$, remove that walker by probability $1-r$; otherwise, add $\lfloor(r-1)\rfloor$ walker(s) with the same value of ${\bm{x}}$ as $i$ and add one more by probability $r-1-\lfloor(r-1)\rfloor$. Modify $N_{\rm w}$ accordingly. 3. 3.

${\rm exp}\{\hat{L}^{\prime}_{\rm Sm}\tau/2\}$:

Repeat step (1)

According to the argument above, we have proposed as a specific form of $V({\bm{x}})$ as

[TABLE]

With this form, the original potential function $U({\bm{x}})$ is straightforwardly devised and its convex/concave structure is reflected. Note that the Laplacian of $U({\bm{x}})$ is then needed to evaluate $L_{\rm rate}({\bm{x}})$ for each walker. We have demonstrated that this setting enables us to track the valley lines of $U({\bm{x}})$ from a known potential minimum to the saddle points for a case of two-dimensional potential surface. In this work, we develop systematic algorithms to make the present scheme applicable to higher dimensions.

3 Scheme of the algorithms

The scheme to search for various possible MEPs in general systems is composed of two steps: First, generate initial points ${\bm{x}}_{0}$ that are connected to some saddle points through the valley lines of the potential surface. Second, evolve the walkers from the initial distribution $q({\bm{x}},0)=\delta({\bm{x}}-{\bm{x}}_{0})$ according to the abovementioned procedure. In the previous two dimensional case [14], the first step was apparently easy, but in general dimensions, where the many-dimensional potential landscape cannot be straightforwardly visualized, we need a special treatment to specify appropriate ${\bm{x}}_{0}$. We therefore develop a method for this task, according to a principle that is obviously applicable to many dimensions. Afterwards, a simple and efficient parallel implementation of the second time evolution step is introduced. The success rate of tracking of the MEPs, defined with the frequency to generate the trajectories to the saddle point without straying from the proper paths, has been shown to depend strongly on the parameter $\delta$ in the biasing potential as $V({\bm{x}})=(1-\delta)U({\bm{x}})$, as well as the number of walkers $N_{\rm w}$. We also explore efficient control of $\delta$ to gain high success rate with relatively small $N_{\rm w}$.

3.1 Generation of initial atomic positions

In many dimensional cases the “valley lines” of the potential can be formally defined as the lines generated by tracking upward the directions of local minimal potential gradient. That is why the walkers on a valley line are expected to relax more slowly to the minimum than those off the lines. This principle is demonstrated by a molecular mechanics simulation [Eq. (1)] for the 2D case with random initial positions and low temperature (see Supplementary movie). During the relaxation, the walkers drift to the stable energy minima in two stages. At first, the walkers run to the valley lines and then drift along the lines to the minimum slowly. Of course, the walker velocities depends on potential surface, but the relative velocities are slower on the valley lines. We assume that this simple principle applies to many-dimensional cases: The slowest walkers in the relaxation from random distribution must be on the reaction path.

A simple algorithm to generate an initial position ${\bm{x}}_{0}$ according to the above principle is to execute usual molecular mechanics simulation with relatively high temperature $T_{\rm high}$ for a while, set a threshold energy $E_{\rm th}$, relax the walkers with the heat bath quenched and specify the walker which comes down to $E_{\rm th}$ last. With this algorithm, however, it is obvious that only one valley line related to the smallest descent valley line can be tracked. To gain the diversity of the MEPs generated in the later time evolution step, we need to consider the different speeds of relaxations along various paths. For this purpose, we define atomic energy $E_{a}(a=1,2,\dots,N_{\rm atom})$ with the interatomic potential $U_{ab}$ by $E_{a}=\sum_{b}U_{ab}$ and characterize the valley lines by the atoms showing largest energy deviations $\Delta E_{a}\equiv E_{a}-E_{a}^{\rm min}$ ($E_{a}^{\rm min}$: value of $E_{a}$ with the optimized atomic configuration). The background consideration of this strategy is as follows: Any reaction starts from the rupture of interatomic bonds. The potential energies between neighboring atoms are increased during the breaking of the bond. The reaction paths are therefore characterized by the combination of atoms showing largest $\Delta E_{a}$.

We have thus reached to the general algorithm to generate the initial positions that are related to diverse reaction paths, as summarized below (Algorithm 1). The algorithm has two stages: the first is to heat the system and the second one is to cool it. Simply speaking, we group the walkers according to the atoms showing $n$th largest $\Delta E_{a}$ and find the most slowly relaxing walkers among those groups.

3.2 Parallel implementation of the time evolution with strict conservation of $N_{\rm w}$

As we explained above, the time evolution operator Eq. (10) consists of the Langevin and the replication/removal steps. A simplistic implementation of the latter suffers from the fluctuating number of walkers $N_{\rm w}$, which is not suitable for efficient parallelization of the code. We therefore adopt an exactly number-conserving copying/removal algorithm in Ref. [25], where this algorithm is named ”[eq]”. Let us here append the algorithm

For the respective walkers (indicated by index $iw$) calculate

[TABLE]

$p_{N_{w}+1}$ is defined as $1$ for convenience. 2. 2.

Generate a random number $d\in[0,1/N_{w})$. 3. 3.

For each walker $iw$, increment $jw$ to find the range so that $\alpha_{iw}\equiv d+(iw-1)/N\in[\tilde{p}_{jw},\tilde{p}_{jw+1})$; the $iw$th walker then inherits the position ${\bf x}^{(jw)}$.

Through this algorithm, in the $N_{w}\rightarrow\infty$ limit the whole distribution function $q({\bf x},t)$ is modified to $q({\bf x},t+\tau)\propto{\rm exp}\{L_{\rm rate}({\bf x}^{(iw)},t)\tau\}q({\bf x},t)$ with its norm $\int d{\bf x}q({\bf x},t)$ kept constant.

We have implemented a flat-MPI code for this algorithm. The pseudocode is as follows:

3.3 Parameter settings for successful tracking

In the previous paper it has been found that, during the time evolution, the center of $q({\bm{x}},t)$ has fluctuations and the path search sometimes failed: the walkers stray normal to the MEP and move up the potential slope to a local maximum. This failure is understandable as the case where a large fraction of walkers are accidentally driven normal to the valley due to the random effect of the thermal fluctuation and copying/removal step. To reduce the probability of such failure with feasible number of walkers $N_{\rm w}$, we here consider appropriate setting of the parameters governing the simulation; the time step divided by the friction constant $dt/\Gamma$, temperature $T$ involved in the thermal fluctuation as well as the rate term, and $\delta$ entering the biasing potential as $V({\bm{x}})=(1-\delta)U({\bm{x}})$.

Temperature $T$ determines the spread of $q({\bm{x}},t)$ and its small spread allows us to describe the whole distribution accurately with a small number of walkers. On the other hand, too small $T$ induces radical ${\bm{x}}$ dependence of the rate term, as it enters $L_{\rm rate}$ in Eq.(10) as $\beta$. This effect can be mitigated by setting $dt/\Gamma$ small enough. A simple principle should therefore applies to $T$ and $dt/\Gamma$: Smaller is better (though they must be compromised not to make the evolution too slow).

The effect of $\delta$ is more complicated [14]. The sign of the potential force in the modified Langevin step [Eq. (11)] is inverted from the original $U({\bm{x}})$ with $\delta<1/2$, which is the condition needed for the center of $q({\bm{x}},t)$ to climb up the slope in a short time. Also, the speed of increase of the center and spread are both proportional to $1/\sqrt{\delta}$, which implies that, with too small $\delta$, the fluctuation of the walker distribution from the ideal $q({\bm{x}},t)$ could be amplified with time. According to our preliminary analysis, the probability of failure is actually increased with a time evolution and tends to be larger with smaller $\delta$.

In order to avoid the fluctuation problem, we set $\delta$ as near to $1/2$ as possible and introduce a “reset and pullback” step by some frequency; namely, ${\bm{x}}$ is set to their sample average for all the walkers to reset the sample distribution to the delta function, and later we execute a short simulation without $V({\bm{x}})$ to make sure that the distribution is well on the valley line.

The whole procedure to track the MEPs from a given initial position, with the reset and pullback steps, is summarized as Algorithm 3. The performance of the current treatment is examined later for the Lennard-Jones cluster with seven particle (Sec. 4.2).

4 Benchmark test for argon cluster

Let us show the performance of the present method with the Lennard-Jones (LJ) clusters with 7-38 atoms. The LJ clusters are well-studied systems and a good target to check our method. We used the model parameters in Ref. [26], which corresponds to the interatomic potential of argon atoms. In the present study, we did not introduce the cutoff length and therefore numerical cost to calculate all the component of the force vector and Laplacian of the potential for each walker is proportional to $N_{\rm atom}^{2}$.

The numbers of the degrees of freedom are then 114 at maximum, and the assembly of the walker represents the distribution function $q({\bm{x}},t)$ in the 114-dimensional space. We visualize the structures of the clusters with the sample average of the position $\langle{\bm{x}}\rangle_{q({\bm{x}},t)}$.

In application of our scheme, we assume the situation where any local minimum structure (reactant) is known but the neighboring minima (products) are unknown. Here we start the simulations with the known minima as summarized in Fig. 1, pentagonal bipyramid for LJ7, icosahedron for LJ13, and truncated octahedron for LJ38, and seek for the possible MEPs connected from those.

The simulation temperature in the MEP tracking process was $T=0.0001$ for LJ7, $0.001$ for LJ${}_{1}3$ and LJ${}_{3}8$, respectively. The time step $dt/\Gamma=0.0004$ for all the simulations. The parameter $\delta$ was adaptively controlled as mentioned in Algorithm 3. The total number of walkers were $N_{w}=3200$ for the simulations in Secs. 4.3 and 4.4

4.1 Parallel efficiency

Here we summarize the scaling performance of our MPI-parallelized code. The benchmark calculation has been performed for the Lennard-Jones clusters with 7 and 38 particles (LJ7 and LJ38). The total number of walkers were set from 3456 to 345600, though tracking of the MEPs can be successful with far fewer walkers as discussed below. The calculations were performed with K computer, which contains 82944 compute nodes in total and each nodes has two 4-core SPARC64 VIIIfx processors.

We measured the computational time for 100 timesteps and show the resulting speedup for various numbers of MPI processes in Fig. 2. There is a small excess of the speedup rate over ideal line in the region $N_{processors}<100$ for $N_{w}=345600$. It could be due to uniform memory distribution in the hard disk space for 1 processor with a huge number of walkers.

In the case of LJ7 the speedup saturated with the number of walkers per process around 300. This relatively quick saturation is probably because of the moderate computational cost to calculate the potential and its derivatives in comparison with the MPI data exchange between processors. In fact the saturation of speedup rate for LJ38 occurs when the number of walkers per process equals 36 and the computational cost to calculate the potential is increased by $(37\cdot 38)/(6\cdot 7)\approx 33.5$. Still, practically this performance was enough in this work for completing the single path tracking utilizing only 16 processors with 32000 time steps in 5 minutes for LJ7 and in no more than several ten minutes for LJ13. For LJ38 we used 80 nodes (640 processes) with 80000 time steps in 7 hours.

4.2 Success rate

To estimate the influence of $\delta$ controlling the strength of the biasing potential $V({\bm{x}})$ on a successful rate of simulation we repeated the simulation starting from a single initial position for LJ7 and LJ13 (Fig. 3) to count the number of successful tracking to the saddle point connected through the valley line. The results are compared with and without the reset and pullback steps and with various ways to control $\delta$. The number of walkers was 16, 160, 240, 320, 384, 480, 560, 800, 1600, 3200 and the number of repeating simulations was 20 for each point. For $N_{w}>3200$ all simulations were successful. It is clear that the most significant parameter in seeking the reaction paths is the delta parameter $\delta$. In the current scheme (Algorithm 3), the $\delta$ adaptive selection is performed by gradually increasing the influence of the biasing potential $V({\bf x})=(1-\delta)\cdot U({\bf x})$

On Fig. 3 there is no successful simulation with $\delta=0.3$ because under strong V(x) all simulations go to the higher order saddle points and the system is destroyed into several parts. If the delta parameter is too near to $1/2$ like 0.4, then the V(x) is too weak and the walkers do not go far up the potential slope. It is interesting that for the small number of walkers $N_{w}<1000$ the system is able to go out of balance due to fluctuations, an amplitude of which becomes higher with the small number of walkers. As a result, the system sometimes drifts up to the saddle point successfully. Of course, when the walker number is too small $N_{w}<400$, the fluctuations become too large and the system is destroyed again.

We also tried to calculate the reaction paths without the additional procedures; reset, pullback and adaptive $\delta$ (top left Fig. 3). For LJ7, $\delta=0.36$ was the minimal value for drifting up on the potential valley, that’s why it has a good rate of successful simulations. The procedures of the reset and pullback are also important only for region $N_{w}<1000$. We need to keep in mind that for successful simulation it is necessary to choose other parameters like time step $dt$, temperature $T$, simulation time $t_{max}$ and etc. From this point of view, short simulations in order to determine the parameter settings are generally essential. Remarkably, the number of walkers required for successful tracking of MEPs does not change so significantly with increasing the system dimension (bottom Fig. 3), For LJ7, LJ13 and LJ38, $N_{w}=3200$–$6400$ was enough.

4.3 Simulation results for argon cluster with 7 atoms

Utilizing the algorithm of seeking the initial atomic positions (Algorithm 1) for LJ7, we have determined 42 initial positions ${\bm{x}}_{0}$ as entrances to the MEPs from one stable state - pentagonal bipyramid. It should be noted only the simulations from 2 initial positions did not reach the saddle point (SP) in 32000 time steps, but the other 40 successfully yielded the paths to the SP. The LJ7 has $C_{5h}$ group symmetry as shown in Fig. 1, that’s why some simulations have the same trajectories. As a result, there are only 21 unique pathways (Fig. 4). We define the unique paths as a way through the unique SP and unique SP by the labels of atoms showing 1st to 3rd largest energy deviations $\Delta E_{a}$ with the corresponding ${\bm{x}}_{0}$ (atomic energy spectrum of SP).

The unique 21 SPs are shown in (Fig. 4). The curves of the potential energy along the paths to other metastable points are shown in Fig. 6, which goes through the a, c and l saddle points whose energy barriers agree with those determined in Ref. [26]. The corresponding reaction paths are visualized in Figs. 6, 8, 8. In the Supplementary movies, we have shown the dynamics of reaction paths during the simulations.

There are a lot of studies of reactions in LJ clusters, for example [27, 28, 29, 30, 31]. Almost all papers are concentrated on the topic of free-energy landscape and their energy minima, using the different methods: molecular dynamics [27, 28], parallel tempering technique with molecular dynamics [29], path-sampling scheme [30] and diffusion Mote-Carlo method [31]. The free-energy landscape is defined in a coarse-grained low-dimensional space with collective order parameter axes like $Q_{4}$ or $Q_{6}$. In our method, on the other hand, we directly generated the MEPs keeping all the atomic degrees of freedom. We also note that, differently from the existing methods using non-empirical potential force like metadynamics [7], we use the original potential $U({\bm{x}})$ to keep the landscape unmodified.

An interesting advantage of our method is possibility to reach high energy barrier as easy as low energy one. In comparison with the previous study [26], where authors found the 4 saddle points from the LJ7 pentagonal bipyramid, our method have found the 21 SPs (Fig. 4). The 21 SPs include the first order SP (a,b,c,l,n,r points in Fig. 4) and high order SPs.

4.4 Simulation results for argon LJ13 and LJ38 clusters

To check how the method behaves with increasing of dimensionality, we performed several simulations for LJ13 and LJ38 clusters. The parallel efficiency is increased with a number of atoms (Fig. 2), but the total number of walkers is still constant around $3200$–$6400$ (bottom of Fig. 3). It seems possible to apply the method to huge systems with several hundred or even several thousands of atoms without increasing the number of walkers so much, and that is our future task. In our opinion the reason for this “good” behavior of method is the fact that only several atoms take part in the reaction while a movement of other atoms is deterministic, with which the volume spanned by $q({\bm{x}},t)$ is kept feasibly small.

Fig. 9 shows the energy evolution of LJ13 cluster for two reaction pathways from icosahedron: the first is to the twist SP, beyond which the same icosahedron is reached (Fig. 11), and the second is to the capped nido icosahedron (Fig. 11). Notably, the previous study [26] didn’t mention the twist SP, despite it is clear to see this structure as unstable and SP for icosahedron - icosahedron transformation with twisted movement (see Supplementary movie). The structure transformation is shown in Fig. 13. The truncated octahedral configuration has very high symmetry and therefore corresponds to a deep energy global minimum, that is why we performed long simulation to reach to the SP: at first the 32000 time steps with a low temperature and, then we increased the $T$ in order to accelerate the simulation.

For small clusters (like LJ7) all reactions occur in accordance with diamond-square-diamond (DSD) process [26], where is usually one bond is destroyed and then another bond forms (Figs. 6,8,8). For LJ13 in the case of twisted reaction, there are 5 DSD processes at the same time (Fig. 11). Figure 11 shows the ”icosahedron - capped nido icosahedron” transformation, which has a vacancy formation on the surface and DSD process induced by a capped atom.

For LJ38 a reaction path is also complicated (see Supplementary movie): a pair of atoms simultaneously and synchronously change their coordinates during the reaction; in this case, two DSD processes occur, and also the second atom breaks two bonds, forming a vacancy (Fig. 13); the first atom forms a vertex, and the second is the base of a small pyramid on the surface of the initial configuration. In addition, during the reaction path, we observed a collective displacement of all atoms towards the atoms forming the pyramid, which reached its maximum at the saddle point and disappeared in a stable state. The amplitude of atomic displacement depends on the distance from atom to reaction couple of atoms. Using collective variable such detailed process is impossible to detect.

5 Conclusion

In this work we showed the concept, mathematical model and algorithmic steps of our method. The method is based on the Master equation (modified Smoluchowski equation), which is related to molecular mechanics described by the Langevin equation with weighted sampling. The model allows us to find the probability distribution through the biasing distribution, which is calculated under biasing potential. To solve the Master equation we use the stochastic walker algorithm with parallelized number-conserving copying/deletion. We verified the scaling performance of our MPI-parallelized code with a benchmark of LJ7 and LJ13 clusters. Using parallel calculation on K-computer we have shown that the speedup effectiveness increases with the dimension of the system. The computational cost is increased due to the calculation of the potential function and diagonal elements of Hessian matrix like $N_{atoms}^{2}\cdot N_{w}$, but $N_{w}$ is constant for 7-38 atoms.

We have designed the algorithm to seek the initial positions as entrances to the MEPs and to track the MEPs efficiently (adaptive delta parameter, reset and pullback) procedures in order to keep the walkers on the MEPs. The successful performance of tracking the MEPs in the LJ cluster systems with various number of particles with the same number of walkers suggests that larger systems could be feasible with the present simulation scheme, paving a path to the simulation of reaction mechanics free from artificial potential or a priori collective coordinates.

Acknowledgment This research was supported by MEXT as Exploratory Challenge on Post-K computer (Frontiers of Basic Science: Challenging the Limits). This research used computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID:hp160257, hp170244, hp180184).

6 Refences