An algorithm for estimating kinetic parameters of atomistic rare events using finite-time temperature programmed molecular dynamics trajectories☆
Introduction
For a large number of material systems, the underlying dynamics can be described in terms of transitions from one long-lived configuration (or state) to another [1], [2], [3], [4], [5]. In order to build an accurate picture of the dynamical evolution of such materials, one requires a detailed understanding of the states of the system and the kinetic pathways from these states [6], [7], [8], [9], [10], [11]. In this paper, we focus on the temperature program molecular dynamics (TPMD) method [12], [13], [14] for performing this task.
TPMD is an efficient approach for seeking thermally-activated kinetic pathways that are rare at timescales accessible with MD (usually MD can access nano- to micro-seconds). In TPMD, transitions are sought with MD by sampling a number of pre-selected temperatures. A temperature program specifies the duration for each temperature stage (, being the number of temperature stages) as the temperature is increased in a stepwise manner from to values as described by Eq. (1). See Fig. 1 for an example of a temperature program. Once a sufficiently high temperature is reached, transitions along the thermally activated pathways from the state of interest become more frequent so that they can be witnessed with MD. The main assumption in TPMD is that the random transition events are first-order processes. The TPMD method provides a framework to analyze transitions from the state as witnessed in several hundreds to thousands of independent MD trajectories, and statistically estimate the rate constants on the basis of the transitions observed, e.g., using maximum likelihood estimate. Rate constants can be estimated at low temperatures where transitions are rare. This results in a significant computational speed-up over standard MD calculations. TPMD has been successfully applied with metals [15], semiconductor [14], ionic materials [16], diffusion processes in presence of solvent [17], and even with biomolecular systems [18], [19]. Temperatures beyond are avoided since the basic assumptions of TPMD, i.e., state in question is long-lived and that transitions are first-order process, may be violated at high temperatures (more discussion on this is provided later). In most trajectories, a transition is witnessed in an early temperature stage as long as is kept reasonably large in Eq. (1). To lower the computational cost, one typically selects .
Despite these advantages, there is scope for improvement in the TPMD algorithm. Previous TPMD implementations required a transition with each TPMD trajectory. However, the probability that transition may not materialize in few MD trajectories (out of several hundreds to thousands) till time can be significant particularly when activation barriers are large. To address this issue, the MD trajectory is allowed to continue at temperature with the hope that soon a transition may occur (see Fig. 1a). A priori it is not known when a transition may happen, and on-the-fly detection of transitions is mandatory. Implementation of the TPMD algorithm with standard MD packages is not straightforward because of the requirement of on-the-fly detection of transitions. Moreover, there is always a chance that with a handful of trajectories one will witness a transition at very long times. This affects the overall use of the TPMD method with standard MD packages.
In this paper, we introduce a new variation of the TPMD algorithm, called the finite time TPMD (FT-TPMD) algorithm, that can be used with standard MD packages, e.g., LAMMPS [20], NAMD [21], Gromacs [22], GULP [23], DL-POLY [24], in a straightforward manner. In FT-TPMD, trajectories are of a fixed duration, namely, . On-the-fly detection of transition is no longer mandatory, and identifying transitions by post-analysis is made possible. Trajectories that do not yield transitions contain useful kinetic information, and a new feature is that during analysis such trajectories are combined with other trajectories that do yield transitions. Rate constants and Arrhenius parameters are estimated from a collection of trajectories even when a sizeable percentage of MD trajectories have not yielded a transition. Implementing FT-TPMD with an on-the-fly transition detection scheme will result in better performance than with post-processing. Most importantly, we find that FT-TPMD can provide better performance than previous TPMD algorithms in some cases, while equivalent performance in others. Therefore, the use of FT-TPMD is recommended.
Section 2 provides an overview of the FT-TPMD method. Section 3 provides validation for FT-TPMD. FT-TPMD and previous TPMD algorithms are compared in Section 4. Section 5 provides important details regarding implementation of the FT-TPMD algorithm. In Section 6, we demonstrate the application of FT-TPMD to the study of Si-vacancy diffusion. Conclusions are listed in Section 7.
Section snippets
TPMD using finite time trajectories
TPMD can be used when the underlying potential energy surface may be smooth (e.g., metal nanoparticle in vacuum [14]) or corrugated (e.g., metal surface in contact with a solvent [17] or biomolecular systems [19]). In both cases, it is required that a long-lived state, which corresponds to a minimum in the free energy surface, can be defined in terms of all-atom coordinates or a sub-set of coordinates. Additionally, TPMD has been applied to examples where states have been lumped into a
Validation of the FT-TPMD algorithm
To validate the FT-TPMD analysis approach, we employ a small kinetic network model with the rate parameters being known a priori. A collection of FPTs is first generated. Next, the rate parameters are estimated without invoking the Arrhenius approximation using Eq. (3), as well as with the Arrhenius approximation using Eqs. (6)–(7).
In the network model of Fig. 3a, states and are accessed from state with rate constants and , respectively. The Arrhenius expression for the rates is
Comparison between finite time TPMD and previous TPMD algorithm
Similar to Section 3, we use a prototype kinetic network model. Two examples are considered. See inset of Fig. 4a and d. The goal is to compare the performance of previous TPMD and FT-TPMD approach when a given number of trajectories are performed. Note that in the analysis of Fig. 4, it is assumed that transitions are detected on-the-fly.
When the activation barriers for kinetic pathways from state 0 are small, such that transitions happen well before the time is reached, almost all
How to perform MD calculations with a temperature program
An MD simulation in the canonical ensemble is performed as follows:
Step 1. Thermalization step: Equations in Section 2 require that the system is in thermal equilibrium at all times in the MD calculation. A short MD calculation is performed with the starting configuration in state at temperature . Initial atom velocities are randomized. Usually, the system is thermalized within few picoseconds. We use the Nosé–Hoover thermostat implemented in LAMMPS [20]. At the end of the thermal
Application of FT-TPMD to Si vacancy diffusion
Kinetic pathways and associated rate parameters are determined using the new TPMD algorithm for a prototype system. Diffusion of Si-vacancy in bulk Si is studied. Based on the observations made advantages of the TPMD method will be discussed in the next section.
Vacancy-assisted diffusion mechanism in silicon (Si) has been extensively studied because of its importance to the semiconductor devices [30], [31]. Here we study the mechanism using TPMD. A periodic crystalline silicon system comprising
Conclusions
In general, the TPMD method shares many advantages with other dynamical rare-event simulation techniques, viz., kinetic pathways can be discovered without any prior knowledge, and dynamical corrections are included in the rate estimate. In addition, previously several advantages of the TPMD method were demonstrated [12], [13], [14]. Orders-of-magnitude acceleration over MD is achieved while largely retaining the simplicity of MD calculations. Implementation of the temperature program
Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Saurabh Shivpuje, Manish Kumawat, Abhijit Chatterjee
Acknowledgment
AC acknowledges support from Science and Engineering Research Board, India , Department of Science and Technology Grant Nos. EMR/2017/001520 and MTR/2019/000909.
References (35)
- et al.
Comput. Mater. Sci.
(2001) - et al.
Surf. Sci.
(2018) J. Comput. Phys.
(1995)- et al.
Comput. Phys. Comm.
(1995) J. Phys. Chem. Solids
(1957)- et al.
Comput. Mater. Sci.
(1998) - et al.
Nucl. Instrum. Methods Phys. Res. B
(1995) - et al.
J. Phys. Chem.
(1996) Stochastic Processes in Physics and Chemistry
(1992)Handbook of Stochastic Methods
(1985)
Theory of Rate Processes
Reaction Rate Theory and Rare Events
Annu. Rev. Mater. Res.
J. Chem. Phys.
J. Chem. Phys.
J. Mater. Res.
MRS Commun.
Cited by (2)
TPMD toolkit: A toolkit for studying rate processes using molecular dynamics trajectories and performing temperature programmed molecular dynamics calculations
2022, Computer Physics CommunicationsCitation Excerpt :This paper provides an implementation of the algorithms discussed in Ref. [16] to accomplish these tasks. In other words, Ref. [16] and this paper complement one another. TPMD offers several advantages including orders of magnitude computational speed-up over standard MD, ability to handle corrugated energy landscapes, it incorporates dynamical corrections [16], no prior knowledge of kinetic pathways and reaction coordinates is required, and uncertainty in the rate estimates can be calculated [19].
- ☆
The review of this paper was arranged by Prof. D.P. Landau.