Simulations of the Gyroid Phase: 3d periodic soft matter from repulsive particles
Our model system consists of uniaxial particles with a tapered shape, the degree of tapering is controlled by the parameter k? (fig. 1).
Figure 1
Our model system consists of uniaxial particles with a tapered shape, the degree of tapering is controlled by the parameter k? (fig. 1).
Our model system consists of uniaxial particles with a tapered shape, the degree of tapering is controlled by the parameter k? (fig. 1).
Figure 1
Our model system consists of uniaxial particles with a tapered shape, the degree of tapering is controlled by the parameter k? (fig. 1).
Figure 1
The interaction between a pair of particles is defined by a so called soft repulsive potential (black graph on figure 2a). This is obtained by shifting the Lennard-Jones potential (blue graph on figure 2a) by the well depth e and cutting it off at an interparticle distance r0. Thus the particles interact by repulsive forces only which increase steeply as the particles come into close contact. As a result the particles behave in a similar fashion to hard particles whose potential consists simply of an infinite step function (red line on fig 2a) at the contact distance s. The contact function defines the separation at which a pair of particles come into close contact for any arbitrary orientations and relative positions as shown in figure 2b.
Figure 2a
Figure 2b
Our model system as a whole consists of a box of several thousand which at the start of a simulation are arranged on a lattice (fig 3a) and given an initial Maxwellian velocity distribution. The molecular dynamics (MD) technque is used to advance the coordinates of the particles forward in time in a stepwise fashion. It essentially involves feeding the current particle coordinates along with the net forces they experience due to their neighbours into a discretised form of Newton’s equations of motion. The output is the new set of coordinates for the next time step.
Figures 3a, 3b and 3c
Periodic boundary conditions are used to ‘fake’ a bulk system in much the same way as a replicated unit cell is used to generate a crystal structure which is conceptually infinite in extent. So if for instance a particle happens to move across the right hand boundary of the simulation box (see figure 3b) its coordinates are adjusted so that it re-enters the box through the opposite left hand face. Also the minimum image convention is applied so that a particle on, say, the left hand side of the box will interact with replicas of particles on the right hand side.
Most of the simulations we carry out involve gradually compressing the simulation box to increase the number density ?*=N/V, where N is the number of particles and V is the box volume. By carrying out many such simulations for a range of particle shapes we map out the density-shape diagram for tapered particles of fixed aspect ratio (fig 4a).
Figure 4a
Fred Bames had previously carried out Monte Carlo simulations on systems of hard tapered particles and mapped out their phase diagram (fig 4b). Note that sometimes it is necessary to employ anisotropic rescaling of the box dimensions i.e. the three box dimensions are allowed to extend/contract by different amounts whilst maintaining constant volume. This allows the box to adapt to the periodicity of the phase forming inside it.
Figure 4b
However several important questions were left unanswered, in particular the nature of the central region of the phase diagram was ambiguous. It is characterised by locally ordered regions but it was unclear whether this was a true equilibrium state or some metastable glassy phase.
Molecular dynamics simulations of soft repulsive particles reveal a similar phase diagram (fig 5) to that obtained from the Monte Carlo simulations of hard particles with the same shape. In addition the MD data revealed that the domain ordered region was in fact fluid also, intriguingly, the k?=3.8 system appeared to exhibit periodicity.
Figure 5
However the size of the repeating structures in the k?=3.8 system is comparable with the size of the simulation box, in cases like this it is possible that the apparently regular structure observed in the particle arrangements is not a characteristic of the system itself but is an artifact of the periodic boundary conditions described earlier.
Therefore to investigate the system further a configuration consisting of 1250 particles was replicated eight times to generate one comprising 10000. This system was then run on in a box of fixed dimensions in the NVE ensemble (constant particle number, constant volume, constant energy) – this is probably about as simple a molecular dynamics simulation as it is possible to perform.
Figure 6a
The periodicity engendered by replication of the smaller system persisted for some time, however after a million time steps or so hexagonal ordering emmerged on two pairs of opposite faces of the simulation box.
Figure 6b
Evidently this system had, purely as a result of the shape of its constituent particles, formed a highly periodic structure.
However it was not clear from looking at the the faces of the box or even cross sections through it what precisely the structure was. To discover what was in the box we used cluster analysis. In this method a cluster is defined as a set of coordinates each of which is lies within a certain arbitrary cut off radius of at least one other member of the group. Any particle outside this cutoff distance must belong to a separate cluster.
Note that close inspection of one of the ‘blooms’ (fig 7a) that make up the hexagonal pattern on the face of the box in (fig 6b), reveals that the structure seems to consist of interconnected bilayers. In general the rounded end of a particle always sits near the rounded end of other particles. Therefore to improve our chances of dintinguishing separate clusters, the particle coordinates were subjected to a translation along their long axes by just over half the particle length (fig 7b).
Figures 7a and 7b
Thanks to this measure we were able to distinguish two more or less identical clusters of translated particle coordinates, moreover the clusters formed an interpenetrating double network structure as shown in figure 8a.
Figure 8a
Figure 8b
(reproduced from Zeng, G Ungar and M Imperor-Clerc, Nature Materials 4, 562 (2005)
This constitutes a bicontinuous cubic phase which we believe to be the Ia3d gyroid (see fig 8b). Analysis of the structure is ongoing and includes a collaboration with a mathematician, Gerd Scroeder, who is attempting to map our simulated configurations onto the gyroid minimal surface (the surface that divides the two interpenetrating networks).
One of his preliminary results is shown in figure 9.
Figure 9
We would also like to study the processes by which this phase forms. In our initial simulations of this system, the cubic phase was obtained by compression of an isotropic phase. We anticipate that just prior to the isotropic-cubic transition precursor structures are present, identifying and studying these precursors may lead to a deeper understanding of why tapered particles form the cubic gyroid phase. Further simulations have also revealed that the same cubic phase can emmerge from decompression of a smectic bilayer phase via pore formation in the bilayers (fig 10). Bilayer to cubic phase transitions occur in cells but they are very difficult to interogate experimentally. Our simulated smectic-cubic transition therefore may provide valuable clues for biologists to understand similar processes in living systems.
Figure 10