\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}

Basic question: What does one need to get particles to self-assemble into a finite-size, monodisperse, ordered object? e.g. virus capsid, but
beyond to think generally

Review what's been done to begin to answer this question.

Specific Qs to this paper:

(1) Beyond ref[previous paper] to look at dependence of difficulty on target geometry [nb restrict targets to platonic solids i.e. each particle
in identical environment and one type of patch].

(2) Interactions directional, but non-torsional, i.e. although the model is general it is more appropriate to patch colloids, not proteins

Review where people are at with making patchy colloids.

Here, we focus on platonic solids.

\section{Methods}

\subsection{Model}

We make use of a minimal model, designed to contain only the essential features required for targeted self-assembly, while allowing for efficient
simulation. The model consists of spherical particles patterned with attractive patches. They are described by a modified Lennard-Jones potential,
in which the repulsive part of the potential is isotropic as usual, but the attractive part is anisotropic and depends on the alignment of patches
on interacting particles. Specifically, the potential is described by
\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 $V_{\rm LJ}$ is the Lennard-Jones potential, given by
\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}
$V_{\rm ang}$ is an angular modulation factor, which depends on the orientations of the patches on the two interacting particles, as well as
the direction of the vector joining them. Specifically,
\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_{\rm ang}$ gives the width 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.

A particular particle design is specified simply by a set of unit vectors describing the positions of the attractive patches.
For each of our target structures there is an associated particle design, in which the patches are placed such that they
point directly at the neighbouring particles in the target structure.

\begin{figure}
\includegraphics[width=7.4cm]{figures/shapePics.eps}
\caption{\label{fig:MainHeatPlots}(Colour Online)}
\end{figure}

\subsection{Simulation methods}

\subsubsection{Use of Metropolis Monte Carlo for studying dynamics}

Emphasise reproduces Brownian motion. Get Refs from Jon if required.

Address possible flaws with lack of cluster moves. Not important.

\subsubsection{Umbrella sampling}

One/two dimensional. What order parameters used, which most successful. How calculated if necessary.

\subsubsection{Hamiltonian exchange}

Dodecahedron was pathologically hard, so needed extra method.

\section{Results}

\subsection{Overview}

Because of the relative simplicity of our model, we are able to perform large numbers of simulations, and
hence map out the assembly behaviour over a range of parameter space. Fig.\ \ref{fig:MainHeatPlots} shows the
yields of each of the platonic solids in simulations over a range of values of the temperature $T$ and the patch width $\sigma$.

%%
Note that throughout, we refer to the plots as ``the tetrahedral plot'' and so forth.

\begin{figure}
\includegraphics[width=7.4cm]{figures/noTorYieldPlots.eps}
\caption{\label{fig:MainHeatPlots}(Colour Online) The numbers of
(a) tetrahedra, (b) octahedra, (c) cubes, and (d) icosahedra formed after $80\,000$ MC cycles
as a function of the patch width $\sigma$ (measured
in radians) and the temperature for systems of 1200 particles at a number density of
$0.15\,\sigma_{\rm LJ}^{-3}$. No plot is included for dodecahedra, because no dodecahedra
were succesfully assembled under any conditions. The white lines show the
thermodynamic transition temperature $T_c$ for the transition from
a gas of clusters to a gas of monomers, calculated for a one-cluster system (see Section \ref{sec:Thermo}).}
\end{figure}

Significant differences between the plots are apparent. Tetrahedra,
octahedra and icosahedra assemble readily over wide ranges of parameter space. By contrast, cubes are more difficult to assemble and are only
formed in a smaller region of parameter space, while dodecahedra never form at all. We shall return to the reasons for these differences later.
However, the same physical principles apply in each case, and as a result the
general shapes of the plots are similar.

%Both similarities and differences between the plots are apparent. We will first explain the underlying reasons for the general
%shapes of the plots, which also explain their similarities. We will then proceed to discuss the differences between the results
%for the different target structures.

The region of successful assembly is determined by a number of constraints, both thermodynamic and kinetic. To begin with the first
of the thermodynamic effects, we see that at high
temperature $T$ the stable state becomes a gas of monomers and small clusters, and target clusters are no longer formed. We
designate the approximate temperature of this transition $T_c$, and $T_c$ lines are marked on Fig.\ \ref{fig:MainHeatPlots}
(the precise definition we have used to obtain $T_c$ is discussed in Section\ \ref{thermo}). A second thermodynamic constraint
arises at high patch width $\sigma$. At high $\sigma$ the patches are so wide that the target clusters cease to be the stable
state of the system. Each patch becomes capable of interacting with more than one neighbour, and so the system can both lower its
energy and raise its entropy by forming large, unstructured liquid-like droplets.

Both of these effects can be seen in Fig.\ \ref{fig:SizeHeatPlots}, which shows the average cluster size for each target structure,
again as a function of $T$ and $\sigma$. At high $T$ the average cluster size approaches unity, signifying a monomer gas, while
at high $\sigma$ the existence of a liquid phase results in clusters containing essentially all the particles in the simulation.

\begin{figure}
\includegraphics[width=7.4cm]{figures/noTorSizePlots.eps}
\caption{\label{fig:SizeHeatPlots}(Colour Online) The mean cluster size (averaged over particles) of systems
of particles designed to form (a) tetrahedra, (b) octahedra, (c) cubes, (d) icosahedra, and (e) dodecahedra, for the same
simulations as Fig.\ \ref{fig:MainHeatPlots}. Again, the white lines show the
thermodynamic transition temperature $T_c$ for the transition from
a gas of clusters to a gas of monomers, calculated for a one-cluster system (see Section \ref{sec:Thermo}).}
\end{figure}

Kinetic constraints become important at low values of $T$ and $\sigma$. At low $T$ the system becomes unable to escape from incorrect
configurations. Since it very likely that some incorrect bonds will be formed this results in very low yields, and the system instead
forms glassy kinetic aggregates, often forming a diffuse network spanning the simulation box. At low values of $\sigma$ the patches
are very narrow. Particles will only rarely be sufficiently well aligned to feel attractive interactions from their neighbours, and
growth is suppressed.

Sandwiched in the middle of these regions of inhibited assembly is a region of parameter space in which the target is both thermodynamically
stable and kinetically accessible; in this region good yields of the target structure are obtained.

These effects are summarised in Fig.\ \ref{fig:SchematicPhaseDiagram}, which shows schematically the different regions of parameter space and the
lines which separate them.

\begin{figure}
\includegraphics[width=8.4cm]{figures/SchematicPhases.eps}
\caption{\label{fig:SchematicPhases}A schematic diagram illustrating the regions of parameter space in which different
behaviours are observed, and the interfaces between them. $T_c$ is the approximate midpoint of the transition between the cluster
gas and the monomer gas (see Section\ \ref{thermo}). $T_{lv}$ and $T_{cl}$ are the liquid-vapour phase transition lines for the
monomer and cluster gases respectively. $T^{agg}_g$ is the approximate temperature below which aggregates become unable to rearrange
to form target clusters on the timescale of the simulations, and instead become trapped in a glass-like state. $\sigma_k$ is the approximate
value of $\sigma$ below which the patches are so narrow that little assembly or aggregation is able to take place on the timescale of the
simulations.}
\end{figure}

We will now examine some of the above effects in more detail.

\subsection{The cluster-monomer gas transition}
\label{sec:TcTransition}

The plots in Fig.\ \ref{fig:MainHeatPlots} show that $T_c$, the temperature above which the stable state is dominated by monomers rather than target
clusters, increases with patch width $\sigma$ regardless of the target structure. We can explain this by considering the free energies of the monomer
and cluster gases (we consider a state dominated by target clusters to be a gas, because the completed clusters have all patches pointing inwards and
feel only very weak attractive forces). At $T_c$, $A_m$ and $A_c$, the free energies of the monomer and cluster gases respectively, are approximately
equal. $A_m$ is largely independent of $\sigma$, as is the configurational energy per cluster in the cluster gas, $V_c \approx \frac{m \epsilon}{2}$,
where $m$ is the number of patches per particle. However, the entropy $S_c$ of the cluster gas increases with
$\sigma$ - at higher $\sigma$ the clusters will be free to undergo larger vibrations, leading to a higher vibrational entropy.

While this is a generic feature, $T_c$ otherwise shows strong target dependence. Fig.\ \ref{fig:TargetCvs} shows the heat capacity
$C_v$ as a function of temperature for each of the target structures at $\sigma = 0.45$. The peaks in the $C_v$ curves correspond to the
cluster-monomer gas transition, and we define $T_c$ as the temperature at which $C_v$ is a maximum. It is important to note that these plots
and the resulting values of $T_c$ were obtained from small systems containing only sufficient particles to form one target cluster. Hence
we expect significant finite size effects, and the values for large systems may be somewhat different.

There are two main factors determining the value of $T_c$. Firstly, the number of patches per particle, $m$ is important, since the more patches
are present the stabler the cluster gas is. The value of $n$, the number of particles per cluster, has the converse effect. When $N$ particles
assemble into $N/n$ clusters there is a reduction in the effective number of particles in the system, and hence in the translational entropy. Hence,
an estimate for $T_c$ can be made using

\begin{equation}T_c = \frac{\Delta E}{\Delta S} \propto \frac{m}{1-\frac{1}{n}} \end{equation}

%############################
(b) Target dependence

\begin{figure}
\includegraphics[width=8.4cm]{figures/allShapeCvsAnnotated.eps}
\caption{\label{fig:SizeHeatPlots}(Colour Online)}
\end{figure}

$T_c$ strongly depends on target.

2 factors:

1) $m$ - no. of particles, $V_c$ roughly equals $\frac{\epsilon}{2}$

2) Loss of translational entropy, $N$ particles -> $N/n$ clusters, i.e. there is a reduction in the effective number of particles.

We can use a very simple model to predict the order of transitions, neglecting vibrations and the mass dependence of 'effective' particles.

\begin{equation}T_c = \frac{\Delta E}{\Delta S} \propto \frac{m}{1-\frac{1}{m}}\end{equation}

This correctly predicts the order of the transitions, if not their exact relative values.

Order: [Make a little chart here perhaps?] Icos, oct, tet, cube, dodec.

\begin{center}
  \begin{tabular}{ l | c || r | }
    \hline
    Icosahedra \\ \hline
    Octahedra \\ \hline
    Tetrahedra \\ \hline
    Cubes \\ \hline
    Dodecahedra \\
    \hline
  \end{tabular}
\end{center}

%##############################

Note that the $C_v$ plot for dodecahedra exhibits a distinct shoulder, which turns out to be very important. We shall return to this
in Section\ \ref{sec}.

Lower values of $T_c$ make assembly more difficult, by decreasing the potential range of $T$ over which assembly might occur. It is noteworthy
that cubes and dodecahedra have the lowest values of $T_c$, and are the hardest of our targets to assemble.

\subsection{The cluster-liquid transition}
\label{sec:Aggregates}

For sufficiently large values of $\sigma$ we see formation of aggregates associated with liquid-vapour coexistence. We can see this for all target
structures on Fig.\ \ref{fig:SizeHeatPlots}. We also see aggregates at low $T$ and with narrower patches. The aggregates formed in these two regions
are in fact strongly related, with all properties varying smoothly between them. This is especially clear for the dodecahedral case,
where there is no region of cluster formation separating the two regions of aggregate formation. All of the aggregates are associated with
liquid-vapour separation.

Nevertheless, it is natural to divide the aggregates into two regions. In the first, at high $\sigma$, the aggregates are the stable state of the system.
In general these thermodynamic aggregates have some liquid-like characteristics, being fairly compact and undergoing constant rearrangement. The second
group of aggregates, those found at low temperature and lower $\sigma$, are metastable, with the target cluster being the true stable state. These kinetic
aggregates display more glassy properties, often forming extended networks which may percolate.

All of the aggregates tend to have irregular forms. Their surface tensions tend to be low, because so long as all patches are pointing inwards there is very
little surface energy. As a result there is little drive towards forming spherical droplets. At low values of
$\sigma$ the surface tension falls even further. This, combined with a reduced ability to rearrange, leads to the formation of increasingly irregular
and ramified aggregates at low $\sigma$. The energy of the aggregates tends to become less negative as well, since increasingly glassy systems are less
able to optimise their configurations to find lower energies. This effect is visible in Fig.\ \ref{fig:EnergyHeatPlots}, which shows the final energies
of simulations over a range of $T$ and $\sigma$ for each of the target structures.

\begin{figure}
\includegraphics[width=7.4cm]{figures/noTorEnergyPlots.eps}
\caption{\label{fig:EnergyHeatPlots}(Colour Online)}
\end{figure}

Aggregate structures also have a strong dependence on the target. The local structure within the aggregates tends to have a strong similarity to that
in the target cluster, so that very icosahedron-like groups are visible in the aggregates of icosahedron-forming particles, and so on. The notable
exception to this is the dodecahedron, where the wide splay angle allows the formation of a wide variety of structures, resulting in highly ramified
structures containing many hexagons and larger rings as well as pentagons. The size of aggregates may also depend on the target. In the case of tetrahedra,
and to some extent octahedra, the average size of the aggregates falls off strongly with decreasing $\sigma$, as seen in Fig.\ \ref{fig:SizeHeatPlots}.
The narrow splay angle promotes high curvature and hence small clusters, although at higher $\sigma$ the patch positions have less control over aggregate
structure and larger aggregates are still seen.  

The extent of the regions in which aggregates are found is important, as it sets boundaries on the region in which the target structure may successfully
be formed. Here we will deal with the factors affecting the region in which liquid aggregates are stable, while in Section \ref{sec:KineticAggs}
we will return to consider those determining the extent of the kinetic aggregate region.

It is useful to define some terms to describe the position of the transition from the cluster phase to the liquid at high $\sigma$. We define
$\sigma_c^{max}$ and $\sigma_l^{min}$ as the maximum value of $\sigma$ at which the cluster phase is found (which will be at zero temperature) and the
minimum value of $\sigma$ at which the liquid phase is found respectively. These two points are connected by the transition line $T_{cl}$. $\sigma_c^{max}$
depends entirely on the energies of the two states since it is found at zero $T$, while the gradient of $T_{cl}$ depends on the entropy difference between
the two states.

Fig.\ \ref{fig:EnergyHeatPlots} shows that the dodecahedra- and cube-forming particles are the most able to satisfy their patchy interactions in liquid
droplets, and attain the greatest bonding energies. The low number of patches per particle $m$ and the wide angles between patches result in fewer
constraints of the liquid structure, making it easier to optimise patch-patch interactions. We can see a clear
correlation between the angles between patches and the energy of the liquid. For the tetrahedra, octahedra and icosahedra, with angles between patches of
$60^\circ$, the particles would be required to pack very closely together in order to satisfy all their bonds, and in fact we see relatively low
bonding energies
per patch. As we increase the angle though, to $90^\circ$ for cubes and then $108^\circ$ for dodecahedra, the energy per patch steadily increases.

The entropy of the liquid aggregates might be expected to be correlated with the energy. When there are fewer constraints it is not only easier to satisfy
patch-patch interactions in a disordered configuration, but there are more ways of doing so, leading to a higher configurational entropy. We would therefore
expect particularly high aggregate entropies for the dodecahedra and cubes. The size of the clusters is also important since, on forming an aggregate of
a given size, more translational degrees of freedom are lost by smaller clusters.

The shallower the gradient of $T_{cl}$ the greater the corresponding entropy differences, with negative gradients indicating that the aggregate is more
entropically favourable. Unfortunately there is no clear approach for finding precise values for $T_{cl}$, but we can get an impression from
Fig.\ \ref{fig:SizeHeatPlots}. From those plots we can see that $T_{cl}$ has a negative gradient for octahedra, icosahedra and probably cubes,
indicating that the liquid has a higher entropy than the cluster gas, and hence that the disorder in the liquid more than compensates for the loss of
translational degrees of freedom on forming a large aggregate. The gradient is less clear for tetrahedra and dodecahedra, although we would expect more
closely similar entropies for tetrahedra, as the smallest of the targets. Later (see Sec.\ \ref{sec:Dodecahedra}) we show that $T_{cl}$ has a far shallower
slope for dodecahedra, such that $\sigma_l^{min}$ lies towards the left of the diagram. This is consistent with our expectations based on the dodecahedron's
low aggregate energy and large cluster size.

%% To conclude...

\subsubsection{3) $A_c < A_l, A_c < A_{g, m?}$, i.e. target most stable}

Er, what?

\subsection{Mechanisms of assembly and misassembly}

Here we consider the mechanics of assembly under three different temperature regimes in which the target clusters are stable. Note that this discussion
assumes a moderate value of $\sigma$; at very low $\sigma$ particles rarely feel each others' attractions, while at high $\sigma$ the stable state
is a liquid rather than a cluster gas.

\subsubsection{$T_c T > T_l$}

In this regime the target cluster is the only structure stable with respect to the monomer gas. Therefore assembly occurs by addition of monomers. There
are no competing states.

Note that as this region is close to $T_c$, the yield may be constrained by the finite width of the cluster-monomer transition, i.e. the equilibrium
yield of clusters may be significantly less than $100\%$.

In this region the speed of assembly is related to the size of the nucleation barrier and to the value of $\sigma$, with the dynamics becoming slower
with decreasing $\sigma$ due to the comparative rarity of bonding events with narrow patches.

The lack of competition in this regime raises the question of why no dodecahedra are formed by this mechanism. In fact, this region is only available
when $\sigma < \sigma^{min}_l$, and for dodecahedra $\sigma^{min}_l$ is very low (around $number$). For $\sigma < \sigma^{min}_l$, assembly is inhibited
by slow dynamics and a very low value of $T_c$, such that there is in fact no region of parameter space where reasonably fast assembly of dodecahedra
can occur.

\subsubsection{$T^{agg}_g < T < T_l$}

In this regime both the target clusters and large liquid aggregates are stable with respect to the monomer gas, leading to competition. However,
mobility in the liquid droplets allows them to rearrange to form target clusters. Once complete clusters are formed all their patches point inwards,
such that they experience very little attractive interaction with other particles. They are therefore able to ``bud off'', separating from the rest
of the liquid droplet in which they formed.

[Figure: Average cluster size against time for each target structure, at optimal conditions]

The optimal conditions for assembly are found in this regime, at a point in parameter space where the assembly mechanisms of monomer addition and liquid 
aggregation followed by budding both operate. Fig.\ \ref{fig:SizesVTime} shows the average cluster size against time under optimal assembly conditions for
each of the target structures.

\subsubsection{$T < T^{agg}_g$}
\label{sec:KineticAggs}

$T^{agg}_g$ is the temperature below which large clusters become unable to rearrange to form target over the timescale of the simulation (and is therefore
dependent on timescale). It is reminiscent of a glass transition temperature. Below $T^{agg}_g$ generally only large aggregates are formed, often forming
extended ramified networks which may percolate.

$T^{agg}_g$ is dependent on target geometry for a number of reasons. The first factor is the degree of order in the aggregates.
Fig\ \ref{fig:aggregatePics} shows kinetic aggregates for target structures of cubes, icosahedra and dodecahedra. In the icosahedral case,
the aggregates themselves contain a considerable degree of icosahedron-like structure, and only a relatively small amount of rearrangement is
required for the aggregates to form complete clusters. For the cubic case a greater degree of rearrangement is required, while in the dodecahedral
case the aggregates have very little dodecahedron-like structure at all, and would require almost total rearrangement to form the target clusters. 

A second factor is the size of the aggregates. In general a smaller splay angle in the assembling particles [make sure
have defined this!] tends to lead to a higher surface curvature and hence to smaller aggregates.
In Fig.\ \ref{fig:SizeHeatPlots} we can see that the low-$T$ aggregates formed with tetrahedral and octahedra target structures are much smaller than for
the other shapes. We would expect that smaller aggregates are easier to rearrange, [and indeed for the tetrahedron especially successful assembly is
observed down to low $T$. *This is a bit problematic because the octahedral assembly doesn't go particularly low, and the tet assembly is explained
by a later point, chance assembly]

A final factor is the size of the target structures themselves. Less rearrangement is needed to form a smaller target. Further, for smaller targets,
the growth of aggregates is less likely to deviate significantly from the correct assembly pathway. For tetrahedra, a some degree of successful
assembly is observed even at very low $T$ - simply by a chance, a significant proportion of groups of four particles coming together will assemble
into tetrahedra. However, even for octahedra this effect becomes negligible.

%% Consequences

\subsubsection{General features}

\begin{figure}
\includegraphics[width=8.4cm]{figures/allShapeYieldsWTime.eps}
\caption{\label{fig:AllShapeYieldsWTime}(Colour Online) The numbers of
(a) tetrahedra, (b) octahedra, (c) cubes, and (d) icosahedra formed after different numbers of simulation steps, as a function of temperature $T$,
in simulations of systems of 1200 particles at $\sigma=0.45$ and a number density of
$0.15\,\sigma_{\rm LJ}^{-3}$.}
\end{figure}

Fig.\ \ref{fig:AllShapeYieldsWTime} shows the yields of each of the target structures as a function of $T$ and time, at $\sigma = 0.45$. These plots show a
number of interesting features.

Firstly, in each case the maximum yield is obtained at high temperatures, close to $T_c$. Here the dynamics are fast, and since liquid-like aggregates
are not stable with respect to the monomer gas at these temperatures, the clusters face no competition. However, at very high $T$ the yield is limited
to the thermodynamic yield, which may be significantly less than $100\%$ around $T_c$.

At lower temperatures the yield falls off, until around $T^{agg}_g$ it becomes severely limited. For target structures of tetrahedra, cubes and icosahedra
rather sharp cut-offs are seen at low $T$, indicating approximate values for $T^{agg}_g$. In the case of the tetrahedra the yield does not fall to zero
at low $T$, because as discussed above some yield of tetrahedra is expected at low temperature by chance assembly with no rearrangement required.
In contrast, the plots for octahedra show an approximately linear decay in yields with decreasing temperature. The reasons for this are unclear.

\subsubsection{The difficulty of assembling dodecahedra}

Of the target structures examined in this paper, the dodecahedron is exceptional in the difficulty of its assembly. While the other challenging target, 
the cube, shows a restricted
region in which assembly occurs, it is nevertheless easy to obtain high yields by conducting long simulations with carefully chosen parameters. By contrast
no dodecahedra are ever formed under any set of conditions. Indeed, as we will discuss later, even strongly biased simulations which make use of umbrella
sampling to direct the system towards forming a dodecahedron are not successful. This apparently anomalous behaviour motivates more detailed study.

In general, because of the low energy of the target structures, there exists for each target at moderate $\sigma$ a temperature range over which the
target is the only stable structure with regard to the monomer gas, so that aggregates can be formed only transiently and cannot compete with the target.
One might expect that in this region at least, dodecahedra should be able to assemble. However, this is in fact not the case, since over most of the
range of $\sigma$ this region does not exist. Fig.\ \ref{fig:AllShapeYieldsWTime} shows a plot of the heat capacity $C_v$ along with a plot of the
average cluster size for a 

 However, as was mentioned in Section \ref{sec:Aggregates}, the
entropy of the dodecahedral aggregates is very favourable compared to the dodecahedra themselves, leading to a very shallow slope in $T_{cl}$. Although
we cannot precisely measure $T_{cl}$, Fig.\ \ref{fig:DodecaTclAndTc} gives some indication of its position. The higher line shows the temperature at
which the average cluster size in the box has a value of $10.5$ (the midpoint of $1$ and $20$), which approximates $T_{cl}$, while the lower line
is $T_c$. The two lines do not meet until around $\sigma = %%what?$. So, for all values of $\sigma$ higher than this, there is no region in which only
dodecahedra are stable. At lower values of $\sigma$ the assembly process is so hindered by the narrow patches that assembly is unable to occur.
%% Mention Fig.\ \ref{fig:AllShapeYieldsWTime} first.

\begin{figure}
\includegraphics[width=8.4cm]{figures/dodecaCvAndAvSize.eps}
\caption{\label{fig:AllShapeYieldsWTime}(Colour Online)}
\end{figure}

We expected that by reducing system to 20 particles, would restrict liquid state such that $A_{20 (l)} << A_{20 (target)}$, but no.

[Figure: $T_c$ and $T$ where av size is 10.5 as functions of $\sigma$]

[Figure: Picture of 20-atom liquid cluster]

First figure shows that liquid is stable to higher $T$ than cluster over a wide range of $\sigma$. The point where the two lines join provides
and upper bound on $\sigma^{min}_l$ (we expect relatively strong finite size effects on $A_{liquid}$, although not as bad as might generally
be expected for such a small liquid droplet because of the low surface tension).

Q: Extend $T_c$ measurements for other sizes into liquid region to see when molten clusters appear => handle on $T_{cl}$? (Perhaps unfair as greater
finite size effects for $n$ small. Or perhaps do it for roughly 20 particles in each case, i.e. 5 tetrahedra, 3 octahedra and so on; would need to
modify bond count order parameter so it was the total bond count rather than the number of bonds in the largest cluster.)

%% Non-energy landscape points:
% Dodecahedra really don't form well
% Unlike other shapes, there is no region where only clusters are stable. Hence assembly must always go via rearrangement of an aggregate.
% At higher temperatures, the aggregate is more stable anyway
% At lower temperatures, rearrangement is very slow
% Even slow rearrangement aside, succesful rearrangement seems really difficult
% This is probably because of geometric factors:
% - Many incorrect shapes are possible, while satisfying bonds well. e.g. interpenetration
% - Pentagons are easily deformed into hexagons or larger polygons without bond-angle problems (true for other shapes?)

In order to more fully understand the difficulty of assembling dodecahedra we attempted to obtain equilibrium thermodynamic data for a system of
20 particles, sufficient to assemble one dodecahedron. This was achieved only with considerable difficulty - indeed, the pathological difficulty
of forming a single dodecahedron even using biasing techniques underlines the extreme unlikeliness of ever forming one under ordinary dynamics, and
implies the presence of a large kinetic barrier to assembly.

\begin{figure}
\includegraphics[width=8.4cm]{figures/landscapes.eps}
\caption{\label{fig:Landscapes}(Colour Online)}
\end{figure}

Success was eventually achieved by combining two-dimensional umbrella sampling (using the number of bonds and the number of pentagons in the largest
cluster as order parameters) with Hamiltonian exchange. Even with heavily biased umbrella sampling no dodecahedra were ever formed using the model
potential that we have used throughout the rest of this paper, and it was found necessary to couple the simulation to an array of simulations which made use
of an additional constraint in the potential to encourage assembly (see Section \ref{sec:HamiltonianExchange}). The resultant thermodynamic data enabled
the plotting of the free energy landscapes shown in Fig.\ \ref{fig:Landscapes}. Note that these simulations were for systems large enough to form only a
single cluster (12 or 20 particles), and we were not able to obtain data for larger systems. As a result we expect finite size effects to have a significant
effect, and the details of the plots will be different than for bulk systems. However, we expect the main features of the landscapes will be unaffected.

The first pair of plots, Fig.\ \ref{fig:Landscapes}(a) and (b), show the free energy as a function of the number of particles and the number of correct
polygons (triangles for the icosahedron and pentagons for the dodecahedron) in the largest cluster in the system. In Fig.\ \ref{fig:Landscapes}(a), showing
the landscape for an icosahedron, we can see clear free energy minima in two regions, representing the monomer gas and the icosahedron. Note that the
monomer gas will typically contain at least one cluster of size three or four, which explains the position of the minimum. Between the two minima is a
region of higher free energy, representing a transition state. Once the transition state is passed, the shape of the landscape directs the system
towards forming an icosahedron, facilitating assembly (note that the high free energy state corresponding to a triangle count of 19 is simply represents
the fact that on forming the final bond, a nascent icosahedron transitions in one step from having 18 to 20 triangles, and does not represent a barrier).

By contrast in Fig.\ \ref{fig:Landscapes}(b), showing the landscape for a dodecahedron, there is little structure. There is no barrier to the formation
of a 20-particle cluster, so the system will tend to rapidly aggregate, and will then experience little drive to form a dodecahedron. Indeed, if the system
follows the steepest downhill path it will tend to form a single aggregate containing few pentagons, so that significant rearrangement would be needed
to form a dodecahedron.

The second pair of plots, Fig.\\ref{fig:Landscapes}(c) and (d), show the free energy as a function of two different parameters, the configurational energy
and a radial disorder parameter. The latter is in fact calculated as the standard deviation in the distance of particles from the centre of mass, and
has a value of zero for particles lying on a spherical shell (such as the two target structures). Again, the plot for the icosahedron shows two valleys,
correspondingto a monomer gas and to an icosahedron, separated by a transition state of higher free energy. Importantly, we see that in order to
obtain lower energies the system is forced to steadily reduce its maximum radial disorder, such that it is directed into the icosahedral state. We also
see that the free energy is a strong function of the radial order parameter.

By contrast, the free energy of the dodecahedral system appears to be almost entirely a function of the configurational energy, with very little
dependence on the radial disorder. Energies close to theat of the global minimum can be reached by highly disordered states. As a result the system is
undirected, and experiences very little pressure towards forming an ordered cluster.

However, neither of the dodecahedral plots fully explain the situation. From these diagrams one might expect the formation of dodecahedra to be rare,
but not essentially impossible. It appears that there is some large kinetic barrier preventing the formation of dodecahedra, which is not visible in these
diagrams due to the choice of axes. Despite some effort, we were not able to find a combination of order parameters which would show this barrier
explicitly.

%% A little conclusion?

\section{Conclusions}

- Yes, strong target dependence.

- Suggests limits to size and complexity of targets of monodisperse self-assembly without more constraints in the potential.

- \emph{Problem} - as $n$ increases, competition between target and disordered aggregates increasingly favours aggregates.

- In addition to that generic size dependence, there are specific geometry dependences. Easier if number of interactions is
  greater e.g. icos with $m = 5$ is fine, while with $m = 3$ cubes struggle and dodecahedra are impossible.
(Shells rather than cages are better).

Ways round this?

- Increase specificity of potential

- In particular if have torsional potential i.e.\ interactions which are specific both in direction and relative orientation, massively reduce
the number of alternative configurations.

- Protein-protein interactions have this kind of specificity

- Examine in accompanying paper

\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}
