Efficient and accurate determination of lattice-vacancy diffusion coefficients via non equilibrium ab initio molecular dynamics

We revisit the color-diffusion algorithm [Aeberhard et al., Phys. Rev. Lett. 108, 095901 (2012)] in non equilibrium ab initio molecular dynamics (NE-AIMD) and propose a simple efficient approach for the estimation of monovacancy jump rates in crystalline solids at temperatures well below melting. Color-diffusion applied to monovacancy migration entails that one lattice atom (colored atom) is accelerated toward the neighboring defect site by an external constant force F. Considering bcc molybdenum between 1000 and 2800 K as a model system, NE-AIMD results show that the colored-atom jump rate $k_{\mathrm{NE}}$ increases exponentially with the force intensity $F$, up to $F$ values far beyond the linear-fitting regime employed previously. Using a simple model, we derive an analytical expression which reproduces the observed $k_{\mathrm{NE}}\left(F\right)$ dependence on $F$. Equilibrium rates extrapolated by NE-AIMD results are in excellent agreement with those of unconstrained dynamics. The gain in computational efficiency achieved with our approach increases rapidly with decreasing temperatures and reaches a factor of 4 orders of magnitude at the lowest temperature considered in the present study.

Figure 1

Ratio of the NE-AIMD to AIMD jump rate $k_{\mathrm{NE}}/k_{\mathrm{E}}$ obtained as a function of $F$ at 2000 K. The inset illustrates the linear-fitting range, previously used to extrapolate equilibrium rates [14].

Figure 2

Model used to clarify the colored-atom jump rate dependence on the intensity $F$ of an applied force seen in Fig. 1. The unperturbed ($F=0$) potential-energy profile $E_{\mathrm{LS}0}$ along the diffusion path $x$ is approximated by a sinusoidal curve. With increasing force intensities $F$, the colored-atom equilibrium position $x_{\mathrm{eq}}$ and the transition-state coordinate $x_{\mathrm{TS}}$ move toward each other, whereas the effective jump activation energy $E_a$ decreases monotonically. ${x_{\mathrm{eq}}}_0,x_{\mathrm{TS}0}$, and ${E_a}_0$ are zero-field quantities. $x_{\mathrm{vac}}$ indicates the vacancy position.

Figure 3

Numerical evaluation of accelerated-rates $k$($F$) as a function of the intensity $F$ of the force applied to the colored-atom based on the model illustrated in Fig. 2, and fit of $k$($F$) up to $F^{\prime }\cong 0.75F_{\max }$ [Eqs. (1) and (2)].

Figure 4

Non equilibrium jump rates $\ln \left[k_{\mathrm{NE}}\left(F,T\right)\right]$ obtained as a function of $F$ for $T=1600$, 1800, and 2000 K. The inset illustrates the agreement between equilibrium rates $k_{\mathrm{NE}\to \mathrm{E}}$ extrapolated by non equilibrium $k_{\mathrm{NE}}$ results and equilibrium rates $k_{\mathrm{E}}$ obtained by nonaccelerated AIMD. $\ln \left(k_{\mathrm{E}}\right)$ is shown with corresponding error bars. The force limit used for $k_{\mathrm{NE}}$-result-fitting, as determined in our model calculations (Fig. 3), is $F^{\prime }\cong 0.75F_{\max }$.

Figure 5

Convergence of extrapolated equilibrium rates $k_{\mathrm{NE}\to \mathrm{E}}$, normalized by AIMD rate $k_{\mathrm{E}}$ values with the number of non equilibrium $k_{\mathrm{NE}}$ interpolation points. $k_{\mathrm{NE}\to \mathrm{E}}/k_{\mathrm{E}}$ ratios are obtained for $f=0.7$ and 0.8 eV/Å (two interpolated rates) and $f$ ranging between 0.3 and 0.8 eV/Å (six interpolated rates). The vertical scale is logarithmic (base 10). Shaded red, green, and blue regions with limits marked by arrows (on the right) correspond to $k_{\mathrm{E}}$ error bars at 2000, 1800, and 1600 K, respectively.

Figure 6

Comparison between theoretical and experimental diffusion coefficients $D$ [36]. The error bars on $D_{\mathrm{E}}\left(T\right)[\propto k_{\mathrm{E}}\left(T\right)c_V\left(T\right)]$ and $D_{\mathrm{NE}\to \mathrm{E}}\left(T\right)[\propto k_{\mathrm{NE}\to \mathrm{E}}\left(T\right)c_V\left(T\right)]$ values account only for the experimental uncertainty on equilibrium vacancy concentrations $c_V\left(T\right)$. The table in the inset lists the gain in computational efficiency $t_{\mathrm{E}}/t_{\mathrm{NE}\to \mathrm{E}}$ achieved by our approach, where $t_{\mathrm{E}}$ and $t_{\mathrm{NE}\to \mathrm{E}}$ are approximate simulation times required to obtain well-converged equilibrium vacancy jump rates as a function of $T$ in nonaccelerated AIMD and accelerated NE-AIMD simulations, respectively.