BoltzmaNN: Predicting effective pair potentials and equations of state using neural networks

TL;DR

This paper employs neural networks to accurately predict equations of state from pair potentials and to derive effective pair potentials from radial distribution functions, improving inverse design and coarse-graining in molecular simulations.

Contribution

It introduces neural network models that outperform traditional methods in predicting equations of state and effective pair potentials, incorporating force information for enhanced accuracy and transferability.

Findings

Neural networks outperform low-density analytic estimates for virial coefficients.

Including force data improves potential predictions, making them smoother and more accurate.

The models are effective for inverse design and coarse-graining applications.

Abstract

Neural networks (NNs) are employed to predict equations of state from a given isotropic pair potential using the virial expansion of the pressure. The NNs are trained with data from molecular dynamics simulations of monoatomic gases and liquids, sampled in the $NVT$ ensemble at various densities. We find that the NNs provide much more accurate results compared to the analytic low-density limit estimate of the second virial coefficient. Further, we design and train NNs for computing (effective) pair potentials from radial pair distribution functions, $g(r)$, a task which is often performed for inverse design and coarse-graining. Providing the NNs with additional information on the forces greatly improves the accuracy of the predictions, since more correlations are taken into account; the predicted potentials become smoother, are significantly closer to the target potentials, and are more transferable as a result.

Tables

Table 1. Table 1: List of abbreviations for the optional parameters of the networks.

Abbreviation	$\overline{U}\left(r\right)$	$\exp \left[-\beta U\left(r\right)\right]$	force	$P_0$
x	✗	✗	✗	✗
c	✓	✗	✗	✗
e	✗	✓	✗	✗
f	✗	✗	✓	✗
i	✗	✗	✗	✓
cf	✓	✗	✓	✗
ci	✓	✗	✗	✓
ef	✗	✓	✓	✗
ei	✗	✓	✗	✓
fi	✗	✗	✓	✓
cfi	✓	✗	✓	✓
efi	✗	✓	✓	✓

Table 2. Table 2: List of abbreviations for the optional parameters of the networks.

Abbreviation	$U_0\left(r\right)$ input	force in loss	$\exp \left[-\beta \hat{U}\left(r\right)\right]$ output
x	✗	✗	✗
b	✓	✗	✗
e	✗	✗	✓
f	✗	✓	✗
be	✓	✗	✓
bf	✓	✓	✗
ef	✗	✓	✓
bef	✓	✓	✓

Equations

$$l_{\mathrm{SE}}(y,\hat{y})$$

$$l_{\mathrm{AE}}(y,\hat{y})$$

$$l_{\mathrm{LE}}(y,\hat{y})$$

$$L_{\mathrm{MSE}}(\mathbf{y},\hat{\mathbf{y}})$$

$$L_U=L_{\mathrm{MLE}}+{\alpha }_{\Delta }{L}_{\Delta }+\sum _{k=1}^4{L}_k{\alpha }^k$$

$$L_{\Delta }=\frac{1}{d-2}\sum _{i=1}^{d-2}{l}_{\mathrm{SE}}\left(y_{i+2}-2y_{i+1}+y_i,{\hat{y}}_{i+2}-2{\hat{y}}_{i+1}+{\hat{y}}_i\right).$$

$$L_k=\frac{1}{d-k}\sum _{i=1}^{d-k}{l}_{\mathrm{SE}}\left(y_{i+k}-y_i,{\hat{y}}_{i+k}-{\hat{y}}_i\right).$$

$$L_k=$$

$$\displaystyle-2\left(y_{i+k}-\hat{y}_{i+k}\right)\left(y_{i}-\hat{y}_{i}\right)\big{]}.$$

$$B_2(T)\approx -2\pi {\int }_0^{\infty }f(r,T)r^2\mathrm{d}r,$$

$$f(r,T)=\exp \left[-\beta U(r)\right]-1.$$

$$\beta P_0=\rho -2\pi {\rho }^2{\int }_0^{\infty }f(r,T)r^2\mathrm{d}r.$$

$$U_0(r):=-k_{\mathrm{B}}T\ln \left[g(r)\right],$$

Full text

1 Introduction

Understanding and predicting the relationship between the (macroscopic) properties of a material and its (microscopic) building blocks is one of the key challenges in materials research and physics. One important goal in statistical physics is the accurate prediction of the phase behavior on the basis of the (effective) pair potential $U(r)$ acting between the particles. According to van der Waals’ theorem of corresponding states, all simple fluids obey the same reduced equation of state (EOS), if the thermodynamic variables are rescaled by their value at the critical point. However, this law only applies for systems with conformal pair potentials, i.e., when the potentials can be fully superimposed by adjusting the interaction strength and particle diameter which is rarely the case in practice. Noro and Frenkel extended this principle by including the reduced second virial coefficient for quantifying the effective range of the attraction.1 This extended approach can provide accurate predictions for pair potentials which are characterized by attractive interactions with ranges much smaller than the particle size,1, 2 but it is expected to fail for more complex pair potentials which, e.g., include a repulsive barrier. Thus, alternative prediction tools are highly desirable, especially given that a large number of (effective) pair potentials in soft matter are bounded or have repulsive barriers.3, 4, 5, 6, 7

Progress in this field has wide implications, not just in terms of our fundamental understanding, but also due to the large number of potential technological applications. Various mechanical, optical, and electronic material properties critically depend on the degree of ordering of their atomic or (macro)molecular constituents. In the conventional forward design approach, the development of new materials typically begins with designing candidate building blocks that are expected to lead to the desired properties. Then these candidates are created, tested, and, if necessary, modified, until the compound with the desired properties has been identified. However, this iterative procedure can quickly become rather time- and resource consuming. Therefore, good initial candidates are required to achieve convergence in a reasonable time frame. Due to these inherent issues, there has been a recent paradigm shift towards the inverse design process, where the building blocks are inferred from the desired target properties. This pathway has been explored for a range of soft materials,8, 9 including athermal granular media,10 colloids,11, 12 and block copolymers.13, 14 In recent years, informatics-driven approaches have gained popularity that utilize machine learning algorithms on large databases to identify previously unrecognized patterns and to predict new candidate materials.15, 16, 17, 18, 19, 20, 21, 22

This inverse design process is strongly related to the problem of top-down coarse-graining, where the goal is to derive effective (pair) potentials for a system with a reduced number of degrees of freedom, while preserving selected target properties of the original fine-grained representation. Several methods have been developed for this task, including Reverse Monte Carlo (RMC),23, 24, 25 Iterative Boltzmann Inversion (IBI),26 and simulated annealing-based optimization.27 These techniques have been employed successfully for, e.g., developing effective pairwise potentials from experimental structure measurements28, 29, 30 and for coarse-graining atomistic simulations.31, 26, 32, 33 Achieving transferability and representability of such coarse-grained models is, however, a key challenge,34, 35 given that the multi-body potential of mean force is usually approximated by an effective pair potential.23, 24, 25, 34, 26, 27, 36, 35 Further, although the IBI scheme should in principle provide a unique solution for a given $g(r)$,37 in practice however, convergence of this iterative procedure is not guaranteed.35 Various strategies have been devised for improving these and other aspects of coarse-graining, such as the addition of thermodynamic constraints,26, 38, 39 and the development of improved methods is an active field of research.40, 41, 42, 43, 44

In this work, we employ artificial neural networks (NNs) for predicting the EOS from a given isotropic pair potential $U(r)$ using the virial expansion of the pressure. Further, we design and train NNs for computing (effective) pair potentials from a given radial pair distribution function, $g(r)$. The training and test data for the NNs are generated from Molecular Dynamics (MD) simulations of monoatomic gases and liquids in the canonical ensemble. The NNs developed in this work as well as example scripts are available online.45 The rest of this manuscript is organized as follows. In Sec. 2, we provide a brief summary of the numerical methods we used as well as how the data were generated. In Sec.  3, we present our results, where Sec. 3.1 focuses on using NNs for determining the EOS from a given pair potential, and Sec. 3.2 discusses how NNs can be used for computing (effective) pair potentials from the radial pair distribution function. Section  4 contains our conclusions and a brief outlook

2 Methodology

2.1 Design of Neural Networks

We used three different NN architectures in this work, i.e., fully connected dense NNs (DN), convolutional NNs (CN), and U-Nets46 (UN). The DNs consist of layers of neurons, where all neurons between two subsequent layers are connected to each other. Operations in DNs are limited to simple matrix-matrix and matrix-vector-multiplications, combined with (nonlinear) activation functions to break linearity. The CNs use convolutional operations which, compared to a DN, drastically reduce the number of parameters necessary by considering only local correlations between points. Therefore, CNs are useful for strongly correlated data, as is typically the case in image processing but also in our problem. In UNs, the input is first processed using convolutional layers then upsampled again and concatenated to a former stage of the processing, hence extracting features and combining them with the original input for further processing.

The accuracy of an NN depends strongly on the information that is provided to it, e.g., the average particle number density and/or the force. Also, the format in which this information is represented plays an important role. For instance, one could use $\exp\left[-\beta U(r)\right]$ rather than $U(r)$, which may be favorable for pair potentials with a strong repulsion combined with an attraction of significantly smaller strength. A detailed study on the effect of these parameters is provided in Sec. 3.

Another key aspect of the NNs is the loss function, which essentially controls the properties that should be optimized by the NN. The loss functions used in this work were combinations of three basic loss functions, the squared error (SE), the absolute error (AE), and the logcosh error (LE). For a given pair of target and prediction, $(y,\hat{y})$, the basic losses are

[TABLE]

The loss of the entire output vector is then determined as the mean over the losses of the individual nodes. For instance, the mean squared error (MSE) is given by

[TABLE]

with $\mathbf{y},\mathbf{\hat{y}}\in\mathbb{R}^{d}$ and dimensionality $d$ of the output vector (e.g. $d=50$ for predicting pair potentials). Analogous expressions are used for $L_{\rm MAE}$ and $L_{\rm MLE}$. We indicate the average loss taken over multiple output vectors through angular brackets, $\left\langle\dots\right\rangle$.

Expressions like Eq. (4) measure the difference between two single points of $\mathbf{y}$ and $\mathbf{\hat{y}}$, but do not capture the correlations between neighboring points. To ensure that the predicted potentials vary smoothly with $r$, it is therefore helpful to include a Laplace-like term as well as loss terms correlating the differences of the potential at different distances

[TABLE]

with weight $\alpha_{\Delta}$ chosen as $\alpha_{\Delta}=2$. The terms $L_{k}$ can optionally be multiplied by a factor $\alpha\leq 1$ to reduce their contribution with increasing distance ($\alpha=1$ was used in this work). The term $L_{\Delta}$ is the discretized Laplace term, which in the case of $l_{\rm SE}$ reads

[TABLE]

The loss terms $L_{k}$ in Eq. (5) are given by the loss function applied to the difference of potential values being $k$ indices apart. For example, with $l_{\rm SE}$ this term reads

[TABLE]

To see how this definition naturally introduces spatial correlations, one can rewrite $L_{k}$ from Eq. (7) as follows

[TABLE]

In this representation of $L_{k}$, it is clear that the last term in Eq. (8) introduces correlations between values at different distances.

All NNs have been constructed and trained using Tensorflow v. 1.13.1.47 The networks were trained using an Adam (adaptive moment estimation) optimizer48 with a learning rate of 0.001. In order to find a suitable NN topology, an extensive grid search was conducted for the three network architectures, where we systematically studied the influence of the loss function as well as the number and shape of the layers, i.e., the width or number of filters and kernel size. During this grid search, the NNs were trained for $10^{4}$ epochs without any optional parameters (see Sec. 3 below). For each architecture, we chose the NN with the lowest $\left\langle L_{\rm MSE}\right\rangle$ and $\left\langle L_{\rm MAE}\right\rangle$ evaluated through 4-fold cross-validation, and trained these NNs then for $2\times 10^{5}$ epochs with optional parameters. For the final benchmarks, the selected NNs were trained for $5\times 10^{5}$ epochs using all data from the training and validation sets.

2.2 Generation of data

In order to generate training and test data for the NNs, we performed a series of MD simulations of monoatomic gases and liquids in the canonical ensemble, and determined the resulting pressures $P$ and radial pair distribution functions $g(r)$. If not stated otherwise explicitly, the particle diameter $\sigma$ and the thermal energy $k_{\rm B}T$ are used as the units of length and energy, respectively. Pair potentials $U(r)$ were generated as spline functions with cutoff distances $r_{\rm cut}\in(0,5\,\sigma]$. The number of base points was chosen randomly in the range $n\in[6,10]$, and the base points were distributed randomly at distances in the range $[0,r_{\rm cut}]$ with magnitudes in the interval $[0,15\,k_{\rm B}T]$. Then a smoothing spline function was applied to connect these points. The degree of the spline function was randomly drawn between $k=2$ and $k=5$ to allow for a large variety of shapes [note that $k\geq 2$ is required to guarantee differentiability of $U(r)$]. To ensure physically meaningful pair potentials, the splines were fitted with the constraint $F(r_{\rm cut})=-\frac{\partial}{\partial r}U(r)|_{r_{\rm cut}}=0$. Further, the potentials were shifted so that $U(r_{\rm cut})=0$. Hard-core interactions were included in selected pair potentials by adding the Weeks-Chandler-Andersen (WCA) potential49 with diameter $\sigma$ and interaction strength $\varepsilon=k_{\rm B}T$.

Using this procedure, 657 potentials were generated, where 344 of them had an additional hard-core contribution. For the MD simulations, tabulated potentials were created with 200 evenly distributed points in the range $[0,r_{\rm cut}]$. For the optimization of the NNs, these tables were then further downsampled to 50 points, evenly distributed in the interval $[0,5\,\sigma]$. Figure 1(a) shows a selection of potentials which were generated using our procedure, while the resulting radial pair distribution functions are plotted in Fig. 1(b,c).

The MD simulations were conducted in a cubic simulation box with edge length $L=40\,\sigma$ and periodic boundary conditions in all directions. The temperature was fixed at $T=1$ using a Langevin thermostat, and the equations of motions were solved using the standard Velocity-Verlet algorithm with a time step of $\Delta t=0.005\tau_{\rm MD}$, with intrinsic MD unit of time $\tau_{\rm MD}=\sqrt{m\sigma^{2}/(k_{B}T)}$. Simulations were conducted at four different particle number densities, i.e., $\rho=0.125\,\sigma^{-3}$ ($N=8000$), $\rho=0.244\,\sigma^{-3}$ ($N=15625$), $\rho=0.512\,\sigma^{-3}$ ($N=32768$), and $\rho=1.0\,\sigma^{-3}$ ($N=64000$). Each simulation was first run for $5\times 10^{6}$ MD steps for equilibration, and then for additional $5\times 10^{6}$ MD steps to sample $P$ and $g(r)$ [see Fig. 1(b,c) for examples]. For the discretization of the calculated $g(r)$ we chose 200 points that were evenly distributed in the interval $[0,L/2]$. In total, 1678 simulations were conducted, with 529 containing potentials with hard-core repulsion. Simulations were discarded from our analysis, which did not reach equilibrium in the allotted time or formed heterogeneous structures [see ESI for more information], leading to a final number of 891 data points. Of those, 790 data points were used as training and validation sets, while the remaining 101 were reserved for testing only. All simulations were performed on GPUs using the HOOMD-blue software package (v. 2.4.2).50, 51

3 Results

3.1 Predicting the equation of state

We applied NNs for predicting the pressure $P$ at a given density from the employed pair potential via three different strategies:

The pair potential is directly mapped to the pressure, $U(r)\mapsto\beta P$; 2. 2.

The pair potential is first mapped to an effective, density-dependent second virial coefficient $U(r)\mapsto B_{2}^{*}$ and then to the pressure, $\beta P=\rho+\rho^{2}B_{2}^{*}$; 3. 3.

The pair potential is first mapped to a set of $n$ virial coefficients $U(r)\mapsto\mathcal{B}_{n}\coloneqq\{B_{i}\;|\;2\leq i\leq 5\}$ and then to the pressure, $\beta P=\rho+\sum_{B_{i}\in\mathcal{B}_{n}}B_{i}\rho^{i}$.

Figure 2 shows schematic representations of the NNs used for predicting $P$ from a given pair potential $U(r)$ and density $\rho$. For the strategies $U(r)\mapsto\beta P$ and $U(r)\mapsto B_{2}^{*}$, the density was included as an additional input after convolutional operations, so that the density-dependence is learned by the NN. The calculation of $P$ for the network calculating $\mathcal{B}_{n}$ is achieved via layers calculating the scalar product or sum as indicated in the schematic. These layers do not contain any adjustable parameters, so that the learned virial coefficients $\mathcal{B}_{n}$ are density independent, while still allowing for the loss function to be applied to the density-dependent pressure $P$.

The only required input parameters of our NNs are $U(r)$ and $\rho$. In principle, one could directly use the tabulated pair potentials $U(r)$ as an input of the NNs. In practice, however, this might cause numerical issues for strongly repulsive potentials, since a further increase of the interaction strength beyond a certain threshold has little or no effect. Therefore, it is sensible to reduce the search space and minimize the risk of overfitting by limiting the maximum repulsion of $U(r)$ during the optimization of the NNs. The range of the pair potential $U(r)$ was optionally constrained either by representing it in exponential space, i.e., $\exp\left[-\beta U(r)\right]$, or by using a clipped potential, $\bar{U}(r)=\max\left[U(r),U_{\rm cut}\right]\forall r$ with $U_{\rm cut}=20\,k_{\rm B}T$. The choice of $U_{\rm cut}$ is somewhat arbitrary but our tests indicated that the specific value does not significantly alter the outcome as long as $U_{\rm cut}$ separates soft and hard potentials.

Furthermore, we tested how additional information can improve the prediction accuracy. The force $F(r)$ was optionally included in the NNs, by applying the same operations as for $U(r)$ and then concatenating the outputs for further processing in fully dense layers (see Fig. 2). If $\bar{U}(r)$ was used as the input, then the force was calculated before clipping to avoid unphysical discontinuities. If the input was $\exp\left[-\beta U(r)\right]$, then $-\frac{\partial}{\partial r}\exp\left[-\beta U(r)\right]$ was used as the force input. Note that including the derivative almost doubles the number of nodes in the network. For the NNs learning $B_{2}^{*}$ or $P$ directly, we optionally provided also a pressure estimate $P_{0}$, calculated using the analytic second virial coefficient $B_{2}$. In a homogeneous system with isotropic interactions, $B_{2}$ is readily available from $U(r)$52

[TABLE]

with Mayer $f$-function

[TABLE]

Thus for $\rho\to 0$, the pressure $P_{0}$ can be directly computed from the pair potential $U(r)$

[TABLE]

In what follows, we will use the following naming convention for the NNs: For a given architecture X and given options y, we name the network X/y, e.g., for a UN with force information we use UN/f. The abbreviations for the different combinations are summarized in Table 1. Hence, we investigated 90 different combinations of network architectures and optional parameters. It is clear that optimizing the hyperparameters (see Sec. 2.1) for each individual case is computationally infeasible, and therefore we optimized only the DN/x, CN/x, and UN/x networks. Then, those optimized hyperparameters are used for all other networks of the same architecture to systematically investigate the effect of the additional information provided to the NN as well as the representation of the input and output data. In most cases, the variation of the hyperparameters resulted in changes of the prediction accuracy on the order of the variation between the folds, thus typically less than the effect of the parameters listed in Table 1. Only in selected cases, we found that the NNs were too small to capture the relevant details or too large to cause overfitting.

Figure (3) shows the mean relative absolute error between the predicted and target pressures, $\left\langle L_{\rm MRAE}\right\rangle=\left\langle\left|P-\hat{P}\right|/P\right\rangle$, for all investigated NNs. The errorbars in Fig. 3 were calculated as the standard error of the mean over all four folds. For comparison, the mean relative absolute error between the target pressure $P$ and the pressure $P_{0}$ estimated via $B_{2}$ was $\left\langle L_{\rm MRAE}\right\rangle\approx 1.4$. In most cases, the NNs provide significantly more accurate predictions compared to $P_{0}$, which is reasonable given that Eq. (11) is strictly valid only in the limit $\rho\to 0$. Further, the NNs that directly map $U(r)\mapsto\beta P$ perform worse than the NNs which first predict $B_{2}^{*}$ or $\mathcal{B}_{n}$ and then compute $P$.

Using either the clipped potential $\bar{U}(r)$ or $\exp[-\beta U(r)]$ as the input instead of $U(r)$ typically leads to better results and reduces overfitting in most of the cases. In contrast, including the forces did not improve the performance notably, but was even detrimental in some cases, as can be seen, for example, by comparing the performance of CN/fi and CN/i predicting $B_{2}^{*}$ [see Fig. 3(e)]. We surmise that this behavior occurs because including the force does not provide any crucial additional information compared to the pair potential, while increasing the complexity of the network architecture and thus increasing the risk of overfitting. Providing $P_{0}$ sometimes improved performance, but the effect was not as pronounced compared to the other parameters. For the UNs predicting $B_{2}^{*}$ [Fig. 3(h)], including $P_{0}$ resulted in a consistently lower accuracy, which was likely due to overfitting as indicated by the small loss on the training set compared to the loss in the validation set throughout.

Based on the comparison shown in Fig. 3, we chose UN/c predicting $\mathcal{B}_{n}$ as the final network because of the following reasons: The pressures predicted by this NN were among the most accurate with small mean relative absolute errors $\left\langle L_{\rm MRAE}\right\rangle\approx 0.06\pm 0.01$ and $\approx 0.11\pm 0.01$ for the training and test set, respectively. Further, the mapping $U(r)\mapsto\mathcal{B}_{n}$ should also generalize better to densities which were not part of the training set, as the predicted $\mathcal{B}_{n}$ are explicitly density independent, with the only density dependence of $\hat{P}$ coming from the virial expansion. Further, we used clipped potentials $\bar{U}(r)$, which guaranteed that the values provided to the NN are bounded, leading to better generalizability and less overfitting. Forces were not included in the final network, because they did not always improve the prediction accuracy, but increased the model complexity and the risk of overfitting.

For the final benchmark, all data (except those in the test set) were used for training the NN. Figure 4 shows the pressure $\hat{P}$ predicted by the UN/c network vs. the target pressure $P$. Here, we have also included the low-density estimate $P_{0}$ [see Eq. (11)]. The predicted pressures $\hat{P}$ are very close to the target values $P$, with coefficients of determination $R^{2}\approx 1.00$ for both the training and test set. In contrast, $P_{0}$ deviated strongly from $P$, especially for the systems at high density and/or with $P<0$. This rather poor agreement with the target pressure is reflected in the coefficients of determination, which was $R^{2}\approx 0.172$.

3.2 Potential prediction from pair distribution function

In this section we will consider NNs (see Fig. 5) for predicting (effective) pair potentials $U(r)$ from radial pair distribution functions $g(r)$ at a given, known particle number density $\rho$. We either used $g(r)$ as input, or computed the (effective) pair potential in the low density limit

[TABLE]

and used $U_{0}(r)$ as the input. The reasoning behind this approach is to provide the NNs with a physically informed estimate which is valid for $\rho\to 0$, so that the NNs only need to learn perturbations to this solution at higher densities. Networks using $U_{0}(r)$ as input are indicated with a “b” for Boltzmann inversion. Optionally, the NNs could optimize and output $\exp[-\beta U(r)]$ instead of $U(r)$ (indicated by “e”), and/or include the force (indicated by “f”) in the loss function through Eq. (7).

The loss was always computed between the predicted, $\hat{U}(r)$, and target pair potential, $U(r)$. To reduce the search space and avoid overfitting, we disregarded parts of $U(r)$ that exceeded $U_{\rm cut}=20k_{\rm B}T$, effectively yielding a sliced loss function. Hence, the predicted $U(r)$ were not optimized in regions where $U(r)>U_{\rm cut}$, allowing the NNs to set arbitrary values for the corresponding ranges.

Methods such as RMC or IBI instead optimize $\hat{g}(r)$ against the reference $g(r)$ by iteratively performing MD simulations with the predicted potential $\hat{U}(r)$ (see ESI for more details). Such an approach is, however, prohibitively time consuming as it would require additional MD simulations at each step of the optimization procedure. Further, minimizing the loss between $\hat{g}(r)$ and $g(r)$ is infeasible from a technical point of view, because the training of the NNs is based on gradient descent, which requires the (unknown) mapping $U(r)\mapsto g(r)$. Thus, the difference between $g(r)$ and $\hat{g}(r)$ is not a viable input for the loss function for training the NNs, but it can be used as a final performance benchmark in selected instances (see Fig. 10 below).

Again, we tested the DN, CN and UN architectures, and systematically analyzed the influence of providing additional information and different representations of the input and output data. Figure 6 shows $\left\langle L_{\rm MSE}\right\rangle$ and $\left\langle L_{\rm MAE}\right\rangle$ between the predicted and the target potentials for all investigated NNs and parameter combinations listed in Table 2. The errorbars were obtained through a 4-fold cross validation. The dotted lines correspond to the MSE (blue) and MAE (green) between $U(r)$ and $U_{0}(r)$ [note that we employed here the same clipping as for $\hat{U}(r)$].

As can be seen from Fig. 6, almost all investigated NNs provide more accurate predictions compared to $U_{0}(r)$. The predicted potentials $\hat{U}(r)$ from the NNs optimizing $\exp[-\beta\hat{U}(r)]$ have consistently larger $\left\langle L_{\rm MSE}\right\rangle$ and $\left\langle L_{\rm MAE}\right\rangle$ compared to the NNs which directly output $\hat{U}(r)$. This discrepancy likely originates from the non-linearity of the transformation, which effectively modifies the importance of different regions of the potential. Figure 6 also shows that including the force in the loss function during the training stage improves the accuracy of the NNs (see Sec. 2.1). Further, the predicted potentials were significantly smoother when the force was included, as shown for one selected example in Fig. 7. Finally, we find that inputting $U_{0}(r)$ instead of $g(r)$ did not improve the prediction significantly. Comparing the average losses $\left\langle L_{\rm MSE}\right\rangle$ and $\left\langle L_{\rm MAE}\right\rangle$ of the different network architectures for a given parameter set, we see that the CNs were most accurate in most cases for both the training and validation sets. In particular, the CN/f network had the smallest average losses and also small variations between the different folds. For final benchmarking, we therefore focuses on the CN/f network, which was trained again with all data except those in the test set.

To evaluate the accuracy of this NN in more detail, we analyzed the probability density function of $L_{\rm MAE}$ for the training and test set, $p(L_{\rm MAE})$, and compared it with $p(L_{\rm MAE})$ obtained from using $U_{0}(r)$ (see Fig. 8). In both cases, $p(L_{\rm MAE})$ can be fitted by an exponential decay $\propto\exp-(L_{\rm MAE}/\delta)$, with $\delta$ quantifying the width of the distribution. For the CN/f network, $p(L_{\rm MAE})$ had a rather narrow distribution with $\delta=0.02\,k_{\rm B}T$ and $\delta=0.03\,k_{\rm B}T$ for the training and test set, respectively, whereas $p(L_{\rm MAE})$ was considerably wider for $U_{0}(r)$ with $\delta=0.5\,k_{\rm B}T$ (the asymptotic standard error of the fits was about $10\,\%$ in all cases). It should also be noted that the estimated $U_{0}(r)$ had several outliers with very large deviations up to $L_{\rm MAE}\approx 8\,k_{\rm B}T$.

At this point, we also want to discuss the transferability of the predicted potentials $\hat{U}(r)$. Given that $g(r)$ can vary significantly for the same $U(r)$ at different particle number densities $\rho$, it is important to check whether $\hat{U}(r)$ remains independent of $\rho$. Figure 9 shows a selected target potential $U(r)$ from the test set and the resulting $g(r)$ at four different $\rho$, highlighting the increasing deviations from ideal-like behavior [$g(r)=1$] with increasing density. We have also included in Fig. 9(a) the predictions of $\hat{U}(r)$ from the CN/f network as well as $U_{0}(r)$ at different $\rho$. Indeed, the predictions $\hat{U}(r)$ are rather close to $U(r)$, and show only minor variations with $\rho$ (as should be the case). In contrast, the estimate $U_{0}(r)$ provides a passable approximation of $\hat{U}(r)$ only at the lowest density $\rho=0.125\,\sigma^{-3}$, and the agreement with $U(r)$ gets significantly worse with increasing $\rho$.

To check whether the predicted $\hat{U}(r)$ lead to similar pair distribution functions $\hat{g}(r)$ as the original $g(r)$, we ran MD simulations with the predicted potential $\hat{U}(r)$ for selected cases from the training and test set. To test the performance of simple Boltzmann inversion, we also ran MD simulations with $U_{0}(r)$ and calculated the corresponding $g_{0}(r)$. Figure 10 shows the comparison between $g(r)$, $\hat{g}(r)$ and $g_{0}(r)$, where we have also included the employed pair potentials $U(r)$, $\hat{U}(r)$, and $U_{0}(r)$ as insets of the corresponding plots. In all considered cases, $\hat{g}(r)$ is much closer to $g(r)$ than $g_{0}(r)$, especially for distances $r\gtrsim\sigma$. The systems shown in Fig. (10)(a,b) have been run with pair potentials from the training set at densities $\rho=0.244\,\sigma^{-3}$ and $0.512\,\sigma^{-3}$, respectively. These cases are rather interesting, because $U(r)$ and $\hat{U}(r)$ look almost identical, but the resulting $g(r)$ and $\hat{g}(r)$ show some deviations for $r\lesssim\sigma$. The estimated potentials $U_{0}(r)$ differ drastically from $U(r)$, but the resulting $g_{0}(r)$ for $\rho=0.244\,\sigma^{-3}$ is relatively close to the target one. At higher density $\rho=0.512\,\sigma^{-3}$, however, the simulations performed with $U_{0}(r)$ lead to a distinct clustering of particles, while the simulations with $U(r)$ and $\hat{U}(r)$ lead to an almost flat, gas-like $g(r)$ with weak oscillations. Thus, in this case, $U_{0}(r)$ would not even be a good starting point for iterative optimization methods like IBI. Figure (10)(c,d) shows results for cases from the test set (same systems as in Fig. 7) at two different densities. At low density $\rho=0.244\,\sigma^{-3}$, the agreement between $g(r)$ and $\hat{g}(r)$ is almost perfect, while $g_{0}(r)$ replicates the qualitative correct trends. For the denser systems $\rho=0.512\,\sigma^{-3}$, both $\hat{g}(r)$ and $g_{0}(r)$ exhibit a qualitatively similar deviation from the target radial distribution function, but the predicted $\hat{U}(r)$ is much closer to the target pair potential $U(r)$ compared to $U_{0}(r)$.

As a final test of our methodology, we used the predicted pair potentials $\hat{U}(r)$ from the CN/f networks as an input of the UN/c networks developed in Sec. 3.1 to calculate the pressure $\hat{P}$ from the radial pair distribution function $g(r)$. Figure 11 shows $\hat{P}$ vs. $P$, demonstrating rather good agreement for the training set ($R^{2}\approx 0.971$). The predictions for the test set have a significantly lower accuracy ($R^{2}\approx 0.461$), but the overall trends are captured in most cases (excluding the outlier at $P\approx 80\,k_{\rm B}T/\sigma^{3}$ from the test set leads to $R^{2}\approx 0.802$). For comparison, we have also included the pressure $P_{0}$, which we estimated by computing $U_{0}(r)$ from $g(r)$. These pressures are almost always much smaller than the target pressure, leading to $R^{2}\approx 0.0192$.

4 Conclusions

Artificial neural networks (NNs) were developed for predicting the equation of state for a given pair potential $U(r)$ and density $\rho$, and for predicting (effective) pair potentials from structural information based on the radial distribution function, $g(r)$. We investigated how the representation of the input and output data as well as additional information (e.g., the forces) influence the prediction accuracy. For both tasks, a key preprocessing step was to limit the input and/or output range by capping the pair potential $U(r)$ so that it could not reach arbitrarily large values. For predicting the pressure $P$, NNs were tested that directly map the input pair potential $U(r)$ to $P$, and NNs that first predict (effective) virial coefficients from $U(r)$ which are then used to compute $P$. The latter strategies resulted in rather accurate predictions, which were much closer to the target pressure ($R^{2}\approx 1.00$) compared to the virial expansion derived in the low-density limit. For predicting (effective) pair potentials, the accuracy of the predictions improved significantly when also the derivatives of the potentials were provided in the loss, since more correlations were taken into account. As a result, the predicted potentials became smoother and were significantly closer to the target potentials compared to results obtained through simple Boltzmann inversion.

There is, however, still room for future improvements: the potentials calculated from our NNs reproduce structural data with lower accuracy compared to potentials obtained from iterative methods, such as Iterative Boltzmann Inversion, which are constructed to reproduce $g(r)$. Further, the systems used in this work only covers a subset of all conceivable cases. Nevertheless, our NN approach is still a viable option for obtaining reasonably accurate initial estimates, which can then be optimized further using other methods until the desired accuracy is achieved. Such a combined approach could drastically cut down the computational cost and development time for coarse-graining applications, and it could also be a useful tool for inverse problems in materials discovery in soft matter.

Supporting Information

Plots of all potentials and radial distribution functions used for training and testing; Additional information on theoretical background; Schematics of final network architectures

Acknowledgements

We thank M. R. Khadilkar for fruitful discussions. This work was funded by the German Research Foundation (DFG) through project number 233630050 - TRR 146. AN further acknowledges financial support provided by the DFG through project number NI 1487/2-1 and NI 1487/2-2. Computing time was granted on the supercomputer Mogon at Johannes Gutenberg University Mainz (www.hpc.uni-mainz.de).