{'fraction_non_alphanumeric': 0.0207179400448184, 'fraction_numerical': 0.006299945034036616, 'mean_word_length': 5.156168662155127, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ESPEI for efficient thermodynamic database development, modification, and uncertainty quantification: application to Cu-Mg

The software package ESPEI has been developed for efficient evaluation of thermodynamic model parameters within the CALPHAD method. ESPEI uses a linear fitting strategy to parameterize Gibbs energy functions of single phases based on their thermochemical data and refine the model parameters using phase equilibrium data through Bayesian optimization within a Markov Chain Monte Carlo machine learning approach. In this paper, the methodologies employed in ESPEI are discussed in detail and demonstrated for the Cu-Mg system down to 0 K using unary descriptions based on segmented regression. The model parameter uncertainties are quantified and propagated to the Gibbs energy functions.


Introduction
Thermodynamics is the foundation of science that concerns the state of a system under given external conditions. The CALPHAD method describes the thermodynamic behavior of a system by a semi-empirical parametric model of the Gibbs energy of each phase in the system.
Developing models and evaluating model parameters is central to the CALPHAD method, but, with the exception of some recent attempts to study the influence of model parameterization in CALPHAD models computationally, the state of the art user tools [1,2] for fitting CALPHAD model parameters have remained relatively stagnant over the years.
Following the publication of the SGTE description of the unary systems [3], significant effort has been made to develop multicomponent databases based on descriptions of the common constituent subsystems. This has led to the convergence of the CALPHAD community on certain key unary and binary descriptions, which adds considerable inertia to updating the descriptions with the most recent scientific understanding. For example, new unary descriptions are being developed that extend down to 0 K using the Einstein [4] or Debye [5] models, include better descriptions of the high temperature behavior of the heat capacity for solid phases, or that take into account variability in the reported data by outlier detection [6] or automatic weighting schemes. In addition, only recently has the topologically close-packed s phase been modeled with five sublattices [7][8][9] corresponding to the five Wyckoff positions in that phase, but nonetheless there are still new assessments published where the s phase is described with fewer sublattices to stay compatible with previous and future assessments [10,11].
With the rapid increase of available computational data and the advances in physicallybased CALPHAD models for Gibbs energies and other physical properties [12], the current methods for optimizing model parameters to data cannot support the growing need to assess and 3 compare different types of models to data and to maintain the existing databases when new data must be incorporated. Furthermore, a key aspect of communicating the efficacy of newly developed models and databases to consumers of CALPHAD databases is to robustly quantify the uncertainty of the model with respect to the underlying data. Therefore, it is necessary to develop a new generation of software tools that can support not only model development, database development and maintenance, but also to quantitatively determine and report the confidence in the optimized models and derived properties to database users.
Here we introduce ESPEI (Extensible Self-optimizing Phase Equilibria Infrastructure) as an open-source, Python-based application for CALPHAD database development and uncertainty quantification. The rest of this paper will develop the theory behind the two fitting modules of ESPEI, single-phase parameter selection and multi-phase Markov Chain Monte Carlo (MCMC) optimization and uncertainty quantification. Finally, ESPEI will be applied to select model parameters and optimize the parameters in the Cu-Mg binary system, propagating the parameter uncertainties to the Gibbs energy functions and the phase equilibria.

Theory
ESPEI implements two steps of model parameter evaluation: generation and MCMC optimization. The parameter generation step uses [13] experimental and first-principles data describing the derivatives of the Gibbs free energy to parameterize the Gibbs energies of each individual phase within the Compound Energy Formalism (CEF) [14], giving the user a complete thermodynamic database based only on thermochemical data. Experimental thermochemical data for virtually all real alloy systems are too sparse to fully describe the Gibbs energies of the phases and are often unable to access the energetics of metastable configurations defined within 4 the CEF. In practice, the CALPHAD method requires that available experimental data is augmented with thermochemical data from estimates, empirical models [15], machine learning models [16], or first-principles calculations [17], all of which only give approximate energies.
Therefore, assessments based on thermochemical data alone cannot reproduce the true Gibbs energies of each phase and, by extension, the phase diagram and other thermodynamic properties. For this reason, the CALPHAD method has modeled the Gibbs energies by fitting to all available thermodynamic data, including both thermochemical and phase equilibrium data, since its inception. The phase equilibrium data require the Gibbs energies of all phases to be defined, so parameters must be optimized iteratively to be self-consistent by the modeler. The MCMC optimization implemented in ESPEI fulfills the role of iterative, self-consistent optimization of Gibbs energy parameters considering all data simultaneously. The following sections will discuss each fitting step in detail.

Parameter selection from thermochemical data
Implementing parameter selection with only thermochemical data allows the Gibbs energy models of each phases to be evaluated directly from the data, therefore models for each phase can be generated independently. The construction of the Gibbs energy within the CEF involves the extrapolation of the low order terms in composition into the interaction space, so all endmembers must be fit first, followed by the excess parameters from low order to high order.
Interaction parameters are evaluated with candidate models from zeroth order to fourth order interactions, fitting all the temperature dependent terms in a Redlich-Kister polynomial for the zeroth order interaction before fitting all the Redlich-Kister terms for the first order interaction and so on. Endmember energies and interactions are evaluated in the following sequence, using a 5 hypothetical three sublattice ternary solution phase as a prototype for a phase with v-sublattices and j-components. , and so on. The parameters that will be considered are limited by the input thermochemical data and how the derivatives of the Gibbs energy can describe the parameters. In general, candidate models that could not physically describe the underlying Gibbs energy function are not generated. For example, isothermal enthalpy of formation data cannot describe high order terms of the Gibbs energy, and thus they will not be considered when only enthalpy data is used.
In this polynomial, all the terms are linear with respect to the model coefficients, so the candidate models can be fit to input thermochemical data by a linear least squares optimization under the assumption that the residuals follow a Gaussian distribution. A tunable hyperparameter 6 for L2 regularization that penalizes the magnitude of the parameters can be adjusted by the user to prevent overfitting, resulting in the objective function given in Eq. 2, where is the vector of coefficients, the matrix of polynomial terms, the data to be fit and the hyperparameter controlling the mixing between the fit to the data and the magnitude of the parameters.
For each fitting step, the parameters for are the candidate models are evaluated by minimizing Eq. 2, then the residual sum of squares (RSS) between the evaluated parameters and the data are used to score and compare the models within the corrected Akaike information criterion (AICc) [18]. The model with the lowest score is the optimal combination of model fitness and complexity. The AICc is a modified version of the AIC that avoids overparameterization when the data is sparse, which is often the case for thermochemical data.
The AICc score is calculated by Eq. 3, where is the number of parameters, the number of samples, and RSS is the residual sum of squares between the model predictions and the data. In this way, models for each phase with a reasonable fit/complexity tradeoff can be built up one parameter at a time.

Parameter optimization to phase equilibria data
The CEF uses sublattices to describe the chemical behavior of occupancy and substation of elements on sublattices, typically representing the Wyckoff positions in crystalline materials.
The endmembers that the CEF introduces into the Gibbs energy model are often not stable or cannot be characterized completely experimentally due to the high dimensionality of the sublattice composition space that is not directly accessible by experiments. Measurements of the

Eq. 3
Gibbs energy derivatives in CEF models cannot be obtained completely from experiments in the general case, requiring the use of thermodynamic data indirectly related to the Gibbs energies of the phases to refine the parameters generated from thermochemical data to construct a selfconsistent set of parameters that reproduce the known data for a given system. In many practical binary and ternary systems of interest, a typical CALPHAD description might have on the order of tens to hundreds of parameters to assess. Current optimization methods require modelers to restrict the degrees of freedom to only a handful of simultaneous parameters and optimize the restricted sub-problems iteratively.
Bayesian optimization is a method of optimizing model parameters by Bayes Theorem [19], obtaining the posterior probability, ( | ), of the model parameters, , conditioned on the data, , in Eq. 4. The prior, ( ), contains the domain knowledge of the modeler in the probability distribution of each parameter. The likelihood, ( | ), is the probability that the data is described by a set of parameters, and the evidence, ( ), is the probability of the data marginalized over the possible parameters. In CALPHAD databases, evaluating a closed form solution for these quantities over all model parameters is not computationally achievable, so the posterior probability is determined numerically using MCMC.
( | ) = Eq. 4 ESPEI uses MCMC to perform a Bayesian optimization of all model parameters simultaneously. In principle, MCMC optimization in ESPEI can be performed for arbitrarily sized multicomponent, multiphase systems with any number of degrees of freedom, however this is not computationally feasible in practice and it is expected that the extrapolation of binary and ternary CALPHAD parameters into multicomponent space precludes the need to model all the 8 degrees of freedom in a multicomponent system simultaneously. However, a challenge in applying MCMC to optimize parameters in CALPHAD databases is that most MCMC samplers assume that model parameters are uncorrelated to efficiently explore parameter space [19]. The parameters in CALPHAD models for each phase are correlated to each other because increasing the value one parameter and decreasing another can give the same Gibbs energy for any given set of conditions. This challenge is addressed in ESPEI by using an ensemble sampler, as introduced by Goodman and Weare [20]. Normal priors follow a Gaussian distribution, usually centered on the parameter. Normal priors are useful for informing a centrality of each parameter in the prior, but not explicitly limiting the parameter bounds. Triangular priors are particularly useful for modeling in ESPEI, as they allow the modeler to inform the centrality of the parameters, e.g. from parameter selection, but also to enforce bounds where the prior probability is zero outside of those bounds, which can be used to ensure that parameters stay the same sign or within the same order of magnitude, for example.
The hyperparameters to generate the prior distributions are specified in ESPEI input files. A 10 detailed example and further explanation are included in the ESPEI software documentation [22].
Three main types of data are considered by the likelihood function defined in ESPEI: single phase thermochemical data of the temperature derivatives of the Gibbs energy, activity data, and multi-phase equilibria data. To ensure consistent weighting of error from different data types, the likelihood for each type of data is normalized by the standard deviation of the error for that data type. For all data types, the error is assumed to follow a normal distribution. ESPEI provides default values for the standard deviations of each type of data, however users can modify the values by adding a weight for each type of data or for each individual dataset. The thermodynamic quantities for each type of data are determined in ESPEI by using pycalphad [23] as the thermodynamic calculation engine.
The error due to thermochemical data is calculated by comparing the difference between the predicted values and the measured or calculated values. Since these values are calculated or measured directly, the standard deviations of the error are relatively well defined. The default standard deviations for enthalpy, entropy, and heat capacity data are chosen to be 500 J/mol, 0.2 J/K-mol, and 0.2 J/K-mol, respectively. To fit activity data, the measured activity is converted to chemical potentials as suggested by Lukas et al. [24]. The error is calculated by the difference between calculated and measured chemical potentials. Although activities are not measured directly, the experimental errors are mathematically related to the activity and the experimental errors reported in the literature can be propagated to the chemical potentials, ESPEI chooses a standard deviation in the chemical potential due to measured activity data of 500 J/mol.
Thermodynamic assessments are judged, to the first order, by the agreement with the measured phase equilibria on the phase diagram, so any complete thermodynamic assessment 11 program must consider optimizing to phase equilibria data. The phase equilibria depend on solving for the global minimum energy for many phases and arbitrary models of the Gibbs energy, so a solution mapping the model parameters to the calculated equilibria is not well defined. Several methods exist for determining the error for phase equilibria. The direct approach is to calculate the phase equilibria and determine the error in the phase boundary due to a different in temperature or composition. However, a drawback of the direct approach is that the phase equilibria of interest are not always stable and a near-optimal solution must already be known to properly calculate the error. Since MCMC optimization considers many degrees of freedom and the initial database may not be close to a globally optimal solution, a more generic approach is warranted. The method implemented by ESPEI is similar to the rough search method implemented in the PanOptimizer software [2]. Given a set of measured phase constitutions in equilibrium (e.g. from EPMA) and a set of candidate parameters for a model the multiphase error can be calculated as follows (illustrated schematically in Figure 1

Cu-Mg system
The Cu-Mg system exhibits features found in most complex CALPHAD databases including terminal phases with solubility, a stoichiometric compound, and an intermetallic compound with solubility. This binary system is a technologically important system for bulk metallic glass (BMG) systems [25] that are produced using molds that are actively cooled, for example with liquid N 2 and for Al-based alloys [26]. CALPHAD models that extend down to low temperatures can be used to describe the energetics for phase precipitation in potential BMG alloys. Because ESPEI can rapidly select and fit different CALPHAD models, it is ideal for comparing different models and unary reference states. Here ESPEI will be used to develop a description for the Cu-Mg system that extends down to 0 K using segmented regression models for the Cu and Mg 13 unary descriptions [5]. The equilibrium Cu-Mg system has four stable solid phases: fcc and hcp solutions, an intermetallic compound, CuMg 2 , and a C15-type Cu 2 Mg Laves phase with considerable solubility. The Cu-Mg system has previously been modeled by several times [27,28], most recently by Zhou et al. [29] [30]. Interested readers are directed to those publications for a thorough review of the existing data.
Recently, first-principles calculations for the terminal solution phases in the Cu-Mg system were calculated by Gao et al. [31]. These data are also considered in our parameter selection and optimization. The latest programming interface, command line usage, and additional documentation including code examples, a changelog of ESPEI versions and links to the development repository can be found at https://espei.org [22]. The documentation and user manual for ESPEI version 0.6.2 is included in the supplemental material.

Unary reference data
The segmented regression (SR) model [5] is one of the alternative physically-based formulations to the recently used SGTE description of reference data for pure elements. The SR model has been developed for the description of temperature dependence in the heat capacities down to 0K considering relevant physical effects and it consists of three terms: is the vector of unknown model parameters to be estimated from experimental and computational data. Here, the main contribution to the heat capacity due to phonon vibration is represented by the well-known Debye model. The second term in Eq. 6 is a so-called bent-cable model. This term is introduced for the decomposition of linearly-dependent physical effects at low and high temperatures. The last term describes the magnetic contributions to the heat capacity. This term is not considered and equal to zero for non-magnetic elements.
Since the analytical formulation of the Debye model cannot be implemented directly in the TDB format, a non-series approximation of the Debye model base on the weighted linear combination of Einstein functions for heat capacity has been developed to overcome this issue and it has been successfully applied for the thermodynamic re-evaluation of the Cr-Nb system [32] and agrees well with available experimental data over the entire temperature range.
The Gibbs energy function for reference data is derived based on the well-known thermodynamic relationship between heat capacity and Gibbs energy, The heat capacities fit for Cu and Mg, presented in Figure 2, and the Gibbs energies, presented in

ESPEI optimization
Parameters were generated for the Cu-Mg system based on the measured and calculated single phase thermochemical data. Experimental enthalpy of formation for the CuMg 2 [33,34] and Laves C15 [33,34] phases and experimental enthalpy mixing for the liquid phase [35] were used. First-principles calculations for the mixing energies for the FCC_A1 and HCP_A3 phases were calculated by Shin [36] and Gao et al. [37], while the endmembers of the Laves phase and CuMg 2 phase were calculated by Zhou et al. [29]. In total, 15 degrees of freedom were generated for the endmembers and binary interactions for the five phases. The resulting phase diagram is shown in Figure 4. This phase diagram, produced only from the thermochemical data, shows the qualitative behavior of the measured phase equilibria, even though those data were not included in the optimization. The parameters that produced this phase diagram were used as input for the MCMC optimization step with triangular priors that spanned from ±0.5 . All the available thermochemical, activity and phase equilibrium data were considered, and all parameters fit to the data simultaneously. Figure 5 shows the phase diagram after 500 MCMC iterations and shows excellent agreement with the phase equilibria data. To correctly quantify uncertainty, each chain must be converged. MCMC convergence cannot be proved for arbitrary simulations, however several criteria exist that are indicators of convergence. The application of these indicators is the same in ESPEI as any MCMC simulation and are discussed in detail elsewhere [19].
The liquid parameters from the Markov chain constructed in the MCMC optimization are shown in the corner plot in Figure 6a. The histograms for each parameter are plotted along the diagonal of the corner plot and the covariance between each pair of parameters is plotted below the diagonal. The histograms of each parameter are not required to follow any distribution because the probability is the product of the likelihood and the prior, which can have varying shape in parameter space. For parameters in CALPHAD models, it is expected that some parameters will exhibit covariance, particularly for parameters in the same phase and between phases in equilibrium. Too many parameters that are correlated are an indication that fewer parameters are required to describe the Gibbs energies that the data require. Uncertainty from the MCMC chains and probabilities were propagated to the Gibbs energy of mixing of the liquid phase in Figure 6b. To propagate the uncertainty, 95% of the parameter values with the highest probability density were used to calculate the Gibbs energies and phase diagrams, forming the 95% highest density intervals (HDI) for each type of data. The mean Gibbs energy is plotted as the line, and the 95% HDI is plotted as the shaded region in Figure 6b. The uncertainty in the Gibbs energies of the liquid region are smaller where liquid is in stable two phase equilibrium on the phase diagram and where the different sets of experiments agree. Near the congruent melting of the Laves phase, there is some disagreement in the data, which leads to a region of relatively larger uncertainty. The maximum uncertainty is on the order of 1 kJ/mol, which is consistent with what would be expected from experiments. Since the unary descriptions are fixed in this optimization, there is no uncertainty bound at either pure Cu or pure Mg will always be 0. Fitting unary parameters in ESPEI would enable the uncertainty quantification. More details on the statistical approach of propagating parameter uncertainty, especially to phase diagrams, will be published in a separate article [38].

Conclusion
ESPEI allows thermodynamic databases to be rapidly generated using thermochemical data, then refined considering all thermodynamic data for a system simultaneously with little user interaction. The parameter selection and MCMC optimization steps in ESPEI make it the first CALPHAD-based software that can generate model parameters, adjust the parameters, and quantify uncertainty for arbitrarily sized multicomponent systems. The two-step optimization approach allows existing and newly created databases to be updated when new experimental or calculated data becomes available. The Cu-Mg system was assessed using ESPEI and a new method for depicting phase diagram uncertainty was presented. ESPEI is developed in the open with the intent that it become adopted as a platform for the development and assessment of Gibbs energy models and other property models.

Acknowledgements
The authors thank Jiong Wang, ShunLi Shang, and Yi Wang for their help in testing ESPEI.