A simple local expression for the prefactor in transition state theory

TL;DR

This paper introduces a computationally efficient method to calculate the frequency prefactor in transition state theory by focusing on a localized active region and applying a correction, reducing computational costs.

Contribution

The authors develop a simple partial phonon DOS approach for the prefactor calculation, avoiding full phonon calculations and enabling easier integration into existing simulation tools.

Findings

The method accurately predicts the prefactor with reduced computational effort.

Convergence improves with larger active regions, balancing accuracy and cost.

Open source implementation available for LAMMPS.

Abstract

We present a simple and accurate computational technique to determine the frequency prefactor in harmonic transition state theory without necessitating full phonon density of states (DOS) calculations. The atoms in the system are partitioned into an "active region", where the kinetic process takes place, and an "environment" surrounding the active region. It is shown that the prefactor can be obtained by a partial phonon DOS calculation of the active region with a simple correction term accounting for the environment, under reasonable assumptions regarding atomic interactions. Convergence with respect to the size of the active region is investigated for different systems, as well as the reduction in computational costs when compared to full phonon DOS calculation. Additionally, we provide an open source implementation of the algorithm that can be added as an extension to LAMMPS software.

Equations

$$\Gamma ={\nu }^{\ast }\exp \left(\frac{-\Delta E}{k_{B}T}\right),$$

$${\nu }^{\ast }=\frac{\prod _{i=1}^{3N}{\nu }_i^B}{\prod _{i=2}^{3N}{\nu }_i^A},$$

$$H^{A}=\left[\begin{array}[]{cc}A&C^{T}\\ C&E\end{array}\right]\text{ and }H^{B}=\left[\begin{array}[]{cc}B&D^{T}\\ D&E\end{array}\right],$$

$$\det H^A$$

$$\det H^B$$

$${\nu }^{\ast }$$

$${\left({\nu }_1^A\right)}^2=\underset{\Vert u\Vert =1}{\min }{u}^T{M}^{-1/2}{\tilde{H}}^A{M}^{-1/2}u=\underset{\Vert M^{1/2}v\Vert =1}{\min }{v}^T{\tilde{H}}^{A}v,$$

$$\det\left(\left[\begin{array}[]{cc}I&0\\ X&L\end{array}\right]\right)=\det\left(L\right),$$

$$\det \left(H_A\right)$$

Full text

Transition state theory, Diffusion prefactor, Effective jump frequency, Atomistic simulation, First-principles calculations

I Introduction

Kinetic processes play a major role in structure formation, microstructure development and materials properties, and diffusion is the underlying phenomenon in many of kinetic processes. The significant progress of atomic level simulations based on first-principles approaches and empirical interatomic potentials paved the way towards determining diffusivities by the use of computer simulations. Additionally, atomic level simulations provide fundamental insight into the atomic mechanisms in diffusion processes, and without relying on fitting to experimental data, they have predictive power for tuning of diffusion rate properties.

A number of previous studies have used electronic structure calculations, model potential or a combination of both to obtain the diffusion coefficient from atomic level simulations. In general, interstitial or point defect jump frequency, which is the dominating mechanics for diffusion in most elemental crystals, is expressed by the equation

[TABLE]

where $k_{B}$ is Boltzmann’s constant, $T$ is temperature, $\Delta E$ represents the potential energy change needed to carry a defect from an initial equilibrium position to a transition state, and ${\nu}^{*}$ is the effective frequency associated with the vibration of the defect along the transition path. The atomic migration energy barrier calculation is more straightforward and can be obtained by using static total energy calculations. On the other hand, calculating the jump frequency needs more effort and depends on the change in the entropy associated with lattice vibration ($\Delta S$), which is due to the constriction or expansion of the diffusing atom’s path during diffusion. One approach is to evaluate the diffusion rate directly from an analysis of the mean squared displacement of the diffusing atom during a molecular dynamics simulation Tang et al. (1997); Mendelev and Mishin (2009). While this method can deliver the pre-exponential factor in Eq. 1, it is only practical at temperatures close to the melting point. Otherwise the diffusion events are too rare to be observed in the molecular dynamics simulation. Milman et al. combined ab initio molecular dynamics and thermodynamic integration to obtain the diffusion entropy change ($\Delta S$) for interstitial diffusion in silicon Milman et al. (1993). Other studies Frank et al. (1996); Sandberg, Magyari-Köpe, and Mattsson (2002); Mantina et al. (2008, 2009) used the harmonic approximation within the transition state theory (TST) and expressed the effective jump frequency as the ratio of the product of normal vibration frequencies of the initial state of atomic migration to that of the non-imaginary normal frequencies of the transition state Vineyard (1957),

[TABLE]

where $\nu^{B}_{i}$ and $\nu^{A}_{i}$ are the normal vibrational frequencies (with $\nu_{1}^{A}$ denoting the omitted imaginary frequency associated with the unstable mode) in the initial basin and in the activated state, respectively, for a system of $N$ atoms. (This expression assumes that the center of mass of the system is held by a fictitious spring. This has the effect of assigning a fixed stiffness to the rigid translation modes whose associated frequencies are then the same in the basin and in the activated states and thus cancel out in the prefactor expression.)

The common approach to calculate the normal vibrational frequencies in Eq. 2 has been the use of direct force-constant calculation Wei and Chou (1992), where the computation cost scales cubically with the number of atoms in the system. This fact limited the previous atomic level simulations of lattice diffusion in solids based on first-principles calculations to system with at most a hundred of atoms Frank et al. (1996); Sandberg, Magyari-Köpe, and Mattsson (2002); Mantina et al. (2008, 2009), which make it difficult to study more complex processes, such as grain boundary diffusion. Relying only on empirical potential models in order to investigate diffusion processes in larger system (thousands of atoms) may not always be a viable alternative, due to complex multi-body bonding behavior. Deploying semi-empirical interatomic potential models by fitting to zero-temperature ab initio data has been performed by Mendelev et.al Mendelev and Mishin (2009), but it requires the task of developing a potential model for various systems. Huang et al. Huang, Voter, and Perez (2013) developed an ingenious linearly scalable approach by approximating the phonon DOS by the kernel polynomial method as an expansion in terms of Chebyshev polynomials to reduce the computational cost of the direct force-constant method.

Here, we use an approach that facilitates transition rate calculations in large-scale systems by replacing the computation and diagonalization of the Hessian of the whole system by the formally equivalent calculation of a much smaller local effective Hessian. For this purpose, we assume that the entire system of atoms can be partitioned into an “active region”, where the kinetic process takes place, and an “environment” surrounding the active region, whose elastic behavior is independent of the atomic geometry of the active region. Our approach is perhaps most related, in its mathematical formalism, to the one of Ref. Binder et al., 2015, although their primary aim is to handle coarse-grained systems rather than fully atomistic simulations. Therefore, our approach entails different methods for calculating the various input quantities needed for the rate calculations.

In this paper, we benchmark this local force-constant approach by applying it to a number of different kinetic processes using realistic energy models (ab initio hamiltonian, ReaxFF or embedded atom potentials): self-diffusion in bulk elemental phases of silver and aluminum, lithium interstitial diffusion in diamond cubic phase of silicon, and oxygen adsorbate surface diffusion on Pt(111) slab.

II Theory

The partitioning of the system into an “active region”and an “environment” (see Fig. 1) implies that the system’s Hessian (after suitable re-ordering of the atoms) can be partitioned in block form as well:

[TABLE]

where $H^{A}$ and $H^{B}$ denote the Hessian in the activated state and in a stable basin, respectively. The block $A$ and $B$ denote the respective Hessians of the active region while block $E$ is the Hessian of the atoms in the environment, assumed to be independent of the configuration in the active region. The $C$ and $D$ blocks represent the coupling between these two regions and can differ in the basin and activated state.

As noted in earlier work,Binder et al. (2015) a standard block-matrix determinant identity implies that:

[TABLE]

(a result shown in Appendix A for completeness). After introducing the mass matrix $M$, we can then write the prefactor as:

[TABLE]

where $\nu_{1}^{A}$ is the imaginary frequency associated with the unstable mode in the activated state. Terms of the form $B-D^{T}E^{-1}D$ and $A-C^{T}E^{-1}C$ can be interpreted as the local effective Hessian of the active region, calculated while allowing the atoms in the environment to relax until they experience no force (rather than being held fixed). These steps closely follow earlier work,Binder et al. (2015) except for the specific way the terms $\left(\nu_{1}^{A}\right)^{2}$, $D^{T}E^{-1}D$ and $C^{T}E^{-1}C$ will be computed in what follows.

In Ref. Binder et al., 2015, the authors seek to obtain a coarse-grained description of the system in which the environment is removed but accounted for via a correction term $-D^{T}E^{-1}D$, which is computed once and then kept fixed throughout the simulation (thus requiring $C=D$). They also aim to approximate $\left(\nu_{1}^{A}\right)^{2}$ only using the coarse-grained representation of the system. In contrast, we do not seek to entirely remove the atoms of the environment, which offers two advantages. First, it allows us circumvent the requirement that $\left(\nu_{1}^{A}\right)^{2}$ be well-approximated in a coarse-grained system. Instead, we simply obtain $\left(\nu_{1}^{A}\right)^{2}$ from

[TABLE]

a quantity that is closely related to what is already needed to find the saddle point via, say, the Nudged Elastic Band Henkelman and Jónsson (2000); Henkelman, Uberuaga, and Jónsson (2000); Sheppard et al. (2012) or the dimer Henkelman and Jónsson (1999) method, again without requiring the calculation of the full Hessian.

A second advantage is to allow us to relax the constraint $C=D$. Relaxing this constraint does not take significantly more computational time because one can use the relaxed environment geometry obtained in the calculation of, say, $H^{B}$ as the starting point for the calculation of $H^{A}$. In a system where $C=D$ is a good approximation, the second calculation will converge much faster.

The approach of Ref. Binder et al., 2015 is most useful when the same partitioning between the active region and the environment can be used for all the transitions of interest while ours is most useful when the transition events can occur anywhere in the system so that the partitioning must be transition-dependent.

In the next sections, we show that the prefactor rates obtained by the above approach agree well with direct force-constant (full Hessian) calculation results, even for small sizes of the active region. Therefore, the use of the aforementioned approach significantly reduces the computation cost by limiting the number of degrees of freedom in force-constant calculations, while delivering accurate results comparable to direct force-constant calculation results. The additional cost of relaxing the atoms in the environment is much less that the burden of full Hessian calculations thanks to the high efficiency of modern conjugate gradient methods in high dimensions. Additionally, the cubically scaling cost of diagonalization of force-constant, which is the computational bottleneck in interatomic potential models, is significantly reduced by narrowing force-constant calculations efforts to atoms in the small active region.

III Computational Details

To test the applicability of the local force-constant method, we employ first-principles calculations in the framework of density functional theory (DFT) for vacancy hopping in bulk aluminum, for which there exist previous ab initio calculations Mantina et al. (2008). However, the finite size effects cannot be excluded in the system sizes that are computationally accessible via DFT, i.e. at most tens of atoms. Therefore, we test the applicability of the model to systems with thousands of atoms by employing interatomic potential models for vacancy hopping in bulk silver, lithium diffusion in diamond cubic silicon, and oxygen surface diffusion on Pt(111) slab. To this end, we have implemented a C++ module in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) framework to add the functionality of local force-constant calculations for any interatomic potential. The source code of this module is open access and is available on the web Kadkhodaei (2018). To invoke this new functionality, one has to simply put our C++ module code in the LAMMPS src directory and re-build it according to the LAMMPS website instructions lam (2018).

For our first-principles calculations, we used the projector augmented wave (PAW) method Blöchl (1994); Kresse and Joubert (1999), as implemented in highly efficient Vienna Ab-initio Simulation Package (VASP) Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996, 1996), with energy cutoff of 300 eV. The generalized gradient approximation (GGA) Perdew (1992) is used for exchange-correlation functional. The Brillouin zone integration is performed using a Monkhorst-Pack k-mesh of $11\times 11\times 11$ for a 31-supercell of fcc aluminum. Tests indicate that this choice of first-principles calculation parameters ensures an energy accuracy of 0.1 meV/atom.

For our interatomic potential calculations, we use an embedded atom (EAM) potential taken from Ref. Williams, Mishin, and Hamilton (2006) to describe atomic interactions for silver and a reactive force field (ReaxFF) with parameters taken from Ref. Fan et al. (2013) to describe interatomic interactions of lithium and silicon. For the ReaxFF potential, we use the parameters derived in Ref. Fan et al. (2013), where they are optimized against a DFT training set consisting of bulk c-$\text{Li}_{\text{x}}$Si and a-$\text{Li}_{\text{x}}$Si phases. A ReaxFF potential with optimized parameters to describe Pt-O system, taken from Ref. Fantauzzi et al. (2014), is used for oxygen diffusion on Pt(111) surface.

The number of atoms in our examples are limited to a couple of thousands atoms for the purpose of comparing our results against the direct force-constant calculation results, while the usage of our method is most advantageous where the number of atoms are in the order of millions and direct force-constant calculations are infeasible.

IV Results

IV.1 Vacancy diffusion

We first demonstrate the local force-constant method by computing the effective jump frequency (pre-exponential factor), $\nu^{*}$, for vacancy hopping in bulk aluminum using DFT calculations. The minimum and saddle point configurations are obtained using the nudged elastic band method as implemented in VASP Henkelman and Jónsson (2000); Henkelman, Uberuaga, and Jónsson (2000). The energy barrier we obtained is 0.52, and when a correction value of 0.05 for the “internal surface” effect Carling et al. (2000) is added to it, the final value of 0.57 eV is in agreement with the vacancy migration energy of 0.57 eV reported in Ref.Mantina et al. (2008). Fig. 2 represents the vacancy migration and its energy profile in fcc aluminum. The atomic configurations are visualized using the VESTA software Momma and Izumi (2011). Having the minimum and saddle point configurations, we calculated the prefactor for different number of atoms in the active region employing the proposed method, as illustrated in Fig. 3. Active region is defined as a sphere originated at the position of migrating atom in the saddle point configuration, as visualized by OVITO software Stukowski (2010) and represented in Fig. 3. A linear increase in the radius of the active region results in a cubic increase in the number of atoms inside the active region. As shown in this figure, our results are compared with a crude local approximation of force-constant, where the effect of atoms in the environment is not accounted for. In other words, in this crude approximation the force-constant calculation is conducted for atoms in the active region while other atoms in the system are fixed. We call this approach the “fixed environment” method, as compared to the proposed method of “relaxed environment”. The results of the “relaxed environment” method indicate a better convergence compared to the “fixed environment” local force-constant calculation. The calculated prefactor values are in perfect agreement with the value reported in Ref. Mantina et al. (2008) that is obtained using the direct force-constant calculation.

Next, we demonstrate our method by computing the effective jump frequency for vacancy hopping in a larger system, a $7\times 7\times 7$ supercell of the conventional cubic cell of fcc silver with a vacancy consisting of 1371 atoms, using the EAM potential model. We used nudged elastic band method, as implemented in LAMMPS Nakano (2008); Maras et al. (2016), to obtain the minimum and transition state configurations. The initial and saddle point configurations for vacancy diffusion in Ag are visualized using VESTA software Momma and Izumi (2011) and are illustrated in Fig. 4. The energy barrier obtained from our calculations is 0.66 eV, in good agreement with the value of 0.65 eV reported in Ref. Huang, Voter, and Perez (2013). As presented in Fig. 4, we calculated the prefactor rate using our method of local force-constant calculation, where the atoms in the environment are relaxed during force-constant calculation of active region atoms. We compared our method with the case of fixing atoms in the environment during force-constant calculations. This comparison indicates that the “fixed environment” approach results happen, in this case, to converge faster due to fortuitous cancellation of errors. As discussed in Supplementary Information, while the partial Hessian calculated in the “fixed environment” approach ($A$ and $B$ in Eq. 3) does not include the effect of atoms in environment, the ratio of partial Hessian determinants ($\frac{\det(B)}{\det(A)}$) happens to be in a good agreement with effective Hessian determinants ratio ($\frac{\det(\tilde{H}^{B})}{\det(\tilde{H}^{A})}$) (see Supplementary Information for a detailed discussion). This case is rather exceptional, however, and the remaining examples we consider below indicate that “relaxed environment” approach more typically performs better. As shown in Fig. 4, even for small sizes of active region, the error magnitude is still fairly small (less than a factor 2), which is typically smaller than errors arising from other sources (such the use of empirical potential models). This fact enables the use of our method for small sizes of active region, which significantly reduces the computational costs involved in vibration frequency calculations.

IV.2 Interstitial diffusion

Here, we demonstrate the effective jump frequency for diffusion of a single Li atom through a diamond cubic supercell of 1728 Si atoms, corresponding to a $6\times 6\times 6$ supercell of 8-atom unitcell, using the ReaxFF potential model. The thermodynamically favorable site for Li atom insertion in cubic diamond Si is a tetrahedral site ($T_{d}$) at $b_{0}$ distance away from a Si atom in the opposite direction of Si nearest neighbor, where $b_{0}$ is the bond length in a cubic Si structure. The diffusion process consists of a Li atom moving between two adjacent interstitial $T_{d}$ sites passing through the hexagonal (Hex) interstitial position. Our calculation results in a lattice constant of 5.46 $\AA$ for c-Si (and a corresponding $b_{0}$ value of 2.37 $\AA$), in good agreement with values reported in Ref. Tritsaris et al. (2012); Ostadhossein et al. (2015), using both DFT and ReaxFF calculations. We calculated the energy of Li atom in $T_{d}$ and Hex positions and the energy difference is 0.63 eV, in good agreement with previous energy barrier values reported in the literature Tritsaris et al. (2012); Ostadhossein et al. (2015). Li interstitial at $T_{d}$ and Hex positions are indicated in Fig. 5 using VESTA. The prefactor rates are calculated with our method and compared with the “fixed environment” approach and the direct force constant calculation results. As indicated in Fig. 5, the discrepancy between our results and full Hessian calculation results remains relevantly small, even for small number of atoms in the active region, confirming the validity of our approach and its effectiveness on reducing the computational cost. We also observe that the “fixed environment” results are very similar to the “relaxed environment” results, slightly underperforming relative to the “relaxed environment” approach.

IV.3 Surface diffusion

We finally represent the effective jump frequency rates for oxygen adsorbate diffusion on platinum (111) surface. We use a 7-layer thick Pt(111) surface with each layer consisting of $4\times 4$ atoms, including 112 platinum atoms in total. This simulation cell size is chosen since there exist previous studies for the diffusion profile for this system size using ReaxFF Fantauzzi et al. (2014). The high symmetry adsorption site on a fcc(111) surface include the fcc, hcp, top and bridge sits, as indicated in Fig. 6. We calculate the prefactor rates for the bridge-diffusion, which consists of an adsorbate moving from an fcc site to the nearby hcp site, passing through the bridge site. To isolate the oxygen adsorbate from the adjacent simulation cells, we included a 30$\AA$ vacuum region in the simulation cell. We conducted nudged elastic band calculations in LAMMPS to obtain the bridge-diffusion path for oxygen surface diffusion on Pt(111). The diffusion energy profile along with the initial, transition and final states are represented in Fig. 6, visualized by OVITO. The calculated forward (fcc-to-hcp) and backward (hcp-to-fcc) energy barriers are 0.65 eV and 0.48 eV, in good agreement with values reported in Ref. Fantauzzi et al. (2014). The prefactor rates are then calculated for different number of atom in the active region, which is a sphere originated at bridge site position. Our results indicate very good convergence to the full Hessian calculation results and the range of error remains as small as 20%. The “fixed environment“ results slightly underperform relative to our approach, although both methods give reliable results..

V Discussion

Our results bring out two key messages. First, one can rigorously define the harmonic TST prefactor in terms of local effective Hessians, where the effect of the environment enters via a simple correction term associated with static elastic relaxations of the environment. Second, we find, through numerical simulations in a wide range of systems, that the corrections due to relaxing the environment tend to be similar in the stable and saddle point configurations, thus justifying an even simpler “fixed environment” approach. This partial error cancellation (that occur beyond what our theoretical model predicts) is analyzed in more detail the Supplementary Information.

These observations open the way to simple, local approaches to computing the TST prefactor. Although we do not explore this possibility in detail here, our observations also suggest efficient hybrid schemes using three different regions: (i) an “active” region whose Hessian is explicitly calculated; (ii) an intermediate region surrounding it where atoms are relaxed during the computation of the effective Hessian of the active region; and (iii) an environment where all atoms are fixed.

VI Conclusion

In this paper, we examined an approach based on a local force-constant calculation in the vicinity of where a kinetic process takes place to obtain the prefactor in transition state theory. While the standard force constant calculation procedure is performed on an “active region”, the overall effect of “environment” atoms is accounted for by relaxing them during this procedure. We demonstrated that the effective jump frequencies obtain using the proposed method rapidly converges to the full force-constant calculation results for a number of diffusion processes, validating the theory employed herein. We also found empirically that a faster and easier local force-constant approach, where the overall effect of atoms in the “environment” is disregarded by fixing them during the local force-constant process, results in very similar prefactor rates due to the cancellation of errors. The proposed approach significantly reduce the computation efforts involved in force-constant calculation by restricting the number of atoms to the vicinity of the kinetic event.

VII Supplementary Information

See supplementary information for a detailed discussion on the cancellation of errors in frequency prefactor calculation for the “fixed environment” method.

VIII Acknowledgments

This work is supported by the Office of Naval Research under grant N00014-17-1-2202 and by Brown University through the use of the facilities of its Center for Computation and Visualization. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.

Appendix A Simple derivation of a block matrix determinant identity.

First observe that

[TABLE]

where $I$ is the identity matrix, [math] is a matrix of zeros, while $X$ and $L$ are conformable arbitrary matrices. We can then write (in the notation of Section II)

[TABLE]

Similarly, $\det\left(H_{B}\right)=\det\left(B-D^{T}E^{-1}D\right)\det\left(E\right).$