DOI: 10.1039/C7ME00081B
(Paper)
Mol. Syst. Des. Eng., 2018, Advance Article

Huikuan Chao and
Robert A. Riggleman*

Department of Chemical and Biomolecular Engineering, University of Pennsylvania, Philadelphia, PA, USA. E-mail: rrig@seas.upenn.edu

Received
29th August 2017
, Accepted 15th November 2017

First published on 15th November 2017

Two dimensional nanoparticle lattices can exhibit unique optical, electrical, and chemical properties giving rise to emerging applications for photovoltaic conversion, electronics and optical devices. In many applications, it is useful to be able to control the particle spacing, the crystal lattice formed, and the local composition of the lattice by co-locating nanoparticles of varying chemistry. However, control over all of these variables requires exquisite control over the interparticle interactions, and a large number of degrees of freedom affect the interactions. Achieving a particular structure by design requires solving the inverse-design problem, where one must optimize the chemistry to meet the structure or property that is desired. In recent years, a variety of examples have shown that one can finely control the interactions between nanoparticles through the use of polymers grafted onto the nanoparticle surface and by controlling the grafting density and the distribution of molecular weights on the nanoparticle surface. In this work, we take the first steps on solving the inverse design problem using an approach that explicitly accounts for the chemistry on the surfaces of the particles. Using two-dimensional hybrid particle/field theory calculations and an evolutionary design strategy, we design polymer grafted nanoparticles that self-assemble into targeted square, honeycomb, and kagome lattices. We optimize both the length and grafting density of the polymers grafted to the nanoparticles, and we show that our design strategies are stable over a range of nanoparticle concentrations. Finally, we show that three-body interactions are critical for stabilizing targeted structures.

## Design, System, ApplicationIn this work, we use an evolutionary design strategy to design polymer-grafted nanoparticles for targeted self-assembly. While the work is fundamental in nature, it could have large impacts in the self-assembly of nano- and colloidal particles for optical applications. There is a clear molecular design component to this work since we apply the evolutionary optimization strategy to the chemistry of the molecules on the nanoparticles' surfaces to achieve a targeted structure. |

One popular strategy to tackle the design challenge is to use inverse design approaches.^{6–14} In this approach, a specific ordered structure is designated as the target. Then an optimization approach is developed to vary the intermolecular potential between the particles such that the particles with the optimized potential will obtain the targeted structure at equilibrium. Pioneering studies by Torquato, Stillinger, and co-workers established an optimization framework that searches for isotropic, pair-wise potentials where the objective function is defined as the difference between the ground state potential energies for the target structure and other candidate structures.^{6} Their approach was successful in searching for isotropic pair-wised potentials for particles self-assembled into two dimensional and three dimensional lattices.^{7,8} More recently, Truskett and coworkers have utilized simulated annealing techniques to optimize the potential such that the stable density range over which the target structure exhibits the lowest ground state potential energy against other candidate structures becomes maximized,^{9} and the form of the potential was restricted to purely-repulsive, convex forms. The method has been demonstrated in potentials that stabilize both two and three dimensional particle lattices against finite variations in the particle density.^{13,14} However, one challenge that remains outstanding in this field is providing direction for achieving a specific interaction between particles in an experimentally-accessible manner. In other words, it remains unknown what specific chemistry should be used to induce a precise interparticle interaction.

Grafting polymers to the surfaces of particles has long been used to stabilize colloidal and nano-particles against aggregation by effectively rendering the interparticle interactions repulsive.^{15} For nanoparticles with monodisperse grafting, previous studies have shown that at high grafting density the potential of mean force (PMF) is purely repulsive when long grafted polymers are wetted by short matrix polymers and screen depletion forces between the nanoparticles. In contrast, the PMF can become attractive when the grafted polymer is shorter than the matrix polymer and dewets the brush layer.^{16,17} When the grafting density becomes sparse, studies show that the interactions between the nanoparticles become anisotropic due to the open surface areas uncovered from the grafted polymer.^{18,19} The anisotropy can be further tuned by changing the grafted polymer length and density giving rise to a series of nanoparticle assemblies.^{20–22} Jayaraman and coworkers have shown how polydispersity in the molecular weight of the grafted polymer would affect the PMF between polymer-grafted nanoparticles in a polymer matrix that is chemically identical to the brush. By combining polymer reference interaction site model (PRISM) theory and Monte Carlo simulation, they found at high grafting density condition the short range repulsion is enhanced by the crowing of monomers from both long and short grafted polymers while the attractive well in the intermediate range can be suppressed by the steric repulsion between the long grafted polymers.^{23,24} The enhanced repulsion between nanoparticles in this grafting strategy has been demonstrated by experimental studies, where inorganic nanoparticles with bimodal PDMS grafts show strong miscibility in polymer matrices with high molecular weight PDMS.^{25,26} Varying the graft architecture and using either semi-flexible polymers, mixed polymer brushes, Janus-grafted polymers, or block polymer brushes are other alternative routes for tuning the interparticle interactions.^{17,27–31}

Evolutionary design optimization schemes, in particular the covariance matrix adaptive evolutionary strategy (CMA-ES),^{32,33} have recently emerged as a route for solving complex optimization problems in soft matter research. CMA-ES has been used in both the design of granular particles with tailored mechanical properties^{10} and in the guided pattern design of self-assembling block copolymers.^{11,12} In both applications, the method has been demonstrated to have efficient convergence despite the complicated function that is optimized. In this method, variables (denoted as a vector, ) to be optimized are updated iteratively via a stochastic process mimicking natural evolution. In an optimization trial denoted as j, mutations of variable _{j} are generated. The quality of each mutation is quantified by a predefined fitness function G_{j}(). Then, the values of _{j+1} for the next generation are determined both upon the covariance matrix of the fitness function of the current mutations and the evolution path of the previous generations. This process is then iterated until convergence, which is typically taken when the fitness function becomes smaller than some threshold value.

In the manuscript, we employed CMA-ES as the optimization method to establish an inverse design framework of grafted nanoparticles based upon a polymer field theory model. In our approach, rather than tuning the non-bonded interactions between the nanoparticles to achieve a particular self-assembled structure, we change the length and number of grafted polymers on the surface of the particle, which in turn modifies the interactions between the nanoparticles. In this first demonstration, the CMA-ES method is applied to a series of two-dimensional lattices by searching for optimized grafting designs that minimize the free energy of square, honeycomb, or kagome lattice unit cells. These optimized designs are then validated using hybrid particle field theory (HPF) simulations where a number of nanoparticles with the optimized designs are sampled and shown to self-assemble into the targeted lattice. After validating our approach, we analyze the role that each grafted chain length plays in the assembly, and we demonstrate that the many-body interactions play a key role in the assembly of the nanoparticles into the targeted structures.

(1) |

The mass of each monomer is distributed about its center using a unit Gaussian function, , where a is the smearing size of the Gaussian function. In a similar manner, the mass of a single NP is distributed about its center according to a complimentary errors function, , where the density decays from bulk value, ρ_{0}, within the particle radius, R_{P}, to zero away from the particle surface over a width controlled by ξ. The grafting sites are distributed homogeneously over the nanoparticle, and the distribution of grafting sites is given by

(2) |

(3) |

The NPs are immersed in an implicit θ-solvent to the polymer monomers, which allows us to neglect the pair-wise non-bonded interactions between monomers, and the only non-bonded interactions to be considered are the interparticle interactions and the monomer-particle interactions. We assume that the only interactions between the polymer monomers and the nanoparticles are the excluded volume interactions, which is defined as a convolution between the density distribution functions of the two species,^{35,36}

(4) |

(5) |

The particle based partition function of the grafted nanoparticle system has the form,

(6) |

To develop the hybrid particle-field theory model, we closely follow the particle-to-field transformation^{37} to decouple pair-wise interactions between monomers and NPs and transform eqn (8) into the field-based form as

(7) |

Here is the effective Hamiltonian with the form,

(8) |

(9) |

(10) |

Here, _{g}(r;[w]) is the density operator of the grafted chains, which can be evaluated as the following,

(12) |

(13) |

q_{gi}(r,0) = e^{−w(r)},
| (14) |

(15) |

We note that the NP coordinates are retained through the particle-to-field transformation and remain explicit in the effective Hamiltonian. This makes the gNP model in the family of HPF simulation approaches^{34} and we denote it the hSCFT model of polymer-grafted nanoparticles. We take b as the unit length in our simulations, and all energies are scaled by k_{b}T.

During the optimization phase where we implement the CMA-ES algorithm (described below), we perform our calculations in a unit cell with fixed particles that are not allowed to move. Thus, evaluating eqn (6) for this unit cell gives us direct access to the free energy of our system with a fixed nanoparticle lattice.

(16) |

Thus, the total inter-particle core potential is

(17) |

The equation of motion for the ith NP is given as

(18) |

(19) |

The parameter ζ in the above equations is set to unity and a time step of Δt = 2 × 10^{−3} is used here. We would like to note that by tracking the inter-particle forces from U_{nn} we do not observe overlap between nanoparticle cores in all the HPF simulations performed here. This is likely due to the repulsion provided by the grafted polymer chains. We note that the force, , experienced by each particle contains contributions from both inter-particle and particle-monomer interactions due to chain grafting.^{34}

In our validation simulations, unlike standard self-consistent field theory calculations, since the nanoparticle positions are evolving direct access to the exact free energy of the system is no long available. Therefore, we resort to an approximation of the free energy. The expression for [w, ρ] given above in eqn (8) contains contributions from the enthalpic interactions between the particles and the polymers as well as the entropy of the grafted chains. The entropy associated with the vibrations of the nanoparticles around their equilibrium positions is absent, and this entropic contribution is approximated by using the classical form of the quasi-harmonic approximation

(20) |

D_{ab} = 〈u_{a}(t)u_{b}(t)〉,
| (21) |

v_{i} = λ_{i}^{−1/2},
| (22) |

(23) |

With the goal of maintaining as simple of a grafted polymer design as possible, we restrict the design to bidisperse grafted nanoparticles; _{j} consists of two chain lengths and two grafting densities. The population size (number of mutations for the set of _{j}) is taken as seven times the number of variables we are optimizing, giving a parameter population size of 28. This choice is based on suggestions from previous studies^{11,32} and our own tests for convergence. The larger one of the two chains lengths is initialized randomly from a Gaussian with a mean of (2R_{P})^{2} (so the average end-to-end distance is on the order of the particle diameter) and a variance of (R_{P}/2)^{2}. The smaller chain length is initialized to one half the longer chain length, and both grafting densities are initialized to 2πR_{P}.

At each iteration in the optimization procedure, a number of mutations of the current generation of _{j} are made, we evaluate our fitness function G_{j}, and we call the CMA-ES routine to update the parameters that comprise _{a,j}. The implementation here uses the CMA-ES developed by N. Hansen and co-workers,^{32,33} which is freely available as a Python package. In each trial, the CMA-ES routine continues until the parameters no longer change and the change of G_{j} within a tolerance of ΔG_{j}/n_{ρ}n = 10^{−4}k_{B}T or the global fitness function G_{j} is less than −5k_{B}T per nanoparticle in the unit cell. The total optimization process is then repeated to convergence for fifty different trials, j ∈ [1, 50], and the solution with the lowest G_{j} is taken as the converged solution and tested using the Brownian dynamics HPF described above.

We next tested whether the converged solutions were able to spontaneously form the targeted structure by performing HPF simulations with 100 NPs with the grafting design from the solutions. Here, the assembled NP configurations from two types of HPF simulations are presented. The first type coming from simulations with initial NP configurations biased to square lattices (Fig. 2a), exhibits a defect-free square lattice pattern for the NPs with design SQ1 over all the entire range of ρ. For the other two sample converged designs SQ2 and SQ3, the assembled structure become less ordered for the highest density test, where ρ = 0.042. Similar results can be observed in the results from simulations with random initial NP configurations (Fig. 2b), though generally speaking more defects are present. The assembled structures from NPs with the design from SQ1 also exhibit square lattice pattern at ρ = 0.042, although the defects in the form of grain boundaries begin to appear in the pattern. However, for the NPs with the designs from SQ2 and SQ3, the square lattice pattern can only be observed in the lowest ρ value at 0.034. It is interesting that the defective systems, SQ2 and SQ3, adopt local motifs with five-particle symmetries, though long-range order is clearly lacking.

Fig. 2 Snapshots from HPF simulations on the three selected designs in Fig. 1a where the simulations were initialized with (a) biased and (b) random initial configurations, respectively. The grafted chains are not shown for clarity. |

Fig. 3 Optimized design for NPs forming SQ, HC, and KG lattices over the range, ρ ∈ [0.034, 0.042] (main figure) and ρ = 0.038 (inset). |

The above designs are tested using HPF simulations with multiple NPs in the system. The assembled configurations shown in Fig. 4a were produced from HPF simulations with initial configuration biased to the target lattice. The NPs using SQ and HC design are found to assemble into defect-free square and honeycomb lattices at all three values of ρ, while defects are observed in the KG lattice. If the initial condition in HPF simulations is changed to a random configuration, as presented in Fig. 4b the NPs using the SQ design still predominantly show a square lattice, though defects and grain boundaries do emerge. However, the NPs using the HC design exhibit localized defects at small and intermediate values of ρ. For the NPs using the KG design, we do not observe well-ordered KG structures at all the values of ρ. Both the KG and HC lattices do provide an open structure, which could be advantageous for 2D assemblies.

To quantify the ordering of the above shown assembled structures, we construct a lattice order parameter S_{L} based on the distribution of angles formed by three adjacent NPs in the structures. Briefly, we calculate the ratio of the number of angles ϕ distributed within a small range (δϕ = ±5 degree) about the characteristic angles in the lattices relative to the total number of angles; the ratio of the number centered about their characteristic angles to the total defines S_{L}. For example, square lattice has two characteristic angles at 90 and 180 degrees such that angles distributed in [85 95] and [175 185] will be accounted to the numerator of the ratio. This ratio has the limits of 1 when the NPs assemble into a perfect lattice and 2δϕ/180 when the angle distribution is completely random. The analyzed S_{L} values from the above structures are presented in Fig. 5. Amongst all the types of NP designs, the structure coming from biased initial configurations exhibit higher values of S_{L} over the range of ρ compared to these from the random initial configurations. This trend is also observed for the structures coming from NPs using designs optimized just at ρ = 0.038. For NPs with SQ design, the S_{L} in the assembled structure coming from random initial configurations reaches highest value at ρ = 0.038, in the HC design S_{L} decreases as ρ increases, and for the KG design S_{L} is approximately constant.

These results demonstrate that the biased states exhibit fewer defects, and we next provide evidence that these configurations are lower free energy states. To approach this, the free energy per NP of the assemble structures is quantified from three contributions: the entropic energy from the grafted chains, H_{g}, the enthalpic energy from NP-grafting chain interactions, H_{w}, and the vibrational free energy from NP configurations, H_{q}. Fig. 6 plots the per NP free energy, F, of the assembled structures as a function of ρ. We observe that for all lattices considered here the free energy has a smaller value from the simulations with biased initial configuration over the range of ρ. As shown in the insets to the figures here, we also find that the majority of free energy contribution arises from the entropy of the grafted chain, H_{g}. Based on this set of results we would argue that the assembled structures from the biased initial configurations are energetically more favored states, and the entropy of the grafted polymers plays a dominant role in driving the system into each configuration.

We next examine the NP–NP interactions through the two-body potential of mean force (w_{2}(r)) between isolated NPs. As shown in Fig. 8, all the potentials of mean force exhibit a purely convex, repulsive nature when the NP–NP distance, r, is below the equilibrium distance in the lattice, r_{0}. Interestingly, we observe that w_{2}(r) remains finite even at distances where r is approximately two times r_{0}. The importance of these long-range effects becomes even more pronounced when we investigate the multi-body interactions. To illustrate this point, the total three-body potential of mean force, w_{3}(r), is calculated. The results show that w_{3}(r) exhibits a finite strength beyond 2r_{0} for each optimized design. Specifically, for the w_{3}(r) from the HC and the KG designs, the interaction range can exceed 4r_{0} (Fig. 9a). We further isolate the three body component of the potential by subtracting off the two-body potential between the three particles involved, u_{3}(r) = w_{3}(r) − 2w_{2}(r) − w_{2}(r_{0}); the result is shown in Fig. 9b. It is interesting to observe that the potential shows a non-monotonic change when r increases across the equilibrium distance r_{0}. In particular, u_{3}(r) from SQ design exhibits a peak intensity right at the position where the three particles would form an equilateral triangle. This local maximum in the potential indicates that the three body interaction destabilizes triangle structure in the assembly of NPs, which likely leads to the stabilization of the square lattice.

We take the methods and results presented here as a promising first step towards designing the chemistry of nanoparticles that could be designed experimentally for targeted self-assembly, but there are several assumptions that need to be relaxed in future work. The first one is the definition of uniform grafting density in the hSCFT NP model, which might become an unphysical assumption when the grafting density is small and the discrete nature of the grafting sites becomes important. It is well-known that grafted particles at low graft densities exhibit anisotropic interactions and assemblies.^{18} One strategy for addressing this issue can be to adopt the dynamic mean field framework introduced by the authors^{39} where the grafting sites on the NP surface are treated explicitly, though this loses direct access to the free energy. In addition, the polymer chain conformations are treated entirely at the mean-field level in a θ solvent, and testing the fluctuation effects and the influence of solvent quality and the effects of other enthalpic interactions remains an outstanding issue.

- H. A. Atwater and A. Polman, Nat. Mater., 2010, 9, 205–213 CrossRef CAS PubMed.
- M. J. Hore and R. J. Composto, Curr. Opin. Chem. Eng., 2013, 2, 95–102 CrossRef.
- M. J. A. Hore and R. J. Composto, Macromolecules, 2014, 47, 875–887 CrossRef CAS.
- Y. Yan, S. C. Warren, P. Fuller and B. A. Grzybowski, Nat. Nanotechnol., 2016, 11, 603–608 CrossRef CAS PubMed.
- B. Rasin, H. Chao, G. Jiang, D. Wang, R. A. Riggleman and R. J. Composto, Soft Matter, 2016, 12, 2177–2185 RSC.
- M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2005, 95, 228301 CrossRef PubMed.
- M. Rechtsman, F. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 011406 CrossRef PubMed.
- É. Marcotte, F. H. Stillinger and S. Torquato, Soft Matter, 2011, 7, 2332 RSC.
- A. Jain, J. R. Errington and T. M. Truskett, Soft Matter, 2013, 9, 3866 RSC.
- M. Z. Miskin and H. M. Jaeger, Nat. Mater., 2013, 12, 326–331 CrossRef CAS PubMed.
- J. Qin, G. S. Khaira, Y. Su, G. P. Garner, M. Miskin, H. M. Jaeger and J. J. de Pablo, Soft Matter, 2013, 9, 11467 RSC.
- G. S. Khaira, J. Qin, G. P. Garner, S. Xiong, L. Wan, R. Ruiz, H. M. Jaeger, P. F. Nealey and J. J. de Pablo, ACS Macro Lett., 2014, 3, 747–752 CrossRef CAS.
- A. Jain, J. R. Errington and T. M. Truskett, Phys. Rev. X, 2014, 4, 031049 Search PubMed.
- W. D. Piñeros, M. Baldea and T. M. Truskett, J. Chem. Phys., 2016, 144, 084502 CrossRef PubMed.
- S. T. Milner, T. Witten and M. Cates, Macromolecules, 1988, 21, 2610–2619 CrossRef CAS.
- D. Meng, S. K. Kumar, J. M. D. Lane and G. S. Grest, Soft Matter, 2012, 8, 5002 RSC.
- D. M. Trombly and V. Ganesan, J. Chem. Phys., 2010, 133, 154904 CrossRef PubMed.
- P. Akcora, H. Liu, S. K. Kumar, J. Moll, Y. Li, B. C. Benicewicz, L. S. Schadler, D. Acehan, A. Z. Panagiotopoulos, V. Pryamitsyn, V. Ganesan, J. Ilavsky, P. Thiyagarajan, R. H. Colby and J. F. Douglas, Nat. Mater., 2009, 8, 354–359 CrossRef CAS PubMed.
- S. K. Kumar, N. Jouault, B. Benicewicz and T. Neely, Macromolecules, 2013, 46, 3199–3214 CrossRef CAS.
- A. Jayaraman and K. S. Schweizer, Macromolecules, 2008, 41, 9430–9438 CrossRef CAS.
- V. Pryamtisyn, V. Ganesan, A. Z. Panagiotopoulos, H. Liu and S. K. Kumar, J. Chem. Phys., 2009, 131, 221102 CrossRef PubMed.
- N. A. Mahynski and A. Z. Panagiotopoulos, J. Chem. Phys., 2015, 142, 074901 CrossRef PubMed.
- T. B. Martin, P. M. Dodd and A. Jayaraman, Phys. Rev. Lett., 2013, 110, 018301 CrossRef PubMed.
- T. B. Martin and A. Jayaraman, Macromolecules, 2013, 46, 9144–9150 CrossRef CAS.
- Y. Li, P. Tao, A. Viswanath, B. C. Benicewicz and L. S. Schadler, Langmuir, 2013, 29, 1211–1220 CrossRef CAS PubMed.
- Y. Li, T. M. Krentz, L. Wang, B. C. Benicewicz and L. S. Schadler, ACS Appl. Mater. Interfaces, 2014, 6, 6005–6021 CAS.
- J. Koski, H. Chao and R. A. Riggleman, Chem. Commun., 2015, 51, 5440–5443 RSC.
- B. Lin, T. B. Martin and A. Jayaraman, ACS Macro Lett., 2014, 3, 628–632 CrossRef CAS.
- N. Nair and A. Jayaraman, Macromolecules, 2010, 43, 8251–8263 CrossRef CAS.
- A. Walther, K. Matussek and A. H. E. Müller, ACS Nano, 2008, 2, 1167–1178 CrossRef CAS PubMed.
- C. E. Estridge and A. Jayaraman, ACS Macro Lett., 2015, 4, 155–159 CrossRef CAS.
- N. Hansen, S. D. Müller and P. Koumoutsakos, Evol. Comput., 2003, 11, 1–18 CrossRef PubMed.
- N. Hansen, in Towards a new evolutionary computation. Advances on estimation of distribution algorithms, ed. J. Lozano, P. Larranaga, I. Inza and E. Bengoetxea, Springer, 2006, pp. 75–102 Search PubMed.
- S. Sides, B. Kim, E. Kramer and G. Fredrickson, Phys. Rev. Lett., 2006, 96, 250601 CrossRef PubMed.
- J. Koski, H. Chao and R. A. Riggleman, J. Chem. Phys., 2013, 139, 244911 CrossRef PubMed.
- M. C. Villet and G. H. Fredrickson, J. Chem. Phys., 2014, 141, 224115 CrossRef PubMed.
- G. H. Fredrickson, The Equilibrium Theory of Inhomogeneous Polymers, Oxford University Press, New York, 2006 Search PubMed.
- D. A. McQuarrie, Statistical Mechanics, University Science Books, Sausalito, CA, 2000 Search PubMed.
- H. Chao, J. Koski and R. A. Riggleman, Soft Matter, 2017, 13, 239–249 RSC.

This journal is © The Royal Society of Chemistry 2018 |