Upscaling from ab initio to Molecular Dynamics:

Interatomic Potentials and Hybrid Methods

Scientific background

The so-called ``up-scaling'' of materials simulations, from *ab initio*
Quantum Mechanics, to classical molecular dynamics, to mesoscopic approaches
such as phase field or discrete dislocation dynamics, to continuum-level
simulations, presents a number of challenging problems which bring together
a variety of fields in science and engineering. Such Multiscale Modeling
approaches are currently being developed for a wide variety of materials,
including soft matter (biomaterials, polymers, etc.), molecular crystals
(such as high explosives), ceramics, and metals. A key step
in ''bridging the scales'' from electrons to devices is the ability to
perform large scale atomistic simulations (Molecular Dynamics or Monte
Carlo) of the fundamental processes that govern the behavior of the material
with an accurate and computationally efficient description of the atomic
interactions. Despite the enormous advances in *ab initio* methods
such calculations are too computationally demanding to study directly most
of the processes relevant for the modeling of mechanical or electrical
properties, chemistry, etc. Thus, the need for accurate interatomic potentials
is evident.

The most common approach to obtain accurate atomic interactions is fitting
a classical potential to a variety of data coming from

experiments and/or *ab initio *simulations. In order to increase
the accuracy and range of validity of the potentials, this data should
sample a wide variety of atomic conformations (equations of state in a
wide pressure range for different structures, defect energies, etc.). *Ab
initio* methods are the most appropriate means of obtaining such fundamental
information. This ``global'' procedure allows the design of a unique
potential, which can reproduce the properties of a material (plasticity,
phase transitions, etc.) with good accuracy under various conditions [1].

An alternative approach considers, instead, the classical potential
as valid only in a particular localized region of the (P,T) phase

diagram. Using a set of configurations that sample the phase
space region of interest, the parameters of a classical potential are fitted
to reproduce the force on each atom, evaluated from high-quality *ab
initio* simulations. Although this potential is in general non-transferable
to other thermodynamic conditions, there is a gain in the accuracy of the
classical simulations results near the specified thermodynamical state.
This ``optimal potential'' method has been successively employed to investigate
the melting curve and shock Hugoniot of iron [2], and the shock Hugoniot
of tin [3]. An important, but largely unaddressed, question regarding this
approach is how the parameters entering the potential function vary across
the phase diagram, particulary phase boundaries.

Both of these approaches have been successfully applied to the study
of metals, including phase transitions and high pressure

properties, such as those found under shock loading conditions. On
the other hand, additional complications arise when covalent species, such
as the organic compounds prevalent in chemical and biochemical systems,
are considered.

A purely classical description for chemical reactions, the reactive
empirical bond order (REBO) potential, has been developed and greatly extended
by White, Brenner, and their collaborators for model explosives [4] and
hydrocarbons [5,6]. The effect of the bond order term is included is to
make the strength of chemical bonds depend on their local environments,
thus allowing bond breaking and forming. A classical potential has been
developed for explosives based on *ab initio* calculations [7]. Recently,
even more complex potentials (including long-range electrostatic interactions
with environment-dependent charges and bond order dependent valence terms)
were developed and applied for hydrocarbons, energetic materials such and
RDX and HMX (containing N and O in addition to C and H), metal oxides and
metals [8].

However, as these empirical potentials become progressively more complex,
it becomes apparent that more computational effort than necessary is being
applied to ``uninteresting'' regions which are merely spectators, while
still leaving questions about the accuracy in the chemically reactive regions
undergoing bond breaking and forming or large strains. It seems worthwhile
to explore hybrid methods which couple *ab initio* and classical molecular
dynamics techniques; applying high-level treatments in chemically active
regions while using simpler classical potentials for the less important
spectator atoms [9,10]. If desired, the classical potential may even
be modified on-the-fly during the simulation, as in the optimal potential
technique. Although very accurate, these methods are still expensive as
far as CPU time is concerned (dominated by the *ab initio* computation,
even for small reactive regions).