Table Of ContentJanuary31,2012 2:3 PhilosophicalMagazine prim
Philosophical Magazine
Vol. 89, No. 34-36, 1-21 December 2009, 3371–3391
Composition-dependent interatomic potentials:
A systematic approach to modelling multicomponent alloys
2
1
B. Sadigha, P. Erharta, A. Stukowskib, and A. Caroa
0
2 a Condensed Matter and Materials Division, Lawrence Livermore
National Laboratory, Livermore, CA; b Institut fu¨r Materialwissenschaft,
n
a Technische Universita¨t Darmstadt, Germany
J (Received 15 July 2009; final version received 27 August 2009)
0
3
We propose a simple scheme to construct composition-dependent interatomic potentials for
multicomponentsystemsthatwhensuperposedontothepotentialsforthepureelementscan
]
i reproducenotonlytheheatofmixingofthesolidsolutionintheentireconcentrationrangebut
c alsotheenergeticsofawiderrangeofconfigurationsincludingintermetallicphases.Weshow
s thatanexpansioninclusterinteractionsprovidesawaytosystematicallyincreasetheaccuracy
- of the model, and that it is straightforward to generalise this procedure to multicomponent
l
r systems.Concentration-dependent interatomicpotentials canbebuiltuponalmostanytype
t ofpotentialforthepureelementsincludingembeddedatommethod(EAM),modifiedEAM,
m
bond-order,andStillinger-Webertypepotentials.Ingeneral,composition-dependentN-body
. termsinthetotal energy leadtoexplicitN+1-bodyforces,whichpotentially rendersthem
t computationally expensive. We present an algorithm that overcomes this problem and that
a
canspeed upthe calculation of the forces forcomposition-dependent pairpotentials insuch
m
a way as to make them computationally comparable in efficiency and scaling behaviour to
- standard EAM potentials. We also discuss the implementation in Monte-Carlo simulations.
d Finally,weexemplarilyreviewthecomposition-dependent EAMmodelfortheFe–Crsystem
n [PRL95,075702,(2005)].
o
c Keywords:empiricalpotentials; multicomponentalloys;concentrated alloys;computer
[ simulations;moleculardynamics;MonteCarlo;composi-tiondependent interatomic
potentials;clusterinteractions
1
v
4
8
2 1. Introduction
6
.
1 Twenty-five years ago, the Finnis-Sinclair many body potential [1], the Embedded
0
Atom Model of Daw and Baskes [2], the Glue model of Ercolessi and Parrinello [3],
2
andtheeffective mediumtheoryduetoPuska,NieminenandNorskov [4,5]marked
1
: the birthday of modern atomic scale computational materials science, enabling
v
i computersimulationsatthemultimillionatomscaletobecomearoutineinmodern
X
materials science research. This family of many body potentials share in common
r thefactthattheexpressionforthetotal energyhas nonlinear contributions ofpair
a
functions, removing in this way the limitations of the pair potential formulation to
describe realistic elastic properties.
Alloys and compounds, where the thermodynamic information is of relevance, is
oneofthemainfieldsinwhichthesepotentials have beenapplied.Intheearly days
of manybodypotentials themainalloy propertyfittedwas theheatof solutionof a
single impurity [6], i.e. the dilute limit of the heat of formation (HOF) of the alloy.
However, when these potentials are applied to concentrated alloys the predictions
are usually uncontrolled; they work well for systems with a mixing enthalpy that is
nearly symmetric and positive over the entire concentration range, as for example
in the cases of Fe–Cu [7, 8], or Au–Ni [6, 9, 10].
Alloys which show a strong asymmetry or even a sign inversion in the HOF such
ISSN:1478-6435print/ISSN1478-6443online
(cid:13)c 2009Taylor&Francis
DOI:10.1080/14786430903292373
http://www.informaworld.com
January31,2012 2:3 PhilosophicalMagazine prim
3372 B. Sadigh, P. Erhart, A. Stukowski, and A. Caro
as Fe–Cr or Pd–Ni are beyond the scope of standard many body potential models,
and there is not yet a unique methodology suitable for their description. Similar
limitations apply with respect to systems with a negative HOF which feature in-
termetallic phases. Frequently, such systems require different parametrisations for
different phases, as in the case of Ni–Al with the B2 phase on one hand [11], and
the γ and γ′ phases on the other [12].
Two schemes have been developed to deal with these shortcomings in the case
of Fe–Cr which displays an inversion in the HOF as a function of concentration,
namely the composition-dependent embedded atom method (CD-EAM) [13] and
the two-band model (2BM) [14]. For neither one of these schemes, it is obvious
how it can be extended to more than two components.
Theobjectiveofthispaperistodevelopaframeworkforconstructinginteratomic
potential models for multicomponent alloys based on an expansion in clusters of
increasing size that can be practically implemented and systematically improved.
Our methodology allows to describe systems with arbitrary heat of mixing curves
andincludesintermetallicphasesinasystematicandphysicallymeaningfulfashion.
Thereby, we overcome the most important disadvantages of current alloy potential
schemes and provide a framework for systems of arbitrary complexity.
In our methodology the interatomic interactions are modified by composition-
dependent functions. This introduces a dependence on the environment which is
somewhat reminiscent of the bond-order potential (BOP) scheme developed by
Abell and Tersoff [15–17]. In this formalism the attractive pair potential is scaled
by a (usually) angular dependent function (the “bond-order”) which describes the
local structure. Thereby, it is possible to distinguish different lattice structures
(face-centred cubic, body-centred cubic, cubic diamond etc.) and also to stabilise
structureswithlowpackingdensitysuchasdiamondorzincblendelattices.(Infact,
theBOPformalismhasbeensuccessfullyappliedtomodelalloyssuchasFe–Ptthat
featureintermetallicphaseswithdifferentlatticestructures[18]).Thecomposition-
dependent interatomic potential (CDIP) scheme introduced in the present work
and the BOP formalism thus both include explicit environment-dependent terms.
However,intheCDIPapproachthisenvironment-dependenceisusedtodistinguish
different chemical motifs while in the BOP scheme it is used to identify different
structural motifs.
This paper is organised as follows: In Sect. 2.1 we introduce the basic termi-
nology and present a systematic approach to fitting potentials for binary systems.
Section 2.2 describes how by including higher order terms it is possible to fit e.g.,
intermetallic phases. In Sect. 3.1 a series expansion is developed which generalises
the concepts introduced in the previous sections and which is used in Sect. 3.2
to obtain explicit expressions for a ternary system. The efficient computation of
forces is discussed in Sect. 4 and an optimal implementation in Monte-Carlo simu-
lations is the subject of Sect. 5. Finally, as an example, the composition-dependent
embedded atom method potential for Fe–Cr is reviewed in Sect. 6.
2. Binary Systems
2.1. Pair Potentials
Forthesakeofclarity ofthefollowingexposition,weassumeEAMmodelsthrough-
outthispaper.Itisimportanttostressthattheformalismtobedevelopedhereafter
can beappliedtoanypotential modelforthepureelements includingmodifiedem-
bedded atom method (MEAM) [19, 20], bond-order [15–17], and Stillinger-Weber
type [21] potentials.
January31,2012 2:3 PhilosophicalMagazine prim
Composition-Dependent Interatomic Potentials 3373
Considerasingle-componentsystemofatomsA,whoseinteractionsaredescribed
by the EAM model,
1
E = U (ρ )+ φ (r ) with ρ = ρ(r ). (1)
A A i 2 A ij i ij
i i j j6=i
X XX X
The first term in Eq. (1) contains the embedding function U (ρ ), which is a
A i
nonlinear function of the local electron density ρ around atom i. It accounts for
i
cohesion due to band formation in the solid state and is constructed to reproduce
the equation of state of system A.Thesecond term represents the remainder of the
interaction energy. It can be interpreted as the effective screened Coulomb inter-
action between pairs of ions in A. The EAM formalism can capture the energetics
associated with densityfluctuations inthelattice andhasbeensuccessfullyapplied
for modelling the formation of crystal defects such as vacancies, interstitials and
their clusters.
Consider now a binary system, where the pure phases are described by EAM
potentials. It can be shown that the total energy expression for this type of poten-
tials is invariant under certain scaling operations [22]. This “effective pair format”
can be used to rescale the two EAM potentials, e.g. such that at the equilibrium
volume for a certain lattice the electron density is 1, to ensure their compatibility.
One part of the total energy of the two-component system can be written as the
superposition of the respective embedding terms and effective pair interactions:
1
E = U ρA+µ ρB + φ (r ) (2)
0 A i A(B) i 2 A ij
i∈A i∈Aj∈A
X (cid:0) (cid:1) XX
1
+ U ρB +µ ρA + φ (r ),
B i B(A) i 2 B ij
i∈B i∈Bj∈B
X (cid:0) (cid:1) XX
where
ρS = ρS(r ). (3)
i ij
j∈S,j6=i
X
NotethatabovewehavenotyetaddedanyexplicitA−B interactions.Equation(2)
is a strict superposition of the interatomic potentials for the pure elements with
the only caveat that the electron density of the A (B) species in the embedding
functionofaB (A)particleisscaledwithaparameterµ inordertoaccountfor
B(A)
the different local electron densities. Thereby, two EAM models can be calibrated
with respect to each other. More elaborate schemes are possible, e.g. one can treat
µ and µ as free parameters. Here for the sake of simplicity, we restrict ourselves
A B
to normalised electron densities.
Starting from a parametrisation for E , we now devise a practical scheme for
0
systematically improving the interaction model. Let us denote the truemany-body
energy functional of the binary system by E . Our goal is to construct an inter-
t
atomic potential model for the difference energy functional ∆E(0) = E −E . We
t 0
begin with the two dilute limits. Consider a lattice of A particles and substitute
the atom residing in the i-th site with a B atom. Let us now assume that ∆E(0)
for this configuration can be satisfactorily represented by a pair potential between
January31,2012 2:3 PhilosophicalMagazine prim
3374 B. Sadigh, P. Erhart, A. Stukowski, and A. Caro
the A−B pairs. In this limit ∆E(0) can thus be written as
∆E(0)(A-rich)= VA (r ). (4)
AB ij
j∈A
X
(There is only one sum in this expression since we are dealing with a configuration
that contains only one B atom). A similar expression is obtained for the B-rich
limit
∆E(0)(B-rich) = VB (r ). (5)
AB ij
j∈B
X
Since we do not require the pair potential models for the two dilute limits to
coincide with each other, an interpolation is needed which preserves the energetics
of the impurities. The main objective of the present paper is to devise such an
interpolation scheme. The simplest ansatz for such an expression is
∆E(0) = xAVA (r )+ xBVB (r ) (6)
ij AB ij ij AB ij
i∈Aj∈B i∈Aj∈B
XX XX
Above, xS denotes the concentration of species S in the neighbourhood of an
ij
A −B pair residing on the i and j sites. Ideally, we require this quantity to be
easy to calculate and to be insensitive to the local density and topology, in other
words itshouldseparatechemistry fromstructure.Inany case, xS has torepresent
ij
an average over the neighbourhood of both centres i and j. Before we derive the
expression for xS,it is instructive to discuss the correspondingone-centre quantity
ij
xS. It describes the local concentration of species S around atom i. A simple way
i
to determine xS is to choose a local density function σ(r ) and then to evaluate
i ij
the following expression
xS = (j∈S,j6=i)σ(rij) = σSi , (7)
i σ(r ) σ
P j6=i ij i
P
which is indeed rather insensitive to the local geometry. This is most obvious in
the dilute limits. The local concentration xS at the site of an impurity atom i is
i
either 0(if S is theminority species)or 1 (if S is the majority species)independent
of the local structure. This is, however, strictly true only for the impurity atom.
For the other atoms in the system xS varies between 0 and 1 depending on the
j
distance to the impurity atom. Also for these particles, atomic displacements may
change the value of xS. A total decoupling of chemistry and structure is therefore
j
not possible. The optimal choice for σ(r ) is the function that minimises the effect
ij
of local geometry on xS. Although it is possible to choose different σ-functions for
i
the different species, we do not expect the quality of the final potential to depend
crucially on the choice of σ(r ). In fact, we expect the best choice for σ(r ) to be
ij ij
the simplest one.
It is now straightforward to extend Eq. (7) to define the concentration xS in
ij
the neighbourhood of a pair of atoms residing on sites i and j. To this end, we
first define a quantity xS to represent the concentration of the species S in the
i(j)
January31,2012 2:3 PhilosophicalMagazine prim
Composition-Dependent Interatomic Potentials 3375
cutoff for atom 2
x = 5 x = 5
1 8 1(2) 7 x = 1(x +x )=11=0.786
2 x2 = 78 x2(1) = 67 12 2 1(2) 2(1) 14
1
x = 3 x = 2
1 8 1(2) 7 x = 1(x +x )= 3 =0.214
x = 1 x = 1 12 2 1(2) 2(1) 14
2 8 2(1) 7
cutoff for atom 1
Figure 1. Schematic illustration of the connection between xS and two-centre concentrations xS and
i ij
their computation in a binary alloy according to Eqs. (7) and (8). Here, the cutoff function σ(r) which
appears inEq.(7)isassumedtobeastepfunctionwhichis1forr<rc andzerootherwise.
neighbourhood of atom i excluding atom j:
xS = (k∈S,k6=i,k6=j)σ(rik) = σSi −δ(S,tj)σ(rij) (8)
i(j) σ(r ) σ −σ(r )
P (k6=i,k6=j) ik i ij
P 1−σ(r )/σS
ij i t = S
j
=xS 1−σ(rij)/σi ,
i 1
tj = S
1−σ(r )/σ
ij i
where t denotes the type of atom i, and δ(t ,t ) is 1 if t = t and zero otherwise.
i i j i j
Using this quantity, the two-centre concentration xS can be defined as follows
ij
1
xS = xS +xS (9)
ij 2 i(j) j(i)
(cid:16) (cid:17)
Hence, the two-centre concentration of the species S about the atom pair (i,j)
is the average concentration in the two separate neighbourhoods of sites i and j
excluding both of these atoms. This definition, which is illustrated in Fig. 1, has
the important advantage that the interpolation scheme introduced in Eq. (6) does
not modify the interactions in the dilute limits, since xS is strictly 0 or 1 in the
ij
two limits irrespective of the local structure. Furthermore, it is straightforward to
generalise Eq. (9) to multi-centre concentrations. For example, in the next section,
we will explicitly discuss the construction of interatomic potentials using three-
centre concentrations.
Let us now revisit Eq. (6). As mentioned earlier this is the simplest ansatz for
∆E(0) that can reproduce the energetics of both dilute limits. A more general
expression is
∆E(0) = hA (xA) VA (r )+ hB (xB) VB (r ), (10)
AB ij AB ij AB ij AB ij
i∈Aj∈B i∈Aj∈B
XX XX
where hB(x) and hA are nonlinear functions with the property hB(0) = hA(0) =
A B A B
0 and hB(1) = hA(1) = 1. By fitting these functions to the energetics of the
A B
concentrated alloys, the quality of the interatomic potential model for the binary
can be improved drastically.
In principle, one can stop here and have an interatomic potential model, E +
0
∆E(0), that can reproduce the energetics of the dilute limits as well as the solid
solutionofthebinary.Itis,however,alsopossibletofurtherrefinetheabovemodel.
January31,2012 2:3 PhilosophicalMagazine prim
3376 B. Sadigh, P. Erhart, A. Stukowski, and A. Caro
For this purpose, let us again define a difference energy functional
∆E(1) = E −E −∆E(0), (11)
t 0
and construct an interatomic potential model for the energy functional ∆E(1).
Consider a lattice of A particles and substitute two atoms, say i and j, with B
particles. Assume that ∆E(1) for this configuration can be well represented by
a potential model describing the interaction of the B-B pair with a lattice of A
particles. In this limit we can express ∆E(1) as
∆E(1)(A-rich) = VA (r )+ VA (r ), (12)
BB ij BBA ijk
k
X
where r is shorthand for the three sets of positions of the i, j and k atoms, i.e.
ijk
{r ,r ,r }. In the same way we obtain for the B-rich limit
i j k
∆E(1)(B-rich)= VB (r )+ VB (r ). (13)
AA ij AAB ijk
k
X
Note that ∆E(1) has both a two-body and a three-body component and thus can
be decomposed as follows
∆E(1) = ∆E(1) +∆E(1) . (14)
pair triplet
In the next section we discuss how to incorporate the three-body contribution into
(1)
the interatomic potential model. For now, we only consider ∆E . Following the
pair
same line of arguments that lead to Eq. (10), we obtain the expression
∆E(1) = hA (xA) VA (r )+ hB (xB) VB (r ), (15)
pair BB ij BB ij AA ij AA ij
i∈Bj∈B i∈Aj∈A
XX XX
which reproduces the contributions of the pair terms in the two limits given by
Eqs. (12) and (13). The two non-linear functions have to fulfil the conditions
hB (0) = hA (0) = 0 (16)
AA BB
hB (1) = hA (1) = 1. (17)
AA BB
By fitting the functions hB and hA in the intermediate concentration range to
AA BB
the energetics of the concentrated alloy, one can obtain a further improvement for
the interaction model for the binary system.
2.2. Beyond Pair Potentials
In this section, we show that the formalism introduced in the previous section can
be extended to multi-body interaction potentials, which enables us to capture the
energetics of awider rangeof phasesincludingordered compounds.Inthe previous
section, we outlined a scheme to construct composition-dependent pair potentials
for the potential energy landscape E +∆E(0)+∆E(1). It was also observed that a
0
properformulation of∆E(1) requiresincorporation of explicit three-bodyterms.In
this section we describe how to construct such composition-dependent multi-body
potentials.
January31,2012 2:3 PhilosophicalMagazine prim
Composition-Dependent Interatomic Potentials 3377
cutoff for atom 2 x = 5 x = 5
1 8 1(23) 6
x = 7 x = 6 x = 1(x +x +x )= 113=0.897
2 8 2(13) 7 123 3 1(23) 2(13) 3(12) 126
2
x = 8 x = 7
3 8 3(12) 7
1
x = 3 x = 1
3 1 8 1(23) 6
x = 1 x = 1 x = 1(x +x +x )= 13 =0.103
2 8 2(13) 7 123 3 1(23) 2(13) 3(12) 126
x = 0 x = 0
cutoff for atom 3 3 8 3(12) 7
Figure2. Schematicillustrationofthecomputationofthree-centreconcentrationsinabinaryalloyusing
Eqs.(7)and(21).Here,thecutofffunctionσ(r)whichappearsinEq.(7)isassumedtobeastepfunction
whichis1forr<rc andzerootherwise.
First, we require an interpolation scheme to connect the two limits of the three-
(1)
body term ∆E in Eq. (14). The simplest ansatz for such an expression is
triplet
∆E(1) = xA VA (r )+ xB VB (r ), (18)
triplet ijk BBA ijk ijk AAB ijk
i∈Bj∈Bk∈A i∈Aj∈Ak∈B
XXX XXX
where xS denotes the concentration of species S in the neighbourhood of the
ijk
triplet residing on sites i, j and k. In analogy with the derivation of the two-centre
concentration Eq. (9), we start from the one-centre concentration xS and define
i
the intermediate quantity xS(jk) that represents the concentration centred around
i
atom i excluding atoms j and k
xS = (l∈S,l6=i,l6=j,l6=k)σ(ril) = σSi −δ(S,tj)σ(rij)−δ(S,tk)σ(rik) (19)
i(jk) σ(r ) σ −σ(r )−σ(r )
P (l6=i,l6=j,l6=k) il i ij ik
P1−[δ(S,t )σ(r )+δ(S,t )σ(r )]/σS
= xS j ij k ik i , (20)
i 1−[σ(r )+σ(r )]/σ
ij ik i
and now following the same line of arguments leading to Eq. (9) we define the
three-centre concentration xS as follows
ijk
1
xS = xS +xS +xS . (21)
ijk 3 i(jk) j(ik) k(ij)
(cid:16) (cid:17)
A graphical illustration of the computation of this quantity is given in Fig. 2. The
three-centre concentration of the species S about the triplet (i,j,k) is the average
concentration (excluding the triplet) in three separate neighbourhoods, each of
which is centred at one of the atoms in the triplet. Thanks to this definition xS is
ijk
strictly 0 or 1 in the two dilute limits described in Eqs. (12) and (13), irrespective
of the local structure. Hence, the interpolation scheme in Eq. (18) does not alter
the interactions in Eqs. (12) and (13). Again, as in Eq. (10), we can improve the
simple interpolation scheme in Eq. (18)
∆E(1) = hA (xA )VA (r )+ hB (xB )VB (r ),
triplet BBA ijk BBA ijk AAB ijk AAB ijk
i∈Bj∈Bk∈A i∈Aj∈Ak∈B
XXX XXX
(22)
where hA and hB are non-linear functions that can be fitted to the energetics
BBA AAB
of the concentrated alloys with the boundary conditions
hA (0) = hB (0) = 0 and hA (1) = hB (1) = 1. (23)
BBA AAB BBA AAB
January31,2012 2:3 PhilosophicalMagazine prim
3378 B. Sadigh, P. Erhart, A. Stukowski, and A. Caro
* * *
−embedded −embedded −embedded
cluster of order 2 cluster of order 3 cluster of order 4
( ) ( ) ( 2 )
(C ,C ) (C ,C ,C ) (C ,C2,C )
1 2 1 2 3 1 2 3
Figure3. SchematicillustrationofS-embeddedcolouredclustersoforders2,3,and4inaternaryalloy.
Theshadedregionindicatesthecutoffrangearoundthecentralatommarkedbyanasterisk.
Followingthisschemecomposition-dependentclusterinteractionsofarbitraryorder
can be included in the interatomic potential model. To summarise, to incorporate
cluster interactions of order n, two cluster potentials are constructed, one for the
configuration wherethecluster is embeddedin theA lattice and onefor the config-
uration where the cluster is embedded in the B lattice. Subsequently these limits
are interpolated using the n-centre concentrations. In the next section, we review
this strategy in detail to show that a systematic series expansion in composition-
dependent cluster interactions is possible for general multicomponent systems.
3. Multicomponent Systems
3.1. Series Expansion in Embedded Cluster Interactions
In the first sections of this paper, we have shown how to practically construct
interatomic potentials for binary systems. First, mixed interatomic pair andtriplet
potentials are generated for the dilute limits which are subsequently extended to
arbitrary concentrations by fitting interpolation functions that dependon the local
concentration about the atomic pairs and triplets. The choice of specific potentials
anddiluteconfigurationswasmainlydrivenbyphysicalintuition.Inthissectionwe
show that this procedure can be formalised and generalised to arbitrarily complex
systems with more than two components.
Consider an n-component mixture of N particles that are distinguishable only
through their species. Assign a unique colour to each of the species: {C ,...,C }.
1 n
We define a colour cluster of order m to be a set of m particles with a specific
colour combination. We use the occupation number formalism to identify colour
schemes, i.e. Ck1,...,Ckn , where k is the number of particles in the cluster with
1 n i
colour Ci, and(cid:16) iki = m.(cid:17)For example, a cluster of order 3 consisting of one parti-
cle with the colour C and two particles with the colour C , is denoted by C ,C2 .
1 3 1 3
P
Furthermore, we define an S-embedded colour cluster of order m to be a set of
(cid:0) (cid:1)
m coloured particles embedded in a pure matrix of species S. Three examples of
such S-embedded coloured clusters are shown in Fig. 3. The key idea is that the
potential energy landscape of an alloy can be expanded in the basis set of elemen-
tary interaction potentials each of which is constructed to reproducethe energetics
of a particular embedded colour cluster. The order of an interaction element in
the series is determined by the order of the corresponding colour cluster. By pro-
gressively includinghigher order colour cluster interactions, one can systematically
increase the accuracy of the model.
To recapitulate, we expand the potential energy landscape of multicompo-
nent systems in the basis set of colour cluster interatomic potential functions
January31,2012 2:3 PhilosophicalMagazine prim
Composition-Dependent Interatomic Potentials 3379
cutoff for atom 2 x = 3 x = 3
1 8 1(2) 7 x = 1(x +x )= 7 =0.500
x = 5 x = 4 12 2 1(2) 2(1) 14
2 8 2(1) 7
2
x = 3 x = 2
1 1 8 1(2) 7 x = 1(x +x )= 3 =0.214
x = 1 x = 1 12 2 1(2) 2(1) 14
2 8 2(1) 7
x = 2 x = 2
1 8 1(2) 7 x = 1(x +x )= 4 =0.286
cutoff for atom 1 x = 2 x = 2 12 2 1(2) 2(1) 14
2 8 2(1) 7
Figure 4. Schematic illustration of the connection between xS and two-centre concentrations xS and
i ij
their computation in a ternary alloy according to Eqs. (7) and (8). Here, the cutoff function σ(r) which
appears inEq.(7)isassumedtobeastepfunctionwhichis1forr<rc andzerootherwise.
VS ({r}), where {r} is the real-space configuration of the respective clus-
C1k1...Cnkn
ter. The expansion coefficient for each basis function is the interpolation function
hS (xS), where xS is the local concentration of the species S in the neigh-
C1k1...Cnkn
bourhood of the cluster. One of the innovations in this work is a simple and com-
putationally expeditious way to determine xS which is illustrated for the case of
a ternary alloy in Fig. 4. Formally the total energy expression for an alloy of n
components and N particles can be written as
E = E + ... hS (xS)VS ({r}), (24)
0 C1k1...Cnkn C1k1...Cnkn
Xm Xk1 Xkn XS
Pni=1ki=m
| {z }
where the first sum is over the order of the cluster potentials and the subsequent
sums are over all distinguishable colour combinations of m-size clusters. Each term
in the above expansion can be evaluated as follows
N N
hS (xS)VS ({r}) = ... hS (xS )VS (r ).
C1k1...Cnkn C1k1...Cnkn C1k1...Cnkn i1...im C1k1...Cnkn i1...im
iX1=1 iXm=1
{i1...im}∈{C1k1...Cnkn}
(25)
| {z }
The sums in Eq. (25) are over all possible m-size atom clusters {i ...i } in the
1 m
system with the colour scheme Ck1,...,Ckn .
1 n
Themain advantage of thissc(cid:16)hemeis thatt(cid:17)hebasisfunctionscan beconstructed
sequentially and independentof the interpolation functions. The lower order terms
can beconstructed with noknowledge of thehigherorder terms andthereforeneed
not be reparametrised when higher order cluster potentials are constructed. The
higher order terms in the expansion become progressively smaller. Furthermore,
addition of new terms in the series expansion is not likely to introduce unphysical
behaviour, a problem that plagues most fitting schemes for interatomic potentials.
3.2. Explicit expressions for ternary alloys
In this section we illustrate the formal discussion in the previous section by con-
structinganexpansioninembeddedpairandtripletpotentialsforaternarysystem.
For simplicity we assume the pure elements are described by EAM models. The
extension to larger number of components and higher order cluster potentials will
be obvious. We consider a system of three components A, B and C, and assume
January31,2012 2:3 PhilosophicalMagazine prim
3380 B. Sadigh, P. Erhart, A. Stukowski, and A. Caro
that three composition-dependent pair potentials for the binary systems A−B,
A−C and B−C have already been constructed. Explicitly, the A−B interaction
is given by the following expression
1
Epair = U ρA +µ ρB + hA (xA)φ (r )+hB (xB)VB (r )
A-B A i A(B) i 2 AA ij A ij AA ij AA ij
i∈A i∈Aj∈A
X (cid:0) (cid:1) XX(cid:0) (cid:1)
1
+ U ρB +µ ρA + hB (xB)φ (r )+hA (xA)VA (r )
B i B(A) i 2 BB ij B ij BB ij BB ij
i∈B i∈Bj∈B
X (cid:0) (cid:1) XX(cid:0) (cid:1)
+ hA (xA)VA (r )+hB (xB)VB (r ) . (26)
AB ij AB ij AB ij AB ij
i∈Aj∈B
XX(cid:0) (cid:1)
By now the notation above should be familiar. The interaction potentials for the
two other pairs can be written in analogous fashion.
Now, we can spell out the expansion in embedded pair potentials for the ternary
A−B−C
Epair = U ρA+µ ρB +µ ρC (27)
A−B−C A i A(B) i A(C) i
i∈A
X (cid:0) (cid:1)
+ U ρB +µ ρA+µ ρC
B i B(A) i B(C) i
i∈B
X (cid:0) (cid:1)
+ U ρC +µ ρA+µ ρB
C i C(A) i C(B) i
i∈C
X (cid:0) (cid:1)
1
+ hA (xA)φ (r )+hB (xB)VB (r )+hC (xC)VC (r )
2 AA ij A ij AA ij AA ij AA ij AA ij
i∈Aj∈A
XX(cid:2) (cid:3)
1
+ hB (xB)φ (r )+hA (xA)VA (r )+hC (xC)VC (r )
2 BB ij B ij BB ij BB ij BB ij BB ij
i∈Bj∈B
XX(cid:2) (cid:3)
1
+ hC (xC)φ (r )+hA (xA)VA (r )+hB (xB)VB (r )
2 CC ij C ij CC ij CC ij CC ij CC ij
i∈Cj∈C
XX(cid:2) (cid:3)
+ hA (xA)VA (r )+hB (xB)VB (r )+hC (xC)VC (r )
AB ij AB ij AB ij AB ij AB ij AB ij
i∈Aj∈B
XX(cid:2) (cid:3)
+ hA (xA)VA (r )+hB (xB)VB (r )+hC (xC)VC (r )
AC ij AC ij AC ij AC ij AC ij AC ij
i∈Aj∈C
XX(cid:2) (cid:3)
+ hA (xA)VA (r )+hB (xB)VB (r )+hC (xC)VC (r ) .
BC ij BC ij BC ij BC ij BC ij BC ij
i∈Bj∈C
XX(cid:2) (cid:3)
The only unknowns in the above equation are VC (r ), VB (r ), VA (r ),
AB ij AC ij BC ij
hC (xC), hB (xB) and hA (xA). The potentials VC (r ), VB (r ) and VA (r )
AB ij AC ij BC ij AB ij AC ij BC ij
describe the interaction between pairs of unlike species embedded in pure lattices
ofthethirdspeciesoftheternary.Inanalogy withtheprevioussection,itisreason-
able to expect that we can construct these potentials separately in their respective
dilute limits and subsequently fit the interpolation functions hC (xC), hB (xB),
AB ij AC ij
hA (xA) to the energetics of the concentrated ternary alloys. However, when the
BC ij
number of species increases certain complications can arise that are not present in
the binaries. This is well illustrated in the situation above. We now show that it
is in fact not possible to separately construct the three pair potentials VC (r ),
AB ij
