Controlling the assembly of colloidal particles into specific structures has been a long-term goal of the soft materials community. Much can be learned about the process of self-assembly by examining the early stage assembly into clusters. For the simple case of hard spheres with short-range attractions, the rigid clusters of N particles (where N is small) have been enumerated theoretically and tested experimentally. Less is known, however, about how the free energy landscapes are altered when the inter-particle potential is long-ranged. In this work, we demonstrate how adaptive biasing in molecular simulations may be used to pinpoint shifts in the stability of colloidal clusters as the inter-particle potential is varied. We also discuss the generality of our techniques and strategies for application to related molecular systems.
I. INTRODUCTION
The design of new functional materials may be separated into two complementary challenges: (1) elucidating relationships between structure, chemistry, and function through the study of existing materials and (2) the controllable assembly of materials that explore new structural and chemical spaces. Both pieces are critical: understanding the fundamental mechanism underlying the material assembly helps make fabrication and synthesis methods easier and more reliable, which enables explicit testing to solidify theoretical models of structure–function relations and enables wider and deeper ranges of applications. Numerous types of solid materials have recently been pursued for the unique functional properties enabled by their open, porous structures. For instance, metal-organic frameworks (MOFs), with their potential for high surface area and synthetically controllable porous structure, are well positioned to make an impact in separations, guest catalysis, and drug delivery.1,2 Mechanical metamaterials, such as self-assembled auxetic materials derived from an interconnected network of certain specifically shaped units, present the unique property of negative Poisson’s ratio, making them potential candidates for applications that include protective textiles and smart sensors.3–5 Related colloidal systems, with extensive applications to optical metamaterials,6–8 porous materials,9,10 and other systems of industrial relevance, have been the subject of extensive scientific investigations, particularly for the creation of well-controlled assemblies such as bespoke clusters and open crystals. Interest in the applications of these classes of materials is so high that research has pressed forward into new structural spaces without a full understanding of the nucleation and assembly processes, controlling their formation. Colloidal materials make attractive candidates for the study of these complex behaviors due to their larger size, modest energies, and slower kinetics, each of which facilitates direct observation using a microscope. Model colloidal materials, such as spheres with diameters on the length scale of the visible spectrum,11–14 have great potential for applications in the processing and diffraction of light, within optical filters, in optical switches, and as photonic band gap materials.14–18 Both types of porous materials have the potential for use as mechanical metamaterials, whose novel structures can act to disperse and dissipate energy.
For developing colloidal technologies, specific threedimensional self-assembled structures are desired. For instance, open lattices such as diamond and pyrochlore17,19 provide the most attractive photonic materials,16,17,20–23 but these large-scale three-dimensional assemblies have proven to be difficult to achieve in the laboratory; there is similarly a gap between the laboratory and the production scale for MOFs.24 Although interactions can be tuned to make colloidal diamond or pyrochlore somewhat favorable, neither of these structures is easily accessible to current methods of manufacturing because neither of them is a close-packed lattice. It has been shown that colloidal particles with interactions near to the hard-sphere limit are driven toward close-packed assemblies25–27 although multiple theoretical and simulation efforts over the past decade have been devoted to finding potentials that favor open structures in the thermodynamic equilibrium.28–32 Some interesting examples play off the conflict between short-range attraction and long-range repulsion although these tend to assemble into disordered structures in experiments.33–38 Charged particles may assemble into lattices with larger vacancies than close-packing when different-sized particles are used,39–41 while DNA-functionalized colloids have constructed a large number of different arrangements controlled by the stoichiometry and size asymmetry of multicomponent particle suspensions.42–44 Additional approaches utilize inventive fabrication methods to facilitate specific binding geometries, favoring the assembly of open lattices.19,45
It is important both for applications of colloids as model systems and in creating open structural lattices to understand and control early cluster formation. In this way, staged protocols46 assembling colloidal spheres into anisotropic objects that mimic molecular structures can be driven and analogs of reticular synthesis47,48 can be applied to the design of colloidal materials. Thus, we seek to understand simple accessible methods by which colloidal cluster structures may be tuned to favor specific arrangements. Inspired by work on small-N clusters (where N is the number of colloidal particles within a particular assembly), we examine methods by which one can drive the free energetically favored cluster assembly from octahedral (OCT) to capped trigonal bipyramid (CTBP, a polytetrahedral assembly).49–52 These structures are illustrated in Figs. 1(a) and 1(b). According to the Cambridge Cluster Database,53 the OCT structure minimizes the potential energy although CTBP has fewer symmetry groups than OCT, leading to a higher configurational entropy.50 For very short-range interactions, only pair contacts are relevant, and entropy is a major contribution to the free energy. This results in a dominant configuration of CTBP. Hence, as a function of particle interaction range relative to the particle size, there must be a transition between the CTBP and OCT morphologies. It is unclear where this transition should take place as current theoretical treatments focus on either the free energy of the short-ranged interaction limit50 or simple potential energy minima of the long-ranged interaction limit.53,54 Whether configurations are modified when varying the interaction range in the intermediate parameter space is, likewise, still unknown.
CTBP and OCT configurations are shown in panels (a) and (b), respectively. Panel (c) illustrates the variation of the Morse potential as α changes between the extremes used in this study. Panel (d) shows the switching function curve we utilize to select CTBP and OCT configurations further, plotted alongside a histogram of contact number vs distance for the perfect clusters in (a) and (b).
CTBP and OCT configurations are shown in panels (a) and (b), respectively. Panel (c) illustrates the variation of the Morse potential as α changes between the extremes used in this study. Panel (d) shows the switching function curve we utilize to select CTBP and OCT configurations further, plotted alongside a histogram of contact number vs distance for the perfect clusters in (a) and (b).
Here, we employ enhanced sampling methods to build a robust free energy landscape on which we expect to locate two minima corresponding to each metastable state. Importantly, this provides a method capable of identifying and sampling the relative weights of arbitrary cluster states in systems with a fixed number of particles and arbitrary molecular interactions, since approximate analytical theories are limited in their applicability by the need to assume constrained, approximately polynomial potential energy minima. While our work here focuses on spherical clusters with N = 6 particles, the techniques are thus readily generalizable to arbitrary shapes and numbers of particles, given that appropriate order parameters may be determined either a priori or on-the-fly using machine learning.55–58 More specifically, the free energy landscapes we measure in this work offer a design principle for colloids whereby polytetrahedral structures (such as the CTBP) may be targeted relative to close-packed structures by tuning the particle size and thus influencing the final colloidal assembly.59
II. METHODS
Ultimately, our interest is in the design of bespoke open colloidal lattices. As such, we are interested in how the molecular interactions may be tuned to favor open structures (such as networks of Bernal spirals or Kagome lattices).60–63 The ability to form these phases depends heavily on the structures favored in early aggregation.64 Inspired by the idea that manipulation of the interaction range can induce structural changes, we choose the Morse potential [see Fig. 1(c)],
a highly controllable potential whose range may be tuned relative to the particle size by adjusting the parameter α, with longer-ranged interactions having smaller α values. We choose the reduced units of energy and length by analogy with the Lennard-Jones system, setting kBT = ϵ = 1 as our energy scale and the equilibrium particle distance rc = σ = 1. D0 is chosen so that the potential well depth has a value of 10kBT; this ensures that the probability of breaking a bond is ∼e−10. This disfavors states with multiple simultaneous broken bonds and has the effect of rendering a precise definition of the transition state unimportant for the free energy computations that follow as their probability e−βF becomes negligibly small.
The two metastable structures have the same nearest neighbor pair count (12) although the number of next-nearest neighbor pairs differs. In terms of bond-lengths (here, essentially equivalent to the particle diameters), the OCT has three diagonal contacts at a distance of , while the CTBP has its two next-nearest contacts at a distance of , equal to the diagonal distance across a triangular-bipyramid, and one contact at a slightly larger distance of 5/3.49 This is shown in Table I and Figs. 1(a) and 1(b), where the nearest neighbor of a blue particle is in red and the next-nearest neighbors are in yellow. When approaching the hard-sphere limit, interactions only exist between the nearest neighbors, and according to the cluster free energy, F = U − TS, CTBP becomes the entropically favored structure by a free energy ΔF ≈ 3kBT.50 As the interaction range increases, the potential energy competes more prominently with entropy in the free energy and the more compact OCT shape becomes preferred due to its significantly lower potential energy.
Number of contacts between colloidal particles at different distances for the two N = 6 configurations discussed in the text.
Configuration . | r/r0 . | Ncontacts . |
---|---|---|
OCT | 1 | 12 |
1.414 | 3 | |
CTBP | 1 | 12 |
1.63 | 2 | |
1.67 | 1 |
Configuration . | r/r0 . | Ncontacts . |
---|---|---|
OCT | 1 | 12 |
1.414 | 3 | |
CTBP | 1 | 12 |
1.63 | 2 | |
1.67 | 1 |
Typically, for deep potential wells, a large amount of time is required by classical molecular dynamics to obtain sufficient sampling of both configurations. Therefore, we apply the adaptive biasing force (ABF)65 enhanced sampling method to compute the free energy landscape relevant for clustering. Essentially, this method computes the mean generalized force vector of the system along a finite number of collective variables (CVs) and then biases against this force so that the free energy landscape becomes “flat” along these coordinates. For more details on the algorithm, interested readers are referred to the review by Comer and co-workers.66 Useful CVs are defined as functional combinations of atomic coordinates. Once the ABF method converges upon an estimated mean force, which can be determined by monitoring CV diffusion (which should mimic a thermodynamically “flat” system along the CVs) or error thresholds in the mean force landscape, the free energy landscape can be obtained by integrating these forces over the defined CVs.
In this study, we apply ABF to two CVs that describe the N = 6 colloidal structures. The first is a pairwise function with a rational switching function kernel,
which, as designed, counts twice the number of “second-neighbor” contacts in our clusters. In this equation, N = 6 is the number of particles in the cluster, rij refers to the inter-particle distance between particles i and j, and r0 and d0 are customized parameters for the switching function. Neighbors that are selected by this function are shown in Figs. 1(a), 1(b), and 1(d).
Additionally, we use the squared radius of gyration of the cluster,
which captures the compactness of the cluster. Here, and refer to the position of each particle and the center of mass, respectively. These two CVs are sufficient to separate the two types of clusters in the N = 6 particle system, provided that the switching function parameters are set to capture the difference in second-neighbor contacts (see Table I). We choose to set d0 and r0 to 1.75 and 0.2, respectively, such that Sij is ≈0 at (filtering out second neighbors in the OCT configuration) and is ≈1 in the range [1.6, 1.9] (including neighbors in the CTBP configuration); this is illustrated in Fig. 1(d).
All simulations utilize the open-source software packages SSAGES67 (for advanced sampling) and LAMMPS68 (for molecular dynamics simulation) to compute the free energy of clusters of N = 6 colloidal spheres with Morse potential interactions of varying α. As mentioned previously, all quantities reported in this paper are non-dimensionalized using standard Lennard-Jones units based on σ = 1.0, ϵ = 1.0, and m = 1.0. We utilize a periodic cubic cell with an edge length of 8.0, which is large enough that colloids do not interact with periodic images of the cluster when fully bonded. Simulations are equilibrated at temperature T = 1.0 using a Langevin thermostat with a damping factor of 0.5 for 105 time steps, after which ABF biasing commences. In the ABF sampling, bias on S is accumulated on the interval [0.2, 7.0] with a restoring force applied outside of a buffer zone (with no bias) with a thickness of 0.1 to push the system back toward the region of interest. Similarly, bias on is acquired on the interval [0.45, 0.65] with a restoring force applied to the regions outside of its buffer zone with a thickness of 0.15. The two-dimensional region utilizes a 100 × 100 grid, and the restoring forces are generated by a half-harmonic potential with a spring constant of 100ϵ/σ2, centered at the edges of the buffer zones. For the ABF algorithm, we set the minimum count of the number of visits required to a bin in the general histogram before the full bias is activated to 200.66 Enhanced sampling calculations are performed for another 5 × 106 time steps after an initial equilibration period of 106 time steps to generate the resulting free energy maps. We found that these parameters were sufficient to explore this system’s basins and combined with multiple independent copies of the ABF simulation enabled accurate data capture and assessment of relative error in free energy measurements. Finally, our time step was set to δt = 10−3t0, and the Langevin relaxation time, which favors inertial motion at large values, is set to 0.5 (meaning that particles in our system have a self-diffusion time of 0.5 in reduced units). However, we stress that the free energy surfaces we obtain were confirmed in test simulations to be insensitive to these choices so long as the combination of time step and friction permits a proper exploration of the Morse potential minima.
III. RESULTS AND DISCUSSION
When α is large, the interparticle interaction decays quickly to zero beyond the equilibrium distance rc, and only contact interactions contribute appreciably to the potential energy. Since the numbers of nearest neighbors are the same for both CTBP and OCT, their free energy will differ primarily through the amount of entropy available, which is maximized by CTBP.50 As compact structures will be energetically favored for small α, there must be an α for which the interactions fall into the intermediate range. Note that while the energetic dependence on this parameter is relatively simple, the entropic dependence is not, as the relative vibrational modes are affected both by the increased range of interactions and by conformational limitations.
We explore the role played by the parameter α to bias the morphology toward the two structures, sweeping from α = 5 to α = 15. For each α, the generalized force is gathered during the biased sampling; then, two-dimensional probability distributions (related to the free energy through for a macrostate i) obtained as a function of and S are subsequently integrated to obtain the relative free energy difference of the two states. For each system, a set of 10 independent runs is utilized in the free energy plots. (See Fig. 2 and the supplementary material for free energy landscapes for each of the α studied, alongside error estimates computed for the surfaces.) As shown in Fig. 2, two basins are identified on each free energy surface, indicating the two metastable states of the six-particle cluster. These are identified through S and values with CTBP and OCT: OCT has both lower S and lower (≈0 and ≈0.5, respectively), while CTBP has higher CV values (S ≈ 6 and ), in accordance with the model shown in Table I. Transition states between these basins that have a single broken bond [see Fig. S2(a) in the supplementary material] are observed frequently at large α [Fig. 2(b)] and correspond to the saddle point between the OCT and CTBP basins.
Two-dimensional free energy landscape as a function of the switching function and radius of gyration. The upper panel (a) is for a long range interaction model with α = 5, and the lower panel (b) is for a short range interaction model with α = 15.
Two-dimensional free energy landscape as a function of the switching function and radius of gyration. The upper panel (a) is for a long range interaction model with α = 5, and the lower panel (b) is for a short range interaction model with α = 15.
Other apparently stable minima are observed in some systems; the minimum with S > 6 and in Fig. 2(b) is one such configuration. This basin links to a set of intermediate configurations when the CTBP cluster is fluctuating between isomers, as shown in Fig. S2(b) in the supplementary material; such states were also observed in Ref. 69. These are not taken into account in the subsequent analysis comparing CTBP and OCT, as they correspond to structures that do not maximize contact interactions and thus do not define mechanically stable clusters.50,51 The apparent free energy stability is a consequence of the large amount of entropy available to configurations with broken bonds. This minimum would be quenched in the limit where the pairwise well depth satisfies D0 ≫ kBT.
To further investigate the relationship between the favorable states and the interaction range, we examine the integrated free energy difference
as a function of α. Here, ΔF12 = F2 − F1 and Ωi refers to the regions of phase space, corresponding to the OCT (2) and CTBP (1) clusters. These are defined by the nearby regions of the free energy landscape that are within 5kBT of the minimum. Other choices can be made which will shift the quantitative answer slightly but will not affect qualitative conclusions about stability. The choice of 5kBT was determined by trial and error to be robust against the influence of higher energy metastable basins in large α surfaces not corresponding to stable cluster states, such as those observed in Fig. 2(b). The results shown in Fig. 3 demonstrate that the CTBP is confirmed distinctly to be more favorable in shorter range (large α) interactions, while OCT becomes more favorable as the interaction range increases.
Free energy difference between CTBP (1) and OCT (2) states (F2 − F1) as a function of the α value in the Morse potential. Each ΔF value is averaged over 10 data points. The green dotted line is at ΔF12 = 0 and defines a crossover between the two preferred states. Bars denote the standard deviation of the sample data points.
Free energy difference between CTBP (1) and OCT (2) states (F2 − F1) as a function of the α value in the Morse potential. Each ΔF value is averaged over 10 data points. The green dotted line is at ΔF12 = 0 and defines a crossover between the two preferred states. Bars denote the standard deviation of the sample data points.
The range from α = 5 to α = 15 was selected due to the strong similarity of α = 5 to the atomic Lennard-Jones interaction, indicating the particles that make up the clusters have atomic dimensions. As α increases, the interaction range decreases, which can be viewed as an increase in the corresponding particle size. A Morse potential with α = 14 is reasonable to describe fullerene molecules, for instance.53 As α increases to values near 40, the Morse potential closely mimics the size and interaction range ratios present in colloidal van der Waals interactions, thus making the Morse potential there a good model for attractive colloids at the μm level.
In practice, colloidal interactions can be described as a summation of van der Waals attractions and electrical double-layer repulsion according to the DLVO theory;70 such interactions could be re-illustrated with our Morse potential model with specific α values. In this study, at α ≈ 10, where the crossover between OCT and CTBP is located, a Morse potential mimicking van der Waals interactions would define the interaction of particles that are smaller than fullerenes but larger than atomic sizes. If we would like to achieve such behavior for colloidal clusters, alternative interactions, such as a combination of repulsive surface grafts and depletion due to non-adsorbing polymer, would have to be used;71–73 although the Morse potential is not perfectly mapped onto these interactions, it should be possible to tune the second virial coefficient with an interaction range in the vicinity of 10% [see Fig. 1(c)] to access this transition.74
IV. CONCLUSIONS
We mapped out the two-dimensional free energy surface of six-particle colloid clusters with a variable interaction range by applying biased sampling to next-nearest neighbor configurations and the radius of gyration. We found that the two basins on the free energy surface correspond to the two equilibrium structures of six-particle colloids, polytetrahedral (CTBP) and octahedral (OCT). As colloidal clusters with the number between N = 1 and N = 10 have been enumerated,49–52 the workflow described here can enable the systematic, swift determination of the most-likely assemblies through subtle modifications of these structural CVs. By integrating the configurational probability distribution over these basins, we found that the favored structure differs at different potential ranges and provides a possible route to the design of clusters with a desired geometry by modifying the particle size and the range of dominant interactions.
The utility of the calculations presented here is twofold. First, we demonstrated the utility of the interaction range in tipping the preferences for self-assembly of colloidal clusters. Second, we demonstrated how molecular simulations that comprehensively obtain cluster assembly landscapes can be helpful for targeting the assembly of stable and metastable cluster morphologies. Extension of our results to engineering other types of interactions, including optimizing parameters for patchy particles and open-lattice precursors, is straightforward. This enables the design and optimization of processes to create functional materials from bespoke particle assemblies.
SUPPLEMENTARY MATERIAL
See the supplementary material for free energy surfaces (and error estimates) for each of the α parameters studied and illustrations of observed cluster transition states.
ACKNOWLEDGMENTS
We thank Hythem Sidky and Vikramjit Singh Rathee for fruitful discussions and technical help. S.H., M.J.Q., and J.K.W. were supported in this work by the National Science Foundation Award No. DMR-1751988. The code SSAGES utilized in this work was developed by the MICCoM Center at Argonne National Laboratory under a grant from the Department of Energy, Basic Energy Sciences.