An algorithm for estimating kinetic parameters of atomistic rare events using finite-time temperature programmed molecular dynamics trajectories

https://doi.org/10.1016/j.cpc.2021.107828Get rights and content

Abstract

A new variation of the temperature programmed molecular dynamics (TPMD) algorithm, called finite-time TPMD (FT-TPMD), is presented. In this variation, a collection of independent molecular dynamics (MD) trajectories is generated, such that each trajectory is of a fixed duration of time and a transition from a given state of the system may/may not occur in a trajectory. This eliminates the requirement of a successful transition in each trajectory and on-the-fly detection of transitions needed in previous TPMD implementations. As a result, FT-TPMD can be easily implemented with standard MD packages, such as LAMMPS, for performing TPMD calculations. Compared to previous TPMD algorithm, FT-TPMD is shown to be computationally as efficient in some cases and better in others. In addition, several refinements are discussed that enable systematic detection of rare transition events, classification of states, and calculation of rate constants. A set of rules is devised to identify one or more atoms involved in a transition and disqualify recrossing events.

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 Tα=maxT0+αΔT,Tmax,α=0,1,2,,Nspecifies the duration τ for each temperature stage α (0αN, N+1 being the number of temperature stages) as the temperature is increased in a stepwise manner from T0 to Tmax 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 Tmax 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 30τ100ps.

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 N+1τ can be significant particularly when activation barriers are large. To address this issue, the MD trajectory is allowed to continue at temperature Tmax 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, N+1τ. 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 1 and 2 are accessed from state 0 with rate constants k1 and k2, 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 N+1τ 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 0 at temperature T0. 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)

  • BarkemaG.T. et al.

    Comput. Mater. Sci.

    (2001)
  • ImandiV. et al.

    Surf. Sci.

    (2018)
  • PlimptonS.

    J. Comput. Phys.

    (1995)
  • BerendsenH.J.C. et al.

    Comput. Phys. Comm.

    (1995)
  • VineyardG.H.

    J. Phys. Chem. Solids

    (1957)
  • GilmerG.H. et al.

    Comput. Mater. Sci.

    (1998)
  • GilmerG.H. et al.

    Nucl. Instrum. Methods Phys. Res. B

    (1995)
  • TruhlarD.G. et al.

    J. Phys. Chem.

    (1996)
  • Van KampenN.G.

    Stochastic Processes in Physics and Chemistry

    (1992)
  • GardinerC.W.

    Handbook of Stochastic Methods

    (1985)
  • GlasstoneS. et al.

    Theory of Rate Processes

    (1941)
  • PetersB.

    Reaction Rate Theory and Rare Events

    (2017)
  • VoterA.F. et al.

    Annu. Rev. Mater. Res.

    (2002)
  • ChillS.T. et al.

    J. Chem. Phys.

    (2014)
  • DellagoC. et al.

    J. Chem. Phys.

    (1998)
  • ChatterjeeA.

    J. Mater. Res.

    (2017)
  • ChatterjeeA.

    MRS Commun.

    (2018)
  • 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 Communications
      Citation 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.

    View full text