Table Of ContentAn Overview of Multiscale Simulations of Materials
Gang Lu and Efthimios Kaxiras
Department of Physics and Division of Engineering and Applied Science,
Harvard University, Cambridge, MA 02138
Multiscalemodelingofmaterialpropertieshasemergedasoneofthegrandchallengesinmaterial
4 science and engineering. We provide a comprehensive, though not exhaustive, overview of the
0 current status of multiscale simulations of materials. We categorize the different approaches in the
0 spatialregimeintosequential andconcurrent, andwediscussinsomedetailrepresentativemethods
2 ineach category. Weclassify themultiscale modeling approachesthat dealwith thetemporalscale
intothreedifferentcategories,andwediscussrepresentativemethodspertainingtotheeachofthese
n
categories. Finally, we offer some views on the strength and weakness of the multiscale approaches
a
J thatarediscussed,andtouchuponsomeofthechallengingmultiscalemodelingproblemsthatneed
tobe addressed in thefuture.
7
]
i I. INTRODUCTION beyond) where a constitutive law governs the behavior
c of the physical system, which is viewed as a continuous
s
- medium. In the macroscale, continuum fields such as
l Some of the most fascinating problems in all fields of density, velocity, temperature, displacement and stress
r
t science involve multiple spatial and/or temporal scales: fields, etc are the players. The constitutive lawsareusu-
m
processesthat occur at a certainscale governthe behav-
ally formulated so that they can capture the effects on
. ior of the system across several (usually larger) scales.
t materials properties from lattice defects and microstruc-
a The notion and practice of multiscale modeling can be
tural elements.
m
traced back to the beginning of modern science (see, for
Phenomena at each length scale typically have a corre-
- example, the discussion in1). In many problems of ma- spondingtime scalewhich,incorrespondencetothe four
d
terials science this notion arises quite naturally: the ul-
n length-scalesmentionedabove,rangesroughlyfromfem-
timate microscopic constituents of materials are atoms,
o tosecondstopicoseconds,to nanosecondsto milliseconds
c andtheinteractionsamongthematthemicroscopiclevel and beyond.
[ (of order nanometers and femtoseconds) determine the
Ateachlengthandtime-scale,well-establishedandef-
behavior of the material at the macroscopic scale (of or-
1 ficient computational approaches have been developed
der centimeters andmilliseconds andbeyond), the latter
v over the years to handle the relevant phenomena. To
3 being the scale of interest for technological applications. treat electrons explicitly and accurately at the atomic
7 The idea of performing simulations of materials across scale,methodsknownasQuantumMonteCarlo(QMC)5
0 severalcharacteristiclengthandtimescaleshastherefore
andQuantum Chemistry (QC)6 canbe employed,which
1 obvious appeal as a tool of potentially great impact on
arecomputationallytoodemandingtohandlemorethan
0 technological innovation2,3,4. The advent of ever more
4 a few tens of electrons. Methods based on density func-
powerful computers which can handle such simulations
0 tional theory (DFT) and local density approximation
providesfurtherargumentthatsuchanapproachcanad-
/ (LDA)7,8 inits variousimplementations, while lessaccu-
t dress realistic situations and can be a worthy partner to
a ratethanQMCorQCmethods,canbereadilyappliedto
m the traditional approaches of theory and experiment.
systemscontainingseveralhundredatomsforstaticprop-
- Inthecontextofmaterialssimulations,onecandistin- erties. Dynamical simulations with DFT methods are
d
guish four characteristic length levels: usually limited to time-scales of a few picoseconds. For
n
(1) The atomic scale (∼ 10−9m or a few nanometers) in materialspropertiesthatcanbemodeledreasonablywell
o
c which the electrons are the players, and their quantum- withasmallnumberofatoms(suchasbulkcrystalprop-
: mechanical state dictates the interactions among the erties or point defects), the DFT approach can provide
v
atoms. sufficiently accurate results. Recent progress in linear
i
X (2) The microscopic scale (∼ 10−6m or a few microme- scaling electronic structure methods9 has enabled DFT-
r ters) where atoms are the players and their interactions based calculations to deal with a few thousands atoms
a
canbedescribedbyclassicalinteratomicpotentials(CIP) (corresponding to sizes of a couple of nanometers on a
which encapsulate the effects of bonding between them, side) with adequate accuracy. Finally, the semi-classical
which is mediated by electrons. tight-binding approximation (TBA), although typically
(3) The mesoscopic scale (∼ 10−4m or hundreds of mi- not as accurate as DFT methods, can extend the reach
crometers) where lattice defects such as dislocations, of simulations to a few nanometers in linear size and a
grainboundaries,andothermicrostructuralelementsare few nanoseconds in time-scale for the dynamics10.
the players. Their interactions are usually derived from For material properties at the microscopic scale,
phenomenological theories which encompass the effects Molecular Dynamics (MD) and Monte Carlo (MC) sim-
of interactions between the atoms. ulationsareusuallyperformedemployingCIPwhichcan
(4)The macroscopic scale (∼ 10−2m or centimeters and often be derived from DFT calculations11,12. Although
2
notas accurateas the DFT andTBA methods, the clas- Conceptually, two categories of multiscale simulations
sical simulations are able to provide insight into atomic can be envisioned, sequential and concurrent. The se-
processesinvolvingconsiderablylargersystems,reaching quential methodology attempts to piece together a hier-
up to ∼ 109 atoms13. The time-scale that MD simula- archy of computational approaches in which large-scale
tions based on CIP can reach is up to a microsecond. modelsusethecoarse-grainedrepresentationswithinfor-
Atthemesoscopicscale,theatomicdegreesoffreedom mation obtained from more detailed, smaller-scale mod-
are not explicitly treated, and only larger scale entities els. This sequential modeling approach has proven ef-
are modeled. For example, in what concerns the me- fective in systems where the different scales are weakly
chanical behavior of solids, dislocations are the objects coupled. Thecharacteristicofthesystemsthataresuited
of interest. In treating dislocations, recent progress has forasequentialapproachisthatthelarge-scalevariations
beenconcentratedontheso-calledDislocationDynamics decouple from the small-scale physics, or the large-scale
(DD) approach14,15,16,17 which has come to be regarded variationsappearhomogeneousandquasi-staticfromthe
as one of the most important developments in computa- small-scale point of view. Sequential approaches have
tional materials science and engineering in the past two alsobeenreferredtoasserial,implicitormessage-passing
decades18. Such DD models deal with the kinetics of methods. The vast majority of multiscale simulations
dislocations and can study systems with a few tens of that are actually in use are sequential. Examples of
micronsinsize andwith amaximumstrain∼ 0.5%for a such approaches abound in literature, including almost
strain rate of 10 sec−1 in bcc metals19. all MD simulations whose underlying potentials are de-
rivedfromelectronicstructuretheory23,24,usuallyabini-
Finally, for the macroscopic scale, finite-element (FE)
methods20 are routinely used to examine the large- tio calculations11,12. One frequently mentioned25,26 ex-
ample of sequential multiscale simulations is the work
scale properties of materials considered as an elastic
continuum21. For example, FE methods have been of Clementi et al.27 who used QC calculations to evalu-
ate the interaction of several water molecules; from this
broughttobearonproblemsofstrain-gradientplasticity,
such as geometrically necessary dislocations22. Contin- database, an empirical potential was parameterized for
use in molecular dynamics simulations; the MD simula-
uum Navier-Stokeequations arealso oftenused to study
tions were then used to evaluate the viscosity of water
fluids.
from atomic autocorrelation functions; and finally, the
The challenge in modern simulations of materials sci-
computed viscosity was employedin computational fluid
enceandengineeringisthatrealmaterialsusuallyexhibit
dynamics calculations to predict the tidal circulation in
phenomenaatonescalethatrequireaveryaccurateand
Buzzard’s Bay of Massachusetts.
computationally expensive description, and phenomena
The second category of multiscale simulations con-
at another scale for which a coarser description is satis-
sists of the so-called concurrent, or parallel, or explicit
factory and in fact necessaryto avoidprohibitively large
approaches. These approaches attempt to link meth-
computations. Since none of the methods above alone
ods appropriate at each scale together in a combined
would suffice to describe the entire system, the goal be-
model where the different scales of the system are con-
comestodevelopmodelsthatcombinedifferentmethods
sidered concurrently and communicate with some type
specialized atdifferent scales,effectively distributing the
of hand-shaking procedure. This approach is necessary
computational power where it is needed most. It is the
for systems that are inherently multiscale, that is, sys-
hope that a multiscale approach is the answer to such a
tems whose behavior at each scale depends strongly on
quest, and it is by definition an approach that takes ad-
what happens at the other scales. In contrastto sequen-
vantage of the multiple scales present in a material and
tial approaches, the concurrent simulations are still rel-
builds a unified description by linking the models at the
atively new and only a few models have been developed
differentscales. Fig. 1illustratestheconceptofaunified
to date. In a concurrent simulation, the system is often
multiscaleapproachwhichcanreachthelengthandtime
partitionedintodomainscharacterizedbydifferentscales
scale that individual methods, developed to treat a par-
and physics. The challenge of the concurrent approach
ticularscaleaccurately,failtoachieve. Atthesametime,
lies at the coupling between the different regions treated
the unifiedapproachcanretainthe accuracythatthein-
by different methods, and a successful multiscale model
dividualapproachesprovideintheirrespectivescales,al-
seeks a smooth coupling between these regions.
lowing, for instance, for very high accuracy in particular
regions of the systems as required. As effective theo- Inprinciple,multiscalesimulationscouldbebasedona
ries, multiscale models are also useful for gaining physi- hybrid scheme, using elements from both the sequential
cal insight that might not be apparent from brute force and the concurrent approaches. We will not examine
computations. Specifically, a multiscalemodel canbe an this type of approach in any detail, since it involves no
effective way to facilitate the reduction and the analysis new concepts other than the successful combination of
of data which sometimes can be overwhelming. Overall, elements underlying the other two types of approaches.
the goal of multiscale approaches is to predict the per- There already exist a few review papers and special
formance and behavior of materials across all relevant editions of articles on multiscale simulation of materials
length and time scales, striving to achieve a balance be- inliterature2,3,28,29,30,31,32. Amathematicperspectiveof
tween accuracy, efficiency and realistic description. multiscale modeling and computation is also available33.
3
The present overview does not aim to provide another text of the phase-field approach, where the lower scale
collection of various multiscale techniques, but rather to knowledge is in the form of ab initio free energies, and
identifythecharacteristicfeaturesandclassifymultiscale the coarse-grainedmodel is again a continuum model.
simulationapproachesintorationalcategoriesinrelation
to the problems where they apply. We select a few illus-
trative examples for each category and try to establish
A. Peierls-Nabarro model of dislocations
connectionsbetweentheseapproacheswheneverpossible.
Since almost all interesting material behavior and pro-
cessesaretimedependent, wewilladdressboththe issue Dislocations are central to our understanding of me-
of length-scales and the issue of time-scales integration chanical properties of crystalline solids. In particular,
in materials modeling. thecreationandmotionofdislocationsmediatetheplas-
Thepaperisorganizedasfollows: InSectionIIwedis- tic response of a crystal to external stress. While con-
cuss in detail representative examples of sequential mul- tinuum elasticity theory describes well the long-range
tiscale approaches in the spatial regime. In Section III elastic strain of a dislocation for length scales beyond
wepresentexamplesofconcurrentmultiscaleapproaches, a few lattice spacings, it breaks down in the immedi-
also in the spatial regime . In Section IV we discuss ate vicinity of the dislocation core. There has been a
representative approaches that extend time-scales in dy- great deal of interest in describing accurately the dis-
namical simulations. Section V contains our comments location core structure on an atomic scale because the
and conclusions for the applicability of the various ap- core structure to a large extent dictates the dislocation
proaches. The examples presented in this overview to properties34,35. So far, direct atomistic simulation of
some extent reflect our own research interests and they dislocation properties based on CIP has not been sat-
are by no means exhaustive. Nevertheless, we hope that isfactory because the CIP is not always reliable or may
they give a satisfactory cross-sectionof the current state even not be available for the material of interest, espe-
of the field, and they can serve as inspiration for further cially when the physical system involves several types of
developments in this exciting endeavor. atoms. Onthe other hand, ab initio calculationsarestill
computationally expensive for the study of dislocation
core properties, particularly of dislocation mobility. Re-
II. SEQUENTIAL MULTISCALE APPROACHES cently, a promising approach based on the framework of
the Peierls-Nabarro(P-N) model has attracted consider-
able interest for the study of dislocation core structure
Two ingredients are required in order to construct a
and mobility36,37,38,39,40,41,42,43,44,45,46,47. This approach
successful sequential multiscale model:
when combined with ab initio calculations for the ener-
(i) It is necessary to have a priori and complete knowl-
getics, represents a plausible alternative to the direct ab
edge of the fundamental processesat the lowestscale in-
initio simulations of dislocation properties.
volved. This knowledge or information can then be used
The P-Nmodelis aninherentlymultiscale framework,
for modeling the system at successively coarser scales.
first proposed by Peierls48 and Nabarro49 to incorporate
(ii) It is necessary to have a reliable strategy for en-
the details of a discrete dislocation core into an essen-
compassing the lower-scale information into the coarser
tially continuum framework. Consider a solid with an
scales. This is often accomplished by phenomenological
edge dislocation in the middle (Fig. 2): the solid con-
theories, which contain a few key parameters (these can
tainingthisdislocationisrepresentedbytwoelastichalf-
be functions), the value of which is determined from the
spacesjoined by atomic levelforcesacrosstheir common
information at the lower scale.
interface, known as the glide plane (dashed line). The
This message-passing approach can be performed in se-
goal of the P-N model is to determine the slip distribu-
quenceformultiplelengthscales,asintheexamplecited
intheintroduction27. Thekeyattributeofthesequential tionontheglideplane,whichminimizesthetotalenergy.
The dislocation is characterizedby the slip (relative dis-
approachisthatthesimulationatahigherlevelcritically
placement) distribution
depends on the completeness and the correctness of the
informationgatheredatthe lowerlevel,aswellasthe ef-
ficiency and reliability of the model at the coarser level. f(x)=u(x,0+)−u(x,0−) (1)
Toillustratethistypeofapproach,wewillpresenttwo
examplesofsequentialmultiscaleapproachesinsomede- which is a measure of the misfit across the glide plane;
tail. The first example concerns the modeling of dislo- u(x,0+) and u(x,0−) are the displacement of the half-
cation properties in the context of the Peierls-Nabarro spaces at position x immediately above and below the
(P-N)phenomenologicalmodel,wherethelowerscalein- glide plane. The total energy of the dislocated solid
formationisintheformoftheso-calledgeneralizedstack- includes two contributions: (1) the nonlinear potential
ingfaultenergysurface(alsoreferredtoastheγ-surface), energy due to the atomistic interaction across the glide
andthe coarse-grainedmodelisaphenomenologicalcon- plane, and (2) the elastic energy stored in the two half-
tinuum description. The second example concerns the spaces associated with the presence of the dislocation.
modeling of coherent phase transformations in the con- Bothenergiesarefunctionalsoftheslipdistributionf(x).
4
Specifically,thenonlinearmisfitenergycanbewrittenas To resolve these problems, a so-called Semidiscrete
VariationalP-N(SVPN)modelwasrecentlydeveloped43,
∞
U = γ(f(x))dx, (2) which allows for the first time the study of narrowdislo-
misfit
cations,asituationthatthestandardP-Nmodelcannot
Z−∞
handle. Within this approach, the equilibrium structure
where γ(f) is the generalized stacking fault energy sur-
ofadislocationisobtainedbyminimizingthedislocation
face (the γ-surface) introduced by Vitek50. The non-
energy functional
linear interplanar γ-surface can, in general, be deter-
mined from atomistic calculations. For systems where U =U +U +U +Kb2lnL, (5)
disl elastic misfit stress
CIP are not available or not reliable (for instance, it
is exceedingly difficult to derive reliable potentials for where
multi-component alloys), ab initio calculations can be
1
performed to obtain the γ-surface. On the other hand, U = χ [K (ρ(1)ρ(1)+ρ(2)ρ(2))+K ρ(3)ρ(3)],
elastic 2 ij e i j i j s i j
the elastic energy of the dislocation can be calculated i,j
X
reasonably from elasticity theory: the dislocation may (6)
be thought of as a continuous distribution of infinitesi-
mal dislocations whose Burgersvectors integrate to that U = ∆xγ (f ), (7)
misfit 3 i
oftheoriginaldislocation51. Therefore,theelasticenergy i
X
of the original dislocation is just the sum of the elastic
energy due to all the infinitesimal dislocations (from the
x2−x2
superpositionprinciple oflinearelasticitytheory),which U =− i i−1(ρ(l)τ(l)), (8)
stress 2 i i
can be written as
i,l
X
µ L df(x)df(x′)
Uelastic = 2π(1−ν) dx dx′ln|x−x′| dx dx′ with respect to the dislocationmisfit density. Here, ρ(i1),
Z Z (3) ρ(i2) andρ(i3) aretheedge,verticalandscrewcomponents
where µ, and ν are the shear modulus and Poisson’s ra- ofthegeneralinterplanarmisfitdensityatthei-thnodal
tio, respectively. L is an inconsequential constant intro- point, and γ3(fi) is the corresponding three-dimensional
duced as a large-distance cutoff for the computation of γ-surface. Thecomponentsoftheappliedstressinteract-
thelogarithmicinteractionenergy52. NotethattheBurg- ingwiththeρ(i1),ρ(i2) andρ(i3),areτ(1) =σ21,τ(2) =σ22
ers vector of each infinitesimal dislocation is the local and τ(3) =σ , respectively. K, K and K are the pre-
23 e s
gradient in the slip distribution, df(x)/dx. The gradient logarithmicelasticenergyfactors52. Thedislocationden-
of f(x) is called dislocation (misfit) density, denoted by sityatthei-thnodalpointisρ =(f −f )/(x −x )
i i i−1 i i−1
ρ(x). Since the P-N model requires that atomistic infor- and χ is the elastic energy kernel43.
ij
mation(embodied in the γ-surface)be incorporatedinto The firstterm in the energyfunctional, U is now
elastic
a coarse-grainedcontinuum framework,it is a sequential discretized in order to be consistent with the discretized
multiscale strategy. Thus the successful application of misfit energy, which makes the total energy functional
the method depends on the reliability of both γ-surface variational. Another modification in this approach is
andtheunderlyingelasticitytheorywhichisthebasisfor that the nonlinear misfit potential in the energy func-
the formulation of the phenomenological theory. tional, U , is a function of all three components of
misfit
In the current formulation, the total energy is a func- the nodal misfit, f(x ). Namely, in addition to the misfit
i
tionalofmisfitdistributionf(x)orequivalentlyρ(x),and alongtheBurgersvector,lateralandevenverticalmisfits
itisinvariantwithrespecttoarbitrarytranslationofρ(x) across the glide plane are also included. This allows the
and f(x). In order to regain the lattice discreteness, the treatmentofstraightdislocationsofarbitraryorientation
integrationoftheγ-energyinEq. (2)wasdiscretizedand in arbitrary glide planes. Furthermore, because the mis-
replacedbyalatticesumintheoriginalP-Nformulation, fit vector f(x ) is allowed to change during the process
i
of dislocation translation,the energy barrier (referredto
∞
as the Peierls barrier) can be significantly lowered com-
U = γ(f )∆x, (4)
misfit i pared to its corresponding value from a rigid transla-
i=−∞
X tion. Theresponseofadislocationtoanappliedstressis
with x the reference position and ∆x the average spac- achieved by minimization of the energy functional with
i
ing of the atomic rows in the lattice. This procedure, respecttoρ atthegivenvalueoftheappliedstress,τ(l).
i i
however,is inconsistentwith evaluation ofelastic energy An instability is reachedwhen anoptimalsolutionfor ρ
i
[Eq. (3)] as a continuous integral. Therefore the total no longer exists, which is manifested numerically by the
energy is not variational. Furthermore in the original P- failure of the minimization procedure to converge. The
Nmodel, the shape ofthe solutionf(x) is assumedto be Peierlsstressisdefinedasthecriticalvalueoftheapplied
invariant during dislocation translation, a problem that stress which gives rise to this instability.
isalsoassociatedwiththenon-variationalformulationof The SVPN model has been applied to study vari-
the total energy. ous interesting material problems related to dislocation
5
phenomena45,46,47. One study involved the understand- continuum theories1. The same strategy can also be ap-
ing of hydrogen-enhanced local plasticity (HELP) in Al. plied to the study of fracture and dislocation nucleation
HELP is regarded as one of three general mechanisms from a crack tip55.
responsible for H embrittlement of metals53. There has It is interesting to note that the analysis of γ-surface
been overwhelming experimental evidence in support of can provide a qualitative understanding of even more
HELP,butatheoreticalfoundationwaslacking. Inorder complex mechanical properties of materials. For exam-
to gain understanding of the physics behind the HELP ple, Rice and coworkers59 formulated powerful criteria
mechanism, Lu et al. carried out ab initio calculations for the brittle behavior of materials, by extending the
for the γ-surface of Al with H impurities placed at the Peierlsanalysisto geometriesinvolvingcracks. Basedon
interstitial sites45. The γ-surface for both pure Al and this framework, Waghmare et al.56,57 were able to pre-
the Al+H systems is shown in Fig. 3. Comparing the dict which alloying elements could improve the ductility
twoγ-surfaces,onefindsanoverallreductioninγ energy ofMoSi byanalyzingtheab initiodeterminedγ-surface
2
inthepresenceofH,whichisattributedtothechangeof of the alloys, and comparing the changes induced by al-
atomicbondingacrossthe glideplane,fromcovalent-like loyingtokeyfeaturesoftheγ-surfaceversusthe changes
to ionic-like54. induced to the surface energy γ . Remarkably, certain
s
Thecorepropertiesoffourdifferentdislocations,screw predictions of this relatively simple theoretical modeling
(0◦), 30◦, 60◦ and edge (90◦) have been studied using were borne out by subsequent experiments58.
the SVPN model combined with the ab initio deter- We have devoted some attention to the description of
mined γ-surface. It was found that the Peierls stress the P-N model and its implementation using ab initio
for these dislocations is reduced by more than an order γ-surfaces, because it is an ideal case of a sequential
ofmagnitudeinthe presenceofH45,whichiscompatible multiscale model: it consists of a well motivated phe-
with the experimental findings that support the HELP nomenologicalframework,within which the set of atom-
mechanism53. Moreover, in order to address the exper- istically derived quantities is well defined and complete
imental observation for H trapping at dislocation cores (in this case the γ-surface). In this sense, it fulfills all
and H-induced slip planarity, the H binding energy to the requirements for a coherent and complete multiscale
the dislocation cores was calculated45. These calcula- model. Therearenodoubtlimitationstoit,arisingfrom
tions showedthatthere is strongbinding betweenH and therangeofvalidityofthe phenomenologicaltheory,but
dislocationcores,thatis,Hisattracted(trapped)todis- within this range there are no other ambiguities in con-
location cores which lowers the core energies. More im- structing the multiscale model. Perhaps, its successes,
portantly,the binding energywasfoundto be a function someofwhichwepresentedabove,areowedtothiscom-
ofdislocationcharacter,withtheedgedislocationhaving plete character of the model.
the greatest and the screw dislocation having the lowest
binding energy. For a mixed dislocation, the binding en-
ergyincreaseswiththeamountofedgecomponentofthe B. Phase-field model of coherent phase
Burgers vector. These results suggest that in the pres- transformations
ence of H, it costs more energy for an edge dislocation
to transform into a screw dislocation in order to cross- The structure-properties paradigm is one of the prin-
slip,sincetheedgedislocationhasalmosttwicethebind- cipal pillars in materials science. The term “structure”
ing energy of the screw dislocation45. In the same vein, here refers to structures at many different scales,includ-
it costs more energy for a mixed dislocation to transfer ing the atomic scale geometry determined by the crys-
its edge component to a screw component for cross-slip. talline arrangement of atoms, the structure of the de-
Therefore the cross-slip process is suppressed due to the fects that exist in a material, and the structure that
presence of H, and slip is confined to the primary glide emerges as a result of the organization of these defects
plane, exhibiting the experimentally observed slip pla- into what is referred to as microstructure. Among these
narity. structures, the microstructure on the scale of microme-
A similar approach was applied to the study of the ters is often directly tied to the mechanical properties
vacancy lubrication effect on dislocation motion in Al. of materials, and has therefore attracted great interest
From this analysis, it was shown that the role of vacan- both in terms of scientific understanding and practical
cies is crucial in reconciling the results of Peierls stress applications60,61,62,63,64.
measuredfromdifferentexperimentaltechniques46. Very Recently, a powerful sequential multiscale approach
recently,amultiple planeP-Nmodelhasbeendeveloped has been put forward for modeling the precipitate
tostudydislocationphenomenainvolvingmorethanone microstructure and its evolution in multicomponent
glide planes, such as dislocation constriction and cross- alloys65,66,materialswhichappearinmanytechnological
slip47. FinallyweshouldpointoutthattheP-Nmodelis applications. The approach is based on the continuum
justoneexampleofmoregeneralcohesivesurfacemodels phase-fieldmodelwhosedrivingforces(freeenergies)are
that are built upon the idea of limiting all constitutive obtained from combined ab initio calculations and the
nonlinearitytocertainprivilegedinterfaces,whilethere- mixed-space cluster expansion technique. One interest-
mainder of the material is treatedvia more conventional ing application of this approach concerned the study of
6
precipitationoftheθ′(Al Cu)phaseinCu-Alalloysdur- for accurately determining these quantities is of critical
2
ing thermal aging66. importancetopredictthemicrostructureevolutionofin-
In the phase-field multiscale approach, the nature of terest. In particular, if the quantities can be determined
phase transformationas wellas the microstructuresthat from ab initio calculations, the goal of an “ab initio”
are produced are described by a set of continuous order- modeling of alloy microstructure evolution would be, to
parameter fields67,68. The temporal microstructure evo- a great extent, achieved73,74.
lution is obtained from solving field kinetic equations
Since direct ab initio calculations of free energies are
that govern the time-dependence of the spatially inho-
either impractical or impossible with currently available
mogeneous order-parameter fields. Within the diffuse-
computational power, a useful method has been devel-
interface description, the thermodynamics of a phase
opedtoextendtheabinitioenergeticstothermodynamic
transformation and the accompanying microstructure
properties of alloy systems with hundreds of thousands
evolution are modeled by a free energy that is a func- of atoms75, referred to as the mixed-space cluster ex-
tion of the order-parameter field, or phase field. For a
pansion (CE). In this scheme, energetics from ab initio
structural transformation, the total free energy can be
calculationsfor anumber ofsmallunitcell(∼ 10atoms)
written as:
structuresaremappedontoageneralizedIsing-likemodel
for a particular lattice type, involving 2-, 3-, and 4-
F =F +F +F (9)
tot bulk inter elast body interactions plus coherency strain energies76. The
Hamiltoniancanbeincorporatedintomixed-spaceMonte
where F is the bulk free energy, F is the interfa-
bulk inter Carlo simulations of N ∼ 105 atoms, effectively allow-
cialfreeenergy,andF is the coherencyelasticstrain
elast ing one to explore the complexity of 2N configurational
energyarisingfromthe lattice accommodationalongthe
space. As demonstrated by Vaithyanathan et al., the
coherentinterfacesinamicrostructure. Foramicrostruc-
bulk free energy can be obtained from Monte Carlo sim-
turedescribedbyacompositionfieldcandasetofstruc-
ulations coupled with thermodynamic integration tech-
tural order-parameters, η , ..., η , the first two terms of
1 p
niques. The precipitate/matrix interfacial free energies
Eq. (9) are given by
canbe determinedfromsimilarMonte Carlosimulations
α orfromlowtemperatureexpansiontechniques. Theelas-
F +F = {f[c(r),η (r)]+ |∇c(r)|2 (10)
bulk inter p tic strain energies are of precisely the same form as the
2
ZV coherencystrainenergyusedtogeneratethemixed-space
1
+ β (p)∇ η (r)∇ η (r)}dV CE.Hence,fromacombinationofabinitiocalculations,a
ij i p j p
2
p mixed-spaceCEapproach,andMonteCarlosimulations,
X
one can obtain all the driving forces needed as input to
where f(c,ηp) is the local free energy density69 and α the continuum phase-field model. The incorporation of
andβij(p)arethegradientenergycoefficientswhichcon- theseenergeticproperties,obtainedfromatomistics,into
trol the width of the diffuse interface. The elastic strain a continuum microstructural model represents a bridge
energyisobtainedfromelasticitytheoryusingthehomo- between these two length scales, and opens the path to-
geneous modulus approximation70. With the total free ward predictive modeling of microstructures and their
energy of an inhomogeneous system written as a func- evolution.
tion of order-parameter fields, the temporal evolution of
Toillustratetheuseofthemethod,wementionbriefly
microstructuresduringaphasetransformationcanbeob-
the work of Vaithyanathan et al. who studied the prob-
tainedby solvingthe coupledCahn-Hilliardequationfor
lem of precipitation of the θ′ (Al Cu) phase in Cu-Al al-
a conserved field c, and the time-dependent Ginzburg- 2
loys. The free energy of the θ′ phase is obtained fromab
Landau equation for a non-conserved field η 71,72:
p initiocalculationsoftheenergyatT =0K,coupledwith
thecalculatedvibrationalentropyofthisphase. Thebulk
∂c ∂f
=M∇2 −α∇2c (11) freeenergiesofmatrixandprecipitatephasesarethenfit
∂t ∂c
(cid:20) (cid:21) to the local free energy as a function of order-parameter
fields in the phase-field model. T =0 K interfacial ener-
∂η δF gies are determined from supercell calculations, both for
p tot
=−L (12)
∂t p δη the coherent interface and for the incoherent interface.
p
The anisotropy of these interfacial energies is large and
whereM is relatedtoatommobilityandL is the relax- has been incorporated in the phase-field model. Elastic
p
ation constant associated with the order-parameter η . energy calculations for the coherent strain of Al/Al Cu
p 2
As the above equations illustrate, the continuum phase- (θ′) and the calculated lattice parameters of each phase
field methodology depends on three input energies: (1) determine the elastic driving force in this system. Hav-
bulkfreeenergiesofsolidsolutionandprecipitatephases, ing determined all the necessary thermodynamic input,
(2) precipitate-matrix interfacial free energies, and (3) Vaithyanathan et al. were able, for the first time, to
precipitate/matrix lattice elastic energies. Experimen- clarify the physical contributions responsible for the ob-
tal determination of these quantities can be difficult and served morphologyof θ′ precipitate microstructure. The
problematic. Therefore a physically motivated method agreementbetweenthecalculatedandexperimentallyob-
7
served microstructure of θ′ in the Al-Cu alloys was ex- equation for ϕ is
cellent, confirming the validity of the approach.
∂ϕ
Althoughthephase-fieldmodelisableto predictcom- +v·∇ϕ=0 (13)
∂t
plex microstructure evolution during phase transforma-
tions, it requires as input phenomenological thermo- Growthisnaturallydescribedbythesmoothevolutionof
dynamic and kinetic parameters. For binary systems, ϕdeterminedbythis differentialequation. Inthe caseof
ab initio calculations can provide these parameters for multilayergrowth,theboundariesΓ (t)oftheislandsare
k
the phase-field model, but it is unrealistic to assume definedasthesetofspatialpointsxforwhichϕ(x,t)=k
that such calculations can be used to determine all the fork =0,1,2,.... The evolutionofthe level-setfunction
thermodynamic informationfor systems beyond ternary. ϕ can be obtained by numerically solving Eq. (13) using
Therefore semi-empirical methods, such as CALPHAD non-oscillatory methods86. The key parameters enter-
(calculated phase-diagram) will remain a useful tool in ing the model are diffusion constants (the terrace and
such an endeavor77,78,79. island-edge diffusion constants) which can in principle
be supplied from atomistic calculations, through the fol-
lowing procedure (see Fig. 4): first, the atomistic pro-
cesses are identified which are responsible for terrace or
C. Other sequential approaches island-edgediffusionandtheirenergeticsareanalyzedus-
ing atomistic (possibly ab initio) calculations; next, the
energy barriers for the atomistic processes are incorpo-
KineticMonteCarlo(KMC)simulations,coupledwith
rated in a KMC model which provides the means for
atomistically determined kinetic energy barriers, repre-
coarse-grainingtheatomisticdegreesoffreedomtoafew
sentapowerfulclassofsequentialmultiscaleapproaches.
mesoscopic degrees of freedom describing the evolution
For example, a large body of research has been car-
ofsurfacefeatures (the islandstepedges); finally,the re-
ried out for surface growth phenomena with KMC sim-
sults of the KMC model are coarse-grained to provide
ulations whose kinetic energy barrier parameters for
the input to the level-set equations, that is, they define
relevant elemental processes are supplied by ab initio the values of the boundary velocity v which depends on
calculations80,81. In an altogether different field, Cai
the local surface morphology. The coarse-graining be-
et al. have used KMC method to study dislocation
tween scales eliminates degrees of freedom that are not
motion in Si based on the well-established double-kink
essential, making the passage to the next scale feasible.
mechanism82. In their approach, the dislocation is rep-
For example, in the illustration shown in Fig. 4, the
resented by a connected set of straight line segments
smallest step width in the KMC scale corresponds to a
which move as the cumulative effect of a large number
two-atom wide region at the microscopic scale, a situ-
of kink nucleation and migration processes. The rate of
ation that is relevant to the Si(100) surface and possi-
these processesis calculatedfromtransitionstate theory
bly to other semiconductor surfaces (such as III-V com-
with the transition energy barrier having contributions
poundsurfaces). Inthesecases,surfaceatomstendtobe
from both atomistically determined energetics (double-
boundtodimerpairs,whichistheessentialunitthatde-
kink formation and migration energy) and elastic inter-
termines the step structure, even though the underlying
actions with other dislocation segments as well as from
dynamics may be determined by the motion of individ-
the externally applied stress.
ual atoms. Thus, the KMC simulation need only take
An example of a multiscale approach, in which KMC intoaccountstructuresconsistingofdimerunits,the dy-
is a key component, employs the so-called level-set namics of which determine the step-edge motion needed
method83,84 forthe largest(macroscopic)scale. Thisap- for the level-set simulation. The middle terrace in Fig.
proachisparticularlywellsuitedforthestudyofepitaxial 4(b) is shown as a grid of squares, each representing a
growth,asubjectofgreatimportanceinmicroelectronics four-atomclusterandbeing the minimalunit relevantto
andoptoelectronicsapplications. Inthelevel-setmethod, stepmotionatthe KMC scaleinthis example,assuming
growthisdescribedbycreationandsubsequentmotionof that only steps of width equal to two atoms in each di-
island boundaries. The model treats the growing film as rection are stable. The level-set method is a manifestly
acontinuuminthelateraldirection,butretainsatomistic multiscale approach, combining information from three
discretenessin the growthdirection. Inthe lateraldirec- different regimes (atomistic, mesoscopic and continuum)
tion,continuumequationsrepresentingthefieldvariables into a neatly integrated scheme. Recently, the level-set
can be coupled to growth through island evolution, by method has also been applied to study dislocation dy-
solving the appropriate boundary-value problem for the namics in alloys87.
field and using local values of this field to determine the Yet another sequential multiscale approach has been
velocity of the island boundaries. The central idea be- successfully applied to the study of crystal plasticity.
hind the level-set method85 is that any boundary curve This is the DD method mentioned earlier, incorporating
Γ, such as a step or the boundary of an island, can be dislocation motion at the macroscopic scale, the mech-
representedas the set of values ϕ=0 (the level-set)of a anism ultimately responsible for crystal plasticity. In
smoothfunction ϕ. For a givenboundary velocityv, the order to predict the mechanical properties of materials
8
using DD simulations, a connection between micro-to- To address these challenges, Abraham, Broughton,
meso scales must be established because dislocation in- Bernstein and Kaxiras developed a concurrent multi-
teractions at close range (when the cores intersect, for scale modeling approach that dynamically couples dif-
instance), are totally beyond the reach of continuum ferent length scales25,26. This multiscale methodology
models. Along these lines, Bulatov et al. were able to aimsatlinkingthelengthscalesrangingfromtheatomic
study dislocation reactions and plasticity in fcc metals88 scale, treated with a quantum-mechanical tight-binding
thatcomparewellwithdeformationexperiments,byinte- approximation method, through the microscale, treated
grating the local rules derived from atomistic simulation via the classical molecular dynamics method, and fi-
of dislocation core interactions into the DD simulations. nally to the mesoscale/macroscale treated via the finite
The same idea has been further explored by Rhee et al. element method in the context of continuum elastic-
in a study of the stage I stress-strain behavior of bcc ity. These authors applied this unified approach,termed
single crystals19. macroscopic-atomistic-ab initio dynamics (MAAD), to
the study of the dynamical fracture process in Si, a typ-
ical brittle material. In traditional studies of fracture,
III. CONCURRENT MULTISCALE only the continuum mechanics level (employing, for in-
APPROACHES stance,the FEmethod) isusuallyinvokedtoaccountfor
the macroscopic behavior. But since there is no natu-
Broadly speaking, a concurrentmultiscale approachis ralsmall-lengthcutoffpresentinthecontinuummechan-
moregeneralinscopethanitssequentialcounterpartbe- ics approach, any important aspect of the atomic-scale
cause the concurrent approach does not rely on any as- mechanisms for fracture is completely missed. This can
sumptions (in the form of a particular coarse-graining be remedied by introducing classical MD to the simula-
model) pertaining to a particular physical problem. As tions. In particular, the MAAD approach employed the
a consequence, a successful concurrent approach can be Stillinger-Weber90 interatomic empirical potential for Si
usedtostudymanydifferentproblems. Forexample,dis- to perform MD calculations at the atomistic level, for a
location core properties, grain boundary structure and largeregionofthematerialnearacracktip. However,the
crack propagation could all be modeled individually or treatment of formation and breaking of covalent bonds
collectively by the same concurrent approach, as long atthe atomic scale is not reliable with any empirical po-
as it incorporates all the relevant features at each level. tential, because bonds between atoms are an essentially
Whatisprobablymostappealing,however,isthatacon- quantum mechanical phenomenon arising from the shar-
current approach does not require a priori knowledge of ing of valence electrons. On the other hand, small devi-
the physical quantities or processes of interest. Thus, ationsfromidealbonding arrangementscanbe captured
concurrent approaches are particularly useful to explore accurately by empirical potentials, because they are to
problems about which little is known at the atomistic first approximationharmonic, a feature that is easily in-
level and its connection to larger scales, and to discover corporated in empirical descriptions of the interaction
new phenomena. We discuss below three instances of between atoms. Therefore, it was deemed necessary to
concurrentapproachesinsomedetail, andmentionsome include aquantummechanicalapproachinto the simula-
additional examples more briefly. tions for a small region in the immediate neighborhood
of the crack tip, where bond breaking is prevalent dur-
ing fracture,while further awayfromthis regionthe em-
pirical potential description is adequate. The particular
A. Macroscopic Atomistic Ab initio Dynamics
methodology chosen to model the immediate neighbor-
hood of the crack-tip, a semi-empirical nonorthogonal
Fracturedynamicsisoneofthemostchallengingprob-
tight-binding scheme91, describes well the bulk, amor-
lems in materials science and solid mechanics. Despite
phous, and surfaces properties of Si. Fig. 5 shows the
nearlyacenturyofstudy,severalimportantissuesremain
spatial decomposition of the computational cell into five
unsolved. Inparticular,thereislittlefundamentalunder-
different dynamic regions of the simulation: the contin-
standing of the brittle to ductile transition as a function
uum FE region at the far-field where the atomic dis-
oftemperatureinmaterials;thereisstillnodefinitiveex-
placements and strain gradients are small; the atomistic
planation of how fracture stress is transmitted through
MD region around the crack with large strain gradients
plastic zones at crack tips; and there is no complete un-
but with no bond breaking; the quantum mechanical re-
derstanding of the disagreement between theory and ex-
gion(labelled TB because ofthe use ofthe tight-binding
periment regarding the limiting speed of crack propaga-
method) right at the crack tip where atomic bonds are
tion. These difficulties stem from the fact that fracture
being broken and formed; the FE/MD hand-shaking re-
phenomena are governed by processes occurring over a
gion; and the MD/TB hand-shaking region. The total
widerangeoflengthscalesthatareallconnected,andall
Hamiltonian, H for the entire system was written as:
contribute to the total fracture energy89. In particular, tot
the physics on different length scales interacts dynami-
H = H ({u,u˙}∈FE)+
cally, therefore a sequential coupling scheme would not tot FE
be adequate for the study of fracture dynamics. H ({r,r˙}∈MD)+
MD
9
H ({r,r˙}∈TB)+ words, the TB terminating atoms are fictitious mono-
TB
H ({u,u˙,r,r˙}∈FE/MD)+ valent atoms forming covalent bonds with the strength
FE/MD
andlength ofbulk Si bonds. These fictitious atoms were
H ({r,r˙}∈MD/TB)
MD/TB called “silogens”: they behave mechanically just like Si,
but chemically like H. The TB Hamiltonian including
ThedegreesoffreedomoftheHamiltonianareatomicpo-
silicon-siliconand silicon-silogenmatrix elements is then
sitionsrandvelocitiesr˙ fortheTBandMDregions,and
diagonalized to obtain electronic energies and wavefunc-
displacements u and their time rates of change u˙ for the
tions, from which the total energy can be computed.
FEregions. Equationsofmotionforalltherelevantvari-
Thus, at the perimeter of the MD/TB region, there are
ables in the system are obtained by taking appropriate
silogens sitting directly on top of the atoms of the MD
derivatives of this Hamiltonian. All variables can then
region,whichare shownas the smallerred circles ontop
be updated in lock-step as a function of time using the
of green circles in Fig. 5. On one side of the TB/MD
same integrator. Thus the entire time history of the sys-
interface,the bonds to an atomarederivedfromthe TB
tem may be obtained numerically given an appropriate
Hamiltonian, and are shown as shaded regions in Fig. 5,
set of initial conditions. Following trajectories dictated
to indicate the electronic distribution responsible for the
bythisHamiltonianleadstoevolutionofthesystemwith
formation of the covalent bonds. On the other side of
conservedtotalenergy,whichensuresnumericalstability.
theinterface,the bondsarederivedfromtheinteratomic
The individual approaches at each level (FE, MD and
potential of the MD simulation. The MD atoms of the
TB) are well established and tested methods. What was
interface have a full complement of neighbors, including
much more important in this study was the seamless
neighborswhosepositionsaredeterminedbythedynam-
hand-shaking of the different methods at the interface
ics of atoms in the TB region; these are shown as small
of the respective domains, namely the hand-shaking al-
greencirclesontopoftheredcirclesinFig. 5. Asbefore,
gorithms between FE and MD regions and between the
the TB code updates atomic positions in lock-step with
MD and TB regions. We present here the main ideas
its FE and MD counterparts.
behind the coupling of the different regions.
FE/MD coupling: To achievethe FE/MD hand-shaking, The MAAD approachwasemployedto study the brit-
the FE mesh spacing is scaled down to atomic dimen- tle fracture of Si in a geometry containing a small crack
sions at the interface of the two regions. In Fig. 5, (notch) within an otherwise perfect solid, with the ex-
the FE nodes are indicated as small open circles con- posednotchfaceinthe(100)planeandthenotchpointed
nected by thin lines. Moving away from the FE/MD in the h010i direction. The system consisted of 258,048
regionanddeep into the continuum, one canexpand the mesh points in each FE region, 1,032,192 atoms in the
mesh size. In this way, the atomistic simulation is em- MD region, and approximately 280 unique atoms in the
bedded in a large continuum solid, indicated by a green- TB region (for computational reasons, the entire region
colored region in Fig. 5(a). FE cells contributing fully modeledbytheTBmethodwasbrokenintosmaller,par-
totheoverallHamiltonian(unitweight)aremarkedwith tially overlapping regions, each assigned to a different
thinsolidlines,whilecellscontributingtothehand-shake processor in a parallel implementation). The lengths of
Hamiltonian(halfweight)arerepresentedbythindashed the MD region are 10.9 ˚A for the slab thickness along
lines. Interactions between the atoms on the MD side, thefrontofthecrack,3649˚A intheprimarydirectionof
which are representedby an interatomic potential, carry propagation, and 521 ˚A in the direction of pull (before
full weight when fully inside the MD region (thick solid pulling). Periodic boundary conditions were imposed at
lines joining neighboring atoms) and half weight (thick the slab faces normal to the direction of the crack prop-
dashed lines) when they cross the boundary, with one of agation (along the front of the crack). The wall-clock
theneighborseffectivelyrepresentedbyanodeintheFE time for a TB force update was 1.5 s, that for the MD
region. The FE/MD interface is chosen to be far from update was 1.8 s, and that for the FE update was 0.7 s.
the fracture region. Hence, the atoms of the MD region The TB region was relocated after every 10 time steps
and the displacements of the FE lattice can be unam- toensurethatitremainsatthe verytipofthepropagat-
biguously mapped onto one another. The assignment of ingcrack. Thecomputationalslabwasinitializedatzero
weights equal to unity within each region and equal to temperature, and a constant strain rate was imposed on
one half at the interface is arbitrary and can be general- the outermostFE boundaries defining the opposing hor-
ized by the introduction of a smooth step function. izontal faces of the slab. Furthermore, a linear velocity
MD/TB coupling: At this interface, the atoms treated gradientwasapplied within the slab, which results in an
quantum mechanically are shown in red while those increasinginternalstrainwithtime. Itwasobservedthat
treated classically are shown in green. The dangling theSisolidfailedinbrittlefashionatthenotchtipwhen
bonds at the edge of the TB region are terminated with the material is stretched by ∼ 1.5%. The limiting speed
pseudo-hydrogen atoms. The Hamiltonian matrix ele- ofcrackpropagationwasfoundtobe85%oftheRayleigh
ments ofthese pseudo-hydrogenatoms arecarefullycon- speedwiththechosencomputationalsetup. Inthecourse
structed to tie off a single Si bond and to ensure the of the simulation, the straight-ahead brittle cleavage of
absence of any charge transfer when that atom is placed the Si slab left behind a rough surface, with increasing
in a position commensurate with the Si lattice. In other roughening as a function of crack distance. Based on
10
these results, the authors suggested that the roughening B. Quasicontinuum model
surface is due to the spawning of dislocations with low
mobility on the time-scale of the crack motion.
One observation from many large-scale atomistic sim-
Ageneralproblemassociatedwithdomaindecomposi- ulations is that only a small subset of atomic degrees
tion, as in the MAAD simulations,is the spuriousreflec- of freedom do anything interesting. The great major-
tion of elastic waves (phonons) at the domain bound- ity of the atoms behave in a way that could be de-
aries due to the changes in system description across scribedby effective continuummodels like elasticity the-
the boundaries. For example, such effects have been ob- ory. The computation and storage of the uninteresting
servedintheatomisticmodelingofdislocationmotion92, degrees of freedom - necessary for a fully atomistic cal-
crack propagation93,94,95,96, and energetic particle-solid culation - consume a large proportion of computational
collisions97,98, all of which involved some domain cou- resources. This observationcalls for novelmultiscale ap-
pling scheme. Since the MAAD method involves do- proacheswhichcanreducethenumberofdegreesoffree-
main decomposition into the TB, MD and FE regions, dom in atomic simulations101,102. One such approach
thequalityofcouplingbetweendifferentregionsneedsto proposed by Tadmor, Ortiz and Phillips is particularly
be examined. In a subsequent paper, the same authors promising and has yielded considerable success in many
reported that there was no visible reflection of phonons applications103. This concurrent multiscale approach is
at the FE/MD interface, and no obvious discontinuities calledthequasicontinuummethod,whichseamlesslycou-
attheMD/TBinterface26. Thus,inthisschemethecou- ples the atomistic and continuum realms. The chief ob-
plingbetweenthevariousdomainsisindeedperformedin jective of the approach is to systematically coarsen the
a seamless manner, closely mimicking the actual behav- atomistic description by the judicious introduction of
ior of the physical system under investigation. Overall, kinematic constraints. These kinematic constraints are
the MAAD approach represents the state of the art of selectedanddesignedsoastopreservefullatomisticres-
current multiscale simulation strategies. It is a finite- olution where required - for example, in the vicinity of
temperature, dynamic and parallel algorithm which, at latticedefects-andtotreatcollectivelylargenumbersof
least as far as general computational aspects are con- atomsinregionswherethedeformationfieldvariesslowly
cerned, is applicable to any type of material. on the scale of the lattice. Variants of the quasicontin-
Ongoing efforts are exploring the possibility of apply- uum model havebeen developedandapplied in different
ing the MAAD strategyto study chemicaleffects onme- situations103,104,105,106,107,108,109,110,111,112,113,114. The
chanical properties of metallic alloys, such as impurity essential building blocks of the static quasicontinuum
effects ondislocationmotion,cracknucleationandprop- modelare: (1)theconstrainedminimizationoftheatom-
agation in various metals. There is an important qual- isticenergyofthesolid;(2)theuseofsummationrulesto
itative difference between such systems and the study compute the effective equilibrium equations; and (3) the
of brittle fracture of Si mentioned above: the nature of use ofadaptationcriteriainorderto tailorthe computa-
bonds in metallic systems is very different from the sim- tionalmeshtothe structureofthedeformationfield. An
plecovalentbondsinSi. Thismakesnecessarythedevel- extension of the method to finite-temperature has also
opment of a different way of coupling the quantum me- been proposed115.
chanicaltotheclassicalatomisticregion,becauseitisno The quasicontinuum model182 starts from a conven-
longerfeasibletoterminatethebondsattheboundaryof tional atomistic description, which computes the energy
the quantum region by simply saturating them with fic- of the solid as a function of the atomic positions. The
titious atoms like the silogens. In such endeavors, other configurationspace ofthe solid is then reduced to a sub-
more efficient and versatile quantum mechanical formu- setofrepresentativeatoms. Thepositionsoftheremain-
lations are desirable. One candidate is the linear scal- ingatomsareobtainedbypiecewiselinearinterpolations
ing real-space kinetic energy functional method99. This oftherepresentativeatomcoordinates,muchinthesame
method approximates the non-interacting kinetic energy manner as displacement fields are constructed in the FE
of DFT as a functional of electron density, and elec- method. The effective equilibrium equations are then
tronic wave-functions are thus eliminated from calcula- obtained by minimizing the potential energy of the solid
tions, and therefore the method is termed as orbital-free over the reduced configuration space. A direct calcula-
density functional theory (OFDFT). As a consequence, tion of the total energy in principle requires the evalua-
no diagonalization of the electronic Hamiltonian and no tion of sums that are extended over the full collection of
sampling of reciprocal space are necessary, making the atoms, namely,
method computationally efficient100. In particular, the
explicit real-spacefeature of this approachmakes it nat- N
urally suitable for domain coupling within the MAAD E = E , (14)
tot i
framework. Whileeffortstoconstructafullyfunctioning i=1
X
scheme along these lines are continuing, we believe this
is a promising method with great potential for applica- where N is the total number of atoms in the solid. The
tions in metallic systems, which are difficult to handle full sums may be avoidedby the introduction of approx-
with other techniques. imatesummationrules. Forexample,the latticequadra-
