Calculations of excess free energies of precipitates via direct thermodynamic integration across phase boundaries
aa r X i v : . [ c ond - m a t . m t r l - s c i ] S e p Calculations of excess free energies of precipitatesvia direct thermodynamic integration across phase boundaries
Babak Sadigh ∗ and Paul Erhart
1, 2, † Lawrence Livermore National Laboratory, Condensed Matter and Materials Division, Livermore, California, 94551, USA Chalmers University of Technology, Department of Applied Physics, Gothenburg, Sweden (Dated: June 12, 2018)We describe a technique for constraining macroscopic fluctuations in thermodynamic variableswell-suited for Monte Carlo (MC) simulations of multiphase equilibria. In particular for multi-component systems this amounts to a statistical ensemble that implements constraints on boththe average composition as well as its fluctuations. The variance-constrained semi-grandcanonical(VC-SGC) ensemble allows for MC simulations, in which single-phase systems can be reversiblyswitched into multiphase equilibria allowing the calculation of excess free energies of precipitates ofcomplex shapes by thermodynamic integration. The basic features as well as the scaling and conver-gence properties of this technique are demonstrated by application to an Ising model. Finally, theVC-SGC MC simulation technique is used to calculate α/α ′ interface free energies in Fe–Cr alloysas a function of orientation and temperature taking into account configurational, vibrational, andstructural degrees of freedom. PACS numbers: 02.70.Tt 05.70.Np 68.35.bd 81.30.Mh
I. INTRODUCTION
Multiphase systems are ubiquitous in nature. Sincetheir properties are crucially dependent on the interfacesbetween the phases involved interface free energies con-stitute crucial input to kinetic and phase field models. Their calculation, however, often constitutes a compli-cated task.
Here, an extension of conventional ther-mostatistical ensembles. Is described that enables calcu-lation of excess free energies via thermodynamic integra-tion across phase boundaries. The approach presented inthis paper can be viewed as a generalization of the ex-tended Gaussian ensemble technique to multicomponentsystems.
An equilibrium ensemble in statistical mechanics is acollection of macroscopic systems that have been pre-pared in the same thermodynamic state. This can beaccomplished by constraining one of each pair of con-jugate thermodynamic variables that represent gener-alized displacements and forces, e.g. volume–pressure,entropy–temperature, or particle number–chemical po-tential. In single-phase regions of the phase diagram,there is a one-to-one correspondence between general-ized displacements and forces. This allows for straightfor-ward computation of free energy differences by thermo-dynamic integration. In contrast, in two-phase regions ofthe phase diagram, the mapping from generalized forcesto generalized displacements becomes multi-valued. Forinstance the melting temperature maps to two distinctentropies, namely those of the solid S s and the liquid S l , respectively. Multiphase equilibria can be obtainedby constraining the entropy to values S s < S < S l . Ifthe conjugate forces are observables of the ensemble thisenables thermodynamic integration across phase bound-aries. Conveniently, two important conjugate forces,namely temperature and pressure, are readily availablein the microcanonical ( N V E ) ensemble. In contrast the chemical potential is not an observable of ensembles thatcontrol the number of particles or chemical compositionin multicomponent systems.In this paper a technique is described for construct-ing ensembles, which constrain fluctuations of generalizeddisplacements, that can be viewed as an extension andgeneralization of the Gaussian ensemble technique.
Such variance-constrained (VC) ensembles provide di-rect access to multiphase regions of the phase diagramwhile allowing for observation of both generalized forcesand displacements. Thereby they enable thermody-namic integration across phase boundaries and compu-tation of free energies for compositions inside miscibil-ity gaps. In this paper this technique is employed tostudy precipitation in an Ising model system and a bi-nary Fe–Cr alloy. Through Monte Carlo simulations inthe VC semi-grandcanonical (SGC) ensemble chemicalpotential–concentration isotherms in the shape of thevan-der-Waals loop are obtained. In this fashion, onecan extract interface free energies as a function of tem-perature and orientation. It is also demonstrated thatinterface free energies can be calculated by monitoringthe growth and/or shrinkage of compact precipitates.The paper is organized as follows. After recapitu-lating key features of the related canonical and semi-grandcanonical ensembles, the next section introducesthe VC-SGC ensemble and discusses strategies for sam-pling this ensemble using Monte Carlo (MC) simula-tions. The computation of interface free energies isdemonstrated using a simple first-neighbor Ising modelin Sect. III A, followed by an analysis of convergence andscaling properties in Sect. III B. Finally, in Sect. IV theVC-SGC MC approach is applied to extract interface freeenergies in Fe–Cr alloys as a function of temperatureand orientation while taking into account both config-urational and vibrational degrees of freedom. The paperis concluded in Sect. V.
II. ENSEMBLES FOR MULTI-COMPONENTSYSTEMS
We start this section with a review of the canonicaland semi-grandcanonical ensembles before deriving thevariance-constrained semi-grandcanonical (VC-SGC) en-semble. We conclude by discussing practical aspects per-taining to sampling the VC-SGC ensemble using MonteCarlo simulations.
A. Canonical and semi-grandcanonical ensembles
Let us start by a brief discussion of multiphase equi-libria in multicomponent systems. For the sake of sim-plicity and clarity, we restrict ourselves to an immisciblebinary alloy. The generalization to systems with an ar-bitrary number of species is straightforward. Considera system of N particles confined in a box of volume V where each particle carries a spin of value 0 or 1. Anarbitrary configuration is specified by ( x , σ ), where x isa 3 N -dimensional vector describing the position of everyparticle in the box, and σ is an N -dimensional spin vec-tor. The number of spin 1 particles is n = P Ni =1 σ i , andtheir concentration is c = n/N . We denote the energyof a configuration by ˆ U ( x , σ ). The canonical partitionfunction for this system at temperature T reads Z C ( c, N ) = 1 n !( N − n )! Z exp h − β ˆ U ( x , σ ) i d x , (1)where β = 1 /k B T and N = { N, V, T } is the set of inde-pendent thermodynamic variables. Using the identity1 n !( N − n )! = 1 N ! X σ ··· σ N δ N X i =1 σ i − n ! , (2)we can rewrite the partition function in terms of an effec-tive potential where the compositional degrees of freedomhave been integrated out: Z C ( c, N ) = 1 N ! Z exp [ − βF C ( x ; c, N )] d x . (3) F C ( x ; c, N ) = − k B T ln X σ ··· σ N δ N X i =1 σ i − n ! × exp h − β ˆ U ( x , σ ) i . (4) F C ( x ; c, N ) is the effective potential energy landscape forthe structural degrees of freedom of the multicomponentsystem at temperature T . In this way, we separate thecompositional and topological degrees of freedom. In thefollowing, we focus our discussion on the contribution ofthe compositional degrees of freedom to the free energy.We thus consider a frozen lattice of particles residingat x . For brevity we drop all references to x and denotethe potential energy of a particular spin configuration as ˆ U ( σ ). The SGC ensemble represents a set of configura-tions that sample the spin degrees of freedom accordingto the Boltzmann distribution while also allowing the to-tal concentration to vary. To tune the average concen-tration an external chemical potential is applied to thesystem, which corresponds to modifying the potential en-ergy function as followsˆ U S ( σ ; ∆ µ ) = ˆ U ( σ ) + ∆ µN ˆ c ( σ ) (5)ˆ c ( σ ) = N X i =1 σ i /N, (6)where ∆ µ is the relative chemical potential that controlsthe average concentration and the free energy is given by F S (∆ µ, N ) = − k B T ln X σ ··· σ N exp h − β ˆ U S ( σ ; ∆ µ ) i . (7)Inserting Eq. (5) into the above equation, the SGC par-tition function defined as Z S = exp ( − F S /k B T ) can beexpressed in terms of the canonical free energy as Z S (∆ µ, N ) = Z exp [ − β ( F C ( c, N ) + ∆ µN c )] dc. (8)The integrand in the above equation can be used to de-fine a concentration distribution function that is peakedaround the average concentration h ˆ c i S . The condition ofzero derivative at h ˆ c i S yields the well-known thermody-namic relation ∆ µ = − N ∂F C ∂c ( h ˆ c i S , N ) . (9) B. The variance constrained semi-grandcanonicalensemble
Since for immiscible systems, one value of ∆ µ mapsto several compositions inside the miscibility gap stablemultiphase coexistence cannot be established in the SGCensemble. To overcome this restriction we modify theSGC ensemble in such a way as to control concentrationfluctuations inside the miscibility gap. This is most eas-ily accomplished by adding a constraint that fixes theensemble-averaged squared concentration (cid:10) ˆ c (cid:11) . We callthis the variance-constrained semi-grand canonical (VC-SGC) ensemble.This approach is analogous to umbrella sampling whencalculating free energy barriers for reactions in molecularsystems. In this method, given a particular reactioncoordinate ξ , a harmonic external potential u b (cid:0) ξ ; ξ (cid:1) = K/ (cid:0) ξ − ξ (cid:1) is applied to bias the system toward posi-tions ξ while at the same time the fluctuations are re-strained by the force constant K . In the same spirit wecan write the potential energy function of the VC-SGCensemble asˆ U V ( σ ; φ, κ ) = ˆ U ( σ ) + κ ( N ˆ c ( σ ) + φ/ κ ) , (10)where φ and κ are now thermodynamic variables thatcontrol the average concentration as well as its fluctua-tions. The harmonic external potential is parametrizedsuch that ˆ U V → ˆ U S except for a constant shift when-ever κ → φ → ∆ µ . κ can be given a physicalintrepretation as the generalized force that controls con-centration fluctuations. This can be realized by puttingthe system in contact with a finite reservoir. Hencethe VC-SGC ensemble, like the closely related extendedGaussian ensemble, is only applicable to finite systems.As all fluctuations vanish in the thermodynamic limit, sodoes κ . It should thus diminish as the system size grows.This poses no difficulty in our study of interfacial free en-ergies since they also tend to zero in the thermodynamiclimit.Analogously to the SGC ensemble [see Eqs. (7) and (8)]we can express the VC-SGC partition function in termsof the canonical free energy as Z V ( φ, κ, N ) = Z exp (cid:2) − β (cid:0) F C ( c, N ) + ∆ U b V (cid:1)(cid:3) dc, (11)∆ U b V = κ ( N c + φ/ κ ) , (12)where ∆ U b V corresponds to the harmonic bias potentialof the umbrella sampling approach. The integrand inEq. (11) describes the probability distribution of theglobal concentration in the VC-SGC ensemble. It is apeaked function around the average concentration h ˆ c i V .The condition of zero derivative at h ˆ c i V yields the fol-lowing relation between φ and κ and the canonical freeenergy for the VC-SGC ensemble φ + 2 N κ h ˆ c i V = − N − ∂F C /∂c ( h ˆ c i V , N ) . (13)It is important to note that the right-hand side ofEq. (13) is the concentration derivative of the canoni-cal free energy at h ˆ c i V , which is a simulation observable.Hence no more unbiasing is necessary when employingthe above technique for free energy integration. How-ever, for very small system sizes where the condition ofzero derivative at h ˆ c i V may not hold, the simple proce-dure outlined above must be modified. The standardapproach in such cases is to extract the probability dis-tribution of the system as a function of concentration viahistogram methods and sampling bias must be cor-rected for as in the usual implementations of umbrellasampling. .Notwithstanding, the key advantage of the VC-SGCensemble is that while the chemical driving force (right-hand side of Eq. (13)) is an ensemble observable, theconstraint on concentration fluctuations simultaneouslyallows multiphase equilibria to be stabilized. This is sim-pler and more efficient than the state-of-the-art methodsfor calculating concentration-derivatives of the free en-ergy within the canonical ensemble, which measure thework required to perform a gradual/instantaneous trans-mutation of one species into the other while keeping theparticle number fix. While quite clever path-samplingalgorithms with optimized estimators for calculating this work have been devised, for the particular applicationdiscussed in this paper, the variance-constrained ensem-ble presents a simpler and more efficient route, since thefree energy derivative can be obtained as an observable ofthe equilibrium ensemble, and no calculation of externalwork is required. Furthermore as discussed in detail in aprevious publication, the VC-SGC MC method allowsfor faster convergence to equilibrium in multiphase re-gions of the phase diagram than the canonical ensemble.Finally the VC-SGC ensemble is also very well suited forefficient parallel MC algorithms that enable simulationsof very large systems containing millions of particles. C. Efficient sampling of the VC-SGC ensemble
When carrying out MC simulations in the SGC en-semble according to Eq. (9), the parameter ∆ µ is di-rectly proportional to the free energy derivative ∂F/∂c ,which renders choosing suitable parameters straightfor-ward. We will now demonstrate that choosing the rele-vant range of values for VC-SGC simulations is just assimple as in the case of the SGC ensemble.In the previous section the VC-SGC ensemble was de-rived in terms of the parameters φ and κ , which led toexpressions that very much resembled the SGC ensemble.In practice it turns out that it is more convenient to usethe substitutions κ = ¯ κ/N and φ = ¯ κ ¯ φ . The VC-SGCpotential defined above in Eq. (10) then readsˆ U V ( σ ; ¯ φ, ¯ κ ) = ˆ U ( σ ) + ¯ κN (cid:0) c + ¯ φ/ (cid:1) (14)and the associated expression for the first derivative ofthe free energy [compare Eq. (13)] is¯ κ (cid:0) ¯ φ + 2 h ˆ c i V (cid:1) = − N − ∂F C /∂c ( h ˆ c i V , N ) . (15)The above transformations effectively decouple the av-erage condtraint from the variance constraint parame-ter, which simplifies parameter selection for VC-SGC MCsimulations as will be seen in the following.It is instructive to first consider a symmetric free en-ergy profile, in which case ∂F/∂c = 0 at the solubilitylimits where F is minimal as well as around 50% where F is maximal. From Eq. (15) it follows that for non-zero¯ κ one can thus install a concentration of 50% by choosing¯ φ = −
1. Furthermore one obtains h ˆ c i V → φ ≈ − h ˆ c i V → φ ≈
0. These observations suggest asimple protocal for an efficient sampling of the full con-centration range: First choose a constant value for ¯ κ andthen vary ¯ φ in the range − . < ∼ ¯ φ < ∼ .
2. In our ex-perience this approach is very transferable and we haveused it successfully in exactly this fashion for various sys-tems including for example the (very asymmetric) Fe–Crsystem as described in Sect. IV.In the following section it will be shown among overthings that the results for the free energy derivative areinsensitive to the variance constraint parameter ¯ κ over a T e m pe r a t u r e F r ee ene r g y de r i v a t i v e ¶ F / ¶ c SGCVC−SGCfit of g prec −2 0 2 flatinterface(a) sphericalprecipitate ¶D F xc / ¶ c (cid:181) g prec c −1/3 F r ee ene r g y F Concentration (%) D F xc (cid:181) g flat (b) D F xc (cid:181) g prec c FIG. 1. (Color online) (a) Derivative of free energy with re-spect to concentration at T = 1 . k B for system size 8 × × k B . (b) Free energy profileobtained by integration of data in (a). wide range of values. This provides for a lot of freedomin choosing hatκ .It is straightforward to formulate a Monte Carloalgorithm for sampling the VC-SGC ensemble, wheretransmutation trial moves comprise1. selecting a particle at random2. flipping its spin (type) and3. computing the energy change ∆ E and the concen-tration change ∆ c These trial moves are then accepted with probability A V = min (cid:8) , exp (cid:2) − β (cid:0) ∆ E + ¯ κN c (cid:0) ¯ φ + ∆ c + 2 c (cid:1)(cid:1)(cid:3)(cid:9) , (16)which satisfies detailed balance. III. BASIC FEATURES OF THE VC-SGC MCTECHNIQUEA. Interface free energies from directthermodynamic integration
In this section we illustrate the utility of the VC-SGCensemble by studying phase segregation in a simple Ising I n t e r f a c e f r eeene r g y Temperature ( k B )(100)(110)(111)spherical FIG. 2. (Color online) Interface free energies as a function oftemperature and orientation for the Ising model described inthe text. The inset shows a snapshot of a { } interface froma VC-SGC MC simulation close to 50% composition (picturerendered with ovito , Ref. 17). model. We will calculate excess free energies of interfacesbetween two coexisting phases by direct thermodynamicintegration from the single-phase region of the phase dia-gram into the miscibility gap. The Ising model is definedon a body-centered cubic (BCC) lattice with its Hamil-tonian parametrized as follows H = − X ij S i S j , (17)where the summation runs over first-nearest neighborpairs and S i ∈ {− , } . The phase diagram of this systemis shown in the inset of Fig. 1(a).Our first goal in this section is to show that the VC-SGC ensemble, unlike the SGC ensemble, allows for cal-culation of free energy derivatives in both single and mul-tiphase regions of the phase diagram. To this end, sim-ulations were carried out using cells with basis vectorsoriented along [100], [010], and [001]. In the VC-SGCsimulations the constraint ¯ φ (controlling average concen-tration) was varied from − .
05 to 0.05 in steps of 0.01at a constant variance constraint of ¯ κ = 100. In caseof SGC simulations the chemical potential difference ∆ µ was changed between − φ (or ∆ µ ) the system was equilibrated for 5 × MCsweeps, after which statistics were gathered for 15 × MC sweeps.The thus obtained relations between free energyderivative ∂F/∂c and concentration c are shown inFig. 1(a) for a temperature of 1 . k B . The gap in theblue line in Fig. 1(a) corresponds to the miscibility gap.This reflects the inability of the SGC ensemble in sta-bilizing multiphase equilibria. The two concentrationsdelineating the boundary of the miscibility gap, i.e. 1%and 99% in Fig. 1, are commonly referred to as binodals.Using the VC-SGC MC method F ′ ( c ) can be calcu-lated for all concentrations as shown by the yellow linein Fig. 1(a). Three regimes can be distinguished: Atconcentrations between the pure phases and the binodals(regime I), single-phase equilibria are obtained and the n z fi
10 9 8 7 6 5 4 3 M a x i m u m o f e xc e ss f r ee ene r g y Inverse linear dimension 1/ n z T = 1.0 k B T = 1.5 k B T = 2.0 k B T = 2.5 k B FIG. 3. (Color online) Scaling of maximum of excess freeenergy [compare Fig. 1(b)] with inverse linear dimension forthe Ising model described in the text. Simulation cell vectorsare oriented along h i directions inducing the formation of { } interfaces. SGC-MC results are reproduced. In this regime, F ′ ( c )has a positive slope. This behavior is preserved for com-positions between binodal and critical nucleation (ex-trema in ∂F/∂c in Fig. 1(a)), which only support for-mation of subcritical minority-phase clusters (regime II).In this regime there is a finite thermodynamic drivingforce for phase segregation but due to small sizes of theminority-phase clusters the interfacial free energy dom-inates and leads to their overall formation free energiesbeing positive. Finally at concentrations beyond criticalnucleation (regime III), supercritical clusters are formedand F ′ ( c ) is a monotonically decreasing function. In thisregime the excess free energy increases due to the in-crease in the interface area of the cluster. Note that thevan-der-Waals loop obtained in Fig. 1(a) is the equilib-rium solution under the constraint of finite system size.In particular, the concentration range spanned by regimeII diminishes as the simulation cell increases, and even-tually vanishes in the thermodynamic limit.Studing the shapes of the supercritical clusters inregime III, we find several subregimes corresponding todifferent types of precipitates. For concentrations imme-diately beyond critical nucleation, minority-phase clus-ters are compact. The excess free energy (proportionalto the interface area) as a function of concentration canthus be written as∆ F xc ( c ) = ξ [ V cell ( c − c )] / γ prec . (18)Here, ξ is a shape factor, which equals (36 π ) for spheri-cal precipitates, and V cell denotes the volume of the simu-lation cell. The bold dashed green lines in Fig. 1 indicatefits to Eq. (18) and its concentration derivative. The ef-fective interface free energies calculated in this way areshown as a function of temperature in Fig. 2.At still higher concentrations, the precipitates becomeso large that spherical clusters can no longer be containedinside the finite-sized simulation cell and a transition to S t anda r d de v i a t i on o f c on c en t r a t i on ( % ) Variance constraint − k x −1/2 · ·
8, 20%4 · ·
8, 50%8 · ·
8, 20%8 · ·
8, 50%(a) 10 Number of atoms N − k = 10 − k = 100 − k = 1000(b) FIG. 4. (Color online) Standard deviation of concentrationat a temperature of T = 2 k B as a function of (a) variance con-straint parameter and (b) system size. Data in (b) measuredat 50% composition. The black dashed line indicates scalingwith 1 / √ x . a cylindrical shape induced by periodic boundary condi-tions is observed. While the free energy profile recordedin this region is of limited physical interest per se , it pro-vides the continuous derivative needed for integration ofthe full excess free energy.Finally in the concentration range around 50% inFig. 1(a), the excess free energy assumes a constant valuedue to formation of superlattices consisting of alternatingslabs of the two phases, which follow the periodicity ofthe simulation cell. The excess free energy at these con-centrations corresponds to the free energy cost associatedwith two flat interfaces with the free energy density givenby γ flat = ∆ F xc /A, (19)where A is the cross-sectional area of the computationalcell. By changing the geometry of the simulation cell itis possible to change A and thereby interface orientation.To see this remember that during the course of a simula-tion the system will establish the interface configurationthat minimizes ∆ F xc not γ flat . Since the cross section A for different orientations can be changed by varying shapeand size of the simulation cell, it is possible to stabilizedifferent interface orientations as exemplified in Fig. 2.To obtain for example the free energy density of { } interfaces, we used an orthorhombic simulation cell withlattice vectors oriented along [1¯10], [¯1¯12], and [111] and2 × × F c for { } and { } interfaces is larger thanfor { } interfaces. Similarly shortening the dimensionsperpendicular to [1¯10] and extending the cell parallel tothis direction stabilizes a { } interface. In this fashionit is possible to not only extract the temperature but alsothe orientation dependence of the interface free energy,which leads to the data displayed in Fig. 2. To verify the reliability of the thus computed interfa-cial free energies using the VC-SGC MC technique, wehave also studied the dependence of the excess free en-ergies of the two-phase equilibria inside the miscibilitygap on the size of the simulation cell. Since the excessfree energy is dominated by the interfacial free energy itshould scale inversely with the longest linear dimensionof the simulation cell. This behavior is indeed observedas shown in Fig. 3. B. Convergence and parameter sensitivity
In this section we study the convergence properties ofVC-SGC MC simulations and the sensitivity of the cal-culated free energies to the choice of the variance con-straint parameter. Figure 4 shows the dependence of thestandard deviation of concentration on (a) variance con-straint parameter and (b) system size. In both cases thestandard deviation scales inversely with the square rootof the respective parameter. To understand this behav-ior consider again Eq. (11). The standard deviation ofconcentration in the VC-SGC ensemble can be expressedas the second moment of its probability distribution asfollows s V = 1 Z V Z c exp (cid:2) − β (cid:0) F C ( c, N ) + ∆ U b V (cid:1)(cid:3) dc − h ˆ c i V . (20)For large enough system sizes, the ensemble prob-ability distribution becomes a normal distribution p A/π exp (cid:0) − A ( c − B ) (cid:1) , with its mean B = h ˆ c i V , andits width A inversely proportional to the standard devi-ation: s V = 1 / A . Now the Gaussian width A is relatedto the second derivative of the distribution taken at theaverage concentration. We thus obtain the following re-lation for s V βs V = ∂ F C ( h ˆ c i V , N ) ∂c + 2¯ κN (21)This is a very important relation that clearly predictsthe dependence of the standard deviation of concentra-tion in the VC-SGC ensemble on the variance constraintparameter ¯ κ as well as on the system size N as observedin Figs. 4(a) and 4(b). It also provides a clear picture ofwhy and under which conditions the VC-SGC ensemblecan stabilize multiphase equilibria. Equation (21) illus-trates that whenever ¯ κ → s V diverges for values of h ˆ c i V , for which the second deriva-tive of the free energy function F C ( c, N ) becomes nega-tive. For these compositions, stable single-phase equilib-ria cannot be achieved. Introducing a non-zero varianceconstraint ¯ κ leads to finite s V . However, in order for ¯ κ to stabilize compositions inside the miscibility, it has tobe chosen large enough to make the right-hand side ofEq. (21) positive. Hence for practical calculations, the −1 0 1 F r ee ene r g y de r i v a t i v e ¶ F / ¶ c T = 24 · · − k = 10 − k = 100 − k = 1000(a) 0 2 4 6 0 20 40 60 80 100 S D V o f c on c en t r a t i on ( % ) Concentration (%) k = 5 k = 10 k = 20 k = 100 k = 1000(b) FIG. 5. (Color online) Effect of variance constraint parame-ter ¯ κ on (a) free energy derivative and (b) standard deviationof concentration at T = 2 k B for system size 4 × × best choice for ¯ κ is a value slightly larger than the max-imum of N | F ′′ C ( c, N ) | inside the miscibility gap. Ac-cording to Fig. 4(a) this threshold is about ¯ κ = 20 for asystem with 4 × × × × κ → κ → ∞ and/or N → ∞ concentration fluctuations go to zeroand thus approach the canonical ensemble.Figure 5(a) illustrates the important fact that the freeenergy derivative obtained via Eq. (13) is indeed unaf-fected by the strength of the variance constraint over awide range of parameter values. The lower limit of theacceptable parameter range is set by the divergence ofthe standard deviation that was described in the previ-ous paragraph and that is shown as a function of concen-tration in Fig. 5(b). The upper limit, however, is due todecreasing acceptance probability of the MC trial moves,which as illustrated in Fig. 6 diminishes only slightly with¯ κ over several orders of magnitude before dropping morerapidly. As a result for large values of ¯ κ one can no longergather sufficient statistics to compute meanigful thermo-dynamic averages. The range within which ¯ κ can be se-lected nonetheless extends at least over two to three or-ders of magnitude depending on system size. Since large¯ κ values lead to a reduction of the acceptance probabil-ity, it is, however, recommended to choose ¯ κ as small as
0 5 10 15 20 25 0 20 40 60 80 100 A cc ep t an c e P r obab ili t y ( % ) Concentration (%) system size4 · · k = 10 k = 100 k = 1000(a) 2 5 20 50 10 10 A cc ep t an c e P r obab ili t y ( % ) Variance constraint parameter − k · · · · · · FIG. 6. (Color online) Effect of variance constraint parame-ter ¯ κ on (a) free energy derivative and (b) standard deviationof concentration at T = 2 k B for system size 4 × × possible while maintaining an acceptable standard devi-ation. Figure 4 suggests that a reasonable value is about1%. IV. INTERFACE FREE ENERGIES IN THEFE–CR ALLOY SYSTEM
Up until now we have considered a simple Ising model.It can be considered as the most simple representativeof alloy cluster expansions, which are widely usedto study complex multi-component alloys. The VC-SGCMC technique is, however, equally applicable to empiricalpotential models and then enables computation of excessfree energies taking into account not only configurationalbut also structural and vibrational degress of freedom. To illustrate this point, we now consider the calculationof excess free energies in BCC the Fe–Cr binary alloy sys-tem as described by the concentration-dependent embed-ded atom method.
The latter reproduces the phasediagram below the Curie temperature featuring both sub-stantial solubility of Cr in Fe ( α -phase) and a pronouncedmiscibility gap as shown in the inset of Fig. 7(a).Simulations were carried out using orthorhombic cellscontaining between 384 and 480 atoms depending onsample orientation. The variance constraint parameterwas set to ¯ κ =100 eV/atom and ¯ φ was varied between − . ¶ F / ¶ c ( m e V / a t o m ) fit of g prec ¶ F / ¶ c (SGC) ¶ F / ¶ c (VC−SGC)−300 0 300 flatinterface(a) a ’−precipitate a −precipitate F ( m e V / a t o m )
0 30 60 excess free energy(b) D F xc ( m e V / a t o m ) Cr concentration (%)
0 10 20 30 0 20 40 60 80 100 D F xc (cid:181) g flat (c) D F xc (cid:181) g prec c D F xc (cid:181) g prec c T e m pe r a t u r e ( K ) a + a ’ FIG. 7. (a) Derivative of free energy with respect to con-centration obtained from MC simulations in the SGC andVC-SGC ensembles at 1000 K. The inset shows the calcu-lated phase diagram, which is in excellent agreement with thepublished phase diagram of the potential used in the presentwork. (b) Free energy profile obtained by integration of thedata in (a) and (c) the corresponding excess free energy. Thecentral dark gray area indicates the region within which flatinterfaces are obtained while the dark gray areas to the leftand right show the regions in which compact precipitates areobtained. thermal vibrations were sampled using displacement MCtrial moves, while volume trial moves were employed tosample zero pressure conditions. The system was evolvedfor 6,000 MC steps at each value of ¯ φ . Statistics weregathered over the last 5,000 steps.First we mapped out ∂F/∂c using SGC MC simula-tions as shown for a temperature of 1000 K by the boldblue line in Fig. 7(a). The gap in the blue line corre-sponds to the miscibility gap, which again reflects the in-ability of the SGC ensemble to stabilize multiphase equi-libria. Using the VC-SGC MC method we then calcu-lated ∂F/∂c over the entire concentration range as shownby the yellow line in Fig. 7(a). The free energy profileobtained in this fashion again coincides with the SGCresults in the single phase regions and inside the misci-
100 200 300 400 400 600 800 1000 1200 1400 1600 I n t e r f a c e f r ee ene r g y ( m J / m ) Temperature (K) {100}{110}{111}spherical precipitate
FIG. 8. Interface free energies for Fe–Cr model alloy for dif-ferent interface orientations as well as spherical precipitates. bility gap shows the same features that were discussed inSect. III A for the Ising model.Once again it is possible to extract free energy densi-ties corresponding to spherical and flat interfaces fromthe excess free energy curves, as shown in Fig. 7(c). Theresults of this analysis are summarized in Fig. 8. Notethat in this case the computation of zero-K interface en-ergies is significantly more complicated than in the caseof the Ising model. The asymmetric phase diagram fea-tures a large solubility of Cr in α -Fe that remains finite astemperature goes to zero and is accompanied by short-range ordering. We have therefore not attempted toconstruct a zero-K α/α ′ interface.Two features in Fig. 8 are particularly noteworthy incomparison with Fig. 2. The ordering of the interfaceorientations is different with { } interfaces now beinglowest in energy, and the anisotropy is much smaller.It is furthermore instructive to separate the differentcontributions to the interface free energy. First, by re-peating the calculations for { } interfaces without dis-placement trials moves it is possible to separate vibra-tional and configurational degrees of freedom. The resultof this analysis is shown in Fig. 9(a), which illustrates themagnitude of the vibrational contribution to the interfacefree energy. Since the simulations also provide averageinternal energies one can furthermore separate explicitlythe interface internal energies and entropies. As shownin Fig. 9(b) and (c) the internal energies are hardly af-fected if vibrations and atomic relaxations (other thanvolume changes) are suppressed. Instead their contribu-tions show up as a practically temperature independentcontribution to the entropy. Figures 9(a) and (b) finallyshow that interface free and internal energies extrapolateto a value of about 372 ± at zero K. V. CONCLUSIONS
In this paper we have introduced a free energy integra-tion technique along compositional degrees of freedomthat allows for efficient calculations of free energies of
200 300 400 F r ee ene r g y ( m J / m ) w/o displacementswith displacementsextrapolation(a) 400 500 600 E ne r g y ( m J / m ) (b) 0.0 0.1 0.2 0 200 400 600 800 1000 1200 1400 1600 E n t r op y ( m J / m / K ) Temperature (K) (c)
FIG. 9. (a) Free energy, (b) internal energy, and (c) entropyof { } interface in Fe–Cr model alloy. The dashed linesrepresent fits to second degree polynomials. multiphase systems with various degrees of heterogene-ity. The construction of the underlying VC-SGC ensem-ble represents a generalization of the extended Gaussianensemble technique to multicomponent systems. To illustrate the basic features of the VC-SGC MCtechnique as well as aspects such as scaling and con-vergence, an extensive characterization of a first-nearestneighbor Ising model was carried out. In this way it wasshown that the VC-SGC MC technique allows ( i ) sta-bilizing arbitrary concentrations inside (and outside) themiscibility gap and ( ii ) simultaneously computing deriva-tives of the free energy with respect to concentration.Efficient sampling of the entire concentration range isstraightforward since in essence it merely requires settingthe variance parameter ¯ κ , which can be chosen rather de-liberately.Since VC-SGC MC simulations provide free energyderivatives as a continuous function of composition, ther-modynamic integration can be carried out along the con-centration axis, which allows determination of interfacefree energy densities. Specifically by monitoring the freeenergy derivative as precipitate size and shape evolve, itis possible to extract the excess free energies of compactprecipitates. By pushing the system close to 50-50 phasecomposition while varying simulation cell shape and size,one can also derive the interface free energies of flat in-terfaces as a function of orientation.To demonstrate the applicability of the VC-SGC ap-proach in more general cases, we finally considered Fe–Cralloys as described by an empirical potential scheme. Theinterface free energies for this system were determinedtaking into account not only the configurational but alsovibrational and structural degrees of freedom. Our re-sults indicate a strong temperature but weak orientationdependence of the interface free energies, suggesting thatboth α and α ′ precipitates should adopt a spherical shapeat practically all relevant temperatures.As described in Ref. 16, the VC-SGC MC ensembleleans itself to efficient parallelization, enabling simula-tions of precipitation in multi-million atom samples. Itshould finally be stressed that our approach is also di-rectly applicable for obtaining interface free energies and nucleation barriers from cluster-expansions, whichare widely applied to study complex multi-componentalloys.
ACKNOWLEDGMENTS
We would like to thank G. Gilmer at LLNL and M.Ath`enes at CEA-Saclay for helpful discussions. LawrenceLivermore National Laboratory is operated by LawrenceLivermore National Security, LLC, for the U.S. DOE-NNSA under Contract DE-AC52-07NA27344. P.E. hasbeen partly funded by the Swedish Research Coun-cil. Computer time allocations by NERSC at LawrenceBerkeley National Laboratory and by the Swedish Na-tional Infractstructure for Computing are gratefully ac-knowledged. ∗ [email protected] † [email protected] Q. Bronchart, Y. Le Bouar, and A. Finel, Phys. Rev. Lett. , 015702 (2008) J. Q. Broughton and G. H. Gilmer, J. Chem. Phys. ,5759 (1986) R. L. Davidchack and B. B. Laird, J. Chem. Phys. ,7651 (2003) J. J. Hoyt, M. Asta, and A. Karma, Phys. Rev. Lett. ,5530 (2001) Y. Mishin, M. Asta, and J. Li, Acta Mater. , 1117 (2010) S. Angioletti-Uberti, M. Ceriotti, P. D. Lee, and M. W.Finnis, Phys. Rev. B , 125416 (2010) J. H. Hetherington, J. Low Temperature Phys. , 145(1987) R. S. Johal, A. Planes, and E. Vives, Phys. Rev. E ,056113 (2003) G. M. Torrie and J. P. Valleau, J. Comp. Phys. , 187(1977) J. K¨astner and W. Thiel, J. Chem. Phys. , 144104(2005) A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 2635 (1988) A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 1195 (1989) B. Roux, Comp. Phys. Comm. , 275 (Sep. 1995) B. M. Dickson, F. Legoll, T. Leli´evre, G. Stoltz, andP. Fleurat-Lessard, J. Phys. Chem. B , 5823 (2010) G. Adjanor, M. Ath`enes, and J. M. Rodgers, J. Chem.Phys. , 044127 (2011) B. Sadigh, P. Erhart, A. Stukowski, A. Caro, E. Martinez, and L. Zepeda-Ruiz, Phys. Rev. B , 184203 (2012) A. Stukowski, Model. Simul. Mater. Sci. Eng. , 015012(2010) J. M. Sanchez, F. Ducastelle, and D. Gratias, Physica ,334 (1984) M. Asta, V. Ozolins, and C. Woodward, JOM , 16 (2001) In general, the two phases that are separated by the in-terface can exhibit size mismatch leading to an additionalstrain contribution to the interface (free) energy.
Thiscontribution can be added to the present approach withoutfurther complications. In the case of the Fe–Cr system thesize mismatch is miniscule and the strain contribution isneglible. A. Caro, D. A. Crowson, and M. Caro, Phys. Rev. Lett. , 075702 (2005) A. Stukowski, B. Sadigh, P. Erhart, and A. Caro, Model.Simul. Mater. Sci. Eng. , 075005 (2009) G. Bonny, P. Erhart, A. Caro, R. C. Pasianot, L. Malerba,and M. Caro, Model. Simul. Mater. Sci. Eng. , 025006(2009) P. Erhart, A. Caro, M. Serrano de Caro, and B. Sadigh,Phys. Rev. B , 134206 (2008) Volume change trial moves were still included in order toallow the system to relax lattice strain. P. Erhart, J. Marian, and B. Sadigh, “Equilibrium struc-ture, shape and orientation relations of BCC and 9R Cu-precipitates in Fe from atomistic simulations,” (2012), tobe submitted Y. Mishin, Acta Mater.52