# Square-triangle direct simulation

###### tags: `SoftQC` `self-assembly`

---

---

### Owners (the only one with the permission to edit the main text)

Frank, Giuseppe, Etienne

---

---

## Background

Chiara is trying to perform self assembly of nano platelets of square and triangle shape.

We want to have some simple simulations to compare with.

The first idea is to consider hard squares and triangles in 2D. Then, we could add some edge to edge attraction between the polygons to mimick possible ligand interactions.

## Bibliography

("Simulation square triangle" on Scholar gives quite a lot of stuff.)

* B. Rubinstein and S. I. Ben-Abraham, Mater. Sci. Eng. A **294–296**, 418–420 (2000)
A toy model for self assembly of squares and triangles. Authors use ad-hoc rules for the assembly that are not really justified. They find that favoring square-triangle contacts prevents demixing and results in qualitatively more satisfying structures.
* X. Ye et al., Nat. Chem. **5**, 466–473 (2013)

* J. de Graaf et al., J. Chem. Phys., **137**, 214101 (2012)
Implementation of the separating axes method for overlap detection.

## MC code

I recycle my standard Monte Carlo base to add support for hard convex polygons.

### Overlap detection

Separating axes method is used to detect overlaps between polygons.

To detect overlaps, we consider a set of axis, and project each polygon onto those axis. If the intervals of the projections on one of the axis are disjoint, the polygons do not overlap.

For each type of particle, the normalised projection axis vectors in the reference frame of the particle, the angle of the projection axis with respect to horizontal, and the extent of the projection of the particle onto those axis are constants and thus computed once and for all at the beginning of the run.

The projection interval of a square or a triangle with an arbitrary angle projected onto an axis can be computed analytically.
Therefore, to detect overlap on a given axis of particle i:
* Get the projection interval of the particle i onto its axis (constant, pre-calculated)
* Rotate $\vec{r_{ij}}$ into the frame of particle i
* Project $\vec{r_{ij}}$ onto the axis ($p_{ij} = \vec{r_{ij}} . \vec{ax}$)
* Compute the projection of particle j $[p_{lo}, p_{hi}]$ onto the axis using analytical formula
* Check whether $[p_{lo} + p_{ij}, p_{hi} + p_{ij}]$ and the projection interval of particle i on the axis are disjoint.

The extent of the projection of particle j onto the axis can also be computed by projecting the vertices, but this is much slower.

### Optimisations

* Inline various functions (trial and benchmark).
* For each particle type, store the radii of the inscribed and outscribed circles. Before running the full-fledge overlap detection, check that the two particles are separated by more than the sum of their inscribed radii, ans less than the sum of their outscribed radii. The speed gain is system and density dependent, but can be very significant (about a factor 2 for 50/50 squares and triangles of edge length 1 at pressure 10).
* For each particle, store a rotation matrix $R_{-\theta}$ that allows to rotate $\vec{r_{ij}}$ in the particle's frame. This matrix only needs to be updated upon particle rotation, which saves cos and sin evaluations.

### Optimal rejection rate target

I did a rough analysis of the optimal rejection rate targets.

#### Translations and rotations

The following two figures display the MSD at last time divided by the computational time of the production run, as a function of rejection rate targets.
The system is composed of 100 hard squares and 100 hard triangles of equal side length. For all simulations, an equilibration run of $5\times 10^4$ MC cycles is first performed, during which the system is gradually (linear ramp) compressed to the target pressure, and amplitudes of the MC steps are optimised to match the target rejection rates. Then a production run of the same length is performed, and the MSD recorded. At high pressures, these simulations are too short to equilibrate the system, but this is just a crude analysis to comfort intuition.
The first figure is obtained with a system at pressure $10 k_B T /\sigma^2$, and the second at $25 k_B T / \sigma^2$.


Without much surprise, the usual rationale for hard disks/spheres still holds for those hard polygons: rejection rate of about $0.7$ are optimal. From this very crude study, it seems that pressure does not change the location of the sweet spot a lot, however the optimal parameters concentrate in a narrower region at high pressures.

#### Box scaling

We compress a system of 100 squares and 100 triangles of equal edge length by linearly increasing pressure from $5k_BT\sigma^{-2}$ to $10k_BT\sigma^{-2}$ with steps of $1k_BT\sigma^{-2}$, during 10^8 MC steps.

The two following figures display the MSD at last time divided by the simulation time, and the volume as a function of timestep respectively.


Again, the usual rational applies. A rejection rate of 0.9 seems optimal for box scaling moves.

## Litterature

Looking for EOS data to validate the code, I glimpsed at the world of 2D phase transitions in hard polygon systems. It's quite fascinating. Here are a few papers that caught my attention:

* Anderson et al., *Shape and Symmetry Determine Two-Dimensional Melting Transitions of Hard Regular Polygons*, Phys. Rev. X **7**, 21001 (2017)
* K.W. Wojciechowski and D. Frenkel, *Tetratic phase in the planar hard square system ?*, Comput Methods Sci Technol **10**, 235-255 (2004)
* Gantapara et al, *A novel chiral phase of achiral hard triangles and an entropy-driven demixing of enantiomers*, Soft Matter **11**, 8684-8691 (2015)

## First simulations (shooting blindly)

I have started some simulations at compositions $x_{sq} = 0.5$ and $x_{sq} = 0.302$ (QC composition). Since I don't know the freezing pressure, I compressed the systems to relatively high pressures, hopping to cross the transition.
During compression, the pressure is either increased linearly or exponentially.

### Linear vs exponential compression

[18/08/2021]

I tried increasing the pressure in a linear of exponential fashion for a system of 500 squares and 500 triangles, from reduced pressure $5$ to $50$. The compression runs last for $10^6$ MC cycles in both cases. Simulations of the same system were also performed with a compression runs of $10^7$ MC cycles.
During linear compression, the reduced pressure is increased by 1 every $2.2 \times 10^4$ MC cyles.
During exponential compression, the reduced pressure is multiplied by $1.01$ every $4.31 \times 10^3$ MC cycles.

The following figures display the MSD and density as a function of MC step during the compression run.

||||
|:-:|:-:|:-:|
|MSD during compression run of $10^6$ MC cycles.| Density during compression run of $10^6 MC cycles.| Density during compression run of $10^7$ MC cycles.|

Linear compression yields larger densities that exponential compression. Consistently, since the density is lower in the exponentially compressed system, particles move more freely and the MSD is larger.

This could allow for a better exploration of the possible structures in exponentially compressed systems, and hence produce better equilibrated crystals. However, in this particular case, the two final structures seem of comparable quality.

|||
|:-:|:-:|
|||
|Resulting structure for exponential compression after $10^6$ MC cycles (top) and $10^7$ MC cycles (bottom).| Resulting structure for linear compression after $10^6$ MC cycles (top) and $10^7$ MC cycles (bottom).|

The linear compression yields slightly higher densities, without a cost in structure quality, therefore I use it in the following. Note that this might be very much protocol-dependent (I did not play with the exponential compression factor for instance). We know that both compression schemes are "good enough", so this should not be a matter of worry yet.

### Simulation results at QC composition

[18/08/2021]

The following figures are the results of simulations at the QC composition, linearly compressed from reduced pressure 5 to 30, with steps of 1. I tried gradually increasing the compression run duration to improve the structure quality.

|||
|:-:|:-:|:-:|
|Compression run of $10^7$ MC cycles (~ a few hours of simulation).|Compression run of $10^8$ MC cycles (~ 1 day of simulation).|Compression run of $10^9$ MC cycles (~ 10 days of simulation).|

With the slowest compression, it some random tiling regions appear. It is not clear whether they span the whole box or not, and carry (quasi) long range order.
I tried to compute structure factors. I am not quite sure about the appropriate procedure. So far, I located a diffraction center at the centre of each tile.
Computing the structure factor for squares alone, triangles alone, or both, give fluid like results. Alternatively, one could locate the diffraction centers at the vertices of the tiles, but this would result in arbitrarily close points. Could that be a problem ?

||||
|:-:|:-:|:-:|
|Structure factor for diffraction centres located at the centre of all tiles after $10^7$ MC cycles of compression.|Structure factor for diffraction centres located at the centre of all tiles after $10^8$ MC cycles of compression.|Structure factor for diffraction centres located at the centre of all tiles after $10^9$ MC cycles of compression.|

Maybe looking at simpler quantities, such as the orientation distribution of the times, could provide better insights.

Although the eyes of faith can spot some QC-inspiring local strctures, the structures obtained thus far do not deserve the name of quasicrystals.

## Equation of state measurements

[18/08/2021]

Annealing often helps to improve structures. One need to know the location of the freezing transition to compress and expend the system multiple times in its vincinity.

I computed the equation of state at the QC composition to try locate the freezing transition.
The simulations start with an equilibration run at the target pressure (quench) of $10^7$ MC cycles. Then, a production run of the same length is performed.

The color of the points signals the equilibration status, as detected by the slope method (see miscelaneous page). Red points are probably not well equilibrated.


No kink is clearly visible in the equation of state. This might be caused by the poor equilibration of quite a lot of simulation points. Considering how complicated melting transitions are in 2D (in particular for hard particles, not to mention mixtures of them... expected to form QC ! Interesting question though.), I am not even sure that looking for a kink in the EOS is a good idea.

Improvement ideas:
* Use a proper compression instead of a quench.
* Run for longer to equilibrate better.