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.

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 science1,2 for its control over the material properties of suspensions.3–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–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–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 nanocomposites28 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–34 These short-range attractions may lead to percolation which has implications for dynamical arrest35–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 rod-like particles32 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–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 bodies49 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.

Cylindrical particles were modeled with the superquadric equation,57 

xax2ϵ2+yay2ϵ2ϵ2ϵ1+zaz2ϵ1=1,
(1)

where x, y and z are the Cartesian coordinates, ϵ1 and ϵ2 parameters determine the curvature, and the ai 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 ax = ay, ϵ2 = 1 are imposed, such that the length of the cylinder, L = 2az, which extends in the z direction in the body-fixed frame of reference, and the diameter of the cylinder, D = 2ax. 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 = Uh + Ua + Ue, which includes a hard-particle steric repulsion, Uh, a short-range attractive interaction, Ua, and a screened electrostatic repulsion, Ue. These contributions are described below.

FIG. 1.

A pair of cylinders with L/D = 3, ϵ1 = 0.1, θi = θj = ψ = π/2 and r = 2D.

FIG. 1.

A pair of cylinders with L/D = 3, ϵ1 = 0.1, θi = θj = ψ = π/2 and r = 2D.

Close modal

The hard-particle interaction between two particles, Uh, is given by

Uh(r,θi,θj,ψ)=r<rh(θi,θj,ψ)0rrh(θi,θj,ψ).
(2)

where rh(θi, θj, ψ) is the hard center separation distance at contact, which is computed numerically using Equation (1) as described in  Appendix B.

The short-ranged attractive interaction was modeled with the pair-wise implicit-depletant potential,

Ua(r,θi,θj,ψ)ϵ=ΔVex(r,θi,θj,ψ)ΔVexm
(3)

where ΔVex is the excluded volume overlap (for a hard sphere of radius Rg = 0.04D) between two cylinders, i and j, and ΔVexm is the maximum of ΔVex 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, rc, 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 force16 induced by a dilute solution of particles with a packing fraction of ϕd. The potential energy is given by Ua=ΔVex(r,θi,θj,ψ)43πRg3ϕdkBT, where kB 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 Ua. The algorithm used to compute ΔVex(r, θi, θj, ψ) is described in  Appendix C. The excluded volume of a particle was approximated by Equation (1) with the size parameters aiai + Rg, which is accurate for axRg. For axRg, 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 axRg.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 particles43) may be accurately modeled with this short ranged potential by applying extended corresponding states.60,61

Finally, we include an electrostatic repulsive term, Ue, which we assume to have a Yukawa form given by

Ue(r,θi,θj,ψ)ϵ=Dreκ[rrh(θi,θj,ψ)],
(4)

where the parameter κD = 104 was chosen such that Ue 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 = 104, 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, rh(θi, θj, ψ), and the attractive interaction at contact, Ua(rh, θ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, ΔVexm 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 ΔVexm 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.

FIG. 2.

The total potential energy between a pair of cylinders with aspect ratio of L/D = 3 and short-ranged attractions, Rg/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.

FIG. 2.

The total potential energy between a pair of cylinders with aspect ratio of L/D = 3 and short-ranged attractions, Rg/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.

Close modal
FIG. 3.

Example configurations of pairs of cylinders with aspect ratio of L/D = 3 and the following relative positions: (a) cos θi = cos θj = 0, and cos ψ = ±1, resulting in rh = D and Ua/ϵ = −1 (b) cos θi = cos θj = 0.2, and cos ψ = −1, resulting in rh = 1.02062D and Ua/ϵ = −0.952 (c) cos θi = cos θj = 1, resulting in rh = 3D and Ua/ϵ = −0.965 and (d) cos θi = 0.98, cos θj = 0.2, and cos ψ = 1, resulting in rh = 2.04156D and Ua/ϵ = −0.464.

FIG. 3.

Example configurations of pairs of cylinders with aspect ratio of L/D = 3 and the following relative positions: (a) cos θi = cos θj = 0, and cos ψ = ±1, resulting in rh = D and Ua/ϵ = −1 (b) cos θi = cos θj = 0.2, and cos ψ = −1, resulting in rh = 1.02062D and Ua/ϵ = −0.952 (c) cos θi = cos θj = 1, resulting in rh = 3D and Ua/ϵ = −0.965 and (d) cos θi = 0.98, cos θj = 0.2, and cos ψ = 1, resulting in rh = 2.04156D and Ua/ϵ = −0.464.

Close modal
FIG. 4.

(top) The center separation distance of a pair of cylinders at hard contact, rh, for L/D = 3 and (top left) cos ψ = −1, or (top right) cos θj = 0.2. Letter-coded labels refer to structures shown in Figure 3. (bottom) The short-range attractive interaction, Ua, at the center separation distance of a pair of cylinders at hard contact, rh, for L/D = 3 and (bottom left) cos ψ = −1, or (bottom right) cos θj = 0.2.

FIG. 4.

(top) The center separation distance of a pair of cylinders at hard contact, rh, for L/D = 3 and (top left) cos ψ = −1, or (top right) cos θj = 0.2. Letter-coded labels refer to structures shown in Figure 3. (bottom) The short-range attractive interaction, Ua, at the center separation distance of a pair of cylinders at hard contact, rh, for L/D = 3 and (bottom left) cos ψ = −1, or (bottom right) cos θj = 0.2.

Close modal
FIG. 5.

The maximum excluded volume overlap, ΔVexm, as a function of the aspect ratio, L/D. The solid line is a linear fit for aspect ratios L/D ≥ 3, with a resulting slope of 0.021 and zero intercept on the ordinate axis, and the dashed line is a linear fit for aspect ratios L/D ≤ 2, which intersects at L/D = 2.9.

FIG. 5.

The maximum excluded volume overlap, ΔVexm, as a function of the aspect ratio, L/D. The solid line is a linear fit for aspect ratios L/D ≥ 3, with a resulting slope of 0.021 and zero intercept on the ordinate axis, and the dashed line is a linear fit for aspect ratios L/D ≤ 2, which intersects at L/D = 2.9.

Close modal

The interaction potential used here requires the calculation of two computationally expensive quantities, rh and ΔVex. The latter is significantly more expensive than the former, and it is not computationally viable to calculate ΔVex on-the-fly. It follows that any effort to speed up simulations of these types of systems should focus on speeding up ΔVex. Thus, we have adopted to tabulate these interactions as described in  Appendix A. Even during the tabulation process, the time to calculate rh is negligible to that of ΔVex. The multi-dimensional tabulation of the anisotropic potential accelerated the simulation by pre-computing 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.

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, B2, for solids of revolution is given by62,63

B2(βϵ)=140πdθi0πdθj02πdψ0f(r,θi,θj,ψ,βϵ)r2sinθisinθjdr,
(5)
f(r,θi,θj,ψ,βϵ)=eβU(r,θi,θj,ψ,ϵ)1,
(6)

where β = 1/kBT 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, Uh, and the continuous potential, Ua + Ue, such that B2=B2h+B2ae.

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 L2 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 B2 with respect to the aspect ratio follows a similar trend to ΔVexm (see Figure 5).

FIG. 6.

The second virial coefficient, B2, of hard cylinders (βϵ = 0) as a function of the aspect ratio, L/D. The line is B2fit=0.773+2.22L/D+0.629(L/D)2. The error bars for this Figure, and all remaining figures, show the standard deviation of 3 independent simulations. Error bars are smaller than the size of the symbols and are not visible in this Figure.

FIG. 6.

The second virial coefficient, B2, of hard cylinders (βϵ = 0) as a function of the aspect ratio, L/D. The line is B2fit=0.773+2.22L/D+0.629(L/D)2. The error bars for this Figure, and all remaining figures, show the standard deviation of 3 independent simulations. Error bars are smaller than the size of the symbols and are not visible in this Figure.

Close modal
FIG. 7.

The second virial coefficient, B2, as a function of the interaction strength, βϵ, for integer value aspect ratios in the range L/D ∈ [1, 8].

FIG. 7.

The second virial coefficient, B2, as a function of the interaction strength, βϵ, for integer value aspect ratios in the range L/D ∈ [1, 8].

Close modal

The theta solvent condition, B2(βϵθ) = 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 B2 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.

FIG. 8.

The theta solvent condition, B2(βϵθ) = 0, is shown as a function of aspect ratio, L/D.

FIG. 8.

The theta solvent condition, B2(βϵθ) = 0, is shown as a function of aspect ratio, L/D.

Close modal
TABLE I.

Properties of cylinders as a function of the aspect ratio.

L/DVp/D3ΔVexm/D3βϵθ
0.779585 0.06712 15.96(6) 
1.55918 0.06388 15.31(6) 
2.33877 0.06378 14.45(6) 
3.11837 0.08448 16.15(8) 
3.89797 0.10535 16.8(1) 
4.67758 0.12603 17.3(1) 
5.45717 0.14671 17.4(1) 
6.23672 0.16739 17.6(2) 
L/DVp/D3ΔVexm/D3βϵθ
0.779585 0.06712 15.96(6) 
1.55918 0.06388 15.31(6) 
2.33877 0.06378 14.45(6) 
3.11837 0.08448 16.15(8) 
3.89797 0.10535 16.8(1) 
4.67758 0.12603 17.3(1) 
5.45717 0.14671 17.4(1) 
6.23672 0.16739 17.6(2) 

The theta solvent condition, B2(βϵθ) = 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, Rg. 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 

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 ensemble65,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, rh(θi, θj, ψ) < r < rc(θ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 Np/5, where Np 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 106 trials, to yield approximately 25 % acceptance of the trial move.

TABLE II.

Monte Carlo trials and relative weights for the probability of selection.

trialweight
single-particle translation or rotation 
cluster translation or rotation 1/5Np 
geometric cluster algorithm 1/Np 
βϵ change 1/100 
trialweight
single-particle translation or rotation 
cluster translation or rotation 1/5Np 
geometric cluster algorithm 1/Np 
βϵ change 1/100 

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 Np particles in a cubic domain of edge length lb/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, lnf, was initially set to unity, and was multiplied by 0.5 whenever the flatness criteria of 80 % was met. See Appendix A of Ref. 72 for implementation details of WL. The simulations were then equilibrated with 2.5 × 108 Monte Carlo trial moves, or 10 update factor reductions (i.e., lnf < 10−6). Following this equilibration, quantities of interest were averaged, and configurations were stored every 106 trials until at least 14 Wang-Landau flatness conditions were met, although some simulations reached over 37. Each simulation consisted of (3 - 50) × 109 Monte Carlo trials, depending on the conditions, and were run for well over a month of computer time. Expanded ensemble simulations were conducted for a range of βϵ in the interval [0.001ΔVexm/(4πRg3/3),0.201ΔVexm/(4πRg3/3)], using increments of δβϵ=0.004ΔVexm/(4πRg3/3). 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 lb/D = 10 and 15 to test the system size dependence of the calculated properties discussed below.

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

I(q)=Vρ(r)eiqrdr2Ω.
(7)

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 FFTW375 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.

FIG. 9.

A comparison between scattering intensity from experiment43,52 for (black dotted line) USAXS and (red dashed line) USANS at 15 °C, L/D = 4 and ϕ = 0.1125, and simulations shown by the solid lines for varying values of attractive strength, βϵ. The experimental scattering intensity was manually shifted by a constant value for comparison.

FIG. 9.

A comparison between scattering intensity from experiment43,52 for (black dotted line) USAXS and (red dashed line) USANS at 15 °C, L/D = 4 and ϕ = 0.1125, and simulations shown by the solid lines for varying values of attractive strength, βϵ. The experimental scattering intensity was manually shifted by a constant value for comparison.

Close modal

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π/lb, where lb 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, βϵ, Seff(4π/lb) is less than unity, but increases with increasing βϵ. The value of βϵ at which Seff(4π/lb) = 1 tends to increase with increasing packing fraction, ϕ. The increase of Seff(4π/lb) 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 Seff results for cylinders with other aspect ratios are summarized in Figure 12, where we plot the attraction strength corresponding to Seff(4π/lb) = 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.

FIG. 10.

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. 10.

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, ϕ.

Close modal
FIG. 11.

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

FIG. 11.

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

Close modal
FIG. 12.

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

FIG. 12.

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

Close modal

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, nc, 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% probability40 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 ⟨nc⟩ = 2.4,82 which has been shown previously to indicate gelation and dynamical arrest in isotropic model systems.39,42 The coordination number, nc, was computed as the average number of particles, denoted by i, whose exclusion volume overlaps with a central common particle, denoted by ji (e.g., when rh(θi, θj, ψ) < r < rc(θ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 nc ⪅ 2.4, the coordination number tends to increase with volume fraction; however, this trend reverses when nc ⪆ 3. Figure 14 shows the distributions of nc at conditions where the average, ⟨nc⟩ ≈ 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 nc = 2.4, precisely. An example configuration of the cylinders with ⟨nc⟩ = 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 ⟨nc⟩ = 2.4, increases with increasing aspect ratio, L/D. In addition, the attractive strength corresponding to ⟨nc⟩ = 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.

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.

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.

Close modal
FIG. 14.

The probability distribution of coordination number, nc, at the value of βϵ nearest an average coordination of 2.4, for (blue square) L/D = 4, ϕ = 0.1125, βϵ = 15.44, where nc = 2.86 (orange x) L/D = 4, ϕ = 0.2125, βϵ = 12.92, where nc = 2.6 and (green +) L/D = 8, ϕ = 0.1125, βϵ = 20.6, where nc = 2.58.

FIG. 14.

The probability distribution of coordination number, nc, at the value of βϵ nearest an average coordination of 2.4, for (blue square) L/D = 4, ϕ = 0.1125, βϵ = 15.44, where nc = 2.86 (orange x) L/D = 4, ϕ = 0.2125, βϵ = 12.92, where nc = 2.6 and (green +) L/D = 8, ϕ = 0.1125, βϵ = 20.6, where nc = 2.58.

Close modal
FIG. 15.

Cylinders colored on a red-green-blue gradient according to the number of cylinders in a cluster with L/D = 4, ϕ = 0.1125 and βϵ = 15.44. The average coordination number is 2.4. The periodic boundaries are shown in gray.

FIG. 15.

Cylinders colored on a red-green-blue gradient according to the number of cylinders in a cluster with L/D = 4, ϕ = 0.1125 and βϵ = 15.44. The average coordination number is 2.4. The periodic boundaries are shown in gray.

Close modal
FIG. 16.

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. 16.

The loci of points where the coordination number is 2.4, for a given interaction strength, βϵ, particle volume fraction, ϕ, and aspect ratio, L/D.

Close modal

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, rh(θi, θj, ψ) < r < rc(θ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.

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.

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.

Close modal
FIG. 18.

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. 18.

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.

Close modal

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, rh(θi, θj, ψ) < r < rc(θ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.

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. 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.

Close modal

In addition, the nematic order parameter was computed as the largest eigenvalue of the tensor Q = (2N)−1N(3uiuiI).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.

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.

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.

Close modal

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, B22 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 end-end 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.

H.W.H. acknowledges William P. Krekelberg for discussions about the tabulated potentiel and Yun Liu for the polydispersity smearing in  Appendix E. R. P. M. acknowledges support of the National Institute of Standards and Technology, U.S. Department of Commerce, in providing the neutron research facilities under agreements 70NANB10H256 and 70NANB12H239. Access to the BT5-USANS beamline at NIST was provided by the Center for High Resolution Neutron Scattering under Agreement No. DMR1508249. X-ray scattering was collected at the 9-ID-C USAXS beamline at the Advanced Photon Source, a U.S. DOE Office of Science User Facility operated by Argonne National Laboratory under Contract No. DEAC02-06CH11357. This work benefited from the use of the SasView application, originally developed under NSF award DMR-0520547. SasView contains code developed with funding from the European Union’s Horizon 2020 research and innovation programme under the SINE2020 project, grant agreement No 654000.

The tabular potential was constructed to obtain, rh(θi, θj, ψ), rc(θi, θj, ψ) and Ua(r, θi, θj, ψ) efficiently during the course of the simulation, as described in Section II, where rc 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 nz + 1 values of z=rrhrcrh, 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 Ua(rh)/ϵ was 4 × 10−5, and these squared differences did not vary greatly with orientation. During the course of the simulations, values of rh, rc and Ua are obtained for arbitrary orientation by multi-dimensional linear interpolation. Example two-dimensional slices of the tabular potential are shown in Figure 4. In addition, with nz = 4, the tabular potential stores in less than 1 gigabyte of memory. This smaller value of nz 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.

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, rh(θ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 (ns + 1, ns + 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, ns = 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 (ax = ay = az = 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 ns. Convergence was assessed by the quantity σr=iN=57046(1/N)(rh(i)1)2, which is the standard deviation from the theoretical result. As shown in Figure 21, the contact function calculation with ns = 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, ax = ay = az = 0.5 and ϵ1 = ϵ2 = 1 averaged over 57 046 orientations.

FIG. 21.

Convergence of the number of surface points for a sphere, ax = ay = az = 0.5 and ϵ1 = ϵ2 = 1 averaged over 57 046 orientations.

Close modal

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 dv. 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 nt random test points. A test point is inside a shape described by Equation (1) if the “inside outside” function, which is the left hand side of Equation (1), is less than or equal to unity. For each test point, if the point is inside both shapes, then there is overlap at that point. The volume of overlap within the voxel is determined by the number of test points which are inside the overlap. For each voxel that is computed by Monte Carlo integration and found to have overlaps, all neighbors of that voxel are then considered for Monte Carlo integration recursively. This algorithm applies to pairs of concave shapes which are ensured to have a single contiguous region of overlap.

The parameters for the overlap volume are now compared with theory for the spherical case. The accuracy of the overlap volume is determined by the average density of test points, here given by nt/dv3. The number of test points in a voxel, nt, is chosen to be the smallest value with no appreciable overhead cost for the flood-fill algorithm. For example, if nt is too small (e.g., nt = 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 nt is too large, for fixed point density, then the flood fill algorithm is not fully utilized. The side length of the voxel, dv, 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 Vo = π(2Rd)2(d + 4R)/12 for two equal-sized spheres of radius R = 0.54D and separation distance d = D. Thus, convergence may be assessed by the quantity σV=iN=57046(1/N)((V(i)Vo)/Vo)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, ns. This plateau occurs when the relative precision of the contact function due to ns becomes higher than the precision due to dv. For the simulations of the cylinders, dv = 0.05D and nt = 512, which has a similar point density as dv = 0.03D and nt = 128 for the sphere case, where overlap volumes were computed within an average precision lower than 0.01 percent.

FIG. 22.

Convergence of the voxel side length, dv with nt = 128 for a sphere, ax = ay = az = 0.54 and ϵ1 = ϵ2 = 1 for (red +) ns = 50 and (black x) ns = 100.

FIG. 22.

Convergence of the voxel side length, dv with nt = 128 for a sphere, ax = ay = az = 0.54 and ϵ1 = ϵ2 = 1 for (red +) ns = 50 and (black x) ns = 100.

Close modal

The volume of a single particle, Vp, was computed by the same algorithm as described here, with nt = 105 and dv = 0.06D. The resulting values of Vp, reported in Table I were utilized to determine volume fractions. Note that, with ϵ1 = 0.1, the volumes of the cylinders are within 0.74 % of the expected cylindrical volume of Vp/D3=LDπ4.

The tabulated interaction described in  Appendix A can be used directly to calculate the second virial coefficient, B2, 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πfsinθdθ=20π/2fsinθdθ,
(D1)
02πfdψ=20πfdψ.
(D2)

The tabular potential described in  Appendix A may be conveniently integrated according to this equation with some changes of variables. In particular,

0π/2fsinθdθ=01fd(cosθ),
(D3)
rhrcf(r,)r2dr=01f(z,)r2(rcrh)dz,
(D4)

where z=rrhrcrh and rc 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, Uh, and the continuous potential, Ua + Ue, such that B2=B2h+B2ae. For the hard particle contribution,

0fh(r,θi,θj,ψ,βϵ)r2dr=13rh(θi,θj,ψ)3,
(D5)

where fh=eβUh1. Thus, the second virial coefficient contribution from the hard particle potential, Bh, is obtained by substitution of Equations (6), (7), (D1)–(D3), and (D5) into Equation (5),

B2h=2301d(cosθi)01d(cosθj)0πrh(θi,θj,ψ)3dψ.
(D6)

The second virial coefficient contribution from the continuous interaction, B2ae, is obtained by substitution of Equations (6), (7), (D1)–(D4) into Equation (5),

B2ae(βϵ)=201d(cosθi)01d(cosθj)0πdψ01f(z,θi,θj,ψ,βϵ)r2(rcrh)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 nez = 100 is sufficient, the second virial coefficient was computed for a number of interpolated values in z, nez, as shown in Figure 23. Direct integration of the table is expected to be more efficient than Monte Carlo integration (e.g. Ref. 48) because a very small portion of configuration space contributes significantly, especially at relatively high values of βϵ and L/D.

FIG. 23.

The second virial coefficient of a cylinder with L/D = 3 as a function of interaction strength, βϵ, for the following values of nez: (dashed blue line) 1, (dotted green line) 10, (solid red line) 100 and (black points) 1000.

FIG. 23.

The second virial coefficient of a cylinder with L/D = 3 as a function of interaction strength, βϵ, for the following values of nez: (dashed blue line) 1, (dotted green line) 10, (solid red line) 100 and (black points) 1000.

Close modal

For slit aperture smearing,

Is(qi)j=iN1qj+12qi2qj2qi2I(qj)
(E1)
j=iN1WijI(qj)
(E2)

where Wij = 0 when qj+1qi > Δq = 0.117Å−1. See the SASVIEW 4.1.2 manual77 for derivation of the above equation in the smearing functions section. Using a diameter of 2610Å for comparison to the experimental system,43 ΔqD = 305.37.

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

Is(qi)=012πσq2exp(DiDj)22σq2I(qj)dDj,
(E3)

where qjDj = qiD. In order to implement this polydispersity smearing on the simulation data, the integration variable is changed from Dj to qj 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 = ΔqDj, where Δq = 0.135. Note that we normalize by Dj instead of D so that the discretized Gaussian kernel width does not depend on qj. The working equation is given by

Is(qi)j=iN1exp(qiqj)22(qiΔq)2qiqjI(qj).
(E4)

In order to apply the smearing kernel at the minimum and maximum q values, I(qj) was extended assuming it was constant.

1.
S. C.
Glotzer
and
M. J.
Solomon
,
Nat. Mater.
6
,
557
(
2007
).
2.
A. V.
Pethukhov
,
R.
Tuinier
, and
G. J.
Vroege
,
Curr. Opin. Colloid Interface Sci.
30
,
54
(
2017
).
3.
E.
Bianchi
,
J.
Largo
,
P.
Tartaglia
,
E.
Zaccarelli
, and
F.
Sciortino
,
Phys. Rev. Lett.
97
,
168301
(
2006
).
4.
P.
Akcora
,
H.
Liu
,
S. K.
Kumar
,
J.
Moll
,
Y.
Li
,
B. C.
Benicewicz
,
L. S.
Schadler
,
D.
Acehan
,
A. Z.
Panagiotopoulos
,
V.
Pryamtisyn
,
V.
Ganesan
,
J.
Ilavsky
,
P.
Thiyagarajan
,
R. H.
Colby
, and
J. F.
Douglas
,
Nat. Mater.
8
,
354
(
2009
).
5.
K.
Bian
,
J. J.
Choi
,
A.
Kaushik
,
P.
Clancy
,
D.-M.
Smilgies
, and
T.
Hanrath
,
ACS Nano
5
,
2815
(
2011
).
6.
M. N.
O’Brien
,
M.
Girard
,
H.-X.
Lin
,
J. A.
Millan
,
M.
Olvera de la Cruz
,
B.
Lee
, and
C. A.
Mirkin
,
Proc. Natl. Acad. Sci. USA
113
,
10485
(
2016
).
7.
W.
Liu
,
N. A.
Mahynski
,
O.
Gang
,
A. Z.
Panagiotopoulos
, and
S. K.
Kumar
,
ACS Nano
11
,
4950
(
2017
).
8.
Y.
Tian
,
T.
Wang
,
W.
Liu
,
H. L.
Xin
,
H.
Li
,
Y.
Ke
,
W. M.
Shih
, and
O.
Gang
,
Nat. Nanotechnol.
10
,
637
(
2015
).
9.
Y.
Tian
,
Y.
Zhang
,
T.
Wang
,
H. L.
Xin
,
H.
Li
, and
O.
Gang
,
Nat. Mater.
15
,
654
(
2016
).
10.
M.
Adams
,
Z.
Dogic
,
S. L.
Keller
, and
S.
Fraden
,
Nature
393
,
349
(
1998
).
11.
Z.
Dogic
and
S.
Fraden
,
Curr. Opin. Colloid Interface Sci.
11
,
47
(
2006
).
12.
L.
Onsager
,
Ann. N. Y. Acad. Sci.
51
,
627
(
1949
).
13.
R.
Zwanzig
,
J. Chem. Phys.
39
,
1714
(
1963
).
15.
A.
Cohen
,
S.
Dorosz
,
A.
Schofield
,
T.
Schilling
, and
E.
Sloutskin
,
Phys. Rev. Lett.
116
,
098001
(
2016
).
16.
S.
Asakura
and
F.
Oosawa
,
J Chem. Phys.
22
,
1255
(
1954
).
17.
S.
Jungblut
,
K.
Binder
, and
T.
Schilling
,
J. Phys.: Condens. Matter
20
,
404223
(
2008
).
18.
D.
Frenkel
,
H. N. W.
Lekkerkerker
, and
A.
Stroobants
,
Nature
332
,
822
(
1988
).
19.
J. A. C.
Veerman
and
D.
Frenkel
,
Phys. Rev. A
41
,
3237
(
1990
).
20.
P.
Bolhuis
and
D.
Frenkel
,
J. Chem. Phys.
106
,
666
(
1997
).
21.
P. G.
Bolhuis
,
A.
Stroobants
,
D.
Frenkel
, and
H. N. W.
Lekkerkerker
,
J. Chem. Phys.
107
,
1551
(
1997
).
22.
S. V.
Savenko
and
M.
Dijkstra
,
J. Chem. Phys.
124
,
234902
(
2006
).
23.
S.
Jungblut
,
R.
Tuinier
,
K.
Binder
, and
T.
Schilling
,
J. Chem. Phys.
127
,
244909
(
2007
).
24.
S.
Jungblut
,
K.
Binder
, and
T.
Schilling
,
Comput. Phys. Commun.
179
,
13
(
2008
).
25.
A.
Widmer-Cooper
and
P.
Geissler
,
Nano Lett.
14
,
57
(
2014
).
26.
C.
Ferreiro-Córdova
and
J. S.
van Duijneveldt
,
J. Chem. Eng. Data
59
,
3055
(
2014
).
27.
D.
Marenduzzo
,
K.
Finan
, and
P. R.
Cook
,
J. Cell Biol.
175
,
681
(
2006
).
28.
S. K.
Kumar
,
V.
Ganesan
, and
R. A.
Riggleman
,
J. Chem. Phys.
147
(
2017
).
29.
H. H.
Wensink
,
G. J.
Vroege
, and
H. N. W.
Lekkerkerker
,
J. Chem. Phys.
115
,
7319
(
2001
).
30.
L.
Wu
,
A.
Malijevsky
,
G.
Jackson
,
E. A.
Mueller
, and
C.
Avendaño
,
J. Chem. Phys.
143
,
044906
(
2015
).
31.
E.
Zaccarelli
,
J. Phys.: Condens. Matter
19
,
323101
(
2007
).
32.
M. J.
Solomon
and
P. T.
Spicer
,
Soft Matter
6
,
1391
(
2010
).
33.
C. E.
Alvarez
and
S. H. L.
Klapp
,
Soft Matter
8
,
7480
(
2012
).
34.
N.
Kazem
,
C.
Majidi
, and
C. E.
Maloney
,
Soft Matter
11
,
7877
(
2015
).
35.
R.
Zhang
and
K. S.
Schweizer
,
Phys. Rev. E
83
,
060502
(
2011
).
36.
M.
Tripathy
and
K. S.
Schweizer
,
Phys. Rev. E
83
,
041406
(
2011
).
37.
M.
Tripathy
and
K. S.
Schweizer
,
Phys. Rev. E
83
,
041407
(
2011
).
38.
R.
Jadrich
and
K. S.
Schweizer
,
J. Chem. Phys.
135
,
234902
(
2011
).
39.
N. E.
Valadez-Pérez
,
Y.
Liu
,
A. P. R.
Eberle
,
N. J.
Wagner
, and
R.
Castañeda-Priego
,
Phys. Rev. E
88
,
060302
(
2013
).
40.
T.
Schilling
,
S.
Jungblut
, and
M. A.
Miller
,
Phys. Rev. Lett.
98
,
108303
(
2007
).
41.
T.
Schilling
,
S.
Dorosz
,
M.
Radu
,
M.
Mathew
,
S.
Jungblut
, and
K.
Binder
,
Eur. Phys. J. Spec. Top.
222
,
3039
(
2013
).
42.
A. P. R.
Eberle
,
N. J.
Wagner
, and
R.
Castañeda-Priego
,
Phys. Rev. Lett.
106
,
105704
(
2011
).
43.
R. P.
Murphy
,
K.
Hong
, and
N. J.
Wagner
,
Langmuir
32
,
8424
(
2016
).
44.
R. P.
Murphy
,
K.
Hong
, and
N. J.
Wagner
,
J. Colloid Interface Sci.
501
,
45
(
2017
).
45.
R.
Ni
,
A. P.
Gantapara
,
J.
de Graaf
,
R.
van Roij
, and
M.
Dijkstra
,
Soft Matter
8
,
8826
(
2012
).
46.
L.
Rossi
,
V.
Soni
,
D. J.
Ashton
,
D. J.
Pine
,
A. P.
Philipse
,
P. M.
Chaikin
,
M.
Dijkstra
,
S.
Sacanna
, and
W. T. M.
Irvine
,
Proc. Natl. Acad. Sci. USA
112
,
5286
(
2015
).
47.
J.-M.
Meijer
,
A.
Pal
,
S.
Ouhajji
,
H. N. W.
Lekkerkerker
,
A. P.
Philipse
, and
A. V.
Pethukhov
,
Nat. Commun.
8
,
14352
(
2016
).
48.
H. W.
Hatch
,
W. P.
Krekelberg
,
S. D.
Hudson
, and
V. K.
Shen
,
J. Chem. Phys.
144
,
194902
(
2016
).
49.
C.
De Michele
,
J. Comput. Phys.
229
,
3276
(
2010
).
50.
C.
De Michele
,
T.
Bellini
, and
F.
Sciortino
,
Macromolecules
45
,
1090
(
2012
).
51.
I.
Torres-Diaz
and
M. A.
Bevan
,
Langmuir
33
,
4356
(
2017
).
52.
R. P.
Murphy
,
H. W.
Hatch
,
N. A.
Mahynski
,
V. K.
Shen
, and
N. J.
Wagner
, under review (
2018
).
53.
H. N. W.
Lekkerkerker
,
Physica A
244
,
227
(
1997
).
54.
C.
Vega
and
D.
Frenkel
,
Mol. Phys.
67
,
633
(
1989
).
55.
C.
Vega
and
S.
Lago
,
Comput. Chem.
18
,
55
(
1994
).
56.
A. G.
Orellana
,
E.
Romani
, and
C. D.
Michele
,
Eur. Phys. J. E
41
,
51
(
2018
).
57.
A. H.
Barr
,
IEEE Computer Graphics and Applications
1
,
11
(
1981
).
58.
J.
Glaser
,
A. S.
Karas
, and
S. C.
Glotzer
,
J. Chem. Phys.
143
,
184110
(
2015
).
60.
M. G.
Noro
and
D.
Frenkel
,
J. Chem. Phys.
113
,
2941
(
2000
).
61.
H. W.
Hatch
,
S.-Y.
Yang
,
J.
Mittal
, and
V. K.
Shen
,
Soft Matter
12
,
4170
(
2016
).
62.
J.
Corner
,
Proc. R. Soc. A
192
,
275
(
1948
).
63.
C.
Law
,
D. J.
Ashton
,
N. B.
Wilding
, and
R. L.
Jack
,
J. Chem. Phys.
145
,
084907
(
2016
).
64.
E.
Herold
,
R.
Hellmann
, and
J.
Wagner
,
J. Chem. Phys.
147
,
204102
(
2017
).
65.
K. S.
Rane
,
V.
Kumar
, and
J. R.
Errington
,
J. Chem. Phys.
135
,
234102
(
2011
).
66.
K. S.
Rane
,
S.
Murali
, and
J. R.
Errington
,
J. Chem. Theory Comput.
9
,
2552
(
2013
).
67.
H. W.
Hatch
,
N. A.
Mahynski
, and
V. K.
Shen
,
J. Res. Natl. Inst. Stan
123
(
2018
).
68.
S.
Whitelam
and
P. L.
Geissler
,
J. Chem. Phys.
127
,
154101
(
2007
).
69.
J.
Liu
and
E.
Luijten
,
Phys. Rev. Lett.
92
,
035504
(
2004
).
70.
J.
Liu
and
E.
Luijten
,
Phys. Rev. E
71
,
066701
(
2005
).
71.
D. W.
Sinkovits
,
S. A.
Barr
, and
E.
Luijten
,
J. Chem. Phys.
136
,
144111
(
2012
).
72.
V. K.
Shen
and
D. W.
Siderius
,
J. Chem. Phys.
140
,
244106
(
2014
).
73.
K.
Schmidt-Rohr
,
J. Appl. Crystallogr.
40
,
16
(
2007
).
74.
D. G.
Greene
,
D. V.
Ferraro
,
A. M.
Lenhoff
, and
N. J.
Wagner
,
J. Appl. Crystallogr.
49
,
1734
(
2016
).
75.
M.
Frigo
and
S. G.
Johnson
,
Proc. IEEE
93
,
216
(
2005
).
76.
J. G.
Barker
and
J. S.
Pedersen
,
J. Appl. Crystallogr.
28
,
105
(
1995
).
77.
M.
Doucet
,
J. H.
Cho
,
G.
Alina
,
J.
Bakker
,
W.
Bouwman
,
P.
Butler
,
K.
Campbell
,
M.
Gonzales
,
R.
Heenan
,
A.
Jackson
,
P.
Juhas
,
S.
King
,
P.
Kienzle
,
J.
Krzywon
,
A.
Markvardsen
,
T.
Nielsen
,
L.
O’Driscoll
,
W.
Potrzebowski
,
R.
Ferraz Leal
,
T.
Richter
,
P.
Rozycko
,
T.
Snow
, and
A.
Washington
,
Zenodo
(
2017
).
78.
P. D.
Godfrin
,
N. E.
Valadez-Pérez
,
R.
Castañeda-Priego
,
N. J.
Wagner
, and
Y.
Liu
,
Soft Matter
10
,
5061
(
2014
).
79.
R. B.
Jadrich
,
J. A.
Bollinger
,
K. P.
Johnston
, and
T. M.
Truskett
,
Phys. Rev. E
91
,
042312
(
2015
).
80.
M. A.
Miller
and
D.
Frenkel
,
J. Chem. Phys.
121
,
535
(
2004
).
81.
S.
Varga
,
E.
Meneses-Júarez
, and
G.
Odriozola
,
J. Chem. Phys.
140
,
134905
(
2014
).
82.
H.
He
and
M. F.
Thorpe
,
Phys. Rev. Lett.
54
,
2107
(
1985
).
83.
J. S.
Horner
,
M. J.
Armstrong
,
N. J.
Wagner
, and
A. N.
Beris
,
J. Rheol.
62
,
577
(
2018
).