Entropy and heat capacity calculations of simulated crystal – hexatic smectic-B – smectic-A liquid-crystal phase transitions

Entropy and heat capacity are calculated in phase sequence of crystal, hexatic smectic-B, and smectic-A liquid crystals through constant pressure and temperature molecular-dynamics simulations of parallel soft spherocylinders. The transition from crystal to hexatic smectic-B phase is continuous while the transition to smectic-A phase is first order. The dependence of the phase sequence against the molecular shape anisotropy is investigated and there exists a triple point at a rather small anisotropy. Hopping diffusion of molecules is observed in the hexatic smectic-B phase.

Entropy S is a fundamental quantity in physics.In molecular simulations, S has been estimated through statistical physics approach ͓1͔ calculating S =−k B ͚ p B ln p B .The value of Boltzmann probability p B depends not only on each configuration but also on the entire ensemble through the partition function, thus the quality of sampling is always crucial.A different approach to calculate S through the absorption or production of heat in the system, i.e., thermodynamic approach, has been proposed recently ͓2͔.Precise calculation of the change in heat is possible by this method due to the structure of the Hamiltonian.
Phases with hexatic order, characterized by quasi-longrange sixfold bond orientational order ͑BOO͒, have been extensively studied in smectic phases ͑which are stacks of liquid layers͒ of liquid crystal compounds ͓3͔ and are relevant to many other systems including superfluids, superconductors, and liquid membrane to mention a few.Despite of the work accumulated to date, many aspects remain to be clarified in relation to the theory of defect-mediated melting transition ͓4͔.The hexatic smectic-B ͑HexB͒ phase can be considered as consisting of stacked two-dimensional ͑2D͒ hexatic layers, however, the correlations between the layers are also very important since three-dimensional ͑3D͒ sixfold BOOs are observed in these systems.The effect of the most basic characteristic of liquid crystal materials, i.e., anisotropy of the molecular shape, is investigated in systems of hardcore spherocylinders ͓5͔ as well as soft spherocylinders with rotation ͓6͔.However, in neither of these systems the HexB phase appears.By hindering the rotation of the long axis of soft spherocylinders, it has been shown that excluded volume effect alone can give rise to the HexB phase ͓7͔, however due to the limitation of the simulation methods, the investigation was quite limited.The system of parallel soft spherocylinders is the most primitive model to give rise to the HexB phase which is consistent with the fact that the orientational order of the long molecular axis is generally very high in HexB phase.
By thermodynamic approach, the relative entropy and heat capacity are calculated through the phase sequence of crystal ͑Cry͒, HexB, and smectic A ͑SmA͒ and are used to define the phase-transition temperatures.We obtain a phase diagram under constant pressure and temperature.In the SmA phase, no BOO exists.Not only the nature of the firstorder transition of HexB-SmA but also the continuous nature of Cry-HexB phase transition is clearly identified.In former studies, the identification of the Cry-HexB phase-transition boundaries has been done through analysis of bond orientational correlation functions ͓8͔ or dual criterion of molecular hopping rate and bond order ͓9͔.We also show that dynamical heterogeneity usually discussed in the context of metastable states, such as glass or supercooled liquid ͓10͔, also appears in HexB phase.
The model particles are parallel soft spherocylinders resembling that of Kihara ͓11͔ which express the interaction between anisotropic molecules through minimum distance between them.However, instead of using interaction with both repulsive and attractive parts, we use purely repulsive pairwise potential ͑just as Weeks et al. potential ͓12͔ is the purely repulsive part of the Lennard-Jones potential͒ where r 0 =2 1/6 D and R ij is the shortest distance between lines ͑where the potential is infinite͒ representing the long axis of the two spherocylinders, i.e., where ͑x ij , y ij , z ij ͒ is the relative position of the center of mass of the spherocylinders i and j, and L is the length of the line representing the long axis of the spherocylinder which is fixed parallel to the z axis.This corresponds to taking the moment of inertia as infinite.Translation of the spherocylinders is free.Periodic boundary conditions are applied in x, y, and z directions.Reduced units where length, energy, and mass are measured in D, ⑀, and m ͑where m is the mass of a spherocylinder͒, respectively, are used throughout this work.We mainly report results for systems with number of spherocylinders N = 1344 which form six layers.Calculations with N = 1170, 2340, 4680, and 9360 were also conducted to * aoki@icfd.co.jp check the system size effects.Time step of calculation is ␦t = 1.0ϫ 10 −4 for all the runs.Data for different temperatures are obtained by independent runs starting from identical initial configurations for each L. Each data is an average value for time duration ⌬t = 100 ͑i.e., 1 ϫ 10 6 time steps͒ at the end of t Ն 1100 calculation unless otherwise stated.The masses of the barostat and thermostat are, respectively, M = 0.01 and K = 20 000 for most of the calculations.
The molecular-dynamics ͑MD͒ method we use was designed to simulate soft matter effectively ͓2͔ by using an explicit symplectic integrator ͓13͔.The barostat is designed to conform to the changes in both shape and bulk elasticities by effectively dealing the anisotropic fluctuations ͓14͔.The thermostat is a Hamiltonian version of Nosé-Hoover dynamics ͓15͔.The entropy is calculated from the change in heat ⌬Q, i.e.,

⌬S = ⌬Q
where is the time scaling factor of the Poincaré time transformation and W is the work ͓2͔.Especially the function ͑1−͒ / T is useful in identifying the nature of the phase transition, since W is a monotonically slow varying function quite small compared to the value of the original Hamiltonian itself.The change in the value of ͑initially set to =1͒ is reflecting the absorption or production of heat in the system.Phase diagram under constant pressure P = 1000 of parallel spherocylinders with various lengths L is given in Fig. 1.Triple point of Cry, HexB, and SmA phases appears at 1.0 Ͻ L Ͻ 1.05.This length is close to the value where HexB phase has been previously reported ͓7͔.As the length of the spherocylinder gets longer, Cry-HexB transition temperature T CB decreases.This tendency coincides with that of crystalhexatic transition temperature of membranes confined in two parallel plates ͑see Fig. 1.6 in ͓4͔͒.The transition temperature to SmA phase ͑independent whether it is from Cry or HexB phase͒ is linearly dependent on the molecular length L. The linear dependence of crystal-smectic phase-transition temperature on molecular length has also been observed in soft parallel spherocylinders interacting via inverse power potential ͓16͔.
In   is evident that the ͑1−͒ / T term ͓Fig.3͑a͔͒ shows the characteristics of the phase transitions clearly, i.e., the continuous nature at T CB and the first-order nature at T BA .Especially, the usefulness of function ͑1−͒ / T to determine T CB is clearly demonstrated since there exists an extremum at T CB whereas there is only a kink in the entropy curve.In the case of spherical soft particles, the function ͑1−͒ / T has been almost proportional to the ⌬S / N curve near the crystalisotropic liquid phase transition ͓2͔.For shorter parallel spherocylinders where crystal melts directly to SmA phase ͑L Յ 1͒, we observed a similar behavior in ͑1−͒ / T.
The differential dH / dT of enthalpy with respect to temperature is shown in Fig. 4.This quantity corresponds to the curve obtained by differential scanning calorimetry ͑DSC͒ in experiments.However, the thermodynamic properties reported here are free from effects of cooling or heating rates since independently equilibrated calculations are conducted at each temperature.
Figure 5 shows the molar heat capacity C P calculated for the same system as Figs. 2 and 3.The value of C P is obtained via the fluctuation of the relative entropy ⌬S ͓19͔, i.e., C P ϰ ͑͗S 2 ͘ − ͗S͘ 2 ͒, thus expressed in arbitrary units.In Fig. 5, not only peaks at T CB but also a step in the value at T BA can be discerned.The double peak at T CB may result from fluctuation induced force near the critical point confined in periodic boundary conditions ͓20͔.Anomalous fluctuation at low temperatures, as reported in ͓21͔, was also observed in our simulations.To get reasonable values of second-ordered thermodynamic quantities such as C P required slow relaxation at lower temperatures with weak coupling to the thermostat, especially in the crystal phase ͑t Ն 4000 in simulation time using mass of thermostats kinetic energy K = 20 000 ͓2͔͒.
When the system was forced to relax fast by strong coupling to the thermostat ͑through K =20͒, reasonable average values of first-order thermodynamic quantities were obtained but resulted in larger values of fluctuation which do not decay, thus, not leading to reasonable values of second-order thermodynamic quantities.It has been reported that strong coupling to the thermostat results in cutting off fluctuations of low frequencies ͓22͔, which must influence the second-order thermodynamic quantities.These results clearly show that it is important to relax the system with weak coupling to the thermostat to get proper values of second-order thermodynamic quantities.
Figure 6 shows the sixfold bond orientational order parameter where j runs over all molecules in the system, k runs over all the nearest-neighbor molecules of j, kj is the angle between the fixed x axis and the projection of the line connecting the center of mass of particles k and j on the xy plane, and N b denotes the total number of bonds in the system.The bracket ͗ ͘ denotes the time average of duration ⌬t =10 ͑1 ϫ 10 5 time steps͒.Note that C 6 contains the information of all layers in the system.The inset shows how the system size affects the value of C 6 .It can be seen that the number of layers in the  Although the specific volume and C 6 vary continuously at T CB , the dynamics drastically change.At temperatures above T CB , hopping dynamics is observed.In HexB phase, diffusion of the particles inside each layer ͑in xy plane͒ increases as the temperature rises and shows a high value of diffusion at temperatures near T BA .Fluctuations of the layer positions are also observed in the HexB phase, however, diffusion of molecules among different layers ͑z direction͒ only occurs at temperatures much higher than T BA .In Fig. 8, positions of all the molecules ͑projected on xy plane͒ relative to that at a certain time are plotted.If the molecule does not move for the period of time, a dot will be plotted on ͑0,0͒ for such diagrams.In the crystal phase, dots are crowded near ͑0,0͒ since molecules only fluctuate around their equilibrium positions.In HexB phase, hopping to other vacant sites occurs as shown in Fig. 8 and dynamical heterogeneity is clearly observed.There are patches of hexagonal lattice which are out of position with each other.In addition to the molecular hopping, the hexatic patches themselves diffuse in the long run which is the reason for not having many points near ͑0,0͒.Since Fig. 8 contains data of all the molecules in six layers, it also shows that correlation between the layers exists.The hopping dynamics in HexB phase was also observed in Monte Carlo simulations and was used as a criterion to define the boundary between Cry and HexB phases ͓9͔ along with the value of C 6 .From our simulation, it is clear that at temperatures close to T CB ͑identified from relative entropy͒, the hopping rate is negligibly small and has intermittent characteristics, thus difficult to get reliable statistics even in the longest run.We also have seen in the previous paragraph that the values of C 6 cannot be a reliable quantity either to predict the transition temperatures, especially T CB .
MD simulations have been conducted for a model of liquid crystals and the nature of Cry-HexB-SmA phase transition has been investigated.The transition temperatures were identified by relative entropy and heat capacity.In the HexB phase, hopping dynamics has been observed.It has been shown that the high-precision symplectic integrator designed for simulating soft matter is effective in investigating not only the dynamics but also the thermodynamics of the system.
FIG. 1. Phase diagram of N = 1344 parallel soft spherocylinders at P = 1000 for various lengths L. Transition temperatures are identified from function ͑1−͒ / T consisting the relative entropy for Cry-HexB phase transition ͑᭺͒ and transition to SmA phase ͑b͒.
FIG. 5. Molar heat capacity C P in arbitrary units ͑obtained from fluctuation of relative entropy͒ vs temperature T for N = 1344 parallel soft spherocylinders with L =2 at P = 1000.Inset shows C P around T CB .