Monte Carlo simulation of cylinders with short-range attractions

Cylindrical or rod-like particles are promising materials for the applications of fillers in nanocomposite materials and additives to control rheological properties of colloidal suspensions. Recent advances in particle synthesis allows for cylinders to be manufactured with short-ranged attractions to study the gelation as a function of packing fraction, aspect ratio and attraction strength. In order to aid in the analysis of small-angle scattering experiments of rod-like particles, computer simulation methods were used to model these particles with specialized Monte Carlo algorithms and tabular superquadric potentials. The attractive interaction between neighboring rods increases with the amount of locally-accessible surface area, thus leading to patchy-like interactions. We characterize the clustering and percolation of cylinders as the attractive interaction increases from the homogenous fluid at relatively low attraction strength, for a variety of aspect ratios and packing fractions. Comparisons with the experimental scattering results are also presented, which are in agreement.


I. INTRODUCTION
Anisotropic shapes are ubiquitous in nature, often conferring unique adaptations over more symmetric counterparts. They manifest on a variety of length scales, which range from the internal structure of the atom, to the domains within proteins, to the internal organization of cells, and beyond. Anisotropy has become well appreciated in the field of colloidal science 1,2 for its control over the material properties of suspensions. [3][4][5][6][7] By tuning the particle shape and taking advantage of anisotropy, the self-assembly of colloidal particles can be carefully controlled.
Perhaps the simplest perturbation conferring anisotropy onto an otherwise isotropic object, is the elongation of a sphere into a spherocylinder. These rod-like particles occur in many materials such as cements, paints, pharmaceuticals and consumer products. They have been shown to be facile building blocks in the assembly of nanoscale polyhedrons, 8,9 and are also observed in viruses. 10,11 Historically, the study of athermal, hard rods has been the subject of great interest because the anisotropic shape of the rods can lead to orientational ordering transitions on the basis of entropy alone due to its transfer between rotational and translational modes. [12][13][14] As the aspect ratio of the rods increases, their modes progressively become decoupled, which facilitates the formation of a variety of morphologies. 15 For instance, interactions between grafted ligands can introduce greater complexity, and even be used to confer directional interactions. 7 The addition of small polymers or other colloids to suspensions of rods can alternatively introduce depletion effects, which are also entropic in nature. 16,17 The latter instance has been the source of much theoretical interest, [17][18][19][20][21][22][23][24][25][26] as depletion effects are estimated to play a significant role in conditions where rods occur naturally, as in cellular environments, 27 and in the manufacturing of nanocomposites 28 and liquid crystals. 1,29,30 Furthermore, understanding depletion effects may act as a surrogate for a plethora of other short-ranged interactions, regardless of their cause. [31][32][33][34] These short-range attractions may lead to percolation which has implications for dynamical arrest [35][36][37][38][39] and electrical conductivity. 40,41 Although the connection between rigidity percolation, coordination number and gelation has been studied for isotropic particles, 42 the phase diagram for rodlike particles 32 possesses an additional degree of freedom to explore relative to spherical particles (i.e., the aspect ratio). The focus of this work is to study clustering and percolation over a variety of attraction strengths, volume fractions and aspect ratios.
Owing to the symmetry of the spherical caps at their termini, it is generally simpler to both simulate and synthesize colloidal spherocylinders rather than cylinders. Recent advances in experimental synthesis, however, have led to the ability to produce silica-based colloidal rods that have one spherical end while the other is flat. 43,44 Although the termini of a cylinder are both flat, this intermediate structure reflects the growing necessity for theoretical and simulation methods that can be used to study the behavior of these systems.
To this end, we performed computer simulations of cylindrical particles with flat-ends and short-ranged attractions. We achieved efficient simulation of flat-end cylinders (not spherocylinders) by developing a numerical scheme to compute the interaction between convex superquadric solids of revolution. The surfaces of superquadric solids may be described analytically, which has been exploited in the past to study a subset of convex superquadrics known as superballs that can smoothly interpolate from a cube to a sphere and an octahedron. [45][46][47] The addition of short-ranged attractions due to depletion has also been investigated recently. 48 Superquadric shapes also have been simulated as hard rigid bodies 49 and with attractive patches. 50 Similarly, analytical solutions for the interaction potential between convex superquadric shapes and a planar wall under the influence of depletion have also been investigated. 51 In this computational work, we study cylinders with an attraction range that is approximately 4% of the cylinder diameter, which have also been the focus of recent experimental studies. 43,52 Thus, these systems represent an opportunity to compare experiment and simulation. It is also worth noting for systems with such short interaction ranges that the liquid state is expected to be metastable with respect to the solid. 53 In addition, this study is limited to packing fractions of 22.5% or less, which is below the nematic and smectic regions of the phase diagram of hard spherocylinders. In this work, we study the clustering and percolation as a function of aspect ratio, packing fraction and attraction strength.
Although analytical methods have been developed to simulate hard cylinders, 54,55 including those with patchy interactions requiring iterative procedures, 50,56 we use a tabulated potential to speed up the simulation of superquadrics with attractions due to the required calculation of excluded volume overlaps between pairs of cylinders. This excluded volume overlap is the most computationally intensive contribution to the interaction potential. The tabular potential allows the interactions to be computed only one time for a given relative orientation and position before the simulations begin. With this procedure, the expensive calculation of the interaction between pairs of anisotropic particles is replaced by a query of the stored table. The relative efficiency of the tabular potential with respect to analytical pair wise calculations depends upon the specific model and the resolution of the table. Here, we first validate the simulations using the tabulated potential by calculating virial coefficients and scattering profiles. The virial coefficients by themselves yield a cross over from end-end dominant attractions to side-side dominant attractions as a function aspect ratio, and scattering profiles are compared with experiment. We then use the simulations to study metrics which signal percolation and dynamic arrest (e.g., cluster percolation probability, coordination number and orientational order), by simulating colloidal cylinders over a wide range of aspect ratios, attraction strengths, and packing fractions.
The outline of this paper is as follows. In Section II, we describe the model for the cylindrical particles and the interaction potential. To explore the effects of aspect ratio on the interaction potential, virial coefficient calculations are presented in Section III. The Monte Carlo methods, including Wang-Landau sampling and expanded ensemble approaches, used to simulate fluids composed of cylindrical particles are described in Section IV. In Section V, we compare simulated and measured scattering profiles to further test the model and simulation method. We then compare various metrics used to identify clustering and percolation in Section VI. Finally, conclusions are made in Section VII and the tabular potential for the computational modeling approach is detailed in the Appendixes.

II. MODELS
Cylindrical particles were modeled with the superquadric equation, 57 x a x 2 ϵ 2 + y a y 2 ϵ 2 where x, y and z are the Cartesian coordinates, ϵ 1 and ϵ 2 parameters determine the curvature, and the a i parameters, with i = x, y, z, determine the maximum extent of the shape from the center, in the given dimension. To model cylinders, the constraints a x = a y , ϵ 2 = 1 are imposed, such that the length of the cylinder, L = 2a z , which extends in the z direction in the body-fixed frame of reference, and the diameter of the cylinder, D = 2a x . The parameter ϵ 1 = 0.1 was chosen such that the ends of the cylinders were flat. An example is shown in Figure 1. As shown in the figure, the relative position and orientation of two solids of revolution, i and j, are given by the center separation distance, r, and three angles, θ i , θ j , and ψ. The first two angles, θ i and θ j , are the angles formed between the center separation vector and the axis of revolution of a particles i and j. The third angle, ψ, is the dihedral angle formed by the three vectors described above (i.e., ψ is the angle between two planes formed by the center separation vector and the axes of revolution of each particle). The potential energy, U, between a pair of cylindrical particles, i and j, is U = U h + U a + U e , which includes a hard-particle steric repulsion, U h , a short-range attractive interaction, U a , and a screened electrostatic repulsion, U e . These contributions are described below.
The short-ranged attractive interaction was modeled with the pair-wise implicit-depletant potential, where ΔV ex is the excluded volume overlap (for a hard sphere of radius R g = 0.04D) between two cylinders, i and j, and ΔV ex m is the maximum of ΔV ex over all non-overlapping positions and orientations between a pair of rods. The parameter ϵ is the scale of the interaction strength. Note that the potential goes to zero at the interaction center separation cut off distance, r c , when the edges of the excluded volumes of the particles touch. A physical interpretation of Equation (3) may be obtained by considering the Asakura-Oosawa depletion force 16 induced by a dilute solution of particles with a packing fraction of ϕ d . The potential energy is given by U a = − Δ V ex (r, θ i , θ j , ψ) 4 3 πR g 3 ϕ d k B T , where k B is the Boltzmann constant and T is the temperature. Equation (3) is then obtained when ϵ is equal to the maximum of the absolute value of U a . The algorithm used to compute ΔV ex (r, θ i , θ j , ψ) is described in Appendix C. The excluded volume of a particle was approximated by Equation (1) with the size parameters a i → a i + R g , which is accurate for a x ⪢ R g . For a x ⪢ R g , the implicit depletant potential is more computationally efficient than explicitly simulating the depletant molecules, although there are alternative methods. 58 With this implicit depletion treatment, many-body effects can arise when the excluded volume of more than two particles overlap, 59 but these effects are expected to be small for a x ⪢ R g . 48 This model is amenable to studying systems with attractions that are physically distinct from the depletant interaction, but are also similarly short-ranged. For example, van der Waals attractions in colloidal systems (e.g., recently synthesized cylindrical particles 43 ) may be accurately modeled with this short ranged potential by applying extended corresponding states. 60,61 Finally, we include an electrostatic repulsive term, U e , which we assume to have a Yukawa form given by where the parameter κD = 10 4 was chosen such that U e only contributes over an extremely short distance of approximately 5D × 10 −4 . Note that this is an approximation of the electrostatic interaction, where a more rigorous treatment might require accounting for the Gaussian curvature by using the Derjaguin approximation. 51 We note that simulations were also performed without Eq. (4) and we did not see a difference in the computed properties. This term is relatively insignificant for the results of this study due to the size of the parameter κD = 10 4 , relative to the strength of the attractive interactions, which we will discuss quantitatively in the following paragraph. But we used this term to make the model more amenable to molecular dynamics simulations with continuous potentials.
Examples of the potential energy as a function of center-center distance for four different relative orientations are shown in Figure 2. Additionally, examples of the distance between centers at contact, r h (θ i , θ j , ψ), and the attractive interaction at contact, U a (r h , θ i , θ j , ψ), are provided in Figure 4, with sample configurations shown in Figure 3. The configuration where the cylinders are parallel and the axes of symmetry lie in a plane is given by cos θ i = cos θ j = 0 and cos ψ = 1. As shown in Figure 4 and illustrated in Figure 3(a), this configuration is highly favorable. Also note, in the context of the relative insignificance of Eq. (4) discussed in the previous paragraph, the softness introduced by this term is visually imperceptible and the potential well reaches very nearly a value of −ϵ because the contribution of Eq. (4) is relatively small compared to the contribution of Eq. (3). The cylinders may also tilt slightly, as shown in Figure 3(b) and still interact favorably. In addition, because the ends of the cylinders are relatively flat, the attraction of the flat ends, as illustrated in Figure 3(c), is nearly as strong as the parallel configuration shown in Figure   3(a), for an aspect ratio of about L/D « 3. As shown in Figure 5, ΔV ex m increases linearly with L/D for aspect ratios of L/D ≥ 3, because the overlap volume in the parallel configuration grows with the length of the cylinder. Note that the slight decrease in ΔV ex m . shown in Figure   5 for L/D ≤ 3 is controlled by ϵ 1 = 0.1. The curvature at the ends of the cylinders gradually increases with length, which leads to a slight decrease in end-end attraction as length increases.
The interaction potential used here requires the calculation of two computationally expensive quantities, r h and ΔV ex . The latter is significantly more expensive than the former, and it is not computationally viable to calculate ΔV ex on-the-fly. It follows that any effort to speed up simulations of these types of systems should focus on speeding up ΔV ex . Thus, we have adopted to tabulate these interactions as described in Appendix A. Even during the tabulation process, the time to calculate r h is negligible to that of ΔV ex . The multidimensional tabulation of the anisotropic potential accelerated the simulation by precomputing these computationally expensive quantities described above over a range of relative positions and orientations, and then interpolating from the stored values during the course of the simulation.

III. SECOND VIRIAL COEFFICIENT
The second virial coefficient is a useful measure of the relative balance of repulsive and attractions. It is also the central quantity of interest in applying extended corresponding states to compare results obtained from different computational models or experimental measurements. The second virial coefficient, B 2 , for solids of revolution is given by 62,63 where β = 1/k B T and the definitions of the relative orientation variables are described in Section II and Figure 1. In practice, the virial coefficients were computed via Eq. (5) by numerically integrating the tabular potential, as described in Appendix D. The second virial coefficient was split into contributions from the hard particle potential, U h , and the continuous potential, U a + U e , such that B 2 = B 2 The second virial coefficients for hard cylinders with no attractions (i.e., βϵ = 0) are shown in Figure 6. These hard-particle virial coefficients are expected to increase with the volume of the particle, here determined by L. Comparing with theoretical calculations for hard spherocylinders, 64 there is expected to be a L 2 dependence, as appears to be the case in Figure 6. Second virial coefficients for cylinders with various aspect ratios as a function of the interaction strength, βϵ, are shown in Figure 7. As the attractive interaction is increased, the virial coefficients become negative. Note that the dependence of B 2 with respect to the aspect ratio follows a similar trend to ΔV ex m (see Figure 5).
The theta solvent condition, B 2 (βϵ θ ) = 0, is a convenient parameter to characterize the various shapes, and is shown in Figure 8 and Table I. In practice, the theta solvent condition was obtained by computing B 2 over a large range of values of βϵ with a spacing of 5 × 10 −4 . The reported value for βϵ θ is the average of the following two values of βϵ: the maximum βϵ for which B(βϵ) > 0, and the minimum βϵ for which B(βϵ) < 0. The reported error for βϵ θ is the difference between these two values.
The theta solvent condition, B 2 (βϵ θe ) = 0, may be used to identify a cross over from end-end dominant attractions to side-side dominant attractions at approximately L/D = 3, as seen in both Figure 5 and 8. For aspect ratios L/D < 3, βϵ θ decreases because the end-end contact interaction dominates energetically, and the end-end contact depends slightly on L/D, as discussed previously in Section II for Figure 5. For aspect ratios L/D ≥ 3, as the aspect ratio increases, βϵ θ increases toward an asymptotic value. The location of this cross over depends on both the shape of the end of the cylinder, and also the interaction range, R g . For the remainder of this work, we focus on cylindrical shapes with 3 ≤ L/D ≤ 8, because a thorough study of particle shapes with end-end dominant interactions and the effect of end-shape is beyond the scope of this study. In addition, the focus of this work is on aspect ratios of cylindrical particles that were recently synthesized. 43

IV. WANG-LANDAU MONTE CARLO SIMULATIONS IN AN EXPANDED ENSEMBLE
Computer simulations employing specialized techniques to sample highly-attractive and short-ranged interactions are used to investigate packing, non-equilibrium gelation and macroscopic phase separation of the cylindrical particles. Wang-Landau (WL) sampling is a flat-histogram method used to obtain the probability distribution function of some specified order parameter. By setting the order parameter to βϵ, the expanded ensemble 65,66 effectively functions as parallel tempering aided by the flat-histogram methods to enhance sampling of transitions with large energy barriers. Note that in this work, the attraction strength, βϵ, is trivially related to the temperature, and therefore the expanded ensemble in βϵ is analogous to the temperature expanded ensemble. Wang-Landau Monte Carlo with an expanded ensemble in attraction strength enhances sampling of transitions between macroscopic phases and microscopic structural changes. 48 In addition, specialized Monte Carlo trials are utilized to enhance sampling of clusters which are expected to form due to the strongly attractive interactions. These simulations were conducted using the FEASST simulation package. 67 The following Monte Carlo trials were employed. Translations and rotations of particles were attempted with equal probability. For expanded ensemble simulations, βϵ, was increased or decreased by a fixed amount, ± δβϵ, subject to Metropolis acceptance criteria.
Collective trial moves were also implemented to facilitate convergence in systems with short-ranged attractions that self-assemble. 46,68 This included rigid-body translations and rotations of clusters, where clusters were defined as all particles with excluded volume overlap, r h (θ i , θ j , ψ) < r < r c (θ i , θ j , ψ), with at least one other particle in the cluster, obtained via a recursive flood-fill algorithm. To obey detailed balance, cluster moves resulting in a particle joining a different cluster were rejected.
The geometric cluster algorithm (GCA) was also used. 46,69,70 The GCA is a rejection-free algorithm that collectively moves particles, and results in better sampling of clusters of particles than traditional single particle moves. The algorithm proceeds as follows. A particle and a pivot point in space are randomly selected, and the particle is reflected about the pivot. All other particles which interact with the pivoted particle, in both the old and newly pivoted positions, are then attempted to be pivoted with a probability related to the pair interaction energy between the two particles. Each attempted pivot was carried out recursively until all the interacting particles were attempted to be pivoted. To avoid inefficient moves involving most of the particles in the system, the pivot point was confined to a cubic box centered on the first randomly selected particle. The size of this bounding cubic box was tuned during the course of the simulation in order to obtain an average target number of particles involved in a pivot, set to N p /5, where N p is the maximum number of particles in the simulation. While the rigid cluster moves could not create or destroy clusters due to detailed balance, the GCA does not suffer from this limitation. For anisotropic particles, pivots about a point, as implemented in this work, result in reflections of the particle orientation, and cannot sample all particle orientations without other Monte Carlo trials (e.g., rigid body single-particle and cluster rotations). Note that this is not a deficiency of the GCA method, because particle reflections about a plane and line may be used to sample arbitrary orientations. 71 The weights for the probability of selecting each trial type are provided in Table II. For each Monte Carlo trial that involved movement of particles, the parameter associated with the maximum change was optimized, via a 5 % change every 10 6 trials, to yield approximately 25 % acceptance of the trial move.
Simulations were conducted for aspect ratios of L/D = [3,4,5,6,7,8] and 15 different volume fractions in the interval ϕ = [0.05, 0.225], with a spacing of 0.0125. The initial configuration was generated by grand-canonical insertion of N p particles in a cubic domain of edge length l b /D = 20 with periodic boundary conditions in order to obtain a given particle volume fraction, φ, rounded to the nearest particle number. Note that volume fractions were computed using particle volumes given in Table I, which were calculated as described in Appendix C. The Wang-Landau update factor, Inf, was initially set to unity, and In order to estimate the standard deviation, all simulations were conducted with 3 identical, independent replicas with different random number seeds. Note that these simulations were well equilibrated and reproducible due in large part to the expanded ensemble flat histogram sampling. Over the course of each simulation, which consisted of 3 -50 billion Monte Carlo trials, we observed that the system was able to sample fluidly over the range of attraction strengths, indicating that freezing did not occur. Additional simulations were conducted for l b /D = 10 and 15 to test the system size dependence of the calculated properties discussed below.

V. SMALL-ANGLE SCATTERING INTENSITY AND STRUCTURE FACTOR
Small-angle scattering experiments are able to probe short-to-intermediate length scale structures. As such, it is useful to compare the simulations of our model with experiments in order to validate the model and also aid in the interpretation of the experimental results. In order to compare the model and simulations with small-angle X-ray or neutron scattering experiments, as related to the recent experimental synthesis of cylinders and scattering measurements, 43 the scattering intensity, I(q), was obtained as the three-dimensional numerical Fourier transform of the scattering density, ρ(r), averaged over all orientations, Ω, 73,74 The scattering density was discretized on a grid with 256 elements along each dimension. For each particle, the center point was used to initialize a flood-fill algorithm to fill the spatial grid with scatters if the center point of the grid was inside the particle. The software package FFTW3 75 was utilized to compute the discrete three-dimensional Fourier transform. Note that discrete Fourier transforms of ρ(r) naturally take into account the periodicity of the simulation domain. The complex conjugate was then averaged over many uncorrelated configurations stored during the course of the simulations. Orientational averaging of the intensity, I(q), was performed with channel sharing. 73 The effective structure factor was obtained by dividing the intensity by the number of particles and by the form factor, which mimics traditional experimental practices. The form factor was computed as the intensity from a single-particle simulation with the same spatial grid resolution.
A comparison between the scattering intensity of experiments and simulations is shown in Figure 9 for an aspect ratio of L/D = 4 and a volume fraction of ϕ = 0.11 for experiments and ϕ = 0.1125 for simulations. To facilitate comparison with experimental scattering intensity results, 43 smearing functions were utilized to represent the aperture of the instrument and the polydispersity, 76,77 as described in Appendix E. There is deviation between the experiments and simulations at the lowest qD value calculated, but this is expected due to finite size effects of the simulation periodic boundary conditions. In addition, small differences observed for qD values above an approximate value of 7 arise from polydispersity of the rods, which were only qualitatively accounted for by smearing the intensity of the monodisperse simulations. In addition to polydispersity, other deviations between experiment and simulation may be due to the approximate modeling of the attractive interactions, which are due to the polymer brush in experiment. This comparison allows us to qualitatively relate the attractive strength, βϵ, of the simulations to the temperature in the experiments.
The effective structure factor is shown in Figure 10 for an aspect ratio, L/D = 4, and various attraction strengths, βϵ, and volume fractions, ϕ. As the attractive strength increases, well defined peaks appear in the structure factor, and the low-qD values increase. Structure factor values which exceed unity at the low-q values (i.e., the largest length scales) are typically associated with macroscopic phase separation, cluster formation or aggregation. 78,79 Here, the smallest q value attainable is qD = 4π/l b , where l b is the side length of the cubic simulation box. The structure factor for the smallest qD values are summarized in Figure 11 for L/D = 4. At low attraction strength, βϵ, S eff (4π/l b ) is less than unity, but increases with increasing βϵ. The value of βϵ at which S eff (4π/1 b ) = 1 tends to increase with increasing packing fraction, φ. The increase of S eff (4π/1 b ) upon decreasing density could be an indication that the fluid is metastable with respect to the solid. We do not impose a constraint which forbids the formation of the solid, as previously used in computations with the adhesive hard sphere model. 80 This is why we do not report on the critical temperature as a function of aspect ratio, as has been studied previously for ellipsoids with an order of magnitude longer range attraction. 81 As noted previously, our flat-histogram simulations sampled a wide range of attraction strengths without any difficulties, consistent with the avoidance of the stable solid phase. The low-q S eff results for cylinders with other aspect ratios are summarized in Figure 12, where we plot the attraction strength corresponding to S eff (4π/1 b ) = 1 for a cylinder of given L/D and packing fraction, φ. An interesting trend that can be gleaned from Figure 12 is that cylinders with smaller aspect ratios tend to exhibit low-q structure factors greater than unity more readily than higher aspect ratio cylinders. Analysis of the peaks in the structure factor may be aided by computation of clustering and orientational order, which is presented in the following section.

VI. CLUSTERS, ORIENTATIONAL ORDER AND PERCOLATION
In the final results section of this manuscript, a detailed analysis of the clustering, percolation and orientational order of the attractive cylinders is presented. In particular, two different measures of clustering will be shown. The first is the average coordination number, n c , which is a local quantity that depends on the nearest neighbors. The second is the probability of connectivity percolation, which is a global quantity that encompasses the entire simulation box, and has been previously shown to be relatively system size independent at the point of 50% probability 40 and verified for a few of our simulations. In addition, a local and global measure of orientational order will also be shown. The local measure is the angle between axes of revolution of nearest neighbors, and the global measure is the nematic order parameter.
Rigidity percolation has been proposed to occur at 〈n c 〉 = 2.4, 82 which has been shown previously to indicate gelation and dynamical arrest in isotropic model systems. 39,42 The coordination number, n c , was computed as the average number of particles, denoted by i, whose exclusion volume overlaps with a central common particle, denoted by j ≠ i (e.g., when r h (θ i , θ j , ψ) < r < r c (θ i , θ j , ψ)). As shown in Figure 13 for L/D = 4, the coordination number increases with increasing attractive strength, βϵ, as expected. For low attraction strengths, when n c ≲ 2.4, the coordination number tends to increase with volume fraction; however, this trend reverses when n c ≳ 3. Figure 14 shows the distributions of n c at conditions where the average, 〈n c 〉 ≈ 2.4. Because the overall shape of these distributions resemble those of previously published spherical particles, 39 the underlying physics of rigidity percolation in the different systems could be similar. However, a more detailed comparison with the distributions in the literature is complicated by the difficulty in finding the conditions where n c = 2.4, precisely. An example configuration of the cylinders with 〈n c 〉 = 2.4 is shown in Figure 15 for L/D = 4, ϕ = 0.1125 and βϵ = 15.44. Figure 16 shows that the attraction strength, βϵ, at which 〈n c 〉 = 2.4, increases with increasing aspect ratio, L/D.
In addition, the attractive strength corresponding to 〈n c 〉 = 2.4 tends to decrease with increasing volume fraction. Also notice that a coordination number of 2.4 can be found at all packing fractions, even when the packing fraction is too low to allow percolation. This is due to the formation of clusters. The coordination number criterion for rigidity percolation should be used with some caution for freely diffusing particle systems.
Percolation of clusters also serves as an important metric for relating structural quantities to dynamical arrest. Clusters of particles were defined as all particles with overlap of their excluded volumes, r h (θ i , θ j , ψ) < r < r c (θ i , θ j , ψ), with at least one other particle in the cluster, obtained via a recursive flood-fill algorithm. A cluster was considered percolated if the particles in the cluster were connected to their own periodic images. 40 The percolation of clusters of cylindrical particles was investigated as a function of the interaction strength, βϵ, the particle volume fraction, φ, and aspect ratio, L/D. In Figure 17, we plot the percolation probability for L/D = 4. As the volume fraction, ϕ, and interaction strength, βϵ, increase, the probability of percolation increases. The point at which there is a 50 % probability of percolation was shown to be system size independent in Ref. 40 and verified for a few of our simulations (not shown). Figure 18 shows that the attraction strength, βϵ, corresponding to 50 % percolation probability generally increases with aspect ratio, and decreases with packing fraction. There is a minimum packing fraction below which percolation does not occur.
The orientation of the rods was investigated with two different order parameters. One measure is the angle between the axes of revolution, denoted by θ, between neighboring particles with excluded volume overlap, r h (θ i , θ j , ψ) < r < r c (θ i , θ j , ψ). Thus, this is a measure of local orientational order. The absolute value of the cosine of this angle, | cos θ|, is unity for parallel configurations, and the ensemble average of this quantity, 〈| cos θ|〉, has an average expectation value of 0.5 for random orientation. 74 As shown in Figure 19, the local orientational order increases with attractive strength, βϵ, and slightly decreases with packing fraction, φ. This trend for 〈| cos θ|〉 is very similar to that observed for the coordination number in Figure 13, because similar physical mechanisms dictate the behavior of both of these local ordering measures. For L/D = 4, 〈| cos θ|〉 does not reach unity even for high attractive strengths, βϵ. Although the parallel orientation is strongly favored, a perpendicularly-oriented cylinder may interact favorably with a cluster of parallel cylinders which present an approximately flat face.
In addition, the nematic order parameter was computed as the largest eigenvalue of the tensor Q = (2N) −1 Σ N (3u i u i −I). 40 The nematic order parameter is unity when rods are fully aligned, and zero in a perfectly isotropic fluid. Figure 20 shows that the nematic order parameter rarely reaches a value above 0.5. Nematic order decreases with packing fraction because the glassy and gel-like configurations are less likely to sample the states with global orientational order. Furthermore, for larger aspect ratios, particles are more likely to orient globally.

VII. CONCLUSIONS
The structural properties of cylindrical particles with short-range attractions were investigated using Wang-Landau Monte Carlo computer simulations over a wide range of attraction strengths, volume fractions and aspect ratios. The second virial coefficients were calculated for different aspect ratios and attraction strengths. Interestingly, B 22 exhibits a crossover from end-end dominant attractions to side-side dominant attractions as a function of L/D. The scattering intensity compared well with experiment, and the low wave-number effective structure factor provides a useful measure for the structure of the cylinders. The conditions where percolating clusters form were investigated with two independent measures: connectivity of the clusters across periodic boundaries and also 2.4 coordination number, which has been related to rigidity percolation. 82 The orientations of the rods were also investigated with two different measures. First, the relative orientation of the axes of symmetry of neighboring particles, which is a local measure of orientational order, was reported. Then, the nematic order parameter was also shown, and this is a global measure of orientational order. These measures suggest the rods may possess orientational order at length scales corresponding to nearest-neighbors, but the global orientational order of the entire system rarely, if ever, reaches a full nematic phase. Instead, small crystallites with local order pack together into more disordered configurations.
This manuscript lays the groundwork for more detailed comparisons with experimental results of the gelation line in the temperature-packing fraction projection of the phase diagram obtained from recently synthesized cylindrical particles. 43,44 In particular, a method has been developed to collapse this gelation line with respect to aspect ratio. 52 This method involves fitting the structure factor to a model for an adhesive hard sphere in order to quantify the temperature of gelation for general potentials and models.
Other future work includes a more thorough analysis of particles with L/D < 3, where endend interactions become dominant. For example, one interesting phenomena observed at L/D = 1 was that the particles stacked into larger cylinders which resembled the rouleaux of red blood cells. 83 Finally, the methodology for efficiently simulating particles according to the superquadric equation will be used to investigate the phase behavior of a variety of different anisotropic-shaped particles in future studies.

APPENDIX A:: TABULAR POTENTIAL
The tabular potential was constructed to obtain, r h (θ i , θ j , ψ), r c (θ i , θ j , ψ) and U a (r, θ i , θ j , ψ) efficiently during the course of the simulation, as described in Section II, where r c is the center separation distance at which the edges of the excluded volumes of the particles touch. The construction of the tabular potential is similar to that described in the Appendix of Ref. 48 for two dimensional particles. In this work, the algorithm to find the contact distance, as described in Appendix B, also returned the surface point which was closest to the contact.
This surface point was subsequently used to obtain the initial voxel of the recursive flood-fill algorithm to compute the overlap volume, as described in Appendix C, over the range of n z + 1 values of z = r − r ℎ r c − r ℎ , z ∈ [0, 1].
In this work, the tabulated angles were evenly spaced over cos θ i , cos θ j and cos ψ, in order to avoid the arccos operation during the simulation. Because the shapes are symmetric about the plane which possesses the center point and is perpendicular to the axis of symmetry, negative values of cos θ did not need to be tabulated, and cos θ ∈ [0, 1] over a range of k + 1 evenly-spaced values. No symmetry operation on ψ was necessary because cos ψ ∈ [−1, 1] over a range of 2k + 1 evenly-spaced values. Nearly half of the possible combinations of cos θ i and cos θ j may be ignored by swapping i and j if cos θ i > cos θ j , which is also required for unique computation of the potential energy due to the randomness in Monte Carlo integration. In this work, k = 150, which was shown to be sufficient in Ref. 48 and resulted in 6 863 101 tabulated orientations. As a measure of the expected error in the tabular potential due to Monte Carlo integration, two tables with identical parameters, and L/D = 3, but different random number seeds were compared. The square root of the sum of the squared difference of U a (r h )/ϵ was 4 × 10 −5 , and these squared differences did not vary greatly with orientation. During the course of the simulations, values of r h , r c and U a are obtained for arbitrary orientation by multi-dimensional linear interpolation. Example twodimensional slices of the tabular potential are shown in Figure 4. In addition, with n z = 4, the tabular potential stores in less than 1 gigabyte of memory. This smaller value of n z was used to capture the qualitative shape of the implicit-depletant potential, but more precision was not necessary when one intends to use extended corresponding states to compare with experiments, by matching second virial coefficients between models and experiments.

APPENDIX B:: HARD CONTACT
The following algorithm used to obtain the distance between centers at hard contact is less efficient than alternative algorithms. 54,55 But the efficiency of this algorithm is not an issue for two reasons. First, this algorithm is not utilized during the simulations because the tabular potential described in Appendix A requires these contact distances to be computed only once before the simulations begin. Second, the bottleneck in the interaction between the particles is the attractions, and not the hard contact.
The distance between the centers at hard contact is obtained by decorating the surface of the particles with many points, and checking whether those surface points overlap with the other particle. For a given relative orientation, r h (θ i , θ j , ψ) is obtained by incrementally moving a particle center by an amount δ until the overlap status changes, and then halving δ and changing its sign, until δ < 10 −15 . Surface points are obtained by solving the implicit equation, Equation (1) in evenly-spaced grids for two of the three independent variables for the following three cases: the x-y plane, the x-z plane, and the y-z plane. Each of these grids was of size (n s + 1, n s + 1) and spanned the range of the respective size parameter "a" for the relevant dimension. In addition, the parametric equations with longitude and latitude parameters and a similar-sized grid, was also used to obtain a fourth set of surface points. Combining these four cases ensures that there are no large gaps in coverage of the surface.
In this work, n s = 100, which resulted in 64 838 points for the cylinders, and 57 046 for a sphere, because the edges of the Cartesian planes may have no solution for Equation (1).
The number of surface points was tested by the convergence of the contact function of two spheres of unit diameter (a x = a y = a z = D/2 and ϵ 1 = ϵ 2 = 1) for 57 046 orientations.
Because each sphere was represented by a mesh of surface points, each orientation may deviate from the theoretical center separation distance of unity, and these deviations diminish as the number of surface points increase, as shown in Appendix D. The square of the deviations for each orientation was summed and plotted as a function of n s . Convergence was assessed by the quantity σ r = ∑ n N = 57046 (1/N)(r ℎ (i) − 1) 2 , which is the standard deviation from the theoretical result. As shown in Figure 21, the contact function calculation with n s = 100 represents a sphere within a standard deviation of approximately 3.4D × 10 −5 .

FIG. 21.
Convergence of the number of surface points for a sphere, a x = a y = a z = 0.5 and ϵ 1 = ϵ 2 = 1 averaged over 57 046 orientations.

APPENDIX C:: OVERLAP VOLUME
The following algorithm was used to compute the overlap of the excluded volumes of two particles during the tabulation of the potential before the simulations were run, as described in Appendix A. Because this algorithm is not used during the course of the simulations, the efficiency of the algorithm is not the simulation bottleneck as long as the tabular potentials do not take too long to compute for a given pair of shapes. Here we used a combination of a recursive flood-fill algorithm and Monte Carlo integration. In this algorithm, space is divided into a grid of voxels with side length of d v . To begin, a voxel is selected which is known to contain an overlap. Determining this voxel is described more in the following Appendix A on construction of the tabular potential. This selected voxel, and all neighbors of this selected voxel are then subject to Monte Carlo integration with n t random test points.
A test point is inside a shape described by Equation (1)  . The number of test points in a voxel, n t , is chosen to be the smallest value with no appreciable overhead cost for the flood-fill algorithm. For example, if n t is too small (e.g., n t = 1) then the algorithm slows down, for a given average point density, due to increased memory to store the voxels for the flood-fill algorithm. But if n t is too large, for fixed point density, then the flood fill algorithm is not fully utilized. The side length of the voxel, d v , was tested against the overlap of two contacting spheres with diameter of 1.08 and a center separation of 1 averaged over 57 046 realizations of random numbers. Theoretically, the overlap volume is known to be given by V o = π(2R − d) 2 (d + 4R)/12 for two equalsized spheres of radius R = 0.54D and separation distance d = D. Thus, convergence may be assessed by the quantity σ V = ∑ i 2 , which may be thought of as the percent standard deviation. As shown in Figure 22, σ V reaches a plateau that depends on the number of surface points, n s . This plateau occurs when the relative precision of the contact function due to n s becomes higher than the precision due to d v . For the simulations of the cylinders, d v = 0.05D and n t = 512, which has a similar point density as d v = 0.03D and n t = 128 for the sphere case, where overlap volumes were computed within an average precision lower than 0.01 percent.
The volume of a single particle, V p , was computed by the same algorithm as described here, with n t = 10 5 and d v = 0.06D. The resulting values of V p , reported in Table I were  Convergence of the voxel side length, d v with n t = 128 for a sphere, a x = a y = a z = 0.54 and ϵ 1 = ϵ 2 = 1 for (red +) n s = 50 and (black x) n s = 100.

VIRIAL COEFFICIENT
The tabulated interaction described in Appendix A can be used directly to calculate the second virial coefficient, B 2 , using Equation (5). Note that the tabulated potential applies to solids of revolution possessing an axis of rotational symmetry and a plane of symmetry at the center point and perpendicular to the axis of symmetry. In other words, we use the following 0 π f sin θdθ = 2 0 π/2 f sin θdθ, The tabular potential described in Appendix A may be conveniently integrated according to this equation with some changes of variables. In particular, 0 π/2 f sin θdθ = 0 1 fd(cos θ), r ℎ r c f(r, …)r 2 dr = 0 1 f(z, …)r 2 (r c − r ℎ )dz, where z = r − r ℎ r c − r ℎ and r c is the center separation "cut off" distance at which the edges of the excluded volumes of the particles touch. Note that no change of variable for ψ was used to avoid numerical instabilities due to a singularity. In addition, the second virial coefficient was split into contributions from the hard particle potential, U h , and the continuous potential, U a + U e , such that B 2 = B 2 ℎ + B 2 ae . For the hard particle contribution, 0 ∞ f ℎ (r, θ i , θ j , ψ, βϵ)r 2 dr = − 1 3 r ℎ (θ i , θ j , ψ) 3 , where f h = e −θU h − 1. Thus, the second virial coefficient contribution from the hard particle potential, B h , is obtained by substitution of Equations (6), (7), (D1)-(D3), and (D5) into Equation (5), The second virial coefficient contribution from the continuous interaction, B 2 ae , is obtained by substitution of Equations (6), (7), (D1)-(D4) into Equation (5), B 2 ae (βϵ) = − 2 0 1 d(cos θ i ) 0 1 d(cos θ j ) 0 π dψ 0 1 f(z, θ i , θ j , ψ, βϵ)r 2 (r c − r ℎ )dz . (D7) The numerical integration was then performed using a multi-dimensional trapezoidal rule, which is expected to converge well for periodic functions. Because z is tabulated with lower resolution than θ and ψ, the resolution of z used to perform the integration was increased by a factor of 100 with respect to the resolution of the tabular potential. To show that n ez = 100 is sufficient, the second virial coefficient was computed for a number of interpolated values in z, n ez , as shown in Figure 23.

APPENDIX E:: SMEARING OF THE SCATTERING INTENSITY
For slit aperture smearing, For polydispersity smearing, we assume that the diameter of the cylinder follows a Gaussian function and that the scattering intensity scales with the diameter of the cylinder. This assumes that the length and the diameter are correlated and follow the same distribution and that the polydispersity has no effect on the relative packing of all the cylinders. Following these assumptions, the smeared intensity is given by where q j D j = q i D. In order to implement this polydispersity smearing on the simulation data, the integration variable is changed from D j to q j and converted to a summation. In addition, the Gaussian width parameter, σ q , was truncated at 3σ q . Finally, the Gaussian width parameter was normalized such that σ q = Δ q D j , where Δ q = 0.135. Note that we normalize by D j instead of D so that the discretized Gaussian kernel width does not depend on q j . The working equation is given by In order to apply the smearing kernel at the minimum and maximum q values, I(q j ) was extended assuming it was constant.

FIG. 2.
The total potential energy between a pair of cylinders with aspect ratio of L/D = 3 and shortranged attractions, R g /D = 0.04, as a function of the separation distance of the centers for a few orientations. Letter-coded labels refer to structures shown in Figure 3.  The second virial coefficient, B 2 , as a function of the interaction strength, βϵ, for integer value aspect ratios in the range L/D ∈ [1,8].

FIG. 9.
A comparison between scattering intensity from experiment 43 The effective structure factor is shown as a function of qD for an aspect ratio, L/D = 4, and various attraction strengths, βϵ, and volume fractions, φ.

FIG. 11.
The value of the effective structure factor at the smallest accessible value of qD = 4π/l b is shown for various attraction strengths, βϵ, volume fraction, ϕ, and L/D = 4. In addition, the dashed black line shows where S eff = 1.

FIG. 12.
The loci of points where S eff (qD = 4π/l b ) = 1 is shown for a given aspect ratio, L/D, attraction strength, βϵ, and volume fraction, ϕ.

FIG. 13.
The coordination number for a given interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D = 4. The black, dashed line shows a coordination number of 2.4. The loci of points where the coordination number is 2.4, for a given interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D.

FIG. 17.
The probability that there is a percolated cluster for a given interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D = 4. The black, dashed line shows the 50 % probability.  The loci of points when there is a 50 % probability that a cluster is percolated, for a given interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D.

FIG. 19.
The ensemble average of the absolute value of the cosine of the angle between axes of revolution of rods with surfaces within a distance 0.08D as a function of interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D.

FIG. 20.
The nematic order parameter is shown at the highest attractive interaction strength investigated as a function of particle volume fraction, ϕ, and aspect ratio, L/D. Properties of cylinders as a function of the aspect ratio.