Computing energy barriers for rare events from hybrid quantum/classical simulations through the virtual work principle

Hybrid quantum/classical techniques can flexibly couple ab initio simulations to an empirical or elastic medium to model materials systems that cannot be contained in small periodic supercells. However, due to electronic nonlocality, a total energy cannot be defined, meaning energy barriers cannot be calculated. We provide a general solution using the principle of virtual work in a modified nudged elastic band algorithm. Our method enables ab initio calculations of the kink formation energy for 〈100〉 edge dislocations in molybdenum and lattice trapping barriers to brittle fracture in silicon.


I. INTRODUCTION
The two-way chemomechanical coupling of chemical and elastic fields creates inextricably multiscale problems with a simultaneous requirement for chemical accuracy and large system sizes. Density functional theory (DFT) has been shown to have excellent predictive power [1], but its typically high O(N 3 ) computational cost limits its application to problems with fewer than around 1000 atoms [2]. This problem is particularly acute for crystal defects such as dislocation lines [3], grain boundaries [4], and cracks [5], which all possess a long-range elastic field that can rarely be contained in small periodic supercells without unrealistically strong image interactions or strain gradients. While linear elastic corrections have successfully removed finite-size effects for small point defect clusters [6] and screw dislocation dipoles in bcc metals [7], in the majority of cases crystal defects require very large supercells which even O(N) first-principles approaches [8] cannot readily accommodate, especially in metallic systems. Furthermore, complex processes such as dislocation emission or thermally activated crack growth occur on time scales that are far too slow for direct dynamical simulations at the ab initio level. As a result it is necessary to determine rare event rates using transition state theory [9], for which the ability to calculate energy barriers is essential, using, e.g., the nudged elastic band (NEB) method [10].
Large systems can be accurately modeled by combining a local QM description with classical models using hybrid multiscale approaches [11] such as the quantum mechanics/molecular mechanics (QM/MM) [12] and "learn on the *  fly" (LOTF) schemes [13]. The LOTF approach has been used extensively to perform ab initio molecular dynamics [5] but is limited to dynamical simulations and cannot compute energy barriers. Energy-based QM/MM schemes for metals developed by Lu and co-workers [12] have been applied to energy pathways [14,15] but do not provide seamless coupling for materials systems [11] and require the definition of system-specific interaction potentials between the QM and MM regions, restricting the generality of the approach. Flexible boundary DFT calculations couple a fully quantum mechanical simulation to an infinite continuum through a lattice Green's function (LGF) [3,[16][17][18][19]. These methods are ideally suited to crystal defect calculations, as the heavily deformed defect core is treated quantum mechanically while the weakly deformed elastic field is captured in the bulk region for comparatively negligible computational cost [20]. However, while elastic embedding methods allow complex local chemical effects to be modeled, they cannot include thermal or entropic effects and rely on the existence of analytical elastic solutions not readily available for complex threedimensional problems. Recent work to numerically compute the lattice Green's function of large-scale defects extends the applicability of the approach [19], but it remains restricted to structural optimization and does not yet allow energy barriers or temperature effects to be modeled. Similarly, the QM coupled atomistic/discrete dislocation method (CADD) approach [21] couples a DFT region directly with a finite-element model but also cannot be used to compute energy barriers.
While hybrid and flexible boundary DFT calculations have been successfully applied to treat a wide range of crystal defects, they suffer from a well-known limitation: due to the nonlocality of the electronic energy an "energy-per-ion" cannot be defined in the QM region, meaning that the total system energy, which in principle should be a sum of classical and quantum contributions, cannot be defined [20]. (We note that energy differences could in principle be calculated by relying on cancellation of errors near boundaries; however, this uncontrolled assumption would have to be tested on a case-by-case basis.) As a result, important and highly desirable quantities such as migration barriers and segregation energies have long been considered inaccessible.
In this paper we detail a general solution to the problem of extracting energy barriers from hybrid simulation schemes without a total energy function. We exploit the fact that ionic forces in both the classical and quantum region are well defined and localized, allowing us to apply the principle of virtual work to construct energy barriers for a given configurational pathway. Combining this principle with the nudged elastic band routine for finding minimum energy pathways allows the calculation of energy barriers in systems much larger than can be treated in periodic DFT supercells. This is related to using thermodynamic integration to reconstruct free energy profiles in biochemical QM/MM methods [22,23], with the key difference that here we target zero temperature potential energies since entropic effects are comparatively small in hard condensed-matter systems. We demonstrate our method on two problems typically considered inaccessible to ab initio methods: kink formation on 100 edge dislocations in Mo and lattice trapping barriers to brittle fracture in Si.

II. HYBRID SIMULATION SCHEME
A prototypical hybrid simulation scheme is shown in Fig. 1. To provide correct forces on atoms in the QM region, at each force call a DFT calculation is performed which contains the QM region, a surrounding "buffer" region, and a vacuum layer to remove periodic image effects. The presence of free surfaces in the DFT supercell induces electronic (though not elastic) surface states, whose effects must be contained within the buffer region, which in practice determines the required buffer width. For insulators dangling bonds are created whose effects can be suppressed through hydrogen bond termination, while in metals a charge dipole is induced with decaying Friedel oscillations [24]. As the buffer region is treated in DFT only to provide correct forces in the QM region, forces on atoms in the buffer (and bulk) region are given by the classical force field, following the "abrupt force mixing" coupling scheme, which gives accurate forces throughout the overall QM/MM system, in contrast to other handshaking methods that typically incur large force errors close to the QM/MM interface [11].
Here, DFT calculations are performed on clusters composed of QM and buffer atoms surrounded by vacuum. Alternative embedding approaches have been proposed for metals that use periodic QM calculations surrounded by bulklike regions instead of vacuum [16,25,26]. For example, Woodward [16] modeled dislocation cores in periodic DFT cells by incorporating a domain boundary at the edge of the cell. This is appropriate where the embedding region is bulklike; however, the topology of dislocations and cracks of interest here restricts the general applicability of such an approach. We therefore used mixed boundary conditions, with periodicity retained in the direction along dislocation lines and crack fronts, and vacuum added in the other two directions.
For hybrid simulation schemes to produce accurate results, the quantum/classical transition region should typically be only weakly deformed by the presence of the defect, such that an interatomic potential with identical elastic properties and lattice constants (B cl ,a cl ) to the DFT system (B qm ,a qm ) would give an identical mechanical response. However, while modern interatomic potentials typically reproduce DFT elastic properties well, the agreement is not perfect; as a result, the atomic positions used for the classical calculation must be scaled by a factor α = a qm /a cl such that atoms in a perfect bulk lattice are fully relaxed in both systems. In addition, using the classical and quantum bulk moduli B cl and B qm to represent the elastic properties of each medium, the classical atomic forces are scaled by a factor αβ, where β = B qm /α 3 B cl . A derivation of this scaling is given in the Appendix. Beyond elastic matching, there are no further QM/MM interaction terms to be calibrated, unlike for energy-based QM/MM schemes where an interaction potential describing the energetic coupling between QM and MM regions must be specified [14].
Our hybrid force mixing implementation was performed in the Atomic Simulation Environment [27], using LAMMPS [28] to generate classical interatomic forces and VASP [29] to perform DFT simulations using projected augmented wave pseudopotentials [30]. To test the force mixing scheme and buffer size we first considered a perfect fcc lattice of aluminum using an embedded atom method (EAM) interatomic potential by Liu et al. [31]. The QM region was a cube of 13 atoms, with a buffer region of width w containing all atoms within a distance w from an atom in the QM region. In this instance the DFT system is a free cluster, meaning only a Ŵ-point calculation is required, with a plane-wave cutoff of 320 eV. As there should be no residual forces on atoms in a perfect lattice configuration, we measured the total magnitude of atomic forces on all atoms in the QM region with buffer size. As shown in Fig. 1, convergence was achieved for a buffer width of 6.5Å, or around three atomic layers, with the total residual atomic force in the QM region being around 10 −3 eV/Å, well below the tolerance of at most 10 −2 eV/Å per atom used during structural minimization.

III. VIRTUAL WORK PRINCIPLE
The virtual work principle states that the energy E(r) required to traverse a pathway U(r) ∈ R 3N , r ∈ [0,1] in configuration space is given by where F(r) ≡ F(U(r)) ∈ R 3N is the force vector for a given configuration U(r). When the force is a gradient field of some energy function V (X) (for which only the spatial gradient, namely, the force, can be calculated in hybrid simulations), it is simple to show that E(r) = V (U(r)) − V (U(0)). We have implemented the virtual work principle in a modified nudged elastic-band-constrained minimization routine [10], evaluating U(r) from a splined set of (possibly unconverged) NEB knots, using (1) to extract energy differences along the pathway. In the NEB routine an energy functional is only required to define the climbing image and, in some variations of the method, to determine the finite difference scheme used to construct pathway tangents. As a result, we first run iterations with no climbing image defined, until a certain tolerance in the maximum force component perpendicular to the pathway is reached, then use (1) to define energy differences along the pathway to identify a climbing image. Typically, a larger number of knots are required as compared to standard NEB calculations to ensure the splined configuration is as smooth as possible; we have found 10-15 knots to be adequate for all the systems considered here. In Fig. 1 we demonstrate an implementation of this method for the migration of a vacancy in fcc aluminum, treated in the hybrid scheme using only a single Ŵ-point calculation, as the DFT region is again a free cluster. We also performed the same simulations in pure DFT using a 3 × 3 × 3 supercell of 107 atoms and a 7 × 7 × 7 k-point grid. In the latter case we are able to extract the total energy and therefore compare the accuracy of our method. As can be seen in Fig. 1, it is clear that the virtual work energy landscape as calculated in hybrid and DFT and the energy landscape as extracted from the total DFT energy are in extremely good agreement, demonstrating the convergence of the hybrid scheme and the validity of the virtual work principle. In contrast to the LOTF scheme, in which the QM region can be moved during a dynamical simulation, here we use the same set of QM atoms for all knots along the NEB path.

IV. MIGRATION OF AN 100 (010) EDGE DISLOCATION IN MO
100 (010) edge dislocations in bcc metals are known to migrate through a double-kink mechanism [32] and play an important role in irradiation damage of bcc metals, forming the core of 100 prismatic dislocation loops [33]. Edge dislocations possess a strong, long-ranged deformation field which, unlike 1/2 111 screw dislocation dipoles [7,34], cannot be contained in periodic DFT supercells. As a result, flexible boundary DFT calculations [3] or the hybrid methods presented here must be used to capture the long-range elastic field. NEB calculations have been successfully applied to calculate the Peierls barrier to rigid dislocation motion in a wide variety of materials [7,20,34,35]. However, in order to correctly calculate a Peierls stress [7,35,36], care must be taken to accurately determine the dislocation core position as a function of the NEB coordinate r, captured through some remapping function x dislo (r). In the present setting the complication of finding a suitable x dislo (r) does not arise, as we focus on the zero stress double-kink formation energy, which controls the thermally activated diffusion of 100 prismatic dislocation loops [37]. By the chain rule one can demonstrate that the maximum energy difference E = max r [− r 0 dr ′ ∂ r ′ U · F(r ′ )] obtained in the virtual work expression (1) is invariant under the substitution x dislo (r) and thus does not affect our results. We note that existing methods [7,20] to calculate the Peierls stress through the determination of a suitable function x dislo (r) can be applied in postprocessing, without modification, to the NEB pathways produced using our approach.
An 100 (010) edge dislocation dipole of length b = |a[100]| was formed in a square supercell such that the dislocation dipoles are separated by half the supercell height, with one dislocation migrating by a[001]. The system was relaxed using a recently developed modified embedded atom method (MEAM) potential by Park et al. [38], which includes an angular dependence to capture the highly directional bonding of bcc metals. We find the MEAM migration barrier converges with increasing system size and dipole separation, as shown in Fig. 2; this size convergence was confirmed in calculations with a single dislocation in a cylindrical supercell, the outermost atoms fixed to the displacements predicted by anisotropic elasticity theory [39]. The size convergence of the migration barrier can also be investigated by considering the local work done by an atom j along the migration pathway where u j ,f j ∈ R 3 are the per-atom values of U,F. With a saddle point r = r s , the locality of the total work can be probed by summing all values of W j (r s ) less than a distance d from the dislocation core, shown in Fig. 4(b). While it is clear that the immediate core region gives the dominant contribution to the migration barrier, the far field is essential to give a convergence result, which is only accessible to the hybrid simulation technique presented here. The final system used for our hybrid simulations consisted of around 10 000 atoms, far too large for a purely ab initio treatment.
In the hybrid simulations illustrated in Fig. 3, the QM region was defined to contain three atomic planes around the joint initial and final positions of the moving dislocation, with the surrounding buffer region constructed as before. Although the DFT simulation has free surfaces normal to the dislocation line, the supercell remains periodic along the line direction, meaning that we must introduce k points in one dimension [2]. Figure 4 shows the result of NEB calculations to determine the Peierls barrier of the dislocation using a variety of buffer widths and total number of k points. Unlike the vacancy and pure bulk systems, where a buffer size of around three atomic planes was required, for the dislocation system we require a buffer of at least five atomic planes leading to DFT clusters containing around 400 atoms, which we attribute to the much greater degree of deformation caused by the dislocation and the lower atomic density of the bcc structure. Nevertheless, across the range of buffer sizes and k points we find a variation from the final converged value of around 10%. Significantly, the value from our hybrid simulations is around 5 times smaller than that found using the MEAM potential, demonstrating the importance of using ab initio forces to treat highly deformed defect cores.

A. Kink formation energy
As the dislocations under study migrate through a kink mechanism, we have estimated the kink formation energy through careful parametrization of the well-known Frenkel-Kontorova (FK) model [32,40]. In the FK model, a regular array of N nodes of spacing a along the dislocation line direction are free to move in the perpendicular glide direction with positions (ia,x i ), i ∈ [0,N − 1]. The nodes are coupled by a line tension of strength Ŵ and a Peierls potential V P (x) = V P (x + b) of period b, which is given in the current setting by the Peierls potential V (r) shown in Fig. 4 using the linear relation The FK model has been successfully applied to calculate the kink formation energy on screw dislocations in bcc metals [32,40], though in the current setting we have found it necessary to allow an additional position dependence in Ŵ(r), giving a total FK system energy where V (r i ) is the migration potential shown in Fig. 4. To determine Ŵ(r) and thus the kink formation energy, we will conjoin two copies of the dislocation core configurations with the same core positions r and calculate the restoring force between when the core positions differ by a small quantity δ. Explicitly, a line profile r i = r + (i − N/2 + 1/2)δ can be formed in atomistic simulations (with N = 2) by conjoining two NEB configurations U(r) ∈ R 3N and U(r + δ) ∈ R 3N  along the dislocation line direction to give an expanded system U ext (r,δ) ∈ R 6N , being a dislocation line twice the original length. In this extended system we can calculate the force to perturb the relative core positions by δ by projecting the force from atomistic simulations F[U ext (r,δ)] ∈ R 6N against the tangent (∂/∂δ)U ext (r,δ) ∈ R 6N , yielding a restoring force The same line profile r i = r + (i − N/2 + 1/2)δ (where again N = 2) can also be constructed in the FK model, which yields a restoring force to order δ of As we have already calculated V (r) through a spline interpolation we can readily calculate the derivatives ∂ r V (r) and ∂ 2 r V (r) and thus can determine Ŵ(r) by setting f FK (r,δ) = f (r,δ), which yields We emphasize that while to extract accurate Peierls stresses a function x dislo (r) which correctly extracts the "true" dislocation position is required, the formation energies calculated using the virtual work technique detailed here are independent of the choice of x dislo (r). We have performed these calculations using both MEAM and hybrid forces to evaluate f (r,δ), yielding a calculation of Ŵ which we show in lower portion of Fig. 4. The FK model (3) can then be used to simulate a much longer dislocation line to obtain a kink formation energy, which in the pure MEAM case can be directly compared to the kink formation energy in molecular statics [32]. This technique yields a kink formation energy of 1.12 eV, which is closely approximated by the MEAM FK kink formation energy of 1.09 eV. Due to the lower line tension and Peierls barrier found in our hybrid simulations, we find a much lower kink formation energy of 0.54 eV, just less than half the MEAM value. It is interesting to note that similar calculations [40] on 111 screw dislocations in bcc Mo, using periodic DFT supercells, find a kink formation energy of 0.52 eV, meaning that both slip systems have similar activation energies for plastic flow [41].

V. BRITTLE CRACK GROWTH IN SILICON
As a final example, we carry out a direct ab initio calculation of the discrete lattice trapping barriers [42,43] to brittle crack growth in silicon in the (110)[110] cleavage system (Fig. 5). We modeled an 8112 atom system with dimensions 297 × 97.9 × 5.43Å 3 , periodic along the crack front direction and with clamped top and bottom edges and applied strains corresponding to strain energy release rates of 5.0, 5.5, and 6.0 J/m 2 (above the Griffith load of 3.44 J/m 2 computed from the relaxed DFT surface energy [44]). An initial configuration with N = 70 broken bonds along the crack line was relaxed using the force-based hybrid scheme using the Stillinger-Weber potential [45] for the MM region and DFT with the Perdew-Burke-Ernzerhof exchange-correlation functional for a QM region containing 32 atoms centered on the crack tip plus a buffer radius of 6Å. The corrugated reconstruction of the (110) surface leads to a slightly blunted crack tip, with an alternating up-down structure that means the next stable minimum occurs with N + 2 broken bonds. Applying the virtual work NEB approach with (N,N + 2) end fracture system, involving a blunt-sharp-blunt tip reconstruction with two bonds opening simultaneously. Inset images show near-tip region for the initial, transition, and final states following hybrid relaxation of the path, with undercoordinated atoms shown in green. Upper-right inset gives dependence of energy barrier on strain energy release rate, with fracture becoming easier as load increases. points identifies a minimum energy path where two diagonally oriented bonds cleave simultaneously, passing through a sharptip transition state. The lattice trapping barrier decreases as the strain energy release rate is increased, predicting thermally activated crack growth rates similar to earlier work where DFT barriers could only be estimated from cluster calculations [44], but with the tip-blunting reconstruction indicating slow crack growth remains important for larger strain energy release rates than previously thought.

VI. CONCLUSION
In summary, we have proposed a method to compute energy barriers for activated processes that combines DFT and classical interatomic potentials in materials systems where strong bonds cross the interface between QM and MM regions. The method has been used to perform an ab initio calculation of the Peierls barrier for 100 (010) edge dislocations in Mo and to identify a novel crack advance mechanism in Si. The method is expected to be generally applicable to any system where localized chemical processes are driven by long-range elastic fields. For example, the technique could be applied to provide ab initio mechanistic insight into the dynamics of three-dimensional crack fronts, where fracture proceeds through kink formation and advance [44], or to provide a QM-based analog of the Rice-Thomson criterion for the transition from brittle cleavage to dislocation emission [46,47]. All data created during this research are openly available from the University of Warwick Research Archive Portal (WRAP) at http://wrap.warwick.ac.uk/91747.

APPENDIX: DERIVATION OF SCALING LAWS
We wish to define a position scaling α and energy scaling β on the classical system to match the bulk moduli and lattice constant of the quantum system. We define a potential energy function E(X), and then a scaled function E ′ (X) = βE(αX).
The corresponding force in the original coordinate system is The equilibrium lattice constant changes from a 0 to a ′ 0 , and the equilibrium cell volume changes from V 0 to V ′ 0 144102-5 according to

The scaled bulk modulus is
Thus if we want to match a target volume V ′ 0 and bulk modulus B ′ we should use where a 0 and a ′ 0 are the lattice constants before and after rescaling. For quantum/classical force mixing, where we label the quantum region as 1 and the classical region as 2, the aim is to rescale the classical region to match the quantum lattice constant a 1 and bulk modulus B 1 , so we have where a 2 and B 2 are the unmodified classical lattice constant and bulk modulus, respectively. The force scaling is thus αβ = B 1 /α 2 B 2 , as given in the main text.