Consistent Integration of Experimental and Ab Initio Data into Effective

Subscriber access provided by UNIVERSITY OF THE SUNSHINE COAST

Article

Consistent integration of experimental and ab initio data into effective physical models Lukas Vlcek, Rama K. Vasudevan, Stephen Jesse, and Sergei V. Kalinin J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00114 • Publication Date (Web): 11 Sep 2017 Downloaded from http://pubs.acs.org on September 12, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Consistent integration of experimental and ab initio data into effective physical models Lukas Vlcek*,1 Rama K. Vasudevan,2 Stephen Jesse,2 and Sergei V. Kalinin2 1

Joint Institute for Computational Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6173, United States

2

Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6496, United States

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paidup, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Corresponding Author *E-mail address: [email protected]

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 48

Abstract: We describe and test theoretical principles for consistent integration of experimental and ab initio data from diverse sources into a single statistical mechanical model. The approach is based on the recently introduced concept of statistical distance between partition functions, uses a simple vector algebra formalism to describe measurement outcomes and coarse-graining operations, and takes advantage of thermodynamic perturbation expressions for fast exploration of the model parameter space. The methodology is demonstrated on a combination of thermodynamic, structural, spectroscopic, and imaging pseudo-experimental data along with ab initio-type trajectories, which are incorporated into models describing the behavior of a nearcritical fluid, liquid water, thin-film mixed oxides, and binary alloys. We evaluate how different target data constrain the model parameters and how the uncertainty associated with incomplete target information and limited sampling of the system’s phase space might influence the choice of optimal parameters.

Keywords: statistical distance, optimization, force field, coarse-graining, measurement


2

Page 3 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction One of the most critical tasks in constructing a model of a physical system is the objective evaluation of its quality by comparing its predictions with given target data. The optimization of model parameters is then achieved by minimizing a loss function embodying the appropriate quality metric. In the case of statistical mechanical models, such as atomistic and coarse-grained force fields, the optimization is performed by modifying particle interactions so that the model system matches the target data, which may include microscopic information from ab initio calculations as well as the results of macroscopic experiments. A number of approaches have been proposed and loss functions constructed to handle diverse sources of information, such as force matching for ab initio trajectories,1-3 relative entropy minimization for coarse-grained models,4,5 the inverse Boltzmann method for pair distribution functions,6,7 or the least squares and chi-squared metrics for thermodynamic and dynamic properties.8-10 Several approaches have been designed to combine diverse sources of information in loss functions targeting simultaneously several different properties that provide joint constraints on the model.11-16 Most notably, these include multicriteria optimization using the Pareto approach by Stobener et al.,14 and a variational approach combining structural and thermodynamic data into coarse-grained models by Dunn and Noid.15 The main limitations of the currently available methods stem from the fact that they are not generally applicable to all force field functional forms (e.g., the energy and force matching techniques can diverge for hard body potentials) or types of reference data (e.g., the Boltzmann inversion principles are not directly applicable to thermodynamic data). Moreover, a


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 48

combination of several optimality criteria for different target properties into a single optimization loss function may be inconsistent, while the relative weights assigned to individual terms of the total loss function may lack theoretical justification. A typical example of such inconsistency is the use of the force matching loss function, which is asymmetric with respect to sampling from model and target system trajectories, in combination with symmetric least squares expressions for macroscopic properties.1,17 To overcome these issues and derive a well-justified form of the loss function from fundamental statistical mechanical concepts, we have recently proposed an alternative approach based on the minimization of a metric known in physics as statistical distance,17 which had been used to quantify distinguishability between two sets of experimental outcomes.18 The approach bears similarity to minimization of the difference between partition functions proposed by Stillinger in early 1970s,19,20 except here we deal with the normalized probability distributions of system states. At the highest resolution level, the loss function quantifies the differences between the microstate probability distributions of the target and model systems and, through systematic coarse-graining, can arrive at the appropriate expressions for macroscopic structural and thermodynamic data.17 We demonstrated that the resulting loss functions are robust, applicable to any type of model functional forms, and can incorporate any kind of observable target data in the same unified way. More recently, we have proposed a simple way to include information from interatomic forces (and higher energy derivatives) obtained from ab initio MD simulations to this optimization loss function, increasing thus significantly the amount of target data available for force field development.21 Moreover, by matching overall microstate probability distributions, rather than select sharply defined properties, the minimization of statistical distance naturally optimizes entropy and fluctuation properties of all orders, and therefore enhances the


4

Page 5 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


transferability of the model over a range of thermodynamic conditions. One important question left open in our previous study was how to consistently and systematically combine different sources of information in a single loss function without resorting to arbitrary weighting of individual contributions. In this article we will introduce and discuss a theoretical framework that enables seamless integration of data from arbitrary sources into a model representation of a target system. For that purpose in Section 2 we review the theoretical background and formalism underlying the model optimization principles, and in Section 3 describe our parameter space exploration technique and the overall optimization procedure. Then, in Section 4, we demonstrate the methodology on examples of bulk and interfacial systems involving the use of different types of target data. Finally, we conclude with a review of the basic principles underlying the proposed approach, discuss its role in general problems of data integration, and outline its current limitations and future directions. 2. Theoretical background 2.1. Basic relations We will restrict the following discussion to physical systems that can be described by semiclassical statistical mechanics, in which the quantum nature of real systems is only reflected in the fundamental resolution limit given by the Heisenberg uncertainty principle. All measurable properties A of such system P are fully determined by the probability distribution over its k distinguishable microstates as, k

A = Ai

= ∑ Ai pi , p

(1)

i=1

where

denotes the ensemble average over all states of system P, Ai is the contribution to


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 48

property A from microstate i, and pi is the probability that a random sample (measurement) will find the system in this i-microstate. In the special case of equilibrium thermodynamic systems, the microstate probabilities follow the Boltzmann statistics according to pi =

exp −β u p,i  Zp

,

(2)

where up,i is the energy corresponding to microstate i, β = 1 kBT , kB is the Boltzmann constant, and T is absolute temperature. The normalizing constant Zp is the partition function defined as k

Z p = ∑ exp −β u p,i 

,

(3)

i=1

with the summation running over all k microstates. In our previous study,17 we showed the convenience of representing k distinguishable microstates of a system as state vectors forming the orthonormal basis of a k-dimensional realvalued Hilbert space H Sk . A particular microstate probability distribution P can be then represented as a unit vector in the positive orthant of this space (Figure 1).18,22


6

Page 7 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. The state space H S3 of a classical three-state system. Red arrows: unit vectors representing individual states; green arrow: a unit vector determining the microstate probability distribution of system P with coordinates given as the square roots of microstate probabilities.

2.2. Representation of measurement outcomes While a statistical mechanical system is defined by its microstate probabilities, in practice any of its properties can only be determined by measurements that are always limited by resolution and finite-size sampling. To properly account for these inherent characteristics of target data, we will briefly discuss the mathematical representation of a measurement. We will assume that it is deterministic, i.e., repeated sampling of the same microstate will result in the same outcome, and that there are no experiment or simulation-specific errors. Let us first consider a measurement consisting of n samples that probe a k-state system at the microstate resolution. Each sample of this system obtained at a distinct time t can be represented as a projection from the system’s Hilbert space H Sk to the k-dimensional Hilbert space of measurement outcomes, H tk , with each possible distinguishable outcome forming an orthonormal basis vector. The outcome of the whole measurement is then a vector summing the n individual sample vectors in a space constructed as the direct sum of n H tk spaces (Figure 2a), n

H k,n = ⊕ H tk

.

(4)

t=1

Each dimension of the combined space can be indexed by a pair of integers: a ‘time’ coordinate t identifying a specific sample, and a ‘space’ coordinate i identifying a detector (histogram bin). As the system is always in one of its states, there is exactly one of the orthogonal vectors in each subspace H tk set to unity, with the others equal to zero. The total sum of n orthogonal unit


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

vectors is thus a vector of length

Page 8 of 48

n uniquely identifying a particular outcome of a

measurement.

Figure 2. (a) A schematic representation of a four-sample measurement performed on a threestate system. Red arrows: unit vectors corresponding to microstates detected at times t1 through t4. (b). A series of two coarse-grained measurements e1 and e2, each consisting of two samples. Red arrows: coarse-grained states detected in the two experiments. (c) Representation of twosample measurements e1 and e2 without distinguishing the sample order.


8

Page 9 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3. Coarse-graining in space and time For measurements at lower resolutions, two or more outcomes become indistinguishable and will be registered as a coarse-grained result. In our formalism, such a measurement will correspond to a projection from H k,n to a lower-dimensional Hilbert space (Figure 2b,c). This operation preserves the total length of the vector summing all microstate vectors from the subspace, but loses information about its orientation in the subspace. In general, coarse-graining applies whenever we are unable (or unwilling) to distinguish between different measurement outcomes. One type of coarse-graining we may consider is merging ‘spatial’ dimensions distinguishing system states at a given time. The outcome of such measurement e distinguishing ke coarsegrained states at time t can be represented in Hilbert space H e,tke (Figure 2b). These may include coarse-graining resulting from limited experimental resolution with respect the closeness of microstate energies, momenta, volumes, spatial coordinates, etc. A special subtype of such coarse-graining, which is important in computational applications, is merging microstates along selected degrees of freedom, effectively thus eliminating these degrees and potentially speeding up a computer simulation. While such coarse-graining can be treated in exactly the same way within the proposed framework, there are additional considerations related to the choice of the degrees of freedom to be eliminated and other technical issues, which are beyond the scope of this presentation, and will be dealt with elsewhere. Another type of coarse-graining can result from losing time information about the order of samples. Even though the H k,n representation applies to any kind of dynamic system, we often deal with systems whose microstate probability distributions do not evolve in time. In this case we may not be interested in the exact order of samples and simply record the total counts over a period of time. This loss of time information corresponds to coarsening over n time coordinates


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 48

for each ‘space’ coordinate i, into a k-dimensional space H eke , as illustrated in Figure 2c. If the samples can be considered statistically independent, the probabilities of different outcomes of such static experiments are described by the multinomial distribution23 k k piyi (5) MD ( y1 yk , n, p1 pk ) = n!∏ ; ∑ yi = n ; ∑ pi = 1 , i=1 yi ! i i where k is the total number of histogram bins, yi is the number of counts in bin i in a k

measurement, pi is the probability of an individual count being registered in bin i. The relative frequency of counts yi n is the maximum likelihood estimate of pi, and as n → ∞ , yi n k converges to pi. The endpoints of all vectors in H representing the results of n-sample

measurements without time information lie on the surface of a sphere with radius

n , and the

vectors corresponding to relative frequencies lie on a sphere with a unit radius. Every point in the positive orthant lying on this k-1 dimensional spherical surface then represents a probability distribution with coordinates given as

pi (Figure 1).

2.4. Statistical distance as a measure of system similarity As shown by Rao,24,25 the probability space of multinomial distributions, Eq. (5), can be equipped with a natural Riemannian metric equal to the Fisher information matrix, with distances between pairs of points representing probability distributions measured by Rao’s distance, whose square can be written as

DR2 = 4ns 2 .

(6)

Here s is statistical distance defined as a geodesic on the spherical probability surface (Figure 3a),26,27


10

Page 11 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


 k  s = arccos ∑ pi qi  ,  i=1 

(7)

where qi is the probability of microstate i in system Q. Distance s was introduced to physics by Wootters,18 who recognized its usefulness as a measure of distinguishability between outcomes of physical measurements (sometimes s is called Wootter’s distance), and further studied by Braunstein and Caves,22 who showed its relation to classical and quantum relations of uncertainty. A critically important property, which distinguishes it from other metrics, such as relative entropy or chi-squared distance, and makes it uniquely suited for the task of model development, i.e., approximation rather than estimation, is the built-in variance equalizing transformation, which ensures that all data have the same statistical weight,17,28 and can be thus easily combined and compared. Another essential property, which it shares with most sensible loss functions, is the invariance of its optimum with respect to the number of independent subsystems.17 For instance, the interaction parameters of a model of a single harmonic oscillator should be identical to those of a system of N such oscillators (e.g., Einstein crystal). The argument in Eq. (7), known as the Bhattacharyya coefficient, CB, is a scalar product between the probability distribution vectors P and Q (Figure 3b). In the special case of an equilibrium canonical system following the Boltzmann statistics, CB attains a special form, k

CB = ∑ pi qi = exp − β ( ∆ui − ∆F ) 2 = exp + β ( ∆ui − ∆F ) 2 i=1

q

p

(8)

where ∆ui = uiP − uiQ is the difference between the energies of microstate i in systems P and Q, while ∆F = −kBT ln ( Z p Zq ) is the Helmholtz free energy difference between the two systems. We note that the probabilities pi used in practical calculations will usually be derived as maximum likelihood estimates based on relative frequencies, yi n , obtained from the statistical


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 48

analysis of available data. The finite resolution of most experimental data limits our ability to distinguish individual microstates. Following the geometric representation, Figure 2 and Eq. (4), the results of such measurements will be represented as vectors in ke-dimensional Hilbert space H eke obtained by coarsening projections from H k . Since projected distances are always equal or less than the original ones, minimization of statistical distances between the microstate distributions of two systems reliably leads to the convergence of all their coarse-grained experimental results. We have used these convenient properties as a basis for rigorous force field development methodology with well-defined relations between system properties at different resolutions levels.17

Figure 3. (a) The probability space of trinomial distributions with a statistical distance s between distributions (i.e., statistical mechanical systems) P and Q. (b) The relations between statistical distance s, the Bhattacharyya coefficient CB, and the Hellinger distance DH on the example of binomial distributions.


12

Page 13 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.5. Combining and comparing series of measurements When dealing with a single source of reference data, the total number n of samples is a factor in Eq. (6) that does not affect the location of the minimum of DR and can be therefore omitted from the loss function. In our previous study,17 we discussed several metrics with monotone relations to statistical distance, such as the Bhattacharyya and Hellinger distances, leading to identical optima. However, when multiple experiments with different amounts of sampling and different resolutions are combined, their respective statistical uncertainties depending on the number of samples n need to be taken into account. To derive the appropriate form and weighting factors for a loss function combining diverse target data, we will return to the original description of the full space H k,n as a direct sum of single-sample spaces H tk . A series of E experiments can be thought of as composed of the blocks of ne individual samples, with each block presenting different projections of the full space into a lower-dimensional space H eke ,ne , defined by the resolution and sampling of experiment e (Figure 2b). The resulting space of the series of experiments is then a direct sum of the individual subspaces E  ne  E H k,n = ⊕  ⊕ H e,tke e  = ⊕ H eke ,ne . e=1  te =1  e=1

(9)

In case the individual experiments do not distinguish the order of samples (Figure 2c), coarsegraining in time leads to the space with reduced dimensionality defined as, E

H k,n = ⊕ H eke .

(10)

e=1

Consequently, for the sets of such space and time coarse-grained measurements, model optimization is achieved by minimizing Rao’s distance between vectors defined in space with E

G = ∑ ke dimensions according to, e=1


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 48

E

D = 4∑ ne se2 , 2 R

(11)

e=1

where se is statistical distance on the space of experiment e. Rather than working with absolute sample numbers, which may not be readily available, it may be more convenient to use relative statistical weights of each experiment and define a normalized statistical distance between a target system and its model based on E experiments as, E

ne 2 se . e=1 n

S =∑ 2 E

(12)

2.6. Data-specific forms of statistical distance To use Eq. (12) as a loss function for practical applications, it is necessary to express the abstract distances se2 in forms tailored for given target data. This is achieved by expressing state probabilities pi and qi in Eq. (7) using the appropriate relations, such as the Boltzmann relation, Eq. (2), for equilibrium systems, normalized histograms for various statistics of structural features, or by values implied by thermodynamic fluctuation relations. We have shown before how structural and thermodynamic data can be used within the statistical distance framework,17 but for completeness here we summarize expressions for data appearing in the examples given in Section 4. For the sake of convenience, rather than se2 we will present expressions for the experiment-specific Bhattacharyya coefficients CB,e , related straightforwardly to the loss function as,

se2 = arccos2 (CB,e ) .

(13)

Similarity between two ab initio-type trajectories, i.e., sequences of particle configurations and corresponding configurational energies, which may originate from quantum chemical or


14

Page 15 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


classical simulations, can be quantified using the following estimator based on Eq. (8), M

CB,AI

1 p = ∑ exp− β ( ∆ui − ∆F ) 2 M p i=1

(14)

where MP is the number of samples generated in an equilibrium simulation of the target system P. A simple extension of Eq. (14) incorporating atomic forces has been presented recently.21

Similarity between atom-atom pair correlations in the target and model systems can be evaluated by measuring the statistical distance between histograms of pair distances collected from the two systems. Denoting the relative frequencies of counts within histogram bins l as hp,l and hq,l , for systems P and Q, respectively, we can write the Bhattacharya coefficient as k

CB,g = ∑ hp,l hq,l ,

(15)

l=1 k

with normalization condition

∑h l=1

k

p,l

= ∑ hq,l = 1. The same expression can be used for histograms l=1

of arbitrary properties, such as the statistics of local configurations collected from microscopic 2 , or other structural descriptors discussed in Section 4.2. images, sIM

While the measurement of thermodynamic and mechanical properties is not typically carried out as series of individual samples, we can still use statistical mechanical relations between the observed properties and the underlying microstate probability distributions, and interpret the data in terms of sampling. To address the problems studied in Section 4, we include the examples of matching (pseudo-) experimental densities, compressibilities, and pressures. Using the fluctuation relations,29 we can reconstruct particle number distributions of a fluid in volume V at a given temperature T and chemical potential µ from known average densities and compressibilities. The corresponding Bhattacharyya coefficient CB,µ for two Gaussian


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 48

distributions approximating the particle number statistics in target P and model Q, can be written as,

CB,µ

Here

N

2  N − N  1  p q = exp − 2 2 2 2 ( ∆N ) p + ( ∆N ) q  4 ( ∆N ) p + ( ∆N ) q  . 2   2 ρ p ρq χ T , p χ T ,q ( ρ p − ρq )  − V = exp  4kBT ρ p2 χ T , p + ρq2 χ T ,q  ρ p2 χ T , p + ρq2 χ T ,q  

(

2 ( ∆N ) p ( ∆N ) q

)

(16)

denotes the average number of particles in the simulation box, with the

p

corresponding variance calculated as

( ∆N ) p = 2

N2

2

p

− N p . The Gaussian shape of the

distributions is derived from the fluctuation relations of statistical mechanics29 with mean and variance related to density ρ p and isothermal compressibility χ T , p as ρ p = N

( ∆N ) p 2

N

p

p

V and

= kB T ρ p χ T , p .

Similarly, instantaneous pressure distributions of the target and model systems can be modeled as the Gaussian functions. In case we have no a priori knowledge of the target system compressibility, such as in the example in Section 4.1.2, we can get around such lack of target information by approximating the unknown compressibility value by that of the model itself. Consequently, we can minimize the distance between two Gaussian distributions with the same variance, for which the pressure-based Bhattacharyya coefficient assumes the form

CB,P

2  Pp − Pq )  ( 1  = exp −  4 2 ( ∆P ) 2   q 

(17)

where Pp and Pq are the average pressures of the target and model systems, respectively, and

( ∆P )q is the variance of the model system instantaneous pressure distribution. 2


16

Page 17 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In practical applications the total sample numbers and related statistical weights of individual experiments may not be easily known, and must be therefore estimated. In some cases we can assume that the statistical uncertainty of an experiment is very small, which will make the particular contribution to the overall loss function a nearly rigid constraint on model optimization. Additional experiment-specific uncertainties, such as noise in real data, will typically soften such constraints. To avoid these additional influences, which might obscure the inherent sampling uncertainty of statistical mechanical systems, in this study we use pseudoexperimental or pseudo-ab initio data generated in simulations, and postpone the discussion of these important, but extraneous issues.

3. Technical details 3.1. Exploring model parameter spaces via thermodynamic perturbation relations A search for the optimal set of model parameters can be efficiently accomplished using the perturbation technique outlined and applied in our earlier publications.17,30-32 The approach uses the thermodynamic perturbation relations of statistical mechanics to predict changes of arbitrary properties A of system P as a function of interaction parameters. This is achieved by sampling a reference system Q and reweighting its contributions Ai according to:

Ai

k k  pi  = A p = A   qi = Ai exp −β ( ∆ui − ∆F ) ∑ ∑ i i i p  qi  i=1 i=1

q

,

(18)

where A can either depend explicitly on the interaction parameters, as in the case of configurational energy, or can be independent of interaction parameters, as for volume and structural properties. In the case of structural properties, such as pair distribution functions, each quantity Ai is an array containing contributions to the pair distance histograms from configuration i.


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 48

The sampling of the reference system can be obtained either from an ab initio trajectory of the target system or from simulations using a preliminary model with suboptimal parameters. The predicted properties can be then compared with the target data using the loss function introduced in the previous sections. The accuracy of predictions for the property A based on Eq. (18) depends both on the overlap of the two systems’ probability distributions and on sufficient sampling of the reference system configurational space. Since the best models are characterized by large overlaps with the target microstate distribution, the perturbation relations using the target system trajectories as a reference can be expected to produce very accurate predictions. A combination of several reference simulations using the multistate Bennett’s acceptance ratio method (MBAR)33 is possible, but we have found that simple importance sampling technique, Eq. (18), is sufficiently accurate and efficient. The key to avoiding the need for repeated reference system simulations for each new parameter combination is to express the microstate energy difference, ∆ui , in a form that separates the model parameters from quantities defined solely by the reference system configurations (examples given in the following paragraph). After the quantities associated with the reference system are collected from a simulation, they can be combined with arbitrary sets of model parameters at the comparably negligible cost of a few array operations rather than that of a full new simulation. Subsequent reweighting according to Eq. (18) can be then performed to predict property A of the model and the corresponding statistical distance from the target data. To illustrate the theoretical framework and methodology introduced in this study, we choose model systems with common forms of interaction potential functions. In Section 4.1, we will optimize molecular models with van der Waals interactions defined by the Lennard-Jones (LJ) potential and electrostatic interactions defined by point charges,


18

Page 19 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


 12  6  σ σ 1 q j qk A C E u ( rjk ) = 4ε   −    + = 12 − 6 +  rjk   rjk   4πε0 rjk rjk rjk rjk   where

(19)

σ and ε are the size and energy parameters of the LJ potential, qj is the magnitude of a

point charge on interaction site j, and rjk is the distance between atoms j and k. The configurational energy of the system in microstate i is then given as,

ui ( A,C, E ) = AaiA − CaiC + EaiE where

(20)

aiA = ∑∑1 rjk12 ; aiC = ∑ ∑1 rjk6 ; aiE = ∑∑1 rjk j k> j

j k> j

j k> j

Here numbers aiA , aiC , and aiE are collected from the reference system simulation for all pairs of particle types in each configuration i, and subsequently multiplied by the desired model parameters A, C, and E. The energy difference ∆ui appearing in Eq. (18) is then obtained by subtracting the known reference system energy uiR from the predicted model energy given by Eq. (20) as

∆ui = ui ( A, C, E ) − uiR

(21)

In section 4.2 we will deal with coarse-grained models of crystalline materials approximated by the Ising-type lattice potentials with pairwise-additive interactions between particles. The configurational energy of a two-component system described by this potential is given as, 2

2

ui = ∑∑ wXY aiXY ,

(22)

X=1 Y =1

where X and Y denote the types of the two components, and the quantities aiXY are collected from the simulations of the reference system for each (X,Y) combination as,

aiXY = ∑δ XY ,

(23)

NN


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 48

with the summations running through all nearest neighbor (NN) pairs of lattice sites; the quantity

δ XY is equal to either one, if the particles are of types X and Y, or zero otherwise. The models used in Sections 4.2 will combine several of the terms defined in Eq. (22) to describe interactions between sites in different geometrical arrangements (e.g., next-nearest pairs).

3.2. Optimization procedure 3.2.1. Choice of target data The choice of a particular target dataset depends on the purpose of the model and the demands of its future user. The case may be that the model is supposed to reproduce all available data, or that its user wishes to apply it in a specific range of conditions and therefore only certain data may be used as a target. Since the relations introduced in this study only prescribe how the data should be used consistently to produce the best possible model, the task of data selection ultimately lies outside the scope of the current work. As a general guideline though, the different sources of target data should be minimally correlated. It is also advisable that ab initio data with small number of samples nAI be supplemented by multiple accurate experimental constraints with large ne as illustrated in Section 4.1. In this case it is critical that the microscopic and macroscopic data are combined consistently using Eq. (12) so that they do not drive the model parameters to incompatible optima. It is clear from the direct link of a model with target data that its parameterization will depend on the conditions under which the target system was measured. However, even though the optimal model will depend on such conditions (e.g., temperature), a properly designed model will also be transferable outside the range of the original conditions, a quality that naturally follows from the use of statistical distance.17


20

Page 21 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2.2. Assigning statistical weights (sample numbers) As with the choice of target data, the assignment of relative sample numbers ne/n in Eq. (12) is dictated by the purpose of the model. If we want the model to match an outcome of a concrete set of actual measurements, the sample numbers will be equal to those used in the target experiment. Most of the time, however, a model is designed for general use by different researchers without a single particular experimental setup in mind. In this case we can interpret the relative sample number ne/n as the probability that a model will be used to predict the outcome of experiment e. For instance, models are often designed for simulations of fluids and phase equilibria over a certain range of conditions. We can then use target thermodynamic data measured over a range of conditions and assign weights according to the expected use, such as emphasize region near the critical point or ambient conditions, or we may choose to set all weights equal. Ideally, all models should be accompanied by a ‘certificate’ that would list the target data, their relative weights, and the statistical distance SE to this dataset. 3.2.3. Design of the preliminary model The main goal when designing a preliminary model for generating a reference system trajectory to be used in the perturbation technique (Section 3.1) is for it to have a reasonable overlap with the target system microstate probability distribution. While the final result of the optimization does not depend on this choice, a good starting model will speed up the procedure. This task is heuristic and depends on our prior knowledge about particle interactions within the target system. In case target system ab initio trajectory is available, it can be used for the reference system input in place of the preliminary model. In this case maximum overlap is guaranteed and model parameters can be optimized using Eq. (14). If other researchers have earlier developed a model of the target system, it may save effort to


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 48

build on their work and use such a model for reference simulations. If none of the above is available, the preliminary model parameters can be estimated by analyzing the target data. For instance, the statistics of the nearest neighbor pairs on a lattice can be used for a rough estimate of the nearest-neighbor interactions using the Boltzmann relation w = −kBT ln ( pAA pAB ) , where pAA and pAB are the probabilities of finding the nearest neighbor pairs of like and unlike particles,

respectively, and w is the difference between interactions of such pairs. 3.2.4. Perturbation search for optimal parameters The optimization algorithm can be summarized as follows: (i)

Collect atomic coordinates and configurational energies along with any structural or thermodynamic property of interest from a reference simulation.

(ii)

Use atomic coordinates to calculate quantities aiX and bi,Xj for each configuration i according to Eqs. (20) or (23).

(iii)

Multiply the arrays of aiX and bi,Xj by trial force field parameters to determine the model configurational energies for each configuration i according to Eqs. (20) or (22).

(iv)

Use these model energies along with the known configurational energies of the target system to calculate the configurational energy differences, Eq. (21).

(v)

Rescale the reference system property of interest to obtain the corresponding model property using Eq. (18).

(vi)

Using the rescaled properties and target system properties, calculate the loss function values using Eqs. (15)-(17). If target system ab initio-type trajectory is used as a reference system trajectory, use Eq. (14) to calculate the loss function.

(vii)

Combine all partial loss functions weighted by the target data statistical weights (sample numbers) using Eq. (12) to determine the overall loss function SE2 .


22

Page 23 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(viii) If the minimum of SE2 is reached, output the optimal parameters and end the procedure. Otherwise, modify the force field parameters using a search algorithm, such as the conjugate gradient method,34 and return to (iii). The resulting model with optimized parameters can be used as a new reference system and, if desired, additional refinement can be attempted by repeating the whole procedure starting from (i). The perturbation search is currently implemented as a set of Python scripts using SciPy libraries for the conjugate gradient minimization.

4. Results and discussion To illustrate the versatility of the principles introduced in Section 2 and to demonstrate the accuracy and efficiency of the perturbation technique described in Section 3, we investigated a diverse set of modeling problems which reflect our interests in the development of accurate atomistic force fields, as well as in the methodology of data integration and design of experiments. These include the development of effective force fields based on the combination of microscopic ab initio-type data with macroscopic thermodynamic and structural constraints in Section 4.1, and the design of simple coarse-grained models of solid materials based on structural information from imaging techniques in Section 4.2. Rather than using real experimental and ab initio input, we generated pseudo-experimental and pseudo-ab initio data using well-defined interaction potentials for the target system to eliminate experiment-specific errors, allowing us thus unbiased evaluation of the procedure and testing of model identifiability.

4.1. Accurate atomistic force fields 4.1.1. The equation of state and compressibility of a supercritical fluid


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 48

Target and model systems

In the first example we consider supercritical argon over a range of chemical potentials, with two sources of target data available: (i) an ab initio–type molecular dynamics (MD) trajectory providing a series of energy-configuration pairs, and (ii) a set of fluid densities and compressibilities at m = 20 different chemical potentials obtained from target system simulations. The target system interactions were represented by the Buckingham potential,35 uij ( rij ) = A exp −Brij  − C rij6

(24)

with the values of parameters A, B, and C given in the caption of Table 1. The interactions of the model system were represented by the Lennard-Jones potential, Eq. (19), whose

σ and ε

parameters were to be optimized. All target and model simulations were performed at the constant supercritical temperature of 160 K in a cubic simulation box with side L=24 Å and the interaction cut-off distance of 12 Å. The ab initio-type trajectory was generated in canonical MD simulations36 with N = 123 over 100 ps with 1 fs time steps, and provided nAI = 100 independent samples. The average densities and compressibilities were generated in grandcanonical Monte Carlo (GCMC) simulations37 carried out at selected chemical potentials (see Figure 4) over 5.106 individual translational and insertion-deletion moves, providing nµ = 5000 samples, which were used to calculate fluid densities and compressibilities. Results 2 Following the procedure described in Section 3.2, we first minimized statistical distance sAI

between the model and target ab initio trajectories using Eq. (14). The fast parameter search (~2 seconds on a single computer node) resulted in the parameters listed in Table 1. The comparison


24

Page 25 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


between the target and model-predicted densities and compressibilities over a range of chemical potentials in Figure 4 shows nearly quantitative agreement between the two data sets, with noticeable discrepancies seen only in the sensitive area near the critical region and for compressibility at the highest densities. For illustration purposes, the statistical distance landscape as a function of the model parameters is displayed in Figure 5a.

Figure 4. Argon density (a) and compressibility (b) at 160 K as obtained for the target system (black) and different models defined in Table 1.

2 For subsequent optimizations we used the sAI -derived model as a reference system. Using

reference data generated in GCMC simulations, we optimized the LJ parameters by minimizing statistical distance sµ2 , Eq. (16), based on density and compressibility values for each considered chemical potential µ .Model parameters resulting from simultaneous minimization of the set of m distances sµ2 for all chemical potentials are listed in Table 1, and the landscape of the

combined distance s

2 EOS

1 m 2 = ∑ sµ as a function of model parameters is shown in Figure 5b. m µ =1


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 48

The total loss function for the given set of ab initio and thermodynamic data was then obtained as a linear combination of individual squared distances weighted by the amount of sampling, 2 2 SE2 = ( nAI sAI + mnµ sEOS ) ( nAI + mnµ )

.

(25)

2 2 As evident from Table 1, all sAI , sEOS , and SE2 predict very similar optimal parameters,

highlighting the internal consistency of the present approach. Due to the much larger statistical weight resulting from better sampling of the thermodynamic properties, the total SE2 nearly 2 coincides with sEOS . The parameters are also nearly identical to those found in our earlier

study,17 in which the microstate probability distribution was probed by 106 samples.

Figure 5. The landscapes of statistical distances in the parameter spaces of the argon LJ parameters

σ and ε based on (a) ab initio data, sAI2 , and (b) on thermodynamic data from fluid

2 densities and compressibilities, sEOS . The set of optimal parameters is indicated by a white cross.

(The white regions in (b) are areas of high statistical uncertainty resulting from imperfect sampling.)


26

Page 27 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For comparison, we included in Table 1 the LJ parameters obtained earlier using the basic form of the force matching technique.17 The equation of state predicted by the force matchingbased model plotted in Figure 4, clearly shows the limited accuracy over the range of thermodynamic properties.

Table 1. The optimized Lennard-Jones parameters based on different loss functions. The target system parameters defined in Eq. (24) were set to A = 1017742.0 kJ/mol, B = 3.663 Å-1, and C = 6142.585 Å6 kJ/mol.35 Loss function

σ [Å]

ε [kJ/mol]

2 sAI

3.395

1.100241

2 sEOS

3.3951

1.0821

SE2

3.3951

1.0821

DB [Ref 17]

3.39574

1.100324

DFM [Ref 17]

3.41823

0.938115

4.1.2. Liquid water Target and model systems

For the second illustration we consider liquid water at ambient conditions with target data combining (i) microscopic information from an ab initio-type trajectory, (ii) atom-atom pair distribution functions, and (iii) average pressure generated in target system simulations. The target system interactions were represented by the Gaussian charge polarizable model (GCPM) of water38 characterized by rigid geometry, Gaussian diffuse charges located at the hydrogens and the auxiliary site M, a single polarizable site located at the molecular center of


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 48

mass, and a modified Buckingham potential, Eq. (24), interaction between oxygen sites.39 For the model system we chose the TIP4P-type rigid geometry with three atomic sites and one massless auxiliary site M located on the H-O-H bisector at a distance lOM from the oxygen site; we set lOM = 0.1546 Å, used for the TIP4P/2005 model.40 The interaction potential between individual sites has the general form of Eq. (19), and we aimed optimized parameters

σ and ε

of the LJ potential defining the van der Waals interactions between oxygen atoms, along with a parameter defining the magnitude of the point charge on the hydrogen site H, qH. All data were generated in canonical MD simulations at constant temperature T=298 K and density ρ = 0.998 g/cm3. The ab initio-type trajectories of the target system were obtained from 100 ps simulations in a system with 108 GCPM water molecules, providing the total of nAI = 103 samples. Oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen pair distribution functions were obtained from 1 ns MD simulations at the same density with 256 GCPM molecules and providing ng=nP=104 samples. Model MD simulations were performed with 108 molecules over a time span of 10 ns, and provided 105 samples for the perturbation technique optimization. Results

As in the previous example, the first step in the optimization was to minimize the ab initio2 based statistical distance sAI according to Eq. (14). The resulting parameters are listed in Table

2, and the model predictions of the atom-atom pair correlation functions and instantaneous pressure distributions are shown in Figure 6. 2 In the second step we used the sAI -based model as a reference system for optimization based 2 2 2 on matching pair distribution functions according to distance sg2 = ( sgOO + sgOH + sgHH ) 3, Eq.

(15). Even though the resulting optimal parameters listed in Table 2 differ significantly from those usually assigned to the TIP4P-type models, the predicted pair distribution functions shown


28

Page 29 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Figures 6a-6c are reproduced very accurately. However, the predicted pressure (Figure 6d) deviates significantly from the target value of 100 kPa, suggesting that the information contained in the pair distribution functions is not enough to constrain the model thermodynamic properties.

Figure 6. (a)-(c): Oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen pair distribution functions for the GCPM water target system (black), and TIP4P-type models optimized based on 2 ab initio data, sAI , (red), pair distribution functions, sg2 , (blue), and the combination of ab initio,

pair distribution, and pressure data, SE2 , (green). For comparison, the predictions of the TIP4P/2005 model are included (purple). (d) Pressure distributions for different models following the same color-coding.


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 48

To correct this deficiency, we complemented the target data with the desired pressure of 100 kPa, and performed further optimization using an additional constraint in the form of pressurebased statistical distance sP2 , Eq. (17). The resulting combined statistical distance has the form, 2 2 2 2 SE2 = nAI sAI + ng ( sg,OO + sg,OH + sg,HH ) + nP sP2  ( nAI + 3ng + nP )

,

(26)

where, for the purpose of this study, we assigned the same statistical weights to the structural and thermodynamic data, i.e., nP = ng . In general, the weights nP and ng will reflect specific target data and may not be equal. Minimization of SE2 in a single perturbation iteration led to the parameters listed in Table 2, which are close to those of TIP4P/200540 and other similar models. The predicted pair distribution functions and instantaneous pressure distributions are shown in Figure 6. It can be seen that adjusting the pressure of the model led to a slight deterioration of the structural properties, manifested as the increased heights of the first peaks of gOO and gOH. Such compromises can be expected as a result of the limited flexibility of the model’s functional form. To gain better insight into the constraints brought by different types of target data, we plotted the corresponding statistical distances as functions of

σ and ε LJ parameters in Figure 7. These

plots were generated using the perturbation expressions, Eq. (18), based on data from the same reference simulation used for parameter optimization while keeping the point charges fixed to the optimal value qH = 0.5415 e, Table 2. It is notable, that the landscapes of all statistical distances shown in plots (a)-(d) exhibit similar banana-shaped valleys indicating a range of nearoptimal

σ - ε combinations. However, the landscape of sP2 based on pressure matching shown in

2 plot (c) differs significantly from those of sAI and sg2 by constraining the acceptable parameter


30

Page 31 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


values to smaller

σ and larger ε , which possibly indicates higher sensitivity of pressure to the

slope of the repulsive part of the LJ potential. The strong pressure constraints on the parameter values then heavily influence the landscape of the overall statistical distance SE2 shown in plot (d), indicating that complementing structural data with pressure information is highly beneficial.

Figure 7. The landscapes of statistical distances in the space of

σ and ε LJ parameters based on

2 (a) ab initio data, sAI , (b) pair distribution functions, sg2 , (c) pressure, sP2 , and (d) combination of

all target data, SE2 . The optimal parameters for each partial loss function are indicated by a red cross and the optimal parameters for the combined loss function are indicated by a white cross.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 48

Table 2. The Lennard-Jones parameters and hydrogen partial charges of simple non-polarizable models based on different loss functions, shown along with the TIP4P/200540 model parameters. Loss function/model

q H [ e]

σ [Å]

ε [kJ/mol]

2 sAI

0.5378

3.2159

0.6268

sg2

0.5438

4.1023

0.0272

SE2

0.5415

3.1222

0.8217

TIP4P/200540

0.5564

3.1589

0.7749

4.2. Simple coarse-grained models In this section we illustrate the proposed optimization framework on a class of problems most relevant for solid-state materials science, and our main interest is in modeling effective driving forces behind structures observed in imaging and scattering data. In ways fully analogous to those discussed in Section 4.1., we can use the statistics of various structural features to refine statistical mechanical models. However, there is an important difference when dealing with solid-state imaging data as opposed to the thermodynamic or structural properties of fluids, which stems from the limited number of independent samples (images) of the system states. One way to get around this constraint is to analyze the range of correlations between different structural features, and then treat uncorrelated local structures as independent samples, even though they are in fact a part of the same state. In the following two examples we will be dealing with crystalline systems that can be easily mapped onto a lattice, which makes it natural to approximate the systems by simple and computationally efficient Ising-type (or lattice gas) models. Our primary goal is thus finding


32

Page 33 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


effective pair-wise interactions between neighboring lattice sites that would be consistent with the observed statistical distributions of structural features.

4.2.1. Elemental segregation in mixed oxide thin films Target and model systems

In this example we aim to assess the driving forces behind metal atom segregation in thin mixed metal

oxide films,

such

as

LaxCa1-xMnO3

(LCMO) which

may exhibit

colossal

magnetoresistance.41 Such information would allow us to predict elemental segregation beneath the surface that determines the materials properties, but which would be otherwise inaccessible. The target data mimic the inputs from (i) electron microscopy images providing structural information about the surface layer, and (ii) from angle-resolved X-ray photoelectron spectroscopy (XPS) supplying concentration profiles across the interfacial layers. In this case of an LCMO-like material we only consider the distribution of two metal atom types: Ca2+ (C) and La3+ (L), whose relative concentration is kept constant at x = 0.5, while the remaining species forming the crystal lattice are considered fixed. In addition, the thin-film particles interact with the underlying substrate (S), as illustrated in Figure 8a. The target system properties were generated in simulations with a simple lattice nearest-neighbor model according to the general form of Eq. (22), i.e.,

 1 ui = wCC ∑δ CC + wCL ∑δ CL + wLL ∑δ LL  + wCS ∑δ CS + wLS ∑δ LS , 2  NN NN NN NN NN

(27).

where the individual sums run over all pairs of nearest neighbor (NN) lattice sites (sharing the face on a simple cubic lattice), while δ XY =1 if the NN particles are of types X-Y and otherwise

δ XY = 0 . The target system energy interaction parameters wCC, wCL, and wLL (Table 3) reflect the


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 48

expected repulsive interactions between metal cations, and additional interaction parameters wCS and wLS were introduced to account for metal atom interactions with the substrate. The target canonical Monte Carlo (MC) simulations were carried out in a rectangular box with two periodic dimensions LX = LY = 52 and non-periodic dimension LZ = 20 (Figure 8a) at reduced temperature kBT/ε = 0.35, where ε is the energy unit scaling interaction parameters. The total number of n=105 independent samples was collected from each target system simulation consisting of 107 individual MC steps, and used to calculate the statistics of six local surface configurations shown in Figure 8c. The same samples were also used to calculate the relative concentrations of C and L particles in individual layers. An individual MC step was performed by swapping the chemical identities of randomly chosen nearest neighbor particles. For model systems we employed the same simple nearest-neighbor interactions between particles C and L as for the target system, Eq. (27), to test model identifiability and convergence of the perturbation technique. All optimization simulations were performed in a box of dimensions LX = LY = 6 and LZ = 20 over 7x108 MC steps, and followed the same simulation protocol as for the target system.


34

Page 35 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. (a) The structure of the thin-film ‘mixed oxide’ target system. Particle C (blue), particle L (red), substrate S particles (green), and the interface with vacuum V (gray); (b) An example of a surface configuration of the target system; (c) Six types of local configurations identified in the surface layer that were used in model optimization.

Results

In the first step, we used a reference system with all interactions set to zero, and optimized 2 the interaction parameters by minimizing statistical distance sIM based Eq. (15) using the

statistics of six surface configurations (Figure 8c) obtained from imaging type of data (Figure 8b). The perturbation scan found a good match with the target properties (Figure 9a) after a single iteration, with the resulting parameters listed in Table 3. While the probability of different surface configurations is reproduced well, the concentration profile across the layers, which was


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 48

not included in the target data, is only accurate near the surface, with the accuracy lost at the substrate interface (Figure 9b). The second perturbation iteration did not further improve the match of the model parameters with those of the target (Table 3), suggesting that the statistical uncertainty limits further optimization. This lack of constraints for the substrate interaction 2 in parameter wLS is clearly manifested in the shape of the landscape of statistical distance sIM

Figure 10a, which exhibits a valley that is extended in the wLS dimension, indicating that a wide range of substrate interaction parameters can reproduce the surface configurations well.

Figure 9. (a) Histograms of local surface configurations C1-C6 (see Figure 8c) obtained from 2 2 2 the target (black), sIM (green) and SE2 = ( sIM + sPR ) 2 (blue) systems. (b) The concentration

profiles of particle C in the same three systems.

Inclusion of the concentration profile statistics across the thin-film layers through distance 2 2 2 based on Eq. (15) into a combined loss function SE2 = ( sIM + sPR sPR ) 2 resulted in radical

improvement of the predicted concentration profile (Figure 9c) and better constrains on the substrate interaction parameter wLS (Table 3). This is reflected in the plot of statistical distance


36

Page 37 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


landscape, Figure 10b, which exhibits a valley that is now constrained also in the wLS dimension.

Figure 10. The landscapes of statistical distances in the space of wLL and wLS parameters based 2 on (a) surface structure data, sIM , and (b) combined surface structure and concentration profile

(

2 2 2 data, SE = sIM + sPR

)

2 . The red and white crosses indicate the optimized and target system

parameters, respectively.

2 Table 3. The interaction parameters of the target system and two models based on sIM and SE2

loss functions, which were optimized in one and two perturbation iterations. Parameters obtained after the second iteration are given in parentheses. Parameters wCC and wSC were set to zero, as these only define the arbitrary zero of the total potential energy. Energy is expressed in reduced units relative to ε.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

wCL

wLL

wLS

Target

0.2

0.44

0.4

2 sIM

0.219

0.449

0.226

(0.211)

(0.445)

(0.230)

0.201

0.434

0.382

(0.200)

(0.419)

(0.380)

2 2 SE2 = ( sIM + sPR )2

Page 38 of 48

4.2.2. Structure of binary alloys Target system

The microstructure of metal alloys has a profound effect on their macroscopic mechanical properties. It is therefore of great importance to be able to explore how the structure changes in response to changing conditions of their preparation and processing. Such exploration can be performed using a simple lattice model with its effective interactions optimized to match the available statistics of observed atomic structures. As a target system we chose an example from an interesting class of AB and A3B type binary alloys exhibiting long-period lamellar modulated structures.42,43 The target data consisted of three-dimensional concentration profiles, such as those that can be obtained from atom probe tomography (APT),44 which is routinely used to study metal alloys.45,46 The target system interactions were represented by the axial next-nearest neighbor Ising (ANNNI) model commonly used to study alloys as well as magnetic and ferroelectric systems.47,48 The configurational energy of these systems is defined by three types of interactions according to the general form of Eq. (22) as

ui = J 0

∑δ NN,x,y

AB

+ J1 ∑ δ AB + J 2 NN,z

∑δ

AB

(28)

NNN,z


38

Page 39 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where J0, J1, and J2 are interaction energy parameters between the pairs of atoms occupying specific nearest and next-nearest neighbor lattice sites as illustrated in Figure 11a; δ AB =1 for pairs of sites of the same type (AA or BB) and δ AB = −1 for pairs of sites of different type (AB). For the target system we set the interaction parameters to J0 = 1.0, J1 = 1.0, J2 = -0.6 (energy is expressed in reduced units relative to ε). Below the critical temperature this setting results in the formation of modulated phases,48 as shown in Figure 12a. The pseudo-experimental target data were generated by MC simulations in a periodic box of dimensions 10x10x40 over 106 individual steps at reduced temperature kBT/ε = 3.5. A single MC step consisted of a trial switching of the type (A or B) of a randomly chosen particle. Randomly selected configurations, such as the one shown in Figure 12a, were used to calculate the statistics of selected structural descriptors that were to be reproduced by a model. The descriptors included (i) the concentration profile along the axial direction (perpendicular to the lamellas), i.e., the probability of finding a particle of type A in a specified layer along this direction, Figure 12b, and (ii) the probabilities of finding pairs of particles of the same type in the nearest and next-nearest positions (corresponding to J0 and J2 positions in Figure 11b) as a function of the axial coordinate. Models

We tested two models against the target ANNNI system: (i) the ANNNI model itself to investigate identifiability, and (ii) the Ising model with isotropic competing interactions (IICI).49,50 The Hamiltonian of the latter model has the same general form defined by Eq. (28), but is symmetric in all three dimensions, with the first summation running over six nearest neighbor pairs, the second summation running over 12 neighbors sharing the cubic edge, and the third summation runs over the six next-nearest neighbor pairs (see Figure 11b). The IICI model is also known to form lamellar structures for certain interaction parameter combinations.49,51 The


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 48

simulations of the model systems followed the same protocol as the target system simulations, but were performed over 107 steps in a periodic simulation box of dimensions 6x6x17, which is compatible with the observed period of the modulated phase oscillations of the target system.

Figure 11. A schematic representation of pair interactions in the ANNNI (a) and IICI (b) models. The lattice site represented by the black square interacts with different types of neighbors according to the indicated interaction parameters J0, J1, and J2. Results

For this analysis, the loss function consisted of the sum of statistical distances for all ne structural descriptors, where the individual descriptor is the probability of finding a particle of type A in a specific layer along the z-axis, or finding a pair of A particles in the nearest or nextnearest positions along the axial and lateral directions, i.e., ne

SE2 = ∑ sl2 ne

,

(29)

l=1

where sl2 are statistical distances between the binomial probability distributions (i.e., the probability of a particle or a pair of particles being present or not) of the individual descriptors. The long-range correlations exhibited by the target system, which lead to the formation of


40

Page 41 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


modulated phases, present a significant challenge for the perturbation technique as a method that relies on the occurrence of a significant configurational space overlap between the reference and target systems. It is clear that starting from a simulation producing completely random configurations as in the previous example will not yield such an overlap; consequently, it is unlikely that an initial guess of the interaction parameters will lead to the formation of the correct modulated phase. In fact it is known that the ANNNI systems can form a large variety of longrange structures sensitive to small changes in the interaction parameters.49,51-53 For this reason we proceeded in steps, starting from a preliminary reference system with all pair interaction parameters set to zero, but with all particles subject to an external sinusoidal potential reproducing the observed concentration profiles along the axial direction. Subsequently, we gradually turned the external potential off in several steps and optimized the pair interaction parameters to maintain the axial concentration profile while trying to match the short-range nearest and next-nearest neighbor correlations of the target system. This procedure is similar to the calculation of crystal free energy by using the Einstein crystal as a reference, while gradually turning on particle-particle interactions.54 In the present case, however, the inter-particle interactions were not known beforehand, and therefore the perturbation technique had to identify them. Following this procedure, we were able to identify the original ANNNI interaction energy parameters (J0, J1, and J2) to high accuracy in six iterations (Table 4). Using the same procedure, we were able to optimize the IICI model parameters (Table 4) in 11 iterations. The resulting model reproduces the average modulated structure, Figure 12b, but exhibits larger structural fluctuations. Given that a wider range of parameters were able to reproduce the target data with similar accuracy, it may be possible to look for additional constraints by extending the set of target properties, should that be needed or desired.


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 48

Figure 12. (a) An example of a lamellar modulated structure generated by the target system representing a binary alloy with particles A shown in green and particles B shown in blue. (b) The average concentrations of particles A in layers along the axial z-direction for the target system (black) and the IICI model (green). The ANNNI model predictions are nearly identical to the target system (not shown). Solid lines indicate the sinus function fit of the concentration profiles. (c) The probability of finding pairs of A particles in the nearest (±1) and next-nearest (±2) positions in lateral (x) and axial (z) directions. This type of histograms was collected for each layer; the presented example is for layer 17.


42

Page 43 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4. Interaction parameters for the target system, Eq. (28), and optimized parameters of the ANNNI and IICI models. J0

J1

J2

Target (ANNNI)

-1.0

-1.0

0.6

SE2 (ANNNI)

-0.999

-1.006

0.601

SE2 (IICI)

-3.568

0.652

0.514

5. Conclusions We have presented a theoretical framework for consistent integration of multiple sources of experimental data and ab initio results into statistical mechanical models. The guiding principle underlying the methodology is the acknowledgement that measurement results, along with their unavoidable resolution and sampling characteristics, are the only way to gain information about physical systems. By decomposing any measurement into its atomic components, i.e., individual samples, we are able to express arbitrary results in terms of vectors in a Hilbert space and associate operations such as coarse-graining with simple vector algebra operations. This representation naturally leads to the model optimization loss function based on a metric known as statistical distance, which can be seen as a measure of distinguishability between measurement outcomes expressed as series of samples. This general approach incorporates target data with their full statistical characteristics and avoids introducing arbitrary weighting parameters. We have demonstrated that the model optimization methodology is applicable to both fluid and solid systems, and to experimental properties as diverse as electron microscopy images and the


43


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 48

fluid equation of state (in fact any measurable property). When combined with a physically sensible functional form of the particle interaction potential, the approach is able to generate predictive models, which may not be the case with some other popular force field development approaches. The proposed framework can not only be seen as a means for designing interaction models for atomistic and coarse-grained simulations, but also as a tool for integrating diverse data pertaining to a physical system in a compact form ready for making predictions, testing hypothesis, and designing experiments. Even though our current examples have dealt only with static properties, an extension to dynamic systems is straightforward, as indicated in Section 2.3. One of the most important issues that need to be addressed in applications to real experimental data, which were intentionally omitted in this study, is the role of noise generated by particular experimental setups. Special considerations are also needed when dealing with quantum uncertainty that introduces a probabilistic aspect in the sampling process, such as in diffraction experiments. In these cases extensions to the present formalism will need to include the additional layers of uncertainty. However, these issues are beyond the scope of the present article and will be dealt with systematically in upcoming studies.

ACKNOWLEDGMENTS Research on solid-state materials modeling was supported (L.V., R.KV, S.J., S.V.K) by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Research on fluid model design was sponsored (L.V.) by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division.


44

Page 45 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES (1) Ercolessi, F.; Adams, J. B. Interatomic Potentials from 1st-Principles Calculations - the Force-Matching Method. Europhys. Lett. 1994, 26, 583-588. (2) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. Effective Force Fields for Condensed Phase Systems from Ab Initio Molecular Dynamics Simulation: A New Method for Force-Matching. J. Chem. Phys. 2004, 120, 10896-10913. (3) Noid, W. G.; Chu, J.-W.; Ayton, G. S.; Krishna, V.; Izvekov, S.; Voth, G. A.; Das, A.; Andersen, H. C. The Multiscale Coarse-Graining Method. I. A Rigorous Bridge between Atomistic and Coarse-Grained Models. J. Chem. Phys. 2008, 128, 244114. (4) Shell, M. S. Systematic Coarse-Graining of Potential Energy Landscapes and Dynamics in Liquids. J. Chem. Phys. 2012, 137, 084503. (5) Dama, J. F.; Sinitskiy, A. V.; McCullagh, M.; Weare, J.; Roux, B.; Dinner, A. R.; Voth, G. A. The Theory of Ultra-Coarse-Graining. 1. General Principles. J. Chem. Theory Comput. 2013, 9, 2466-2480. (6) Moore, T. C.; Iacovella, C. R.; McCabe, C. Derivation of Coarse-Grained Potentials Via Multistate Iterative Boltzmann Inversion. J. Chem. Phys. 2014, 140, 10, 224104. (7) Reith, D.; Putz, M.; Muller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624-1636. (8) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020-9041. (9) Borreguero, J. M.; Lynch, V. E. Molecular Dynamics Force-Field Refinement against Quasi-Elastic Neutron Scattering Data. J. Chem. Theory Comput. 2016, 12, 9-17. (10) Bai, P.; Tsapatsis, M.; Siepmann, J. I. Trappe-Zeo: Transferable Potentials for Phase Equilibria Force Field for All-Silica Zeolites. J. Phys. Chem. C 2013, 117, 24375-24387. (11) Wang, L. P.; Martinez, T. J.; Pande, V. S. Building Force Fields: An Automatic, Systematic, and Reproducible Approach. J. Phys. Chem. Lett. 2014, 5, 1885-1891. (12) Hulsmann, M.; Reith, D. Spagrow-a Derivative-Free Optimization Scheme for Intermolecular Force Field Parameters Based on Sparse Grid Methods. Entropy 2013, 15, 36403687. (13) Stobener, K.; Klein, P.; Horsch, M.; Kufer, K.; Hasse, H. Parametrization of TwoCenter Lennard-Jones Plus Point-Quadrupole Force Field Models by Multicriteria Optimization. Fluid Phase Equilib. 2016, 411, 33-42. (14) Stobener, K.; Klein, P.; Reiser, S.; Horsch, M.; Kufer, K. H.; Hasse, H. Multicriteria Optimization of Molecular Force Fields by Pareto Approach. Fluid Phase Equilib. 2014, 373, 100-108. (15) Dunn, N. J. H.; Noid, W. G. Bottom-up Coarse-Grained Models That Accurately Describe the Structure, Pressure, and Compressibility of Molecular Liquids. J. Chem. Phys. 2015, 143, 16, 243148. (16) Mayne, C. G.; Saam, J.; Schulten, K.; Tajkhorshid, E.; Gumbart, J. C. Rapid Parameterization of Small Molecules Using the Force Field Toolkit. J. Comput. Chem. 2013, 34, 2757-2770.


45


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 48

(17) Vlcek, L.; Chialvo, A. A. Rigorous Force Field Optimization Principles Based on Statistical Distance Minimization. J. Chem. Phys. 2015, 143, 15. (18) Wootters, W. K. Statistical Distance and Hilbert Space. Phys. Rev. D 1981, 23, 357-362. (19) Stillinger, F. H. Effective Pair Interactions in Liquids - Water. J. Phys. Chem. 1970, 74, 3677-3687. (20) Stillinger, F. H. Density Expansions for Effective Pair Potentials. J. Chem. Phys. 1972, 57, 1780-1787. (21) Vlcek, L.; Sun, W.; Kent, P. R. C. Combining Configurational Energies and Forces for Molecular Force Field Development. J. Chem. Phys. 2017. (22) Braunstein, S. L.; Caves, C. M. Statistical Distance and the Geometry of Quantum States. Phys. Rev. Lett. 1994, 72, 3439-3443. (23) Rao, C. R. Maximum Likelihood Estimation for the Multinomial Distribution. Sankhya 1957, 18, 139-148. (24) Atkinson, C.; Mitchell, A. F. S. Rao's Distance Measure. Sankhya 1981, 43, 345365. (25) Rao, C. R. Information and the Accuracy Attainable in the Estimation of Statistical Parameters. Bull. Calcutta Math. Soc. 1945, 37, 81-91. (26) Kass, R. E. The Geometry of Asymptotic Inference. Stat. Sci. 1989, 4, 188-234. (27) Bhattacharyya, A. On a Measure of Divergence between Two Statistical Populations Defined by Their Probability Distributions. Bull. Calcutta Math. Soc. 1943, 35, 99110. (28) Thacker, N. A.; Bromiley, P. A.; Trucco, E. In Joint Meeting of the BMVA and the Royal Statistical Society 2004. (29) Landau, L. D.; Lifshitz, E. M. Statistical Physics; Elsevier, 1980; Vol. 1. (30) Chialvo, A. A. Excess Properties of Liquid-Mixtures from Computer Simulation a Coupling Parameter Approach to the Determination of Their Dependence on Molecular Asymmetry. Mol. Phys. 1991, 73, 127-140. (31) Vlcek, L.; Chialvo, A. A.; Cole, D. R. Optimized Unlike-Pair Interactions for Water-Carbon Dioxide Mixtures Described by the Spc/E and Epm2 Models. J. Phys. Chem. B 2011, 115, 8775-8784. (32) Vlcek, L.; Chialvo, A. A. Single-Ion Hydration Thermodynamics from Clusters to Bulk Solutions: Recent Insights from Molecular Modeling. Fluid Phase Equilib. 2015, 10.1016/j.fluid.2015.05.048. (33) Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129, 124105. (34) Hestenes, M.; Stiefel, E. Methods of Conjugate Gradients for Solving Linear Systems. J. Res. Natl. Bur. Stand. 1952, 49, 409-436. (35) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proc. R. Soc. London, A. 1938, 168, 264-283. (36) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1-19. (37) Martin, M. G. towhee.sourceforge.net. (38) Paricaud, P.; Predota, M.; Chialvo, A. A.; Cummings, P. T. From Dimer to Condensed Phases at Extreme Conditions: Accurate Predictions of the Properties of Water by a Gaussian Charge Polarizable Model. J. Chem. Phys. 2005, 122, 244511.


46

Page 47 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(39) Chialvo, A. A.; Vlcek, L. Ewald Summation Approach to Potential Models of Aqueous Electrolytes Involving Gaussian Charges and Induced Dipoles: Formal and Simulation Results. J. Phys. Chem. B 2014, 118, 13658-13670. (40) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: Tip4p/2005. J. Chem. Phys. 2005, 123, 12, 234505. (41) Tselev, A.; Vasudevan, R. K.; Gianfrancesco, A. G.; Qiao, L.; Ganesh, P.; Meyer, T. L.; Lee, H. N.; Biegalski, M. D.; Baddorf, A. P.; Kalinin, S. V. Surface Control of Epitaxial Manganite Films Via Oxygen Pressure. ACS Nano 2015, 9, 4316-4327. (42) Broddin, D.; Vantendeloo, G.; Vanlanduyt, J.; Amelinckx, S.; Portier, R.; Guymont, M.; Loiseau, A. Long-Period Superstructures in Cu-3+/-Xpd. Philos. Mag. A-Phys. Condens. Matter Struct. Defect Mech. Prop. 1986, 54, 395-419. (43) Kulik, J.; Takeda, S.; Defontaine, D. Long Period Superstructures in Ag3mg. Acta Metall. 1987, 35, 1137-1147. (44) Moody, M. P.; Vella, A.; Gerstl, S. S. A.; Bagot, P. A. J. Advances in Atom Probe Tomography Instrumentation: Implications for Materials Research. MRS Bull. 2016, 41, 40-45. (45) Herbig, M.; Raabe, D.; Li, Y. J.; Choi, P.; Zaefferer, S.; Goto, S. Atomic-Scale Quantification of Grain Boundary Segregation in Nanocrystalline Material. Phys. Rev. Lett. 2014, 112, 5, 126103. (46) Jin, K.; Guo, W.; Lu, C. Y.; Ullah, M. W.; Zhang, Y. W.; Weber, W. J.; Wang, L. M.; Poplawsky, J. D.; Bei, H. B. Effects of Fe Concentration on the Ion-Irradiation Induced Defect Evolution and Hardening in Ni-Fe Solid Solution Alloys. Acta Mater. 2016, 121, 365373. (47) Selke, W. The Annni Model-Theoretical Analysis and Experimental Application. Phys. Rep.-Rev. Sec. Phys. Lett. 1988, 170, 213-264. (48) Selke, W.; Fisher, M. E. Monte-Carlo Study of the Spatially Modulated Phase in an Ising-Model. Phys. Rev. B 1979, 20, 257-265. (49) Upton, P.; Yeomans, J. An Ising-Model with Isotropic Competing Interactions and Binary-Alloys. Europhys. Lett. 1988, 5, 575-582. (50) Upton, P.; Yeomans, J. Model for Binary-Alloys - an Ising-Model with Isotropic Competing Interactions. Phys. Rev. B 1989, 40, 479-492. (51) Widom, B. Lattice Model of Microemulsions. J. Chem. Phys. 1986, 84, 69436954. (52) Selke, W.; Duxbury, P. M. The Mean Field-Theory of the 3-Dimensional Annni Model. Z. Phys. B-Condens. Mat. 1984, 57, 49-58. (53) Yeomans, J. The Theory and Application of Axial Ising-Models. Solid State Phys. 1988, 41, 151-200. (54) Frenkel, D.; Smit, B. Understanding Molecular Simulations: From Algorithms to Applications; Academic Press: New York, 2007.


47


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 48

For Table of Contents use only


48

Consistent Integration of Experimental and Ab Initio Data into Effective

Recommend Documents