93, VOLUME NUMBER 17
PHYSICA L R EVIEW LET T ERS
week ending 22 OCTOBER 2004
Learn on the Fly: A Hybrid Classical and Quantum-Mechanical Molecular Dynamics Simulation
Gabor Csa nyi,1 T. Albaret,3 M. C. Payne,1 and A. De Vita 2,3
Cavendish Laboratory, Madingley Road, Cambridge, CB3 0HE, United Kingdom Physics Department, Kings College London, Strand, London WC2R 2LS, United Kingdom 3 INFM-DEMOCRITOS National Simulation Center and Center of Excellence for Nanostructured Materials (CENMAT), University of Trieste, Trieste, Italy (Received 18 May 2004; published 19 October 2004)
2 1
We describe and test a novel molecular dynamics method which combines quantum-mechanical embedding and classical force model optimization into a unied scheme free of the boundary region, and the transferability problems which these techniques, taken separately, involve. The scheme is based on the idea of augmenting a unique, simple parametrized force model by incorporating in it, at run time, the quantum-mechanical information necessary to ensure accurate trajectories. The scheme is tested on a number of silicon systems composed of up to 200 000 atoms.
DOI: 10.1103/PhysRevLett.93.175503 PACS numbers: 62.20.Mk, 02.70.Ns, 31.15.Qg, 71.15.Pd
Many processes in materials involve strongly coupled subregions of the system which are best modeled in a multiscale fashion, i.e., by the concurrent use of different levels of description. As an example, the mechanical failure of brittle solids involves complex bond breaking processes at the advancing crack tip as well as chemically simpler stress concentration phenomena in a much larger region surrounding the tip [1,2]. Using a unique model, e.g., a quantum-mechanical Hamiltonian or a classical force eld, in this or a similar situation is not a viable strategy. On the one hand, performing expensive quantum calculations on the whole system can rarely be afforded, as the system sizes and simulation times treatable in this way are severely limited. Even when possible, such calculations imply a signicant waste of resources, since a simple classical interatomic potential would be sufcient to model most of the system. On the other hand, no classical potential can be used as the only tool to model the entire system unless it can be made accurate enough to describe the whole range of expected chemical reactions. Much effort has recently been devoted in this direction to develop reactive classical force elds capable of high accuracy for some specic materials and classes of reactions [36]. However, it is clear that considerable effort and subtlety are involved in adapting an existing potential [7] or developing a new one each time a new chemical situation is targeted. Indeed, the lack of straightforward transferability is the main disadvantage of the specialized classical potentials techniques, in striking contrast with the application of rst principles techniques, which are computationally intensive but accurate and general. To date, developing a general analytic interatomic potential which is a priori capable of chemical accuracy in all situations still seems an overwhelmingly difcult task [7]. As a result, the use of multiscale modeling appears to be necessary to approach the crack (propagation) problem mentioned above and many others. Consequently, multiHamiltonian schemes have been proposed, which divide 175503-1 0031-9007=04=93(17)=175503(4)$22.50
the system into a quantum region where the more important chemical processes occur and a larger embedding region described by a less computationally intensive model. Matching these different regions mechanically [2,810] and chemically [11,12] is, however, a delicate issue. Perhaps the main practical drawback of these embedding approaches is that dedicated techniques have to be developed and judicious choices made [13] to represent the boundaries between regions. In the present work, we retain the idea that a unique, classical force model should be used to describe the entire system to avoid boundary problems. Also, we aim at constructing a scheme which can be relied upon to be fully transferable, or in all cases no less transferable than some reference quantum model taken, such as an ab initio density functional theory (DFT) scheme. Moreover, we wish to be able to enforce the high accuracy of our reference quantum model where and when we choose to, across the system and during its time evolution. However, we take the view that no feasible analytic form for the classical force model can reproduce simultaneously the reference models accuracy (say, over a small domain of system congurations) and transferability (to all other such domains), since these two properties pose opposite requirements on the force models complexity (measured, for example, by the number of tting parameters). We note that underlying this difculty is a general problem of tting techniques involving parameter optimization, where functional complexity has to be tuned as a compromise between bias and variance of the t [14]. We thus renounce any attempt to develop sophisticated parametrized functional forms with sets of xed-valued optimized parameters. Rather, we select a simple parametrized potential and augment it at run time with the necessary extra information, which is computed on the y by means of accessory quantum calculations. These more expensive quantum calculations are the price we pay to ensure both the transferability and accuracy of 2004 The American Physical Society 175503-1
VOLUME 93, NUMBER 17
PHYSICA L R EVIEW LET T ERS
week ending 22 OCTOBER 2004
electronic-structure calculation in a simple and exible scheme, which is well suited to investigate a range of multiscale processes. These ideas were outlined in [15]. We start by assuming the existence of a parametrized classical potential which captures the general physics of the system studied reasonably well [7,16]. We allow the parameters of the potential to take different values across the system, e.g., with independent parameters ij for each two-body term V2 (ij ; Ri ; Rj ) corresponding to the pair of atoms i and j. The original potential is assumed to accurately describe only processes that are simple from a chemical point of view, but may be taking place throughout the entire system at all times, such as the elastic deformations of the atomic lattice far from an opening crack. We then set out to extend the applicability of the potential to the remaining, more complex processes such as the breaking and forming of chemical bonds which occur during the evolution of the system. The information missing from the classical potential is computed using black box engines [based, e.g., on a DFT or a tight binding (TB) formalism] and then incorporated on the y into the potential by allowing an appropriate subset of its parameters to vary suitably with time. A possible ow structure for the learn on the y (LOTF) hybrid scheme using a force tting [17,18] predictor-corrector approach is as follows: (1) Initialization: start with the physical system in its initial conditions, and with a reasonable choice of initial parameters for the classical potential. (2) Extrapolation: as in standard molecular dynamics (MD), use the potential with xed parameters to generate a small number of time steps of the system trajectory. (3) Testing: in the latest conguration, the local validity of the classical potential is assessed on a site by site basis, and a selected subset of atoms is agged for quantum treatment. (4) Quantum calculations: use any method which provides the desired accuracy to compute the forces on only the selected subset of atoms. (5) Force tting: tune the parameters of the classical potential until they reproduce the accurate forces. (6) Interpolation: return the state of the system to that before the extrapolation and rerun the dynamics, interpolating the potential parameters between the old and the new values. (7) Return to 2. As an example, we consider the diffusion of a vacancy in silicon. The empirical potential used is the StillingerWeber (SW) potential [16]: E
N X X i
1 Bij =rp Aij fcut rij ij 2 cosijk Cijk 2 fcut rij ; rik ; (1)
value (3.77 A). The parameters Aij , Bij , and Cijk are allowed to vary with time. After initialization, in step (2) we extrapolate the dynamics at nite temperature for ve MD steps, using the velocity Verlet algorithm with 1 fs time step. In step (3) we choose to compute accurate (TB) quantum-mechanical forces for the atoms located within 7.5 A of the diffusing defect. A simple way to calculate these forces is to rst construct for each atom a cluster including all its neighbors up to some given cluster radius, here set to 7 A. To chemically terminate the cluster, the Si atoms bound to cluster atoms and located outside the cluster radius are replaced by hydrogen atoms. These are positioned along the direction of each substituted Si-Si bond at the Si-H equilibrium distance from the cluster atoms. The Hellman-Feynman force on the required atom is then computed after direct diagonalization of the cluster system. The forces on the central atom of each cluster (cf. Figure 1), are shown to converge quickly with cluster radius and to be relatively insensitive to the precise termination strategy [19]. We note that since only these forces are kept from the TB calculations, this implies that the system trajectory will not be affected by ne details in representing the boundary. In step (4) the parameters of the classical potential are optimized by minimizing the squared difference between forces given by the classical model and the accurate forces. For atoms that do not have a quantum-mechanical force calculated on them, the current classical force is used as the target of the optimization. Tests show that local parameters related to atoms which are further than a few neighbor shells outside the QM region do not change appreciably with the t. The optimization [step (5)] is therefore limited to the two-body and three-body terms of the potential (1) in volving atoms in a spherical region of 12 A radius centered on the diffusing defect. About 30 conjugate gradient steps in parameter space are sufcient to perform a converged t, maximum the residual deviation from the target forces being typically less than 0:01 eV=A. The LOTF value for the Si vacancy diffusivity [20] obtained from a 1 ns constant temperature (CT) simulation using the above scheme at T 1400 K with a 215 atom periodic system is 2:1 0:4 106 cm2 =s [19]. This compares well with the value calculated using the TB scheme of 2:3 0:2 106 cm2 =s, obtained using the model in [21] [the Stillinger-Weber potential gives 3:4 1:0 105 cm2 =s]. We report in Fig. 2 results for the vacancy diffusivity in Si as a function of temperature, obtained by
Atoms selected for quantum treatment (defect at 3)
1 2 3 4 5
j2hi;ji
X
j;k2hi;j;ki
where the rst sum is over the atoms; the second and third sums are over pairs and triplets of atoms which include atom i. The cutoff functions fcut go to zero exponentially when the interatomic distance reaches a chosen cutoff 175503-2
QM Cluster 1
F1 QM Cluster 5
F5
FIG. 1 (color). Schematic of the QM zone centered on the defect in position 3, and of two of the clusters used by the QM black box to compute the forces used in the t (see text).
175503-2
VOLUME 93, NUMBER 17
4
PHYSICA L R EVIEW LET T ERS
Kwon et al.
90 =
week ending 22 OCTOBER 2004
log10 D [cm /sec]
4.5 StillingerWeber 5
Angle Parameters, C 0.5 0 0.5 5 4 3
2
5.5 Lenosky et al. 6
0.7
0.8
0.9 1000/T [K]
1
Time [ps]
2 1 0 Bulk Angles Surface Angles
FIG. 2 (color). Arrhenius plot of the silicon vacancy diffusivity obtained by the present scheme tted on the y on different TB models (red), compared with the results of fully TB simulations (black) and xed-parameter SW potential (blue).
FIG. 3 (color). Time evolution of angle parameters near a Si[100] reconstructed surface. A representative surface angle is marked in the inset with a black arc. The bulk angle parameters are all for angles centered on atoms of the third atomic layer.
CT simulations using the xed-parameter SW potential, two different TB schemes [22,23], and our hybrid scheme. To be able to perform several fully quantum-mechanical simulations we used a 64 Si atom cubic cell. As the cell is so small, here we compute the exact TB forces to be used in the LOTF t by direct diagonalization of the periodic system Hamiltonian. We note that the large difference between the results obtained by tting the scheme on the two different quantum models is due to the different predictions of these models in accurate quantitative computations, reproduced by the present method. This emphasizes the fact that the present scheme can at best be expected to reproduce the results of the QM model that it is given [24], but can in no way improve its accuracy. To further elucidate how the scheme adapts to the local environment, we report in Fig. 3 the time evolution of the parameters Cijk at a dimerized Si[100] surface at room temperature. In the bulk these angle parameters remain at all times close to their 1=3 equilibrium value (left). However, on the reconstructed surface, the equilibrium angle on the lower side of the buckled dimers (inset) is lowered to almost 90 . The corresponding parameters learn this by switching to a value close to zero, and ip back and forth between zero and 1=3 as the buckling direction varies with time (right). Moving to a problem where a fully quantum approach would be practically unfeasible, we simulate the gliding motion at 900 K of an opposing pair of 30 partial dislocations in Si, using a 4536 atom unit cell. We ag for QM treatment all atoms within 7.0 A of the dangling bonds (undercoordinated atoms) formed during dislocation motion, corresponding to two disjoint QM regions of 100 120 atoms, which follow the diffusing kink of each partial dislocation. On a single 1200 MHz Pentium III processor, the hybrid simulation takes 34 CPU minutes per MD step using a cluster radius of 7 A. For comparison, a MD step takes 11 h using a state-of-the-art linear scaling tight binding code [25], while a conventional 175503-3
direct diagonalization would be approximately an order of magnitude slower. Improving on the accuracy of the classical potentials is crucial to simulate this extended defect. LOTF simulations on this system reveal a bond topology for the free energy minimum at 900 K which differs from that obtained in the same conditions (and at 0 K [26]) from the SW model, with a square of atoms formed adjacent to an antiphase defect (red atom in Fig. 4), and a different defect migration pathway. While the TB model used [23] may still not capture every relevant feature of the phenomenon under study, these results indicate the need to enforce electronic-structure precision. Indeed, in dislocation dynamics the inaccurate canonical phase space sampling provided by classical potentials is regarded as the main source of discrepancy between theory and experiment [26].
FIG. 4 (color). Equilibrium structure of the Si 30 partial dislocation left kink from a 900 K hybrid simulation (only the glide plane atoms are shown). One undercoordinated atom (red) and its neighbors (blue) are agged for QM treatment. The gray atoms further away complete the set of atoms for which potential optimization is performed. The parameters corresponding to interactions between the yellow atoms in the external region are not updated. Empty circles represent the stable 0 K structure [26].
175503-3
VOLUME 93, NUMBER 17
PHYSICA L R EVIEW LET T ERS
720 A [110] [112]
week ending 22 OCTOBER 2004
thank Stefan Goedecker for providing the code that was used to benchmark the linear scaling tight binding method.
FIG. 5 (color). Propagation of the 111110 crack in silicon, using the LOTF scheme and the SIESTA code as quantum engine. Brittle fracture propagation on a (111) cleavage plane is correctly predicted. Note the 2 1 Pandey reconstruction just appearing on the upper opening surface.
We next test the scheme on a much larger system containing approximately 200 000 atoms, to simulate brittle failure in silicon under uniaxial tensile stress. About 300 atoms near the crack tip are treated quantum mechanically using the SIESTA DFT package [27] as the QM...