Modeling thermophysical properties of glasses

Metal oxide glasses are important in various industries because their properties can be tailored to meet application-specific requirements. However, there are few rigorous modeling tools for predicting thermomechanical properties of these materials with acceptable accuracy and speed, yet these properties can play a critical role in material design. In this article, a general multi-scale modeling framework based on Monte Carlo simulation and a cubic equation of state for predicting thermomechanical properties is presented. There are two novel and fundamental aspects of this work: (1) characterization of glass transition and softening temperatures as adjacent saddle points on the heat capacity versus temperature curve, and (2) a new moving boundary equation of state that accounts for structure and ‘soft’ repulsion. In addition, modeling capabilities are demonstrated by comparing thermomechanical properties of a pure B2O3 glass and PbO–B2O3 glass predicted by the equation of state to experimental data. Finally, this work provides a rigorous approach to estimating thermophysical properties for the purpose of guiding experimental work directed at tailoring thermomechanical properties of glasses to fit applications.

1. It uses the Gibbs-Helmholtz equation to constrain the energy parameter. ∂V 2 ) Tc , on the critical isotherm equal to zero. This results in a fixed boundary condition for the energy parameter given by a T c , p c = 0.42748 R 2 T 2.5 c p c , which is not applicable to glasses and metal oxide glasses because many of these materials decompose before ever reaching a critical point.
TlnT, www.nature.com/scientificreports/ 3. Monte Carlo simulations are used to gather only pure component molecular length scale information (i.e., internal energies of departure, U D ) over ranges of temperature and pressure and U D s are placed in lookup tables for use in bulk phase modeling and simulation. 4. Monte Carlo simulations need only be run once because U D lookup tables can be used over and over again.
Thus, the GHC framework is computationally efficient. 5. Mixing rules are used to predict properties of mixtures. 6. The parameters, a and b , are never regressed to experimental data. Experimental data is only used to validate predictions by the GHC equation.
Equation (1) has also been used for vapors, hexagonal ice, aqueous electrolytes, and gas hydrates. The objective of this article is to develop a thermodynamically rigorous modeling framework that addresses the open modeling challenges associated with the design of glassy materials. A general multi-scale equation of state framework for modeling and simulating glassy materials is presented that includes the use of Monte Carlo simulation to gather molecular length scale information and a novel moving boundary equation of state for calculating bulk thermomechanical properties of pure glasses and glassy mixtures. The numerical methodology and simulation results for a pure glass and binary glass mixture are compared to experimental data taken from the open literature where available. This is followed by a discussion of results and the conclusions of this work.

Materials and methods
Modeling glasses and their properties with an equation of state. Unfortunately, Eq. (1) breaks down for glassy materials (e.g., asphaltenes, ionic liquids, polymers, pure glass, metal oxide glasses, etc.) because these materials decompose prior to reaching a critical point. This, in turn, makes the boundary condition for the energy parameter in Eq. (1) undefined. To circumvent this Lucia and Gow 32,33 modified the expression for the energy parameter for pure glasses as follows: where the quantity a 0 T g , p g in Eq. (3) denotes the boundary condition for the energy parameter for pure glassy materials. Lucia and Gow 32 proposed and used an iterative numerical procedure to estimate a 0 T g , p g that minimizes the error between U D from Monte Carlo simulation and U D calculated from fundamental thermodynamics. However, this numerical procedure is computationally expensive and does not consider 'soft' repulsive forces 34 . Therefore, in this work, we propose the following expression for the boundary condition for the energy parameter for glassy materials Note that Eq. (4) has the following properties: 1. It is a closed form, temperature-dependent expression for the attraction parameter boundary condition, making it a moving boundary condition. 2. The first term on the right, −b G U DG T g , is fixed and replaces the computationally expensive numerical procedure described in Lucia and Gow 32 .
3. The second term on the right, b G R lnT ln2 (T − T g ) , is an ideal gas-like contribution that increases the boundary condition as temperature increases. Equation (5) is a closed form and computationally efficient expression for the attraction parameter that accounts for soft repulsive forces through the inclusion of internal energy and entropy of departure terms.
For mixtures, the molecular co-volume parameter, b M , is computed using The moving boundary condition for mixtures, a 0M , and the temperature dependent attraction parameter, a G (T), are computed using one fluid theory (i.e., by replacing b G with b GM , T g with T gM , and U D with U DG M in Eqs. 4 and 5), which gives the following moving boundary condition for a 0M : and the following expression for a GM (T): Molecular length scale information from Monte Carlo simulation. Accurate molecular length scale information is extremely important in building a robust bulk phase model for predicting thermomechanical properties. For the multiscale GHC equation, this means accurate pure component internal energies of departure, U D , are required. As noted earlier, this information is gathered a priori using Monte Carlo simulations that are run over ranges of temperature and pressure relevant to the application. Results of the Monte Carlo simulations are then placed in pure component lookup tables and used in bulk phase property determinations. In the multiscale GHC framework, pure glass and pure metal oxide properties tend to be very sensitive to values of internal energies of departure. This suggests that large numbers of Monte Carlo cycles and many sets are required to estimate values of U D .
Force field models. Force field models play a critical role in computing accurate molecular length scale information and many force fields (potentials) have been used for glasses. Table 1 gives a summary of some of the force field models that have been used to model thermomechanical properties for pure glasses and metal oxide glass mixtures using molecular dynamics and Monte Carlo simulation. As can be seen from Table 1, much of the modeling work in the open literature has focused on silica. However, there is no general agreement on the force field used to model glasses or metal oxide glasses.
In this work, we have chosen to use Lennard-Jones plus electrostatic forces to model the internal energy of departure where ε and σ are the well depth of the potential function and zero potential distance parameters respectively, q denotes a partial charge, ǫ 0 is the permittivity, r is the distance between particles (i.e., atoms), i and j are atom indices, and A is the total number of atoms in a molecule.
Force field parameters. Finding or determining reliable values of the parameters, ε , σ , and q , for glasses and metal oxides can, at times, be quite challenging. Many of the parameters published in the literature [e.g., Si and O parameters in Cruz-Chu et al. 35 , Butenuth et al. 36 , Emani et al. 42 , Anagnostopolous et al. 37 ] give values of U D and bulk phase properties that vary all over the map.
Numerical methodology. Before presenting simulation results it is instructive to describe the numerical methodology used in determining thermophysical properties of glass and metal oxide glass, which determines the glass transition and softening temperatures prior to estimating all other desired thermophysical properties. Motivation for the methodology to follow is illustrated using boron trioxide as an example. The zero potential Lennard-Jones parameters, σ , were taken from Tabrizi et al. 43 and the well depth parameters, ε , were taken from Soper 23 and adjusted to remove the EPSR repulsive term. Partial charges were selected to enforce electroneutrality. The B 2 O 3 force field parameters are shown in Table 2. www.nature.com/scientificreports/ NVT Monte Carlo simulations with N = 320 atoms at 1 bar pressure were run to gather molecular information. Smaller ensembles will give inaccurate molecular information at higher temperatures.
Heat capacity. Figure 1 shows a comparison of the Monte Carlo simulation predictions of specific heat capacity, C V , vs. temperature for boron trioxide with experimental data from D' Angelo and Carini 44 , where the predicted heat capacity at constant volume was computed using the fundamental expression ( ∂U ∂T ) V = ∂ ∂T (U ig + U D ) V . For B 2 O 3 glass, N , the number of atoms/molecule, is five and the degrees of freedom, DoF , for a polyatomic, nonlinear ideal gas varies from 6 at low temperature to 3N at high temperature. Therefore, U ig = ( DoF 2 )RT and U D comes from the Monte Carlo simulations. However, the results of the Monte Carlo simulations over-predict the maximum heat capacity (i.e., 2.3 J/g K vs. 2 J/g K) and slightly under-predict the heat capacity (i.e., 1.59 vs. 1.75 J/g K) over the range [600-850 K].
Glass transition and softening temperatures. One of the novel aspects of this work is the conjecture that glass transition and softening temperatures can be estimated using adjacent inflection points on each side of the maximum in C V (T) . Using the numerical procedure described in the Supplemental Material, adjacent saddle points on each side of the maximum in C V (T) give values of T g = 535.27K and T s = 660.78K. Volumetric coefficient of thermal expansion. Volume versus temperature results from the equation of state computations can be fit and then differentiated to give approximations of ( ∂V ∂T ) p so that a prediction of α can be computed from the equation

Pseudo-algorithm.
A general pseudo-algorithm for computing glass transition and softening temperatures along with other thermophysical properties is as follows: 1. Compute all pure component molecular properties (i.e., U D i ′ s , etc.) for the glass under consideration using Monte Carlo simulations over ranges of temperatures and pressures of interest, considering any structural changes that can occur. Typically, four sets of all-atom Monte Carlo simulations with periodic boundary conditions and 100,000 equilibration and 100,000 production cycles are run for each temperature and pressure and averaged. Store the molecular information in pure component lookup tables. Compute the mixture internal energies of departure using Eq. (7). 2. Using internal energies of departure and degrees of freedom compute the total internal energy, U , and heat capacity at constant volume, C V , over the desired range of temperature where U = U ig + U D , U ig = DoF 2 RT , and C V = ( ∂U ∂T ) V . 3. Compute the global maximum of C V and saddle points on each side of the global maximum using the procedure described in the Supplementary Material. The saddle points on each side of the maximum give predicted values of T g and T s , where T s > T g . 4. Given T g , compute the boundary condition for the glass using Eq. (4) (or Eq. 8 for mixtures). 5. Using the multiscale GHC equation of state, calculate the molar volume and all desired volume-dependent properties such as volumetric thermal expansion coefficient, isothermal compressibility, bulk modulus, elastic modulus, etc. over the temperature and pressure ranges of interest.
Steps 1 is a critical step of the modeling since it is well known that glasses can exhibit temperature-dependent changes in structure. For example, it is well known that boron trioxide is comprised of two types of elementary building blocks-BO 3 units and boroxol rings (B 3 O 6 ) and that the fraction of boroxol rings decreases with increasing temperature (i.e., from 0.6 at 500 K in the glassy region to 0.2 at 1900 K in the liquid state). See Fig. 2 and Ferlat 38 .
However, it is possible to incorporate structure in the Monte Carlo simulations by using an initial structure in the simulation box that reflects structure. In the case of B 2 O 3 that would mean initializing the atoms in the simulation box to reflect the ratios of the elementary building blocks BO 3 to B 3 O 6 .
Step 2 and 3 provide a way of directly calculating glass transition and softening temperatures for pure glasses or metal oxide glasses and set the stage for subsequently computing volume and all volume-dependent thermophysical properties in steps 4 and 5. However, it is important to understand that accurate Monte Carlo simulations must be performed in step 1 to get reliable estimates of U D (T) and the C V (T) functionality, particularly at high temperature. This means that several sets of Monte Carlo simulations with a large number of equilibration and production cycles and large ensemble are required. Heat capacity. The heat capacity versus temperature curve shown in Fig. 1a shows that the general trend and magnitude of calculated C V (T) in this work match the experimental C P (T) from D' Angelo and Carini 44 quite well.However, the results of the Monte Carlo simulations over-predict the maximum heat capacity (i.e., 2.3 J/g K vs. 2 J/g K), slightly under-predict the heat capacity (i.e., 1.59 vs. 1.75 J/g K) over the range [600-850 K] and show a peak slightly to the right of the peak given in the experimental data. Remember, the computed results are truly predictions with no regression to any experimental data. Also, note that there is no clear melting point in Fig. 1a or b because the material is glass. This suggests that the force field model for B 2 O 3 consisting of Lennard-Jones + electrostatic forces together with Monte Carlo simulation is capable of predicting the heat capacity of boron trioxide with reasonably good accuracy. In addition, it is important to note that the comparison of calculated values of C V in Fig. 1a with experimental values of C P shown in the Fig. 1b is valid because (a) �H = �U + �(pV ) , (b) the pV term is small, (c) the pressure is low (i.e., 1 bar), and (d) the glass phase is a condensed phase. Thus, U ≈ H is a good approximation, which implies C V ≈ C P .

Glass transition temperature. The saddle point (inflection point) to the left of the maximum in C V
shown in Fig. 1a gives a predicted glass transition temperature, T g , of 535. 26 45 to 571 K in Ref. 46 to as high as 580 K in Ref. 47 .
Softening temperature. Softening temperature is generally described as a transport property and characterized as the temperature at which the viscosity, η, of the glass is 10 7.6 poise (log 10 (η) = 7.6). Experimental viscosity data given in Ref. 48 in Table II, p. 616 shows that log 10 (η) of B 2 O 3 has a value of 7.6 at 375 °C (648 K). Like many properties of boron trioxide, there is limited experimental information for softening temperature and in the few articles containing T s data 48,53 , the reported results vary widely. Nevertheless, the moving boundary equation of state prediction of softening temperature in Fig. 1a is T s = 660.78K , which is in good agreement with the viscosity-based results reported in Table II on p. 616 in Ref. 48 of 375 °C (648 K)-less than 2% error. However, the estimated softening temperature of B 2 O 3 reported on p. 379 in Ref. 53 is much lower-315 °C (588 K) and, in our opinion, is much closer to the glass transition temperature. Table 3 illustrates the volume predictions of the multiscale GHC equation in the glassy region with and without soft repulsive effects. Figure 3 shows a comparison of multiscale GHC equation predictions of B 2 O 3 specific volume with and without soft repulsive forces with experimental volume data for the glassy and melt regions.

Volume predictions.
The comparison of predicted specific volume with experimental data shown in Table 3 and Fig. 3 are in good agreement with the five experimental data sets from the literature, which contain a total of 51 experimental data points and were taken from 10,[48][49][50][51] . The values of the AAD % error for the multiscale GHC equation of state without (i.e., ) and with soft repulsion and temperature-dependent boroxol ring structure (i.e., ) were 1.29 and 0.64% respectively for the temperature range [300, 1600 K] and shows that the inclusion of 'soft' repulsion and structure in the moving boundary condition and in the Monte Carlo simulations lead to improvements in model predictions of specific volume.Also shown in Fig. 3 are the molecular dynamics (MD) simulation  Figure 4 shows plots of the volumetric coefficient of thermal expansion for the moving boundary equation of state, along with the experimental α data reported in Ref. 52 and the coefficient of expansion data reported in Ref. 48 , which appears to have been computed from fits of density versus temperature.     Table I, p. 106 Spaght and Parks 52 provide different sets of results for α for two separate samples of B 2 O 3 glass for partially and carefully annealed samples in which volume was measured using dilatometry and α was approximated using the finite difference expression α = ( 2 V 1 +V 2 )( �V B2O3 T 2 −T 1 ) . However, Spaght and Parks present no volume data or experimental error statistics for the coefficients of expansion-so it is difficult to judge the accuracy of their data. However, for partially annealed and carefully annealed boron trioxide, Spaght and Parks show that α increases with increasing temperature over the range [122, 275 °C] corresponding to temperatures up to the glass transition temperature. Again, it is important to note that these authors present no raw volume data in their paper. Donoghue and Hubbard 53 used interferometry and report a melting point of 450.8 °C (see p. 375) and a softening temperature of 315 °C (see p. 379). However, they only describe the qualitative behavior of α for B 2 O 3 . Fajans and Barber 54 use electronic structure arguments and the experimental data in Ref. 53 and in Ref. 55 respectively to sketch plots of α(T) similar to those given in Ref. 52 for glassy B 2 O 3 and in Ref. 48 for liquid B 2 O 3 . Finally, in Fig. 3, p. 616 in Ref. 48 it is shown that α for boron trioxide decreases monotonically with increasing temperature over the range [400, 1400 °C], which corresponds to temperatures in the melt regime.
The moving boundary equation of state predictions of volumetric coefficient of thermal expansion are in good agreement with experimental α data reported in Refs. 48,52,54 both in temperature trends and in order of magnitude [10 -4 , 10 -5 °C −1 ]. Like the experimental data, for which there is considerable scatter, the moving boundary equation predictions show a rapid rise in α(T) with increasing temperature in the glassy region and a decline in α(T) with increasing temperature in the melt region.

Metal oxide glass: lead monoxide-boron trioxide. Mixtures of lead monoxide and boron trioxide
have been studied by Refs. [56][57][58][59][60][61] . To begin, to calculate glass transition and softening temperatures for mixtures of PbO and B 2 O 3 and generate C VM (x, T) vs. temperature curves similar to Fig. 1a, a force field and force field parameters for PbO are needed. Table 4 gives the Lennard-Jones 6-12 parameters and partial charges for PbO, which were taken from Refs. 62,63 . We first tested the PbO force field parameters by determining internal energies of departure as a function of temperature. We then used those internal energies of departure to calculate the densities of litharge and massicot using the moving boundary equation of state. In each case the calculated predictions were close to values reported in the literature. More specifically, literature values of densities for litharge  www.nature.com/scientificreports/ and massicot range from 9.14-9.53 and 9.56-10 g/cm 3 respectively (e.g., mindat.org) while the moving boundary equation of state predictions were 9.13 g/cm 3 at 300 K for litharge and 9.58 g/cm 3 at 800 K for massicot.
Heat capacity vs. temperature. Figure 5 shows a plot of C VM (x, T) vs. temperature at 1 bar for a mixture of 50 mol% PbO and 50 mol% B 2 O 3 over the temperature range [300, 1500 K], where C VM (x, T) was calculated using the finite difference formula RT, and U D M (x, T) was calculated using NVT Monte Carlo simulations with 25,000 equilibration cycles and 75,000 production cycles. Also, the degrees of freedom, DoF , for the ideal gas contribution to the total internal energy were DoF PbO = 9 because PbO is a triatomic linear molecule and DoF B2O3 = 15 since B 2 O 3 is a polyatomic nonlinear molecule.
The shape of C VM (x, T) vs. temperature for the equimolar mixture of PbO and B 2 O 3 is qualitatively the same as that for pure boron trioxide but is shifted to higher temperatures due to the presence of the network former/ modifier PbO. Also, the peak in the heat capacity of the mixture is less than that of pure B 2 O 3 by about half because the heat capacity of lead monoxide is small (0.2 J/g K) compared to that of pure boron trioxide (2 J/g K).  Table 3, p. 591 in Ref. 66 , there is no volume or density data. There is also limited experimental values of density and linear thermal expansion coefficients measured by gamma densitometry in Table 2, p. 11 in Ref. 67 for three different molar compositions of PbO and B 2 O 3 . Moreover, the data in Ref. 66 and in Ref. 67 are easily reconciled and provide consistent values of linear thermal expansion coefficients. Here the focus for comparisons of modeling results with experimental data will be an equimolar mixture of PbO and B 2 O 3 .

Glass transition and softening temperatures of PbO
From the calculated glass transition temperature, T g = 672.10K , and the linear mixing rules for molecular co-volume and internal energy of departure given by Eqs. (6) and (7) Table 5, in which the volumetric coefficients of thermal expansion were computed analytically from 25 to 300 °C and averaged.
The results in Table 5 show that the moving boundary equation of state gives good estimates of specific volume and average volumetric coefficients of expansion for a 50-50 mol% mixture of PbO and B 2 O 3 . The error in the volume prediction at 800 °C is 4.77% while the average value of the volumetric coefficient of expansion, while higher than the experimental values, is certainly the right order of magnitude.  Table 1, p. 369 in Ref. 64 and 658 K reported in Fig. 1, p. 749 in Ref. 65 . In contrast, a glass transition temperature of T g = 678K for an equimolar mixture of PbO and B 2 O 3 is reported in Ref. 58  For pure boron trioxide, the equation of state predictions of volume was within 0.64% and volumetric coefficients of thermal expansion were in the same range as those reported in the literature. Equation of state predictions of the same thermophysical properties for an equimolar mixture of lead monoxide and boron trioxide were also good but comparisons were limited to very few available experimental data. It is also important for the reader to understand that there is limited experimental data for metal oxide glasses in the open literature and, in many cases, there is considerable disagreement among the experimental data. Moreover, because none of the papers used for comparison in our manuscript presents error bars for the experimental data or gives any statistical information regarding experimental errors and with such wide variation in the experimental data for various properties, it is difficult to identify the source of the disagreement between experiments and modeling. In our opinion, what is important from a modeling and simulation perspective is for the modeling framework to give results that (1) fall within the range of published experimental data (i.e., have the right order of magnitude) and (2) show the correct temperature trends. However, this is not to imply that there are not errors introduced in modeling. While computations should not be expected to match experimental data for different properties to the same level of accuracy, the proposed computational framework can introduce errors due to errors in the (1) Monte Carlo simulations of internal energies of departure, and (2) the mixing rules for mixture internal energies of departure and molecular co-volume.
In conclusion, it was shown that glass transition and softening temperatures can be estimated by computing adjacent saddle points on C V (T) vs. temperature, where C V (T) is computed from finite difference derivatives of the total internal energy. Internal energies of departure from Monte Carlo simulations needed to compute C V (T) were used within the multiscale GHC equation of state framework and provided good estimates of volume and volumetric coefficients of thermal expansion.
Finally, the modeling and simulation framework presented in this work can be readily extended to other metal oxide glasses and has the potential to guide the experimental design of these materials.

Data availability
Input and output data from all numerical simulations and plots presented in this paper are available by contacting the corresponding author.