\documentclass[twocolumn,floatfix,showpacs,preprintnumbers,amsmath,amssymb,prb,superscriptaddress]{revtex4}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math

\begin{document}

\title{Self-assembly of monodisperse clusters: Dependence on target geometry}

\author{Alex W.~Wilber}
\affiliation{Physical and Theoretical Chemistry Laboratory, University of
  Oxford, South Parks Road, Oxford OX1 3QZ, United Kingdom}
\author{Jonathan P.~K.~Doye}
\thanks{Authors for correspondence}
\affiliation{Physical and Theoretical Chemistry Laboratory, University of
  Oxford, South Parks Road, Oxford OX1 3QZ, United Kingdom}
\author{Ard A.~Louis$^*$}
\affiliation{Theoretical Physics, University of Oxford, Keble Road, Oxford OX1 3NP, United Kingdom}

\date{\today}

\begin{abstract}
We systematically make interesting shapes.
\end{abstract}

\pacs{Don't know}

\maketitle

\section{Introduction}

%Background
Self-assembly is very important in biology, with viruses as a classic example, but it is not very well understood.
There have been a number of theoretical studies conducted (Zlotnick etc) which have shed some light on the process,
but frequently they have had to make certain assumptions about the mechanisms of assembly.
Other simulation studies (Rapaport, Chandler, Van Workum, Brooks) also contribute some understanding.

In previous work (last paper) we considered a simple model system of spheres with patchy attractions assembling
to form icosahedra. *Advantages of this approach*
The success with which assembly occurred was examined as a function of temperature and the width
of the patches. At least two
mechanisms of assembly were found to be at work. The first was a simple nucleation mechanism in which icosahedra
were built up from monomers and smaller clusters. This mechanism was found to be most effective close to the maximum
temperature at which icosahedra were stable, since at these temperatures kinetic traps are not so important. A second
mechanism involved the formation of larger clusters which then rearranged to produce icosahedra. This mechanism was
more prevalent at higher patch widths and lower temperatures, where large clusters were more easily formed. At the
optimal conditions for assembly both mechanisms made a contribution.

- Did we find anything else out worth talking about?

\subsection{Target structures (probably remove this header)}

We now seek to extend this study from icosahedra to the other platonic solids (tetrahedra, octahedra, cubes and dodecahedra),
in order to investigate how general the findings are, and in what ways the shape of the target structure affects the ease
and mechanisms of assembly. The platonic solids represent relatively simple targets for assembly because of their high symmetry
- all particles and patches are equivalent. However, they are DIVERSE.

<PICTURES OF TARGET STRUCTURES>
We chose to study the self-assembly of the platonic solids, since they represent interesting but attainable targets for assembly.
Fig.\ \ref{fig:TargetStructures} shows the particle designs and target structures for each. The particles are designed such that
patches are pointing directly at neighbouring particles if assembled into the target structure, so that the target structure
represents an energy minimum.

\section{Methods}

\subsection{Model}

We make use of the same patchy sphere model as in previous work \cite{me}. The interaction potential consists of an isotropic
repulsion, plus a patchy attraction:

\begin{equation}
\label{eq:potential}
V_{ij}({\mathbf r_{ij}},{\mathbf \Omega_i},{\mathbf \Omega_j})=\left\{
    \begin{array}{ll}
       V_{\rm LJ}(r_{ij}) & r<\sigma_{\rm LJ} \\
       V_{\rm LJ}(r_{ij})
       V_{\rm ang}({\mathbf {\hat r}_{ij}},{\mathbf \Omega_i},{\mathbf \Omega_j})
                       & r\ge \sigma_{\rm LJ}, \end{array} \right.
\end{equation}
where $r_{ij}$ is the interparticle vector and ${\mathbf \Omega_i}$ is the orientation of particle $i$.
$V_{\rm LJ}$ is the Lennard-Jones potential,
\begin{equation}\label{eqn:LJ} V_{\rm LJ}(r) = 4\epsilon\left[ \left( \frac{\sigma_{LJ}}{r}
    \right)^{12} - \left( \frac{\sigma_{LJ}}{r} \right)^{6} \right], \end{equation}
and the attraction is modulated by an orientationally dependent term,
\begin{eqnarray}
V_{\rm ang}({\mathbf {\hat r}_{ij}},{\mathbf \Omega_i},{\mathbf \Omega_j})&=&
G_{ij}({\mathbf {\hat r}_{ij}},{\mathbf \Omega_i})
G_{ji}({\mathbf {\hat r}_{ji}},{\mathbf \Omega_j})\\
\label{eqn:AngMod}
G_{ij}({\mathbf {\hat r}_{ij}},{\mathbf \Omega_i})&=&
\exp\left(-{\theta_{k_{\rm min}ij}^2\over 2\sigma^2}\right).
\end{eqnarray}
where $\sigma$
is the standard deviation of the Gaussian,
$\theta_{kij}$ is the angle between patch vector $k$ on particle $i$
and the interparticle vector $\mathbf r_{ij}$,
and $k_{\rm min}$ is the patch that minimizes the magnitude of this angle.
Hence, only the patches on each particle that are closest to the interparticle
axis interact with each other.

When using this potential in our Monte Carlo
simulations, we truncate and shift the potential using a cutoff distance of
3\,$\sigma_{LJ}$ in order to improve the computational efficiency. To maintain the
isotropic nature of the repulsion, we also adjust the distance at which we start
to include the angular modulation so that it corresponds to where the cut-and-shifted
Lennard-Jones potential passes through zero.

\subsection{Monte Carlo simulation}
\label{sec:MonteCarlo}

The bulk of our simulations consist of Metropolis Monte Carlo (MC) simulations in the canonical ensemble, using periodic boundary conditions.
The allowed move types are small single-particle translations and rotations. The translational moves are randomly chosen from a small cube centred
on the selected particle. Rotational moves make use of a quaternion description of the particle's orientation; the proposed quaternion
is given by the renormalised sum of the current quaternion and a smaller, randomly generated 4-vector.

This method does not directly simulate the physical trajectories of the particles, and is more often used to obtain thermodynamic averages of the
properties of a system. However, the restriction to small moves results in the system approximating diffusive dynamics. This method has the advantage
of being very efficient, since derivatives of the potential (which are expensive to calculate) are not required.

As the above Monte Carlo simulations do not generally reach equilibrium, we make use of umbrella sampling to obtain equilibrium properties. Under
this scheme instead of choosing the acceptance probabilities using a Boltzmann distribution, the modified distribution 
$\exp\left[-\beta\left(\mathcal{V}+w(Q)\right)\right]$ is used, where $Q$ is an order parameter
and $w(Q)$ is a weighting function. The weighting function is, in general,  designed to artificially reduce the free energy barriers between different
(meta)stable states of a system. Canonical averages are simple to obtain from a simulation of such a non-Boltzmann (nB) ensemble using the expression
\begin{equation} \left<B\right>_{NVT} = \left<B \exp \left[ \beta w(Q) \right] \right>_{\mathrm{nB}}, \end{equation}
where $B$ is some generic property of the system.

In our case, we wish to locate the temperature at which the system becomes most stable as a collection of completed target clusters, 
and so we need to overcome the free energy barrier for the formation of this structure. As the
completed clusters only very weakly interact with each other, it is sufficient, to a first approximation,
just to consider a system consisting of the required number of particles to form one complete cluster.
For an order parameter we simply use the number of particles in the largest cluster present. We wish to choose
the weight function so that the system spends an approximately equal amount of time at each value of $Q$. To
achieve this we use an iterative scheme in which $w(Q)$ is updated at the end of each of a series of simulations.

\subsection{Umbrella sampling}

\subsection{Parallel tempering?}

\subsection{Cluster growing simulation}
\label{sec:ClusterGrowing}

Large kinetic and thermodynamic aggregates play an important role in our self assembly simulations, but
it can be difficult to obtain detailed information about them. We therefore set about designing a simple
simulation method to examine them more closely. In particular we were interested in whether aggregates tend to
have unbonded patches on their surfaces and in what number, and whether particles in larger aggregates were able
to arrange themselves to form effective bonds through most of their patches.

A simulation starts with a single particle with the desired set of patches. A second particle is then added, such
that patches on the first and second particles line up to form a strong bond. The system then undergoes a short
MC simulation of $1000$ MC cycles at low temperature (we used $T = 0.03$) to allow the system to find nearby low-energy 
states - although with only two
particles this will have little effect. Following the MC simulation each patch in the system is considered in turn,
to see whether another particle could be added at that site or whether the locality of the patch is obstructed by
nearby particles, or in fact occupied by a bonding particle. One of the unobstructed patches is chosen at random
and a third particle is placed at that site, followed again by a low temperature MC simulation and enumeration of
the free patches, and so on. The number of free patches and the overall energy of the system is recorded as a function
of the number of particles in the aggregate.

\section{Results}

\subsection{Overview}

%In order to investigate the self-assembly behaviour of our system we performed arrays of Monte Carlo (MC) simulations
%in the canonical ensemble. <POSSIBLY INSERT A TYPICAL TRAJECTORY OR TWO HERE, ALONG WITH AN EXPLANATION OF HOW ASSEMBLY
%PROCEEDS>.

%Because of the relative simplicity of our model the simulations are inexpensive, and so we are able to map out the
%behaviour of the system in detail as a function of temperature and the width of the patches. Figs.\ \ref{fig:MainPlots}
%and \ref{fig:SizePlots} summarize the results of a large number of simulations at different combinations of $\sigma$
%and $T$, for each of the different target structures. Fig.\ \ref{fig:MainPlots} shows the percentage yield of the target
%structure for each simulation, while Fig.\ \ref{fig:SizePlots} shows the average size of the clusters present at the end
%of each simulation. Note that there is no plot included for dodecahedra in Fig.\ \ref{fig:MainPlots}, since no correctly
%formed dodecahedra were observed in any of the simulations.

%It can be seen from Fig.\ \ref{fig:MainPlots} that for each of the target structures there is a temperature window within
%which successful assembly takes place. At high temperatures the target structures cease to be stable, and the particles
%primarily remain as monomers with some transient small clusters present. At low temperatures, by contrast, interactions between
%particles become very hard to break and have long lifetimes. During assembly it is very likely that erroneous bonds will
%be formed, which are not compatible with the target structure. For successful assembly to occur it is therefore important
%that the temperature is high enough for such erroneous bonds to be broken, and hence at low temperatures the system will suffer from
%kinetic trapping and will form flawed structures with long lifetimes.

%Fig.\ \ref{fig:MainPlots} also shows that there is a window of patch widths within which successful assembly occurs. Low patch
%widths inhibit assembly through both kinetic and thermodynamic effects. From the kinetic standpoint, narrow patches make it difficult
%for particles to find each other, and bonds are only rarely formed. The restricted vibrations of the resultant bonded structures also
%lead to a greater entropy penalty on formation, so that bonds formed from narrow patches are thermodynamically less stable. High patch
%widths also inhibit successful assembly, since the target structure eventually ceases to be the most stable configuration and more compact,
%liquid-like structures dominate. Fig.\ \ref{fig:gmin} shows the minimum energies of clusters of icosahedron-forming particles as a function
%of the cluster size, and for different patch widths $\sigma$. The zero axis represents the energy per particle of correctly formed icosahedra.
%It can be seen that for narrow patch widths, the optimal structures are simply multiple correctly formed icosahedra, but beyond a critical
%patch width of $\sigma=0.6$ other structures have a lower energy for $N_c>12$ [have not explained $N_c$, or indeed any of this well].

%<gmin plot from Jon>

%Conclusion: You can only form correct clusters for small values of $\sigma$. For larger
%values liquid-like droplets become stable and cluster formation is no longer favourable.

\subsection{Overview}

\begin{figure}
\includegraphics[width=8.2cm]{figures/mainPlots.eps}
\caption{\label{fig:MainPlots}(Colour Online)
  Main heat maps.}
\end{figure}

\begin{figure}
\includegraphics[width=8.2cm]{figures/sizePlots.eps}
\caption{\label{fig:SizePlots}(Colour Online)
  Average size heat maps.}
\end{figure}

Because of the relative simplicity of our model the simulations are inexpensive, so we are able to map out the
behaviour of the system by performing a large number of simulations. Heat maps showing the yields of target clusters for the different target structures
as a function of temperature $T$ and patch width $\sigma$ are shown in Fig.\ \ref{fig:MainPlots}, while the average sizes
of the clusters formed in each simulation are shown in Fig.\ \ref{fig:SizePlots}. The same general
trends are observed throughout. In each case there is a region of relatively successful assembly, bordered on three sides by regions in which
distinct competing structures are formed instead. At high temperatures the target clusters cease to be stable, and the system consists of an equilibrium
mixture of monomers and small clusters. At low temperatures, by contrast, the system tends to become kinetically trapped. Large aggregates are formed
and are unable to rearrange to form the lower energy clustered state. Finally, at high patch widths liquid-like aggregates become the stable state.
Because of the wide patches the particles are able to form bonds even if they are not correctly oriented, while the relative disorder of the liquid
droplets results in a favourable entropy compared to the clusters. In the following sections we consider the transitions to each of these competing
states in turn, and examine how they differ for the various target structures.

\subsection{Transition to a monomer gas}

The temperature above which target clusters are unstable is labelled $T_{clust}$, and is superimposed on the plots in Fig.\ \ref{fig:MainPlots}.
These values were obtained by performing umbrella sampling simulations to produce equilibrium heat capacity plots, with the plotted data points
in Fig.\ \ref{fig:MainPlots} corresponding to the peaks in the heat capacities. Examples of heat capacity plots, given for $\sigma = 0.4$, are shown
in Fig.\ \ref{fig:HeatCapacities}.

It is intuitively clear that $T_{\mathrm{clust}}$ will increase with the number of bonds per particle (since the clusters will have greater
stability), and decrease with cluster size (since larger clusters pay a greater price in entropy). These factors largely cancel out when comparing
tetrahedra, octahedra and icosahedra, resulting in similar values of $T_{\mathrm{clust}}$ for each at a given value of $\sigma$. However, the values
of $T_{clust}$ for cubes and dodecahedra are markedly lower, since they only have three bonds per particle but are much larger clusters than tetrahedra.

[Why the graphs are wider for smaller clusters: finite size effect]

[Back-of-the-envelope calculation?]

\subsection{Kinetic trapping transition}

The kinetic trapping transition represents an approximate temperature $T_{agg}$ for a given value of $\sigma$ below which
the system is largely unable to escape incorrect configurations, and does not produce a significant yield of the
target clusters over the timescale of the simulation. Below the transition the system tends to produce large
glassy networks of particles, often spanning the simulation box and involving almost all of the particles in the
simulation. The character of the networks varies with the patch width $\sigma$. At low $\sigma$ the ``strings''
of the network tend to be only a single particle across, but are highly interconnected. At higher $\sigma$ the
particles are able to rearrange much more and form thicker strings reminiscent of worm-like micelles \ref{}.
[CHECK THIS IS STILL TRUE FOR NON-ICOSAHEDRAL SHAPES] Eventually at very high $\sigma$ there is a smooth transition 
to the more globular aggregates characteristic of the liquid phase.

The value of the kinetic transition temperature $T_{agg}$ arises from two other characteristic temperatures,
which we shall denote $T_{agg}^{meta}$ and $T^{agg}_{trap}$. $T_{agg}^{meta}$ is the temperature below which aggregates become metastable. Above this
temperature the formation of aggregates cannot take place, since the monomer gas is more stable. Note that $T_{agg}^{meta}$
must always lie below $T_c$ for values of $\sigma$ at which the target clusters are the stable structure, giving at least
a narrow temperature range in which the target clusters are the only aggregates which can form. In general $T_{agg}^{meta}$
will increase with $\sigma$, since aggregates are stabilised both energetically and entropically by higher patch widths.

The second relevant temperature, $T^{agg}_{trap}$, represents the approximate temperature below which the ``incorrect''
bonds in aggregates can be considered to have very long lifetimes, and the system becomes trapped. In general $T^{agg}_{trap}$
decreases with $\sigma$, since wider patches reduce the energy barriers for rearrangements.

The system forms large kinetic aggregates only when the temperature is less than both $T_{agg}^{meta}$ and $T^{agg}_{trap}$, as shown
schematically in Fig.\ \ref{fig:TAggSketch}. This explains the shape of the kinetic aggregate region seen in Figs.\ \ref{fig:MainPlots}
and \ref{fig:SizePlots}.

-- Pictures of aggregates
-- Structures, dangling bonds and energies. Geometric effects on difficulty in rearranging to correct structure.
-- More structural things: Formation of chains, chains of the interlocking type, interpenetration etc
-- Effect of splay angle
-- Ring counting

It is noticeable from Fig.\ \ref{fig:SizePlots} that the sizes of the kinetic aggregates formed by
particles designed for different target structures vary wildly. In order to obtain a better understanding of this
phenomenon 100 cluster growth simulations, as described in Section \ref{sec:ClusterGrowing}, were performed for each
of the five particle designs. The average number of unblocked patches and the average energy per pair of patches are shown 
in Fig.\ \ref{fig:ClusterGrowing}, along with the population of clusters which have succeeded in growing to a given size without
running out of free patches. Note that when a cluster reached zero free patches and could be grown no further it
contributed a zero to the average number of free patches for all higher particle numbers, but was not included in the
average energy. [MAYBE INCLUDE THIS IN THE CAPTION?]

THIS IS NOW WRONG ---------
Fig. \ref{fig:ClusterGrowing} shows that clusters of tetrahedron- and octahedron-forming particles tend to become `closed',
with no remaining dangling patches, before they grow very large. This explains the fact that no large aggregates are observed
in the low-$T$ low-$\sigma$ region of the corresponding heat maps in Fig.\ \ref{fig:MainSizes}. For cubes, icosahedra and
dodecahedra however, the number of free patches continues growing with aggregate size, and very few aggregates become closed,
hence the very large aggregates seen in Fig.\ \ref{fig:MainSizes}.
---------------

AGAIN NEEDS FIXING ----------------
The reasons for these differences are not entirely clear. One important factor is what we term the 'splay angle' of the patches -
that is, the angle at which each of the patches lies from the symmetry axis of the particle. The splay angles for the different
particles are shown in Table\ \ref{tab:SplayAngles}. Some correlation is apparent, since the two designs which tend to form
closed clusters (tetrahedron- and octahedron-forming particles) also have the lowest splay angles. One can imagine that a small splay
angle tends to create greater curvature in the surface of an aggregate and encourages closure. However, this is clearly not the whole story,
since octahedron-forming particles have the greatest tendency to form closed particles but do not have the smallest splay angle. Presumably
the specific behaviour of each design arises from the detail of the geometry of the particles and how they can pack while maximising bonding energies.
--------------------

More on this:

- Tetrahedra and octahedra form chains   ----------- CHECK!!

- Tetrahedra appear to be more likely to have forks in the chains. Large aggregates tend to have many forks.

- Cubes form a more compact structure, though with some linearity

- Icosahedra form very compact structures

- Dodecahedra form interlinking rings. Wide angles mean the structure is less dense, but it is compact in the sense of not consisting of chains or sheets.

\begin{figure}
\begin{tabular}{|r|r|}
\hline
Target structure&Splay angle (degrees)\\
\hline
Tetrahedron&35.3\\
Octahedron&45.1\\
Cube&54.7\\
Icosahedron&58.3\\
Dodecahedron&69.1\\
\hline
\end{tabular}
\caption{\label{tab:SplayAngles}Splay angles for particles designed to form the various target structures.}
\end{figure}

PICTURES OF KINETIC AGGREGATES

\subsection{Cluster-liquid transition}

[put gmin plot in here, if it wasn't in earlier]

At high values of $\sigma$ the stable state of the system is no longer a gas of clusters with the target structure, but rather
one or more liquid-like aggregates. High $\sigma$ allows strong interactions even if patches are poorly aligned, and particles may form
strong attractions with more neighbours than they have patches. Liquid droplets have higher entropy than ordered clusters, so the value of
$\sigma$ above which the liquids are stable decreases with increasing temperature.

The structure of the aggregates formed at high $\sigma$ is very different to that of the kinetic aggregates formed at lower $\sigma$ and
low temperature. The kinetic aggregates tend to form extended networks, and the total bonding energy of the system is relatively small. The
thermodynamic liquid-like aggregates, by contrast, tend to be more compact, with a total bonding energy competing with or even exceeding
that found in the target clusters (i.e.\ $\epsilon/2$ per bond). However, although more compact, the liquid-like aggregates are often very
far from spherical. If all of the patches on the particles are pointing inwards the aggregate has a very small surface tension, and therefore
has little drive to adopt a spherical configuration.

- I get the impression that the cube liquid is relatively stable in comparison to the clusters
  since it appears to infringe on the cluster region at quite low values of $\sigma$.

- Perhaps generally cluster-liquid transition moves to left as a function of size (seems reasonable on entropy grounds), with exception of
icosahedra, which only bond really poorly in liquid state. [do they?]

-- Importance of splay angle
-- Short range similarity to cluster structure (ring counts, adf, rdf?)

PICTURES OF THERMODYNAMIC AGGREGATES [should be combined with pictures of kinetic aggregates?]

\subsection{On formation by chance}

- Tetrahedra can form even at very low patch widths and temperatures, just by chance.
  Similarly at high patch widths and temperatures. Basically anywhere except where the stable
  state is monomers or liquids. (what is going on in the top right?)  

- From the cluster growing experiments it appears that a little over 10\% of stepwise-growing
  clusters will form tetrahedra.

\subsection{Energy landscapes}

\subsection{Summary of consequences for different shapes}

- Tetrahedra form well - good Tclust, low tendency to form large aggregates, dynamically simple & quick to assemble, may even assemble
by chance under otherwise unfavourable conditions.

- Octahedra are similar only somewhat less so

- Cubes are badly squeezed and are sad. They have a great tendency to form large aggregates which will tend to mop up free particles and
inhibit assembly. The liquid is stable and competes strongly with the clusters. Tclust is low because of the low number of bonds and reasonable
size. Perhaps squares are bad because of pentagons?

- Icosahedra recover considerably because of their large number of bonds (high Tclust) and poor energy in the liquid state (transition to liquid moved
to the right). Kinetic aggregates do cause thorough blocking in the region in which they exist, though.

- Dodecahedra are a disaster since (a) Very low Tclust, (b) large kinetic aggregates, and (c) presumably low-sigma cluster-liquid transition. I expect Tclust is the main issue though. While there is presumably some temperature window in which clusters would form, it may be very narrow, and formation may be very slow.

\subsection{Mechanisms of assembly/Dynamics}

-- Yield vs time (with thermo lines) - why shape differences?
-- Ard plots
-- Does budding happen for other shapes?
-- Geometric effects on dynamics? (``are dodecahedra hard to form because they're clumsy?'')

WANT TO INCLUDE PICTURE OF BUDDING OFF HERE

\begin{figure}
\includegraphics[width=8.2cm]{figures/tetraYields.eps}
\includegraphics[width=8.2cm]{figures/octaYields.eps}
\includegraphics[width=8.2cm]{figures/cubeYields.eps}
\includegraphics[width=8.2cm]{figures/icosaYields.eps}
\caption{\label{fig:YieldPlots}(Colour Online)
  Yields of tetrahedra, octahedra, cubes and icosahedra.}
\end{figure}

\begin{figure}
\includegraphics[width=8.2cm]{figures/tetraSizeWTime.eps}
\includegraphics[width=8.2cm]{figures/octaSizeWTime.eps}
%\includegraphics[width=8.2cm]{figures/cubeSizeWTime.eps}
\caption{\label{fig:SizeWTimePlots}(Colour Online)
  Average cluster size as a function of time for tetrahedra, octahedra and cubes.}
\end{figure}

\begin{figure}
\includegraphics[width=7.0cm]{figures/tetraArd.eps}
\includegraphics[width=7.0cm]{figures/octaArd.eps}
\includegraphics[width=7.0cm]{figures/cubeArd.eps}
\caption{\label{fig:ArdPlots}(Colour Online)
  Ard plots for tetrahedra, octahedra and cubes.}
\end{figure}

\subsection{Noise}

ARE WE GOING TO INCLUDE NOISE?

\section{Conclusions}

Overall, I think I like shapes even more than when I started.

\begin{acknowledgments}
The authors are grateful for financial support from the EPSRC and the Royal Society.

Are we still grateful to the Royal Society?
\end{acknowledgments}

\bibliography{references}

\end{document}
