I. INTRODUCTION
Given the availability of parallel computing resources, it has become more valuable than ever to design atomistic simulation algorithms that are well suited for distributed computing clusters. When it comes to Molecular Dynamics (MD) simulations, algorithm development has mostly focused on improving weak scaling through spatial decomposition. These efforts have been very effective, allowing current researchers to routinely simulate billions of atoms or more.Reference Kadau, Germann and Lomdahl1,Reference Germann and Kadau2 Unfortunately, spatial decomposition has done very little to extend the maximum temporal reach of MD beyond a few microseconds. This is because MD requires the sequential integration of Newton’s equations, and a femtosecond-scale time step is typically needed to ensure numerical stability.
The Accelerated Molecular Dynamics (AMD) approach to the MD time scale problem is to produce a state-wise trajectory at a faster rate than direct MD. It does so by leveraging the separation of time scales between the numerical integration time step and state-to-state transitions. This approach is particularly well suited for applications where this separation is large, like point-defect evolution in crystalline materials. The performance of an AMD algorithm is often expressed in terms of a computational boost factor, given by WCTMD/WCTAMD, where WCTMD and WCTAMD refer to the required wall clock time (WCT) for direct MD and AMD, respectively.
Since the late 1990’s, three generalizations of the AMD approach have been established: Hyperdynamics,Reference Voter3 Parallel Replica Dynamics (PRD),Reference Voter4,Reference Perez, Uberuaga and Voter5 and Temperature-Accelerated Dynamics (TAD),Reference Sorensen and Voter6 each containing a variety of sub-methods. In all cases, an MD-based approach is used to identify an accurate state-to-state trajectory at an accelerated rate. Of the three general approaches, PRD is the most accurate and the only one relying completely on computational parallelism to obtain a boost. Although PRD offers both parallel scaling and accuracy, other approaches usually require less cpu time to achieve the same boost (but potentially longer wall clock times). By contrast, when the dynamics of a system is dominated by activated transitions that are well described by Harmonic Transition State Theory (HTST),Reference Vineyard7 the TAD method can be used on a single processor to achieve an exponential speedup over MD. For the silver adatom system discussed below, PRD would require more than 103 processors to produce the same boost as serial TAD.
Although TAD is a powerful method for the study of HTST-abiding systems, the serial boost factor can be further improved by introducing computational parallelism. For example, the Speculatively-Parallel TAD (SpecTAD) method improves upon the performance of TAD by concurrently executing high-temperature MD within a dynamically generated list of expected future states.Reference Zamora, Uberuaga, Perez and Voter8 This speculation-based technique, which will be investigated further in this work, has been shown to improve the TAD boost by orders of magnitude and has already been used to illuminate new details of defect clustering behavior in complex oxides.Reference Zamora, Voter, Perriot, Perez and Uberuaga9 With the obvious successes of parallel speculation in mind, the current work is intended to demonstrate that other means of parallelism will often lead to a more efficient use of computational resources.
Here, we present two new parallel TAD methods, both intended to utilize replica-based parallelism to accelerate the execution of high-temperature MD within a TAD-based algorithm. First, we present the Parallel Replica TAD (ReplicaTAD) method, which is designed to utilize every available cpu resource within the current state of the official TAD trajectory. Since all concurrent MD studies are executed in the same state, this parallelism approach can be understood as pure replication; similar to the combination of TAD and PRD. We also present the hybrid SpecTAD method, which is designed to combine the advantages of both speculation and replication into a single algorithm. Using both theoretical predictions and real simulations, we compare all three parallel approaches and ultimately motivate further development of an adaptive hybrid approach.
For all algorithms described in this work, a stochastic Langevin thermostat is assumed to provide temperature control during MD execution.Reference Allen and Tildesley10 All transition checks are to be performed by minimizing the potential energy of the system with respect to the configuration of atoms, X atoms, using a force-based minimization procedure. This transition-checking procedure will be referred to here as quenching.
II. TEMPERATURE-ACCELERATED DYNAMICS
Before discussing approaches to parallelism, we first review the basics of the serial TAD algorithm. This brief introduction is not intended to provide the reader with a comprehensive understanding of TAD but to provide a necessary foundation for the parallel algorithms investigated in this work. For a more thorough discussion, we refer the reader to the original paperReference Sorensen and Voter6 and a recent review article.Reference Zamora, Uberuaga, Perez and Voter8
Consider the 3N + 1-dimensional potential energy surface of a material system: it contains a large number of basins that correspond to different conformations of the system. These basins are often also referred to as states. For any AMD method to perform well, these basins should be separated by barriers that are high compared to k BT, where T is the temperature of interest and k B is the Boltzmann constant. TAD further relies on the HTST assumption that the transition rate from a given state, along some pathway j, is given by
$${k_{\rm{j}}} = {\nu _{{\rm{0j}}}}{e^{\left( {{{ - {E_{\rm{j}}}} \over {{k_{\rm{B}}}T}}} \right)}}\quad ,$$where ν0j is the temperature-independent pre-exponential factor and E j is the transition energy barrier. Note that the transition rate depends on the inverse exponential of the temperature. This allows for the construction of a hypothetical trajectory at a low temperature of interest, T Low, by observing transitions at an exponentially accelerated rate at T High.
As illustrated in Fig. 1(a), TAD can be expressed as a recursive function, where each call to TAD is tied to a specific official state visit. Here, we use the word official to refer to a state that is visited by the hypothetical low-temperature trajectory. Each of these visits will correspond to a minimum-energy configuration, X atoms, and the time when the official low-temperature trajectory began the visit, t Low. Beginning with these initial conditions, we first execute a block of MD to thermalize the system at T High [line 10 in Fig. 1(a)]. Although X atoms is allowed to evolve during this MD block, the official high-temperature time, t High, is left unchanged. This thermalization stage is repeated until no transition is observed during the prescribed thermalization time.
(a) Simple pseudocode description of the TAD algorithm, formulated as a recursive function. (b) The algorithm to be executed by each speculation process in SpecTAD. These processes must wait until an initial state (
$X_i^{atom}$) is provided by the master process.
Once thermalization is complete, and a so-called hot point is prepared, we evolve the trajectory by executing a block of MD and advancing t High [line 14 in Fig. 1(a)]. This stage is repeated until a transition is detected, or the high temperature stop time, t Stop, is reached. When a transition is detected, the exact transition time, t High, is drawn from a uniform distribution of times over the preceding MD block. This high-temperature transition time is then mapped onto a hypothetical low-temperature time using
$$t_{\rm{j}}^{{\rm{low}}} = t_{\rm{j}}^{{\rm{High}}}{e^{ - {E_{\rm{j}}}\left( {{1 \over {{k_{\rm{B}}}{T_{{\rm{Low}}}}}} - {1 \over {{k_{\rm{B}}}{T_{{\rm{High}}}}}}} \right)}}\quad .$$This calculation depends on an accurate definition of the barrier height for the j pathway: E j = E jsaddle − E(X Min), which is determined using the nudged elastic band (NEB) methodReference Jonsson, Mills, Jacobsen, Berne, Ciccotti and Coker11,Reference Henkelman, Uberuaga and Jónsson12 [line 18 in Fig. 1(a)]. Once a transition has been detected and its corresponding t Low has been defined, the previously described stages must be repeated until t High reaches t Stop. At this point, it would be naive to accept the first transition observed at high-temperature because it is possible that a later high-temperature transition will map onto an even earlier low-temperature time. This means that t Stop would need to be equal to t Low to accept an event with certainty. However, to obtain an AMD boost over MD, t Stop is typically defined by
$${t^{{\rm{Stop}}}} = {{\ln \left( {{1 / \delta }} \right)} \over {{\nu _{{\rm{Min}}}}}}{\left( {{{{\nu _{{\rm{Min}}}}{t^{{\rm{Low}}}}} \over {\ln \left( {{1 / \delta }} \right)}}} \right)^{{{{T_{{\rm{Low}}}}} / {{T_{{\rm{High}}}}}}}}\quad ,$$so that the winning transition is accepted with a desired confidence level. Here, δ defines the acceptable level of uncertainty and νMin is the assumed lower bound on the pre-exponential factor.
Before moving on to parallelism, it is important to recognize that the interatomic force calculations (the force calls) account for most of the execution time of any AMD method. Even for a classical interatomic potential, the total run time of a serial TAD simulation is approximately equal to the product of the total number of force calls and the mean execution time of a single force call. As illustrated in the cartoon diagram in Fig. 2(a), this means that both quenching and NEB calculations can account for a significant fraction of the total execution time when many high-temperature transitions are detected before one is accepted.
The hypothetical run time behaviors of (a) TAD, (b) SpecTAD, and (c) ReplicaTAD. The diagrams are intended to be educational aids and are not drawn to scale. In all diagrams: The solid gray rectangles depict many separate blocks of basin-constrained MD, with transition-check quenches not shown between each block. Blue rectangles depict an MD block, and corresponding quench, during which a transition was detected. Orange blocks depict the NEB calculations for new pathways. Green blocks depict communication.
III. PARALLEL ALGORITHMS
Starting with the serial TAD algorithm described in Sec. II, several opportunities for parallel execution are relatively obvious. For applications in which the force calls are computationally expensive, for example, a typical parallel decomposition routine is likely to provide a significant speedup. This approach does not require any modification to the serial TAD algorithm and is highly recommended whenever it is efficient to do so. In fact, the code used for this work allows for spatial decomposition through a library binding to the popular LAMMPS engine.Reference Plimpton13 Although this type of parallelism was not advantageous for the small physical test system used in this work, one should assume that any of the presented algorithms can (and should) be implemented with parallel force call capabilities.
In addition to spatial decomposition, there is also a clear opportunity to accelerate the saddle-point calculations. Since the NEB is performed by minimizing the total energy of an entire chain of configuration images, it is relatively straightforward to parallelize this calculation by decomposing the chain into groups of images. Although this approach should improve performance, the remaining serial parts of the TAD algorithm are likely to prevent the overall speedup from scaling very far. In contrast to the spatial decomposition of force calls, the code used for this work does not include parallel NEB calculations, as the implementation has been left for future work.
Although there are a variety of ways to use parallelism to improve the performance of TAD, we will focus on speculation and replication algorithms for the remainder of this paper.
A. Speculative parallelism
The basic concept of SpecTAD is that for many physical systems, it is typical for a winning transition to be observed by a TAD process long before its corresponding t Stop is actually reached. When there are multiple processors available, an idle core can (and should) be used to immediately begin executing a separate TAD process in what is expected to be the next official state (state B, e.g.). Recall that an official state is one that is visited by the hypothetical low-temperature trajectory. This early execution in state B is speculative in nature because there is still a possibility that the parent TAD process will replace B as the next official state.
An ideal pattern of SpecTAD parallelism is illustrated in Fig. 2(b) for the simulation of a three-transition trajectory with three speculative worker processes. This shows a special case because the first transition observed at high temperature always corresponds to the transition that is ultimately accepted. With this simplification in mind, the diagram still accurately depicts the intended overlap of work between sequential states. This overlap implies that pure speculation can use approximately 
${{{{\bar{w}}_{{\rm{stop}}}}} / {{{\bar{w}}_{{\rm{winner}}}}}}$ cores efficiently, where the 
$\bar{w}$ quantities represent average WCT values and 
${\bar{w}_{{\rm{winner}}}}$ is the average WCT needed to observe the winning transition in each state. Figure 2(b) makes it clear that the number of cores used at any given time in a SpecTAD simulation cannot exceed the total number of official states accepted during the entire simulation. For the simple three-transition example in Fig. 2(b), this means that the parallel speedup must be less than three. When the number of accepted events (N accept) is much larger than the number of available speculation processes (N spec), the total run time can be approximated by
$${\rm{WC}}{{\rm{T}}_{{\rm{total}}}} = \max \left( {{{{{\bar{w}}_{{\rm{winner}}}}} \over {{{\bar{w}}_{{\rm{stop}}}}}},{1 \over {{N_{{\rm{spec}}}}}}} \right){\bar{w}_{{\rm{stop}}}}{N_{{\rm{accept}}}}\quad .$$The SpecTAD implementation used in this work requires two types of processes: a single master process and a predefined number of speculation processes. The purpose of the master process is to build the official low-temperature trajectory and to manage the mapping of speculated work onto the various speculation processes. To understand the SpecTAD algorithm, it is most important to understand the behavior of a single speculation process, whose primary responsibilities are to wait for messages and to execute the SpecTAD function detailed in Fig. 1(b).
A typical SpecTAD simulation begins with the master process sending the initial geometry, X atoms, to an idle speculation process. The speculation will then follow a procedure that is very similar to the serial TAD algorithm. The most important difference is the call to a Speculate function [line 24 of Fig. 1(b)], during which the process must alert the master that it has replaced its best guess for the next official state. Since the SpecTAD function is not recursive, and the speculation process is externally managed, it also becomes necessary for the speculation to alert the master when it has officially accepted the winning transition (line 27). Since it is often the case the first transition observed at high-temperature is not accepted, the master must also be able to interrupt speculation processes, and any dependent child processes, to reassign work. That is, when a speculation process is executing the SpecTAD function in a state that is no longer in the list of expected future states, the master must be able to inform the speculation process to wait for a different state to explore. This implementation detail was accomplished for this work by including message probes that are excluded from Fig. 1(b) for simplicity.
B. Replica parallelism
The goal of ReplicaTAD is to use all available resources to accumulate high-temperature MD time as quickly as possible within the TAD procedure. Like SpecTAD, the ReplicaTAD implementation used in this work uses two types of processes: a single master process and a predefined number of replica processes. Using this structure, the master process is effectively responsible for executing the serial TAD procedure [Fig. 1(a)], with the modification that all MD-based work (steps 8 through 15) is off-loaded to replica processes. The procedure followed by the master can be framed as the recursive ReplicaTAD function defined in Fig. 3.
The recursive ReplicaTAD function executed by the master process in the ReplicaTAD method. Note that the master process must communicate with the replica processes. The inset (blue) shows the procedure executed by each replica process.
While the master process is responsible for executing most of the TAD procedures, the replica processes are used to discover transitions. The idea is very similar to the conventional PRD method. However, in the case of ReplicaTAD, a slightly different procedure is used to leverage the fact that the high-temperature MD is constrained to a single energy basin until the total sum of MD time (across all replicas) has reached t Stop. In contrast to PRD, this means that a replica process does not need to stop after each transition is detected. Instead, the replica processes can follow the more efficient procedure detailed in Fig. 3. Note that, in order for the master process to deduce an accurate t High for each transition, the replica processes must record both the WCT and the local MD time whenever an MD block is either initiated or finished. These pairs of values are stored in a history so that, given the WCT of a reference transition on a different process, the replica has a good estimate of the local MD time when that transition occurred.
A simplified diagram of the ReplicaTAD procedure is illustrated in Fig. 2(c) for a single-transition trajectory with three replica processes. The diagram shows that the master process is responsible for all NEB calculations, meaning that the execution time for an ideal ReplicaTAD simulation with infinite replicas would be approximately equal to the time required to thermalize the system in each state and to calculate all the necessary saddle-point energies. This also implies that a ReplicaTAD simulation using only one replica process (and one master) will achieve better performance than Serial TAD because the NEB and MD work will be performed concurrently. The expected parallel speedup of the MD and transition checks for N rep replica processes is given by
$${s_{{\rm{rep}}}} = {{1 + {n_{{\rm{quench}}}} + {n_{{\rm{therm}}}}} \over {{{\left( {1 + {n_{{\rm{quench}}}}} \right)} / {\left( {{N_{{\rm{rep}}}} + {n_{{\rm{therm}}}}} \right)}}}}\quad .$$Here, n therm represents the expected fractional number of thermalization force calls needed for every MD step and n quench represents the expected fractional number of transition-check force calls needed for every MD step. Since the master process will always be responsible for NEB calculations, the total ReplicaTAD run time can be approximated by
$${\rm{WC}}{{\rm{T}}_{{\rm{total}}}} = \left[ {\max \left( {{{\bar{w}}_{{\rm{NEB}}}}{N_{{\rm{NEB}}}},{{{{\bar{w}}_{{\rm{stop}}}} - {{\bar{w}}_{{\rm{NEB}}}}{N_{{\rm{NEB}}}}} \over {{s_{{\rm{rep}}}}}}} \right) + {{\bar{w}}_{{\rm{therm}}}}} \right]{N_{{\rm{accept}}}}\quad .$$Although ReplicaTAD modifies the traditional PRD procedure to avoid unnecessary synchronization costs, it is important to note that this asynchronous behavior requires special attention to implement. For example, there is no guarantee that the master process will receive transition messages from replica processes in the order that the transitions actually occurred. Therefore, it is possible for the master to reduce a value of t High, and then subsequently receive a transition message corresponding to an even earlier high-temperature time. This can be a problem if the master decides that t High > t Stop because it may still be possible for an out-of-order message to replace the winning transition. For this reason, it is necessary for a replica process to inform the master of any outstanding messages whenever it is asked to participate in a reduction operation.
Another important characteristic of ReplicaTAD is that the error can be very sensitive to the number of MD steps between sequential transition checks (the MD block size). When a replica process detects a transition, it defines the WCT of the transition by drawing a random number from a uniform distribution over the last MD block. Theoretically, this value should be drawn from an exponential distribution, but this would require a good estimate of the expected transition rate (which is rarely known a priori). Since an exponential distribution is not used in practice, the definition of the transition time will always introduce some error that is proportional to the MD block size. In ReplicaTAD, this error is multiplied by the number of replica processes because the transition WCT is used to reduce the total t High across all replicas. Since frequent transitions can clearly reduce the efficiency of the simulation, we find that a more practical way to avoid error is to frequently record the geometry of the system over the course of a single MD block. This way, a more accurate transition time can be refined, without increasing n quench.
One final point about the ReplicaTAD algorithm is that the initial thermalization stage must always be long enough to guarantee that the system is properly prepared at T High. If the system is consistently cold when the thermalization stage is concluded, the observed transition rates will be artificially slow, and the errors will be multiplied by the total number of replica processes. We note that this error is easy to avoid by employing a conservative warm-up period for each new state. It is also important to make sure that the random thermal momenta for an initially minimized state correspond to twice the desired temperature. The common mistake of using T High, rather than 2T High, will require the thermostat to gradually increase the temperature as the system approaches thermal equilibrium because roughly half the initial kinetic energy will be converted into potential energy.
C. Hybrid parallelism
In addition to pure speculation and pure replication, one can also employ a hybrid version of SpecTAD, in which each hybrid speculation process behaves as the master for its own local set of replica processes. The hybrid speculation process follows the HybridSpecTAD procedure detailed in Fig. 4, while the local replica processes still follow the procedure detailed in Fig. 3 inset. The ideal behavior results in each speculation process using many local replicas to reduce the time between NEB calculations. If most of the NEB calls are performed after each winning pathway is observed, then many speculation processes can also be used to overlap the work between sequential states. Although a functioning hybrid SpecTAD code was implemented for this work, we will not present a parametric study of speculation and replication process–count combinations. Instead, we use the hybrid algorithm to demonstrate that neither SpecTAD nor ReplicaTAD is likely to be the optimal choice for a given system, but that a combination of the two will usually be best.
The function executed by the master process within the hybrid SpecTAD method. Note that the replica processes still execute the procedure detailed in the inset algorithm in Fig. 3.
Since the hybrid SpecTAD method aims to combine the advantages of SpecTAD and ReplicaTAD, the theoretical run time is given by
$${\rm{WC}}{{\rm{T}}_{{\rm{total}}}} = \max \left( {{{{{\bar{w}}_{{\rm{winner}}}}} \over {{{\bar{w}}_{{\rm{stop}}}}}},{1 \over {{N_{{\rm{spec}}}}}}} \right)\left[ {\max \left( {{{\bar{w}}_{{\rm{NEB}}}}{N_{{\rm{NEB}}}},{{{{\bar{w}}_{{\rm{stop}}}} - {{\bar{w}}_{{\rm{NEB}}}}{N_{{\rm{NEB}}}}} \over {{s_{{\rm{rep}}}}}} + {{\bar{w}}_{{\rm{therm}}}}} \right)} \right]{N_{{\rm{accept}}}}\quad ,$$where the replica speedup for each speculation process (s rep) is still determined by Eq. (5). Since each speculation process will effectively own N rep replica processes, the simulation will consume a total of N spec(N rep + 1) processors. An important limitation of both ReplicaTAD and hybrid SpecTAD, as they are presented here, is that N rep is always static. When the transition rates remain constant over the course of the simulation, then the optimal number of replica processes may also remain constant. For a very simple physical system, it may be possible to search the space of allowed T High–N rep–N spec combinations to determine an optimal parameter set. However, in most realistic systems, the transition rate will vary from one state to the next, and so the parallel efficiency will also vary. To properly address this limitation, an adaptive resource-management procedure should be used to dynamically allocate how all computational resources are used. An optimal algorithm would most likely resemble that of the Parallel Trajectory Splicing (ParSplice) method.Reference Perez, Cubuk, Waterland and Kaxiras14 Since an adaptive algorithm requires rather significant modifications to the hybrid algorithm presented here, its implementation has been left for future work.
IV. PARALLEL PERFORMANCE
A. Simulation details
All simulations presented in this paper correspond to a simple physical test system containing 301 silver atoms at room temperature (300 K). The geometric configuration corresponds to an isolated adatom on an ideal (100) surface. This system was specifically chosen to match previous investigations of TAD performance.Reference Mniszewski, Junghans, Voter, Perez and Eidenbenz15,Reference Zamora, Voter, Perez, Santhi, Mniszewski, Thulasidasan and Eidenbenz16 The interatomic potential corresponds to an embedded atom method (EAM) formulation,Reference Daw and Baskes17 with the Voter parameterization for silver.Reference Voter and Jacobson18 At room temperature, only two types of activated events are typically accepted during a TAD simulation: (i) An adatom hop with an energy barrier of 0.492 eV and (ii) an adatom exchange with an energy barrier of 0.586 eV. All remaining pathways are rarely accepted, as they correspond to much higher energy barriers (≥0.98 eV).
The presented results also use consistent TAD parameters for accepting transitions, that is, νMin and δ are always set to 5 × 1011 and 0.05, respectively. To ensure that transitions are resolved with sufficient accuracy, transition checks are performed every 500 time steps, with a system snapshot being recorded every 25 steps. Here, each step is 4 fs, meaning that the transition-time resolution is 0.1 ps. To guarantee that the system is properly prepared in each new state, the initial thermalization stage of a minimized geometry is performed for 1500 steps (6 ps). After a transition is detected and a hot point is already available, the thermalization stage is only performed for 500 time steps (2 ps). We note that these block sizes were all chosen conservatively since the intention was to maintain a constant set of algorithm parameters across a range of T High settings and processor counts.
For the chosen physical system, the behavior of an actual TAD run is illustrated for two different high-temperature settings in Fig. 5. Here, the black blocks depict the first 20 accepted transitions, with the unique (ordered) ID of the transition on the x-axis and the total WCT of the simulation on the y-axis. The blue, green, and red blocks represent the same exact run but with the ideal scaling behavior for SpecTAD, ReplicaTAD, and hybrid SpecTAD, respectively. For SpecTAD, this means that the next official state is explored as soon as the pathway leading to it is discovered. For ReplicaTAD, this means that the total time in each state can be approximated by the cost of NEB calculations, plus the initial thermalization stage. For hybrid SpecTAD, both the SpecTAD and ReplicaTAD assumptions hold. Note that the normal and hybrid SpecTAD blocks only represent the WCT needed to make the winning speculation for each state, while the WCT passed between speculating and accepting the winning transition is represented by a line.
The ideal parallel execution of states using SpecTAD, ReplicaTAD, and hybrid SpecTAD at high temperatures of (a) 600 K and (b) 900 K. The TAD curve corresponds to the performance of a real serial simulation, while all other curves correspond to the hypothetical scaling of the same trajectory. For serial TAD (black), and ReplicaTAD (green), each block represents the time needed to accept a single transition. For normal SpecTAD (blue) and hybrid SpecTAD (red), each block represents the time needed to find the winning transition, and the line represents the remaining time needed to reach the ultimate stop time.
B. Results
To compare the performance of SpecTAD and ReplicaTAD, both algorithms were used to simulate 1000 transitions (∼8 ms) in the silver adatom system, using 1 through 118 processors and a T High of 600 K through 900 K. The results of the scaling study are summarized in Figs. 6(a) and 6(b) for SpecTAD and ReplicaTAD, respectively. In both plots, the x-axis corresponds to the total number of processors used for the simulation, while the y-axis corresponds to the average WCT needed to accept each transition over the course of the many-transition simulation. The results are also summarized in Figs. 6(c) and 6(d), where the x-axes correspond to the high-temperature setting, and each curve represents the total processor count. Since the performance is measured in terms of an average run time, it is important to note that the error bars correspond to the standard error of the mean for the specific set of physical processors used to perform the long-time scale simulation. Across all measured runs, the average WCT needed to perform a single force call was approximately 0.72 ms.
Summary of performance results for (a) and (c) SpecTAD and (b) and (d) ReplicaTAD. In (a) and (b), the x-axes correspond to the number of processes used for the simulation, while in (c) and (d), the x-axes correspond to the high-temperature setting. In all plots, the y-axis corresponds to the average WCT needed to accept each of 1000 transitions (∼8 ms) in the silver adatom system. The circular data points represent real simulations, while the dashed lines represent ideal (communication-free) theoretical predictions using Eqs. (4) and (6). In (a), the solid green line represents real hybrid simulations with 2 replica processes belonging to every speculation process. In (b), the solid blue line represents real hybrid simulations with 23 replica processes belonging to each speculation process. In (d), the solid black line represents real hybrid simulations with 11 replica processes and 10 speculation processes.
For the real simulations considered here, every accepted pathway corresponded to an adatom hop or exchange event. This means that a theoretical run-time prediction can treat every official state as identical. In Fig. 6, the dashed-line curves represent simple theoretical predictions that employ this assumption, using Eqs. (4) and (6). All theoretical scaling curves were produced using the itemized force call statistics summarized in Table I. These statistics, which were obtained using serial runs at each of the high temperatures of interest, make it clear that the performance of SpecTAD becomes much better than that of ReplicaTAD when T High = 900 K. At this temperature, each winning transition will be discovered (and its corresponding NEB performed) after completing about 5.33% of the total force calls needed to reach t Stop. This implies that the SpecTAD procedure will perform roughly 94.7% of the total work in parallel. ReplicaTAD, on the other hand, is strongly limited by the serial execution of NEB calculations and can be expected to accelerate less than 87% of the work at 900 K. At T High = 600 K, the story is very different because the NEB calculations account for less than 1% of the total simulation time in serial, while speculation can only hope to perform about 90% of the work in parallel, at best.
Breakdown of the average force-call counts per state visit. Total corresponds to the average total force calls needed for each state visit. N NEB and N trans correspond to the average number of saddle-point searches and transitions detected, respectively.
The measured performance results, shown in Fig. 6, agree reasonably well with Eqs. (4)–(7). However, since the models do not consider communication and synchronization costs, there are some qualitative discrepancies. The most important discrepancy is that the theoretical model for ReplicaTAD ignores master–replica synchronization delays, which leads to a negative error in the run time prediction, an error that is proportional to the number of processes. It is likely that much of this synchronization time can be avoided by implementing a more efficient communication strategy, but it is important to recognize that the theoretical model leads to optimistic predictions. For SpecTAD, there is less communication between processes, although the real performance is clearly worse than the theoretical prediction when SpecTAD is performed on 22 cores at T High ≥ 700 K. For these cases, SpecTAD has a processor count that should be close to the expected scaling limit (note that the theoretical predictions are all flat above 22 cores). Therefore, the performance reduction observed in these runs is most likely due to a delay whenever a speculation process is reassigned to a new state. When there are only a small number of processors available, it is unlikely that a speculation process will be interrupted with a reassignment message. This is because there is already a delay between the time when a pathway is observed and the time when a speculation process becomes available to explore the corresponding state. When there are too many processors available, there are also fewer reassignment penalties because there is always a speculation process available when a new pathway is observed.
In addition to pure speculation and replication, Fig. 6 also shows some results for the hybrid SpecTAD approach detailed in Fig. 4. These results demonstrate that a given SpecTAD or ReplicaTAD simulation can often be improved by using the same number of total cores to perform a combination of speculation and replication. At T High = 900 K, in Fig. 6(a), a hybrid SpecTAD simulation was performed using N rep = 2 replica processes and N spec = 1 − 58 speculation processes. At T High = 700 K, in Fig. 6(a), a hybrid SpecTAD simulation was performed using N rep = 23 replica processes and N spec = 1 − 5 speculation processes. Across a range of T High settings, in Fig. 6(a), a hybrid SpecTAD simulation was performed using N rep = 11 replica processes and N spec = 10 speculation processes. In all three cases, the hybrid simulation led to a clear improvement in performance relative to the best-performing pure algorithm. Some of these hybrid results are also shown as circular markers in Fig. 7(a), alongside (dashed line) theoretical predictions and (solid line) pure SpecTAD and ReplicaTAD results.
(a) A comparison of SpecTAD and ReplicaTAD results (solid lines) with selected hybrid SpecTAD results (circular markers). Dashed lines correspond to the theoretical prediction from Eq. (7). (b) Results of a theoretical sweep over various high-temperatures and N rep–N spec combinations. The open circular marker represents to a real hybrid simulation using the parameters corresponding to the optimal theoretical prediction.
The disagreement between the real and estimated results for the 700 K hybrid simulations shows, once again, that our simple theoretical model does not capture the costs of communication and synchronization. With this limitation in mind, the theoretical model can still be used to perform an approximate sweep of various combinations of T High, N rep, and N spec. Several of the performance curves generated by such a sweep are shown in Fig. 7(b). With the total number of processes being limited to 108, the model tells us that the optimal combination will be N rep = 6 replica processes and N spec = 18 speculation processes at T High = 800 K. While this combination does lead to a better performance than pure speculation or replication, the average WCT to accept each state is almost a factor of two greater than the theoretical estimation. Note that this real data point is shown as a red circle in Fig. 7(b).
V. DISCUSSION
The results presented in this work provide significant motivation for the use of parallelism in TAD. To fully appreciate the advantages of both TAD and its parallel extensions, it is useful to recognize that direct MD would take approximately four weeks to reach the average time required for a single silver adatom transition at room temperature. Using serial TAD with T High = 900 K, the first transition is accepted in approximately four minutes, corresponding to an approximate boost factor of 5 × 103. With speculation introduced, this boost factor further improves to 1 × 105, meaning that the average transition is accepted in less than twelve seconds.
Although the performance of SpecTAD is certainly impressive, the most crucial result of this work is that a mixture of speculation and replication can be used to improve the boost even further. The benefit of trading speculation processes for replication processes becomes even more pronounced as the T High setting is reduced. This is very useful in practice because the behavior of most physical systems tends to depart from HTST at elevated temperatures. For the silver adatom system used here, it is known that the average TAD transition rate will be faster than the HTST prediction whenever T High exceeds ∼500K.Reference Mniszewski, Junghans, Voter, Perez and Eidenbenz15 As the T High setting is varied from 600 to 900 K, the error is expected to vary from approximately 10–30%. Therefore, if the goal of this work was to accurately calculate the adatom diffusion rate at room temperature, we would probably want to choose T High < 600 K. In this case, a speculation-only approach would perform significantly worse than a replication-only approach [as suggested by the dashed black curves in Figs. 6(a) and 6(b)]. Furthermore, the theoretical and measured performances of hybrid SpecTAD tell us that a combination of speculation and replication is likely to get us to the same boost obtainable at T High = 900 K, assuming that sufficient cpu resources are available.
The fact that temperature-acceleration can be traded for parallel processing is an important and exciting result. It tells us that, to a point, we can alternatively think of a ReplicaTAD simulation as temperature-accelerated PRD, that is, we can start with a PRD simulation in mind (T High = T Low), and then we can adjust T High to achieve an optimal balance of acceleration and accuracy. When temperature acceleration is introduced to a PRD simulation in this way, it is likely that the parallel efficiency of the replica processes will be reduced. Although this means that ReplicaTAD will not scale as well as traditional PRD, much of the lost efficiency can be recaptured by the hybrid SpecTAD method. Hybrid parallelism can lead to significantly better scaling because the replica processes become more efficient when they are split between a number of speculative states.
To fully appreciate this last point about hybrid parallelism, it is useful to compare the results plotted in Figs. 6(c) and 6(d). These results suggest that ReplicaTAD is more scalable than SpecTAD when T High is less than 750 K, while the opposite is true when T High is greater than 750 K. However, since we know [from Figs. 6(a) and 6(b)] that both replica and speculation processes will be consistently efficient when there are only ∼10 processes, we can expect a combination of 11 replicas and 10 speculations to perform well across all T High settings of interest. In Fig. 6(d), real performance results for this scenario are shown. These results not only confirm that an intuitive combination of replica and speculation processes is effective but also show that the combination can yield performance well beyond the reach of either pure parallelism approach.
Keeping in mind the obvious advantages of hybrid SpecTAD over the other parallel TAD methods, there are still a few reasons why further algorithm development is motivated by this work. As already mentioned, for realistic systems that exhibit significant state-to-state variation in their activated transition rates, the static-replica approach (used in this work) is likely to limit the parallel scaling. Additionally, the current master–replica synchronization procedure requires conservative parameter choices to avoid multiplicative errors in the total transition rate. With these limitations in mind, the investigation presented here makes it clear that the future development of an adaptive parallel-TAD approach must continue to combine the advantages of both speculation and replication.
Since the silver adatom system is dominated by only two transition pathways, with similar energy barriers (0.492 and 0.586 eV), it is also important to note that the presented comparison between speculation and replication is somewhat limited. In general, the optimal allocation of parallel resources is strongly dependent on the specific distribution of possible transition rates. This is because these rates ultimately determine the 
${{{{\bar{w}}_{{\rm{stop}}}}} / {{{\bar{w}}_{{\rm{winner}}}}}}$ ratio, which in turn corresponds to the efficiency of parallel speculation. For the silver adatom system, the large disparity in energy between the dominant mechanisms, and the more-complex pathways (∼0.6 eV versus ∼1.0+ eV), means that early speculations are especially likely to correspond to winning transitions.
What the simplicity of the chosen physical test system ultimately means for the analysis presented above is that the use of replication is likely to become more advantageous than speculation when the system becomes more complex. For systems with transition rates that are strongly dependent on the current state of the official trajectory, low-energy flickers can certainly degrade the efficiency of pure replication. However, it is also much easier to guarantee that all processes will produce work within an officially visited state. Therefore, in practice, a reasonable approach to the simulation of complex systems is to begin with a replication-heavy allocation of processes and to gradually add speculative character. Since the optimal allocation of resources depends strongly on the exact physical system, this recommended approach represents yet another reason why an adaptive procedure would be valuable.
VI. CONCLUSION
In summary, we have demonstrated that the speculative approach to TAD can often be improved by incorporating replication-based parallelism. To illustrate this point, we have compared three parallel TAD methods: SpecTAD, ReplicaTAD, and hybrid SpecTAD. After reviewing the basics of serial TAD, we discussed the behavior of each parallel algorithm in detail and presented theoretical predictions for their performance. Using a combination of theoretical models and simulation results, we have investigated the capabilities and limitations of each method. For a given number of cpu resources, the results clearly demonstrate that hybrid SpecTAD will outperform both SpecTAD and ReplicaTAD for some N rep–N spec combination. Therefore, our primary conclusion is that optimal performance will typically call for hybrid parallelism.
ACKNOWLEDGMENTS
This work was supported by the United States Department of Energy (U.S. DOE), Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division. LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. DOE, under contract DE-AC52- O6NA25396.
