Optimal control of electrostatic self-assembly of binary monolayers

A simple macroscopic model is used to determine an optimal annealing schedule for self-assembly of binary monolayers of spherical particles. The model assumes that a single rate-controlling mechanism is responsible for the formation of spatially ordered structures and that its rate follows an Arrhenius form. The optimal schedule is derived in an analytical form using classical optimization methods. Molecular dynamics simulations of the self-assembly demonstrate that the proposed schedule outperforms other schedules commonly used for simulated annealing.

In general, the term self-assembly refers to a non-intrusive process that transforms a system of components from a less to a more ordered state. There are numerous examples of selfassembly [1] including colloids [2], suspensions [3] and biological systems [4]. This work is motivated by self-assembly in model systems proposed by Grzybowski et al [5]. Those systems are horizontal monolayers of oppositely charged millimeter-size spherical particles placed inside a large square container. In such systems, the electrostatic interactions alone are not sufficient for arranging randomly placed particles into ordered lattices, as the particles are trapped in meta-stable states. To escape those traps, the container was subjected to vibration in the horizontal plane. With properly chosen vibration parameters it was possible on the one hand to escape meta-stable states and, on the other hand, preserve ordered lattice states. While Grzybowski et al were able to choose vibration parameters leading to self-assembly, their search strategy was heuristic and limited to very basic vibration histories. Our goal here is to develop an optimal strategy for guiding self-assembly in systems similar to those considered by Grzybowski et al. We note that the setup in [5] is not the only way to 2 obtain a square lattice in 2D: for example, it is possible to obtain a square lattice by tailoring the interaction potential in a single component system [6].
The systems considered by Grzybowski et al are very challenging for direct numerical simulations, primarily because they require modeling of friction and contact. Therefore we considered similar but simpler systems governed by basic rules of molecular dynamics (MD). We restricted ourselves to monolayers in which there were the same number of positively and negatively charged particles. Other than the charge sign, the particles were identical spheres of radius a and mass m. A perfectly ordered state for such monolayers is a square lattice. The electrostatic interactions of two particles with charges q i and q j was computed using Coulomb's law, where 0 is the vacuum permittivity and r is the distance between the particles. The container was treated as electrically neutral. Contact interactions were described by a repulsive potential where δ is the penetration distance between two overlapping surfaces. The penalty constant κ = 10 3 was chosen large enough to ensure small penetration distances. The particles agitated by the vibrating container were modeled as a canonical ensemble in which the temperature was maintained with a Nose-Hoover chain thermostat [7]. Throughout the paper, we measure energies in the units q 2 −1 0 a −1 and the temperature θ in the units q 2 −1 0 a −1 k −1 B , where k B is Boltzmann's constant. Time is measured in the units of q −1 √ 0 ma 3 , which is comparable to the period of vibration of two particles.
An optimization problem for a typical system can be stated as follows. Given an initial state of the system, determine a temperature history such that the system reaches a desired state as quickly as possible. Although this is a standard optimal control problem, in our case, we have to deal with a stochastic setting and a large problem size. In this regard, let us observe that we have already replaced a large number of the microscopic variables (velocities) with a single macroscopic variable (temperature). In what follows, we develop a simple macroscopic model that allows us to obtain a closed-form solution for the optimal temperature history. The proposed macroscopic model is based on an energetic description of the system. To this end, we consider the following energies: • U * is the potential energy of the instantaneous configuration.
• U is the potential energy of the instantaneous configuration subjected to a local energy minimization. Following [8], we refer to this configuration of mechanical equilibrium as the inherent structure. Figure 1 illustrates the difference between U * and U .
• V is the ensemble average of U . In the proposed model, we assume that U fully describes the system state and θ is an externally controlled input. At constant temperature, if the system is left to its own devices, it will eventually approach thermodynamic equilibrium, so that U → V . The simplest model describing this process is where τ is the relaxation time at a given temperature. In what follows, we extend the model one step further by assuming that it holds for non-isothermal processes. Under this assumption, we treat τ and V as functions of θ.
The function V (θ ) was obtained by fitting simulation results into the functional form in which the pre-exponential factor γ , the energy β, the melting temperature θ * and the ground state energy are treated as fitting parameters. We chose this functional form based on the equilibrium density of point defects in a perfect crystal at low temperature [9]. For each temperature, we computed U for multiple realizations, and then computed V as the average. The number of realizations was chosen such that at each temperature the standard error in V did not exceed 1.6 × 10 −3 (10% of the difference between V (θ max ) and ). Simulation results (discrete symbols) and the fit (continuous line) for γ = 0.007, β = 0.36, θ * = 0.018 and = −0.40 are shown in figure 2. These data show a transition between low-and high-temperature regimes near θ * . We are mostly interested in the low-temperature regime, because of its importance to self-assembly. In principle, obtaining the correct equilibrium value of V by averaging may not always be feasible due to large relaxation times at lower temperatures. To evaluate the feasibility of this approach in our computations we turn to a discussion of the relaxation time τ . The relaxation time τ (θ ) is chosen in the Arrhenius form where ν and α are dimensionless fitting parameters. These parameters were evaluated for a system of 400 particles using MD simulations at a fixed set of temperatures. In all cases, the initial configuration was obtained by equilibrating the system at temperature θ = 0.025. Then the system was allowed to relax to equilibrium at the desired temperature. The results of the simulations are shown in figure 3. The data were fitted in the range 0.6θ * θ θ * . For these data we obtain the Arrhenius fit with the energy barrier α = 0.043 and frequency ν = 0.012. Note that in this range of temperatures, dimensionless relaxation times estimated from (5) are in the thousands, which is well below the limitation imposed by the computational speed. Therefore it is feasible to compute V via ensemble averaging.
To rationalize the parameters of the Arrhenius fit in terms of a microscopic mechanism, we calculated a collection of energy barriers for diffusion on the surface and in the bulk. These calculations were done with the climbing-image nudged elastic band method to determine the minimum energy pathway and saddle point between minima [10]. Figure 4 shows the lowest energy mechanism for surface diffusion with a barrier of 0.05. In contrast, for bulk diffusion the barrier is in the range of 0.3-0.4. Thus we suggest that the likely rate-controlling mechanism is surface diffusion.  In summary, for n = 400 the model is governed by equations (3) and (4). The model parameters are Using standard methods of optimal control theory [11], one can formulate the following optimization problem. For the state equation the initial condition U (0) = U 0 and the termination condition Here U 0 and U f are prescribed values and t f is the final time to be determined as part of the optimal solution. The Hamiltonian corresponding to the stated optimal control problem is where λ is the co-state variable. It allows us to express the optimality condition for θ as which yields the optimal relationship between U and θ : With this relationship, equation (6) becomes an ordinary differential equation for θ . After obvious manipulations we obtain the optimal schedule with the initial temperature Let us note that logarithmic time-dependency, similar to that in equation (11), commonly arises in optimal scheduling [12]. The final time is expressed in the form with the final temperature It is straightforward to establish that the optimal schedule can be also used to achieve the minimum U f for a fixed t f . This view of optimality is adopted in simulations described next. We compared schedule (11) with four types of schedules commonly used in simulated annealing, which are listed in table 1. Schedule (i) is equivalent to vibrating the container at a fixed intensity and is directly related to the strategies implemented in the experiments of Grzybowski et al. Schedule (ii) results in linear temperature history. Schedule (iii) lowers the energy at a constant rate [13] and (iv) is a constant thermodynamic speed schedule [14]. Schedules (iii) and (iv) require knowledge of V (θ) and, in addition, in the case of (iv), τ (θ). To solve for the schedules, we used equations (4) and (5) for V (θ) and τ (θ), respectively.  Figure 5. Comparison of the performance of schedules (i)-(v). The dashed and the black solid curves represent, respectively, the energy U (θ(t)) (from equations (10) and (11)) and the equilibrium energy V (θ(t)) (from equations (4) and (11)) for schedule (v).
All simulations involved self-assembly of 400 particles. The performance of each schedule was evaluated using the parameter: withṼ = V |θ =θ = −0.38. Initial positions were obtained by equilibrating the system atθ = 0.5 θ * . The termination time t f = 10 000 was chosen such that it was sufficient for the average value of η to drop below 0.1 when using schedule (v). To choose the constant operating temperature for schedule (i), we conducted simulations at three different temperatures θ = {θ 0 , 1 2 (θ 0 + θ f ), θ f }. Here θ 0 and θ f are, respectively, initial and final temperature of schedule (v) for the chosen Figure 6. Evolution of the microstructure for a system of 400 particles. simulation time. We found θ = θ 0 = 0.94θ * to produce the best results based on the optimality condition (7). The initial temperature for the constant rate schedules (ii)-(iv) was θ 0 = θ * . The final temperature θ f = 0.5θ * was determined from equation (5) for τ = t f . This represents a temperature with the characteristic relaxation time comparable to the total simulation time t f . The proportionality constant ξ in (iii) and the thermodynamic speed v s in (iv) were treated as parameters whose values were determined by our choice of θ 0 , θ f and t f .
Simulations results are plotted as log η versus t in figure 5. The results lead to several conclusions. Firstly, schedule (v) outperforms its competitors up to t ≈ 6000, where a crossover with schedules (i) and (iv) occur. Secondly, the macroscopic prediction (10) agrees well with the microscopic simulations only up to t = 2500. Thirdly, self-assembly occurs far from equilibrium, independent of the schedule used. This is evident by comparing V [θ(t)] to all other histories.
To better understand the discrepancy between the macroscopic and microscopic model predictions for t > 2500 we show the evolving microstructure in figure 6. Up to t = 2500, microstructural changes are primarily associated with the elimination of internal defects by particle rearrangements along internal boundaries. By t = 2500 (third insert), most of the internal defects have been eliminated and the particles have been assembled into one ordered cluster (or a few ordered clusters, in some cases). The subsequent transformation of this structure to that of a perfectly ordered one likely involves a different relaxation mechanism. The existence of (at least) two distinct relaxation mechanisms has been previously noted in supercooled liquids and glasses as well as ionic crystals [15]. This helps explain why the constant thermodynamic speed schedule (iv) is able to overcome schedule (v) in figure 5, if it is able to better capture slower relaxation mechanisms at longer times.
It has been previously demonstrated that constant thermodynamic speed schedules are optimal with respect to thermodynamic efficiency, which, in general, is not equivalent to maximum power production [16,17]. This is confirmed by the fact that schedule (v) performs better than (iv) up to t = 6000. This allows us to conjecture that the true optimal schedule can be derived in the framework of optimal control theory but must be adaptive with respect to the change of the rate-controlling mechanism. Another interesting observation from figure 6 is that for the case of macroscopically small systems such as those considered in this paper, capturing one of the relaxation mechanisms is almost sufficient to drive the system into the perfectly ordered state.
In summary, we proposed a simple macroscopic model that provides a useful description of self-assembly in systems considered by Grzybowski et al [5]. Based on the proposed model, an optimal annealing schedule is derived in an analytical form using classical optimization methods. The model predictions agree well with results of microscopic simulations under conditions when the model captures correctly the rate-controlling mechanism. Furthermore, in this regime, the optimal schedule outperforms commonly used simulated annealing schedules.