A novel orientation-dependent potential model for prolate mesogens

J. Chem. Phys. 122, 024908 (2005); https://doi.org/10.1063/1.1830429 122, 024908 © 2005 American Institute of Physics. A novel orientation-dependent potential model for prolate mesogens Cite as: J. Chem. Phys. 122, 024908 (2005); https://doi.org/10.1063/1.1830429 Submitted: 01 December 2003 • Accepted: 19 October 2004 • Published Online: 22 December 2004 B. Martınez-Haya, A. Cuetos, S. Lago, et al. ARTICLES YOU MAY BE INTERESTED IN Columnar phases of discotic spherocylinders The Journal of Chemical Physics 129, 214706 (2008); https://doi.org/10.1063/1.3028539 Columnar phases of discotics with orientation-dependent interactions The Journal of Chemical Physics 131, 074901 (2009); https://doi.org/10.1063/1.3207284 A re-examination of the phase diagram of hard spherocylinders The Journal of Chemical Physics 104, 6755 (1996); https://doi.org/10.1063/1.471343 A novel orientation-dependent potential model for prolate mesogens B. Martı ´nez-Haya,a) A. Cuetos, 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 Postal 1065, 41080 Seville, Spain 共Received 1 December 2003; accepted 19 October 2004; published online 22 December 2004兲 An intermolecular potential is introduced for the study of molecular mesogenic fluids. The model combines distinct features of the well-known Gay-Berne and Kihara potentials by incorporating dispersive interactions dependent on the relative pair orientation to a spherocylinder molecular core. Results of a Monte Carlo simulation study focused on the liquid crystal phases exhibited by the model fluid are presented. For the chosen potential parameters, molecular aspect ratio L*⫽5 and temperatures T*⫽2, 3, and 5, isotropic, nematic, smectic-A, and hexatic phases are found. The location of the phase boundaries as well as the equation of state of the fluid and further thermodynamical and structural parameters are discussed and contrasted to the Kihara fluid. In comparison to this latter fluid, the model induces the formation of ordered liquid crystalline phases at lower packing fractions and it favors, in particular, the appearance of layered hexatic ordering as a consequence of the greater attractive interaction assigned to the parallel side-to-side molecular pair configurations. The results contribute to the evaluation of the role of specific interaction energies in the mesogenic behavior of prolate molecular liquids in dense environments. © 2005 American Institute of Physics. 关DOI: 10.1063/1.1830429兴 I. INTRODUCTION The use of simple models to represent the overall pair interactions in molecular fluids has proven to be a successful strategy to study the behavior of mesogens, as well as to predict and characterize their liquid crystal phases. The aim of such models is to capture the essential aspects of the physics underlaying the mesogenic behavior of the real systems, such as excluded volume effects and dispersive interactions, while keeping reasonable analytical and computational efficiency for theoretical and simulation studies. Fluids of elongated or rodlike molecules are an important class of mesogens with relevant technological and biological applications and, therefore, different models have been introduced in order to explore their properties. For instance, a family of rigid molecular models of ellipsoidal symmetry has been proposed, among which the Gay-Berne 共GB兲fluid ranks as one of the most extensively studied.1 This model extended the pioneering studies of Frenkel and co-workers on the hard ellipsoid fluid,2and was specifically introduced as an improvement of the Gaussian overlap model.3The main feature of the GB model is a fourparameter functionality that controls the aspect ratio of the ellipsoidal core and the anisotropy of the attractive interactions. The phase diagram of the GB fluid has been characterized for a broad range of aspect ratios and interaction parameters and, in particular, nematic and layered smectic and hexatic liquid crystal phases have been reported.4–6 In spite of the success of the Gay-Berne model, more detailed interaction approaches, such as site-site LennardJones chain models, indicate that the actual core of prolate molecules is significantly better reproduced by a spherocylinder core 共i.e., a cylinder of height/diameter aspect ratio L*⫽L/ ␴ , capped at both ends with a hemisphere of the same diameter ␴ 兲.5,7 In fact, a number of fluid models of this latter symmetry have been introduced in the past decades, especially after efficient algorithms were developed to compute the minimum distance between the central rods of such molecules.8Examples of this family of models are the hard spherocylinder 共HSC兲fluid and its square-well 共SWSC兲or soft repulsive 共SRS兲variants,9–14 and the Kihara fluid.15 The Kihara model was introduced as a generalization of the Lennard-Jones fluid for anisotropic molecules, and has been employed in numerous investigations of thermodynamic, structural, and transport properties of fluids of linear molecules.16 Perhaps surprising, it has not been until recently that the ability of the Kihara fluid to form liquid crystals phases has been systematically investigated.17 One drawback that the Kihara model shares with the HSC, SWSC, and SRS models is that it assigns the same interaction energy to all pair orientations, as long as the minimum distance between the molecules remains constant. This is in contrast with the interactions of real systems where the dispersive forces are orientation dependent and, unless specific interactions come into play, in thermotropic fluids they usually tend to be greater for aligned molecular configurations 共e.g., side-to-side parallel pairs兲than for the misaligned ones 共e.g., head-to-tail or T-shaped pairs兲.18,19 This behavior of the interaction energy is qualitatively reproduced by the rigid chain models with multiple interaction sites. However, although multiple-site models of this latter type have been employed to study liquid crystal phases,7,20 their a兲Electronic mail: [email protected] THE JOURNAL OF CHEMICAL PHYSICS 122, 024908 共2005兲 122, 024908-10021-9606/2005/122(2)/024908/8/$22.50 © 2005 American Institute of Physics use is limited by the computational cost associated to the large number of sites required in order to mimic realistic mesogenic molecules. The main idea behind the present work is to correct for the deficiency of the Kihara model commented above by incorporating one the most distinct features of the Gay-Berne fluid, namely, its parametric modulation of the pair orientation dependence of dispersive forces. The model is presented in Sec. II and the remaining of the paper is then devoted to explore, by means of Monte Carlo simulations, the qualitative effects that this feature introduces in the liquid crystal phase diagram of the fluid. II. INTERACTION MODEL In this work, we introduce an interaction model that incorporates the pair orientation dependence of the dispersive interactions of the Gay-Berne potential to the spherocylinder molecular core of the Kihara fluid. In order to achieve this in a straightforward and easily recognizable way, we have built a ‘‘hybrid’’ interaction energy functional by multiplying the Kihara potential by the same orientational prefactor of the Gay-Berne potential. These latter factor depends explicitly on the three-vector correlations between the directors of the given pair of particles (uˆi,uˆj) and a unit vector in the direction of the center-of-mass intermolecular distance vector (rˆij), in contrast to the Kihara interaction which depends on the relative orientation of the molecular pairs only through its influence on the minimum distance between the molecular cores dm(rij,uˆi,uˆj). We will refer to this model as the GayBerne-Kihara potential or, in short, GB-K potential. The model is thus defined by the following expressions: UGB-K共rij,uˆi,uˆj兲⫽ ⑀ GB共rˆij,uˆi,uˆj兲UK共dm兲,共1兲 UK共dm兲⫽4 ⑀ 关共 ␴ /dm兲12⫺共 ␴ /dm兲6兴,共2兲 ⑀ GB共rˆij,uˆi,uˆj兲⫽ ⑀ GO ␯ 共uˆi,uˆj兲 ⑀ ⬘ ␮ 共rˆij,uˆi,uˆj兲,共3兲 ⑀ GO共uˆi,uˆj兲⫽关1⫺ ␹ 2共uˆi•uˆj兲2兴⫺1/2,共4兲 ⑀ ⬘共rˆij,uˆi,uˆj兲⫽1⫺ ␹ ⬘ 2 冋 共rˆij•uˆi⫹rˆij•uˆj兲2 1⫹ ␹ ⬘共uˆi•uˆj兲 ⫹共rˆij•uˆi⫺rˆij•uˆj兲2 1⫺ ␹ ⬘共uˆi•uˆj兲 册 .共5兲 The prefactor ⑀ GB is characterized by the usual GayBerne four-parameter set 共 ␬ , ␬ ⬘, ␮ , and ␯ 兲, with the notation ␹ ⫽( ␬ 2⫺1)/( ␬ 2⫹1) and ␹ ⬘⫽( ␬ ⬘1/ ␮ ⫺1)/( ␬ ⬘1/ ␮ ⫹1). Simple geometrical arguments show that, in our model, the parameter associated to the ellipsoidal aspect ratio ␬ is related to the spherocylinder aspect ratio through ␬ ⫽L*⫹1 共note that L*is the height of the central cylinder, whereas the total length of the spherocylinder, including the end caps, is L*⫹1, always in units of the diameter ␴ 兲. On the other hand, the anisotropy of the dispersive interaction is controlled by ␬ ⬘; for instance, the attractive energy well of a parallel pair of molecules is ␬ ⬘times deeper for a side-toside configuration than for a head-to-tail one. For the present study, we have taken the set of parameters L*⫽5, ␬ ⫽6, ␬ ⬘⫽5, ␮ ⫽2, and ␯ ⫽1, so that the model here employed will be henceforth denoted GB-K共6,5,2,1兲according to the notation GB-K( ␬ , ␬ ⬘, ␮ , ␯ ) which is similar to the one used previously for the Gay-Berne fluid.14 The spherocylinder aspect ratio L*⫽5共and, hence, ␬ ⫽6) was chosen in order to compare the present results to our recent study of the Kihara fluid,17 whereas ␬ ⬘⫽5, ␮ ⫽2, ␯ ⫽1, are the values originally suggested by Gay and Berne1also employed in several previous studies of the liquid phase diagram of fluids of shorter aspect ratios ␬ ⭐4.4–6 It must be noted that the choice of ␬ ⬘, ␮ , ␯ is actually not a trivial task, especially when comparing with real systems, and other sets of values have been proposed for specific thermotropes.14,21 Figure 1 illustrates some of the main features of the GB-K potential model introduced in the previous paragraph. The top panel 关contour plot 共A兲兴 shows the equipotential contours for the interaction of two perfectly parallel GBK共6,5,2,1兲particles 共i.e., uˆi⫽uˆj) whose relative position is defined by the intermolecular distance vector in polar coordinates (rij, ␪ ) in the plane of the figure. As can be readily seen in the figure, the attractive well for the parallel side-toside configuration ( ␪ ⫽0°) is significantly enhanced with respect to the head-to-tail one ( ␪ ⫽90°). The same representaFIG. 1. Top panel: Equipotential energy surfaces for the interaction of two parallel particles interacting through 共a兲the GB-K共6,5,2,1兲potential introduced in this work 共with parameters L*⫽5, ␬ ⫽6, ␬ ⬘⫽5, ␮ ⫽2, and ␯ ⫽1, see text for details兲,共b兲the Kihara potential 共with L*⫽5), and 共c兲the Gay-Berne potential GB共6,5,2,1兲. The position of the pair of particles is described by the polar coordinates (rij , ␪ ), with ␪ ⫽0° for the side-to-side configuration and ␪ ⫽90° for the head-to-tail one. Bottom panel: GBK共6,5,2,1兲pair interaction energy as a function of the minimum distance between the molecular cores dmfor pairs of parallel molecules in side-toside, crossed, head-to-tail, and T-shaped configurations. 024908-2 Martı ´nez-Haya et al. J. Chem. Phys. 122, 024908 (2005) tion for the Kihara fluid 关plot 共B兲兴 obviously yields a uniform well around the molecular core. On the other hand, the similar diagram for the Gay-Berne potential with the same parameters, GB共6,5,2,1兲关plot 共C兲兴, illustrates its characteristic anisotropic well and ellipsoidal core. A closer quantitative representation of the anisotropy of the dispersive interactions in the GB-K共6,5,2,1兲model is provided in the bottom panel of Fig. 1, where the potential energy as a function of the minimum distance dmis represented for pairs of molecules in parallel, crossed, head-totail, and T-shaped configurations. The fivefold deeper well for the parallel configuration with respect to the head-to-tail, as corresponds to ␬ ⬘⫽5, is apparent in the figure. It must be noticed that our formulation of the GB-K model assigns a well depth of unity 共in reduced units兲to the crossed configuration, which must be kept in mind when comparing, at a given temperature, the phase diagram of this fluid to that of the Kihara fluid 共for which a well depth of unity applies to all molecular orientations兲. III. SIMULATION DETAILS We have carried out isothermal-isobaric (N-P-T) ensemble Monte Carlo 共MC兲simulations to study the liquid crystal phase diagram at different temperatures of the GBK共6,5,2,1兲fluid 共i.e., with parameters ␬ ⫽6, and hence L* ⫽5, ␬ ⬘⫽5, ␮ ⫽2, and ␯ ⫽1). For computing efficiency, in the present simulations the interaction was truncated at a distance of dC⫽3 ␴ , which corresponds to center-of-mass distances ranging from rij⫽3 ␴ to rij⫽8 ␴ , depending on the relative pair orientation, and shifted so that the potential energy UGB-K vanishes at the truncation boundary. Hence, the potential actually employed is given by UGB-K共rij,uˆi,uˆj兲 ⫽ 再 ⑀ GB共rˆij,uˆi,uˆj兲关UK共dm兲⫺UK共dC兲兴,dm⭐3 ␴ 0, dm⬎3 ␴ .共6兲 The simulations were run for a system of Np⫽1080 molecules at three reduced temperatures, T*⫽kT/ ⑀ ⫽2, 3, and 5, where kdenotes the Boltzmann constant and ⑀ the well depth for the crossed configuration 共see Fig. 1兲. As discussed throughout the following sections of the paper, at these temperatures the GB-K共6,5,2,1兲fluid presents stable isotropic (I), nematic (N), smectic-A(Sm-A), and hexatic 共Hex兲 phases. The calculation was started at each temperature by compressing the fluid to a state of high density well inside the hexatic region of the phase diagram (P*⫽P ␴ 3/kT ⫽2.4, 2.6, and 3.1, respectively, for T*⫽2, 3, and 5兲. Such state was properly equilibrated over some 106MC cycles before a systematic isothermal expansion of the fluid was performed down to the isotropic phase. Such approximate procedure for the estimation of liquid crystal phase boundaries is reliable for our purposes and has been extensively employed in the past,9,11,22 in particular, for similar fluids of prolate molecules. It must be stressed, however, that a proper study of the phase diagram should include the accurate calculation of the free energy in each of the phases, which is beyond the scope of the present study. Each state was typically equilibrated over 106MC cycles and ensemble averages of the thermodynamic and structural properties of the system were computed over 3 ⫻105cycles. The last molecular arrangement was then 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 a trial change of the box volume. The usual periodic boundary conditions are employed and the acceptance ratios were kept within 30%–40% for the tilt and displacement of the particles, and within 20%–30% for the box volume change. Box volume changes were attempted by randomly changing the length of each side of the box independently, with the restriction that none of them could become shorter than twice the range of the interaction potential, of 8 ␴ in the present case. Interestingly, we found a tendency of the box to spontaneously adopt an anisotropic geometry with side lengths fulfilling Lx⬇Ly⬍Lz共the subscript x denotes the shortest side兲, which is consistent with the geometry chosen for the simulations with fixed box geometry (Lx⬇Ly ⬇Lz/2) of our earlier studies of the SWSC, SRS, and Kihara fluids.11,17 In spite of this latter finding, in order to perform a proper comparison between the GB-K共6,5,2,1兲and the Kihara fluids, we have also recomputed full simulations for the four isotherms T*⫽kT/ ⑀ ⫽1.5, 2, 3, and 5 of the Kihara fluid studied in a previous work17 employing the same methodology described above. The result of these new simulations, which essentially reproduce the data of the earlier work,17 are presented below and include an extended region at low and high density with respect to our previous work in order to overlap with the present study of the GB-K model. The liquid crystalline transitions observed in the expansion of the fluid are characterized by discontinuities in the density and the nematic order parameter, as well as by sudden qualitative changes observed in the different correlation functions g 储 (r 储 ), g⬜(r⬜), g⬜ (n)(r⬜) defined in previous works.6,9,11,12,23 For instance, the function g 储 (r 储 ) represents the projection of the pair distribution function along the nematic director of the fluid and develops a characteristic oscillatory structure when layered phases are formed. On the other hand, g⬜ (0)(r⬜) accounts for the correlation between particles within the same layer and presents long-range structure when solidlike order is present. An illustrative insight into these correlation functions is provided below within the framework of the fluidlike smectic-Aor hexagonally packed hexatic phases found in this work for the GB-K共6,5,2,1兲 fluid. The formation of these latter hexatic phases was inspected, in addition, through the calculation of the bond hexagonal order parameter:22,24 H6⫽ 冓 冏 1 Np兺 j 1 nj兺 (kl)exp共6i ␪ kl兲 冏 冔 ,共7兲 where njis the number of pairs of neighbors within the first in-layer coordination shell of particle j, which is defined as a cylindrical volume of radius 1.6 ␴ and height 1.0 ␴ centered at the center of mass of the molecule. The sum over (kl) applies to all possible pairs of such first coordination shell particles and ␪ kl is the angle between the projection of the in024908-3 Novel potential model for prolate mesogens J. Chem. Phys. 122, 024908 (2005) termolecular vectors rjk and rjl onto the plane perpendicular to the director of molecule j. Hence, H6is so defined as to approach unity for perfect hexatic order within the smectic layers 共with nj⫽6 in-layer nearest neighbors around each particle forming a hexagon兲and tend to zero for disordered fluidlike layers. We employ the generic denomination hexatic throughout the paper, since the limited number of particles of our simulations does not allow to fully characterize this phase as a fluid smectic-Bphase or a solid crystal Bphase, as will be discussed below. IV. RESULTS AND DISCUSSION Figure 2 共top panel兲represents the liquid phase diagram of the GB-K共6,5,2,1兲fluid for the three isotherms T*⫽2, 3, and 5. The pressures and densities of the boundary states at each of the phase transitions observed are listed in Table I. A more detailed information about the isotherms is provided in Figs. 3, 4, and 5 which depict, respectively, the equations of state of the fluid, the average energy per particle, and the nematic and hexatic order parameters for each of the states sampled in the present simulations. The phase diagram and equations of state of the Kihara fluid, recalculated here employing the variable box geometry simulation methodology FIG. 2. Summary of the densities of the boundary states at the isotropicnematic (I-N), isotropic–smectic-A(I-Sm-A), nematic–smectic-A (N–Sm-A), or the smectic-A–hexatic (Sm-A–Hex) transitions obtained in the MC–N-P-Tsimulations of the GB-K共6,5,2,1兲fluid of the present work. The bottom panel shows the phase diagram for the Kihara fluid with the same molecular aspect ratio. These latter data are revised with respect to those reported in a previous work 共Ref. 17兲. See also Tables I and II. TABLE I. Reduced pressures and densities of the boundary states of the liquid crystal phase transitions of the GB-K共6,5,2,1兲fluid considered in this work at temperatures T*⫽2, 3, and 5. Note the reduction of pressure with temperature (P*⫽P ␴ 3/kT) when comparing the data. I–Sm-ASm-A–Hex pI * ␳ I *pSmA * ␳ SmA *pSm-A * ␳ Sm-A *pHex * ␳ Hex * T*⫽2.0 0.70 0.086共1兲0.75 0.104共1兲1.80 0.130共1兲1.90 0.137共1兲 I–NN–Sm-ASm-A–Hex pI * ␳ I *pN * ␳ N *pN * ␳ N *pSm-A * ␳ Sm-A *pSm-A * ␳ Sm-A *pHex * ␳ Hex * T*⫽3.0 1.10 0.100共1兲1.15 0.104共1兲1.40 0.112共1兲1.45 0.123共1兲2.25 0.140共1兲2.30 0.146共1兲 T*⫽5.0 1.30 0.107共1兲1.35 0.110共1兲2.35 0.137共1兲2.40 0.141共1兲2.95 0.151共1兲3.00 0.158共1兲 FIG. 3. Equations of state for the isotherms T*⫽2, 3, and 5 of the GBK共6,5,2,1兲and the Kihara fluids with spherocylinder aspect ratio L*⫽5. The filled circles correspond to the GB-K fluid whereas the open squares are for the Kihara fluid. The corresponding liquid crystal phases (I,N, Sm-A,or Hex兲are indicated next to the data. In the regions where the isotherms of the two fluids overlap but the phases are different, the phase corresponding to the Kihara fluid is indicated in parentheses. The statistical error bars associated to the density values are of the size of the symbols or smaller. Note the reduction of pressure with temperature (P*⫽P ␴ 3/kT) when comparing the data. 024908-4 Martı ´nez-Haya et al. J. Chem. Phys. 122, 024908 (2005) described in the preceding section, are also included in the same figures for direct comparison with the GB-K data. It was found that essentially the same results are obtained for the Kihara fluid with respect to our earlier work,17 except for a⫾0.05–0.10 difference in the reduced pressure of the system in the smectic states of same density. This latter effect is likely to be related to anisotropic stress effects25 induced by the fixed box geometry approach of our previous work that should not be affecting the present results. Other than that, the liquid crystal transitions relevant to the present work occur in the expansion of the Kihara fluid at the same packing fractions, within statistical uncertainty, as found in the fixed box geometry simulations of our earlier work.17 We focus now on the results for the GB-K model. At the three temperatures investigated, the GB-K共6,5,2,1兲fluid presents a stable hexatic phase at sufficiently high densities in which the fluid is internally ordered in smectic layers of hexagonally packed molecules. We leave the characterization of this liquid crystalline phase to the last part of this Section. When decreasing the system pressure, and hence expanding the fluid, the hexatic phase eventually melts to a smectic-A phase in which the layers become fluidlike and the twodimensional positional order within them is lost. A further expansion of the system leads to a transition to either a nematic phase 共for the isotherms T*⫽3 and 5兲or directly to an isotropic phase 共for T*⫽2). This latter phase transition shows a stronger first-order character than the rest of transitions observed, as the fluid rearranges from a relatively dense layered structure to a fully disordered isotropic phase with a much smaller packing fraction. One noticeable aspect of the phase diagram of the GBK共6,5,2,1兲fluid is the substantial reduction in the range of stability of the nematic phase when cooling down the system from T*⫽5toT*⫽3, as the transition from smectic-Ato nematic shifts to lower pressures and densities. In fact, the nematic phase disappears altogether at T*⫽2, as already pointed out above. Hence, the result of the MC simulations suggests the presence of an isotropic–nematic–smectic-A triple point at a temperature of roughly T*⬇2.5 and a density close to ␳ *⬇0.1 共see Fig. 2兲. The comparison of the phase diagram of the GBK共6,5,2,1兲fluid to that of the Kihara fluid in the same range of temperatures should provide relevant information about the relative importance of energy and entropy as ‘‘driving forces’’ of the liquid crystal transitions. As can be seen in Fig. 2, the overall features of the phase diagram of both fluids present many similarities, although it is interesting to find that the I–N–Sm-Atriple point of the Kihara fluid is located at roughly similar temperature but at a higher density FIG. 4. Average pair potential energy (U*⫽U/ ⑀ ) along the isotherms studied in this work for the GB-K共6,5,2,1兲and the Kihara fluids. The filled circles correspond to the GB-K fluid and the open squares to the Kihara fluid. The vertical error bars correspond to one standard deviation of the computed energies. FIG. 5. Top panel: Nematic order parameter S2and bond hexagonal order parameter H6as a function of number density for the GB-K共6,5,2,1兲fluid at reduced temperature T*⫽3. The vertical error bars correspond to one standard deviation of the computed order parameters. Bottom panel: In-layer distribution function in the direction perpendicular to the director g⬜ (0)(r⬜ *) in the smectic-A共short-dash curve兲and hexatic 共long-dash and solid curves兲 phases of the same fluid. The variable r⬜ *⫽r⬜/ ␴ represents the projection of the pair intermolecular distance vector onto the plane perpendicular to the director. Hence, this distribution shows the spatial pair correlation between the particles in the same smectic layer. The reduced pressure of each state is indicated in the legend, whereas the temperature is T*⫽3 in all cases. 024908-5 Novel potential model for prolate mesogens J. Chem. Phys. 122, 024908 (2005) than in the GB-K共6,5,2,1兲fluid. There are further significant differences between both phase diagrams which have to be attributed to the anisotropy of the dispersive interactions in the GB-K potential around the molecular core in comparison to the isotropic interactions of the Kihara fluid. A first relevant aspect is that the smectic phase of the Kihara fluid does not present any trace of in-layer hexagonal packing order within the range of densities scoped in the present study. In order to corroborate this observation, the density range covered in our earlier work for the Kihara fluid17 was extended to overlap with the present simulations. In fact, the hexatic states found for the GB-K共6,5,2,1兲fluid melted systematically to smectic-Astates when the GB-K potential was exchanged by the Kihara potential. Another particularity of the GB-K共6,5,2,1兲fluid readily observed when comparing the two panels of Fig. 2 is the systematic decrease of the reduced pressures and densities of the boundary states separating the phase transitions that takes place as the temperature of the fluid cools down from high (T*⫽5) to low (T*⫽2) temperature. This trend, which is especially pronounced for the N–Sm-Atransition, as commented above, is as well apparent for the Sm-A–Hex and I-Ntransitions of the GB-K fluid. In contrast to this behavior, the phase boundaries of the Kihara fluid 共Table II兲are considerably less sensitive to temperature 共with the exception of the N–Sm-Atransition, due to the instability of the nematic phase at low temperature兲. A first overall conclusion that can be drawn from these findings is that the orientation dependence introduced in the GB–K potential and, in particular, the enhanced interaction for aligned pair configurations 共see Fig. 1兲, favors the formation of the liquid crystal phases explored in the present work. In fact, the isotropic phase becomes significantly destabilized in comparison to the Kihara fluid, especially at low temperature. At T*⫽2, for instance, the density of the boundary isotropic state at the I–Sm-Atransition is ␳ *⫽0.086 in the GB-K共6,5,2,1兲fluid, in comparison to the much higher value of ␳ *⫽0.112 in the Kihara fluid 共see Table I兲. A somewhat unexpected finding that to some extent seems to contradict this latter conclusion is that the nematic phase of the GB-K共6,5,2,1兲fluid at the highest temperature investigated, T*⫽5, presents a greater range of stability than in the Kihara fluid. In other words, the N–Sm-Atransition is delayed toward higher densities within the nematic phase in the former fluid in comparison to the latter one. In addition, it is also observed that the change in density at the transition is largely reduced in the GB-K共6,5,2,1兲fluid. A possible explanation for these features is that the GB-K nematic states become actually more compressible than the Kihara ones as a consequence of the slower rise of the short-range repulsive forces for a large fraction of pair orientations. In Fig. 1 it can be appreciated that the repulsive part of the interaction at distances shorter than dm/ ␴ ⫽1共the zero of the potential for all pair orientations兲becomes weaker for molecular configurations with a prefactor ⑀ GB smaller than unity. At the same time, the core repulsion is enhanced in the GB-K共6,5,2,1兲 potential with respect to the Kihara potential for the side-toside parallel configurations typical of the smectic phase. As a consequence of this, the effective volume excluded by the molecules in the smectic layers is greater in the GB-K fluid, so that the entropic constrains work in favor of the more compressible nematic phase 共in which many pair orientations far from the side-to-side one are relevant兲, hence delaying the N–Sm-Atransition. Since the short pair distances at which the repulsive part of the potential dominates are more efficiently accessed by the molecules of the fluid at high temperature, it is not surprising that this effect is only observed in our simulations at T*⫽5, whereas at the lower temperatures, more sensitive to the attractive part of the interaction, the smectic-Aphase is stable down to significantly lower densities in the GB-K共6,5,2,1兲fluid than in the Kihara one.Figure 3 compares the equations of state (P* ⫽P ␴ 3/kT versus ␳ *⫽ ␳␴ 3) along the isotherms T*⫽2, 3, and 5 for the GB-K共6,5,2,1兲and the Kihara fluids. We first of all remark that in all cases the isotherms of both models consistently converge at low density. With increasing temperature the behavior of the fluids is expected to become progressively less sensitive to the details of the pair interaction potential 共especially at low density where the attractive forces are dominant兲and, in fact, at T*⫽5 the good overlap between the equations of state of the GB-K共6,5,2,1兲and the Kihara fluids extends to a substantial region of the nematic TABLE II. Reduced pressures and densities of the boundary states of the liquid crystal phase transitions of the Kihara fluid considered in this work at temperatures T*⫽1.5, 2.0, 3.0, and 5.0. Note the reduction of pressure with temperature (P*⫽P ␴ 3/kT) when comparing the data. The tabulated values differ slightly with respect to those given in a recent work 共Ref. 17兲due to the greater system size and the improved pressure equilibration methodology of the present simulations 共see text for details兲. I–Sm-A pI * ␳ I *pSMA * ␳ SMA * T*⫽1.5 1.75 0.114共1兲1.80 0.132共1兲 T*⫽2.0 1.65 0.112共1兲1.70 0.125共1兲 I–NN–Sm-A pI * ␳ I *pN * ␳ N *pN * ␳ N *pSM-A * ␳ SM-A * T*⫽3.0 1.40 0.108共1兲1.45 0.111共1兲1.70 0.119共1兲1.75 0.130共1兲 T*⫽5.0 1.30 0.107共1兲1.35 0.110共1兲1.75 0.128共1兲1.80 0.135共1兲 024908-6 Martı ´nez-Haya et al. J. Chem. Phys. 122, 024908 (2005) phase. At the sufficiently high density, however, significant differences between both systems eventually arise. At T* ⫽5, for smectic-Astates of similar density, the pressure is significantly greater for the GB-K fluid than for the Kihara fluid. Interestingly, this trend tends to reverse at the lower temperatures, so that at T*⫽2 and 3 the pressure for the Kihara states becomes comparable to that of the GB-K states with the same density. At these lower temperatures, however, the comparison between the equations of state of the GBK共6,5,2,1兲and the Kihara fluids at high density is not as straightforward as for T*⫽5, since the overlap between the liquid crystal phases is greatly reduced. For instance, the accidental coincidence of the equations of state of both fluids is noticeable at T*⫽2 in the high density range ( ␳ * ⬎0.137) where the GB-K共6,5,2,1兲fluid has entered the hexatic phase, whereas the Kihara fluid remains in the Sm-A phase. The reason for finding greater pressures in the GB-K fluid at high temperature is again related to the larger prefactor in the GB-K interaction potential which enhances the short-range repulsion of the side-to-side pair configurations 共at the same time as it increases the depth of the well at larger distances兲. The repulsive wall is felt more efficiently by the molecules at high temperatures with the net effect of increasing the pressure of the system with respect to the Kihara fluid. At lower temperatures, the greater average attractions reduce the pressure of the GB-K system, especially in the smectic-Aphase. From the comparison of the phase diagrams and equations of state of the GB-K共6,5,2,1兲and the Kihara fluids it becomes apparent that, even though steric excluded volume effects constitute the main constrain that drives the liquid crystal transitions of the prolate molecular fluids of the type considered here, the energetic contribution to the free energy has an important role that cannot be neglected and, hence, that an appropriate description of the dispersive forces is required in order to model the properties of real mesogens. Further support for these considerations is provided in Fig. 4, which reveals the qualitatively different evolution in the GB-K and Kihara fluids of the average internal energy per particle along the three isotherms considered in our work. In the Kihara fluid, the energy grows monotonously with density at high temperature (T*⫽3 and 5兲, whereas a weak decreasing trend is observed at lower temperatures (T* ⫽2) in the low density side of the isotropic branch and for the smectic-Astates.17 The fact that the internal energy grows with density at high temperature can be traced back to the progressive relevance of the repulsive part of the pair potential at the relatively high densities considered in our study. In the Kihara fluid, the attractive well only dominates at sufficiently low temperatures or, alternatively, in a diluted density regime, much closer to the ideal gas limit than here considered. Thus, it turns out that in the range of densities relevant for liquid crystal behavior, the isotropic Kihara attraction incorporates little novel qualitative features with respect to the similar behavior also found for the purely repulsive SRS fluid 共basically a truncated Kihara fluid without attractive well兲.11 On the contrary, the specific topology of the attractive interactions in the GB-K fluid leads to a qualitatively different role of energy in the stability of the liquid crystal phases. Especially noticeable is the drop in energy that takes place as the system enters the smectic-Aphase from either the isotropic (T*⫽2) or nematic (T*⫽3 and 5兲 phase. Such energy discontinuity is absent or much weaker in any of the spherocylinder model fluids proposed in the past, such as the Kihara17 or the SRS and SWSC fluids,11 and is a direct consequence of the bias of the GB-K potential favoring specific pair configurations. It follows that the stability of the smectic phases is supported not only by steric effects, but also by the beneficial contribution of the dispersive interactions to the free energy. Interestingly, a significant further drop in energy is also observed in the three isotherms at the Sm-A–hexatic transition. This reinforces the idea that the appearance of hexagonal order in the GB-K共6,5,2,1兲fluid at densities where the Kihara fluid maintains a stable smectic-Aphase is also closely related to the energetic effects induced by the anisotropic nature of the interaction potential. In order to characterize the hexatic phases observed in the GB-K共6,5,2,1兲fluid and clearly discern them from the smectic-Aphases, several diagnostics were applied during the simulations. The transition from the smectic-Ato the hexatic phase implies a sudden change of the hexagonal bond order parameter 关Eq. 共7兲兴, in our case from roughly H6⬇0.2 to H6⬇0.5, as can be seen in the upper panel of Fig. 5 where H6is represented along with the nematic order parameter for the isotherm T*⫽3. The formation of the hexatic phase can also be visualized through the long-range structure appearing abruptly in the in-layer distribution function g⬜ (0)(r⬜)共the pair distribution within the plane of the layer兲. Figure 5 illustrates the qualitative change that this function undergoes at the Sm-A–hexatic transition at T*⫽3, from a smooth liquidlike distribution with only short-range order to a much more structured two-dimensional hexatic crystal-like correlation function with well-defined nearest-neighbor positions over several coordination shells. A controversial topic present in several recent works on the Gay-Berne fluid has been whether the hexatic phases formed in that system are actually liquidlike smectic-B phases or solidlike crystal Bphases on a macroscopic basis.6 We advance that our present results do not serve to close this far from trivial question. At sufficiently high density a stable solid crystalline phase is known to form in spherocylinder core fluids.9In fact, it has been argued that the absence of a well-defined phase transition when compressing and/or cooling down the hexatic Gay-Berne fluid indicates that the phase must be crystalline from the beginning. Further evidence in favor of the crystalline nature of the hexatic phases of the type observed for the Gay-Berne system arise from the calculation of shear viscosities.23 However, the difference between the smectic-Band the crystal Bphases is actually a subtle one for the finite size systems employed in computer simulations. In addition, care must be taken when dealing with such high density states where the mobility of the particles is greatly reduced and metastable states or inefficient sampling of the phase space may affect the apparent structure of the fluid. One further aspect about the smectic-B phases is that the long-range order within the layers is truncated on a macroscopic scale by the presence of a large num024908-7 Novel potential model for prolate mesogens J. Chem. Phys. 122, 024908 (2005) ber of punctual defects, in contrast with the crystal Bphase where the long-range order extends over a macroscopic distance.6We close by noting to this respect that, within the range of densities and temperatures explored in our simulation, the hexatic order parameter remains at values below 0.6 共see Fig. 5兲, that is much smaller than the limiting value of unity, which is an indication of a substantial presence of defects within the layers of the fluid, in principle compatible with a smectic-Bphase. However, the limited size of our simulation box prevents us from making definite statements about the range of the in-layer correlations and we cannot therefore draw conclusions in favor of the smectic-Bor crystal Bcharacter of the observed hexatic phase. V. SUMMARY AND CONCLUSIONS A rigid model potential, referred to as Gay-Berne-Kihara or GB-K potential, has been introduced which is expected to be reliable for the study of molecular mesogenic fluids. The GB-K model features a spherocylinder molecular core dressed with dispersive interactions dependent on the relative pair orientation. The mathematical formulation of the model is compact and combines the functionalities of the wellknown Kihara and Gay-Berne potentials. The results of the Monte Carlo simulation study of the model fluid for a specific set of parameters at temperatures T*⫽2, 3, and 5 presented in this paper show that it is capable to reproduce the isotropic, nematic, smectic-A, and hexatic liquid crystal phases observed in real mesogens. At the three temperatures investigated, the GB-K共6,5,2,1兲fluid presents a stable hexatic phase at sufficiently high densities, characterized by layers of hexagonally packed molecules. When expanding the fluid, the hexatic phase eventually melts to a smectic-Aphase structured in fluidlike disordered layers. A further expansion of the system leads to a transition to either a nematic phase 共at T*⫽3 and 5兲or directly to an isotropic phase 共at T*⫽2). In this latter phase change the fluid undergoes a strongly first-order transition from a dense layered structure to a fully disordered isotropic phase with a substantial change in density. In comparison to the Kihara fluid, the GB-K共6,5,2,1兲 fluid is found to favor the formation of the ordered liquid crystalline phases and, specifically, the appearance of layered hexatic order, at lower packing fractions. This property can be interpreted as being a direct consequence of the greater dispersive interactions assigned in the GB-K fluid to specific pair orientations, such as the parallel side-to-side configuration. In fact, the specific topology of the GB-K potential with respect to the Kihara potential leads to a qualitatively different behavior of the internal energy and of its influence on the stability of the liquid crystal phases. In particular, it is found that the entrance of the fluid in the smectic phases (Sm-Aor hexatic兲is accompanied by a substantial stabilization of the internal energy of the fluid, in contrast to the behavior of other spherocylinder model fluids studied previously, such as the Kihara, or the SRS and SWSC fluids. Hence, it becomes apparent that the energetic contribution to the free energy has an important role in the mesogenic behavior of prolate molecular liquids in dense environments, and that an appropriate description and treatment of the dispersive forces is required in order to model accurately the properties of real mesogens, even at a qualitative level. The influence of the short-range interactions on the internal structure of the molecular fluids has as well been stressed recently in a study of systems composed of linear dipolar molecules.26 The main advantage of the GB-K potential is that it combines, within the intrinsic limitations of the rigid models, a qualitatively more adequate description of both, the pair interactions and the molecular shape of the typical mesogens, in comparison to previous models, whereas it keeps a comparable compactness and numerical efficiency in its formulation. The present study has focused on the presentation of the GB-K model and has stressed the qualitative effects introduced by the anisotropic dispersive interactions of the model in the equation of state and in the phase diagram of the fluid. Future work in our group will be devoted to compare the behavior of the 共spherocylinder兲GB-K fluid to that of the 共ellipsoidal兲Gay-Berne fluid, so that the relevance of the exact shape of the molecular core at supercritical and subcritical temperatures will be exposed. In addition, the implementation of the GB-K model to oblate 共disk-like兲mesogens will be explored. ACKNOWLEDGMENTS The authors acknowledge support from the Spanish Direccio ´n General de Investigacio ´n Cientı ´fica y Te ´cnica 共Groups No. BQU2001-3615-C02兲and Plan Andaluz de Investigacio ´n共Grant Nos. FQM-205 and FQM-319兲. 1J. G. Gay and B. J. Berne, J. Chem. Phys. 74, 3316 共1981兲. 2M. P. Allen, G. T. Evans, D. Frenkel, and B. M. Mulder, Adv. Chem. Phys. 86,1共1993兲, and references therein. 3B. J. Berne and P. Pechukas, J. Chem. Phys. 56,4213共1972兲. 4C. Zannoni, J. Mater. Chem. 11, 2637 共2001兲. 5L. F. Rull, Physica A 220,113共1995兲, and references therein. 6E. de Miguel and C. Vega, J. Chem. Phys. 117, 6313 共2002兲. 7G. V. Paolini, G. Ciccotti, and M. Ferrario, Mol. Phys. 80,297共1993兲. 8C. Vega and S. Lago, Comput. Chem. 18,55 9S. C. McGrother, D. C. Williamson, and G. Jackson, J. Chem. Phys. 104, 6755 共1996兲, and references therein. 10D. J. Earl, J. Ilnytskyi, and M. Wilson, Mol. Phys. 99, 1719 共2001兲. 11 A. Cuetos, B. Martı ´nez-Haya, L. F. Rull, and S. Lago, J. Chem. Phys. 117, 2934 共2002兲;117, 11405共E兲共2002兲. 12M. S. Al-Barwani and M. P. Allen, Phys. Rev. E 62, 6706 共2000兲. 13F. del Rı ´o, E. A ´valos, R. Espı ´ndola, L. F. Rull, G. Jackson, and S. Lago, Mol. Phys. 100, 2531 共2002兲. 14M. A. Bates and G. R. Luckhurst, J. Chem. Phys. 110,7087共1999兲. 15T. Kihara, Adv. Chem. Phys. 5, 147 共1963兲. 16S. Lago, B. Garzo ´n, S. Calero, and C. Vega, J. Phys. Chem. 101,6763 共1997兲. 17A. Cuetos, B. Martı ´nez-Haya, S. Lago, and L. F. Rull, Phys. Rev. E 68, 011704 共2003兲. 18R. A. Kromhaut and B. Linder, J. Phys. Chem. 99, 16909 共1995兲. 19D. P. Ojha and V. G. K. M. Pisipati, Z. Naturforsch., A: Phys. Sci. 57a, 645 共2002兲. 20J. Ilnytskyi and M. R. Wilson, Comput. Phys. Commun. 134,23共2001兲. 21G. R. Luckhurst and P. S. J. Simmonds, Mol. Phys. 80, 233 共1993兲. 22M. Houssa, L. F. Rull, and S. C. McGrother, J. Chem. Phys. 109, 9529 共1998兲. 23J. T. Brown, M. P. Allen, E. Martı ´n del Rı ´o, and E. de Miguel, Phys. Rev. E57, 6685 共1998兲. 24H. Zewdie, Phys. Rev. E 57, 1793 共1998兲. 25H. Domı ´nguez, E. Velasco, and J. Alejandre, Mol. Phys. 100,273共2002兲. 26S. Lago, S. Lo ´pez-Vidal, B. Garzo ´n, J. A. Mejı ´as, J. A. Anta, and S. Calero, Phys. Rev. E 68, 021201 共2003兲. 024908-8 Martı ´nez-Haya et al. J. Chem. Phys. 122, 024908 (2005)