# Marijke's project

###### tags: `bachelors projects`

---

### Owners (the only one with the permission to edit the main text, others can comment)

Marijke, Laura, Marjolein

---

# ModSim Exercises
## Estimating $\pi$
### Direct sampling
Using direct sampling I found the following values.

### Markov chain
Using the Markov chain method I found the following values.

## Random Walk
### 2D, no PBC
A simple random walk in 2D without boundary conditions

### 2D, with PBC
Again in 2D, but now with periodic boundary conditions. The box is 1x1.
The mean square displacement is roughly linear.

### 3D, with PBC
Now in 3D, with periodic boundary conditions. The box is 1x1x1.

## Monte Carlo Simulation of (Pseudo) Hard Spheres in the NVT Ensemble
### Hard Spheres
Maximum volume fractions
For cubic:
$$\frac{8 \cdot \frac{1}{8} \cdot \frac{4}{3} \pi r^3}{(2r)^3} = \frac{\pi}{6} \approx 0.5235$$
For FCC:
$$ \frac{(8 \cdot \frac{1}{8} + 6 \frac{1}{2}) \pi r^3}{(\sqrt{8}r)^3} = \frac{\pi}{3\sqrt{2}} \approx 0.7405$$

I finished the read_data() and move_particles() function of the provided code.

Starting from fcc.dat, mc_steps = 10000, output_steps = 100, diameter = 1.0, delta = 0.1:
packing fraction = 0.6 (move acceptance around 10%):

packing fraction = 0.5 (move acceptance around 28%-30%):

packing fraction = 0.4 (move acceptance 43%-45%):

I added a function that allows volume changes, so it can do simulations in the NPT ensemble.

To make the code faster, I implemented cell lists. A function make_cell_list() makes a cell lists. The number of cells in every direction is the length of the box divided by the cutoff radius, rounded down. This ensures that the length of the cell is at least the cutoff radius.
The function update_cell_list updates the cell lists when a particle moves to a different cell. It does this by removing the particle from the old cell list (by moving all the particle numbers after this particle to the index before it and reducing the length by 1) and adds the particle to the new cell list (at the end).

### Pseudo-Hard Spheres (WCA)

I rewrote the code so it uses a Lennard-Jones potential with a cutoff (truncated and shifted). With a cutoff of $2^\frac{1}{6}$, this is the WCA potential.

I wrote functions to calculate the virial pressure and energy of the system. I checked the function of the virial pressure with NPT simulations.

For a WCA potential in the NVT ensemble, I ran simulations for N=256 with 20.000 mc steps (29.000 for $\rho = 0.7$, as it took longer to reach equilibrium), measuring the average (virial) pressure and energy during the last 19.000 steps and got the following results.


I also did 3 runs for N=2048 to see if it made a difference (it does not).

# NVT simulations of the Modified Lennard-Jones potential (Broughton-Gilmer)
The Broughton-Gilmer potential is given by (Davidchack & Laird)
$$ \phi_{BG} =
\begin{cases}
4\epsilon \left[ (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^{6} \right] + C_1
& r \leq 2.3\sigma \\
C_2 (\frac{\sigma}{r})^{12} + C_3 (\frac{\sigma}{r})^{6}
+ C_4 (\frac{r}{\sigma})^2 + C_5
& 2.3\sigma < r < 2.5\sigma \\
0
& 2.5 \sigma \leq r
\end{cases} $$
where
$$\begin{align*}
C_1 &= 0.016132\epsilon \\
C_2 &= 3136.6\epsilon \\
C_3 &= -68.069\epsilon \\
C_4 &= -0.083312\epsilon \\
C_5 &= 0.74689\epsilon.
\end{align*}$$

According to Broughton-Gilmer (1986), the triple point is at $T=0.617$.

I ran simulations of a cubic box with $N=864$ particles at $T=0.3$, $T=0.617$ and $T=1.0$. The paper by Davidchack and Laird (2003) has a table with coexistence conditions for $T=0.617$ and $T=1.0$ (and $T=1.5$), so I plotted them too.


Close up of coexistence region with also the coexistence pressure:

Snapshot of long box (6144 particles) with $\rho_{init} = 0.89$ at $T=0.617$ and then increasing the $z$-direction by a factor 1.04.

Snapshot of long box (6144 particles) with $\rho_{init} = 0.96$ at $T=0.3$ and then increasing the $z$-direction by a factor 2.00.

## Theory