Importance of feature construction in machine learning for phase transitions

Machine learning (ML) has become an important tool in the study of phase transitions in molecular simulations.¹ Supervised and un-supervised methods have been successfully applied to a variety of systems.^2–12 While the majority of studies have focused on lattice models such as the 2D Ising model,^8,13–18 there have also been supervised machine learning studies of complex fluids, such as polymer blends¹⁹ and polymers in ionic liquids.⁷ In this work, we use un-supervised machine learning methods to study the phase behavior of two continuous-space binary mixtures.

Machine learning methods can be broadly classified into two categories: supervised and un-supervised. In supervised methods, such as neural networks and support vector machines, the algorithm is trained on known data and used to predict the behavior of an unknown evaluation set. In un-supervised methods, a method is used to decrease the dimensionality of the input data, and correlations (or clustering) in the reduced dimensions are used to provide insight into the behavior. The objective is to discover patterns in the data without any prior training.^2,13 The advantage of un-supervised methods for phase behavior is that prior knowledge of the existence of a phase transition is not required. Furthermore, they can be useful for complex fluids where molecular simulation of the phase behavior using conventional methods (required for a training set in supervised methods) is computationally intensive. In this work, we use un-supervised machine learning methods with two dimensionality reduction techniques, namely, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE), to study the phase behavior of off-lattice models.

There have been some studies of continuous-space methods using un-supervised ML methods. Jadrich et al. used PCA to study the phase behavior of 2D model (continuous-space) particle systems, such as the RandOrg model, hard ellipses, and non-additive, hard sphere Widom–Rowlinson (WR) mixture, as well as 3D hard sphere and patchy colloid mixtures.^3,4 They emphasized the importance of “feature” construction in the application of PCA to particle-based systems. The feature vector is the data that is input to the ML method. In spin models, such as the 2D Ising model, the raw configuration, i.e., the values of all the spins, is usually an excellent choice as the input feature vector. For particle-based systems, the positions of the particles prove to be uninformative. Inter-particle distances are more useful, since these are invariant to rotation of the simulation cell.^3,4,20

In this work we study the phase behavior of two 3D off-lattice systems showing liquid–liquid phase behavior, a binary Lennard-Jones (LJ) mixture and the Widom–Rowlinson (WR) mixture, using PCA and t-SNE techniques. We find that the construction of the feature vector is crucial to the ability of the ML method to predict the phase behavior. A distance-based feature provides good estimates for the critical point of the LJ mixture, but not for the WR mixture. We propose an affinity-based feature, where the elements are +1 if a particle has predominantly neighbors of the same species and −1 otherwise. We find that this feature provides an excellent description of the phase behavior for both models.

The two models we investigate are an LJ mixture and the WR mixture.²¹ In the symmetric LJ mixture, the particles interact via a truncated and force shifted 12-6 LJ potential, i.e., the potential of interaction between a particle of species A and a particle of species B, Φ_AB, is given by

where

σ_AA = σ_BB = σ_AB = σ, r_c = 2.5σ, and ϵ_AA = ϵ_BB = 2ϵ_AB = ϵ. Since the attraction between like species (A–A and B–B) is stronger than that of unlike (A–B) species, this system phase separates upon cooling. We define a reduced temperature by T* = k_BT/ϵ, where k_B is Boltzmann’s constant and T is the temperature, and a reduced density ρ = Nσ³/V, where N is the total number of particles, N_A + N_B = N, where N_A and N_B denote the number of A type and B type particles, respectively. V is the volume of the simulation box.

In the WR mixture, there is no interaction between particles of the same species, and particles of unlike species (A-B) interact via a hard sphere interaction, with diameter σ_HS. This system phase separates when the density is increased. We define a reduced density $ρ=NσHS3/V$⁠. All simulations are performed for a mole fraction x = 0.5 which, by symmetry, is the critical composition of both mixtures.

Configurations of the WR mixture are obtained from Monte-Carlo simulations in the NVT (number of particles, volume, and temperature fixed) ensemble. Initial configurations are obtained by placing N particles in a cubic cell with linear dimension L. The system is evolved via single particle moves, where a randomly chosen particle is moved at a random angle with displacements uniform on (0, σ_HS). The system is equilibrated for 10⁶ steps per particle. The production run consists of 10⁷ attempted moves per particle from which 1000 configurations are extracted. A total of 51 number densities are explored, ranging from $ρσHS3$ = 0.5–1. We study six system sizes with N = 1024, 2048, 3072, 4096, 8192, and 16 384.

Configurations of the binary LJ mixture are obtained via molecular dynamics (MD) simulations in the canonical (NVT) ensemble using the LAMMPS package.²² Initial configurations are obtained by placing N particles in a cubic cell with linear dimension L = V^1/3, such that ρ = 1. The system is propagated using a velocity verlet integrator, with a time step of Δτ* = 0.004, where τ* is the reduced time (⁠$τ*=tσ2ϵ/m$, where m is the mass of a particle). A Nose–Hoover thermostat is used to maintain the desired simulation temperature.^23,24 The initial configuration is equilibrated for 2 × 10⁶ time steps (8000 τ*) and properties are averaged over 10⁷ time steps (40 000 τ*) from which 1000 configurations are extracted and used. Simulations are performed for temperatures ranging from 1.0 ≤ T* ≤ 1.8 for distance feature matrix calculations and from 0.9 ≤ T* ≤ 2.375 for affinity matrix calculations. For affinity matrix calculations, we study varying system sizes with L/σ = 10, 12, 14, 16, 18, 20, 24, and 30.

The feature vector is a one-dimensional vector that comprises the data from each configuration that is input into the ML algorithm. In spin systems, this is usually a list of the values of all the spins in a configuration. In continuous-space systems, the choice is not obvious. We test two different choices for the feature vector in this work. The first is based on the work of Jadrich et al.,⁴ which we refer to as the distance-based feature. The second is a choice that we propose, inspired by spin systems, and we refer to this as the affinity-based feature.

The distance based feature vector is constructed as follows [see Fig. 1(a)]: For each configuration, we choose a probe particle of one of the species and then calculate the distances from this particles to the N_c nearest neighbors of the same species. We select a particle of species A as a probe and use the distances from N_c nearest neighbors for analysis as described below. If v_i is the distance from the probe to the ith nearest neighbor, the average distance is

We define a new distance $xi=vi−v̄$ and manually re-order the elements such that $x1≤x2≤…≤xNc$⁠. A single component of the feature vector is given by

For the WR mixtures, the particles’ co-ordinates are first normalized by the box length.

The feature matrix X is constructed by stacking all the feature vectors

where N_T is the number of configurations. Since there are N_c neighbors, X is a matrix of dimension N_T × N_c. We average the results over ten probes in the case of the LJ mixtures and N/2 probes in the case of the WR mixture. We perform PCA for each feature matrix constructed from each probe and then derive the mean and standard deviation profile. Then, we average them to get the final results.

The affinity-based feature vector is constructed as follows [see Fig. 1(b)]: For each particle, we calculate the number of particles of each species within a cut-off distance r_cut, which we set equal to half the box length. If the number of particles of the same species around particle i is greater than the number of particles of the other species, then the quantity g_i = 1, otherwise g_i = −1.

For each configuration the feature vector is given by

where N is the total number of particles. The feature vector matrix is an N_T × N matrix,

We optimize the value of r_cut in the affinity-based feature vector by selecting a value that can best distinguish between different structures. Figure 2 depicts the probability distribution function, $Pḡ$⁠, of the average feature vector, $ḡ=1N∑i=1Ngi$⁠, for two values of r_cut. There is a greater distinction between the different values of $ḡ$ for r_cut = 0.5L. We choose a value of r_cut where the variance $σḡ2$ is the largest. For both the WR and LJ mixtures, this occurs for r_cut = 0.5L.

We use Principal Component Analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE)^25,26 to obtain a low-dimensional representation of our data. PCA uses an orthogonal linear mapping of the original high dimensional data into a lower dimensional space. The method we use is standard, but we describe it briefly here for completeness.

The first step is standardization of the feature matrix to give a matrix X_S, where the mean of each column is zero. The covariance matrix

is a positive, semi-definite, symmetric matrix (of dimension N_c × N_c or N × N) with non-negative eigenvalues and orthogonal eigenvectors. The eigenvalues b_n and eigenvectors w_n are obtained by solving

The transformation matrix has the eigenvectors as columns, i.e., $W=w1T,w2T,…,wNT$⁠, where the eigenvectors are ordered so that they correspond to a decreasing order of eigenvalues, b₁ > b₂ > b₃, etc. The principal component matrix is given by

The principal component (PC) axes are the columns of P. Since the eigenvalues are ordered in decreasing magnitude, the first PC will have the largest variance. A dimensionality reduction is achieved by choosing the first few principal components. If b_n is the variance of PC n, called the explained variance, the relative explained variance λ_n is defined as

where M is the dimension of W.

Figure 3 depicts the relative explained variances plotted against the PC index for the LJ mixture and WR models. In both models, λ_n decreases slowly with n with the distance-based feature, but the first index is dominant with the affinity-based feature. The affinity-based feature, therefore, provides a better reduction in dimensionality.

We quantify the phase transition using PCA-derived order parameters (OPs). These OPs are calculated by using only the first two PC scores: in other words, only the first two columns of the final matrix P. The first two columns contain the first two principal components P₁ and P₂ for each state point (temperature or density). One OP is the ensemble average of P₁, denoted $P1$⁠. We normalize $P1$ so that the range is between 0 and 1. The second OP is related to the variance of the principal components. If $σ12$ and $σ22$ are the ensemble averages of the variances of P₁ and P₂, respectively, the OP is defined as $σ12≡σ12+σ22/2$⁠.

The theoretical framework of t-SNE is based on effective minimization of the divergence between the two conditional probabilities (P and Q). Here, for the distance feature matrix, we have a higher dimensional input dataset X_S = $[x1,x2,…xNT]$⁠. For the input, we use the PCA transformed data. The first step is measuring the conditional probability of pairwise distances between the input data points based on a Gaussian distribution,

where

and

Here, σ_i is the variance of the Gaussian centered on each high dimensional data point x_i. This parameter is automatically determined by the algorithm for a given value of user-provided perplexity, which is basically the effective number of neighbors around x_i, mathematically given as

where H(P_i) is Shannon entropy. In order to keep the symmetry between P_i∣j and P_j∣i, t-SNE uses a probability function P_ij, which is the mean of two probabilities

where N is the size of the input dataset.

The second conditional probability (Q) measures the similarity between the high-dimensional data points in lower dimensional space. The low dimensional mapping is represented as Y = $[y1,y2,…yNT]$ and the embedding similarity Q_ij between two low dimensional data points y_i and y_j is calculated using a Student-t kernel given as

where

and

The t-SNE approach minimizes the Kullback–Leibler divergence that is defined using a cost function (C), given as

In order to minimize the difference between two probability distributions P_ij and Q_ij and updating Q_ij iteratively, the cost function is minimized by descending along the gradient

Solving this gradient iteratively allows t-SNE to find a low dimension projection y_i of the corresponding input x_i. Since, here we are projecting our data to a two-dimensional space, the final low-dimensional matrix is of size N_T times two, and the columns give us the first and second components of the projections. These are labeled S₁ and S₂ in the scatter plots. In the case of t-SNE, the order parameters are the same as described in connection with the PCA, although they are denoted as $S1$ and σ₁₂.

All the PCA and tSNE calculations were performed by using the libraries available in scikit-learn.²⁷ 

A scatter plot of the PCA and t-SNE as a function of state point shows that the affinity based feature vector more clearly distinguishes between the phase separated and mixed states, when compared to the distance-based feature, for both the LJ mixture and WR model. Interestingly, there is no significant difference in the ability of the PCA or t-SNE to distinguish the phases. These conclusions are evident in Figs. 4 and 5, which depict the scatter plots for the two models with the distance based (Fig. 4) and affinity-based (Fig. 5) feature vectors, respectively. In the figures, the colors represent the state point going from mixed (red, high temperature in the LJ mixture and low density in the WR model) to phase separated (blue, low temperature in the LJ mixture and high density in the WR model).

A comparison of the PCA with the t-SNE (compare left to right panels in all figures) shows that there is no significant distinction between the ability of the PCA and t-SNE to distinguish mixed and phase separated states, i.e., the degree of clustering of data points is similar. With the distance-based feature, both methods show some clustering in the case of the LJ mixture, but are unable to display a meaningful separation in the case of the WR model. With the affinity-based feature, both methods show a clear clustering, with the mixed phase (red) distinct from the phase separated phase (blue). This suggests that the affinity-based feature vector is a more promising route to extracting the phase behavior of these models.

Often higher order components in the PCA are required to distinguish different structures. This is particularly relevant here because, with the distance-based feature, the explained variances show more than one dominant component (see Fig. 3). The higher components do not provide significant additional information in this case, however. Figure 6, which is a three-dimensional plot of the first three components of the PCA (P₁, P₂ and P₃), does not show a clear distinction between the mixed and separated states.

The order parameters obtained with the affinity-based feature vectors provide an accurate estimate of the critical point of both models. The order parameters with the PCA are depicted in Fig. 7 and those with the t-SNE are in the supplementary material (Fig. S1). The order parameter $P1$ is zero in the mixed phase and increases rapidly in the two-phase region, and σ₁₂ shows a peak in the neighborhood of what could be the critical point. Interestingly, the behavior of σ₁₂ is quite similar to the constant volume specific heat (Fig. S2).

For the smaller systems, the PCA is not able to distinguish between the mixed and phase separated states, for both models, but a clear transition is seen for larger systems. The order parameter $P1$ often shows a gradual increase rather than a rapid increase to non-zero values, and we, therefore, identify the peak in σ₁₂ with the phase transition. For the systems shown in Figs. 7(b) and 7(d), which are the largest system sizes for each case, we estimate a critical temperature of 1.425 for the LJ mixture, which compares well with the literature value 1.423 ± 0.0005.²⁸ Similarly, the transition density for the WR model is 0.75, which compares well with the literature value of 0.762.²⁹ Results for other system sizes are shown in Figs. S3–S8 of the supplementary material.

We estimate the infinite system transition point using finite size scaling. We find the transition temperatures and densities from the position where σ₁₂ is the highest. Uncertainties in the temperature or density are measurement uncertainties, i.e., the simulation intervals of temperature (LJ) or density (WR). The machine learning transition temperatures do not show the usual scaling expected from finite size simulations in the Ising universality class.³⁰ Usually, T_c(L) ∼ T_c,∞ + AL^−1/ν, where T_c,∞ is the critical temperature of the infinite system, ν = 0.629 is the Ising finite size scaling exponent, and A is a positive constant. This implies that T_c(L) > T_c,∞ for finite values of L. This is because, in a finite system, the critical point is reached when the length-scale of fluctuations exceeds the box size, and this happens at a higher temperature for smaller systems.

We find that for the WR model, T_c is essentially independent of L (within uncertainties), and for the LJ mixture, T_c increases with increasing L. The variation of T_c with L^−1/ν for the LJ mixture and WR model are shown in Figs. 8(a) and 8(b), respectively. From a fit to T_c(L) = T_c,∞ + AL^−1/ν, we extract T_c,∞ = 1.43 ± 0.01 for the LJ mixture, and ρ_c,∞ = 0.76 ± 0.01 for the WR mixture. These estimates are in good agreement with the critical point obtained using conventional simulations.

We use unsupervised machine learning methods to obtain the critical point of three dimensional off-lattice systems. The most important result is that the construction of the input feature vector plays a crucial role in the ability of the ML method to predict the critical point. We use two feature vectors, one based on the distances between particles (“distance-based”) and one based on the local concentration (“affinity-based”), and study two systems, a symmetric Lennard–Jones binary (LJ) mixture and a non-additive, hard sphere Widom–Rowlinson (WR) mixture. We find that the distance-based feature vector is not able to distinguish clearly between mixed and two-phase states for the WR mixture, although it does show a distinction for the LJ mixture. The affinity-based feature vector, however, shows a distinct separation of one-phase and two-phase configurations for both models. We use two ML methods, PCA and t-SNE, and find that the choice of the ML algorithm is not important—the key choice is the construction of the feature vector.

We estimate the critical point from the variation of the peak in variance of the first two principal components with temperature and density, and extrapolate to an infinite system using the critical scaling of the Ising universality class. The critical point obtained in this fashion is in good agreement with previous estimates from conventional simulations, and requires standard simulations with no particle insertions or deletions.

Supervised and un-supervised ML methods have been widely used for simple lattice systems and some off-lattice systems. In lattice systems, the choice of the feature vector is clear, namely, the values of the spins. We find that in off-lattice systems, a judicious choice of the input features can play a crucial role in the predictions, and it remains to be seen if the choice of feature is system dependent. We hope this paves the way to study other complex fluids using un-supervised methods.

An interesting direction is to connect the machine learning methods to statistical mechanics and liquid state theory. For the LJ mixture, we find empirically that the order parameter σ₁₂ shows behavior quite similar to the constant volume heat capacity. Although this suggests that there might be a link, the former relies only on the locally averaged concentrations, whereas the latter is a measure of energy fluctuations. Note that ML relies on how configurations at a given state are different from those at other states, and it is, therefore, important to have a range of temperatures (in the LJ case) that span the critical point. In statistical mechanics, on the other hand, one only requires a single state to determine all the properties. In this regard, it is important to note that identifying the peak in σ₁₂ with the critical point is merely an ansatz, and one that must eventually be verified to represent the critical point from a more fundamental standpoint.

See the supplementary material for analysis of the phase behavior using the distance-based feature vector and order parameters with the t-SNE method.

This work was supported by the National Science Foundation through Grant No. CHE-1856595. All simulations presented here were performed using computational resources provided by UW-Madison Department of Chemistry HPC Cluster under NSF Grant No. CHE-0840494.

The authors have no conflicts to disclose.

Inhyuk Jang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Supreet Kaur: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Arun Yethiraj: Conceptualization (equal); Funding acquisition (lead); Investigation (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.