Atomistic algorithms, in which nuclear degrees of freedom are evolved in an underlying potential energy surface (PES), are commonly used to unravel physical phenomena in diverse scientific fields such as materials science, biology, or chemistry, to name a few. This underlying PES may be given by empirical classical interatomic potentials or by fully quantum-mechanical descriptions including electronic degrees of freedom. The main goal of these algorithms is to predict the response of an arrangement of atoms under different external conditions.

Molecular dynamics (MD) simulations evolve the nuclear degrees of freedom according to a set of equations of motion (EOM) fully resolving the atomic vibrations. These EOM involve the forces acting on each atom, which calculation represents an important computational bottleneck of the MD simulations. The more accurate the forces, the larger the computational burden. Also, for the system to be stable in a given PES, a sufficiently small timestep related to the atomic vibrations needs to be applied, which represents a second critical source of inefficiency. However, it often times occurs that the phenomena of interest take place in time scales much longer than atomic vibrations, which control the maximum integration timestep that can be used. This makes straightforward MD inefficient in exploring the PES, such that it might be unable to capture crucial mechanisms in the process to be studied. Typically, MD simulations are limited to sub-microsecond timescales, even at the lowest end of the accuracy spectrum. Overcoming this limitation calls for the development of advanced algorithms to overcome the so-called time scale problem to more efficiently predict long-time evolution or more extensively sample the PES to explore regions in phase space that traditional MD fails to reach.

Two main strategies have been pursued to mitigate the time scale problem following the main computational bottlenecks discussed above: (i) accelerating the calculation of the forces acting on every atom and (ii) bypassing or mitigating the limitations inherent to the short integration timestep by focusing on the slow state-to-state dynamics that control the evolution of materials over long timescales.

In the context of quantum MD (including electronic degrees of freedom), the main strategy so far has been to make the calculation of the forces faster, as accuracy is required to study the problem at hand. As a matter of fact, even in Kohn–Sham density functional theory,^{Reference Kohn and Sham1} which is a widely used mean-field approximation of the quantum many-body problem, the cubic scaling complexity of computing the single electron wave functions makes it prohibitively expensive to consider large-scale calculations. To overcome this limitation, there are many ongoing efforts to develop reduced order scaling approaches for Kohn–Sham DFT.^{Reference Bowler and Miyazaki2,Reference Goedecker3} Another promising approach is the development of orbital-free density functional theory (OF-DFT), where the kinetic energy of non-interacting electrons is modeled as an explicit functional of electron density and the computational complexity scales linearly with system size. The article in this Focus Issue by William Witt^{Reference Witt, del Rio, Dieterich and Carter4} presents a comprehensive overview of OF-DFT, including the development of various kinetic energy density functionals and their accuracy, as well as a wide range of problems where OF-DFT has been adopted.

Another recent development that accelerates the time evolution of a quantum mechanical MD simulation is the Extended Lagrangian Born-Oppeheimer MD. The Extended-Lagrangian Born-Oppenheimer (XLBO) MD formalism has proven extremely efficient for describing the evolution of an ensemble of atoms with quantum accuracy. The XLBO approach relies on the introduction of an auxiliary electronic degree of freedom as a dynamical variable, for which the equations of motion can be integrated without breaking time reversibility, even under approximate self-consistent field convergence. This framework therefore enables accurate energy conserving simulations at a fraction of the cost of regular Born-Oppenheimer simulations.^{Reference Niklasson, Tymczak and Challacombe5–Reference Niklasson and Cawkwell11} While originally formulated for microcanonical trajectories, it has been recently shown^{Reference Martínez, Cawkwell, Voter and Niklasson12} that the XLBO formalism is also compatible with the noise introduced via thermostatting in canonical simulations, maintaining the efficiency of the original formulation.

These developments, among others, enable simulated times larger than their traditional counterparts, which opens the possibility of coupling them with methods to explore state-to-state dynamics like the ones described below. Still, the development of advanced algorithms in this front is crucial to broaden the range of applicability of these methods and gain access to physical problems of interest with unprecedented accuracy.

On the other hand, when the computational expense of the forces is not daunting, methods to overcome the second bottleneck can be applied. In the case of one-ended methods, where the goal is to generate individual trajectories that are representative of the long-time evolution of the system, acceleration techniques tend to fall into one of two broad classes: Accelerated MD (AMD) and kinetic Monte Carlo (KMC) algorithms. The acceleration provided by these methods results in a tradeoff: in exchange for longer times, one has to give up on generating trajectories that are fully continuous in phase space; instead, one only requires that the trajectories be statistically correct at a coarse level. Specifically, these methods produce so-called state-to-state trajectories, where a state typically corresponds to a finite volume in configuration space corresponding to a given metastable configuration of the material. In other words, one might not be able to precisely determine the phase space variables of a system while it explores a given state, but the statistics of jumps from state to state should be preserved. When one is interested in predicting the behavior of materials over timescales that far exceed typical vibrational periods, this is typically all of the information that is needed.

The first class of methods, the AMD techniques, was pioneered by Arthur F. Voter in the late 1990s.^{Reference Voter13–Reference So/rensen and Voter15} The philosophy that underpins AMD is to hang on to the most powerful features of MD: that it can be used to generate statistically correct escape times and pathways without any a priori knowledge of the process. While the inability to do so directly in the conditions of interest is the root cause of the timescale problem, this limitation can be mitigated by instead running dynamics in specially crafted conditions where transitions are faster, and by unbiasing the results back to the target dynamics. This usually entails resorting to various rate theories, but it can also rely on simple limiting behaviors. The main strength of AMD methods is that fairly strong statements can often be made on the quality of the resulting state-to-state dynamics.^{Reference Bris, Lelievre, Luskin and Perez16,Reference Lelièvre17} However, the performance of AMD methods can be variable as it relies on the presence of a large separation of timescales between vibrating within a state and escaping from a state. When this separation is large, AMD can provide considerable acceleration, sometimes into the seconds of simulation time; however, when this separation is limited due to the presence of low barriers, acceleration over MD can be modest.

A contemporary area of interest in the development of AMD methods addresses the question of how to best leverage the massive parallelism provided by current super-computers.^{Reference Perez, Huang and Voter18,Reference Zamora, Perez and Voter19} This is non-trivial, as conventional MD shows notoriously poor strong-scaling, i.e., the increase in simulation throughput at fixed system size quickly saturates as more and more processing elements are added. This is because integrating atomistic EOM is an inherently serial task: one timestep has to be completed before work on the next can start. This means that computation can only be distributed within a given timestep, not across. Of course, more accurate, and hence more computationally intensive, methods can exploit more processing elements efficiently than computationally cheaper empirical potentials; nonetheless, scaling often breaks down before computational resources are fully utilized, especially on leadership-scale computers with hundreds of thousands of cores. This presents an opportunity if these extra resources can be put to good use. Another area of interest is the development of methods that can maintain efficiency even when the separation of timescales between vibrations and transitions is small, either because the system is large, so it contains many local environments that enable transitions at significant rates, or because the system is inherently complex. While a truly general solution to this problem might not exist, there is a great potential to improve the performance in the common case where trajectories are trapped within super-basins, i.e., set of states that are connected by low barriers compared to those that need to be overcome to exit the set. This kind of situation leads to repetitive transitions, which offers opportunities for optimization.^{Reference Perez, Cubuk, Waterland, Kaxiras and Voter20} Both of these aspects are examined in articles found in this Focus Issue, “Long-Time Molecular Dynamics Simulations on Massively-Parallel Platforms: a Comparison of Parallel Replica Dynamics” (p. 813) and “Parallel Trajectory Splicing; Speculation and Replication in Temperature Accelerated Dynamics” (p. 823).

The second class of algorithms focused on the long-term dynamics belong to the kinetic Monte Carlo (KMC) family. The fundamental hypothesis of the KMC algorithm is that the average probability, to first order in time δt, that a particular event α takes place in the next time interval (δt) equals c _αδt, where c _α is a rate constant.^{Reference Bortz, Kalos and Lebowitz21–Reference Gillespie23} This assumption implies that the process studied is Markovian (i.e., the probability of escaping a state in a given time is just a function of the phase space variables at the present time). Gillespie^{Reference Gillespie22,Reference Gillespie23} showed how, from this fundamental hypothesis, a simple algorithm can be derived to follow the state-to-state dynamics of a process characterized by such fundamental relation. The KMC algorithm samples the time evolution of the process with the right probability distribution, giving one realization of the Master Equation. Often times it is the case in physical processes that all the rates for the possible events are not known a priori, or the number of events is so large that makes the solution of the Master Equation daunting. In these cases the KMC poses itself as a feasible alternative for the study of such processes. This is the case, for example, in the context of chemical reactions, kinetic processes of alloys, systems under irradiation, or a combination of such. It is also important to note that the fact that the KMC solves the state-to-state dynamics integrating out the atomic vibrations makes it suitable to study processes ranging a large variety of time scales unlike traditional molecular dynamics models.

The KMC has been extensively used to study chemistry,^{Reference Yang, Sing-Long and Reed24} thermal properties,^{Reference Martínez, Senninger, Fu and Soisson25} surface growth,^{Reference Archarya and Rahman26} or microstructure evolution under irradiation.^{Reference Caturla27–Reference Opplestrup, Bulatov, Gilmer, Kalos and Sadigh36} Different subclasses of algorithms have been developed depending on the specific processes that are to be studied. Object KMC^{Reference Martin-Bragado, Rivera, Valles, Gomez-Selles and Caturla28,Reference Hudson, Dudarev, Caturla and Sutton37,Reference Wen, Ghoniem and Singh38} or Event-driven KMC^{Reference Fu, Torre, Willaime, Bocquet and Barbu39,Reference Oppelstrup, Jefferson, Bulatov and Zepeda-Ruiz40} algorithms, which uniquely model the evolution of certain pre-defined objects, rely on the a priori knowledge of the rates of all possible events. These models are ubiquitous in the study of the behavior of irradiation-created defects. Rigid lattice models (Lattice KMC) assume that the possible events can be described in an underlying lattice and constitute a known set to calculate the rates on-the-fly. These models have been used to study chemistry evolution under thermal and irradiation conditions.^{Reference Soisson32,Reference Soisson and Martin34,Reference Enrique and Bellon41} On the other hand, newer algorithms have been developed that are capable of calculating the event catalog on-the-fly, making no assumptions about what these events might be. The adaptive KMC (AKMC),^{Reference Henkelman and Jónsson42,Reference Vernon, Kenny, Smith and Sanville43} k-ART^{Reference Béland, Brommer, El-Mellouhi, Joly and Mousseau44–Reference Marinica, Willaime and Mousseau46} or the self-evolving atomistic KMC (SEAK-MC)^{Reference Xu and Henkelman47,Reference Xu, Osetsky and Stoller48} are some examples. These algorithms would be able to accurately describe the evolution of a complex process but, in general, are limited in size and physical time compared to the less flexible or accurate previously-mentioned algorithms. It is also worth mentioning that several parallelization schemes have been developed over the years to take advantage of the large-scale parallel architectures commonly used these days.^{Reference Shim and Amar49–Reference Nandipati53}