Liquid crystal behavior of the Kihara fluid

Liquid crystal behavior of the Kihara fluid A. Cuetos, B. Martı ´nez-Haya, and S. Lago Departamento de Ciencias Ambientales, Universidad Pablo de Olavide, 41013 Seville, Spain L. F. Rull Departamento de Fı ´sica Ato ´mica, Molecular y Nuclear, Area de Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado 1065, 41080 Seville, Spain 共Received 28 February 2003; published 15 July 2003兲 The liquid crystal phases of the Kihara fluid have been studied in computer simulations. The work focuses on the isotropic–nematic–smectic-Atriple point region, especially relevant for the understanding of the properties and the design of real mesogens with specific phase diagrams. The Kihara interaction resembles more appropriately than other related models, the shape of elongated polymers and biomolecules, and a closer assertion is provided for the role of the configurational entropy and the dispersive interactions in the behavior of such molecules in dense phases or under macromolecular crowding conditions. DOI: 10.1103/PhysRevE.68.011704 PACS number共s兲: 64.70.Md, 61.30.Cz, 87.15.Aa I. INTRODUCTION A number of computationally efficient pair interaction models have been proposed over the past decades for the study of liquid crystalline mesogens. In particular, fluids composed of elongated 共prolate兲molecules constitute an important class of liquid crystals with relevant industrial and biological applications. Therefore, different models have been used to explore their behavior such as the hard ellipsoid 共HE兲fluid, the Gay-Berne 共GB兲fluid, or the hard spherocylinder 共HSC兲fluid and its square–well 共SWSC兲or soft repulsive 共SRS兲variants 关1–9兴. The Kihara fluid model was introduced in the 1960s 关10兴 as a generalization, for anisotropic molecules, of the spherically symmetric Lennard-Jones fluid 关see Eq. 共1兲兴 and has been thoroughly employed to investigate the thermodynamic, structural, and transport properties of fluids of linear molecules. In particular, this model has been found to be appropriate to resemble the liquid-vapor equilibrium of different solvents and alkanes of moderate elongation 关6兴. The Kihara interaction potential incorporates a soft 共Lennard-Jones type兲 attractive well around the molecular core, which is absent or oversimplified in other models investigated previously, such as the SRS or SWSC fluids. The flat square-well attraction of the SWSC fluid, for instance, allows for small compressions or expansions of the system without change of energy. In addition, the spherocylindrical shape of the molecules in the Kihara model is more realistic for many 共bio兲chemical systems than the ellipsoidal shape imposed by the well-known HE or GB models. For instance, recent investigations have explored the relevance of the spherocylindrical core in the stability of microtubules in solution 关11兴. We present here a computational study of the Kihara fluid focused on the stability of the isotropic 共disordered兲, nematic 共orientational order兲, and smectic 共layering order兲liquid crystal phases at different temperatures. To our knowledge, this is the first time that the liquid crystal behavior of this model fluid has been explored, since all of the previous investigations have been restricted to moderate molecular elongations and/or packing fractions corresponding to thermodynamic states deep within the isotropic phase of this system. We hope that the present results will aid in the assessment of the roles of entropy and of the dispersive interactions in the stability of ordered phases of real mesogens and into the appropriate strategy to model them. II. METHOD We have applied the isothermal-isobaric (NPT) ensemble Monte Carlo 共MC兲technique to evaluate thermodynamic and structural properties of a Kihara fluid with elongation L* ⫽L/ ␴ ⫽5. Within this model, the molecules are represented by rigid cylinders of diameter ␴ with semispherical caps 共i.e., a so-called spherocylinder兲and the pair interaction potential is given by U共r,⍀兲⫽ 再 4 ⑀ 关共 ␴ /dm兲12⫺共 ␴ /dm兲6兴,dm⭐3 ␴ 0, dm⬎3 ␴ .共1兲 The minimum distance between the central axes of the two prolate molecules, dm⫽dm(r,⍀), is a function of their relative orientation ⍀, determined by three independent angles, and of the distance between their centers of mass, r. For the present simulations, the interaction was truncated at a distance dC⫽3 ␴ , which corresponds to center-of-mass distances ranging from r⫽3 ␴ to r⫽8 ␴ , depending on the relative pair orientation. The MC simulations were run over a system of Np ⫽768 molecules at four reduced temperatures, T*⫽kT/ ⑀ ⫽1.5, 2, 3, and 5 (kdenotes the Boltzmann constant and ⑀ denotes the well depth兲and started from an equilibrated high pressure state well inside the smectic A region of the phase diagram. After equilibration (1⫻106MC cycles兲and averaging (3⫻105MC cycles兲, the last molecular arrangement is stored and used as initial configuration for the subsequent run with a smaller system pressure. Each MC cycle consists of Npattempts for random displacements and/or reorientations of the particles, plus an attempt to change the box volume. The usual periodic boundary conditions are employed and the simulation box is anisotropic with fixed ratios PHYSICAL REVIEW E 68, 011704 共2003兲 1063-651X/2003/68共1兲/011704共4兲/$20.00 ©2003 The American Physical Society68 011704-1 Ly/Lx⫽1.15 and Lz/Lx⫽1.94. As in our previous study of the SWSC and SRS fluids 关5兴, the longer dimension along the zdirection is meant to increase the number of possible smectic layers with respect to the cubic box for a given number of molecules. Further methodological details can be found in Ref. 关5兴. The 共first-order兲liquid crystalline transitions observed in the expansion of the fluid are characterized by discontinuities in the density, the nematic order parameter, the bond hexagonal order parameter, and by sudden qualitative changes observed in the different correlation functions 关3,5,12兴. No trace of the hexagonal ordering or long-range correlations characteristic of smectic-Bor solidlike phases were found within the scope of the present work. III. RESULTS The equations of state in reduced units (P*⫽P ␴ 3/kT vs ␳ *⫽ ␳␴ 3) of the L*⫽5 Kihara fluid along the four isotherms investigated are represented in Fig. 1. In addition, Fig. 2 and Table I summarize the coexistence pressures and densities for each of the liquid crystal phase transitions observed. At the highest temperature investigated, T*⫽5, the Kihara fluid presents a stable smectic-Aphase that undergoes a first-order transition to a nematic phase, as the system is expanded to the pressures below P*⫽1.8. A further expansion of the system induces a transition to an isotropic phase at a pressure just below P*⫽1.35. A qualitatively similar behavior is found in the expansion at T*⫽3. However, at this lower temperature the nematic phase presents a reduced range of stability, as the transition from smectic-Ato nematic is delayed to lower pressures and densities in comparison to T*⫽5. Also noticeable is the increased difference in the density of both phases at coexistence. On the other hand, although the subsequent transition to the isotropic phase occurs at a greater reduced pressure in the expansion along the T*⫽3 isotherm than along the T*⫽5 one, the nematicisotropic coexistence densities are very similar at both temperatures. At the lower temperatures, T*⫽1.5 and T*⫽2, the phase diagram of the Kihara fluid changes qualitatively. The most relevant feature is the absence of the nematic phase and, hence, the appearance of a smectic-A–isotropic transition. This kind of phase transition features a greater density asymmetry in the coexistent phases and therefore a stronger firstorder character than the smectic-A–nematic or isotropicnematic transitions. Interestingly, the density difference between the smectic-Aand isotropic phases at coexistence is greater at T*⫽1.5 than at T*⫽2. This behavior has been also found for the Gay-Berne fluid 关2,12,13兴and appears to be characteristic of systems with phase diagrams in which the temperature of the I–N–Sm-Atriple point is greater than the critical temperature, which for this fluid is significantly smaller than T*⫽1关14兴. The coexistence states derived from our simulations for the L*⫽5 Kihara fluid allow for an estimation of the triple-point temperature of roughly T* ⬇2.5, as tentatively indicated in Fig. 2. Figure 1 compares the equations of state and phase transitions at T*⫽5 for the Kihara fluid and its repulsive counterpart, the SRS fluid 关5兴. The purely repulsive SRS interac- FIG. 1. Equations of state for the isotherms T*⫽1.5, 2, 3, and 5 of the Kihara fluid with elongation L*⫽5. The vertical lines denote the coexistence densities at the isotropic-nematic (I-N), isotropic– smectic-A(I–Sm-A) or the nematic–smectic-A(N–Sm-A) transitions of the fluid. The top panel shows the equation of state of the SRS fluid of same elongation at T*⫽5共Ref. 关5兴兲 for comparison. FIG. 2. Summary of the coexistence densities at the isotropicnematic (I-N), isotropic–smectic-A(I–Sm-A), or nematic– smectic-A(N–Sm-A) transitions obtained in the present MC-NPT simulations of the Kihara fluid. The open square represents a tentative extrapolation of the results to the I-N–Sm-Atriple point. CUETOS et al. PHYSICAL REVIEW E 68, 011704 共2003兲 011704-2 tion potential constitutes a generalization of the Weeks- Chandler-Andersen potential 关15兴and is built from a shifted Kihara potential truncated at the well minimum. Due to the upward shift of the repulsive wall with respect to the Kihara one, the effective size of the molecules at a given temperature is larger in the SRS model 关4,5兴. Hence, the SRS fluid is expected to be less compressible and to be affected by entropic excluded volume effects to a greater extent than the Kihara fluid, which explains the shift of both the isotropicnematic and the nematic–smectic-Atransitions toward lower densities in the SRS system with respect to the Kihara one at T*⫽5共see Fig. 1兲. The argumentation is as well consistent with the location of these liquid crystal transitions at even lower densities in the rigid HSC and SWSC fluids 关5兴,ascan be seen in Table I. These latter considerations indicate that, in general terms, the main ‘‘driving force’’ of the isotropic-nematic-smectic transitions in the fluids studied in this work is of entropic nature, the system trying to maximize the volume of phase space accessible to the molecules. The role of the energetic contribution to the free energy is a more subtle, yet relevant one. The most significant effect of the pair interaction potential energy is linked to entropy, since it is the nature of the short-range repulsive interaction that imposes the effective size of the molecules and, thus, the volume excluded by the core of the molecules. Another important qualitative aspect is the role played by the presence of the attractive well in the pair interaction. A first approach to this question is given in our previous work on the SWSC fluid 关5兴, where it was shown that the square-well interaction induces a delay in the I–N–Sm-Atransitions with respect to the purely repulsive HSC fluid 共Table I兲. Thus, the effect of the attractive well affects the transitions in a similar way as a reduction of the effective molecular size would do. Furthermore, the influence of the attractive well should be more relevant at low temperature, but this aspect has not been much explored for fluids of spherocylindrical symmetry 关5兴, even though it is well characterized for the Gay-Berne fluid 关2,9,12,13兴.A glance at the energetic constrains on the liquid crystal phase transitions as a function of temperature is provided in Fig. 3, which shows the evolution of the average potential energy per particle along the four isotherms of the Kihara fluid studied in this work. Interestingly, the qualitative behavior of the energy at T*⫽5 is quite similar to the one found for the SRS fluid at the same temperature 关5兴共also shown in Fig. 3 for direct comparison兲. In both of these cases, the energy grows monotonously with density, except at the nematic–smectic-A transition where a slight decrease in energy is observed. This qualitative trend changes at the lower temperatures, for which a decrease in energy is eventually observed as the density grows within the smectic phase. Hence, it appears that at low temperature the smectic phase is favored by both entropic and energetic constrains. In fact, the increased stabilization of the smectic phase opposes that of the nematic phase, which contributes to the disappearance of this latter phase at T*⭐2. TABLE I. Coexistence pressures and densities of the liquid crystal phase transitions for the Kihara fluid and temperatures investigated in this work. The tabulated densities are those of the consecutive states in the MC simulation belonging to different phases, whereas the coexistence pressure is obtained from the mean value of the pressures of the two boundary states. Similar results for the HSC 关3兴, the SWSC, and the SRS fluids 关5兴are also included for comparison. I–Sm-AI-NN–Sm-A P* ␳ I * ␳ SmA *P* ␳ I * ␳ N *P* ␳ N * ␳ SmA * T*⫽1.5 1.675 0.1117 0.1326 T*⫽2 1.675 0.1126 0.1288 T*⫽3 1.425 0.1069 0.1101 1.725 0.1193 0.1294 T*⫽5 1.325 0.1067 0.1101 1.775 0.1260 0.1348 SRS T*⫽5 1.175 0.0990 0.1027 1.625 0.1191 0.1250 SWSC T*⫽5 1.325 0.0968 0.0993 1.675 0.1099 0.1155 HSC 1.190 0.0914 0.0932 1.540 0.1061 0.1095 FIG. 3. Pair potential energy along the four isotherms studied in this work for the Kihara fluid. The pair energies for the T*⫽5 isotherm of the SRS fluid from Ref. 关5兴are also included for comparison. LIQUID CRYSTAL BEHAVIOR OF THE KIHARA FLUID PHYSICAL REVIEW E 68, 011704 共2003兲 011704-3 A further trend readily noted in Fig. 3 for the Kihara fluid is the global shift of the energy per particle toward more negative values as the temperature is reduced; that is, the attractive contribution to the pair interaction becomes progressively dominant when cooling down the system. The inspection of the radial distribution functions 共not shown兲indicates that, in fact, the average confinement of neighboring particles within the attractive well is more effective at low temperature, whereas at higher temperatures the pairs show a greater degree of penetration into the short-range repulsive region of the potential. IV. SUMMARY AND CONCLUSIONS It can be concluded that the Kihara potential constitutes a valid model to study the phenomenology of mesogenic fluids. The potential interest of the present results is twofold: they aid to the development of molecular modeling strategies for mesogens with specific phase diagrams and provide a closer assertion of the role played by the configurational entropy and the dispersive interactions in the behavior of molecules in dense environments, typical of liquid crystals but also of the macromolecular crowding conditions of many biological systems. As already mentioned, the advantage of the Kihara model is that it incorporates a more realistic description of the molecular shape and of the dispersive pair interactions in comparison to other related models. As shown above, the pair interaction in the Kihara fluid brings the liquid crystal phases to significantly higher densities in comparison to the HSC, the SWSC, and the SRS fluids, due to the combined effect of energetic and excluded volume entropic constrains 关5兴.In particular, the attractive well has a beneficial influence upon the stabilization of the smectic phase with respect to the nematic one and contributes to the destabilization of this latter phase at low temperature. The present results have served to estimate, for the L* ⫽5 fluid, the isotropic-nematic–smectic-Atriple point, a thermodynamical parameter relatively scarce in the literature for any of the model fluids commonly employed. A physical interpretation of the supercritical location of the triple point found here for the Kihara fluid can be grasped from the behavior of peptide chains, which are often modelled by rigid long rods. Recent simulations for polyalanine show that the critical temperature of this peptide is well above its decomposition temperature 关16兴. Since this is a common situation for proteins 共with rare exceptions, such as collagen兲,itmay seem that a supercritical triple point, although realistic for ideal stable molecules, would take place at temperatures unattainable for pure systems. However, elongated peptides and, in particular, polyalanine undergo first-order transitions in solution 关16兴. Our simulations are then relevant within this context, where the effect of solute concentration is analogous to that of an inverse temperature. The more complete picture of the interplay of entropy and energy in the Kihara fluid would be desirable from an extension of the present work to fluids of varying elongation. In addition, one drawback of the Kihara model is that it assigns the same interaction energy to all pair orientations as long as the minimum distance between the molecules remains constant. In real systems, in contrast, the dispersive interactions are typically greater for aligned configurations than for misaligned configurations of the molecular pairs. This kind of behavior is qualitatively reproduced by models composed of rigid chains with Lennard-Jones interaction sites. Although multiple-site models of this latter type have been employed to study liquid crystal phases 关2,17兴, their use is limited by the computational cost associated to the large number of sites required in order to mimic realistic mesogenic molecules. Our group is at present working on an alternative approach based on a refinement of the Kihara model that incorporates an orientational dependence of the attractive well with the aim to explore the qualitative effects that this feature may introduce in the phase diagram of the fluid. ACKNOWLEDGMENT We acknowledge support from the Spanish Direccio ´n General de Investigacio ´n Cientı ´fica y Te ´cnica 共Grant No. BQU2001-3615-C02兲, Instituto de Salud Carlos III 共Grant No. 01/1664兲, and Plan Andaluz de Investigacio ´n共groups FQM-205, FQM-319兲. 关1兴M.P. Allen, G.T. Evans, D. Frenkel, and B.M. Mulder, Adv. Chem. Phys. 86,1共1993兲. 关2兴L.F. Rull, Physica A 220,113共1995兲, and citations therein. 关3兴S.C. McGrother, D.C. Williamson, and G. Jackson, J. Chem. Phys. 104, 6755 共1996兲. 关4兴D.J. Earl, J. Ilnytskyi, and M. Wilson, Mol. Phys. 99, 1719 共2001兲. 关5兴A. Cuetos, B. Martı ´nez-Haya, L.F. Rull, and S. Lago, J. Chem. Phys. 117, 2934 共2002兲;117, 11405 共2002兲. 关6兴S. Lago, B. Garzo ´n, S. Calero, and C. Vega, J. Phys. Chem. 101, 6763 共1997兲. 关7兴M.S. Al-Barwani and M.P. Allen, Phys. Rev. E 62, 6706 共2000兲. 关8兴F. del Rı ´o, E. A ´valos, R. Espı ´ndola, L.F. Rull, G. Jackson, and S. Lago, Mol. Phys. 100, 2531 共2002兲. 关9兴M.A. Bates and G.R. Luckhurst, J. Chem. Phys. 110, 7087 共1999兲. 关10兴T. Kihara, Adv. Chem. Phys. 5, 147 共1963兲. 关11兴M. Caplow and L. Fee, Biochemistry 42, 2122 共2003兲. 关12兴J.T. Brown, M.P. Allen, E. Martı ´n del Rı ´o, and E. de Miguel, Phys. Rev. E 57, 6685 共1998兲. 关13兴E. de Miguel, E. Martı ´n del Rı ´o, J.T. Brown, and M.P. Allen, J. Chem. Phys. 105, 4234 共1996兲. 关14兴B. Martı ´nez-Haya, L.F. Rull, A. Cuetos, and S. Lago, Mol. Phys. 99, 509 共2001兲. 关15兴H.C. Andersen, J.D. Weeks, and D. Chandler, Phys. Rev. A 4, 1597 共1971兲. 关16兴共a兲N.A. Alves and U.H.E. Hansmann, Phys. Rev. Lett. 84, 1836 共2000兲;共b兲Y. Peng, U.H.E. Hansmann, and N.A. Alves, J. Chem. Phys. 118, 2374 共2003兲. 关17兴J. Ilnytskyi and M.R. Wilson, Comput. Phys. Commun. 134, 23 共2001兲. CUETOS et al. PHYSICAL REVIEW E 68, 011704 共2003兲 011704-4