Phase behaviour of hard cylinders
Joyce T. Lopes, Flavio Romano, Eric Grelet, Luis F.M. Franco, Achille Giacometti
PPhase behaviour of hard cylinders
Joyce T. Lopes, a) Flavio Romano,
2, 3, b) Eric Grelet, c) Lu´ıs F. M. Franco, d) and Achille Giacometti
2, 3, e) Universidade Estadual de Campinas, Faculdade de Engenharia Qu´ımica, Departamento de Engenharia de SistemasQu´ımicos Dipartimento di Scienze Molecolari e Nanosistemi, Universit`a Ca’ Foscari di Venezia Campus Scientifico,Edificio Alfa, via Torino 155, 30170 Venezia Mestre, Italy European Centre for Living Technology (ECLT) Ca’ Bottacin, 3911 Dorsoduro Calle Crosera, 30123 Venice,Italy Universit´e de Bordeaux, CNRS, Centre de Recherche Paul-Pascal, 115 Avenue Schweitzer, 33600 Pessac,France (Dated: 23 February 2021)
Using isobaric Monte Carlo simulations, we map out the entire phase diagram of a system of hard cylindricalparticles of length L and diameter D , using an improved algorithm to identify the overlap condition betweentwo cylinders. Both the prolate L/D >
L/D <
L/D >
L/D <
I. INTRODUCTION
After nearly one century since Onsager’s pioneeringprediction that orientational order can be entropicallyinduced for elongated particles , simple models of rod-like objects continue to play a central role in the study ofcolloidal liquid crystals and self-assembly processes .The simplest model for a rod-like molecule is the hardspherocylinder, an object formed by a cylinder of length L capped with two hemispheres of matching diameter D . This shape can be obtained by rolling a sphere ofradius D/ L . The greatadvantage of this model, and the key to its popularity,is the simplicity of the overlap condition between twosuch hard spherocylinders; this condition can be cast ina simple analytical form that can be computed veryefficiently. As early as 1997, Bolhuis and Frenkel per-formed a remarkably detailed study of the phase diagramof this model that is now reckoned as a classic referencein the field. Other similar shapes have also been pro-posed in the literature, including hard ellipsoids , hardhelices , and hard dumbbells .However, there are physically relevant objects whoseshape cannot be represented as hard spherocylinders butrather as hard cylinders. Examples include biologicallyrelevant cases such as viruses and nucleosomes .Hard cylinders of length L and diameter D have alsothe additional advantage of having a natural oblate limit a) Electronic mail: [email protected] b) Electronic mail: fl[email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected] e) Electronic mail: [email protected]
L/D <
1, approaching a disk for
L/D →
0, as well as theprolate limit
L/D > . By contrast theoverlap condition between two cylinders is significantlymore evolved with respect to the spherocylinder case.This notwithstanding, and given the similarity in shape,one might rightfully wonder what are the differences, ifany, in the two phase diagrams. For instance, the phasediagram of hard ellipsoids is different from the phase di-agram of hard spherocylinders, in spite of the significantsimilarities in their shapes. This issue goes far beyond asimple academic problem in view of the strong propensityof nucleosomes and filamentous viruses to form acolumnar phase, whose existence in the phase diagram ofspherocylinders has been ruled out by recent detailed nu-merical simulations for prolate particles even if it werepredicted theoretically for oblate ones.The aim of the present paper is to tackle this issueby performing a detailed analysis of the phase diagramof hard cylinders both in the prolate ( L/D >
1) and inthe oblate (
L/D <
1) cases. While simulations of hardcylinders have been performed in the past , to thebest of our knowledge, our study is the first one provid-ing the complete phase diagram. For this purpose, weperform isobaric Monte Carlo simulations of a system ofhard cylinders in a wide range of aspect ratios
L/D andvolume fractions, using an efficient method for the over-lap test that compares well with existing ones . Thealgorithm, inspired by Ref. 22, is described in the Ap-pendix. By monitoring the appropriate order parametersand correlation functions, we provide the correspondingphase diagram in the
L/D volume fraction plane andcompare it with the corresponding phase diagram of thehard spherocylinders . Particular care has been devoted a r X i v : . [ c ond - m a t . s o f t ] F e b in avoiding possible finite size effects along the lines of arecent similar analysis for spherocylinders .The outline of the paper is as follows. In Section II wedescribed the details of our numerical approach, as wellas the arsenal of tools (order parameters and correlationfunctions) useful to identify all different phases. Sometechnical details have been confined in Appendix A. Sec-tion III will report the main results of the present studywith additional Figures and Tables reported in SectionB. Finally, Section IV reports the key messages of thisstudy and some interesting perspectives for the future. II. SIMULATIONSA. Monte Carlo simulations
Our particles consist of N cylinders/disks of height L ,diameter D , and whose orientations are identified by aunit vector (cid:98) u , as shown in Figure 1 (a). Pressures aremeasured in reduced units P ∗ = P v HC /k B T , and thedensity ρ = N/V is represented by the volume fraction η ≡ N v HC /V , where v HC = LπD /4 is the volume of ahard cylinder (HC). We then performed isobaric (NPT)Monte Carlo (MC) simulations at different aspect ratios L/D both for rods (
L/D >
1) and for disks (
L/D < ≈ . × MCsteps, with additional production runs of 1 . × steps.The typical number of particle was N ≈ N was adjusted depending on L/D tokeep the simulation box roughly cubic. In our NPT sim-ulations, we have used floppy (i.e. shape-adapting) rect-angular computational box, where one axis was randomlyselected and its length was allowed to change, with peri-odic boundary conditions to obtain an isotropic pressure in the prolate L/D >
L/D < . B. Overlap of hard cylinders
The first method for testing overlaps of hard cylinderswas proposed by Allen et al. , and also used by Blaak,Frenkel, and Mulder . An alternative method was re-cently proposed by Orellana, Romani, and De Michele .In the following, we use a refined version of this method outlined below. The overlap of two cylinders can occurin either of the following three ways: disk-rim, rim-rim,and disk-disk (Figure 1). Therefore, to ensure that thecylinders do not overlap, we have to check whether theoverlap occurs in one of those possible configurations, asdetailed in Appendix A.Preliminary simulations were initially performed to as-sess the computational effort of this algorithm comparedwith the hard spherocylinders counterpart. We found thepresent algorithm to be slightly slower - of the order of20% or less, and hence in line with Ref. 22 (a) (b)(c) (d) FIG. 1: (a) Our cylinder model, where L is the height, D is the diameter, (cid:98) u is the unit vector defining theorientation of the cylinder; possible overlapconfigurations between two cylinders (see appendix):(b) rim-rim; (c) rim-disk; (d) disk-disk.We perform the less expensive test first, and progres-sively include additional more expensive ones. So wefirst check whether the two spheres (of diameter L + D )that encompass the cylinders overlap; if they do not, thecylinders cannot overlap. If the encompassing spheresdo overlap, the test is repeated for the spherocylindersenclosing the particles, using the standard algorithm tocalculate the shortest distance between two rods . Onlyif the spherocylinders overlap, the overlap between twocylinders is tested for. See Appendix A for additionaldetails. C. Order parameters
To identify different thermodynamic phases, we relyon information based on global orientational and trans-lational order, i.e., the nematic, smectic and hexatic or-der parameters, on correlation functions such as the ra-dial g ( r ), parallel g (cid:107) ( r (cid:107) ) and perpendicular g ⊥ ( r ⊥ ) dis-tribution functions, as well as on visual inspection ofthe simulation snapshots. Figure 2 displays represen-tative snapshots of all different phases obtained in the L/D = 10 case – all snapshots were obtained with theOvito Software where different colours represent differ-ent orientations of the cylinders. While the isotropic Iphase is both positionally and orientationally disorderedthe nematic N phase is positionally disordered but orien-tationally ordered, and its presence can be inferred mon-itoring the nematic order parameter P . This is obtainedas the largest eigenvalue of the tensor Q αβ = 1 N N (cid:88) i =1 (cid:98) u iα (cid:98) u iβ − δ αβ (1)where α, β = x, y, z . The corresponding eigenvector thengives the main director (cid:98) n .In addition to the orientational order along one pre-ferred direction (cid:98) n , the smectic phase SmA is further char-acterised by a one-dimensional ordering (layering) along (cid:98) n that is best captured by a combination of the radialdistribution function g ( r ) = 1 N ρ πr (cid:42) N (cid:88) i =1 N (cid:88) j (cid:54) = i δ ( r − r ij ) (cid:43) (2)as well as the parallel g (cid:107) ( r (cid:107) ) = 1 N (cid:42) ρL x L y N (cid:88) i N (cid:88) j (cid:54) = i δ ( r (cid:107) − r ij · (cid:98) n ) (cid:43) (3)positional correlation function. Here r i is the center ofmass of the i -th cylinder, and r ij = r j − r i , and r ij = | r ij | .The smectic order parameter (cid:104) τ (cid:105) = (cid:12)(cid:12)(cid:12)(cid:68) e i2 π r · (cid:98) n d (cid:69)(cid:12)(cid:12)(cid:12) (4)also proves convenient. Here r is the position of a parti-cle’s centre of mass and d the optimal layer spacing. Hereand below, (cid:104) . . . (cid:105) is the average over independent config-urations at equilibrium. Then we have (cid:104) τ (cid:105) ≈ (cid:104) τ (cid:105) ≈ (cid:98) n . This is best cap- tured by the perpendicular positional correlation func-tion g ⊥ ( r ⊥ ) = 12 πr ⊥ N (cid:42) ρL z N (cid:88) i N (cid:88) j (cid:54) = i δ ( r ⊥ − | r ij × (cid:98) n | ) (cid:43) (5)positional correlation functions, as well as by the use ofthe hexatic (or bond) order parameter (cid:104) ψ (cid:105) = (cid:42) N (cid:88) j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n ( j ) (cid:88) (cid:104) lm (cid:105) e θ lm (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:43) . (6)Here θ lm is the angle that the projection of the inter-molecular vectors r jl and r jm onto the plane perpendic-ular to the director (cid:98) n , n ( j ) is the number of nearest-neighbours pairs of molecule within a single layer, andthe sum (cid:80) (cid:104) lm (cid:105) is over all possible pairs within the firstcoordination shell. With this definition, (cid:104) ψ (cid:105) ≈ (cid:104) ψ (cid:105) ≈ et al. and refer-ences therein) for additional details.Finally, the cubatic Cub phase corresponds to a long-range orientationally ordered phase without any posi-tional order but with the presence of three equivalent per-pendicular directions. In this phase, the particles formshort stacks of typically few particles with neighbouringstacks tending to be perpendicular to one another alongthe three selected directions. While a suitable order pa-rameter can be devised , visual inspection is usually suf-ficient to unambiguously identify this phase. The detailsof the cubatic Cub phase will be discussed in Section III.TABLE I: Number of particles N used in thesimulations of rods. L/D N L/D N
III. RESULTSA. Cylindrical rods
L/D > We first consider the prolate case, i.e. cylindrical rodswith
L/D >
1. Figure 3 (a) depicts the reduced pres-sure P ∗ as a function of the volume fraction η (i.e., theequation of state) for L/D = 5 and
L/D = 10. Figure3 (b) shows also the corresponding orientational orderparameter P again as a function of the volume fraction η .FIG. 2: Representative snapshots of the thermodynamic phases found for HC for L/D = 10: Isotropic (I), Nematic(N), Smectic A (SmA), and Crystal (X). Reduced corresponding pressures P ∗ are displayed.In the case L/D = 5 (open symbols) the system is inan isotropic phase I until η ≈ .
4, then switches to asmectic SmA phase, and then to a crystal X. The samesequence of phases is also found for the large aspect ratio
L/D = 10 (closed symbols) but with transitions shiftedto lower η , and with the additional presence of a nematicN phase in the region 0 . ≤ η ≤ .
4. The isotropic-nematic transition is signalled by an abrupt jump in thenematic order parameter P and by a discontinuity in theequation of state as shown in Figure 3 (a).Here it is worth to notice that our definition of crys-tal phase X includes the so-called smectic SmB phase,another name often used in this framework , that is asmectic SmA phase with additional in-plane long-rangehexagonal order , thus hardly distinguishable from acrystal phase due to the finite size of the simulations box.Additional insights can be obtained by looking atthe correlation functions at an intermediate aspect ra-tio L/D = 7 – see Figure 14 for the analogue of Figure 3in the case
L/D = 7. Figure 4 presents the correspond-ing radial g ( r ) (a), parallel g (cid:107) ( r (cid:107) ) (b), and perpendic-ular g ⊥ ( r ⊥ ) (c) distribution functions of cylinders with L/D = 7 for increasing pressures.At P ∗ = 3 .
96 (continuous red line) all correlationfunctions are featureless, indicating the presence of anisotropic I phase. As pressure is increased up to P ∗ =4 .
40 (yellow dashed line), the correlation functions donot show any significant change but the nematic orderparameter P (see Figure 14 (b)) shows an abrupt up-swing, signalling the onset of a nematic N phase.At P ∗ = 7 .
70 (green dotted line), both the radial distri-bution function g ( r ) and the parallel correlation function g (cid:107) ( r (cid:107) ) display a clear periodicity consistent with a smec-tic SmA ordering. The absence of regular oscillations inthe perpendicular correlation function g ⊥ ( r ⊥ ) confirmsthe radial liquid-like order of the mesophase, which there-fore does not correspond to the crystal X phase. The lat-ter phase is eventually reached at P ∗ = 9 .
90 (dash-dottedline) as shown by the characteristic periodicities for all directions in the g ⊥ ( r ⊥ ) as well as in g ( r ) and g (cid:107) ( r (cid:107) ).It comes as no surprise that the low- η behaviour ofHard Cylinders is qualitatively similar to the correspond-ing HSC counterpart , with small quantitative differ-ences for the smaller aspect ratio L/D = 5. However,at high pressure and volume fraction, one possible im-portant element of distinction between the two phasediagrams is the presence of a putative columnar phasethat has already been demonstrated not to exist in theHSC counterpart . We explicitly addressed this prob-lem following the method proposed by Dussi, Chiappini,and Dijkstra who suggested that the apparent stabili-sation of a columnar phase in HSCs could be ascribed tofinite size effects when the number of layers is not suf-ficiently high compared to the aspect ratio L/D . Forsake of consistency, we first reproduced the same resultsfound in Ref. 19 for HSCs, and then applied the samemethod to the HC case. We note that the metastabilityof the columnar phase for HSCs was also independentlyconfirmed by Liu and Widmer-Cooper using a differentmethod.The results obtained are presented in Figure 5. Figure5 (a) shows a production run with aspect ratio L/D = 6at packing fraction η = 0 .
6. In this case, both visualinspection and the behaviour of the corresponding cor-relation functions (see solid line in Figure 15 for resultswith N = 675) strongly suggest the presence of a colum-nar phase. However, if the number of particles is doubledalong the director (cid:98) n , the same calculation produces thefinal configuration shown in Figure 5 (b) that can clearlybe classified as smectic SmA (see dotted line in Figure 15for results with N = 1350). This shows that there is nostable columnar phase in HCs as in the HSCs case. Thiseffect is likely to be ascribed to the preference for finitesize domains to arrange locally in columnar structureswhose stability is eventually overwhelmed by long-rangeeffects.A sketch of the final phase diagram for HCs in theplane packing fraction η as a function of the aspect ra- P * η I SmA X (a) P η IN SmA X (b)
FIG. 3: (a) Reduced pressure P ∗ versus cylinder volumefraction η . Open symbols: L/D = 5, closed symbols:
L/D = 10; (b) Nematic order parameter P versusvolume fraction η for both L/D = 10 and
L/D = 5.Same symbols as above. The different symbols andcolours refer to different mesophases, as detailed inFig 2.tio ranging from
L/D = 2 . .Similar to hard spherocylinders, the system exhibitsthe isotropic (I), nematic (N), smectic A (SmA), and g (r) r/D P* = 3.96P* = 4.40P* = 7.70P* = 9.90 (a) g || (r) r/D P* = 3.96P* = 4.40P* = 7.70P* = 9.90 (b) g ⊥ (r) r/D P* = 3.96P* = 4.40P* = 7.70P* = 9.90 (c)
FIG. 4: Distribution functions of cylinders with
L/D = 7 . P ∗ = 3 .
96 continuous red line (I); P ∗ = 4 . P ∗ = 7 .
70 green dotted line(SmA); P ∗ = 9 .
90 dash-dotted line (X). Note that r = | r | in (a), r = | r (cid:107) | in (b), and r = | r ⊥ | in (c).crystalline (X) phases. Not surprisingly, this behaviouris similar to that of HSC but few differences are worthnoticing. (a) Unstable columnar phase(b) Crystalline phase FIG. 5: Equilibrated configuration with aspect ratio
L/D = 6 at packing fraction η = 0 . N = 675initially distributed on two layers;(b) The same resultwith twice the particles exhibits a different structure.As in the case of HSC, no liquid crystal phases are ob-served below a critical aspect ratio L/D ≈
3. This factcan be easily rationalised via Onsager theory , as the ra-tio between the covolume and volume of rods with lower L/D are not sufficiently larger than that of a sphere,and the excluded volume effects then are insufficient topromote an organised orientationally ordered phase. Bycontrast, at sufficiently high densities and aspect ratios,exclude volume effects tend to promote orientational or-der to increase the translational entropy, then minimisingthe free energy.Accordingly, the system is in isotropic phase for anyratio
L/D below a certain packing fraction that decreasesby increasing the rod aspect ratio, as shown in Figure 6. FIG. 6: Computed phase diagram of hard cylinders ofpacking fraction η versus aspect ratio L/D . Visiblephases are isotropic (I), nematic (N), smectic (SmA),and crystal (X). Colour codes are as in Figure 2.Upon increasing η , the first organised phase encounteredis a smectic SmA phase in the range from L/D = 3 .
25 to
L/D = 6, and a nematic N phase above
L/D ≈
6. Thismirrors the HSC case where, however, the smectic SmAphase is limited to a very small range 3 < L/D < η , the system undergoes asmectic SmA to crystal X transition irrespective of theaspect ratio L/D .The sketched phase diagram in Figure 6 prompts theexistence of an isotropic-smectic-solid (I-SmA-X) triplepoint at η ≈ .
55 and
L/D ≈ .
0, and an isotropic-nematic-smectic (I-N-SmA) triple point at η ≈ . L/D = 6 . L/D ≈ . L/D ≈ . . As a result,the nematic N phase stabilises at shorter aspect ratios inthe HSC system when compared to its HC counterpart.It is also interesting to notice that our results are com-patible with both the I-N and the N-SmA being first-order transitions, thus mirroring what is known for HSCfrom the work by Polson and Frenkel . It would be in-teresting to pursue the same analysis carried out by theseauthors in the present case as well. The same considera-tion holds true for interesting analysis of the L/D → ∞
Onsager limit that has been performed for the HSC case that could be also replicated in this case.As a final remark we note that the length L for HSCscorresponds to L + D in case of HCs. This is impor-tant when comparing the corresponding phase diagramsand indeed it rationalizes why the isotropic-smectic I-SmA transition for HCs occurs at L/D ≈ L/D ≈
4. However, the tendency of aflat edge to promote the onset of a smectic SmA phaseappears to be a general feature as also suggested by a re-cent study on hard equilateral triangular prisms, wherethe particles feature also flat sides and the smectic SmAphase is shifted to considerably lower packing fractionsas compared to HSCs. B. Cylindrical disks
L/D < We now tackle the oblate case of cylindrical disks with
L/D <
1. One important advantage of dealing withcylinders is that this limit can be achieved with no solu-tion of continuity, unlike the spherocylinders counterpartwhere this is not possible . Figure 7 depicts the fourdifferent phases that we find in this case: a disorderedisotropic I, a cubatic Cub, a nematic N, and a columnarC phases, as detailed in Table II and illustrated in Figure7.TABLE II: Colours and symbols used to represent diskphases. Colour Phase Notation Symbolred isotropic I circleyellow nematic N trianglepurple cubatic Cub squaresblue columnar C diamond
TABLE III: Number of particles N used in thesimulations of disks. L/D N L/D N
As in this case of prolate cylinders, we performed thesame detailed analysis of the different obtained phasesin terms of correlation functions and order parametersto derive the equation of states. Supplementary FigureS3 reports the reduced pressure P ∗ and the P nematicorder parameter as a function of the packing fraction η for both L/D = 0 . L/D = 0 .
05 as representativeexamples, from which one can obtain the correspond-ing phase diagram of Figure 8 in the volume fraction η aspect ratio L/D plane, that can be contrasted withits prolate counterpart of Figure 6. Here a range from
L/D = 0 .
05 to
L/D = 0 . L/D > . < L/D < . η above ≈ .
4. In the columnar phase thedisks are arranged on a hexagonal lattice in the directionperpendicular to the main director (cid:98) n but their centresof mass are disorderly distributed in space. For smalleraspect ratios 0 . < L/D < .
3, a cubatic Cub phase ap-pears between the I and C phase. In the cubatic phase,the disks tend to assemble in short stacks of about fouror five units, with neighbouring columns perpendicularto each other. This differs from the cubic phase becauseit lacks translational order . At even smaller aspect ra-tio ( L/D ≤ .
1) the cubatic Cub phase is replaced by anematic N phase up to η ≈ . η . At even higher packing fractions η we didobserve the formation of a crystal phase X but the loca-tion of corresponding boundaries would require a specificinvestigation that was not pursued in the present paper.All these transitions can be best inferred by lookingat the correlation functions as shown in Figure 9 (seeFig. 7 for the corresponding snapshots). In the I phase(red continuous line) the radial distribution function g ( r )displays a flat behaviour for r > D , indicating the ab-sence of short-range aggregation, a feature confirmed bythe behaviour of both g (cid:107) ( r (cid:107) ) and g ⊥ ( r ⊥ ). By contrast,the columnar C phase (blue dashed line) displays char-acteristic regular oscillations in g ( r ) and g ⊥ ( r ⊥ ), but thebehaviour of g (cid:107) ( r (cid:107) ) is irregular, indicating the absenceof a one-dimensional ordering along the main director (cid:98) n .Likewise, the radial distribution function g ( r ) of the cu-batic phase (dotted purple line) is quite different fromboth the I and C phases, while the nematic order param-eter P is close to zero for I, and Cub. An evidence ofthe formation of short stacks is the higher peak at shortdistances ( L/D < r/D < L/D ) in the radial distribu-tion function g ( r ) of the cubatic Cub phase (purple linein Figure 9), which is significantly smaller in the the g ( r )of an isotropic phase (red line in Figure 9). Finally, theonset of the nematic N phase (dash-dotted yellow line)is signalled by the significant oscillation of the radial dis-tribution function g ( r ) and by the abrupt upswing in thenematic order parameter P , as shown in supplementarymaterials.At variance with the prolate L/D >
L/D ≈ . η ≈ .
3. The N-Cub-C triple point is approximately located at
L/D ≈ . η ≈ .
4. Finally, the I-Cub-C triple point hasapproximate coordinates
L/D ≈ .
35 and η ≈ . who pre-dicted the existence of a nematic region for flat disksthat becomes progressively narrower as L/D increases.The same authors also predicted a transition from theisotropic phase directly to the columnar C phase, inagreement with our results. At a more quantitativelevel the predicted volume fractions η IN ≈ πL/D forthe isotropic-nematic transition and η NC ≈ . L/D <
L/D and reduced pressure P ∗ arealso reported.At variance of these theoretical findings, our results in-dicate also the existence of a cubatic Cub phase, in agree-ment with the results by Veerman and Frenkel , as wellas by Duncan et al. , in simulations of cut spheres, andby Blaak, Frenkel, and Mulder in simulations of hardcylinders. Where direct comparison with the above twopapers is possible, we find complete agreement betweentheir results and ours.FIG. 8: Phase diagram of the oblate hard cylindricaldisk case L/D < η aspect ratio L/D plane. Different phases are colour coded asdetailed in Figure 7 and Table II.Duncan et al. simulated cut spheres of L/D = 0 . .
15, 0 .
2, 0 .
25 and 0 . L/D = 0 .
1, and that the opposite is true for
L/D ≥ .
15. Figure 8 shows that for HCs a cubatic Cubphase is already present at
L/D = 0 .
11, as the nematic N phase vanishes. The cubatic Cub phase is present until
L/D ≈ . investi-gated a system of HC with L/D = 0 . L/D = 0 . . < L/D < .
3, as shown in Figure 8.
IV. CONCLUSIONS
In this paper we have used isobaric (
N P T ) MonteCarlo simulations to study the phase diagram of a systemof N hard cylinders as a function of their aspect ratio andvolume fraction. To achieve this, we have implementeda new and efficient overlap test for hard cylinders thatcompares well with those existing in the literature .This allows us to study the complete phase diagrams inthe aspect ratio versus volume fraction for both the pro-late L/D >
L/D <
L/D >
1, we find a phase diagramvery similar to the hard spherocylinders counterpart, fea-turing the presence of a nematic N and a smectic SmAphases, in addition to the isotropic I and the crystal Xphases, as well as two I-Sm-X and I-N-SmA triple points.As in the spherocylinder case , we have shown that theappearance of a columnar C phase can be traced back toa finite-size effect and that it disappears for sufficientlylarge systems. Our simulations confirm the lack of exis-tence of a stable columnar C phase, that was neverthelesspredicted by density functional theory .In the oblate case L/D <
1, we identified the pres-ence of a columnar C, a nematic N and a cubatic Cubphase, in agreement with theoretical prediction , as wellas with past numerical simulations of cut spheres andof hard cylinders . In the latter case, we have providedan explanation of the failure of past simulations of iden- g (r) r/D L/D = 0.1, P* = 1.18L/D = 0.05, P* = 1.37L/D = 0.2, P* = 5.50L/D = 0.5, P* = 8.25 (a) g || (r) r/D L/D = 0.1, P* = 1.18L/D = 0.05, P* = 1.37L/D = 0.2, P* = 5.50L/D = 0.5, P* = 8.25 (b) g ⊥ (r) r/D L/D = 0.1, P* = 1.18L/D = 0.05, P* = 1.37L/D = 0.2, P* = 5.50L/D = 0.5, P* = 8.25 (c)
FIG. 9: Distribution functions of hard cylindrical disks.
L/D = 0 . P ∗ = 1 .
18 red continuous line I;
L/D = 0 .
05 and P ∗ = 1 .
37 blue dashed line C;
L/D = 0 . P ∗ = 5 .
50 dotted purple line Cu;
L/D = 0 . P ∗ = 8 .
25 dash-dotted yellow line N.Colour code is outlined in Table II.Here again r = | r | in(a), r = | r (cid:107) | in (b), and r = | r ⊥ | in (c).tifying the cubatic Cub phase that can be ascribed to the too large aspect ratio used in these simulations. Inter-estingly, the phase diagram also includes three I-N-Cub,N-Cub-C, and I-Cub-C triple points.The present study paves the way to tackling more com-plex systems building upon cylindrical shapes that are ofexperimental interest, such as hard cylinders interactingvia a Yukawa tail , as well as hard cylinders with short-range directional attractions . Such investigations areunderway and will be reported elsewhere. ACKNOWLEDGMENTS
We are indebted to Cristiano De Michele and MariaBarbi for useful discussions. The use of the SCSCF mul-tiprocessor cluster at the Universit`a Ca’ Foscari Veneziaand of the LESC cluster at the University of Campinasare gratefully acknowledged. The work was supportedby MIUR PRIN-COFIN2017
Soft Adaptive Networks grant 2017Z55KCW and Galileo Project 2018-39566PG(AG), by S˜ao Paulo Research Foundation (FAPESP)grant no. 2018/02713-8, and by the Coordena¸c˜ao deAperfei¸coamento de Pessoal de N´ıvel Superior - Brasil(CAPES) - Finance Code 001. FR acknowledges IdExBordeaux (France) for financial support. JL gratefullyacknowledges the hospitality of the Ca’ Foscari Univer-sity of Venice where part of this work was carried out.The authors would like to acknowledge the contributionof the Eutopia COST Action CA17139.
DATA AVAILABILITY
The data that supports the findings of this study areavailable within the article and its supplementary mate-rial.
Appendix A: Algorithm to check overlap between twocylindersa. Parallel Cylinders
If two cylinders are parallel, the overlap can occur be-tween disk-disk or rim-rim only, and it can be easilychecked. For each particle pair we define four vectorsstarting from the the vector joining the two centres ofmass, r . We do this by extracting its parallel and per-pendicular component with respect to the director (cid:98) u j ofeach particle, i.e., r (cid:107) = ( r ij · (cid:98) u ) (cid:98) u r (cid:107) = ( r ij · (cid:98) u ) (cid:98) u r ⊥ = r − ( r · (cid:98) u ) (cid:98) u r ⊥ = r − ( r · (cid:98) u ) (cid:98) u (A1)In a parallel configuration, the directors can either bethe same of one the opposite of the other. The overlap0occurs if all the following conditions are satisfied: | r (cid:107) | ≤ L | r (cid:107) | ≤ L | r ⊥ | ≤ D | r ⊥ | ≤ D (A2)When the cylinders are exactly parallel, (cid:98) u = ± (cid:98) u ,and half of the conditions above are redundant since | r ⊥ | = | r ⊥ | and | r (cid:107) | = | r (cid:107) | ; when implementing com-puter code, however, one has to include tolerances andcare must be taken in in handling these conditions con-sistently. b. Rim-rim overlap Since the overlap between spherocylinders is the firsttest that is done, and the rim of a spherocylinder is sim-ilar to the rim of a cylinder, if the rims of two sphero-cylinders do overlap, than the two cylinders will certainlyoverlap as well. Hence, having performed the sphero-cylinder overlap test, we now check if the overlap occursin a rim-rim configuration.To that end, we define the vectors V = − r + λ (cid:98) u and V = r + µ (cid:98) u , where the numbers λ and µ , consis-tently with Ref. 7, identify the points of closest approachbetween the axes of the two cylinders. These values arecalculated using the Vega and Lago’s algorithm , whichwe implement in the spherocylinder overlap test. If thecylinders are in a rim-rim configuration, the two condi-tions below must both be satisfied: | V · (cid:98) u | The orientations of the cylinders are perpendicular tothe planes of the disks. The planes of the two disks in-tersect in a line parallel to (cid:98) u × (cid:98) u . We define P and P as being the points in the intersection line that are closerto the disks centers d and d , respectively, as shown inFigure 11. FIG. 11: Disks of two cylinders.To find P , we minimise ( P − d ) , which is equivalentto minimising | P − d | . The minimisation can be doneby applying Lagrange multipliers with two constraints:( P − d ) · (cid:98) u = 0 (A4a)( P − d ) · (cid:98) u = 0 (A4b)The constraints presented in Equation A4 ensure that P is in a line perpendicular to both (cid:98) u and (cid:98) u . Applyingthe Lagrange multipliers: L = ( P − d ) − λ ( P − d ) · (cid:98) u − µ ( P − d ) · (cid:98) u (A5)From ∇L = 0, one has: P = d + λ (cid:98) u µ (cid:98) u λ = − µ ( (cid:98) u · (cid:98) u ) (A7)Substituting Equations A4b and A7 into A6 yields: µ = − d − d ) · (cid:98) u − ( (cid:98) u · (cid:98) u )) (A8)Replacing Equation A8 into A7: λ = 2 [( d − d ) · (cid:98) u ] · ( (cid:98) u · (cid:98) u )1 − ( (cid:98) u · (cid:98) u ) (A9)Replacing Equations A9 and A8 into A6: P = d + [( d − d ) · (cid:98) u ] · (( (cid:98) u · (cid:98) u ) · (cid:98) u − (cid:98) u )1 − ( (cid:98) u · (cid:98) u ) (A10)1We define d = d − d and ∆ = ( P − d ) , andrewrite Equation A10 as:∆ = ( d · (cid:98) u ) · (( (cid:98) u · (cid:98) u ) − (cid:98) u · (cid:98) u ) + 1)(1 − ( (cid:98) u · (cid:98) u ) ) (A11)Simplifying Equation A11:∆ = ( d · (cid:98) u ) − ( (cid:98) u · (cid:98) u ) (A12)Similarly for disk 2:∆ = ( d · (cid:98) u ) − ( (cid:98) u · (cid:98) u ) (A13)A necessary, but not sufficient, condition for the over-lap to occur is that both ∆ and ∆ have to be less thanthe cylinder radius D/ 2. If this condition is satisfied, theintersection line crosses both disks through segments oflength 2 δ and 2 δ , as presented in Figure 12.FIG. 12: Disks of two cylinders.The expressions to calculate δ and δ are presented inEquation A14. δ = (cid:114) D − ∆ δ = (cid:114) D − ∆ (A14)Finally, an overlap will occur if the condition in thefollowing equation is true: | P − P | = (cid:12)(cid:12)(cid:12)(cid:12) d · ( (cid:98) u × (cid:98) u ) | (cid:98) u × (cid:98) u | (cid:12)(cid:12)(cid:12)(cid:12) ≤ δ + δ (A15) d. Disk-rim overlap Let us take a disk with centre in d j and a cylinder withcentre in r i . We define U i as the point on cylinder i thatis the closest to d j , P d a point on the disk j that is theclosest to cylinder i , P c a point on cylinder i that is theclosest to disk j , φ an angle between (cid:98) w j and d j − P d , (cid:98) v j , (cid:98) u j , an axis system fixed on cylinder j and, finally, φ as an angle between (cid:98) w j and d j − P d . FIG. 13: Disk-rim configuration. U i is obtained from: U i = r i + [( d j − r i ) · (cid:98) u i ] (cid:98) u i (A16)First, we test the following conditions:1. If | d j − U i | > d : there is no overlap2. If | d j − U i | < d/ | d j − r i | > L/ | d j − U i | ≤ d/ | ( d j − r i ) | < L/ j is within cylinder i .Test number 3 is a sufficient, but not necessary, condi-tion for the overlap to occur, since another point can betouching cylinder j even if d j is not within cylinder i .Hence, if condition 3 is not satisfied, we have to find P d , the closest point in disk j to cylinder i .Arbitrary points on the border of disk j ( d ), and onthe line of cylinder i ( c ) are defined as: d = d j + R cos ( φ ) (cid:98) w j + R sin ( φ ) (cid:98) v j (A17a) c = r i + λ (cid:98) u i (A17b)where R ≡ D/ d and c is thus:( d − c ) = ( d j − r ) + R + λ (A18)+ 2 R cos φ (( d j − r i ) · (cid:98) w j )+ 2 R sin φ (( d j r i ) · (cid:98) v j ) − λ (( d j − r i ) · (cid:98) u i ) − λR cos φ ( (cid:98) w j · (cid:98) u i ) − λR sin φ ( (cid:98) v j · (cid:98) u i ) P c and P d are the points that minimise Equation A18,therefore: λ − r cos φ ( (cid:98) w j · (cid:98) u i ) − r sin φ ( (cid:98) v j · (cid:98) u i ) + (A19) − (( d j − r i ) · (cid:98) u i ) = 02sin φ [ λ ( (cid:98) w j · (cid:98) u i ) − (( d j − r i ) · (cid:98) w j )] + (A20) − cos φ [ λ ( (cid:98) v j · (cid:98) u i ) − (( d j − r i ) · (cid:98) v j )] = 0Rewriting Equation A21 gives:sin φ cos φ = λ ( (cid:98) v j · (cid:98) u i ) − (( d j − r i ) · (cid:98) v j ) λ ( (cid:98) w j · (cid:98) u i ) − (( d j − r i ) · (cid:98) w j ) (A21)If the numerator and denominator of Equation A21are taken as the catheti of a triangle, the hypotenuse canthen be found to give the expressions for cos φ and sin φ .Once we have these expressions, they are applied intoEquation A20, resulting in an equation for λ . Since wewere not able to find an analytical solution to the previ-ous equation, a numerical method such as the Newton-Raphson or bisebsection method is used to find λ . Inour code, we combine both methods, running a few stepswith one and a few with the other until convergence isfound to machine precision.Once P d is obtained, we define T = P d − r i , and cal-culate the components of T that are parallel T (cid:107) and per-pendicular T ⊥ to (cid:98) u i . T (cid:107) = ( T · (cid:98) u ) (cid:98) u (A22a) T ⊥ = T − T (cid:107) (A22b)Finally, the overlap only occurs if | T (cid:107) | ≤ L/ | T ⊥ | ≤ D/ Appendix B: Supplementary material1. Figures P * η I NSmA X (a) P η I N SmA X (b) FIG. 14: (a) Reduced pressure P ∗ versus volumefraction η for L/D = 7. (b) Nematic order parameter P versus volume fraction η for L/D = 7.3 g (r) r/D N = 1350N = 675 (a) g || (r) r/D N = 1350N = 675 (b) g ⊥ (r) r/D N = 1350N = 675 (c) FIG. 15: Distribution functions of hard cylinders withaspect ratio L/D = 6 and η = 0 . § (a) g ( r ); (b) g (cid:107) ( r (cid:107) );(c) g ⊥ ( r ⊥ ). Dotted line N = 1350 and P ∗ = 9 . 42 , Solidline N = 675 and P ∗ = 8 . P * η I NCub C (a) P η I CubN C (b) FIG. 16: (a) Reduced pressure P ∗ versus volumefraction η . Open symbols: L/D = 0 . 2, closed symbols: L/D = 0 . 05; (b) Nematic order parameter P versusvolume fraction η for both L/D = 0 . L/D = 0 . 2. Tables TABLE IV: L/D = 0 . P ∗ η S ψ τ Phase0.20 0.066 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± TABLE V: L/D = 0 . P ∗ η S ψ τ Phase0.79 0.194 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± TABLE VI: L/D = 5 . P ∗ η S ψ τ Phase0.79 0.193 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± L/D = 7 . P ∗ η S ψ τ Phase3.96 0.354 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± TABLE VIII: L/D = 10 . P ∗ η S ψ τ Phase0.79 0.165 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± L. Onsager, “The effects of shape on the interaction of colloidalparticles,” Ann. N. Y. Acad. Sci. , 627–659 (1949). T. Gibaud, E. Barry, M. Zakhary, M. Henglin, A. Ward, Y. Yang,C. Berciu, R. Oldenbourg, M. Hagan, D. Nicastro, R. Meyer,and Z. Dogic, “Self-assembly through chiral control of interfacialtension,” Nature , 348 (2012). E. Grelet, “Hard-Rod Behavior in Dense Mesophases of Semiflex-ible and Rigid Charged Viruses,” Phys. Rev. X , 021053 (2014). Y. Liang, Y. Xie, D. Chen, C. Guo, S. Hou, T. Wen, F. Yang,K. Deng, X. Wu, I. I. Smalyukh, and Q. Liu, “OSymmetry con-trol of nanorod superlattice driven by a governing force,” NatureCommun. , 1410 (2017). B. Sung, A. de la Cotte, and E. Grelet, “Chirality-controlledcrystallization via screw dislocations,” Nat. Commun. , 1405(2018). M. P. Allen, G. T. Evans, D. Frenkel, and B. M. Mulder, “Hardconvex body fluids,” Advances in chemical physics , 1–166(1993). C. Vega and S. Lago, “A fast algorithm to evaluate the shortestdistance between rods,” Comput. Chem. , 55–59 (1994). P. Bolhuis and D. Frenkel, “Tracing the phase boundaries of hardspherocylinders,” J. Chem. Phys. , 666–687 (1997). D. Frenkel and B. Mulder, “The hard ellipsoid-of-revolution fluid,” Molecular Physics , 1171–1192 (1985),https://doi.org/10.1080/00268978500101971. E. Frezza, A. Ferrarini, H. B. Kolli, A. Giacometti, and G. Cinac-chi, “The isotropic-to-nematic phase transition in hard helices:Theory and simulation,” J. Chem. Phys. , 164906 (2013). K. Milinkovi´c, M. Dennison, and M. Dijkstra, “Phase diagram ofhard asymmetric dumbbell particles,” Phys. Rev. E , 032128(2013). J. Bernal and I. Fankuchen, “X-RAY AND CRYSTALLO-GRAPHIC STUDIES OF PLANT VIRUS PREPARATIONS,”J. Gen. Physiol. , 111–165 (1941). X. Wen, R. B. Meyer, and D. L. D. Caspar, “Observation ofSmectic-A Ordering in a Solution of Rigid-Rod-Like Particles,”Phys. Rev. Lett. , 2760–2763 (1989). A. Leforestier, A. Bertin, J. Dubochet, K. Richter, N. SartoriBlanc, and F. Livolant, “Expression of chirality in columnarhexagonal phases or DNA and nucleosomes,” Comptes RendusChim. , 229–244 (2008). F. Livolant, S. Mangenot, A. Leforestier, A. Bertin, M. de Frutos,E. Raspaud, and D. Durand, “Are liquid crystalline properties ofnucleosomes involved in chromosome structure and dynamics?”Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. , 2615–2633(2006). G. Odriozola, “Revisiting the phase diagram of hard ellipsoids,” The Journal of chemical physics , 134505 (2012). E. Grelet, “Hexagonal order in crystalline and columnar phasesof hard rods,” Phys. Rev. Lett. , 168301 (2008). E. Grelet and R. Rana, “From soft to hard rod behavior in liquidcrystalline suspensions of sterically stabilized colloidal filamen-tous particles,” Soft Matter , 4621 (2016). S. Dussi, M. Chiappini, and M. Dijkstra, “On the stability andfinite-size effects of a columnar phase in single-component sys-tems of hard-rod-like particles,” Molecular Physics , 2792–2805 (2018). H. Wensink and H. Lekkerkerker, “Phase diagram of hard col-loidal platelets: a theoretical account,” Mol. Phys. , 2111–2118 (2009). R. Blaak, D. Frenkel, and B. M. Mulder, “Do cylinders exhibita cubatic phase?” J. Chem. Phys. , 11652–11659 (1999). A. G. Orellana, E. Romani, and C. De Michele, “Speeding upMonte Carlo simulation of patchy hard cylinders,” Eur. Phys. J.E , 51 (2018). A. Stukowski, “Visualization and analysis of atomistic simulationdata with OVITO-the Open Visualization Tool,” MODELLINGAND SIMULATION IN MATERIALS SCIENCE AND ENGI-NEERING (2010), 10.1088/0965-0393/18/1/015012. H. B. Kolli, G. Cinacchi, A. Ferrarini, and A. Giacometti, “Chi-ral self-assembly of helical particles,” Faraday Discuss. , 171–186 (2016). P. D. Duncan, M. Dennison, A. J. Masters, and M. R. Wilson,“Theory and computer simulation for the cubatic phase of cutspheres,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. ,1–11 (2009). Y. Liu and A. Widmer-Cooper, “A versatile simulation methodfor studying phase behavior and dynamics in colloidal rod androd-polymer suspensions,” The Journal of chemical physics ,244508 (2019). S. C. McGrother, D. C. Williamson, and G. Jackson, “A re-examination of the phase diagram of hard spherocylinders,” J.Chem. Phys. , 6755–6771 (1996). J. M. Polson and D. Frenkel, “First-order nematic-smectic phasetransition for hard spherocylinders in the limit of infinite aspectratio,” Physical Review E , R6260 (1997). M. Marechal, S. Dussi, and M. Dijkstra, “Density functional the-ory and simulations of colloidal triangular prisms,” The Journalof chemical physics , 124905 (2017). J. A. C. Veerman and D. Frenkel, “Phase behavior of disklikehard-core mesogens,” Phys. Rev. A , 5632–5648 (1992). A. Repula, M. Oshima Menegon, C. Wu, P. van der Schoot, andE. Grelet, “Directing liquid crystalline self-organization of rodlikeparticles through tunable attractive single tips,” Phys. Rev. Lett.122