Monte Carlo study of liquid crystal phases of hard and soft spherocylinders

Monte Carlo study of liquid crystal phases of hard and soft spherocylinders A. Cuetos, B. Martınez-Haya, L. F. Rull, and S. Lago Citation: The Journal of Chemical Physics 117, 2934 (2002); View online: https://doi.org/10.1063/1.1491872 View Table of Contents: http://aip.scitation.org/toc/jcp/117/6 Published by the American Institute of Physics Articles you may be interested in A re-examination of the phase diagram of hard spherocylinders The Journal of Chemical Physics 104, 6755 (1996); 10.1063/1.471343 Tracing the phase boundaries of hard spherocylinders The Journal of Chemical Physics 106, 666 (1997); 10.1063/1.473404 Orientational ordering and phase behaviour of binary mixtures of hard spheres and hard spherocylinders The Journal of Chemical Physics 143, 044906 (2015); 10.1063/1.4923291 A new anisotropic soft-core model for the simulation of liquid crystal mesophases The Journal of Chemical Physics 128, 044906 (2008); 10.1063/1.2825292 Columnar phases of discotic spherocylinders The Journal of Chemical Physics 129, 214706 (2008); 10.1063/1.3028539 A novel orientation-dependent potential model for prolate mesogens The Journal of Chemical Physics 122, 024908 (2004); 10.1063/1.1830429 Monte Carlo study of liquid crystal phases of hard and soft spherocylinders A. Cuetos and B. Martı ´nez-Hayaa) Departamento de Ciencias Ambientales, Universidad Pablo de Olavide, 41013 Sevilla, Spain L. F. Rull Departamento de Fı ´sica Ato ´mica, Molecular y Nuclear, Area de Fı ´sica Teo ´rica, Universidad de Sevilla, Apdo. 1065, 41080 Sevilla, Spain S. Lago Departamento de Ciencias Ambientales, Universidad Pablo de Olavide, 41013 Sevilla, Spain 共Received 2 January 2002; accepted 15 May 2002兲 We report on a Monte Carlo study of the liquid crystal phases of two model fluids of linear elongated molecules: 共a兲hard spherocylinders with an attractive square-well 共SWSC兲and 共b兲purely repulsive soft spherocylinders 共SRS兲, in both cases for a length-to-breadth ratio L*⫽5. Monte Carlo simulations in the isothermal–isobaric ensemble have been performed at a reduced temperature T*⫽5 probing thermodynamic states within the isotropic 共I兲, nematic 共N兲, and smectic A 共Sm A兲 regions exhibited by each of the models. In addition, the performance of an entropy criterion to allocate liquid crystalline phase boundaries, recently proposed for the isotropic–nematic transition of the hard spherocylinder 共HSC兲fluid, is successfully tested for the SWSC and the SRS fluids and furthermore extended to the study of the nematic–smectic transition. With respect to the more extensively studied HSC fluid, the introduction of the attractive square well in the SWSC model brings the I–N and N–Sm A transitions to higher pressures and densities. Moreover, the soft repulsive core of the SRS fluid induces a similar but quite more significant shift of both of these phase boundaries toward higher densities. This latter effect is apparently in contrast with very recent studies of the SRS fluid at lower temperatures, but this discrepancy can be traced back to the different effective size of the molecular repulsive core at different temperatures. © 2002 American Institute of Physics. 关DOI: 10.1063/1.1491872兴 I. INTRODUCTION The current understanding of the behavior of liquid crystalline mesogens relies largely on detailed studies of simple molecular fluid models. Over the past decades, a number of models of varying complexity have been developed with the aim of incorporating the essential features of the real systems, while minimizing the expense of their theoretical and numerical treatment. Such expense is mainly imposed by the intrinsic orientational dependence of the pair interaction and the long-range order characteristic of the liquid crystal phases. An important class of such liquid crystal models encompasses fluids of linear molecules interacting via a hard or soft repulsive core of either ellipsoidal or spherocylindrical symmetry, eventually dressed by attractive interactions of different nature.1–5 In particular, the hard spherocylinder fluid 共HSC兲constitutes one of the most extensively studied models and, in spite of its simplicity and purely repulsive interaction, it has been found to display differentiated isotropic, nematic, smectic, and solid phases.1,3 In the present work, we have concentrated on two of the simplest modifications of the HSC model which, perhaps surprisingly, have not deserved as much attention as other elemental liquid crystal models: a fluid of hard spherocylinders with an attractive square-well 共SWSC兲,6–8 and a fluid of soft repulsive spherocylinders 共SRS兲.5,9–11 We hence intend to explore the effect of the simple features introduced in these two models with respect to the HSC fluid, namely the presence of an attractive well or the softness of the short range repulsive interaction, on the liquid crystal phase diagram of the fluid. In addition, unlike the HSC fluid, these models allow for the study of systems at finite temperature in a natural way. One further fundamental aspect of this type of study will necessarily be related to the assessment, at an elementary level, of the roles played by two of the main extensive magnitudes, internal energy and entropy, in driving the liquid crystalline transitions observed in the fluids under study. We will focus, in particular, on the extension of earlier works that have explored the relationship between statistical entropy and ordering transitions in atomic and molecular fluids by resorting to the multiparticle correlation expansion of entropy established by Green and Nettleton12 and later generalized by Lazaridis and co-workers to nonspherical systems.13 In this context, many authors have studied the response to ordering in the fluid of the residual multiparticle entropy, defined as the total contributions to the excess entropy of all correlations involving more than two particles.14–17 Costa et al.17 were the first to employ this methodology to orientational transitions in liquid crystals and found a sudden growth of ⌬sat the isotropic–nematic transition in a HSC fluid. The authors of this latter work imposed a reduced dimensionality in the orientational depena兲Author to whom correspondence should be addressed. Electronic mail: [email protected] JOURNAL OF CHEMICAL PHYSICS VOLUME 117, NUMBER 6 8 AUGUST 2002 29340021-9606/2002/117(6)/2934/13/$19.00 © 2002 American Institute of Physics dence of the radial distribution function and proposed the crossover of ⌬sfrom negative to positive values as an empirical indication for the incipient transition to the nematic phase, a criterion that had previously been applied to the freezing transition in the hard sphere fluid,14 but that showed only moderate success in atomic fluids with attractive interactions.15 Thus, in the present study we have addressed the question of whether the entropy criterion, although contrasted for the HSC fluid, does keep its consistency for the liquid crystalline transitions in the SWSC and SRS fluids. The relevant details of the SWSC and SRS models and of the Monte Carlo simulation methodology employed in this work, including the characterization of the different liquid crystalline phases and the computation of the different structural and thermodynamic magnitudes, are described in Sec. II. The simulation results are presented and discussed thoroughly in Sec. III. Finally, a summary of the main conclusions is drawn in Sec. IV. II. SIMULATION METHOD A. Model fluids and simulation details We have applied the isothermal–isobaric 共NPT兲version of the Monte Carlo technique to evaluate thermodynamic and structural properties of the SWSC and SRS fluids. In the SWSC fluid, the hard core of the molecules comprises a cylinder of diameter ␴ and elongation L*⬅L/ ␴ ⫽5, with hemispherical caps on both ends of the same diameter ␴ . The attractive square well interaction is characterized by a depth ␧and a range ␭⫽1.5 ␴ , and has the same anisotropy of the hard-core. Thus, the pair interaction potential for two such spherocylinders is given by V共r,⍀兲⫽ 再 ⬁dm⭐ ␴ ⫺␧ ␴ ⬍dm⭐␭ 0dm⬎␭ 共1兲 as represented in Fig. 1. The distance between a pair of molecules, dm⫽dm(r,⍀), is defined as the minimum distance between the segments described by the axis of the central cylinder. The value of dmis a function of the relative orientation of the molecules, ⍀, determined by three independent angles, and of the distance between their centers-of-mass, r, as well as of the molecular parameters: diameter and elongation.In the second model, the SRS fluid,5,9–11 the same type of linear molecules interact through a soft and purely repulsive potential of spherocylindrical shape, built from a truncation of a shifted Kihara potential 共see Fig. 1兲which assures continuity of the potential and its first derivative, V共r,⍀兲⫽ 再 4␧关共 ␴ /dm兲12⫺共 ␴ /dm兲6⫹1/4兴dm⭐ 冑 62• ␴ 0dm⬎ 冑 62• ␴ . 共2兲 The MC-NPT simulations were run at reduced temperature T*⬅kT/␧⫽5共kdenotes the Boltzmann constant兲and started at low pressure (P*⬅P• ␴ 3/kT⫽0.1) deep in the isotropic region of the phase diagram, over a system of NP ⫽768 molecules initially arranged on a hexagonal closepacked 共hcp兲lattice. Following in part the procedure of McGrother et al.,3we built the hcp lattice with the 关111兴 direction along the zaxis and a distance between the close packed 关111兴planes scaled with (1⫹L*). The hcp unit cell contains four molecules and we chose a box with 12, 8, and 2 unit cells in the x,y, and zdirections, respectively. In this way, the dimensions of the simulation box along each axis fulfill the ratios Ly/Lx⫽1.15 and Lz/Lx⫽1.94. The longer dimension along the zdirection is meant to increase the number of possible smectic layers within the simulation box 共up to 4 layers in the close-packed configuration兲with respect to the cubic box for a given number of molecules. After equilibration and averaging, the last molecular arrangement is stored and used as initial configuration for the subsequent run with an increased system pressure. Thus, in this way the fluid is compressed through the isotropic–nematic and nematic–smectic transitions. In order to check for possible metastability and hysteric effects, an expansion run, all the way back to the isotropic region, is subsequently performed beginning with the smectic state of highest density in the calculation. As will be discussed below, the expansion run indeed resulted in a better equilibration of the calculations at high density. The 共first-order兲liquid crystalline transitions are located by the discontinuities in the density and the order parameter, aided by the analysis of the different correlation functions. The procedure is similar to that employed by McGrother et al. in their study of the HSC fluid3and very recently also by Earl et al. for the SRS fluid.5The nematic order parameter S⬅ 具 P2(ui"n) 典 is calculated as the ensemble average of the mean value of the second Legendre polynomial, its argument being the dot product between the unitary vector along the molecular axis of each particle, ui, and the 共also unitary兲 nematic director n. For the present study, together with the usual pair distribution function g(r), a relevant orientational distribution function is g2(r)⬅ 具 P2(ui"uj) 典 , defined as the average, for each intermolecular distance, of the second Legendre polynomial with argument the cosine of the relative angle between the axis of particles iand j,ui"uj⫽cos ␪ .In addition, the translational order of the molecules is further examined by means of the functions g 储 (r 储 ) and g⬜(r⬜), FIG. 1. Pair interaction potential energies as a function of the minimum distance, dm, between the molecules for the model fluids studied in this work: square-well spherocylinder fluid 关SWSC, see Eq. 共1兲兴 and soft repulsive spherocylinder 关SRS, see Eq. 共2兲兴. The value of dmdepends on the distance between the centers-of-mass of the molecules, r, and on their orientations as defined by the vectors u1and u2. 2935J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Liquid crystal phases of hard and soft spherocylinders which are the projections of the pair distribution function a distance r 储 parallel and a distance r⬜perpendicular to the nematic director, respectively. As shown by McGrother et al.,3smectic ordering induces well defined oscillations in g 储 (r 储 ), whereas the structure within the smectic planes is reflected in g⬜(r⬜). Furthermore, we checked for smectic B hexatic ordering within the smectic layers. We already advance, however, that only smectic A phases were found in the present study. The thermodynamic 共density, energy兲and structural properties of the system, as well as the different distribution functions are obtained as ensemble averages in the MC simulation. Typically, for each state an initial simulation of up to 2⫻106cycles 共depending on convergence兲is performed to equilibrate the system before averaging for some 1–2⫻105 cycles. Each cycle consists of NPattempts for random displacements and/or reorientations of the particles 共whether either or both types of moves are performed is chosen randomly兲, plus an attempt to change the volume 共actually done by rescaling the molecular diameter, ␴ , in units of box length兲. The maximum tilt angle and displacement of the particles and the maximum volume change are adjusted as to give an acceptance ratio of between 30% and 40%. The usual periodic boundary conditions and minimum image conventions are employed.18 B. Entropy criterion for liquid crystalline transitions As commented in the Introduction, we have explored the correlation between the configurational entropy and the phase diagram of our liquid crystal models. Following, among others, the works of Lazaridis et al.13 and Costa et al.,17 we have focused on the behavior of the different components of the pair entropy, s2, and of the residual entropy, ⌬s⫽sex⫺s2, defined as the multiparticle contribution to the excess entropy, sex , i.e., all correlations involving more than two particles. As shown in these previous works,13 the formal factorization of the correlation function, g(r, ␻ ) ⫽g(r)•g( ␻ 兩 r), in terms of g(r), the radial distribution function, and g( ␻ 兩 r), the conditional distribution function of two molecules with orientations described by ␻ 共which denotes all the relevant orientation angles兲at an intermolecular distance r, leads to a decomposition of the pair entropy in orientational and positional contributions. For the computation of ⌬sin a fluid of linear molecules, Costa et al. proposed a reduced dimensionality of the correlation function, g(r, ␪ ), in which the relative orientation of each pair of particles is described by the angle ␪ between the two molecular axis, already introduced in Sec. IIA. With this assumption, the formal factorization of the correlation function, g(r, ␪ )⫽g(r)•g( ␪ 兩 r), leads to the following expression for the pair contribution to the excess entropy, s2共per particle and in units of the Boltzmann constant兲: s2⫽s2 tr⫹s2 or ,共3兲 s2 tr⫽⫺2 ␲␳ 冕 关g共r兲lng共r兲⫺g共r兲⫹1兴r2dr,共4兲 s2 or⫽4 ␲␳ 冕 g共r兲Sor共r兲r2dr,共5兲 Sor共r兲⫽⫺ 1 4 冕 0 ␲ g共 ␪ 兩 r兲lng共 ␪ 兩 r兲sin ␪ d ␪ .共6兲 In the above expressions 共3兲–共6兲, ␳ ⫽N/Vdenotes the number density and s2 tr and s2 or refer to a formal distinction between translational and orientational contributions to s2; whereas s2 tr is expected to monitor any changes in the translational order of the molecules of the fluid, s2 or should respond to their orientational order. It must be remarked that the formal decomposition of the pair entropy s2⫽s2 tr⫹s2 or ,as well as Eq. 共4兲for s2 tr , are independent of the approximation introduced by the use of the correlation function g(r, ␪ ), whose justification relies upon the comparatively simplified computation of s2 or , as defined in Eq. 共6兲, and its appropriate behavior when applied to linear nematogens. In fact, Costa et al. found a rapid decrease of s2 or towards large negative values in the vicinity of the isotropic–nematic of the HSC fluid and suggested the zero of the residual entropy, ⌬s ⫽sex⫺s2⫽0, as an indicator for the incipient nematic ordering.17 The total excess entropy of the isotropic fluid may be evaluated from the MC equation of state, by means of the exact expression, sex共 ␳ 兲⫽uex T⫺ 冕 0 ␳ 冋 P kT ␳ ⬘⫺1 册 d ␳ ⬘ ␳ ⬘,共7兲 where uex denotes the excess energy 共i.e., the potential energy兲per particle in units of the Boltzmann constant. Hence, the excess entropy is determined by the internal energy and by the integral of (Z⫺1)/ ␳ , where Z⫽P/(kT ␳ ) is the compressibility factor of the fluid. In the present work, we have extended the computation of the excess entropy to the nematic and smectic phases by including the entropy change at the corresponding phase transitions. Since we are dealing with transitions at constant N, P, T, the entropy of transition, ⌬strans , can be readily obtained from the enthalpy of transition, ⌬htrans , as given by the changes in the internal energy, ⌬utrans , and volume ⌬vtrans , ⌬strans⫽⌬htrans T⫽⌬utrans T⫹P*⌬vtrans ,共8兲 where we recall that ⌬strans ,⌬htrans ,⌬utrans , and P*⌬vtrans(⫽P•⌬vtrans /kT) are per particle and in kunits. Thus, for instance, the extension of Eq. 共7兲to the calculation of the excess entropy of a particular state after the I–N phase transition is given by sex共 ␳ 兲⫽uex T⫺ 冕 0 ␳ 1 冋 P kT ␳ ⬘⫺1 册 d ␳ ⬘ ␳ ⬘⫹PI–N *⌬vI–N ⫺ 冕 ␳ 2 ␳ 冋 P kT ␳ ⬘⫺1 册 d ␳ ⬘ ␳ ⬘,共9兲 2936 J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Cuetos et al. where ␳ 1and ␳ 2, denote the coexistence densities of the two phases, represented in our NPT simulations by the densities of the states closest on either side to the phase boundary. The transition pressure is taken as the average of the simulation pressure of both boundary states. Note that in Eq. 共9兲the energetic part of the entropy of transition, ⌬uI–N/T, is effectively included in uex /T. The computation of the first integral in Eq. 共9兲was performed from the Monte Carlo values for (Z⫺1)/ ␳ and their linear extrapolation at low density. This latter extrapolation should provide, at ␳ ⫽0, the second virial coefficient, B2, of the fluid at the title temperature. In fact, as shown below 共Fig. 7兲, the extrapolation to ␳ ⫽0of the computed (Z⫺1)/ ␳ curve for the SWSC system is consistent with the known analytical expression of B2for this model fluid.6 We have employed Eqs. 共3兲–共9兲to extend the study of the behavior of sex ,⌬sand the different components of s2, through both the isotropic–nematic and the nematic–smectic transitions of the SWSC and SRS fluids. III. RESULTS AND DISCUSSION We have performed MC-NPT simulations for the SWSC and SRS fluids of elongation L*⫽5 at a reduced temperature T*⫽5. The thermodynamic data resulting in the compression and expansion simulation runs for this isotherm of both model fluids are listed in Tables I–IV. In addition, the equations of state 共P*vs ␳ *兲, order parameters and internal energies obtained for both models are represented in Figs. 2, 4, and 6. Within the range of pressures covered in our study, P*⬇0.1–2.2, both the SWSC and the SRS fluids exhibit differentiated isotropic, nematic, and smectic A phases. The TABLE I. Isothermal–isobaric Monte Carlo 共MC-NPT兲simulation results for the equation of state of the SWSC fluid of molecular elongation L*⫽L/ ␴ ⫽5, square well depth ␧, and range ␭*⫽␭/ ␴ ⫽1.5 共see Fig. 1兲 at temperature T*⫽kT/␧⫽5. The pressure, fixed in the simulations, is expressed in reduced units either with respect to the molecular diameter, ␴ 3共first column兲, or to the volume of the molecular hard core, vHSC 共second column兲, to allow for a direct comparison with earlier works. U*and Sdenote, respectively, the averages of potential energy per particle and the order parameter of the fluid. The simulations were run by compressing the fluid from the least dense state. The values in brackets denote the statistical uncertainty 共one standard deviation兲 in the last digit. The isotropic 共I兲, nematic 共N兲or smetic A 共Sm A兲phase corresponding for each state is indicated in the last column. P*⫽P• ␴ 3/kT P*•vHSC / ␴ 3 ␳ *⫽ ␳ • ␴ 3 ␩ ⫽ ␳ •vHSC U*⫽U/␧SPhase 0.01 0.045 0.0080共3兲0.035共1兲⫺0.36共2兲0.030共1兲I 0.02 0.089 0.0135共4兲0.060共2兲⫺0.64共3兲0.030共1兲I 0.05 0.22 0.0245共4兲0.109共2兲⫺1.25共4兲0.032共1兲I 0.10 0.45 0.0360共4兲0.160共2兲⫺2.00共4兲0.032共1兲I 0.20 0.89 0.0494共4兲0.220共2兲⫺2.98共5兲0.033共1兲I 0.30 1.34 0.0580共4兲0.258共2兲⫺3.66共5兲0.036共1兲I 0.40 1.78 0.0645共4兲0.287共2兲⫺4.20共5兲0.037共1兲I 0.50 2.23 0.0699共6兲0.311共3兲⫺4.63共6兲0.042共1兲I 0.60 2.67 0.0741共4兲0.330共2兲⫺4.97共5兲0.040共1兲I 0.70 3.12 0.0773共4兲0.344共2兲⫺5.21共5兲0.049共1兲I 0.80 3.56 0.0816共4兲0.363共2兲⫺5.52共5兲0.042共1兲I 0.90 4.01 0.0845共4兲0.376共2兲⫺5.69共5兲0.060共1兲I 1.00 4.45 0.0876共4兲0.390共2兲⫺5.89共5兲0.060共1兲I 1.10 4.90 0.0910共4兲0.405共2兲⫺6.02共5兲0.072共1兲I 1.20 5.34 0.0932共4兲0.415共2兲⫺6.13共5兲0.088共2兲I 1.25 5.56 0.0948共4兲0.422共2兲⫺6.18共5兲0.068共2兲I 1.30 5.79 0.0959共4兲0.427共2兲⫺6.20共5兲0.115共1兲I 1.35 6.01 0.0980共4兲0.436共2兲⫺6.21共4兲0.233共1兲I 1.40 6.23 0.0993共4兲0.442共2兲⫺6.18共6兲0.290共1兲I 1.45 6.45 0.1016共4兲0.452共2兲⫺6.20共5兲0.486共2兲N 1.50 6.68 0.1031共4兲0.459共2兲⫺6.13共5兲0.533共2兲N 1.55 6.90 0.1047共3兲0.466共1兲⫺6.08共4兲0.665共2兲N 1.60 7.12 0.1065共3兲0.474共1兲⫺6.05共4兲0.719共2兲N 1.65 7.34 0.1085共4兲0.483共2兲⫺5.98共4兲0.784共1兲N 1.70 7.57 0.1132共6兲0.504共3兲⫺5.79共4兲0.882共2兲Sm A 1.75 7.79 0.1146共4兲0.510共2兲⫺5.80共5兲0.890共1兲Sm A 1.80 8.01 0.1153共8兲0.513共2兲⫺5.86共5兲0.852共1兲Sm A 1.85 8.23 0.1180共4兲0.525共2兲⫺5.81共5兲0.880共2兲Sm A 1.90 8.46 0.1202共4兲0.535共2兲⫺5.77共5兲0.905共3兲Sm A 1.95 8.68 0.1225共6兲0.545共3兲⫺5.73共5兲0.914共2兲Sm A 2.00 8.90 0.1240共4兲0.552共2兲⫺5.74共5兲0.915共3兲Sm A 2.05 9.12 0.1251共4兲0.557共2兲⫺5.77共5兲0.906共3兲Sm A 2.10 9.35 0.1263共4兲0.562共2兲⫺5.81共6兲0.906共4兲Sm A 2.15 9.57 0.1272共4兲0.566共2兲⫺5.82共5兲0.917共3兲Sm A 2.20 9.79 0.1285共4兲0.572共2兲⫺5.82共5兲0.920共3兲Sm A 2937J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Liquid crystal phases of hard and soft spherocylinders phase transitions involve relatively weak discontinuous changes in density 共see Fig. 6兲and internal energy but are accompanied by appreciable jumps in the order parameter and/or in the structure of the different radial correlation functions, g(r), g2(r), and in g 储 (r 储 ), defined above and depicted in Figs. 3 and 5. As we shall see below, the pair and residual entropies also respond to the orientational and translational ordering of the fluid 共Figs. 7 and 8兲. TABLE II. Same as Table I, but for a set of simulations run by expanding the SWSC fluid from the most dense state (P*⫽2.20). P*⫽P• ␴ 3/kT P*•vHSC / ␴ 3 ␳ *⫽ ␳ • ␴ 3 ␩ ⫽ ␳ •vHSC U*⫽U/␧SPhase 1.00 4.45 0.0881共6兲0.392共3兲⫺5.87共7兲0.085共1兲I 1.10 4.90 0.0905共4兲0.403共2兲⫺6.03共5兲0.084共1兲I 1.20 5.34 0.0930共4兲0.414共2兲⫺6.11共5兲0.146共2兲I 1.25 5.56 0.0953共6兲0.424共3兲⫺6.15共5兲0.175共2兲I 1.30 5.79 0.0968共6兲0.431共3兲⫺6.17共5兲0.215共2兲I 1.35 6.01 0.0993共4兲0.442共2兲⫺6.02共5兲0.557共2兲N 1.40 6.23 0.1007共4兲0.448共2兲⫺6.04共5兲0.609共2兲N 1.45 6.45 0.1027共8兲0.457共4兲⫺5.98共6兲0.691共2兲N 1.50 6.68 0.1043共6兲0.464共3兲⫺5.99共5兲0.716共2兲N 1.55 6.90 0.1058共6兲0.471共3兲⫺5.97共5兲0.772共3兲N 1.60 7.12 0.1078共6兲0.480共3兲⫺6.00共5兲0.789共3兲N 1.65 7.34 0.1099共5兲0.489共3兲⫺5.92共5兲0.824共3兲N 1.70 7.57 0.1155共6兲0.514共3兲⫺5.67共3兲0.907共3兲Sm A 1.75 7.79 0.1171共5兲0.521共3兲⫺5.70共3兲0.913共3兲Sm A 1.80 8.01 0.1189共4兲0.529共2兲⫺5.70共5兲0.911共3兲Sm A 1.90 8.46 0.1218共4兲0.542共2兲⫺5.72共4兲0.915共3兲Sm A 2.00 8.90 0.1243共6兲0.553共3兲⫺5.80共5兲0.915共3兲Sm A 2.10 9.35 0.1263共4兲0.562共2兲⫺5.81共5兲0.926共3兲Sm A 2.20 9.79 0.1285共4兲0.572共2兲⫺5.82共5兲0.920共3兲Sm A TABLE III. Isothermal–isobaric Monte Carlo 共MC-NPT兲simulation results for the equation of state of the SRS fluid 关Eq. 共2兲,Fig.1兴of molecular elongation L*⫽L/ ␴ ⫽5 at temperature T*⫽kT/␧⫽5. The simulations were run by compressing the fluid from the least dense state. The same notation is employed as in Table I. In particular, vHSC denotes the volume of a hard spherocylinder with the same elongation L*⫽5. P*⫽P• ␴ 3/kT P*•vHSC / ␴ 3 ␳ *⫽ ␳ • ␴ 3 ␩ ⫽ ␳ •vHSC U*⫽U/␧SPhase 0.01 0.044 0.0079共4兲0,035共2兲0.12共3兲0.031共1兲I 0.02 0.089 0.0134共4兲0,060共2兲0.23共4兲0.030共1兲I 0.05 0.22 0.0249共4兲0,111共2兲0.48共6兲0.031共1兲I 0.10 0.45 0.0368共4兲0.164共2兲0.83共7兲0.033共1兲I 0.20 0.89 0.0521共6兲0.232共3兲1.4共1兲0.035共1兲I 0.30 1.34 0.0627共6兲0.279共3兲1.9共1兲0.037共1兲I 0.40 1.78 0.0708共6兲0.315共3兲2.4共1兲0.040共1兲I 0.50 2.23 0.0775共6兲0.345共3兲2.9共1兲0.040共1兲I 0.70 3.12 0.0887共6兲0.395共3兲3.7共1兲0.048共1兲I 0.80 3.56 0.0937共6兲0.417共3兲4.1共1兲0.055共1兲I 0.90 4.01 0.0982共6兲0.437共3兲4.5共1兲0.061共1兲I 1.00 4.45 0.1025共6兲0.456共3兲4.8共1兲0.071共1兲I 1.10 4.90 0.1065共6兲0.474共3兲5.2共2兲0.134共2兲I 1.15 5.12 0.1088共8兲0.484共4兲5.4共2兲0.110共3兲I 1.20 5.34 0.1117共6兲0.497共3兲5.5共2兲0.198共2兲I 1.25 5.56 0.1157共6兲0.515共3兲5.7共2兲0.541共2兲N 1.30 5.79 0.1200共8兲0.534共4兲5.8共2兲0.722共2兲N 1.35 6.01 0.1216共8兲0.541共4兲6.0共2兲0.715共2兲N 1.40 6.23 0.1240共8兲0.552共4兲6.2共2兲0.752共2兲N 1.45 6.45 0.1258共8兲0.560共4兲6.3共2兲0.766共2兲N 1.50 6.68 0.1281共8兲0.570共4兲6.5共2兲0.795共2兲N 1.55 6.90 0.1303共8兲0.580共4兲6.6共2兲0.820共3兲N 1.60 7.12 0.1319共8兲0.587共4兲6.8共2兲0.826共4兲N 1.65 7.34 0.1339共8兲0.596共4兲7.0共2兲0.838共3兲N 1.70 7.57 0.1360共8兲0.605共4兲7.1共2兲0.855共3兲N 1.75 7.79 0.1370共8兲0.610共4兲7.3共2兲0.862共3兲N 1.80 8.01 0.1398共8兲0.622共4兲7.4共2兲0.880共4兲N 1.85 8.23 0.145共1兲0.644共5兲7.4共2兲0.915共3兲N 1.90 8.46 0.147共1兲0.655共5兲7.5共2兲0.920共2兲Sm A 1.95 8.68 0.151共1兲0.673共6兲7.5共2兲0.937共2兲Sm A 2.00 8.90 0.153共1兲0.681共5兲7.6共2兲0.940共3兲Sm A 2.10 9.35 0.158共1兲0.702共5兲7.9共2兲0.950共3兲Sm A 2938 J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Cuetos et al. A. Phase diagram of the SWSC fluid We focus now on the discussion of the simulation results for the SWSC fluid. The isotherm obtained by compressing the system from the isotropic phase 共Table I and Fig. 2兲 evolves smoothly with increasing pressure up to P*⬇1.30. In this low pressure interval, both the density and the absolute value of the potential energy become progressively larger, while the order parameter remains small (S⬍0.1) and roughly constant. At P*⫽1.35–1.40, the number density is sufficiently high and close to the incipient isotropic–nematic transition as to induce fluctuations in the system that take the order parameter above 0.2. At the same time, the internal energy approaches an extremum and stabilizes at U* ⬇⫺6.2. In fact, at P*⫽1.45 the order parameter jumps to S⫽0.486 and enters already the nematic phase of the fluid; this is the first state included in the compression run that displays long-range orientational order, as monitored by g2(r). The density, internal energy and order parameter of this first nematic state, as well as the behavior of the thermodynamic variables in the proximity of the transition, are consistent with the result of similar simulations of Williamson and del Rı ´o6performed over a larger number of particles (NP⫽1020). Within the nematic phase, the further increase of pressure and density of the SWSC fluid induces a greater orientational ordering which takes the order parameter to S⬎0.7 for P*⬎1.55. At the same time the potential energy becomes slightly less negative 共i.e., the energy increases兲. The density, energy and order parameter maintain their increasing trend with growing P*, and at P*⫽1.70 the translational ordering characteristic of a smectic A phase develops and persists in the remaining higher pressure states included in our study. No evidence for smectic B or a solid phase order within each of the smectic layers was detected, which is not unexpected, since for the HSC fluid of same elongation the smectic FIG. 2. MC-NPT results for the equation of state 共top兲, order parameter 共middle兲, and internal energy per particle 共bottom兲for the SWSC fluid model. Solid circles and open triangles denote the simulation series run by compressing and expanding the fluid, respectively 共see text兲. The vertical dashed lines indicate the coexistence densities at the isotropic–nematic 共I–N兲and at the nematic–smectic A 共N–Sm A兲phase transitions, obtained in the simulations when expanding the fluid 共see Table V兲. TABLE IV. Same as Table III, but for a set of simulations run by expanding the SRS fluid from the most dense state (P*⫽2.10). P*⫽P• ␴ 3/kT P*•vHSC / ␴ 3 ␳ *⫽ ␳ • ␴ 3 ␩ ⫽ ␳ •vHSC U*⫽U/␧SPhase 1.00 4.45 0.1025共6兲0.456共3兲4.8共1兲0.058共2兲I 1.10 4.90 0.1067共6兲0.475共3兲5.2共2兲0.096共2兲I 1.20 5.34 0.1119共8兲0.498共4兲5.5共2兲0.267共2兲I 1.25 5.56 0.1157共8兲0.515共4兲5.7共2兲0.504共2兲N 1.30 5.79 0.1195共8兲0.532共4兲5.8共2兲0.698共2兲N 1.35 6.01 0.1220共8兲0.543共4兲6.0共2兲0.722共2兲N 1.40 6.23 0.1240共8兲0.552共4兲6.2共2兲0.740共2兲N 1.50 6.68 0.1285共8兲0.572共4兲6.5共2兲0.802共2兲N 1.60 7.12 0.1323共8兲0.589共4兲6.8共2兲0.830共2兲N 1.65 7.34 0.1343共8兲0.598共4兲7.0共2兲0.823共3兲N 1.70 7.57 0.1371共8兲0.610共4兲7.1共2兲0.864共3兲N 1.75 7.79 0.1387共8兲0.617共4兲7.3共2兲0.879共3兲N 1.80 8.01 0.1431共8兲0.637共5兲7.4共2兲0.904共3兲Sm A 1.85 8.23 0.1465共8兲0.652共5兲7.5共2兲0.931共3兲Sm A 1.90 8.46 0.1492共8兲0.664共6兲7.7共2兲0.938共3兲Sm A 2.00 8.90 0.1557共8兲0.693共6兲7.8共2兲0.953共3兲Sm A 2.10 9.35 0.1577共8兲0.702共5兲7.9共2兲0.950共3兲Sm A 2939J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Liquid crystal phases of hard and soft spherocylinders A–solid transition is observed at significantly higher pressure (P*⬇2.5) and density ( ␳ *⬇0.138).3 As mentioned above, in order to assess the stability of the different phases observed in the simulations, an expansion run was performed from the smectic state of highest pressure (P*⫽2.20) back to the isotropic region (P* ⫽1.00). The relevant data concerning the expansion simulation set are listed in Table II and also shown in Fig. 2. As the system pressure is decreased in this expansion run, the fluid undergoes the reverse sequence of phase transitions Sm A–N–I and a moderate hysteresis is observed. The isotropic–nematic and the nematic–smectic A phase boundaries are in this case located within the intervals (P* ⫽1.30–1.35, ␳ *⫽0.0968–0.0993) and (P*⫽1.65–1.70, ␳ *⫽0.1099–0.1155), respectively 共see Table V兲, and, hence, appear slightly shifted toward smaller pressures and/or densities in comparison to the compression run discussed in the precedent paragraphs. Similar effects were observed in the simulations of Williamson and del Rı ´o for the isotropic–nematic transition of this same fluid.6Since the equilibration of an ordered state from an initial disordered configuration is more demanding than the reverse case when performing MC simulations, we conclude that, in the compression run, the states of higher density within the isotropic and nematic phases are actually metastable, i.e., unstable with respect to the more ordered nematic and smectic phases, respectively. Therefore, when discussing the liquid crystalline behavior of the SWSC fluid in the remaining of the paper, we will consider solely the results arising from the simulations of the expansion run. Figure 3 represents the radial functions g(r), g2(r), and g 储 (r 储 )共see above for definitions兲for relevant states of the expansion MC run for the SWSC fluid. It is remarkable, for instance, how g2(r) clearly reflects, for P*⭓1.35, the long range orientational order that differentiates qualitatively the nematic and smectic phases from the isotropic one. On the other hand, the smectic order that characterizes the states at P*⭓1.70 in the expansion run is clearly observed not only in g 储 (r 储 ) but also in the radial distribution function g(r). As can be seen in Fig. 3, with growing density all through the isotropic and nematic phases, g(r) develops progressively more differentiated maxima associated to the first and second neighboring molecules 共at r/ ␴ ⬇1.1 and 2.3, respectively兲. However, after the nematic–smectic A transition 共i.e., from P*⫽1.65 to P*⫽1.70兲, the structure of g(r) changes qualitatively, with a sudden increase of the area of the first, second, and even third, next-neighbor peaks, and the appearance of a depression at r/ ␴ ⬇3.5–6.0, where g(r) stays below unity, as a consequence of the layered structure of the fluid. It must be noted, however, that, due to the weak first-order character of the liquid crystal transitions presently studied, and in spite of the qualitative changes undergone by g(r), g2(r), and g 储 (r 储 ) at the I–N and N–Sm A phase transitions, it is not necessarily straightforward to assign the boundary states for each transition in MC simulations unavoidably performed over a finite number of particles. For instance, the assignment of the N–Sm A transition, based on the observation of an appreciable jump in density 共see Fig. 6兲accompanied by a simultaneous jump in the amplitude of the layering maxima in g 储 (r 储 ) and in the structure of g(r), leaves a weak onset of layering 关⬇1.2 amplitude in g 储 (r 储 )兴in the state of highest density considered nematic in our study 共P*⫽1.65 in the SWSC fluid, P*⫽1.75 in the SRS fluid, see bottom of Figs. 3 and 5兲. This same effect was noted by McGrother et al. in their study of the HSC fluid.3It seems also timely to comment on the apparent more pronounced oscillations of g 储 (r 储 ) at the edge of the simulation cell, whereas the correct structure of this function expected in the smectic phase would be a sequence of maxima of the same height. However, due to the limited number of particles employed in our simulations, in the computation of g 储 (r 储 ) the first neighboring smectic layers are not sampled as efficiently as the central layer, and thus the computed maxima at 兩 r 储 / ␴ 兩 ⬇6.5 are too narrow and overestimated in height. Simulations with a larger number of particles 共and, thus, an increased size of the simulation cell兲, as those of Ref. 3 for the HSC system with L*⫽5, largely correct for this effect. Summarizing, the pressures and coexistence densities for the isotropic–nematic and nematic–smectic A transitions for FIG. 3. Correlation functions for representative states of the SWSC fluid in the present study 共see Sec. II A for definitions兲. The system pressure of each state is indicated next to the corresponding curve. 共Top兲radial distribution function g(r); 共middle兲orientational distribution function g2(r) ⬅ 具 P2(cos ␪ ) 典 ;共bottom兲projection of the pair distribution function a distance r 储 parallel to the nematic director g 储 (r 储 ). Oscillations in this latter function are indicative of layered smectic ordering in the fluid. 2940 J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Cuetos et al. the SWSC fluid, as obtained by averaging the pressures of the two boundary states of each phase transition 共in the expansion MC run兲and from their individual densities, are PI–N *⫽1.325; ␳ I *⫽0.0968, ␳ N *⫽0.0993, and PN–Sm A * ⫽1.675; ␳ N *⫽0.1109, ␳ Sm A *⫽0.1155, respectively 共see Table V兲. The specific role of the square well in the liquid crystal behavior of the SWSC fluid may be assessed by contrasting the present results with the MC-NPT simulations of McGrother et al. for the hard spherocylinder HSC fluid.3We begin by recalling the remarkable stabilization of the internal energy of the SWSC fluid at roughly constant values around U*⫽⫺6.0 at high density ( ␳ *⬎0.09) 共lower panel of Fig. 2兲. Since such stabilization takes place right before the I–N transition, it might be tentatively concluded that the energetic contribution to the free energy of the system ‘‘saturates’’ and has a negligible effect in the liquid crystalline behavior of the fluid. In fact, such a saturation effect may be understood as a direct consequence of the square-well interaction imposed in the model: at high enough densities all nearest-neighbors are at distance closer than the SW range and thus move freely inside the square well without change of energy.19 Under this interpretation it follows that the SWSC fluid should virtually resemble the liquid crystalline behavior of the HSC fluid. However, we find that both the I–N and the N–Sm A phase transitions of the SWSC fluid are delayed toward higher densities and pressures with respect to the HSC fluid. This can be seen in Fig. 6 which depicts in more detail the equation of state for the HSC and the SWSC fluids in the vicinity of the I–N and the N–Sm A transitions. The corresponding transition pressures and densities are compared in Table V. For the HSC system, the I–N and N–Sm A transitions were determined to be within (PI–N *⫽1.19, ␳ I–N *⫽0.0914–0.0932) and (PN–SmA *⫽1.540, ␳ N–Sm A * ⫽0.1061–0.1095), respectively.3Thus, the energetic contribution of the attractive well to the free energy does seem to drive partly these ordering transitions. The presence of the square well imposes configurational constrains that stabilize the isotropic phase with respect to the nematic phase, and also this latter one with respect to the smectic A phase. It follows that the role of the internal energy opposes in this case that of the main driving force of the liquid crystalline transitions; the configurational entropy, a magnitude that is controlled by excluded volume effects 共the only relevant effect for the HSC fluid兲, mainly a competition between the phase space accessible for the rotation and translation of the molecules. We finally note that for arbitrarily high temperatures 共i.e., T→⬁兲 the phase diagram of the SWSC system should tend asymptotically to that of the HSC fluid of same elongation. Hence, with increasing temperature the location of the liquid crystal transitions of the former are expected to shift smoothly toward smaller pressures and densities. B. Phase diagram of the SRS fluid We concentrate now on the simulation results for the SRS fluid. The isotherms obtained when compressing the fluid from the isotropic phase 共Table III and Fig. 4兲or expanding the system back from the smectic phase 共Table IV and Fig. 4兲display a sequence of isotropic–nematic–smectic A phases involving the same qualitative behavior in the relevant structure parameters and correlation functions 共Fig. 5兲 as that observed for the SWSC system. The simulations for the SRS fluid showed a weaker hysteresis than the SWSC fluid, yielding more similar results in the compression and expansion runs. However, some degree of hysteresis is still appreciable and we will refer to the results of the expansion run in the following discussion of the liquid crystalline properties of the SRS fluid. There are quite significant quantitative differences between the SRS system with respect to the SWSC and HSC models. As can be readily observed in Figs. 4 and 6, for a given pressure, the soft core of the SRS interaction brings the fluid to considerably larger densities in comparison to the HSC 共Ref. 3兲and the SWSC fluids. In addition, with respect to these latter systems, the SRS fluid presents I–N and N–Sm A transitions under quite different pressure/ density conditions, namely within the intervals (PI–N *⫽1.20–1.25, ␳ I–N *⫽0.1119–0.1157) and (PN–Sm A * FIG. 4. MC-NPT results for the equation of state 共top兲, order parameter 共middle兲, and internal energy per particle 共bottom兲for the SRS fluid model. The same notation as in Fig. 2 is used. The isotropic–nematic 共I–N兲and the namatic–smectic A 共N–Sm A兲coexistence densities are indicated by vertical dashed lines 共see Table V兲. Note that, in spite of the soft core of the SRS fluid 关Eq. 共2兲兴, the packing fraction given in the upper axis 共for a more direct comparison with earlier works兲is defined with respect to vHSC , the volume of a hard spherocylinder of the same elongation L*⫽5. 2941J. Chem. Phys., Vol. 117, No. 6, 8 August 2002 Liquid crystal phases of hard and soft spherocylinders