Machine learning approaches for structural and thermodynamic properties of a Lennard-Jones fluid

The principal goal of statistical mechanics is to predict a system’s collective behavior from its constituent interactions. In condensed-phase molecular systems, solving this problem is typically intractable analytically, and therefore, computer-assisted methods must be employed.^1,2 Over the last 60 years, molecular simulations, particularly with respect to fluids, have significantly advanced our understanding of how macroscale observables arise from interatomic interactions.^3–8 A major limitation of using computer simulations, however, is that they are typically performed in an ad hoc fashion and as such the predictive power of these models is limited. This is because the observations taken from computer simulations are virtual measurements of a system’s properties under a particular set of physical conditions, and these measurements are, in general, not transferable to other physical conditions. In addition to the lack of transferability in the solutions produced by molecular simulations, these methods can also incur significant computational costs in order to make simple predictions. This is particularly true when the interatomic interactions are treated quantum mechanically, the system is very large, and/or the system must be simulated over a long timescale in order to make meaningful observations.^8–15 

Recently, new machine learning (ML) methods have been developed and applied to solve problems in physics and chemistry that were previously intractable using traditional molecular simulation methods.^16–20 The success of these ML methods has illustrated that many problems in the physical sciences can be solved faster and more accurately using data-driven approaches based on information analysis than through the application of protocols that institute some set of physics principles directly.^17,21–23 Machine learning methods are not without significant shortcomings, however. For example, if not trained and tested on a suitably general dataset, ML programs often fail catastrophically when asked to solve problems that require extrapolation outside of the dataset on which they are trained.^24–29 Despite having usability that is typically constrained to a narrow application window, ML approaches have been successfully applied to solve both new problems in physics and to develop improved solutions to old problems. The latter is the focus of this article in the context of a canonical problem in statistical mechanics—predicting the structural and thermodynamic properties of a Lennard-Jones (LJ) fluid.^30–37 

Previous attempts to model structural properties of the LJ fluid, specifically, the radial distribution function (RDF), have relied on fitting simulation data to empirically motivated functional forms^33–45 or using integral equation methods.^46–53 Classical density functional theories and quasi-continuum theories have also been developed to determine properties of simple fluids.⁵⁴ There are a number of empirical expressions in the literature that can be applied to produce the RDF of a LJ fluid.^33–37 The pioneering expression of Goldman³³ gives accurate predictions for the RDF in certain temperature and density regimes but has a limited range of applicability. The expression by Morsali et al.³⁵ gives accurate predictions for the RDF over a selected range of temperatures and densities, but it does suffer from the significant problem that its functional form is not continuous and therefore results in unphysical predictions. It also fails under low density conditions. Lack of continuity and a limited range of applicability are also problems in the expressions proposed by Bamdad et al. and Matteoli and Mansoori.^34,36 Other than fitting empirical functional forms, another common approach to predict the RDF of simple fluids is to use integral equation methods. However, integral equations have historically not produced acceptable results for the LJ fluid.^46–48 This is, in part, due to the fact that in order to produce accurate predictions, the functional form of the integral equation, i.e., the closure relation, must be modified specifically to treat the LJ system.^49–51,53

Machine learning approaches have been applied to examine structural correlations in a LJ fluid. Moradzadeh and Aluru have trained a neural network to predict what values of the LJ parameters will give a specific form of the RDF. They then applied this deep learning approach to address inverse design problems for coarse-graining applications.¹⁹ They have also trained an autoencoder network to determine the RDF of simple fluids from a limited number of atomistic configurations.⁵⁵ We have previously developed a ML method based on segmented linear regression and multivariate function decomposition that is able to predict the RDF of simple fluids with significantly increased accuracy in comparison to traditional theoretical approaches.²⁰ 

Determining the equation of state (EOS) of a LJ fluid is another long-standing problem, which has attracted significant interest.^{37,38,41–45,56} Although not an exhaustive list, prominent LJ equations have been developed by Nicolas et al.,³⁸ Johnson et al.,⁴¹ Mecke et al.,⁴² and Thol et al.,⁴⁴ among others. In each of these previous studies, the equation of state is constructed by fitting simulation data to an empirically motivated functional form. Comprehensive discussions of the advantages and disadvantages of many LJ equations of state can be found in Refs. 44 and 57. Recently, ML approaches have been applied to predict thermodynamic properties of fluid systems.⁵⁸ To our best knowledge, however, there have been no attempts to understand how ML can be applied to determine the thermodynamic properties of the LJ fluid and if these data-driven approaches provide any advantages compared to standard frameworks.

In this article, we illustrate how ML methods can be applied to determine properties of a LJ fluid with significantly increased accuracy in comparison to traditional human-guided fitting models and integral equation approaches. First, we provide a detailed analysis of the error decreases that can be expected in comparison to traditional theoretical methods when a modified version of the ML methodology developed in Ref. 20 is used to generate the RDF of a LJ fluid. Second, we use the learned RDF to determine the pressure and internal energy of a LJ fluid and compare those values to the values generated by several analytical functions. Finally, we train a kernel ridge regression (KRR) process^59–61 to determine the pressure and internal energy of a LJ fluid. Comparisons are made between the thermodynamic properties generated by each ML method and the values produced by two important EOS expressions for the LJ fluid. We find that the presented ML methods are able to produce properties of a LJ system more accurately in comparison to previous theoretical methods.

In general, the computational advantage of applying ML approaches to generate structural and thermodynamic properties of condensed-phase systems will be particularly significant when batches of data are needed for applications such as coarse-grained model development, analytical function fitting, inverse design approaches, and interfacing between scales in multi-scale physics applications. The most computationally expensive step in the implementation of the ML methods is the generation of training data, which is an upfront investment of computational effort. We have previously shown, however, that similar ML approaches to those applied here can produce significant error reductions in comparison to traditional theoretical methods even when trained on sparse amounts of data.²⁰ After training, the developed procedures produce structural and thermodynamic properties at a negligible computational cost in comparison to using molecular simulations.

The rest of this article is organized as follows: Sec. II A contains the details of the datasets that are used to train the different ML procedures. The technical details of the ML methods are presented in Sec. II B. A discussion of the testing datasets that are used to quantify the error generated by the developed methods is given in Sec. II C. The results of applying ML to predict the RDF, pressure, and internal energy of a LJ fluid are presented in Sec. III. Comparisons are made between the developed ML procedures, the results from integral equation methods, and several analytical models whose parameters are obtained by fitting simulation data to empirical functional forms. Concluding remarks and thoughts on the ramifications of this work are presented in Sec. IV.

We use both home-generated and literature datasets to train our ML models. The data used to train the ML process to produce the LJ RDF was generated using MD simulations of a system of N = 2500 particles interacting through the cut and shifted LJ potential,

with

where σ and ϵ are the standard LJ parameters and r_c is the cutoff distance. The values of the LJ parameters and the masses of the particles were chosen to correspond to argon. A cutoff of r_c = 4.0σ was used, which is a typical cutoff chosen in LJ fluid studies, e.g., by Johnson et al. in Ref. 41. The simulations were performed in the NVT ensemble using the modified impulse Langevin integrator developed by Goga et al. in Ref. 62. A time step of Δt = 1 fs and a relaxation rate of γ = 0.5 ps⁻¹ were used in the Langevin integrator. A neighbor list that was updated every ten time steps was employed to accelerate the sampling time. Standard periodic boundary conditions were used. The propagation of each MD trajectory was performed in three stages: In the initial equilibration stage lasting 5 ps, linear velocity scaling of each particle was used at every time step in order to maintain a constant temperature. In the second phase, which was the primary equilibration phase, the system was propagated for 500 ps using the modified impulse Langevin integrator with no sampling taking place. Finally, in the third phase, the system was propagated and sampled over 500 ps. During this phase, the RDF, pressure, and energy were measured.

Using the results of these simulations, we calculated the RDF g(r^*, ρ^*, T^*) of the system as a function of the reduced LJ parameter’s distance r^* = r/σ, density ρ^* = ρσ³, and temperature T^* = k_BT/ϵ. The RDF was computed using a histogram bin width of 0.01 Å/σ ≈ 0.002 94. Simulations were performed at 629 state points on a grid in the ρ^* × T^* plane. The grid spacing was Δρ^* = 0.05 and ΔT^* = 0.2 in the respective dimensions. These points were mostly constrained to the vapor, liquid, and supercritical regions of the phase diagram, as shown in Fig. 1, although some of the sampled points were outside of the fluid region. The highest temperature and density sampled were T^* = 6 and ρ^* = 1.25, respectively.

We also trained KRR processes to determine the reduced pressure P^* = Pσ³/ϵ and reduced residual internal energy U^* = U/Nϵ of a LJ fluid. The data used to train these processes were generated by Gottschalk using Monte Carlo simulations.⁵⁶ The Gottschalk dataset contains thermodynamic properties of the LJ system from 8374 state points with temperatures T^* ≤ 6.4. Most of these state points are confined to the region shown in blue in Fig. 1. This dataset was used for training because it consisted of a large number of datapoints and it also contains a high density of points in the supercritical region of the phase diagram, a region that is often not well-sampled in other datasets. An extensive review and analysis of various other datasets containing thermodynamic properties of the LJ system can be found in Ref. 57.

The ML method we applied to determine the LJ RDF is a variant of the linear regression (LR) with the multivariate function decomposition approach developed in Ref. 20. Here, we term this approach the LR-RDF method. The MD data used to train this process consisted of feature vectors of the form $rk*,ρi*,Ti*$⁠, which were mapped to labels $g(rk*,ρi*,Ti*)$⁠, where $rk*$ is a particular value of the variable r^* and $pi=ρi*,Ti*$ is a state point in the training data.

The implementation algorithm for this method is as follows: First, take an input state point $p=ρ*,T*$ (the point where the RDF would like to be predicted) and determine the N_neigh = 4 nearest neighbor points in the training data to the input point. The metric used to determine the distance between input point p and training point p_i is the weighted Euclidean distance

with weights $w$₁ = 1 and $w$₂ = 0.035, which were obtained using a grid search hyperparameter optimization. Next, take the discretized training data in the r^* variable, and at each particular value, $r*=rk*$ construct a linear approximation to the RDF,

where ĝ_k is a linear function in ρ^* and T^* with coefficients that are determined from multivariate least squares regression. The linear approximation is based on the assumption that in the local region in ρ^* × T^* space around input point p, the RDF will vary approximately linearly in ρ^*and T^* at each $rk*$⁠. Finally, combine the collection of regression processes at each $rk*$ to construct the total RDF at point p,

where $I$ is an interpolant over the variable r^* between functions in the set G = {ĝ₁, ĝ₂, …, ĝ_K} and K is the index of the maximum r^* value used in the regression. Structural correlations in the LJ fluid occur over length scales in the dimensional r variable commensurate with σ. Therefore, if the spacing of the training data in the r^* = r/σ dimension is ≪1, the choice of the functional form of the interpolant $I$ (linear, polynomial, etc.) will not be significant. Moreover, for most numerical applications, knowledge of the discretized set of functions G will be sufficient, and therefore, interpolation over r^* will not be necessary. Because the core regression process used in the LR-RDF method is two-dimensional linear regression that is fit to four neighboring state points, there is no significant overfitting in the procedure.

We developed an equation of state for the LJ fluid by training KRR processes to predict the thermodynamic properties pressure P^* and internal energy U^*. The training data for these processes consisted of feature vectors (state points) $pi=ρi*,Ti*$⁠, which were mapped to corresponding labels $P*(ρi*,Ti*)$ and $U*(ρi*,Ti*)$⁠. The data used for training was taken from the Gottschalk dataset.⁵⁶ This dataset is well-suited to use for training because of its size and because the cutoff distance in the simulations used to generate the data is large. A proximity-based approach was employed in the KRR training procedures, meaning that for a particular input point p, only local training data was used to determine the pressure and internal energy at that point. Specifically, the training data used to generate P^* and U^* at point p consisted of N_neigh = 16 nearest neighbor points in the training data as quantified through the distance metric in Eq. (3). The number of neighboring points used in the regression process was found through hyperparameter optimization.

The KRR procedure was implemented using the scikit-learn ML package for Python.⁶³ A polynomial kernel

of degree d = 4 was employed to quantify the similarity between state points in the training data. For each regression process, we performed optimization of the c₀ parameter in Eq. (6) using a grid search where the optimization scoring function was the absolute percent error of the corresponding thermodynamic property P^* or U^*. The conditioning factor α in the KRR process (see Ref. 63 and associated documentation) was set to unity. In each regression, we used a value of γ = 100 in Eq. (6), which was found through a grid search. We took γ to be constant because we found that its optimized value would typically not deviate significantly from the given constant value and because doing so decreases the computational time needed to perform the KRR. In comparison, we found that the optimal value of c₀ would change significantly for each regression process, and therefore, a grid search optimization of this parameter was performed during each KRR procedure.

The test data for the ML procedures were obtained using the same simulation protocol described in Sec. II A. We generated two different test sets. The first test set was generated by performing MD simulations at 300 random points uniformly distributed over the region $0.05,1.15×1.0,5.0$ in the ρ^* × T^* plane and measuring the RDF, pressure, and internal energy at each point. In order to account for errors incurred by using the cut and shifted potential in Eq. (1), standard tail corrections were applied to the measured pressure and internal energy values.^1,41 Any points outside of the vapor, fluid, or supercritical regions were removed from the test set, leaving 265 points. The second test set was generated by performing simulations along select isodensity contours (ρ^* = const.) using a spacing of ΔT^* = 0.2 in the T^* dimension. For each contour, the lowest temperature sampled was dictated by the fluid phase boundaries shown in Fig. 1, and the highest temperature sampled was T^* = 5.

Additionally, there are multiple datasets in the literature containing thermodynamic properties of the LJ system at various state points measured using molecular simulation methods. Here, we used the datasets generated by Johnson et al., Stephan et al., and Meier as test sets to further gauge the accuracy of the developed ML methods.^41,57,64 In each of these datasets, we removed any outlier points as identified in the work of Ref. 57, any points that were outside of the fluid region of the phase diagram (the blue region in Fig. 1), and any points with temperature T^* > 6. After this procedure, the trimmed Johnson, Stephan, and Meier datasets consisted of 134, 317, and 269 points, respectively. We chose to use these datasets as test sets because they are some of the largest datasets available for the LJ fluid and also the cutoff distance used in each set is large (r_c ≥ 4.0σ).

Most previous work developing predictive models for the RDF of the LJ fluid has focused on fitting the results of simulations to empirically selected analytical forms or applying integral equation approaches. The main advantage of using empirical models is that they allow fast and often accurate approximation of the RDF of a LJ fluid without the need to perform computationally taxing simulations or numerical procedures. The ML methodology we apply here can produce LJ RDFs at a negligible computational cost, but with improved accuracy in comparison to human-guided procedures and integral equations.

Comparisons between RDFs measured in simulation and RDFs predicted by the ML procedure described in Sec. II B are shown in Fig. 2 for various temperature and density values. Excellent agreement is observed between the ML and MD results. The results for high-density systems with strong structural correlations are shown in Figs. 2(a) and 2(b). In these systems, the ML protocol almost exactly captures the highly structured form of the RDF. As shown in Figs. 2(c) and 2(d), the ML method also generates accurate predictions for the RDF in intermediate- and low-density systems. This is notable because, as discussed in detail below, previous analytical expressions for the LJ RDF typically fail to give accurate predictions at low densities. In all cases, the ML results are in excellent agreement with the simulation results, illustrating the power of using ML to generate structural properties of fluids.

Figure 3 illustrates the root mean square deviation (RMSD) of the RDF,

predicted using various theoretical methods g_theory taken with respect to the RDF measured in simulation g_MD along several isodensity contours. The RMSD was computed at K_c = 1362 points along the variable r^* using an equidistant spacing of ∼0.002 94 between points. This corresponded to a maximum value of r^* ≈ 4 used in the RMSD calculation. We truncated the RMSD calculation at this value to avoid sampling the structureless regime r^* → ∞ where g(r^*, ρ^*, T^*) = 1 as sampling this regime would bias the RMSD toward a lower value. The lower bound for T^* along each contour is dictated by the fluid phase boundaries shown in Fig. 1.

Along the isodensity contour ρ^* = 1.11, shown in Fig. 3(a), the RMSD of the Morsali expression is ∼0.025 while the RMSD of the Goldman expression is ∼0.02. In contrast, the RMSD generated using ML varies from 0.0005 to 0.001. We also tested if approximating the RDF using the data from the nearest neighbor point in the training set yields an accurate estimate for the RDF. The nearest neighbor approximation yields a RMSD of ≈0.01. This illustrates that the accuracy of the ML method arises from the fitting procedure and is not simply a function of the density of training data. The ML approach also gives more accurate results with respect to well-known integral equations. Specifically, the ML RMSD is greater than an order of magnitude lower than solutions to the Percus–Yevick (PY) and hypernetted-chain equations (HNC). The solutions to the integral equations were obtained using the pyPRISM program.⁶⁵ 

The results along the contour ρ^* = 0.825 are shown in Fig. 3(b). The Morsali expression results in a RMSD of ∼0.01, while the Goldman expression produces an RMSD that is ∼0.015 for low T^* but diverges as T^* is increased. The RMSDs produced by the PY and HNC solutions are large for low T^* but decrease to ∼0.015 and 0.04, respectively, as the temperature is increased. The ML result ranges between ∼0.0005 to 0.001, greater than an order of magnitude decrease when compared to the RMSDs obtained using other methods.

Figure 3(c) illustrates the results for ρ^* = 0.6, which was a density included in the training set. We have used this density to illustrate the power of ML in comparison to other methods at a commonly examined density, not to infer any accuracy metrics about the ML procedure itself. The RMSD due to the Morsali expression is ∼0.005. The Goldman expression produces an RMSD of ∼0.015 for T^* < 4 but diverges for T^* > 4. The PY and HNC solutions generate errors that are, in general, greater than the other methods. In comparison, the RMSD due to ML is ∼0.001—a significant error reduction.

The results along the contour ρ^* = 0.47 are shown in Fig. 3(d). The Morsali RMSD is ∼0.02 at low temperatures and decreases to ∼0.004 at high temperatures. The Goldman expression again produces an RMSD result of ∼0.015 for low values of T^* but quickly diverges as T^* is increased. Away from the low-temperature regime, the integral equations result in errors of the order 0.005 and 0.012 for the PY and HNC solutions, respectively. The ML RMSD is ∼0.004 at low temperatures but then decreases quickly to ∼0.0007 as T^* is increased. Typically, we found that using ML to predict the LJ RDF resulted in greater than an order of magnitude decrease in the RMSD in comparison to traditional theoretical methods in high- and intermediate-density regimes.

The RMSD calculated along the isodensity contours ρ^* = 0.25 and ρ^* = 0.125, which both correspond to the low-density regime of the LJ fluid, are shown in Figs. 4(a) and 4(b) respectively. Notice that the RMSD for each density is shown on a log scale. In both cases, the ML protocol developed in this article significantly outperforms both the Morsali and Goldman expressions. In Fig. 4(a), the Morsali RMSD is of the order 1 for low T^* but quickly decreases to 0.1 as T^* is increased, while the Goldman RMSD is of the order 0.1 at lower temperatures but diverges to values greater than 100 for T^* > 4. In comparison, the RMSD calculated using ML is ∼0.002 at low T^* but decreases to ∼0.0005 as the temperature increases. For the case of ρ^* = 0.125, shown in Fig. 4(b), the Morsali expression fails catastrophically and produces an RMSD that is, at best, of the order 10, while the Goldman RMSD is again ∼0.1 for T^* < 4 but diverges for T^* > 4. The failure of the Morsali expression arises because it is a piecewise noncontinuous function, and one of the branches of the expression diverges at low density creating the large deviations from the MD results, using the ML results in an RMSD of ∼0.002 at this density.

In contrast to the empirically motivated functional forms, both the PY and HNC integral equations perform reasonably well at low density. For ρ^* = 0.25, shown in Fig. 4(a), the PY and HNC equations produce RMSDs that are ∼2 to 5 times greater than ML. Therefore, the integral equation RMSDs are large in comparison to the ML method but provide a dramatic improvement over both the Goldman and Morsali expressions. For ρ^* = 0.125, shown in Fig. 4(b), the RMSD produced by both integral equations is approximately equal to that produced using ML. Also shown in Fig. 4(b) is the result given by taking the theoretical RDF to correspond to the infinite dilution limit $g(r*,T*)=e−ϕ*(r*)/T*$⁠, with $ϕ*(r*)/T*=ϕ(r)/kBT=4/T*[(1/r*)12−(1/r*)6]$⁠. Infinite dilution corresponds to the limit in which multi-body effects vanish. The RMSD given by infinite dilution approximation is approximately an order of magnitude greater than the RMSD generated by ML. This illustrates that theoretical procedures such as ML or integral equation methods must be employed to obtain accurate estimates for the RDF in systems with significant multi-body structural correlations. In general, we find that using ML to predict the RDF in low-density LJ systems will typically provide two orders of magnitude or greater decrease in the RMSD when compared to empirical functional forms and will often improve on the results of integral equation solutions.

Shown in Table I is the RMSD generated by each theoretical method calculated over the random test set described in Sec. II C. The first column lists the various theoretical methods. The second column is the RMSD result for each method calculated over every point in the test set. The ML approach reduces the error by factors of ∼20 and 30 in comparison to the PY and HNC equations, respectively. The results are not shown for the Morsali and Goldman expressions because the claimed range of validity of both expressions is ρ^* > 0.35, and the lower density bound in the random test set is ρ^* = 0.05. We therefore cannot assess the error generated by these expressions using the full test set because both expressions return very large and often numerically divergent RMSD values in the low density regime (see Fig. 4). In order to construct a test set to compare the theoretical methods away from the low density regime, we removed any points from the full set with ρ^* < 0.35 and calculated the RMSD generated by each method over this trimmed dataset. These results are shown in the third column of the table. Over the trimmed test set, ML again generates the lowest error. The error from the Morsali expression is the next lowest, roughly 15 times greater than ML. The ML approach reduces the error by respective factors of ∼40 and 60 in comparison to the PY and HNC equations. Note that the RMSD generated by each integral equation increases for the trimmed test set in comparison to the full test set because integral equations perform best in the low-density regime, and therefore, removing the low-density state points increases the average RMSD. The Goldman expression returns divergent RMSD values in the high-temperature regime, as shown in Fig. 3. Therefore, to calculate the average RMSD for the Goldman expression, we also removed any state points from the test set that generated an RMSD > 1. The RMSD value generated by ML over the trimmed test set is approximately a factor 60 lower than RMSD calculated using the Goldman expression with high-temperature state points removed. The RDF can also be estimated using data from the nearest point in the training data set. Over both the full and trimmed test sets, the ML method reduces the RMSD by approximately an order of magnitude in comparison to using the RDF from the nearest neighbor to estimate the RDF at the input point.

Shown in Fig. 5 is the learning curve for the LR-RDF method illustrating how the accuracy of the method varies as a function of the number of points in the training set. When using as few as 100–200 training points, ML reduces the RMSD by approximately an order of magnitude, or more, in comparison to each of the other methods. Moreover, the LR-RDF method outperforms all the other methods when trained on as few as ∼30 points. Approximating the RDF at the input point using data from the nearest neighbor point in the training set yields better results than all the other methods except ML when the number of training points, i.e., the density of training data, is large but yields poor results in comparison with the other methods as the training data becomes sparse. It should be noted that as the density of the training data decreases, a threshold is reached where the ML method ceases to be more accurate than other theoretical methods. This threshold is reached, however, when there are only a few points in the training set.

It is important to note that the different analytical expressions developed to generate the LJ RDF have been fit to data that were obtained using many different simulation protocols.^33,35,36,38 Because of these discrepancies between simulation procedures, we can only make observations about how well any previously developed functions are able to reproduce the data obtained under the specific simulation conditions used here. However, when the results of these prior analytical expressions are compared to our simulation data, we observe RMSD values that are typically similar to or less than those reported in the original articles. This implies that the simulation data obtained in this manuscript is well-representative of the data used to fit these functions. Moreover, based on the observed accuracy of the developed ML methods, we expect that ML will provide substantial improvements when compared to all of the previous empirical functions over a broad range of simulation conditions.

The major bottleneck of applying the present ML method is generating training data in a way that effectively samples the fluid region of the phase diagram while also being computationally efficient. While the comparisons given here illustrate the improvements in predictive accuracy that can be expected when using ML to generate the RDF of a simple LJ fluid (see also the supplementary material of Ref. 20), questions remain with respect to the best way to generate training data and if a grid-based approach will scale well when applying the developed ML method to determine structural properties of more complex systems such as fluid mixtures. These questions will be addressed in our future work.

The ability to accurately predict macroscopic thermodynamic quantities from the function used to generate the RDF is an important test of both the precision of the fitting procedure and the validity of the chosen functional form used in the fitting.^33–37 In simple fluids that are dominated by pairwise interactions, structurally based macroscopic observables can be expressed using standard thermodynamic relations containing the RDF. In the specific case of a LJ fluid, the reduced pressure P^*(ρ^*, T^*) and reduced internal energy U^*(ρ^*, T^*) are related to the RDF through the respective expressions¹ 

and

Therefore, once a functional form for the RDF is known, it provides a direct route to generate the thermodynamic properties of the LJ system.

To further gauge the accuracy of the developed ML protocol, we evaluated Eqs. (8) and (9) numerically using the learned RDF and compared those results to pressure and internal energy values that were measured directly in MD simulations. Calculating macroscopic observables using the predicted RDF and comparing those values to values measured in simulation is a standard method that is typically used to test the accuracy of the specific procedure that is applied to generate the LJ RDF. We applied a long-range correction to the ML RDF by taking g(r^*, ρ^*, T^*) = 1 for values of r^* > L/2, where L is the simulation box length in the MD training data.

Figure 6 illustrates a comparison between the results for P^* and U^* measured directly in simulation and the results calculated by using the ML RDFs in Eqs. (8) and (9). The pressure and internal energy results are shown in Figs. 6(a) and 6(b), respectively. In both figures, the results are shown as a function of T^* along various isodensity contours. The agreement between the two methods is excellent in all cases.

Shown in Fig. 7 is the percent error,

for each thermodynamic property $X∈P*,U*$ calculated along various isodensity contours using the LR-RDF method as well as the Goldman and Morsali expressions. In Eq. (10), X_theory is the value predicted by the particular theoretical method, and X_MD is the value taken from MD simulation. Figures 7(a) and 7(b) show $Perror*$ along the contours ρ^* = 1.021 and ρ^* = 0.88, respectively. The ML procedure generates significantly lower error than the other methods at these high densities. It is interesting to note that the error generated by the Goldman and Morsali expressions is quasi-oscillatory, while the error due to the ML approach presents as noise. As shown in Fig. 7(c), the error reduction produced by using the LR-RDF method is particularly pronounced in the low-density regime. Specifically, for ρ^* = 0.25, the error produced by the Goldman expression is very large, while the Morsali expression produces an error of ≈10%. In contrast, ML produces an error of ≈0.06%. The percent error in the internal energy $Uerror*$ along the contours ρ^* = 1.021 and ρ^* = 0.88 are, respectively, shown in Figs. 7(d) and 7(e). The LR-RDF method outperforms the other expressions, particularly at high temperatures. Figure 7(f) shows the internal energy error for ρ^* = 0.25 where, again, applying ML produces significantly lower error than the other methods.

In general, the absolute error |X_error(ρ^*, T^*)| in each thermodynamic property calculated using the LR-RDF method is typically in the range 0.05%–0.5% across all density regimes. In contrast, both the Morsali and Goldman expressions result in average absolute errors of ∼2% at high densities and greater than 10% at low densities. These results highlight the effectiveness of applying ML to predict a system’s macroscopic properties using the learned RDF.

We also trained a KRR process to produce thermodynamic properties of the LJ system as described in Sec. II B 2. This process is trained directly on thermodynamic properties measured in simulation and therefore does not include information from the RDF. The black dashed lines in Fig. 6 are the results for P^* and U^* generated using KRR along various isodensity contours. In all cases, the agreement between the MD results and the results generated using ML is excellent. In fact, upon detailed examination, it can be observed that the KRR process yields results that are in better agreement with the MD values than those generated using the LR-RDF method.

Figure 8 illustrates a comparison between the percent error in the pressure and internal energy values generated using the KRR method, the LR-RDF method, the Johnson EOS, and the Mecke EOS. The results are shown along various isodensity contours. The errors in the predicted pressure for high-density systems with ρ^* = 1.021 and ρ^* = 0.88 are shown in Figs. 8(a) and 8(b), respectively. The KRR method is the most accurate in both cases, generating an error along each contour that fluctuates around zero. The errors generated by all of the theoretical methods are correlated because the same MD data are used to compute the error for each method. The MD data fluctuates along each isodensity contour due to sampling errors, and so calculating the percent error generated by each theoretical equation of state produces correlated fluctuations. The results for a low-density system with ρ^* = 0.25 are shown in Fig. 8(c). In this case, both the KRR and LR-RDF methods generate errors along the contour that fluctuate about zero, while the Mecke and Johnson expressions yield respective errors of ∼0.25% and 0.5%. The percent error in the internal energy produced by each method for ρ^* = 1.021 and ρ^* = 0.88 is shown in Figs. 8(d) and 8(e). At these high densities, the Johnson EOS is clearly not as accurate as the other methods and predicts values for the internal energy that varies significantly from the MD results. The Mecke EOS and the LR-RDF method yield similar errors at both densities, while the KRR method is again the most accurate. The error in the predicted internal energy for ρ^* = 0.25 is shown in Fig. 8(f). The Mecke and Johnson expressions produce errors that change sign as T^* is increased. Therefore, while the average percent error along the contour is approximately zero, the absolute percent error generated by each expression will be relatively large. The LR-RDF method yields slightly better results than the two analytical EOSs, and the KRR process is, again, the most accurate method. The KRR method and the Mecke and Johnson expressions are fit to two-dimensional input data while the LR-RDF method uses three-dimensional inputs. It is therefore noteworthy that the LR-RDF procedure produces comparable results to the others while being trained on more complex data structures.

Table II shows the results of using various methods to calculate the average absolute percent error in the pressure $|Perror*|$ over four different test sets. The specific methods compared are the KRR and LR-RDF ML methods, the Johnson and Mecke EOS expressions, and the result obtained by approximating the pressure at the input point using the pressure value from the nearest neighbor point in the training set. The examined test sets are the random test set generated here and the datasets of Johnson et al., Stephan et al., and Meier.^41,57,64 Descriptions of each test set can be found in Sec. II C. The combined error averaged over all of the test sets is also shown in the last column of the table. The Mecke expression is widely considered to be the most accurate LJ EOS for generating thermodynamic properties. The KRR process results in significant error decreases in comparison to the other methods, outperforming the Mecke EOS over every test set except the Johnson dataset where the methods yield similar results. Specifically, the KRR procedure reduces the error when compared to the Mecke EOS over the combined test set by a factor of ≈1.7. The LR-RDF method is approximately as accurate as the Johnson EOS and, in general, results in errors in the pressure that are between 2 and 10 times greater than that produced by the KRR method and the Mecke EOS. Because the KRR process is more accurate than the Mecke EOS, we conclude that the KRR process is currently the most accurate method for generating the pressure of the LJ system. The nearest neighbor approximation leads to significant errors, typically close to an order of magnitude greater than the other methods.

The average absolute percent error in the internal energy $|Uerror*|$ generated by each theoretical method calculated over each test set is shown in Table III. In general, the LR-RDF method is more accurate than the Johnson EOS, but not as accurate as the Mecke EOS, while the KRR method reduces the absolute error in comparison to the Mecke EOS. Specifically, the developed KRR approach outperforms the Mecke EOS on every test set except for Johnson et al. data. When the absolute error is averaged over every point in all of the test sets combined, as shown in the last column of the table, the KRR method again reduces the error in the comparison to the Mecke EOS by roughly a factor of 1.7. We therefore conclude that the KRR approach developed here is currently the most accurate theoretical method for determining the internal energy of the LJ system. The nearest neighbor approximation is overall the worst method; however, over the Johnson et al. test data, this method generates more accurate predictions for the pressure than the LR-RDF method and the Johnson EOS and more accurate predictions for the internal energy than the Johnson EOS. This behavior is explained, in part, because the Johnson dataset is the smallest we have examined, and the values in this dataset are sampled from the fewest number of configurations in comparison to the other datasets, leading to uncertainty in the error calculations.

Figure 9 shows learning curves for the KRR process. Specifically, the average absolute percent error in the pressure and internal energy are shown, respectively, in Figs. 9(a) and 9(b) as a function of the number of state points in the training set. The averages are taken over the test set constructed in this manuscript (see Sec. II C). For both of the thermodynamic properties, there are only small increases in the error as the training set is reduced in size from ∼10 000 points to ∼100 points. It is interesting to note that when the KRR process is trained on as few as ≈200 points, it provides more accurate predictions for both the pressure and internal energy in comparison to all of the other methods. The fluctuations in the learning curves arise because the size of the training dataset cannot be reduced in a completely uniform way because some regions of the phase diagram contain a higher density of state points than others in the full Gottschalk training set.

The EOS expressions that have been developed for the LJ system are some of the most accurate due to the importance of the LJ system as the prototypical fluid model and the simplicity of the interatomic LJ potential. It is significant that the ML approaches developed here are more accurate than previous methods/expressions applied to the LJ system, and we therefore expect that in other condensed-phase systems, ML will provide even larger increases in predictive accuracy than those observed here for the LJ fluid. For example, in complex fluid systems, analytical EOS expressions typically do not exist or they provide only a qualitative description of the system’s thermodynamical properties. An advantage of the present ML approaches is that because no assumptions are made about the functional form of the EOS or the RDF, the developed methods are easily transferable to more complex systems. This opens the possibility to apply similar approaches to complex condensed-phase systems that have previously been intractable to treat with traditional methods other than brute-force simulation.

Machine learning protocols have been developed that are able to determine structural and thermodynamic properties of a LJ fluid more accurately in comparison to previous human-guided fitting attempts and integral equation methods. This work illustrates an example case in which ML methods can be used to supersede human efforts in the development of predictive functions that are generated by fitting large datasets. All of the raw training data used in this article as well as the ML programs that are used to compute the RDF, pressure, and internal energy of a LJ fluid have been made publicly available at https://github.com/gtcraven. The primary use of these programs is to avoid the computationally taxing task of performing molecular simulations in order to generate properties of the LJ system. This will be especially useful when a large number of simulations are needed over various parameter values for: interfacing between scales in multi-scale physics codes, analytical function fitting, molecular design, and other applications in which large batches of data are needed over diverse system states.

Applying the developed LR-RDF method was shown to dramatically improve the predictive accuracy of generating the RDF of a LJ fluid when compared to solving integral equations and fitting human-guided functional forms. Typically, we found that using the LR-RDF method reduced the error by an order of magnitude in comparison to traditional theoretical methods. The thermodynamic properties generated using the developed KRR ML method were also able to reduce the predictive accuracy when compared to previous methods, but by a smaller factor, typically of the order 2. The reason that the LR-RDF approach produces a much larger error decrease for structural properties than the KRR approach does for thermodynamic properties is because the dimensionality and complexity of the structural functions is greater, and therefore, previous theoretical methods have historically had difficulty producing accurate results for structural properties of the LJ fluid. Said in more detail, the LJ thermodynamic properties are fit to two-dimensional inputs, and therefore optimization and regression over this data is much simpler than performing regression using the three-dimensional inputs that are used to fit the LJ RDF.

The major disadvantage of applying the presented ML methods is that generating enough training data to sufficiently sample a particular region of a system’s phase diagram can be computational expensive. Moreover, in some systems, this step may be so computationally taxing that it precludes the use of ML to develop predictive models. In fact, it should be noted that in any situation in which (a) data are only needed at a single state point, (b) a simulation methodology exists that can accurately produce properties of the system, and (c) no ML procedure has been previously trained to produce the properties of interest for the system, performing a molecular simulation at this single point will be much more computationally efficient than first generating training data and then implementing ML. As described in detail previously, however, we anticipate that there will be a number of systems and situations in which the present ML approaches will be useful.

We have applied ML to produce particular properties of the LJ fluid; however, the developed protocols should serve as general blueprints that can be used to predict other structural, thermodynamic, transport, and dynamical properties in a multitude of fluid systems, including soft and complex fluids.^66–75 Our future work will address applications of ML in the development of predictive functions for other properties of isotropic fluids. Extending the ML methods developed in this article to the case of anisotropic fluids is an important next step, and work in this direction is currently underway.

The data and computational codes that support the findings of this study are openly available at https://github.com/gtcraven, Refs. 76 and 77.

We acknowledge support from the Los Alamos National Laboratory (LANL) Directed Research and Development funds (LDRD). This research was performed, in part, at the Center for Nonlinear Studies (CNLS) and the Center for Integrated Nanotechnologies (CINT) at LANL. The computing resources used to perform this research were provided by the LANL Institutional Computing Program. We thank Ryan Jadrich and David Rosenberger for insightful discussions.