Simulating the nematic-isotropic phase transition of liquid crystal model via generalized replica-exchange method

The nematic-isotropic (NI) phase transition of 4-cyano-4'-pentylbiphenyl (5CB) was simulated using the generalized replica-exchange method (gREM) based on molecular dynamics simulations. The effective temperature is introduced in gREM, allowing the enhanced sampling of configurations in the unstable region, which is intrinsic to the first-order phase transition. The sampling performance was analyzed with different system sizes and compared with that of the temperature replica-exchange method (tREM). It was observed that gREM is capable of sampling configurations at sufficient replica-exchange acceptance ratios even around the NI transition temperature. A bimodal distribution of the order parameter at the transition region was found, which is in agreement with the mean-field theory. In contrast, tREM is ineffective around the transition temperature owing to the potential energy gap between the nematic and isotropic phases.


I. INTRODUCTION
Highly anisotropic molecules can form a variety of mesophases, including nematic, smectic, and columnar phases between crystalline solids and isotropic liquids. 1 Such molecules exhibiting liquid crystal (LC) phases are referred to as mesogens. One of the common mesogens is 4-cyano-4'pentylbiphenyl (5CB), which undergoes a nematic-isotropic (NI) phase transition at room temperature. An order parameter for characterizing the nematic LC phase is introduced with respect to the director, which is defined as a collection of mesogen molecules. The order parameter reduces with an increase in temperature. Typically, it reduces from 0.6 to 0.4 on approaching to the transition temperature and exhibits a discontinuous drop to zero at the transition. This behavior is regarded as a sign of the first-order phase transition, as demonstrated well by the mean-field theory for the NI transition. [2][3][4] Because the time scale associated with the equilibration of the nematic phase becomes longer with an isotropic starting configuration, it is still challenging to perform molecular dynamics (MD) simulations of the NI phase transition from atomistic levels, even though considerable research has been conducted on the use of all-atom, united-atom (UA), and coarse-grained models. [5][6][7][8][9][10][11][12] In particular, MD simulations of systems exhibiting first-order phase transitions often encounter difficulties in sampling the configurations around the transition temperatures, owing to the presence of unstable states bridging the potential energy gaps between two stable phases.
The temperature replica-exchange method (tREM) (also known as parallel tempering) is a promising sampling method for simulating the phase transitions in various physical and chemical systems. 13,14 The tREM is developed to sample a a) Electronic mail: kk@cheng.es.osaka-u.ac.jp b) Electronic mail: nobuyuki@cheng.es.osaka-u.ac.jp wide range of configurations, where many replicas at different temperatures are simulated in parallel and the configurations between two replicas are exchanged to prevent replicas at lower temperatures from staying in the local minimum state. However, the exchanges of replicas become ineffective near the first-order phase transition temperature, which results in tREM inefficiency. It should be noted that various enhanced sampling methods beyond the original tREM have been developed. 15 Berardi et al. 16 and Kowaguchi et al. 17 attempted to improve the sampling efficiency near the NI phase transition with the Gay-Berne model and anisotropic Lennard-Jones fluids using the multidimensional REM 18 (also known as Hamiltonian REM 19 ) and the isobaric-isothermal REM, 20 respectively.
Recently, Kim et al. proposed the generalized replicaexchange method (gREM), 21 which was designed to efficiently simulate first-order phase transition systems. In gREM, the effective potential and effective temperature are introduced for each replica to ensure the unimodal distribution of the potential energy E of the system, even in the thermodynamically unstable region. The effective potential is derived from the ensemble weight, which is set to the Tsallis form in practical implementations. Further, by utilizing its dependence on the most probable value of E, the effective temperature is connected to the statistical temperature of the system. Sufficient overlaps were achieved for the probability distributions of E for replicas with different effective temperatures, which enabled efficient sampling of the transition region. The gREM has been applied to various systems showing solidliquid 22,23 and vapor-liquid 24,25 transitions.
The isobaric extension of gREM was proposed by Małolepsza and Keyes. [26][27][28][29] In this extension, the enthalpy H is used instead of the potential energy E, which is employed in the canonical case (to be exact, H refers to the sum of E and the product of pressure and volume). The effective temperature for each replica is introduced with respect to H, and owing to the use of H, gREM is appropriate for simulating constantpressure ensembles. The isobaric gREM was also utilized to investigate the liquid and vapor phases of water 30 and the order-disorder phase transitions in lipid bilayer systems. 31 Therefore, it is of interest to assess the applicability of isobaric gREM to the NI phase transition, which is characterized by the orientational ordering. In this study, we performed MD simulations in combination with isobaric gREM for 5CB based on the united-atom (UA) model developed by Tiberio et al. 32 Finite-size effects on the NI transition are important because an equilibration time scale near the phase transition region drastically increases with increasing the number of molecules N. In computer simulations for LC systems, the finite-size effects were indeed examined, for example, using the Lebwohl-Lasher model [33][34][35] and anisotropic Lennard-Jones fluids. 36 We investigated system size effects by the use of gREM with the number of molecules N ranging from 250 to 4000. The sampling performance was also analyzed by comparison with tREM. From sampled configurations, we examine the temperature dependence of the density and order parameter and discuss the temperature and enthalpy relationship.

II. MODEL AND SIMULATION DETAILS
The UA model for 5CB developed by Tiberio et al. was utilized in this study. In this model, the generalized AM-BER force field (GAFF) was modified to reproduce the temperature of the NI phase transition in the experiments, T NI ≈ 308.2 K. 32 The comparison with experimental results of the density and order parameter was also reported in Ref. 32. This UA model has also been employed in previous studies, [37][38][39] which considered the NI phase transition in a consistent manner with Ref. 32. Another UA model was developed by the parametrization of the TraPPE-UA force field for the 5CB molecule. 40 Moreover, the transferability of the coarsegrained model has also been proposed. 41,42 The present work is mainly methodological, and Tiberio et al.'s model was employed since it reflects the chemical reality and reproduces the NI transition temperature well.
The simulated system was composed of 5CB molecules in a cubic box with periodic boundary conditions. The number of molecules was varied as N = 250, 1000, 2000, and 4000. The pressure was set to 1 atm, and the temperatures were ranged between 300 K and 320 K. The NPT ensemble with the Nóse-Hoover thermostat and isotropic Parrinello-Rahman barostat was used with a time step of 2 fs in all the simulations. We performed gREM and tREM using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 43 .
The basic idea of gREM is summarized as follows (see the detail in Ref. 21): The effective temperature T α of replica α (=1, 2, · · · , M) is introduced by an inverse mapping of the effective potential w α , T α = (∂w α /∂E) −1 . w α is related to the ensemble weight W α through the relation, w α = − ln W α . T α is set to intersect the statistical temperature, T s = (∂S /∂E) −1 .
Here, S and E denote the configurational entropy and potential energy in the canonical ensemble, respectively. Furthermore, the parametrization with a linear function for the effec- , provides the ensemble weight that is equivalent with the form of the Tsallis statistics, Here, E 0 represents an arbitrarily chosen potential energy and γ is chosen to have a sufficiently large negative value to ensure that T α (E) intersects T s (E) only once. Note that λ α is a control parameter for determining the distribution of E in replica α for a given slope γ.
The non-Boltzmann ensemble with the weight W α and effective temperature T α enables to sample the thermodynamically unstable states bridging two stable phases. The isobaric gREM is formalized considering the enthalpy H, statistical temperature T s = (∂S /∂H) −1 , and effective temperature T α = (∂w α /∂H) −1 for replica α. 26 In isobaric gREM, the enthalpy dependent effective temperature for replica α is given by where H is the sum of the potential energy and the product of pressure and volume, and H 0 and γ represent the reference enthalpy and the slope of T α (H), respectively. Again, γ and λ α are control parameters of the temperature intercept at a chosen enthalpy H 0 . The acceptance probability of replica exchange between neighboring replicas α (with enthalpy H) and α (with enthalpy H ) is given by In contrast, the acceptance of replica exchange in tREM is given by where β α = 1/T α and β α = 1/T α are inverse temperatures of replica α and α , respectively. The volume is equilibrated with the barostat of 1 atm if the replica exchange is occurred.
In the present simulations, M = 11 replicas were utilized for both gREM and tREM. The lowest and highest temperatures were 300 K and 320 K; hence, temperatures of replicas 1 and M were T 1 = 300 K and T M = 320 K, respectively. The values of enthalpy at the two temperatures (denoted byH 1 and H M , respectively) were quantified from the averages of 10 ns simulations after equilibration over 50 ns using the NPT ensemble. For supplementary comparison, we also simulated isobaric-isothermal REM by Okabe et al. 20 (2) with the pressure and volume pairs (p α , V α ) and (p α , V α ) of replica α and α , respectively. The second term is designed to employ the enthalpic quantity for generation of an isobaric-isothermal ensemble.
For gREM, the slope γ was determined from γ = (T M − T 1 )/(H 1 −H M ). The values ofH 1 ,H M , and γ are listed in Table I. The reference enthalpy H 0 in Eq. (1) was set toH 1 in the present gREM. It is also seen thatH 1 /N andH M /N are insensitive to the system size; therefore, γ varies roughly in inverse proportion to N. The control parameter λ α for replica α is given by In contrast, for tREM, the temperature of replica α is given by T α = T 1 + ∆T (α − 1) with ∆T = 2 K. A 50 ns simulation was run, in which a temperature swap between adjacent replicas were attempted every 10 ps both for both gREM and tREM. We calculated various quantities presented in Sec. III from the last 10 ns simulation containing 1000 configurations, where the stationary state was achieved. Since all the events of temperature swap were tracked, the data can be extracted for each replica in both gREM and tREM.
Note that the most probable enthalpy H * in gREM provides the relation, T s = T α (H * ). 26 Henceforth, T of gREM is represented by the temperature determined from the peak of the enthalpy distribution for each replica from Eq. (1) (see also Fig. 4(a) below). In contrast, T of tREM is represented by the designed temperature of each replica.

III. RESULTS AND DISCUSSION
First, we investigate the temperature dependence of the density. Figure 1 shows the results for gREM (a) and tREM (b) by changing the number of molecules N. The sudden drop of the density is evident at around 305-310 K, particularly for the larger system size with N = 4000. The transition from nematic to isotropic phases occurs in this temperature range, as demonstrated in previous studies using the same UA model. 32,[37][38][39] We next examine the orientational order parameter of the LC nematic phase. The order parameter P 2 is generally expressed by the second-order Legendre polynomial as follows: where θ i represents the angle between the director n and the unit vector of the long axis of molecular i, u i . Here, the director n is the unit vector representing the preferred direction of the local volume. The brackets denote the an ensemble averages. The most widely used method to quantify P 2 from MD simulations is the diagonalization of the order parameter tensor Q. 44,45 The expression of Q is given by with the unit matrix I. The eigenvalues of the matrix Q, λ − < λ 0 < λ + , guarantee λ − + λ 0 + λ + = 0 because Q is traceless. The largest eigenvalue λ + provides the order parameter P 2 , and the corresponding eigenvector is the director n. In practice, P 2 was calculated as −2λ 0 , which was suggested for a better estimation of P 2 , particularly in the isotropic phase. 46 Here, u i in this work refers to the C≡N bond of the cyano group. Note that the dependence of the choice of u i on P 2 was negligible if the inertia axis of the molecule is used. 32 Figure 2 shows the temperature dependence of the order parameter P 2 from gREM (a) and tREM (b) simulations. It is observed that the decrease in temperature T leads to a decrease in P 2 from 0.5-0.6 to nearly zero for all the system sizes that were investigated. The temperature dependence of P 2 exhibits a sharper drop across the NI transition for the larger system with N = 4000. However, P 2 for the smallest system size with N = 250 gradually decreases with increasing T . In general, the discontinuity of the first-order phase transition between two stable states are smeared out in finite systems. 47,48 The N dependence in Fig. 2 is in accordance with the well-known fact that unstable regions of the phase transition become narrow as the system size decreases.
Here, we used the empirical equation, 49,50 for fitting with the order parameter P 2 . We obtained the residual order parameter under the isotropic phase P iso 2 ≈ 0.12, NI transition temperature T NI ≈ 307 K, and pseudo-critical exponent β ≈ 0.20 for the gREM result with N = 4000. The fitting result is shown in Fig. 2(a). Note that T NI ≈ 307 K is slightly lower than the experimental value of 308 K. 51,52 In contrast, the fitting for tREM with N = 4000 becomes uncertain due to the unsampled value of P 2 around the NI transition temperature, as observed in Fig. 2(b). The fittings of P 2 with N = 250 using Eq. (5) are also shown in Fig. 2 (a) and (b). For both gREM and tREM, the NI transition temperature is estimated as T NI ≈ 310 K. This value is in good agreement with the NI temperature reported in Ref. 32 The distribution of the orientational order parameter P 2 is also important to examine the NI transition. 32 Figure 3 shows the results from replicas 1, 5, and 11 for the gREM with N = 4000. The corresponding temperatures T α (H * ) of gREM are 300 K, 307.0 K, and 320 K, respectively. The value of 307.0 K corresponds to T NI , which was determined from the fitting results using Eq. (5). The distributions of P 2 in replicas 1 and 10 are unimodal, whereas two peaks at P 2 ≈ 0.1 and 0.4 are observed in replica 5, resulting in the large fluctuation of P 2 at 307.0 K, as shown in Fig. 2(a). The bimodal distribution of the order parameter P 2 is consistent with the free energy landscape predicted by the mean-field theory for first-order phase transitions 1 . However, the theoretical exponent β = 0.5 is different from the simulation result. In addition, it is interesting to examine pretransitional behaviors in the isotropic phase because the NI transition is relatively weak first-order phase transition. This regard can be characterized by the short-ranged orientational order from the orientational correlation functions G 1 (r) = δ(r − r i j )(u i · u j ) i j and 32 where r i j denotes the distance of center-of-mass between molecules i and j. We also calculated G 1 (r) and G 2 (r) using the gREM with N = 4000 and obtained results consistent with those reported in Ref. 32 (data not shown).
It should be noted that the difference between gREM and tREM becomes significant for the system with N = 4000 around T NI (see Figs. 1 and 2). Although the reduction of density and P 2 with temperature was also reproduced using the tREM, the values of density and P 2 from the tREM and gREM were not in agreement in the transition region. In particular, the averaged P 2 at 306 K (replica 4) of the tREM in Fig. 2(b) is close to the value in the isotropic phase with a large error bar. To elucidate the difference between the tREM and gREM with N = 4000, we quantified the acceptance ratios of the replica exchanges, the results of which are sum-TABLE II. Acceptance ratio of the replica exchanges during the last 10 ns simulation between two replicas α and α + 1 for the gREM and tREM with N = 4000. Note that * indicates that no exchanges were observed during 1000 trials.  marized in Table II. The acceptance ratio in gREM is larger than that in tREM for any pair of replica exchanges. The exchange ratio in gREM reduces between replicas 4 and 5, near T NI ≈ 307 K. Nevertheless, it exceeds 10% and assures the efficiency of the replica-exchange method. This is clear evidence that gREM is an effective simulation method, even near the phase transition temperature for larger system sizes. By contrast, the replica exchanges between replicas 4 (T = 306 K) and 5 (T = 308 K) and between replicas 5 (T = 308 K) and 6 (T = 310 K) are significantly restricted in tREM, which is an effect of the phase transition, as anticipated beforehand. We confirmed that the isobaric-isothermal REM relatively en- hances the replica exchange between neighboring replicas, but the acceptance ratio between replicas 4 and 5 remains small at 0.2%. This is due to the fact that the volume change of the NI transition is small about 4% from 300 K to 320 K (see Fig. 1). The acceptance ratios between neighboring replicas are determined by the degree of overlap in the enthalpy or potential energy distributions in gREM or tREM, respectively. Figure 4 shows the probability distribution of enthalpy in gREM (a) and potential energy in tREM (b) for each replica. The transition region corresponds to replicas 4, 5, and 6. The overlap of the distributions is significant in this region when gREM is employed. However, the overlap becomes scarce when using tREM, leading to a significant reduction in replica-change events. In principle, the acceptance ratios of replica exchanges in tREM can be improved by using more replicas to increase the degree of overlap of the potential energy distribution. However, simulating the tREM with more replicas is impractical in terms of computational costs, particularly for significantly larger system sizes. In contrast, for the smallest system size with N = 250, we confirmed that replica exchanges of tREM remain even at the NI transition region, showing behaviors in the density and P s similar to those of gREM, as seen in Figs. 1 and 2.
Finally, we examine the relationship between the temperature T and most probable enthalpy H * , which is plotted in Fig. 4(a). Here, T of gREM was determined from T s = T α (H * ) using Eq. (1) with the peak value H * of the enthalpy distribution for replica α (see Fig. 4(a)). It has been demonstrated that T becomes flat when crossing the NI transition temperature (T NI ≈ 307 K). For comparison, the temperature as a function of H * from tREM simulations with N = 4000 is also shown in Fig. 5(a). Although T (H * ) of tREM also exhibits a curve similar to that of gREM, the plotted points become sparse, particularly for the NI transition region. The temperature dependence of the heat capacity C p = ∂(H * /N)/∂T with N = 4000 can be further evaluated from a central difference approximation for H * /N as a function of T , as plotted in Fig. 5(b). The sharp peak associated with the NI phase transition was observed at T NI ≈ 307 K for gREM. This peak is also consistent with the drop of the density at the NI transition with N = 4000, as observed in Fig. 1(a). In contrast, the peak of C p becomes unclear for tREM, again indicating the necessity of more replicas to increase the resolution around the NI transition. For smaller system sizes, the inflection of the T − H * /N curve is gradually smeared out, causing smaller peaks of C p for both gREM and tREM (data not shown).

IV. CONCLUSIONS
In this paper, we report on MD simulation results combined with the gREM for the NI phase transition of the 5CB model. We demonstrated that the gREM can be applied to the NI phase transition of LC systems. It is effective in the sense that the acceptance ratios of the replica exchange remain significant at the NI transition temperature even with increasing N. By contrast, the replica exchange of tREM became inefficient, particularly for larger system sizes, although tREM simulated the results similar to those with gREM for N = 250.
The temperature dependence of the density and orientational order parameter P 2 was also examined. In particular, P 2 in the transition region between the nematic and isotropic phases can be sampled well using the gREM with the help of the effective temperature. A sharp drop in P 2 toward the NI temperature was observed, and this was more prominent for larger system sizes. In addition, a sharp peak of heat capacity C p was clearly observed around the NI temperature.
An advantage of gREM is its availability for the smectic phase of 4-octyl-4'-cyanobiphenyl (8CB) using the relevant UA model. 53 Note that a recent work shows the free energy landscape of smectic-nematic phase transition using machine learning technique. 54 It is interesting to examine these complex phase transition behaviors using gREM. Moreover, it is important to investigate the applicability of gREM in all-atom MD simulations of LC systems, for example, self-assembling helical structures 55 and nanochannels, 56 where phase transitions occur continuously over a wide temperature range. Further studies focusing on these aspects are required.