Accelerated molecular-dynamics (MD) simulations based on hyperdynamics (HD) can significantly improve the efficiency of MD simulations of condensed-phase systems that evolve via rare events. However, such simulations are not generally easy to apply since appropriate boosts are usually unknown. In this work, we developed a method called OptiBoost to adjust the value of the boost in HD simulations based on the bond-boost method. We demonstrated the OptiBoost method in simulations on a cosine potential and applied it in three different systems involving Ag diffusion on Ag(100) in vacuum and in ethylene glycol solvent. In all cases, OptiBoost was able to predict safe and effective values of the boost, indicating that the OptiBoost protocol is an effective way to advance the applicability of HD simulations.
The temporal evolution of many materials systems is governed by rare events, where the system spends a relatively long time in one free-energy minimum before escaping and moving on to another one. While molecular dynamics (MD) simulations can, in principle, elucidate the atomic-scale processes and rates of rare events, the practical time scale is limited to the microsecond range. Accelerated MD methods were developed to address this issue.1–8 In this paper, we focus on hyperdynamics (HD)1,2—one accelerated MD technique that has been used in many studies.6,7,9–16 As we will elaborate below, a bias potential is added to the potential energy near the minima in HD. A careful design of the bias potential induces rapid transitions without affecting the relative transition frequency and a weighted time increment allows for long-time simulations that can exceed the microscale by many orders of magnitude.
One HD bias potential is given by the bond-boost method.17,18 In this method, the strain of the bond between atoms i and j is defined as
where Rij is the current distance between atoms i and j, and the distance when the potential energy is at a minimum. A bias potential, or “boost,” is applied to the maximally stretched bond—the bond that is most likely to “break” and result in a transition. For a fixed atom configuration R, the boost energy ΔV of the maximally stretched bond, with ɛmax is given by
where ΔVmax is a parameter that controls the magnitude of the boost and q is a cutoff parameter of the limiting strain at which the bias potential goes to zero. The q parameter should be set such that the boost is zero at transition states. MD trajectories are then run on a “boosted” potential-energy surface (PES) , where the potential energy is given by
There are two relevant times in HD simulations. The first is the MD time tMD, which is given by
where N is the number of MD time steps and Δt is the MD time increment, typically on the order of fs. The physical time is given by
It is evident from Eq. (2) that the values of both q and ΔVmax control the time-boost factor in Eq. (6). In this work, we adopt the value of q = 0.3, which has been shown to be a safe value for many cases in previous studies,16,17 and we focus our efforts on establishing an appropriate value of ΔVmax.
Figure 1 illustrates both the promise and the difficulty in using HD simulations with the bond-boost method. Figure 1 shows the escape rate from HD simulations as a function of ΔVmax run for motion on a one-dimensional cosine potential of the form
The boosted potential in these simulations is given by Eq. (3), and the bias potential is given by Eq. (2) with q = 0.3. The inset to Fig. 1 shows the time-boost factor from Eq. (6) as a function of ΔVmax. For each value of ΔVmax, we ran 500 different trajectories on the cosine potential using Mathematica® to integrate the equation of motion with the Verlet algorithm and the Andersen thermostat. The unit of mass is m0, the unit of time is t0, the unit of length is x0, and the simulations were run at a temperature of . We calculated the rate as the reciprocal of the average physical time, given by Eq. (5), to exit a minimum after entering. The exact rate to move from one minimum to another on the cosine potential is given by
where m is the mass, which we take to be 1.0m0.
In Fig. 1, we see that the simulated rate matches the exact rate given by Eq. (8) until , then the simulated rate begins to deviate from the exact value and the difference between the exact and simulated rate increases with increasing ΔVmax. In the inset to Fig. 1, we see the boost increases with increasing ΔVmax up to . At , the time-boost factor has reached a value of B = 81, meaning that the HD simulation is 81 times faster than regular MD. While the large efficiency afforded by HD simulations is promising, this can come at the expense of accuracy if ΔVmax is too large.
The origins of the deviation of the rate from the exact value in Fig. 1 can be seen from a plot of , given by Eq. (3) and shown in Fig. 2 for the cosine potential with different values of ΔVmax. For , the biased PES retains the same shape as the original potential (ΔVmax = 0). From the inset to Fig. 1, the boost for ΔVmax = 0.2 is around 2.3, which is non-negligible but modest in terms of what can be achieved. When ΔVmax increases to 0.8, the shape of begins to deviate from that of the original potential. However, the shape deviation is not large and the rate on the boosted surface in Fig. 1 is still the same as that on the original surface. From the inset to Fig. 1, the boost for ΔVmax = 0.8 is almost 15, which represents a significant increase in efficiency over ΔVmax = 0.2. When we reach ΔVmax = 2.0, the shape of V* is significantly distorted from the original potential, as there is a maximum where a minimum occurs on the original potential, and sub-minima appear near the transition states. Trajectories become trapped in the sub-minima on V* and they cannot easily access the region of the original minimum. As a result, the rate increases on V* compared to that on the original potential and becomes inaccurate.
Thus, care needs to be taken in choosing the boost parameter for HD simulations based on the bond-boost method. As Figs. 1 and 2 demonstrate, the magnitude of ΔVmax should not far exceed the value of the energy barriers in the system. In the current state of HD simulations, a safe and effective magnitude of ΔVmax is usually unknown and needs to be decided before running the simulation. In this work, we develop a method called OptiBoost to adjust the magnitude of the boost in HD simulations with the bond-boost method.
II. METHODS AND RESULTS
A. Cosine potential
The main ideas behind the OptiBoost method are as follows: (1) The bond-boost potential in Eq. (2) affects the potential energy of a bond on an arbitrary 3N-dimensional PES in a similar way that it affects the potential energy of the cosine potential shown in Fig. 2. (2) Dynamics on , but on the MD time scale, given by Eq. (4), displays a distinct signature of the boost potential. Regarding (1), the bond-boost potential is applied to bonds, whose potential energies as a function of strain are typically harmonic or slightly anharmonic around their minima and have a similar shape to the cosine potential. Regarding (2), the probability P that a trajectory originating in the potential well will escape to a neighboring well within a time tMD is given by
Figure 3 shows a plot of Eq. (9) (dashed lines) as a function of ΔVmax for three different time intervals. From Fig. 3, we see has an approximate sigmoid form in each case. This form can be seen to originate from Eq. (9), which we can write as
If all the rates are small, as they are in this example (and as is generally the case for rare events), and the time intervals tMD are sufficiently short, we have
where we note that . Equation (12) has a sigmoid shape. As we can anticipate from Fig. 2, rMD is a complicated function of ΔVmax, though we expect rMD to increase with increasing ΔVmax and approach a constant value when ΔVmax is sufficiently large. Both features will promote a sigmoid curve for with increasing ΔVmax, though not of the exact form as Eq. (12). Because we will not generally know the exact relationship between rMD and ΔVmax for any arbitrary system, we can write P as
This equation also has a sigmoid shape in which we can handle the fact that using three adjustable parameters (A, B, and C).
To test Eq. (13), we ran a series of HD simulations on the boosted cosine potential for various values of ΔVmax. We calculated the value of P as an average over 1000 trajectories, each with a duration of tMD. Figure 3 shows the results from the simulations as points and a fit of the simulation results to the sigmoid curve in Eq. (13) as lines. It is evident that there is excellent agreement between the exact values of P (dashed lines), the simulation values (points), and the fits to a sigmoid curve (solid lines).
The excellent fit of the sigmoid curve to the P data in Fig. 3 is an opportunity to define an appropriate boost parameter. Namely, reaches an asymptotic value for large ΔVmax because developed a substantial maximum where the original potential had a minimum, resulting in sub-minima that trap the trajectory near the transition state. We can identify the value of ΔVmax where P begins to reach an asymptotic value from the minimum in . This value of ΔVmax should represent an aggressive boost. From this perspective, the inflection point in , where and exhibits a maximum, is a well-defined quantity that should provide a safe estimate of ΔVmax.
Figure 4(a) shows , and Fig. 4(b) shows for the simulations shown in Fig. 3. The maxima in Fig. 4(a) and the minima in Fig. 4(b) are listed in Table I. It is evident from Figs. 1 and 2 that the barrier of the cosine potential is and that HD simulations become inaccurate when ΔVmax becomes much greater than the barrier. From Table I, we see that ΔVmax implied from the maximum in is less than the barrier and the point of minimum is somewhat larger than the barrier. Though the values of ΔVmax vary for different time intervals, they are similar. From the viewpoint of computational efficiency, it is advisable to use the smallest time interval tMD that will yield the sigmoid form—producing P ≅ 0 for ΔVmax around zero and a non-zero plateau for large ΔVmax. From Fig. 3, we see tMD as small as 0.2t0 are acceptable for the cosine potential. Generally, the smallest acceptable time interval depends on the barrier for the transition under consideration.
|t ( .||Maximum in ( .||Minimum in .|
|t ( .||Maximum in ( .||Minimum in .|
From Table I, it is evident that maxima on provide safe and conservative estimates for ΔVmax. In Fig. 1, these values all yield the exact rate, with time-boost factors of B ≈ 15. ΔVmax estimates from minima on also yield exact rates in Fig. 1, with time-boost factors of B ≈ 30—however, the HD rate begins to deviate from the exact value immediately beyond these estimates. Thus, ΔVmax estimates from are safe and those from are bold, but likely still satisfactory. In summary, the OptiBoost method involves (1) obtaining escape probabilities P as a function of ΔVmax from a series of short trajectories; (2) fitting to the sigmoid function in Eq. (13); and (3) delineating the boost range from the maximum in and from the minimum in .
B. Ag diffusion on Ag(100)
To demonstrate that the results for the cosine potential can be observed for other systems, we applied OptiBoost to three different systems involving Ag atom diffusion on Ag(100): (1) an Ag atom in vacuum; (2) an Ag dimer in vacuum; and (3) an Ag atom in the presence of ethylene glycol (EG) solvent. The third case is relevant to understanding the growth of Ag nanocrystals in solution.19–21 We simulated these system using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code,22 version 29 Oct 2020, compiled with the “REPLICA,” “MANYBODY,” “MOLECULE,” “KSPACE,” “USER-MISC,” and “PYTHON” packages. LAMMPS was used as a package to Python, for which we also imported “multiprocessing,” “os,” “random,” “shutil,” and “numpy.” We note that the parallel implementation of HD in LAMMPS has been described by Plimpton recently.16
We used an embedded-atom method (EAM) potential for Ag23 and the metal-organic many-body (MOMB) force field for the Ag–EG and EG–EG interactions.24–27 The size of the Ag substrate in the vacuum systems was 8 × 8 × 6 atoms, and the substrate consisted of 6 × 6 × 6 atoms with 48 EG molecules in the Ag–EG system. Prior to production runs, we performed equilibration in the NPT ensemble using the Nosé–Hoover thermostat and barostat for 10 ns to account for thermal expansion. All systems were equilibrated at 433 K and 1 bar. The molar mass of EG is 62.02 g/mol, and the volume of EG in the system was 17.5 × 17.5 × 13 Å3 after NPT equilibration. Using 48 EG molecules, the density is 1.24 g/ml, a value somewhat higher than the experimental value of 1.1 g/ml.28 Subsequent HD production runs were performed in the NVT ensemble using the Langevin thermostat. The MD time step was 1 fs in all runs.
For the HD simulations, we use q = 0.3 [see Eq. (2)]. There are several specific LAMMPS parameters for HD simulations.16 One is the bond cutoff Dbond, which we take to be 3.32 Å. The bond cutoff is the distance over which LAMMPS defines bonds for the bond-boost method. As the simulation runs, a check is performed every 1000 time steps to determine whether an event has occurred. In the solvent system, a check is also performed every 500 time steps and the observed event rate did not change. The check consists of quenching the system using the “quickmin” method, with dimensionless energy and force tolerances of 0.1, a maximum of 40 iterations, and 50 force evaluations. An event is said to have occurred if the displacement of the new quenched state from the current quenched state is greater than a distance of Devent, which we take to be Devent = 1.2 Å. We note that this value of Devent is ∼70% of the nearest-neighbor distance.
In LAMMPS, our OptiBoost HD simulations are driven by Python, in a protocol where N short trajectories are launched in parallel for each value of ΔVmax. The outcome of each trajectory is either “yes” if an event occurred or “no” if no event occurred. We then determine the event probability as the number of “yes” trajectories (Nyes) divided by the total number of launched trajectories, i.e., . All analyses of (i.e., curve-fitting to the sigmoid function and determination of appropriate values of ΔVmax) are done in Python.
1. Ag diffusion on Ag(100) in vacuum
We first performed HD simulations of a single Ag atom diffusing on Ag(100). Figure 5 shows the physical hopping rate as a function of ΔVmax. Each data point in Fig. 5 is an average over five simulations, with each run covering 40–200 ns of MD time [tMD in Eq. (4)]. The hop rate fluctuates around a value of 0.035 hops/ns for ∆Vmax between 0.0 and 0.4 eV and then decreases for larger ∆Vmax. This trend is opposite to that for the cosine potential in Fig. 1, and the discrepancy can be attributed to the differences in way LAMMPS detects events and the way we checked for events in our Mathematica simulations.
In the Mathematica code for the cosine potential, the locations of the transition states are obvious, and we checked every time step for a transition-state crossing. In LAMMPS, an event check is performed every 1000 time steps and an event is detected when the current quenched state of the atom displaces longer than a distance of Devent from the previous one. The quench process in HD uses the “quickmin” minimization routine without a strict tolerance, which implies the quenched state is not necessarily an energy minimum. When ΔVmax is large, the trajectory is confined near the transition state (as we see in Fig. 2), with a small barrier such that the transition state could be crossed more than once before an event check occurs. Moreover, when a trajectory on the boosted surface is confined far from the minimum on the original potential, a relatively long distance needs to be traveled in “quickmin” to reach the minimum. With loose tolerances, “quickmin” times out, the minimization does not exactly reach the minimum, and the atom displacement is shorter than Devent. An event is not identified, and the net rate decreases. Although this may be seen as an impediment, we note that the HD simulations yield correct rates (i.e., the same rate as for diffusion with ΔVmax = 0) for a wide range of ΔVmax. Moreover, the value of ΔVmax where the boosted rate begins to deviate from the exact rate is around the value of the barrier for hopping in this system, which we earlier estimated from climbing-image nudged-elastic bind (CI-NEB) calculations.29 Additional simulation overhead is required for more frequent event checking and for more thorough optimization, so the loose parameters for these routines are justified for the present system. However, these parameters may not be satisfactory for all systems.
Event probabilities for single-adatom diffusion are shown as a function of ΔVmax for different short time intervals in Fig. 6, along with their fits to the sigmoid curve in Eq. (13). Each data point in Fig. 6 is an average over 200 simulations. Here, we see that P has a similar form to the cosine potential, and it is evident that the sigmoid curve is an excellent fit to the data. As we discussed for the cosine potential, appropriate values of ΔVmax can be estimated from the maximum of and/or the minimum of . Table II lists the optimal values of ΔVmax obtained from these derivatives, along with the time-boost factor B for each of value of ΔVmax.
|t (ps) .||.||.|
|ΔVmax .||B .||ΔVmax .||B .|
|t (ps) .||.||.|
|ΔVmax .||B .||ΔVmax .||B .|
From Fig. 5, we see that the ΔVmax predictions listed in Table II for are safe, with values of less than half the barrier29 (∼0.5 eV) estimated from Climbing-Image Nudged Elastic Band (CI-NEB) calculations.30 The ΔVmax predictions from are also safe and are closer to the CI-NEB barrier. The time-boost factors for these two estimates are B = 25 from and B = 41–43 from . To gauge the efficiency, we required 9 h to obtain the 150 ps curve and the prediction range in Fig. 6. The simulations to acquire P were run in parallel on 33 cores in several batches. For regular MD simulations, it will take 135 wall-clock hours to reach 1 μs in this system, running on four cores of Penn State’s Roar cluster. With the time-boost factor range of 20–40 in Table II, a simulation with OptiBoost will require 12–14 h (9 h for selecting a boost and 3–5 h for HD simulations running on four cores). It is evident that the computational effort has moved from executing the simulation to selecting the boost.
Even with the overhead of selecting an appropriate value of ΔVmax, HD simulations are an order of magnitude faster than regular MD. Here, we note that (1) the efficiency of selecting an appropriate value of ΔVmax could be improved by more extensive or complete parallelization; (2) a boost can be selected “once and for all” to benefit future simulations, dramatically increasing the efficiency. Thus, once the boost(s) have been established for a particular system, HD simulations with the bond-boost method would be more than 20 times faster than regular MD.
2. Ag dimer diffusion on Ag(100) in vacuum
Figure 7 shows the physical hopping rate as a function of ΔVmax for Ag dimer diffusion on Ag(100) in vacuum. Each data point in Fig. 7 is an average of five simulations, with each run covering 40–100 ns of MD time [tMD in Eq. (4)]. The observed dimer event rate fluctuates around 0.14 hops/ns for ΔVmax from 0.0 to 0.20 eV and then decreases. Unlike single-adatom hopping and motion on the cosine potential, several different kinds of events can occur for dimer diffusion and the rate in Fig. 7 reflects all these different events. The most frequent dimer motions are twirling dissociation and linear dissociation. We quantified these two events using CI-NEB calculations,30 as shown in Fig. 8. The smallest-barrier event is backward aggregation from twirling dissociation in Fig. 8(a), which has a barrier of 0.24 eV, and we note that the HD simulations in Fig. 7 become inaccurate when the boost is larger than around 0.2 eV. The forward and backward barriers for linear dissociation, indicated in Fig. 8(b), are both larger than this value.
Event probabilities for dimer diffusion beginning with the intact dimer state are shown as a function of ΔVmax for different short time intervals in Fig. 9 along with their fits to the sigmoid curve in Eq. (13). Each data point in Fig. 9 is an average over 200 simulations. As for the previous cases, the sigmoid curve is an excellent fit to the data. Table III lists the values of ΔVmax obtained from and , along with the time-boost factor B. While the ΔVmax from and are smaller than both barriers to break the dimer from the intact configuration, these values are larger than the smallest barrier for dimer recombination after twirling dissociation in Fig. 8(a) and larger than the barrier for which the physical rate deviates from the exact value in Fig. 7.
|tMD (ps) .||.||.|
|ΔVmax .||B .||ΔVmax .||B .|
|tMD (ps) .||.||.|
|ΔVmax .||B .||ΔVmax .||B .|
It is possible to isolate the recombination event in Fig. 8(a) and determine an appropriate boost for this event. Beginning with the final state for twirling dissociation in Fig. 8(a), we determined event probabilities for recombination of the dissociated dimer. The results are shown in Fig. 10 where we again see an excellent fit to a sigmoid curve. The vertical lines in Fig. 10 indicate the optimal values of ΔVmax obtained from and . Based on the barrier for twirling dissociation in Fig. 8(a) and the barrier for which the physical rate deviates from the exact value in Fig. 7, we can see that these are safe values. Thus, applied to a single event, the bold boost estimates from the OptiBoost method can predict safe values for ΔVmax.
For the dimer simulations, 10 min was required to obtain the 5 ps curve and the prediction range in Fig. 10. The time-boost factors were B = 3 and B = 10 based on conservative and aggressive predictions, respectively. For regular MD simulations, it will take 129 h to reach 1 μs in this system, running on four cores of our local facility. With the time-boost factor range of B = 3–10 achieved in Fig. 10, a 1 μs simulation with OptiBoost will require 12–36 h running on four cores using a boost that is safe.
Thus, HD simulations running with a minimal boost are up to an order of magnitude faster than regular MD. Here, we note that it is the current state of HD simulations to run with one boost that is set at the beginning of the simulation based on the lowest energy barrier. It is evident from the two examples for Ag in vacuum that such a mode of running would be sub-optimal since the boost determined for twirling dimer recombination (0.18 eV) is much less than that for single-atom diffusion (0.40 eV) or dimer dissociation (0.44 eV). Our calculations with the dimer indicate that it is possible to associate a particular boost with a particular local atomic environment. Thus, once the boosts have been established for various environments, HD simulations with the bond-boost method could run with high efficiency.
3. Ag adatom diffusion on Ag(100) in solvent
Special concerns apply to HD simulations of Ag adatom diffusion in the system with EG solvent. Namely, the EG solvent is highly active so we must frequently update the bond list between Ag and EG and computational effort is wasted on quenching and re-bonding the system. If the bond list is not updated frequently, some EG bonded with the adatom will displace further than the bond cutoff and the bond will remain highly strained. This will result in a zero-bond bias and no acceleration (B = 1). Figure 11(a) depicts the typical way to group atoms in our HD simulations of surfaces. In this set-up, the bottom two layers of Ag are fixed and the upper four layers of Ag along with the adatom(s) are simulated using HD. We used NVT simulations for EG solvent.
We ran test HD simulations with this protocol, using ΔVmax = 0.30 eV and updating the bond list every 100 time steps. 20 simulations were run for 100 ps each, and no events were detected. Moreover, 99% of the boosted bonds had zero boost, and only 55% of the computational effort was spent on dynamics. Subsequently, 50 simulations were performed in which the bond list was updated every ten time steps. Although no event was detected, the time-boost factor was B ≅ 10, 40% of the boosted bonds had zero boost, but only 10% of the computational effort was on dynamics. By updating the bond list more frequently, larger boosts could be achieved, but a significant fraction of the computational effort was wasted on quenching and re-bonding.
As an alternative to our standard protocol in Fig. 11(a), we implemented the method in Fig. 11(b). In this protocol, the top Ag layer and the adatom are run with NVT simulations (together with the EG), so that Ag-EG bonds are mostly excluded from the HD bond list. 200 HD simulations were performed for 100 ps each using the modified group method with ΔVmax = 0.30 eV. The HD bond list was updated every 1000 time steps. In this protocol, an event occurred in one simulation out of nine. Only 10% of the boosted bonds had zero boost, and over 90% of computational effort was spent on dynamics. With time-boost factors of B ≈ 80, the modified group method can be applied to systems with a boundary between active and relatively inactive parts.
In long-time simulations based on the method in Fig. 11(b), the observed event rate as a function of ΔVmax is shown in Fig. 12. Each data point in Fig. 12 is an average over five simulations with times ranging from 40 to 200 ns. Interestingly, while single-atom hopping of the adatom on top of the surface occurred in vacuum, exchange diffusion, in which a surface atom and the adatom exchange places, was the preferred mechanism for a single Ag atom in solvent. We were able to determine that this was not a result of using the boosting protocol in Fig. 11(b) because HD simulations of single-atom diffusion in vacuum still exhibited hopping using this protocol. Additionally, the lattice parameters were the same in vacuum and in solvent, so the exchange mechanism in the solvent systems did not result from strain. Comparing Figs. 8 and 12, we see that Ag adatom diffusion in solvent environment is almost an order of magnitude slower than that in a vacuum environment. From Fig. 12, we see the simulated rate becomes inaccurate for ΔVmax = 0.4 eV.
To test the OptiBoost method, we conducted 200 HD simulations each for two different time intervals and various ΔVmax to obtain the event probabilities shown in Fig. 13. In Fig. 13, we see that the event probabilities for diffusion on a surface in EG solvent exhibit the same sigmoid form as the other systems. Using the OptiBoost analysis of and , we find similar boost ranges as we saw previously: For the conservative analysis with , we find a safe boost of ∼0.3 eV for both time intervals. Using the “bold” analysis from , we predict a safe, but aggressive boost of ∼0.4 eV. From the inset to Fig. 12, the time-boost factor for ΔVmax = 0.4 eV is B = 230. Interestingly, even though the boosts are similar for Ag diffusion in a solvent and in a vacuum environment, B is about 40 times greater in the solvent environment. We attribute these differences to differences in the boosting mode for the solvent system [Fig. 11(a) in vacuum vs Fig. 11(b) in solvent], which leads to a greater fraction of non-zero boosts.
Regarding the efficiency of these simulations, 80 h were required to obtain the 150 ps curve and the prediction range for ΔVmax. The time-boost factors were B = 100 and B = 230 based on conservative and aggressive predictions, respectively. For normal MD simulations, it would take 1316 h to reach 1 μs in this system, running on four cores. With the time-boost factor range of 100–230 indicated by Fig. 12, a simulation with OptiBoost will require 86–91 h running on four cores using a safe boost—an acceleration of over an order of magnitude, even with the substantial overhead of establishing an appropriate value of ΔVmax. Without the overhead of determining ΔVmax, the acceleration could be two orders of magnitude.
We developed a method called OptiBoost to adjust the value of ΔVmax in HD simulations based on the bond-boost method. The OptiBoost method involves (1) selecting a value of ΔVmax and running a series of short trajectories to obtain the probability that an event occurred; (2) fitting to a sigmoid function and; (3) delineating the boost range from the maximum in and from the minimum in .
We demonstrated the OptiBoost method in simulations on a cosine potential, where exact rates were known, and exact analysis could be performed. These results showed that the OptiBoost method could determine safe and aggressive boosts. We then implemented OptiBoost in three different simulations involving Ag atom diffusion on Ag(100): (1) an Ag adatom on Ag(100) in vacuum; (2) an Ag dimer on Ag(100) in vacuum; and (3) an Ag adatom on Ag(100) in EG solvent. We implemented the OptiBoost method for the Ag diffusion systems in the LAMMPS code driven by a Python wrapper. All trajectory analyses for the Ag diffusion systems were performed in Python.
For all the Ag diffusion systems, we obtained similar results to the cosine potential, demonstrating the robust nature of OptiBoost. Even though OptiBoost incurs computational overhead, the increases in efficiency enabled by this approach were typically an order of magnitude. Since determination of the event probability is an “embarrassingly parallel” calculation, the overhead associated with this calculation would become negligible on a massively parallel architecture, allowing for more substantial acceleration. Additionally, there are several aspects of this algorithm that could be optimized, including the length of the run to determine the event probability and the range of ΔVmax values. While we did not optimize these aspects here, this would be a worthwhile future goal. Perhaps most importantly, our simulations with the dimer showed that an appropriate boost can be established “once and for all” for each local atomic environment and recalled once the local configuration is re-visited, so that HD simulations can be run with multiple boosts, free of the overhead incurred to obtain the boost. In this mode of operation, HD simulations would be particularly efficient.
This work was funded by the National Science Foundation (Grant No. OAC-1835607).
Conflict of Interest
The authors have no conflicts to disclose.
The data that support the findings of this study are available from the corresponding author upon reasonable request.