Nematic-smectic transition of parallel hard spheroellipsoids Franz J. Vesely Citation: The Journal of Chemical Physics 141, 064109 (2014); doi: 10.1063/1.4892378 View online: http://dx.doi.org/10.1063/1.4892378 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/6?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Nonuniform liquid-crystalline phases of parallel hard rod-shaped particles: From ellipsoids to cylinders J. Chem. Phys. 129, 054907 (2008); 10.1063/1.2958920 Smectic, nematic, and isotropic phases in binary mixtures of thin and thick hard spherocylinders J. Chem. Phys. 124, 234904 (2006); 10.1063/1.2207141 Isotropic-nematic transition of hard rods immersed in random sphere matrices J. Chem. Phys. 121, 12067 (2004); 10.1063/1.1815294 The nematic-isotropic phase transition in semiflexible fused hard-sphere chain fluids J. Chem. Phys. 114, 3314 (2001); 10.1063/1.1340606 The nematic–isotropic phase transition in linear fused hard-sphere chain fluids J. Chem. Phys. 110, 11630 (1999); 10.1063/1.479102
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
THE JOURNAL OF CHEMICAL PHYSICS 141, 064109 (2014)
Nematic-smectic transition of parallel hard spheroellipsoids Franz J. Veselya) Computational Physics Group, Faculty of Physics, University of Vienna, Sensengasse 8, A-1090 Vienna, Austria
(Received 17 June 2014; accepted 25 July 2014; published online 11 August 2014) Spheroellipsoids are truncated ellipsoids with spherical end caps. If gradients are assumed to change smoothly at the junction of body and cap, the truncation height z0 determines the geometry uniquely. The resulting model particle has only two shape parameters, namely, the aspect ratio c/a of the basic ellipsoid and the cutoff z0 /a. These two parameters can be tuned to yield a continuous transformation between a pure ellipsoid and a spherocylinder. Since parallel hard spherocylinders display a nematic-smectic A phase transition, while ellipsoids do not, the influence of the particle shape on the possibility of a smectic phase may be investigated. A density functional analysis is used to detect the dividing line, in the (c/a, z0 /a) plane, between the presence and absence of the N-S transition. Since spheroellipsoids may be useful as generic model particles for anisotropic molecules, we provide a computationally efficient overlap criterion for a pair in a general, non-parallel configuration. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4892378] I. INTRODUCTION
The study of fluids consisting of anisotropic hard particles is now well advanced. The initially surprising richness of structural phases is well understood, and some early presumptions about the necessity of anisotropic attractive interactions have been disproved. In particular, simulation and theoretical studies have established that the hard repulsive interaction is sufficient to explain the existence of smectic phases in some hard rod systems.1, 2 An intriguing problem that still remains is the qualitative difference between the phase behavior of hard prolate ellipsoids and spherocylinders. The latter display a nematicsmectic transition, while the spheroids do not. If the nematic phase is “ideal,” i.e., made up of completely parallel rods, the absence of a smectic phase in ellipsoids is easily explained: scaling the (nematic) z axis by the aspect ratio of the particles transforms them into spheres – which obviously do not tend to arrange themselves in layers, except at crystalline density. Still, we do not know what detail of the particle shape – for example, what amount of non-sphericity – gives rise to layering. In the general case of freely rotating bodies the difference between spherocylinders and ellipsoids is even more pronounced. As Bolhuis and Frenkel have shown,2 spherocylinders can form six different phases, two of which are of the liquid crystal type, while three are solid states of various kinds of spatial and orientational order. On the other hand, prolate spheroids have recently evoked renewed interest due to the discovery of an unexpected high-density SM2 structure.3–6 In this situation, it is desirable to have a model particle that may be “tuned” in a continuous manner such that it will resemble either an ellipsoid or a spherocylinder – or something in between. More generally, a standard model a) [email protected]
0021-9606/2014/141(6)/064109/6/$30.00
for anisotropic molecules is still lacking, and it is difficult to compare simulation or theoretical results for cylinders,1 spherocylinders,2, 7 fused spheres,8 ellipto-cylinders,9 hard Gaussian overlap particles,10 etc. A few-parameter model that can be varied in shape would certainly help to reduce this overabundance of models. A recent attempt to address this problem is the study of Martínez-Raton and Velasco11 who investigated superellipses with varying exponents ε. Starting out from an ellipsoid with ε = 2, one can change the shape of this object such that for large ε it eventually resembles a cylinder with slightly rounded rims. Applying both MC simulation and density functional theory, these authors indeed observed a qualitative change in phase behavior when the exponent exceeded ε ≈ 2.5, which is quite near to the ellipsoid limit. While the specific question for the stability of a smectic phase may be attacked in terms of the superellipse with changeable shape, this geometrical object is probably not suited for the role of a standard model for elongated molecules. The reason is that we have no overlap criterion for pairs of non-parallel superellipses, which is a prerequisite both for MC and for density functional theory. In Sec. V, we will relate such a criterion for general spheroellipsoids with arbitrary positions and orientations.
II. SPHEROELLIPSOIDS
Let a, c denote the short and long semiaxes of an ellipsoid. Truncate the c axis at some distance z0 from its center and complete the body by spherical end segments. Requiring that the tangents change continuously at the junction, we find that the geometrical center of the caps must 2 2 be placed at ±zc with zc = z0 (1 − a /c ); their radius is rc = a 1 − (z02 /c2 )(1 − a 2 /c2 ). The total length of the resulting spheroellipsoid is L = 2l = 2(zc + rc ). The particle volume is v = vbody + 2vcap ,
141, 064109-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-2
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
with aspect ratio l/a = 3 (blue) and with v = vell (c/a = 3) (green) are depicted. For later reference we consider the density at close packing of a system of strictly parallel hard spheroellipsoids. The maximum packing density is ηcp = v/Vcp , where Vcp is the unit cell volume. Assuming a hexagonal configuration, the cell volume may be determined as follows: (a) In a hexagonal configuration of non-truncated ellipsoids the point of contact between a particle centered at the origin and its nearest neighbor in √ the adjacent upper layer has a z coordinate of z∗ = c 2/3. If the cutoff height z0 is larger than z* the geometry of the hexagonal unit cell is the same as in the√pure ellipsoid case, = 8a 2 c/ 2. with a cell volume of Vcp (b) If z0 < z* the contact occurs no more between the ellipsoidal trunks but between the end caps; in this case the unit cell volume changes to Vcp √ = 4a 2 3 (zc + rc2 − 4a 2 /3). III. STABILITY ANALYSIS FIG. 1. Spheroellipsoid: Parameters are a, c, and z0 ; all other quantities follow.
where vbody = 2π a 2 z0 (1 − z02 /3c2 ) and vcap = (2π/3)rc3 − π z(rc2 − (z)2 /3), with z ≡ z0 − zc = z0 a2 /c2 . In the parameter plane of z0 /a vs. c/a all possible spheroellipsoids are situated below the diagonal z0 = c (see Fig. 1). It may be interesting to define a sequence of shapes with fixed aspect ratio l/a. Given c, a, and l, the required value of z0 is given by z0 (c, a, l) = l[1 −
(1/ l 2 − 1/c2 )/(1/a 2 − 1/c2 )].
(1)
If for given (a, c) a certain particle volume v0 is desired, the necessary value of z0 must be determined by solving the equation v(z0 ) = v0 numerically. Figure 2 shows a number of isoaspect and iso-volume lines; also, several spheroellipsoids
A. General relations
The following density functional analysis was first applied to parallel hard cylinders in 1987,1 later to spherocylinders,7, 12 and generalized to mixtures of rods and spheres13–15 and of rods with different lengths or widths.16–18 Recently, Martínez-Raton and Velasco used the method to study aligned superellipsoids.11 Let ρ(r) denote the local number density, with ρ(r) = ρ¯ + δρ(r) and ρ¯ = N/V . At the level of the second virial approximation the excess free energy per particle is given by φex {ρ(r)} ≡ Fex {ρ(r)} /N kT = −(1/2N ) dr dr ρ(r) ρ(r )f (r , r ), where f (r , r ) is the Mayer overlap function for a pair of particles at the given positions. Since we are dealing with hard particles, f = 0 except for the overlap region of the two objects, where f = −1. Due to the local density variation δρ(r) the free energy differs from that of the homogeneous density system.
FIG. 2. Parameter plane of spheroellipsoids. The red diagonal is inhabited by the basic ellipsoids with z0 = c. Blue curves: lines of equal aspect ratio; green curves: lines of equal volume. The bodies on the right correspond to the blue and green points; they are related to the ellipsoid with c/a = 3 in that they have the same aspect ratio (blue) or volume (green); the leftmost bodies are identical to the generic ellipsoid.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-3
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
The footprint of a layered phase will be the appearance of a preferred periodic density variation of the form δρ(r) = ρ¯ a cos kz. Inserting this in φex {ρ(r)} we find for the second order variation of the excess free energy δ2 φex
a2 2 ρ¯ I (k), = 4
+ 2Icap (k), where
Z2 Ibody (k) = 2π A C 1 − 02 C + 4π A2 C
(2)
where I (k) = v (0) dr cos kz is the cosine transform of the exc excluded volume around a particle centered at the origin. The variation of the ideal free energy density is derived from 1 1 Fid {...} = φid {...} ≡ dr [ρ ¯ + ρ(r) ln ρ(r)] . N kT N (3) Again inserting the periodic density profile ρ(z) = ρ¯ (1 + a cos kz) we find for the second-order term in δφ id , 1 2 a2 δ2 φid = dr cos2 kz = a . (4) 2V 4 Combining all this we find for the quadratic term of the free energy variation a2 I (k) a2 [1 + ρI ¯ (k)] = 1+η , (5) δ2 φ = 4 4 v where η = ρv ¯ is the packing fraction. Instability of the spatially homogeneous phase sets in when at some packing density ηs and wave number ks the ideal free energy variation is cancelled by the excess term. Thus, ηs = −v/I (ks ), where ks minimizes I(k); the corresponding wavelength λs ≡ 2π /ks is then the natural period of the smectic lamelles. To account for the higher, so far neglected terms in the virial expansion of the excess free energy we apply the correction introduced by Parsons and Lee.19, 20 Equation (5) is modified by replacing η with the Carnahan-Starling term χ (η) ≡ (η − 3η2 /4)/(1 − η)2 , I (k) a2 1 + χ (η) . (6) δ2 φ = 4 v Requiring that χ (ηs )I (k)/v = −1 we arrive at a quadratic equation for ηs ; writing w ≡ v/I (k) we find √ (2w − 1) − 1 − w . (7) ηs = 2w − 3/2
B. Calculation of I(k)
For parallel spheroellipsoids the excluded volume around one particle is again a spheroellipsoid whose specifications are double the values of the particle parameters: axes A = 2a, C = 2c; truncation half-length Z0 = 2z0 , and similarly for the center position Zc and radius Rc of the cap. The cosine integral I(k) is then made up of a part referring to the truncated ellipsoidal body and to the spherical caps: I(k) = Ibody (k)
2
sin kZ0 kC
sin kZ0 − kZ0 cos kZ0 (kC)3
(8)
and, with ≡ Z0 − Zc , Z +R c c
Icap (k) ≡ π dz Rc2 − (z − Zc )2 cos kz Zc +
= π K1 cos kZc − π kK2 where
K1 ≡
Rc
sin kZc , k
dz Rc2 − z2 cos kz
(9)
(10)
sin kRc − kRc cos kRc sin k − k cos k −2 k3 k3 sin k
(11) − Rc2 − 2 k
=2
k K2 ≡ k
Rc
dz Rc2 − z2 sin kz
(12)
Rc sin kRc − sin k 2 + Rc cos kRc − 2 cos k k
cos kRc − cos k − Rc2 cos kRc − cos k . (13) −2 2 k All terms in these expressions are well-behaved when k → 0. = −2
IV. RESULTS
In the parameter plane of Figure 2, we wish to identify the dividing line between particles that can develop a smectic phase and those that cannot. Since the ellipsoids reside on the red diagonal and the smectogenic spherocylinders are at the far right, we have a general idea where that boundary must be. While stability analysis yields the density at which layering may occur, it does not discern, at least in the simple formulation used here, between smectic and crystalline phases. Experience with all sorts of elongated hard particles shows that the relative packing fraction η* ≡ ηs /ηcp provides a hint: typically the transition from the nematic to a smectic phase occurs at a relative density of 0.35−0.46.21–25 In contrast, if the crystalline phase preempts the smectic, this is usually indicated by a higher relative density at the transition.22, 24 Hence we will set the limit for acceptable transition densities at η* = 0.5. The variation of the transition density with shape may be appreciated from Figure 3. As a typical example we consider a spheroellipsoid with (c/a, z0 /a) = (16.0, 5.071). It is practically indistinguishable from the frequently studied spherocylinder with aspect ratio (L + D)/D = 6. Accordingly, the relative density at the layering instability is 0.422, close to
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-4
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
low this line will be able to form a smectic phase, while the region above this line describes more “ellipsoid-like” particles which have no layered phase. Simulation studies will be needed to test this expectation by probing the vicinity of the boundary line. For example, spheroellipsoids with a fixed axis ratio c/a = 6 should change their phase behavior as the second parameter z0 /a is increased from 2.8 to 3.2, the theoretical boundary being at z0 /a = 2.916. Alternatively, several shapes with the same aspect ratio – say, l/a = 6 as in the previous example – may be tested. V. OVERLAP CRITERION
If spheroellipsoids are to play a role as model particles, there must be a quick procedure to test for pairwise overlap. Here we present such a criterion for a pair of particles with arbitrary orientations and positions. FIG. 3. DFT results for particles with aspect ratio l/a = 6: λ/2l . . . layer distance; ηcp . . . close packing density; ηs . . . packing density at smectic transition (if any). The dotted lines show the limiting values for true spherocylinders. The arrow denotes the parameter value (c/a)crit = 9.70 for which the relative packing density at transition equals 0.5; the second parameter is then z0 /a = 5.21.
the spherocylinder value of 0.386; clearly, such a particle will be smectogenic. When we decrease c/a while maintaining the same aspect ratio – see Eq. (1) – the relative bifurcation density slowly rises until it reaches η* = 0.5 at (c/a, z0 /a) = (9.7, 5.210). Proceeding to even smaller values of c/a we arrive at particle shapes that resemble ellipsoids so much that they will not display a smectic phase. Figure 4 shows the locus of parameter pairs (c/a, z0 /a) for which the relative density at the nematic instability equals 0.5. It is to be expected that particles situated to the right and be-
A. Shape potential
To describe the surface of a general spheroellipsoid A we start from the usual representation of an ellipsoid by the implicit constraint equation σAell (r) ≡ (r − rA ) · A · (r − rA ) − 1 = 0, where r is a point on the ellipsoid surface, rA is the ellipsoid center, and A is the matrix A = R−1 · A0 · R, with R a rotation matrix, and A0 the unrotated ellipsoid matrix. Let cA denote a unit vector along the positive long semiaxis of A; then r± c = rA ± zc cA are the centers of the sphercap ical caps. Their surfaces are simply defined by σA ≡ (r ± ± − rc ) · C · (r − rc ) − 1 = 0 with C = (1/rc2 ) I. The two constraint equations may be combined to describe a spheroellipsoid: σA ≡ σAell if the projection of r − rA onto the long cap axis cA is smaller than z0 , and σA ≡ σA else, σA (r) = (r − rA ) · A · (r − rA ) − 1 if |sA (r)| ≡ | cA · (r − rA ) | < z0 ± = (r − r± c ) · C · (r − rc ) − 1
if |sA (r)| ≥ z0 .
FIG. 4. Estimated dividing line between smectogenic and non-smectogenic spheroellipsoids. For shapes left of the boundary, density functional theory predicts that a layering transition occurs, if at all, at a packing density η* > 0.5. Since smectic phases normally start at lower relative densities, these particles are probably unable to form a smectic. Shapes right of the dividing line resemble spherocylinders and thus are expected to support a smectic phase.
(14)
We note that the function σA (r) for arbitrary points r in space may be regarded as a “shape potential” function, with the equipotential surface σA (r) = 0 describing the spheroellipse surface. By construction, the function σA (r) = 0 is continuous everywhere. This is not true for other equipotential surfaces, which display sharp discontinuities at the joining of ellipsoid trunk and cap. To account for this we redefine the function σA (r) such that all equipotential surfaces cap are continuous. To this end we multiply σA by the factor cap γ ≡ |∇σAell |/|∇σA |, where both gradients are taken at the junction of body and cap, z = z0 . It is easy to see that γ = rc2 /a 2 . The final formula for the shape potential then reads σA (r) = (r − rA ) · A · (r − rA ) − 1 = γ [(r −
r± c )
· C · (r −
r± c )
− 1]
if |sA (r)| < z0 else.
(15)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-5
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
Accordingly, we have for the local gradient gA (r) ≡ ∇σA (r) of the shape potential if sA (r) < z0 gA (r) = 2 A · (r − rA ) = 2 γ C · (r − rc )
else.
(16)
The equipotential surfaces are now smooth; for any point outside of A, σA (r) is positive, and it is negative for points within the spheroellipsoid. Values of σ A vary between −1 at the particle center and → ∞ at large distances. B. Contact function
Given two spheroellipsoids A, B with arbitrary center positions rA,B and ellipsoid matrices A, B, we wish to decide whether the two objects overlap. Drawing on the ideas of Perram and Wertheim26 we define a “contact function” as follows: F (r, λ) ≡ λσA (r) + (1 − λ)σB (r).
(17)
This function is defined everywhere, and for each value of λ [0, 1] there is a unique “bridge point” rb where F (r, λ) is minimal, i.e., λgA (rb ) + (1 − λ)gB (rb ) = 0. Recalling the possible expressions for gA,B we find that the following cases may occur:
FIG. 5. Bridge curve according to Perram and Wertheim. The points rb on the red curve fulfill λgA (rb ) + (1 − λ)gB (rb ) = 0, with λ[0, 1]. The intersection with any of the spheroellipse contours is a “plummet point”; it provides an overlap criterion (see text.)
± rb (λ) = λr± c,A + (1 − λ)rc,B
provides a convenient criterion: the bodies overlap if and only if σB (rp ) is below zero. It should be noted that the plummet points are in general distinct from the “proxy points” on the surfaces of a pair of convex shapes. The latter are defined as those points where the two surfaces are closest to each other. They may be found by requesting the normals to the surfaces A and B to be antiparallel: gA ∝ −gB . The plummet and proxy points coincide only when the two surfaces are in contact.
(18)
provided that |sA (rb )| ≥ z0 and |sB (rb )| ≥ z0 , meaning that the projections of rb onto the long axes of A and B are both outside the ellipsoidal trunks. Otherwise, rb (λ) = [λA + (1 − λ)B ]−1 · [λA · rA + (1 − λ)B · rB ], (19) if |sA (rb )| < z0 and |sB (rb )| < z0 . Finally,
−1 rb (λ) = λA + (1 − λ)(γ /rc2 )I
(20) · λA · rA + (1 − λ) γ /rc2 r± c,B if |sA (rb )| < z0 and |sB (rb )| ≥ z0 . This situation occurs if rb lies “athwart” of the truncated long axis of A but not of B; similarly for interchanged B and A. The bridge points rb (λ) trace a curved path between the centers of A and B; see Figure 5 for an example in two dimensions. Perram and Wertheim proceeded by locating the (singular, by construction) maximum of F (λ) = λσA (rb (λ)) + (1 − λ)σB (rb (λ)). Obviously, the value of Fmax provides an overlap criterion: if it is above zero, then the respective bridge point must be situated outside of both ellipsoids, which is possible only if the ellipsoids are separate; similar considerations can be found for Fmax < 0 (overlap) and = 0 (contact). This method has been used extensively in Monte Carlo simulations of ellipsoidal models. Instead of locating the maximum of F(λ) – which entails a possibly time-consuming extremal search26–28 – we will determine the intersection of the bridge curve with the surface of A. This may be done applying a fast Newton-Raphson pro˜ = 0. It is obvious that the cedure to find λ˜ such that σA (r(λ)) ˜ being that point on the resulting “plummet point” rp = r(λ), surface of A which is “lowest” in the shape potential of B,
C. Application of the overlap criterion
Given spheroellipoids A, B and the bridge curve rb (λ) be˜ is on the tween their centers, we wish to find λ˜ such that rb (λ) ˜ = 0. As mentioned before, the resurface of A, i.e., σA (rb (λ)) spective “plummet point” is by construction that point on A ˜ and which is lowest in the shape potential of B: rp ≡ rb (λ), σB (rp ) = min(r on A) σB (r). Here is the procedure: 1. Given A and B, choose λ(1) = 0.5 and compute r(1) b on the bridge curve. 2. For the Newton-Raphson derivative we find
(1) drb dσA r(1) b (λ) , r · = g D≡ A b (1) dλ dλ λ(1) λ
(21)
where the local tangent vector drb /dλ may be calculated trivially from Eqs. (18) to (20). 3. The next estimate for λ˜ is then λ(2) = λ(1) − σA (r(1) p )/D. 4. Iterate until the change in λ(n) is negligible; set λ˜ = λ(n) ˜ and rp = rb (λ). Having found the plummet point rp in this manner, check for σB (rp ) : 0.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-6
Franz J. Vesely
VI. CONCLUSIONS
Spheroellipsoids promise to serve as a versatile model for a large group of elongated convex particles. By varying the two model parameters c/a nd z0 /a various shapes may be achieved, from simple ellipsoids to almost perfect spherocylinders. The large variation of mesogenic properties over this class of particles may thus be explored on the basis of a single few-parameter model. In particular, by applying Parsons-Lee level density functional theory to systems of parallel spheroellipsoids we have been able to suggest a dividing line between smectogenic and non-smectogenic particle shapes. In order to facilitate future MC simulations as well as density functional theory applied to freely reorienting spheroellipsoids, we have derived a fast overlap detection procedure. The method is based on Perram and Wertheim’s overlap criterion for ellipsoids, which was generalized to spheroellipsoids and accelerated by replacing an extremal search by a Newton-Raphson procedure.
1 B.
Mulder, Phys. Rev. A 35(7), 3095 (1987). Bolhuis and D. Frenkel, J. Chem. Phys. 106, 666 (1997). 3 A. Donev, F. H. Stillinger, P. M. Chaikin, and S. Torquato, Phys. Rev. Lett. 92, 255506 (2004). 4 P. Pfleiderer and T. Schilling, Phys. Rev. E 75, 020402(R) (2007). 5 G. Odriozola, J. Chem. Phys. 136(13), 134505 (2012). 6 G. Bautista-Carbajal, A. Moncho-Jordá, and G. Odriozola, J. Chem. Phys. 138(6), 064501 (2013). 7 A. Poniewierski and R. Holyst, Phys. Rev. Lett. 61(21), 2461 (1988). 2 P.
J. Chem. Phys. 141, 064109 (2014) 8 C.
Vega, C. McBride, and L. G. MacDowell, J. Chem. Phys. 115(9), 4203 (2001). 9 G. T. Evans, Mol. Phys. 76, 1359 (1992). 10 E. de Miguel and E. Martin del Río, J. Chem. Phys. 115, 9072 (2001). 11 Y. Martínez-Raton and E. Velasco, J. Chem. Phys. 129, 054907 (2008). 12 H. Xu, H. N. W. Lekkerkerker, and M. Baus, Europhys. Lett. 17(2), 163 (1992). 13 T. Koda and H. Kimura, J. Phys. Soc. Jpn. 63, 984 (1994). 14 T. Koda, M. Numajiri, and S. Ikeda, J. Phys. Soc. Jpn. 65, 3551 (1996). 15 F. J. Vesely, Mol. Phys. 103(5), 679 (2005). 16 Sz. Varga, E. Velasco, and F. J. Vesely, Mol. Phys. 107(23), 2481 (2009). 17 Sz. Varga, A. Gabor, E. Velasco, L. Mederos, and F. J. Vesely, Mol. Phys. 106(15), 1939 (2008). 18 Z. Dogic, D. Frenkel, and S. Fraden, Phys. Rev. E 62, 3925 (2000). 19 J. D. Parsons, Phys. Rev. A 19, 1225 (1979). 20 S.-D. Lee, J. Chem. Phys. 87(8), 4972 (1987). 21 A. Stroobants, H. N. W. Lekkerkerker, and D. Frenkel, Phys. Rev. Lett. 57(12), 1452 (1986). 22 J. A. C. Veerman and D. Frenkel, Phys. Rev. A 43(8), 4334 (1991). 23 D. Costa, F. Saija, and P. V. Giaquinta, J. Phys. Chem. B 107, 9514 (2003). 24 J. A. Capitán, Y. Martínez-Ratón, and J. A. Cuesta, J. Chem. Phys. 128, 194901 (2008). 25 M. Torikai, J. Stat. Phys. 148, 345–352 (2012). 26 J. W. Perram and M. S. Wertheim, J. Comput. Phys. 58, 409 (1985). 27 To locate F max , Perram and Wertheim proposed the use of a Brent-type extremal search which converges with order 1.3247. Donev et al.28 have suggested to use a Newton-Raphson search for the zero of dF(λ)/dλ. NR converges quadratically, which makes their approach as fast as ours. It should be mentioned that generally it is not recommended to replace a minimum search by a search for the zero of the derivative, but here the functions involved seem to be very well-behaved, and Donev also recommends to use a “safeguarded”NR scheme, meaning that in case of convergence problems the code would automatically switch to the slow but safe interval bisection. In our scheme the location of Fmax is not needed. 28 A. Donev, S. Torquato, and F. H. Stillinger, J. Comput. Phys. 202(2), 737 (2005).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
THE JOURNAL OF CHEMICAL PHYSICS 141, 064109 (2014)
Nematic-smectic transition of parallel hard spheroellipsoids Franz J. Veselya) Computational Physics Group, Faculty of Physics, University of Vienna, Sensengasse 8, A-1090 Vienna, Austria
(Received 17 June 2014; accepted 25 July 2014; published online 11 August 2014) Spheroellipsoids are truncated ellipsoids with spherical end caps. If gradients are assumed to change smoothly at the junction of body and cap, the truncation height z0 determines the geometry uniquely. The resulting model particle has only two shape parameters, namely, the aspect ratio c/a of the basic ellipsoid and the cutoff z0 /a. These two parameters can be tuned to yield a continuous transformation between a pure ellipsoid and a spherocylinder. Since parallel hard spherocylinders display a nematic-smectic A phase transition, while ellipsoids do not, the influence of the particle shape on the possibility of a smectic phase may be investigated. A density functional analysis is used to detect the dividing line, in the (c/a, z0 /a) plane, between the presence and absence of the N-S transition. Since spheroellipsoids may be useful as generic model particles for anisotropic molecules, we provide a computationally efficient overlap criterion for a pair in a general, non-parallel configuration. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4892378] I. INTRODUCTION
The study of fluids consisting of anisotropic hard particles is now well advanced. The initially surprising richness of structural phases is well understood, and some early presumptions about the necessity of anisotropic attractive interactions have been disproved. In particular, simulation and theoretical studies have established that the hard repulsive interaction is sufficient to explain the existence of smectic phases in some hard rod systems.1, 2 An intriguing problem that still remains is the qualitative difference between the phase behavior of hard prolate ellipsoids and spherocylinders. The latter display a nematicsmectic transition, while the spheroids do not. If the nematic phase is “ideal,” i.e., made up of completely parallel rods, the absence of a smectic phase in ellipsoids is easily explained: scaling the (nematic) z axis by the aspect ratio of the particles transforms them into spheres – which obviously do not tend to arrange themselves in layers, except at crystalline density. Still, we do not know what detail of the particle shape – for example, what amount of non-sphericity – gives rise to layering. In the general case of freely rotating bodies the difference between spherocylinders and ellipsoids is even more pronounced. As Bolhuis and Frenkel have shown,2 spherocylinders can form six different phases, two of which are of the liquid crystal type, while three are solid states of various kinds of spatial and orientational order. On the other hand, prolate spheroids have recently evoked renewed interest due to the discovery of an unexpected high-density SM2 structure.3–6 In this situation, it is desirable to have a model particle that may be “tuned” in a continuous manner such that it will resemble either an ellipsoid or a spherocylinder – or something in between. More generally, a standard model a) [email protected]
0021-9606/2014/141(6)/064109/6/$30.00
for anisotropic molecules is still lacking, and it is difficult to compare simulation or theoretical results for cylinders,1 spherocylinders,2, 7 fused spheres,8 ellipto-cylinders,9 hard Gaussian overlap particles,10 etc. A few-parameter model that can be varied in shape would certainly help to reduce this overabundance of models. A recent attempt to address this problem is the study of Martínez-Raton and Velasco11 who investigated superellipses with varying exponents ε. Starting out from an ellipsoid with ε = 2, one can change the shape of this object such that for large ε it eventually resembles a cylinder with slightly rounded rims. Applying both MC simulation and density functional theory, these authors indeed observed a qualitative change in phase behavior when the exponent exceeded ε ≈ 2.5, which is quite near to the ellipsoid limit. While the specific question for the stability of a smectic phase may be attacked in terms of the superellipse with changeable shape, this geometrical object is probably not suited for the role of a standard model for elongated molecules. The reason is that we have no overlap criterion for pairs of non-parallel superellipses, which is a prerequisite both for MC and for density functional theory. In Sec. V, we will relate such a criterion for general spheroellipsoids with arbitrary positions and orientations.
II. SPHEROELLIPSOIDS
Let a, c denote the short and long semiaxes of an ellipsoid. Truncate the c axis at some distance z0 from its center and complete the body by spherical end segments. Requiring that the tangents change continuously at the junction, we find that the geometrical center of the caps must 2 2 be placed at ±zc with zc = z0 (1 − a /c ); their radius is rc = a 1 − (z02 /c2 )(1 − a 2 /c2 ). The total length of the resulting spheroellipsoid is L = 2l = 2(zc + rc ). The particle volume is v = vbody + 2vcap ,
141, 064109-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-2
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
with aspect ratio l/a = 3 (blue) and with v = vell (c/a = 3) (green) are depicted. For later reference we consider the density at close packing of a system of strictly parallel hard spheroellipsoids. The maximum packing density is ηcp = v/Vcp , where Vcp is the unit cell volume. Assuming a hexagonal configuration, the cell volume may be determined as follows: (a) In a hexagonal configuration of non-truncated ellipsoids the point of contact between a particle centered at the origin and its nearest neighbor in √ the adjacent upper layer has a z coordinate of z∗ = c 2/3. If the cutoff height z0 is larger than z* the geometry of the hexagonal unit cell is the same as in the√pure ellipsoid case, = 8a 2 c/ 2. with a cell volume of Vcp (b) If z0 < z* the contact occurs no more between the ellipsoidal trunks but between the end caps; in this case the unit cell volume changes to Vcp √ = 4a 2 3 (zc + rc2 − 4a 2 /3). III. STABILITY ANALYSIS FIG. 1. Spheroellipsoid: Parameters are a, c, and z0 ; all other quantities follow.
where vbody = 2π a 2 z0 (1 − z02 /3c2 ) and vcap = (2π/3)rc3 − π z(rc2 − (z)2 /3), with z ≡ z0 − zc = z0 a2 /c2 . In the parameter plane of z0 /a vs. c/a all possible spheroellipsoids are situated below the diagonal z0 = c (see Fig. 1). It may be interesting to define a sequence of shapes with fixed aspect ratio l/a. Given c, a, and l, the required value of z0 is given by z0 (c, a, l) = l[1 −
(1/ l 2 − 1/c2 )/(1/a 2 − 1/c2 )].
(1)
If for given (a, c) a certain particle volume v0 is desired, the necessary value of z0 must be determined by solving the equation v(z0 ) = v0 numerically. Figure 2 shows a number of isoaspect and iso-volume lines; also, several spheroellipsoids
A. General relations
The following density functional analysis was first applied to parallel hard cylinders in 1987,1 later to spherocylinders,7, 12 and generalized to mixtures of rods and spheres13–15 and of rods with different lengths or widths.16–18 Recently, Martínez-Raton and Velasco used the method to study aligned superellipsoids.11 Let ρ(r) denote the local number density, with ρ(r) = ρ¯ + δρ(r) and ρ¯ = N/V . At the level of the second virial approximation the excess free energy per particle is given by φex {ρ(r)} ≡ Fex {ρ(r)} /N kT = −(1/2N ) dr dr ρ(r) ρ(r )f (r , r ), where f (r , r ) is the Mayer overlap function for a pair of particles at the given positions. Since we are dealing with hard particles, f = 0 except for the overlap region of the two objects, where f = −1. Due to the local density variation δρ(r) the free energy differs from that of the homogeneous density system.
FIG. 2. Parameter plane of spheroellipsoids. The red diagonal is inhabited by the basic ellipsoids with z0 = c. Blue curves: lines of equal aspect ratio; green curves: lines of equal volume. The bodies on the right correspond to the blue and green points; they are related to the ellipsoid with c/a = 3 in that they have the same aspect ratio (blue) or volume (green); the leftmost bodies are identical to the generic ellipsoid.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-3
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
The footprint of a layered phase will be the appearance of a preferred periodic density variation of the form δρ(r) = ρ¯ a cos kz. Inserting this in φex {ρ(r)} we find for the second order variation of the excess free energy δ2 φex
a2 2 ρ¯ I (k), = 4
+ 2Icap (k), where
Z2 Ibody (k) = 2π A C 1 − 02 C + 4π A2 C
(2)
where I (k) = v (0) dr cos kz is the cosine transform of the exc excluded volume around a particle centered at the origin. The variation of the ideal free energy density is derived from 1 1 Fid {...} = φid {...} ≡ dr [ρ ¯ + ρ(r) ln ρ(r)] . N kT N (3) Again inserting the periodic density profile ρ(z) = ρ¯ (1 + a cos kz) we find for the second-order term in δφ id , 1 2 a2 δ2 φid = dr cos2 kz = a . (4) 2V 4 Combining all this we find for the quadratic term of the free energy variation a2 I (k) a2 [1 + ρI ¯ (k)] = 1+η , (5) δ2 φ = 4 4 v where η = ρv ¯ is the packing fraction. Instability of the spatially homogeneous phase sets in when at some packing density ηs and wave number ks the ideal free energy variation is cancelled by the excess term. Thus, ηs = −v/I (ks ), where ks minimizes I(k); the corresponding wavelength λs ≡ 2π /ks is then the natural period of the smectic lamelles. To account for the higher, so far neglected terms in the virial expansion of the excess free energy we apply the correction introduced by Parsons and Lee.19, 20 Equation (5) is modified by replacing η with the Carnahan-Starling term χ (η) ≡ (η − 3η2 /4)/(1 − η)2 , I (k) a2 1 + χ (η) . (6) δ2 φ = 4 v Requiring that χ (ηs )I (k)/v = −1 we arrive at a quadratic equation for ηs ; writing w ≡ v/I (k) we find √ (2w − 1) − 1 − w . (7) ηs = 2w − 3/2
B. Calculation of I(k)
For parallel spheroellipsoids the excluded volume around one particle is again a spheroellipsoid whose specifications are double the values of the particle parameters: axes A = 2a, C = 2c; truncation half-length Z0 = 2z0 , and similarly for the center position Zc and radius Rc of the cap. The cosine integral I(k) is then made up of a part referring to the truncated ellipsoidal body and to the spherical caps: I(k) = Ibody (k)
2
sin kZ0 kC
sin kZ0 − kZ0 cos kZ0 (kC)3
(8)
and, with ≡ Z0 − Zc , Z +R c c
Icap (k) ≡ π dz Rc2 − (z − Zc )2 cos kz Zc +
= π K1 cos kZc − π kK2 where
K1 ≡
Rc
sin kZc , k
dz Rc2 − z2 cos kz
(9)
(10)
sin kRc − kRc cos kRc sin k − k cos k −2 k3 k3 sin k
(11) − Rc2 − 2 k
=2
k K2 ≡ k
Rc
dz Rc2 − z2 sin kz
(12)
Rc sin kRc − sin k 2 + Rc cos kRc − 2 cos k k
cos kRc − cos k − Rc2 cos kRc − cos k . (13) −2 2 k All terms in these expressions are well-behaved when k → 0. = −2
IV. RESULTS
In the parameter plane of Figure 2, we wish to identify the dividing line between particles that can develop a smectic phase and those that cannot. Since the ellipsoids reside on the red diagonal and the smectogenic spherocylinders are at the far right, we have a general idea where that boundary must be. While stability analysis yields the density at which layering may occur, it does not discern, at least in the simple formulation used here, between smectic and crystalline phases. Experience with all sorts of elongated hard particles shows that the relative packing fraction η* ≡ ηs /ηcp provides a hint: typically the transition from the nematic to a smectic phase occurs at a relative density of 0.35−0.46.21–25 In contrast, if the crystalline phase preempts the smectic, this is usually indicated by a higher relative density at the transition.22, 24 Hence we will set the limit for acceptable transition densities at η* = 0.5. The variation of the transition density with shape may be appreciated from Figure 3. As a typical example we consider a spheroellipsoid with (c/a, z0 /a) = (16.0, 5.071). It is practically indistinguishable from the frequently studied spherocylinder with aspect ratio (L + D)/D = 6. Accordingly, the relative density at the layering instability is 0.422, close to
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-4
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
low this line will be able to form a smectic phase, while the region above this line describes more “ellipsoid-like” particles which have no layered phase. Simulation studies will be needed to test this expectation by probing the vicinity of the boundary line. For example, spheroellipsoids with a fixed axis ratio c/a = 6 should change their phase behavior as the second parameter z0 /a is increased from 2.8 to 3.2, the theoretical boundary being at z0 /a = 2.916. Alternatively, several shapes with the same aspect ratio – say, l/a = 6 as in the previous example – may be tested. V. OVERLAP CRITERION
If spheroellipsoids are to play a role as model particles, there must be a quick procedure to test for pairwise overlap. Here we present such a criterion for a pair of particles with arbitrary orientations and positions. FIG. 3. DFT results for particles with aspect ratio l/a = 6: λ/2l . . . layer distance; ηcp . . . close packing density; ηs . . . packing density at smectic transition (if any). The dotted lines show the limiting values for true spherocylinders. The arrow denotes the parameter value (c/a)crit = 9.70 for which the relative packing density at transition equals 0.5; the second parameter is then z0 /a = 5.21.
the spherocylinder value of 0.386; clearly, such a particle will be smectogenic. When we decrease c/a while maintaining the same aspect ratio – see Eq. (1) – the relative bifurcation density slowly rises until it reaches η* = 0.5 at (c/a, z0 /a) = (9.7, 5.210). Proceeding to even smaller values of c/a we arrive at particle shapes that resemble ellipsoids so much that they will not display a smectic phase. Figure 4 shows the locus of parameter pairs (c/a, z0 /a) for which the relative density at the nematic instability equals 0.5. It is to be expected that particles situated to the right and be-
A. Shape potential
To describe the surface of a general spheroellipsoid A we start from the usual representation of an ellipsoid by the implicit constraint equation σAell (r) ≡ (r − rA ) · A · (r − rA ) − 1 = 0, where r is a point on the ellipsoid surface, rA is the ellipsoid center, and A is the matrix A = R−1 · A0 · R, with R a rotation matrix, and A0 the unrotated ellipsoid matrix. Let cA denote a unit vector along the positive long semiaxis of A; then r± c = rA ± zc cA are the centers of the sphercap ical caps. Their surfaces are simply defined by σA ≡ (r ± ± − rc ) · C · (r − rc ) − 1 = 0 with C = (1/rc2 ) I. The two constraint equations may be combined to describe a spheroellipsoid: σA ≡ σAell if the projection of r − rA onto the long cap axis cA is smaller than z0 , and σA ≡ σA else, σA (r) = (r − rA ) · A · (r − rA ) − 1 if |sA (r)| ≡ | cA · (r − rA ) | < z0 ± = (r − r± c ) · C · (r − rc ) − 1
if |sA (r)| ≥ z0 .
FIG. 4. Estimated dividing line between smectogenic and non-smectogenic spheroellipsoids. For shapes left of the boundary, density functional theory predicts that a layering transition occurs, if at all, at a packing density η* > 0.5. Since smectic phases normally start at lower relative densities, these particles are probably unable to form a smectic. Shapes right of the dividing line resemble spherocylinders and thus are expected to support a smectic phase.
(14)
We note that the function σA (r) for arbitrary points r in space may be regarded as a “shape potential” function, with the equipotential surface σA (r) = 0 describing the spheroellipse surface. By construction, the function σA (r) = 0 is continuous everywhere. This is not true for other equipotential surfaces, which display sharp discontinuities at the joining of ellipsoid trunk and cap. To account for this we redefine the function σA (r) such that all equipotential surfaces cap are continuous. To this end we multiply σA by the factor cap γ ≡ |∇σAell |/|∇σA |, where both gradients are taken at the junction of body and cap, z = z0 . It is easy to see that γ = rc2 /a 2 . The final formula for the shape potential then reads σA (r) = (r − rA ) · A · (r − rA ) − 1 = γ [(r −
r± c )
· C · (r −
r± c )
− 1]
if |sA (r)| < z0 else.
(15)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-5
Franz J. Vesely
J. Chem. Phys. 141, 064109 (2014)
Accordingly, we have for the local gradient gA (r) ≡ ∇σA (r) of the shape potential if sA (r) < z0 gA (r) = 2 A · (r − rA ) = 2 γ C · (r − rc )
else.
(16)
The equipotential surfaces are now smooth; for any point outside of A, σA (r) is positive, and it is negative for points within the spheroellipsoid. Values of σ A vary between −1 at the particle center and → ∞ at large distances. B. Contact function
Given two spheroellipsoids A, B with arbitrary center positions rA,B and ellipsoid matrices A, B, we wish to decide whether the two objects overlap. Drawing on the ideas of Perram and Wertheim26 we define a “contact function” as follows: F (r, λ) ≡ λσA (r) + (1 − λ)σB (r).
(17)
This function is defined everywhere, and for each value of λ [0, 1] there is a unique “bridge point” rb where F (r, λ) is minimal, i.e., λgA (rb ) + (1 − λ)gB (rb ) = 0. Recalling the possible expressions for gA,B we find that the following cases may occur:
FIG. 5. Bridge curve according to Perram and Wertheim. The points rb on the red curve fulfill λgA (rb ) + (1 − λ)gB (rb ) = 0, with λ[0, 1]. The intersection with any of the spheroellipse contours is a “plummet point”; it provides an overlap criterion (see text.)
± rb (λ) = λr± c,A + (1 − λ)rc,B
provides a convenient criterion: the bodies overlap if and only if σB (rp ) is below zero. It should be noted that the plummet points are in general distinct from the “proxy points” on the surfaces of a pair of convex shapes. The latter are defined as those points where the two surfaces are closest to each other. They may be found by requesting the normals to the surfaces A and B to be antiparallel: gA ∝ −gB . The plummet and proxy points coincide only when the two surfaces are in contact.
(18)
provided that |sA (rb )| ≥ z0 and |sB (rb )| ≥ z0 , meaning that the projections of rb onto the long axes of A and B are both outside the ellipsoidal trunks. Otherwise, rb (λ) = [λA + (1 − λ)B ]−1 · [λA · rA + (1 − λ)B · rB ], (19) if |sA (rb )| < z0 and |sB (rb )| < z0 . Finally,
−1 rb (λ) = λA + (1 − λ)(γ /rc2 )I
(20) · λA · rA + (1 − λ) γ /rc2 r± c,B if |sA (rb )| < z0 and |sB (rb )| ≥ z0 . This situation occurs if rb lies “athwart” of the truncated long axis of A but not of B; similarly for interchanged B and A. The bridge points rb (λ) trace a curved path between the centers of A and B; see Figure 5 for an example in two dimensions. Perram and Wertheim proceeded by locating the (singular, by construction) maximum of F (λ) = λσA (rb (λ)) + (1 − λ)σB (rb (λ)). Obviously, the value of Fmax provides an overlap criterion: if it is above zero, then the respective bridge point must be situated outside of both ellipsoids, which is possible only if the ellipsoids are separate; similar considerations can be found for Fmax < 0 (overlap) and = 0 (contact). This method has been used extensively in Monte Carlo simulations of ellipsoidal models. Instead of locating the maximum of F(λ) – which entails a possibly time-consuming extremal search26–28 – we will determine the intersection of the bridge curve with the surface of A. This may be done applying a fast Newton-Raphson pro˜ = 0. It is obvious that the cedure to find λ˜ such that σA (r(λ)) ˜ being that point on the resulting “plummet point” rp = r(λ), surface of A which is “lowest” in the shape potential of B,
C. Application of the overlap criterion
Given spheroellipoids A, B and the bridge curve rb (λ) be˜ is on the tween their centers, we wish to find λ˜ such that rb (λ) ˜ = 0. As mentioned before, the resurface of A, i.e., σA (rb (λ)) spective “plummet point” is by construction that point on A ˜ and which is lowest in the shape potential of B: rp ≡ rb (λ), σB (rp ) = min(r on A) σB (r). Here is the procedure: 1. Given A and B, choose λ(1) = 0.5 and compute r(1) b on the bridge curve. 2. For the Newton-Raphson derivative we find
(1) drb dσA r(1) b (λ) , r · = g D≡ A b (1) dλ dλ λ(1) λ
(21)
where the local tangent vector drb /dλ may be calculated trivially from Eqs. (18) to (20). 3. The next estimate for λ˜ is then λ(2) = λ(1) − σA (r(1) p )/D. 4. Iterate until the change in λ(n) is negligible; set λ˜ = λ(n) ˜ and rp = rb (λ). Having found the plummet point rp in this manner, check for σB (rp ) : 0.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
064109-6
Franz J. Vesely
VI. CONCLUSIONS
Spheroellipsoids promise to serve as a versatile model for a large group of elongated convex particles. By varying the two model parameters c/a nd z0 /a various shapes may be achieved, from simple ellipsoids to almost perfect spherocylinders. The large variation of mesogenic properties over this class of particles may thus be explored on the basis of a single few-parameter model. In particular, by applying Parsons-Lee level density functional theory to systems of parallel spheroellipsoids we have been able to suggest a dividing line between smectogenic and non-smectogenic particle shapes. In order to facilitate future MC simulations as well as density functional theory applied to freely reorienting spheroellipsoids, we have derived a fast overlap detection procedure. The method is based on Perram and Wertheim’s overlap criterion for ellipsoids, which was generalized to spheroellipsoids and accelerated by replacing an extremal search by a Newton-Raphson procedure.
1 B.
Mulder, Phys. Rev. A 35(7), 3095 (1987). Bolhuis and D. Frenkel, J. Chem. Phys. 106, 666 (1997). 3 A. Donev, F. H. Stillinger, P. M. Chaikin, and S. Torquato, Phys. Rev. Lett. 92, 255506 (2004). 4 P. Pfleiderer and T. Schilling, Phys. Rev. E 75, 020402(R) (2007). 5 G. Odriozola, J. Chem. Phys. 136(13), 134505 (2012). 6 G. Bautista-Carbajal, A. Moncho-Jordá, and G. Odriozola, J. Chem. Phys. 138(6), 064501 (2013). 7 A. Poniewierski and R. Holyst, Phys. Rev. Lett. 61(21), 2461 (1988). 2 P.
J. Chem. Phys. 141, 064109 (2014) 8 C.
Vega, C. McBride, and L. G. MacDowell, J. Chem. Phys. 115(9), 4203 (2001). 9 G. T. Evans, Mol. Phys. 76, 1359 (1992). 10 E. de Miguel and E. Martin del Río, J. Chem. Phys. 115, 9072 (2001). 11 Y. Martínez-Raton and E. Velasco, J. Chem. Phys. 129, 054907 (2008). 12 H. Xu, H. N. W. Lekkerkerker, and M. Baus, Europhys. Lett. 17(2), 163 (1992). 13 T. Koda and H. Kimura, J. Phys. Soc. Jpn. 63, 984 (1994). 14 T. Koda, M. Numajiri, and S. Ikeda, J. Phys. Soc. Jpn. 65, 3551 (1996). 15 F. J. Vesely, Mol. Phys. 103(5), 679 (2005). 16 Sz. Varga, E. Velasco, and F. J. Vesely, Mol. Phys. 107(23), 2481 (2009). 17 Sz. Varga, A. Gabor, E. Velasco, L. Mederos, and F. J. Vesely, Mol. Phys. 106(15), 1939 (2008). 18 Z. Dogic, D. Frenkel, and S. Fraden, Phys. Rev. E 62, 3925 (2000). 19 J. D. Parsons, Phys. Rev. A 19, 1225 (1979). 20 S.-D. Lee, J. Chem. Phys. 87(8), 4972 (1987). 21 A. Stroobants, H. N. W. Lekkerkerker, and D. Frenkel, Phys. Rev. Lett. 57(12), 1452 (1986). 22 J. A. C. Veerman and D. Frenkel, Phys. Rev. A 43(8), 4334 (1991). 23 D. Costa, F. Saija, and P. V. Giaquinta, J. Phys. Chem. B 107, 9514 (2003). 24 J. A. Capitán, Y. Martínez-Ratón, and J. A. Cuesta, J. Chem. Phys. 128, 194901 (2008). 25 M. Torikai, J. Stat. Phys. 148, 345–352 (2012). 26 J. W. Perram and M. S. Wertheim, J. Comput. Phys. 58, 409 (1985). 27 To locate F max , Perram and Wertheim proposed the use of a Brent-type extremal search which converges with order 1.3247. Donev et al.28 have suggested to use a Newton-Raphson search for the zero of dF(λ)/dλ. NR converges quadratically, which makes their approach as fast as ours. It should be mentioned that generally it is not recommended to replace a minimum search by a search for the zero of the derivative, but here the functions involved seem to be very well-behaved, and Donev also recommends to use a “safeguarded”NR scheme, meaning that in case of convergence problems the code would automatically switch to the slow but safe interval bisection. In our scheme the location of Fmax is not needed. 28 A. Donev, S. Torquato, and F. H. Stillinger, J. Comput. Phys. 202(2), 737 (2005).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.53.24.2 On: Fri, 03 Oct 2014 05:08:26
