Machine Learned Phase Transitions in a System of Anisotropic Particles on a Square Lattice

The area of Machine learning has seen a tremendous growth in the last two decades. Researchers have made great progress in areas like computer vision Sultana et al. (2018), natural language processing Young et al. (2018), medical diagnostics de Bruijne (2016); Ahmed et al. (2019) using deep learning methods LeCun et al. (2015); Schmidhuber (2015). The availability of a large amount of data and higher computing power due to the advancements in hardware have made the growth possible. Machine learning techniques usually perform better with a large amount of data. Data-rich branches of physics like High energy physics and Astronomy thus have employed Machine learning techniques successfully in extracting physical insights from the data BALL and BRUNNER (2010); Baldi et al. (2014); Albertsson et al. (2018); Das Sarma et al. (2019). In recent years Machine learning has made its way into other branches like Condensed matter physics and Statistical physics. Both supervised and unsupervised learning methods have been successfully employed to detect phase transitions and classify different phases Arai et al. (2018); Broecker et al. (2017); Iakovlev et al. (2018); Carrasquilla and Melko (2017); Wetzel (2017); Wang (2016a); Bedolla-Montiel et al. (2020); Liu and van Nieuwenburg (2018); Ponte and Melko (2017); Ch’ng et al. (2017); Boattini et al. (2019); Wang (2016b). Neural networks have shown great potential in identifying new exotic phases and phase transitions even in systems where the order parameter cannot be defined explicitly Li et al. (2018); Lian et al. (2019). Generative neural networks like restricted Boltzmann machines and variational autoencoders have been employed to model physical probability distributions and extract features in spin models D’Angelo and Böttcher (2020); Morningstar and Melko (2018).

Here, we focus on the specific application of machine learning in classifying different phases separated by a phase transition. Usually, in case of a structural phase transition one can construct a suitable order parameter to distinguish between the different phases. Away from the transition point, one may even be able to identify the phases just by looking at the snapshots. However, this distinction gets blurry as one approaches the transition more and more closely (at least for a continuous or a weakly first-order transition). Especially, in the vicinity of a critical point, large fluctuations are expected to interfere. Now the question is whether one can train a statistical classifier to distinguish between such phases just by using the configurations. The classification of the ordered and disordered phases, and the study of the critical properties of Ising model have been extensively investigated in recent times Morningstar and Melko (2018); D’Angelo and Böttcher (2020); Wang (2016b); Mehta et al. (2019); Tanaka and Tomiya (2017). After the success of Machine learning methods in characterizing simpler models like Ising model, they are now being explored in studying more complex systems, like liquid crystals. In experiments, liquid crystals are usually studied using optical imaging Zola et al. (2013). Recently Machine learning techniques are applied to classify phases and predict the physical properties of liquid crystals from the experimental data Terao (2020); Sigaki et al. (2019); Minor et al. (2019). Convolutional neural networks (CNNs) have been shown to successfully classify nematic and isotropic phases in continuum with very high accuracy Sigaki et al. (2020). In this case the isotropic-nematic transition is first-order in nature where the order parameter changes abruptly in the vicinity of the transition point.

In this context, we ask how difficult is the classification task near a continuous transition when compared with the same around a first-order transition. We investigate how different machine learning models perform on this classification task and how the system size affects the performance. We further ask if physics-guided features help in learning the phases. In addition, we examine if the method of Learning by confusion can be extended to the system of rods to determine the critical (without any prior estimate) using the information about the ordered phase van Nieuwenburg et al. (2017); Tan and Jiang (2020). This technique is quite powerful as the usual supervised ML techniques rely on the prior knowledge of the critical point to be trained– this may not always be possible for all complex systems van Nieuwenburg et al. (2017). Understanding these issues can help us extending ML techniques to more complex systems and also tackling systems where the identification of structural order (if exists) is nontrivial, e.g., amorphous systems Bapst et al. (2020) .

To investigate the above questions, we study a model of long hard rods on a two-dimensional square lattice that undergoes isotropic-nematic-disordered phase transitions with increasing density. Note that, systems of hard rectangles with finite width exhibit more complex liquid crystalline phases Kundu and Rajesh (2014). Such hard core lattice gas models are, in general, relevant for understanding self-assembly of nanoparticles Rabani et al. (2003), glass transition Diaz-Sanchez et al. (2002), adsorption of gas molecules on metal substrates Taylor et al. (1985); Bak et al. (1985); Patrykiejew et al. (2000) and entropy-driven phase transitions (realized in liquid crystalline assemblies of various colloidal systems Kuijk et al. (2011, 2012); M. P. B. van Bruggen and Lekkerkerker (1996); Zhao et al. (2011); Barry and Dogic (2010)). In contrast to the continuum model, the high-density phase of the system of hard rods on a lattice is an orientationally disordered phase that remained inaccessible until recently due to large relaxation times Kundu et al. (2012, 2013); Kundu and Rajesh (2014). Using a novel Monte Carlo algorithm with nonlocal moves one can thermalize the system at very high densities to show that it undergoes two continuous transitions with increasing density: first from a low-density isotropic to an intermediate density nematic phase and second, from the nematic phase to a high-density disordered phase Kundu et al. (2013); Kundu and Rajesh (2014). We generate Monte Carlo sampled data for the hard rods system on a square lattice to classify the isotropic (I) and nematic (N) phases around the first I-N critical point using various Machine learning techniques. We show how the system size affects the performance of the models. Note that this classification task is trivial when the I and N phases are far from the critical point. However, closer to the critical point, these two phases are not visually distinguishable (see Fig. 1) as discussed above, making the task much harder. The same is true for the Ising critical point Mehta et al. (2019). Next, by breaking the symmetry between the two different orientations, we induce a first-order phase transition in the system to further show that the classification task gets easier around an abrupt transition. We then train logistic regression and random forest on physical features and compare the results with the models trained on lattice data. Further, we show that one can use simple features rather than the raw snapshots to get better accuracy with simpler ML models.

In the above-mentioned methods of classifying the phases, one must know the critical point in advance to label the data for training. Other methods such as clustering can be used which do not need prior knowledge of the critical point. It is shown that the neural networks can be trained to calculate the critical point with only theoretical ground state snapshots of ordered phase for models like ferromagnetic and anti-ferromagnetic Potts models Tan and Jiang (2020); Tan et al. (2020). Here we employ a similar strategy by training CNN on ordered phase ground state lattice to estimate the critical point.

The rest of the manuscript is organized as follows. We briefly describe the model, the algorithm to generate the equilibrium configurations, and the phenomenology in Sec. II. In Sec. III, we discuss the classification task around the I-N critical point using three machine learning techniques- logistic regression, deep neural networks, and convolutional neural networks that are trained on the I-N transition data. Next, we show the results including the performance of the classification tasks near the second-order and the first-order transitions in Sec. IV. We then show how physical features improve the performance of the models like logistic regression and random forests. In Sec. V we show how the nematic phase information can be used to estimate the critical point of the system of hard rods. Finally, we conclude in Sec. VI.

We consider a system of hard rods of length $k$ on a square lattice of size $L\times L$ with periodic boundary conditions. Each rod occupies $k$ consecutive lattice sites along one lattice direction and thus, can have two possible orientations: horizontal and vertical. The only constraint is that no two rods can overlap. All the configurations with no overlap are equally likely. Each rod is associated with an activity $z=e^{\mu}$ where $\mu$ is chemical potential. The system is treated using grand canonical ensemble where $z$ controls the density (fraction of sites occupied by the rods) of the system.

To simulate the system, we use the algorithm presented in this ref. Kundu and Rajesh (2014). The Monte Carlo algorithm is described below: at each step a row or a column is randomly selected. If a row(column) is selected all the horizontal(vertical) rods in that row are removed. This is the evaporation step. Now we end up with segments of empty sites separated by the forbidden sites due to the vertical rods passing through them. The next step is to reoccupy these empty segments with horizontal rods following correct statistical weights. The problem of occupying the empty row is reduced to filling the empty sites in one dimension. And the probabilities of the new configuration can be calculated exactly Kundu et al. (2013, 2012). This is a deposition step. A Monte Carlo step consists of $2L$ such evaporation-deposition steps. The equilibration is performed for $6\times 10^{5}$ MC steps. And then the snapshots of the system are taken at an interval of $5000$ MC steps to ensure they are uncorrelated.

The system undergoes two continuous phase transitions with increasing density when $k\geq 7$: first from a low-density isotropic phase to a nematic phase and second, from the nematic phase to a high-density disordered phase. We try to classify the isotropic and the nematic phases around the first transition. The critical chemical potential $\mu_{c}$ is determined from the probability distribution $P(Q)$ of the order parameter $Q$ defined as $Q=(n_{h}-n_{v})/(n_{h}+n_{v})$, where $n_{h}$ and $n_{v}$ are the number of horizontal and vertical rods correspondingly (note for the nematic phase $Q>0$ and for the isotropic phase $Q\approx 0$). Below $\mu_{c}$, $P(Q)$ has only a single peak around $Q=0$ corresponding to the isotropic phase and above $\mu_{c}$, $P(Q)$ develops two symmetric peaks corresponding to the nematic phase (see Fig. 2). Further, we study a first-order transition in the system by introducing a variable $\Delta$ (equivalent to the external magnetic field in Ising model) that breaks the symmetry between horizontal and vertical rods. The corresponding chemical potentials for the two types of rods are $\mu_{\rm h}=\mu+\Delta$ and $\mu_{\rm v}=\mu-\Delta$. As $\Delta$ is varied from a negative to a positive value at a given $\mu$, the system undergoes a first-order transition from a horizontal rod rich phase to a vertical rod rich phase (both are nematic). The order parameter $Q$ changes abruptly as shown in Fig. 3. In the following sections, we discuss the classification problem around the I-N criticality and the first-order transition point as mentioned above.

Typical snapshots of the system in the disordered and the nematic phases close to the I-N critical point are shown in Fig. 1. It is evident from the snapshots that they are not visually distinguishable. Hence, we attempt to classify them using Machine learning. The data around the I-N critical point is trained on logistic regression, deep neural network (DNN) and convolutional neural network (CNN). The data are for lattice size $L=98$ and rod length $k=7$. $1500$ snapshots are generated at each $\mu$ value. Equal number of data points are taken on the either side of the critical point $\mu_{c}=0.97$ for training. The snapshots below $\mu_{c}$ as labeled 0 and above $\mu_{c}$ are labeled 1. The data is divided into 3 parts. First the data is divided into 85% train set and 15% test set. The train set is further divided in two sets. The model is trained of 85% of the train set and validated on 15% of the train set.

We compare the accuracy of the phase classification problem around the I-N transition for the three models in Table. 1. It is evident that the CNN outperforms the logistic regression and DNN. This observation is intuitive as CNN is designed to capture the relevant features. Thus, for further study we use CNN with the same architecture and discuss the classification results in detail. We calculate the error bars using Bootstrapping.

To estimate the critical point of the system of hard rods we train neural networks on ground state snapshots of the ordered phase. This method is successfully used to estimate critical points of models like ferromagnetic and anti-ferromagnetic Potts models, 3D classical $O(3)$, 2D XY models Tan and Jiang (2020); Tan et al. (2020). Nematic phase configurations of the system size $182\times 182$ are used as the training set. The data are generated as mentioned in Sec. II by breaking the symmetry between horizontal and vertical rods upon introducing the variable $\Delta$. Snapshots are generated at $\mu=1.50$ and $\Delta=0.10$. At this value of $\Delta$, the system becomes almost fully ordered ($Q\approx 0.98$). All the vertical (horizontal) rods in the horizontal (vertical) rich snapshots are removed to obtain snapshots in which rods are completely aligned in one direction with order parameter $Q$ being $1$. The vertically aligned snapshots are labeled $[1,~0]$ and the horizontally aligned snapshots are labeled $[0,~1]$. Note that the nematic phase has two possible configurations– all horizontal and all vertical. CNN is trained using this data with one convolution layer followed by max pooling layer and two dense layers with softmax activation in the output layer, $L_{2}$ regularization is used to avoid overfitting. The loss, optimizer and activation used in training this model are categorical cross-entropy, adam optimizer, and relu activation. We also train Logistic regression with the physics-guided features as discussed above. In this case we work with $Q$ values $\leq 1$ (almost fully ordered). These trained models are then used to predict the labels of snapshots over a range of $\mu$ values. The norm of the predicted label, $R$, is calculated. The true value of $R$ can vary from 1 to $1/\sqrt{2}$ where $1$ corresponds to the fully ordered nematic phase and $1/\sqrt{2}$ corresponds to the case when the model fails to classify the phase as horizontal or vertical rich nematic phase. As the model is trained only with (almost) fully ordered phases, if one tests the output for a disordered phase, the output vector would ideally be $[0.5,~0.5]$ and the value of $R$ would be $1/\sqrt{2}$. The norm of the predicted labels of ordered states which are fully packed with rods in one direction should ideally be 1. To correct this, the difference between $1$ and the $R$ value of the predicted labels of the fully packed nematic phase, $\delta$ [where $\delta=1-(R_{\rm ver}+R_{\rm hor})/2$ and $R_{\rm ver}$($R_{\rm hor}$) is the norm of the predicted label of the fully packed nematic phase snapshot corresponding to vertical (horizontal) rods] is added to $R$. This $R$ is plotted against $\mu$ values to estimate the critical point. Assuming linearity of $R$ with $\mu$ near criticality, $\mu_{c}$ should be associated with the mid-point and is given by the intersection of the two curves $R+\delta$ and $1+1/\sqrt{2}-R-\delta$. As shown in the Fig 8, the intersection point of $R+\delta$ and $1+1/\sqrt{2}-R-\delta$ is the estimated critical point. The estimated values of $\mu_{c}$ by both the models are in close agreement (within error bars) with the numerical value $\mu_{c}=1.14\pm 0.03$ obtained from the probability distribution of the order parameter (see Fig. 2)

In addition, the modulus of the difference between two elements of the output vector $R$, can be treated as a Machine learned order parameter– the ordered phase will correspond to a value $\approx 1$, and in the isotropic phase it will be $\approx 0$.

In this work, we classify the phases of the system of hard rods on a square lattice. Although the classification task is trivial far from the transition point, sufficiently close to criticality system spanning fluctuations set in, making the phases visually indistinguishable and thus, the classification problem becomes harder. Three machine learning models, logistic regression, deep neural network, and convolutional neural network are trained to classify the phases around the isotropic-nematic transition. CNN has been shown to classify the phases with higher accuracy than the other two methods. We showed that the model performance improves with the increase in the system size– this is attributed to the finite-size effects. We further induce a first-order phase transition into the system using a field variable and the CNN is trained to classify phases around the transition. We demonstrate that classifying phases is easier around a first-order transition than around a second-order transition. Although the classification problem around the first-order transition is not the same as that around the second-order transition, our conclusions are not affected by that. We have also shown that physics-guided features drastically improve the performance of simpler models like logistic regression and random forest. In fact, with such feature engineering, they outperform more complex models like CNN (where the raw snapshots are used as inputs). We then estimate the critical point of the system of hard rods $\mu_{c}$ using only the information about the ordered phase. This estimate is in strong agreement with the value calculated by traditional methods. As discussed above, the Machine learned order parameter can be used to perform the finite size scaling analysis to compute the critical exponents. Our work infers that these methods may further be useful in studying transitions in more complex systems where defining an order parameter is nontrivial cite QCD. Another interesting question that emerges from our work is how to capture the critical correlations near a continuous phase transition using ML techniques. These will be future areas of investigation.