# Rebecca's project

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

---

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

Rebecca, Laura

---

# Bachelor project Chemistry
Subject: Direct coexistence simulations in NVT-ensemble
Start date: 13 november 2023
Report deadline: 25 january 2024
Final: 30 january 2024

---

# Getting started with C
### 2: Estimating pi
2.1: Estimate pi using the direct sampling method.
Results of mean square deviation voor N = 10, 100, 1000, 10000 are 0.014, 0.0013, 0.0002, 0.000014. The mean square deviation thus scales exponentially with N.

2.2: Estimate pi using the Markov Chain method.
A rejection rate of ~0.9 yields the smallest mean square deviation and thus yields the highest precision.

2.3 A random walk on a MxM lattice using the Markov Chain method. Below you see the result for a 5x5 matrix and 100000 steps. All lattice sites are visited more or less the same amount of times.


2.4 The transfer matrix of a 3x3 lattice gives the probability of a step from one position to another:

Correlation time can be calculated using the second eigenvector of the transfer matrix using t=-1/log(e2). The correlation time as a function of M is given in the plot below:

### 3: Random walk
3.2 Simple random walk in 2D


Plot of a simple random walk in 2D with 1000 steps and step length 1.

3.3 Mean square distance from starting point after N steps


Functional dependence of the mean square distance from starting point after N steps for 100 different N and 100 random walks per mean value.

3.4 Mean square distance from starting point after N steps with periodic boundary conditions

For 100 different N, the mean square distance is more or less the same for each number of steps N when periodic boundary conditions are applied.

3.5 Simple random walk in 3D

You cannot generate a random unit vector using angles, because this yiels way more data points at the poles of the sphere than at the rest of the sphere. I created random unit vectors by generating random x,y,z and normalizing the resulting vector.
The functional dependence of the mean square distance from the starting point on N in 3D is similar to that in 2D.

# Monte Carlo simulation of Hard spheres in the NVT-ensemble
### 4.1 Place N hard spheres on cubic lattice

4.2 Maximum volume fraction cubic lattice:

### 4.3 Create FCC structure on cubic lattice
Ideas used for the algorithm to fill the lattice:

Result for 3x3x3 fcc unit cells:

4.4 Maximum volume fraction fcc lattice:

### 4.5 Monte Carlo simulation
Using the code that is given, with added read_file() and move_particle() routines, Monte Carlo simulations of hard spheres in NVT are done.
I would say that the packing fraction at which the fcc configuration melts is roughly 0.5

# Project

For this project I will do simulations on long boxes containing both crystal and fluid. In the end I want to plot the equation of state and the directional forces (z-direction and xy-direction of the box).
To get this data, I need to implement a potential first. This potential then enables me to calculate the pressure. At this point I can make an equation of state. Then I need to generalize the virial theorem for pressure to calculate the pressure tensor, which will give me the directional forces.

### Add WCA potential to Monte Carlo simulation in NVT ensemle

Structure of code:
**read_data()**: reads data from (fcc) file with coordinates for n particles
**set_density()**: changes box dimensions and particle coordinates to match the set density
int **move_particle()** :
- picks random particle to try to move
- calculate particle energy
- move particle + apply periodic boundaries
- calculate particle energy after move
- Metropolis condition for accepting/rejecting moves

**double distance(int part, int i)**: calculates the distance between the particle to analyze and particle i.

**double particle_energy(int part)**: loops over all n particles (except p): calculates distance, calculates energy if distance <= r_cut.
$u(r_{ij} \le 2^{1/6}\sigma)=4\varepsilon[(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^{6}]+\varepsilon$
$u(r_{ij} > 2^{1/6}\sigma)=0$

**write_data(int step)**: writes particle coordinates after every 100 Monte Carlo cycles in .dat-files

**int main()**:
- Monte Carlo cycles: try n moves per cycle
- Every 100 Monte Carlo cycles: calculate acceptance rate, write particle coordinates

### Calculate virial measure of pressure
Add calculation of virial measure of pressure to code using the formula:
$$P = \rho k_BT + \frac{1}{3V}<\sum_{i=1}^{n-1} \sum_{j=i+1}^n f_{ij} \cdot r_{ij}>$$
The force on a particle can be calculated from the energy: $\vec{f}(r_{ij}) = -\nabla u(r_{ij})$
The derivative gives: $-\frac{du_{ij}}{r_{ij}}=\frac{24\varepsilon}{r_{ij}}[2(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^{6}]$
So $P = \rho k_BT + \frac{1}{3V}<\sum_{i=1}^{n-1} \sum_{j=i+1}^n r_{ij} \cdot -\frac{du_{ij}}{r_{ij}}>=\rho k_BT + \frac{1}{3V}<\sum_{i=1}^{n-1} \sum_{j=i+1}^n 24\varepsilon[2(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^{6}]>$

Because the WCA is always in the exact opposite direction of the force, the direction of the force vector is not taken into account.
The pressure is calculated every 100 Monte Carlo cycles. Then, at the end of the run, the average pressure of all those values is calculated.

In code:
**double measure_pressure(int)**: calculates the pressure using the above formula every 100 Monte Carlo cycles.
**write_pressure_data(double pressures[output_cycles])**: writes all pressure values from pressures[output_cycles] array in a .dat file. The pressures are given in eV.

### Equation of state for WCA potential
This equation of state was made using 41 datapoints. For each datapoint I simulated a system of 256 particles for 50000 Monte Carlo cycles.
The pressure values (and energy values) are similar to Laura's data (which was obtained using 2048 particles).

### Stress tensor
For the stress tensor, we use the formula:
$$\mathbf{P}=\rho k_B T\mathbf{I} + \frac{1}{3V} <\sum_{i=1}^{n-1} \sum_{j=i+1}^n \mathbf{f}_{ij} \times \mathbf{r}_{ij}>$$
Every component $\alpha\beta$ in the resulting tensor matrix can be calculated using:
$$P_{\alpha\beta} = \rho k_B T \delta_{\alpha\beta} + \frac{1}{V} <\sum_{i=1}^{n-1} \sum_{j=i+1}^n f_{ij}^{\alpha} r_{ij}^{\beta}>$$

The resulting tensor matrix looks like:
$$\mathbf{P} = \begin{pmatrix}
P_{xx} & P_{xy} & P_{xz}\\
P_{yx} & P_{yy} & P_{yz}\\
P_{zx} & P_{zy} & P_{zz}\\
\end{pmatrix}$$

The values are symmetric, so $P_{xy} = P_{yx}$, $P_{xz} = P_{zx}$, $P_{yz} = P_{zy}$.
The trace of the matrix, $P_{xx} + P_{yy} + P_{zz}$ is equal to $P + 2*\rho k_B T$.

The calculation of the stress tensor is added to the **measure_pressure()** routine, because we calculate the same component for both values: $24\varepsilon[2(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^{6}]$.

### Long box simulations
The simulations for this project will be done using long boxes, containing at least 1000 particles. To initiate a long box, with both a crystalline and fluid phase in the box, I made a long box with FCC-packing and added a bit of extra space in the z-direction so that the crystal can melt into there.
The overall density of the box shoul be within the coexistence region. From Laura's paper it followes that the coexistence region is between $\rho = 0.704$ and $\rho = 0.785$.
(https://pubs-aip-org.proxy.library.uu.nl/aip/jcp/article/134/13/134901/189273/Simulation-of-nucleation-in-almost-hard-sphere)

The long box simulations will be done several times for slightly deformed systems. To deform a system I slightly change the x- and y-dimensions of the box with value $dxy$ so that the lattice constant in this direction changes. The I change the z-dimension of the box so that the volume of the box is the same as before the deformation. Thus, if the lattice parameters in the xy-directions increases, the lattice parameter in the z-direction decreases and vice versa.