# Haadi's Project
###### tags: `bachelor's Projects`

## Project discription
The focus of this project will be elastic constants of hard polyhedral. Negative Poisson ratios are a relatively rare material property that is highly sought after. It refers to materials who when you squeeze them in one direction, instead of expanding in perpendicular directions, they contract. In a recent work (V. Belde, BONZ thesis) it was found – rather surprisingly -- that crystals of hard cubes have this property. In this project we will explore whether this is also the case in other hard polyhedral.

## 1 Benchmark Programming

Prior to delving into simulations of systems will negative Poisson rations, we will begin with working with Monte Carlo methods and their applications to cubic hard sphere systems.

### 1.1 Estimating Pi
We begin by writing a C program using direct sampling. The code runs 20 cycles for various number of trias N = 10, 100, 1000, 10000. The mean square deviation
$$
MSE = \left \langle \left( \frac{N_{hit}}{N} - \left \langle \frac{N_{hit}}{N} \right\rangle \right)^2\right \rangle
$$

is written into a txt file and then plotted in mathematica against N:


Next, using the the Markov Chain method, first beginning at a random position inside the square, selecting two random numbers $∆x$ and $∆y$ from a uniform distribution in [$−δ, δ$]. The new position becomes [$x + ∆x$, $y + ∆y$]. In the case the new position falls outside the square, we reject this trial move and keep the old position. We repeat this procedure until we have generated N trial moves. After 20 cycles at $N = 10000$, varying $δ ∈ [0L, 3L]$ in steps of $0.5 \space \delta$ with $L=1$. Plotting the MSE against $\delta$ in Mathematica:


Next, we plot the rejection rate $N_{fail}/N$ as a function of $N$:


Next, we use the Markov Chain method for a discrete lattice with $M× M$ lattice sites using a neighbour table. For this section, I programmed a random walk by limiting the number generation to 1 step: 

void StartTable(int neighborTable[M][M][4]) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < M; j++) {
            neighborTable[i][j][0] = (i - 1 + M) % M; // Up
            neighborTable[i][j][1] = (i + 1) % M;     // Down
            neighborTable[i][j][2] = (j - 1 + M) % M; // Left
            neighborTable[i][j][3] = (j + 1) % M;     // Right
        }
    }

For a lattice value of $M = 5$, the data seems to agree with the theoretical model of a random walk. The visitation probabilities of each adjacent cite ($p_{ij}$ is nonzero if $|i - j| = 1$ ), `lattice[x][y]/N`, for large cycles show that:

Site (0, 0): 0.0982
Site (0, 1): 0.0535
Site (0, 4): 0.0502
Site (1, 0): 0.0451
Site (1, 1): 0.0964
Site (1, 2): 0.0486
Site (2, 1): 0.0494
Site (2, 2): 0.0996
Site (2, 3): 0.0466
Site (3, 2): 0.0490
Site (3, 3): 0.1061
Site (3, 4): 0.0513
Site (4, 0): 0.0507
Site (4, 3): 0.0525
Site (4, 4): 0.1028

### 1.2 Generation of Hard spheres in Lattices

#### 1.2.1 Cubic Crystal Structure

#### Global Variables:

lattice parameter $L$
number of cycles/spheres $N = L^3$.
sphere diameter $D = 1.000$

#### Functions:

**void apply_periodic_boundary_conditions(double position[3])** adds L to each coordinate and then takes the remainder when divided by L

**void initialize_spheres(double spheres[][3])** takes a 2D array as an argument, where each row represents a sphere, and each column represents the x, y, and z coordinates of that sphere, runs a loops form $0$ to $N$ and populates each according to the format of "xyz.dat" randomly along $L$. Uses a do while loop to check for overlap

#### int main

Algorithm:

Create sphere arrays

Initialization of spheres:
Populate arrays with randomly generate spheres
Pack spheres i.e ensure they do not overlap
Apply Period boundary conditions

Write to outfile

#### Result

Output:

417
0 7.646000
0 7.646000
0 7.646000
4.481758 1.087631 3.333208 1.000000
6.826200 2.898099 6.000387 1.000000
7.512677 3.409715 7.063784 1.000000
...


##### Maximum Volume fraction

The maximum volume fraction you can obtain with a cubic lattice without creating overlaps would be:

$$
\eta = \frac{1 \times \frac{4}{3}\pi r^3}{(2r)^3} = \frac{1}{6}\pi \approx 0.52
$$

### 1.2.3 Face Centered Cubic Crystal Structure

#### Global Variables

Lattice unitcell length $a$

##### Functions

**write_to_file()** standard writing function which writes to "xyz.dat"

#### Main

Algorithm:

Open file
Assign box dimensions $L_x$, $L_y$ and $L_z$(hard-coded)
Calculate $N_x$, $N_y$, and $N_z$ based on corresponding box dimensions divided by $a$
Calculate $N$
Create a unit cell
Paste it throughout the box dimensions:

#### Result
Output

108
4.242641  4.242641  4.242641
0.000000 0.000000 0.000000 1.000000
0.707107 0.707107 0.000000 1.000000
0.707107 0.000000 0.707107 1.000000
0.000000 0.707107 0.707107 1.000000
0.000000 0.000000 1.414214 1.000000
0.707107 0.707107 1.414214 1.000000
0.707107 0.000000 2.121320 1.000000
0.000000 0.707107 2.121320 1.000000
...


## 1.3 Monte Carlo simulation of Hard spheres in the NVT ensemble

#### Global Variables (Apart from those provided)

Lenard-Jones/WCA Parameters: $\sigma$, $\epsilon$, $r_{cut} = r^{1/6}\sigma \approx 1.2225$

Temperature $T$

Boltzmann Constant $k$

#### Functions

**void readcoords()** iterates through input file and assigns values to $N$, $L$, and the coordinate locations of the spheres

**void writecoords()** writes $N$, $L$ and coordinate locations to .dat file

**int overlap(int part)** iterates through entire lattice and checks if the distance between the input particle and any other particle is less than 1

**double wca_potential(double rij)** takes as input the distance between two particles i and j and returns the WCA potential for that distance

**double potential_energy(int part)** calculates the potential energy contribution of a particle due to its neighbouring particles by iterating through the lattice and sees if a given particle in within the cutoff ditance $r_{cut} = r^{1/6}\sigma$ of the particle,if it is, then it is returned

**double total_energy()** sums over all particle contributions and returns half the resulting energy sum

**int move(int part)**

Generates a random displacement from the interval $[-\delta, \delta]$ fy or each dimensions, giving us small displacements $\delta x$, $\delta y$ and $\delta z$.

Assigns current values of sphere coordinates to variables `old_x`, `old_y`, and `old_z`

Calculates total energy of lattice **total_energy()**

Assigns new displacements to current sphere coordinates

Enforces periodic boundary conditions

If overlap detected, reject the move

Calculate new energy and find $\Delta E$

Metropolist Criterion:

If $\Delta E$ < 0 or genrand_real2() < $e^{-\Delta E/kT}$
......then accept move
Else reject move

#### Main

algrithm:

Initialize all functions

iterate over cycles:

Call the move function

Take snapshots at successive cycle points

Print simulation information

### Analytic Calculations for a perfect FCC lattice

The WCA potential $U_{\text{WCA}}(r)$, is given by:

\begin{cases} 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right] - U_{\text{shift}} & \text{if } r < r_{c} \\
0 & \text{if } r \geq r_{c}
\end{cases}

Where $U_{\text{shift}} = 4\epsilon \left[ \left( \frac{\sigma}{r_c} \right)^{12} - \left( \frac{\sigma}{r_c} \right)^6 \right]$. Within the current analysis, $\sigma = 1$, $\epsilon = 1$, $k = 1$ and $T = 1$. From here we can calculate for some values of r:

$U_{\text{WCA}}(0)$ is undefined or infinity
$U_{\text{WCA}}(0.5) = 16328.25$
$U_{\text{WCA}}(1) = 1$
$U_{\text{WCA}}(1.122) = 0$

The total energy for a perfect FCC lattice with N unit cells would be given by:

$$
U_{\text{total}} = \frac{12N}{2} \times U_{\text{WCA}}\left(\frac{a}{\sqrt{2}}\right)
$$

For $r = 0.5$, $U_{\text{total}} \sim 10^5$ and for $r = 1$, $U_{\text{total}}= 648$ for $N = 108$ and $a = \sqrt{2}$, so the simulation total energy should be fluctuating somewheere around these energies (?)

#### Results & Debugging

**Trial 1: With setdensity() enabled**

Input File:

108
4.242641  4.242641  4.242641
0.000000 0.000000 0.000000 1.000000
0.707107 0.707107 0.000000 1.000000
0.707107 0.000000 0.707107 1.000000
0.000000 0.707107 0.707107 1.000000
...

with $\delta = 0.1$
Terminal Output:

Box: (4.24, 4.24, 4.24), Npart: 108
Scaled box: (4.76, 4.76, 4.76), Npart: 108
Total Energy before Simulation: 0.000000
Test TPmoves 108000 Pmoves: 29103
Fraction of accepted particle moves = 0.269472
Total Energy after Simulation: 72.228727

Comments:
Initial energy shouldn't be 0, I'm assuming this might have something to do with setdensity function creating an interparticle distance greater than the cutoff of the function

Energy fluctuations are off by almost exactly order of $10$ to their analytical predictions:


**Trial 2: Same input file but setdensity() is disabled:**

Output:

Box: (4.24, 4.24, 4.24), Npart: 108
Total Energy before Simulation: 647.998853
Test TPmoves 108000 Pmoves: 0
Fraction of accepted particle moves = 0.000000
Total Energy after Simulation: 647.998853
PS C:\Users\Lenovo\Documents\Thesis\C>

Comments: Energy seems closer to analytically calculated value but no moves get accepted, very small and large values of $\delta$ do not change the acceptance rate. I believe this is due to overlap function being implemented in the move function, therefore resulting in all moves getting rejected as the input is a perfect fcc lattice? -> solution to implement overlap in main? Does the overlap count as part of the metropolis method? (Probably not is my guess)

**Trial 3: Intermediate box dimensions between trials 1 and 2**

Input:

108
4.587913  4.587913  4.587913
0.000000 0.000000 0.000000 1.000000
0.707107 0.707107 0.000000 1.000000
0.707107 0.000000 0.707107 1.000000
0.000000 0.707107 0.707107 1.000000
0.000000 0.000000 1.414214 1.000000
0.707107 0.707107 1.414214 1.000000
0.707107 0.000000 2.121320 1.000000
...

Note: $\delta = 0.07$

Output:

Box: (4.59, 4.59, 4.59), Npart: 108
Total Energy before Simulation: 449.999713
Test TPmoves 108000 Pmoves: 26923
Fraction of accepted particle moves = 0.249287
Total Energy after Simulation: 123.347874

Comments: There appears to be a significant drop as the simulation acclimates. I'm under the impression this is expected


Visualisation of the lattice after this:


Seems slightly more irrgular than before but hard to say

**Trial 4: Sanity checking using Felix's input file**
(Thanks Felix)

Identical results and output to trial 2

**Points to discuss:**
1) Role of overlap function in move function
2) Check input/output files
3) Output units for reduced LJ parameters

### Incorporating a Virial Equation for Pressure

Expressing virial equation for pressure in terms of intermolecular force $\sum_{i < j}{\text{f}(r)} = -\nabla u(r)$:

$$
P = \frac{Nk_BT}{V} + \frac{1}{3V} \left< \sum_{j<i}^{N} r_{ij} \times \frac{du(r_{ij})}{dr_{ij}} \right>
$$

Programmatically:
Q) Where would this be implemented? On page 62 of the manuscript it states to do this in the measure() function, would this still be the case?
1. Find gradient of $u(r_{ij})$
2. Take product with sum over $r_{ij}$
3. Implement the $\left< ...\right>$ by summing over cycles
4. Calculate P

Explicitly calculating the derivative of the particle potential, we obtain:

$$
\mathbf{f}_{ij} = f_{ij} \hat{\mathbf{r}}_{ij} = 48\epsilon \left[ \left( \frac{\sigma^{12}}{r_{ij}^{13}} \right) - \left( \frac{\sigma^6}{r_{ij}^7} \right) \right] \hat{\mathbf{r}}_{ij}
$$

For the current computational context, the scalar analogue of the above should suffice (?):

$$
-\frac{du(r_{ij})}{dr_{ij}} = 48\epsilon \left[ \left( \frac{\sigma^{12}}{r_{ij}^{13}} \right) - \left( \frac{\sigma^6}{r_{ij}^7} \right) \right]
$$

When calculating virial pressure in a simulation, according to the virial theorem, we will be computing a time average which would translate to an average over the sampled configurations during the simulation run. If the simulation is long enough and the system reaches thermal equilibrium, this time average can be considered an approximation of the ensemble average. Therefore the importance of the assignment of`initcycles` will be greater now (?).

#### Functions

**compute_force(rij)** which takes as input the distance between two particles and returns the above derivative

**measure()** iterates through the lattice and calculates the force per particle, returns the virial equation for pressure

#### Results

Beginning first with a density of $\rho = \frac{N}{V} = 1.1412$:


Next, for $\rho = 1.04496$:


Last, for $\rho = 0.45141$ (changing $\delta = 0.5$):


Just like energy, there seems to be a staring value which the pressure begins from at cycles = 0 which is also depdendent on $\rho$.

#### Points to discuss

Questions:

1. Graphs and if they are correct
2. Is the energy and pressure dependency on density correct

### Tensor Generalization of the Virial Pressure

To calculate a tensor generalisation of the virial equation for pressure, the equation would involve some object of the form $\left< \sum_{j<i}^{N} \mathbf{r}_{ij} \otimes \mathbf{f}_{ij} \right>$.

Physically what do the components of the two tensors mean?

The tensor analogue of the virial equation for pressure would take the form:

$$
\mathbf{P} = \frac{N k_B T}{V} \mathbf{I} + \frac{1}{V} \left< \sum_{j<i}^{N} \mathbf{r}_{ij} \otimes \mathbf{f}_{ij} \right>
$$

$$
\mathbf{P} = \begin{bmatrix}
P_{xx} & P_{xy} & P_{xz} \\
P_{yx} & P_{yy} & P_{yz} \\
P_{zx} & P_{zy} & P_{zz} \\
\end{bmatrix} \equiv\sigma
$$

If we express the above matrix alternatively as:

$$
P_{\alpha\beta} = \frac{Nk_BT}{V} \delta_{\alpha\beta} + \frac{1}{V} \left< \sum_{j<i}^{N} r_{ij}^\alpha f_{ij}^\beta \right>
$$

then each term component in the tensor product $P'_{\alpha\beta} = \sum_{j<i}^{N} r_{ij}^\alpha f_{ij}^\beta$ can be worked out as:

$$
P'_{xx} = r_x \cdot f_x \\
P'_{xy} = r_x \cdot f_y \\
P'_{xz} = r_x \cdot f_z \\
P'_{yx} = r_y \cdot f_x \\
P'_{yy} = r_y \cdot f_y \\
P'_{yz} = r_y \cdot f_z \\
P'_{zx} = r_z \cdot f_x \\
P'_{zy} = r_z \cdot f_y \\
P'_{zz} = r_z \cdot f_z \\
$$

#### Further Debugging

After debugging an error in the boundary conditions and programming an energy check, the corresponding graphs are:


#### Graphing an Equation of State:

For an initial attempt, I varied the setdensity() function on the 6.25 side input lattice

For $\rho = {1.4, 1.2, 1.3}$, $\delta = 0.07$
for $\rho = {1.0, 0.8,0.6}$, $\delta =0.13$
for $\rho = 0.4$, $\delta = 0.5$
for $\rho = 0.2$, $\delta = 1.3$
for $\rho = 0.15$, $\delta = 2.3$


Note: Here ${T, k_b, \epsilon, \sigma} = 1$

### Computing the Virial Pressure Tensor

#### Global Variables

double tensor_product[3][3] which holds the tensor product $P'_{\alpha\beta} = \sum_{j<i}^{N} r_{ij}^\alpha f_{ij}^\beta$

#### Ammended functions:

**compute_force(rij)** which takes as input the distance between two particles and returns the force as negative derivative of particle potential contribution

**measure()** iterates through the lattice and calculates the force per particle, then assigns the product of the force and component of distance to the appropriate matrix element (worked out explicitly above) then calls **print_tensor_product();**

**print_tensor_product();** Prints the contents of tensor_product[3][3] in terminal and writes it onto a file for external manipulation

#### Results

By taking multiple snapshots of the tensor product $P'_{\alpha\beta} = \sum_{j<i}^{N} r_{ij}^\alpha f_{ij}^\beta$, then taking the piecewise average and dividing all elements of the matrix by $3V$, our virial equation becomes:

$$
\mathbf{P'} = \begin{bmatrix}
74.5777 & 0.00600323 & 0.0889667 \\
0.00600323 & 74.6743 & 0.165513 \\
0.0889667 & 0.165513 & 74.5485 \\
\end{bmatrix}
$$

Adding to this the ideal gas term:

$$
\begin{bmatrix}
\rho & 0 & 0 \\
0 & \rho & 0 \\
0 & 0 & \rho \\
\end{bmatrix}
$$

For our current system $\rho = 1.41045$, the pressure tensor is:

$$
\mathbf{P} = \begin{bmatrix}
75.98815 & 0.00600323 & 0.0889667 \\
0.00600323 & 76.08475 & 0.165513 \\
0.0889667 & 0.165513 & 75.95895 \\
\end{bmatrix}
$$

#### Debugging Energy

After debugging the potential energy function after feedback, the energy per particle of the system, plotted against cycles towards the end of the simulation:


The above has a mean of $\mu = 6.96582$ and a standard deviation of $\sigma = 0.115065$ for a $\rho = 1.4000000000000041$. For Frank's data, the energy per particle for this density is $6.973752766925632$

#### Equation of State

By varying the density and calculating the corresponding pressure, the following plot displays the equation of state alongside Frank's data:


At lower values the simulation tends to match the provided data very accurately, but otherwise is slightly higher consistently

### Elastic Deformation

#### Pure Stretch

In the lattice initialization, a deformation function has been added:

**deform_position(double *x, double *y, double *z, double delta)** takes the positions of each lattice point, and applies the stretch deformation on each axis

The code now initially deforms the lattice, and each hard sphere position as it iterates in a fcc orientation

Next, we focus on a strectch with the corresponding deformation matrix:

$$
\begin{align}
\mathbf{D}_\text{StrechZ} &= \begin{pmatrix} 1 + \delta & 0 & 0 \\ 0 & 1 - \delta & 0 \\
0 & 0 & \frac{1}{1 - \delta^2} \end{pmatrix}
\end{align}
$$

Using $\epsilon = \frac{1}{2} \big(\mathbf{D}^T\mathbf{D} - 1\big)$, we find the corresponding strain
\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\frac{\delta ^2}{2}+\delta & 0 & 0 \\
0 & \frac{\delta ^2}{2}-\delta & 0 \\
0 & 0 & \frac{1}{2 \left(1-\delta ^2\right)^2}-\frac{1}{2} \\
\end{array}
\right)
\end{align}

We can relate the stress tensor of the strained system to the stress tensor before the deformation was applied by using:
\begin{align}
\sigma_{ij}' = \sigma_{ij} + B_{ijkl} \epsilon_{kl} + \mathcal{O}(\epsilon^2)
\end{align}

Where we have $\sigma_{ij}$ from our initial stress tensor calculations, $\sigma_{ij}'$ from calculating the stress tensor in deformation, and $\epsilon_{kl}$ from the above. $B_{ijkl}$ represents the elastic modulus tensor (or stiffness tensor), relating stress and strain in the material. The indices $ijkl$ suggest it's a fourth-order tensor in three-dimensional space, making it possess 81 components. However, due to symmetries, the independent components can be fewer.

>[name=Frank] The text below seems to be a bit confused. For the deformation matrix you show above, $\epsilon_{kl}$ is only nonzero for $k=l$. Take a close look at Eqs. 1.64 and 1.65 in Vera's thesis (and the text around it). You should be able to get $\mu^\prime = (B_{11}-B_{12})/2$ from these simulations, by looking at the diagonal components of the stress tensor. The off-diagonal components will come into play later, when we deform the system by shifting the boundary conditions.

As mentioned in (Velda, 2023) we can obtain $B_{44} = B_{1313} = B_{1331}$ if we take $ij = 13$. By exploiting the symmetry that only when $kl = 13, 31, 33$ the second term on the left-hand side will be non-zero. We know that $\sigma_{ij} = -P \delta_{ij}$, with $\delta_{ij}$ the kronecker delta.
Thus we have
\begin{align}
\sigma_{13}’ = 0 + B_{1313} \bigg(-\frac{\delta}{2}\bigg) + B_{1331} \bigg(-\frac{\delta}{2}\bigg) + B_{1333} \bigg(\frac{\delta^2}{2}\bigg)
\end{align}
Due to the symmetries of the cube (and thus also definitely of the sphere), $B_{1333} = 0$. This leaves us with
\begin{align}
\sigma_{13}' &= -B_{1313} \frac{\delta}{2} - B_{1331} \frac{\delta}{2}\\
&= -B_{44}\delta
\end{align}

>[name=Frank] Note that the cubic symmetry of interest in Vera's thesis is not the symmetry of the cube-shaped particle, but rather of the cubic crystal. You are also looking at a cubic crystal (face-centered cubic) so the same symmetries apply. But not because of the spherical symmetry of the particle: spherical particles can form non-cubic crystals!

Thus, we can find $B_{44}$ by taking taking $\sigma_{13}’$ from the stress tensor we get as a result of the simulation, and dividing it by the $\delta$ we chose.

Next, to get $B_{11}$ and $B_{12}$, we would need to relate relevant components of $B_{ijkll}$ to the bulk modulus $B$ which would require an equation of state:

#### Test Calculations:

For $\rho$ = 1.00000, the present simulation obtains the following pressure tensor:

$$
\begin{align}
\sigma &= \begin{pmatrix} 11.6107 & -0.0080132 & -0.0077402 \\ 0.0080132 & 11.5863 & 0.00337572 \\
-0.0077402 & 0.00337572 & 11.5989 \end{pmatrix}
\end{align}
$$

And the coressponding $\sigma_{ij}'$ for a stretch in the xy direction which preserves density:

$$
\begin{align}
\sigma' &= \begin{pmatrix} 11.4716 & 0.00726063 & -0.015653 \\ 0.00726063 & 11.6222 & 0.0160515 \\ -0.015653 & 0.0160515 & 11.5613 \end{pmatrix}
\end{align}
$$

**Q4:** Is anything changing?
>[name=Frank] Maybe! It looks like your 11 component dropped significantly, for example. You might need some better statistics... Let's talk about using the cluster at the next meeting. How long did the simulation run (in real-world time)?

The corresponding $\epsilon$ is:

$$
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
0.0050125 & 0 & 0 \\
0 & -0.0049875 & 0 \\
0 & 0 & 0.0000999975 \\
\end{array}
\right)
$$

## Determining Elasticity Tensor

For more convenient and effective handelling of deformation, seperate scripts to run each kind of deformation have been added.

### Pure Stretch

readcoords();
    boxx = boxx * (1 + delta);
    boxy = boxy * (1 - delta);
    boxz = boxz * (1.0 / (1 - delta * delta));
    int i;

    for (i = 0 ; i < N; i++)
    {
      r[i][0] = (1 + delta) * r[i][0]; 
      r[i][1] = (1 - delta) * r[i][1];
      r[i][2] = (1.0 / (1 - delta * delta)) * r[i][2];
    } 

    writecoords();

### Shear in $X$

The logical framework of this shear was applied to its $Y$ counterpart:

### Determining $B_{12}$ and $B_{11}$

Following the outline of calculations performed in Belde (2023), we can begin with:
\begin{align}
\mathbf{D}_\text{StrechZ} &= \begin{pmatrix} 1 + \delta & 0 & 0 \\ 0 & 1 - \delta & 0 \\
0 & 0 & \frac{1}{1 - \delta^2} \end{pmatrix}
\end{align}
Using $\epsilon = \frac{1}{2} \big(\mathbf{D}^T\mathbf{D} - 1\big)$, we find the corresponding strain
\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\frac{\delta ^2}{2}+\delta & 0 & 0 \\
0 & \frac{\delta ^2}{2}-\delta & 0 \\
0 & 0 & \frac{1}{2 \left(1-\delta ^2\right)^2}-\frac{1}{2} \\
\end{array}
\right)
\end{align}

Because we only use small $\delta$, we can reduce it to the first order in $\delta$. We then get:
\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\delta & 0 & 0 \\
0 & -\delta & 0 \\
0 & 0 & 0 \\
\end{array}
\right) + \mathcal{O}(\delta^2)
\end{align}

We again use the following relation:
\begin{align}
\sigma_{ij}' = \sigma_{ij} + B_{ijkl}\epsilon_{kl} + \mathcal{O}(\epsilon^2)
\end{align}

From this, we can derive:

\begin{align}
\frac{\text{d}(\sigma_{11}' - P)}{\text{d}\delta} = B_{11} - B_{12} + \mathcal{O}(\delta)
\end{align}

From the elastic constants for isotropic materials, the bulk modulus $B$ can also be defined as:

$B = \frac{1}{3}(B_{12} +2B){12})$

Expressing the bulk modulus in terms of equation of state:

$B = \frac{1}{3}(B_{12} +2B_{12}) = \rho \frac{dP}{d\rho}$

#### Equation of state interlude

Now, to fit the previously constructured equation of state, below are two condidate fits.

First, we use `curve_fit()` in `scipy`:


With the corresponding equation: $P = a\rho^3 +b\rho^2 +c \rho +d$ where $a = 102.09370879310924, b = -154.09055691274017, c = 70.75750354230874, d = -6.988246787631052$

**Q1:** I have a second fit equation I derived in an earlier project:


The fit equation is: $P = d + \frac{(a-d)}{(1 + (\rho/c)^b)}$ where $a = 1.654307, d = 52287540, c = 12.81897, b = 6.095766$.
**Q1.1:** Is this overfit?
**Q1.2:** Despite being more accurate to the data, making further computations with this expression will be taxing I assume. Finding an analytical derivative is possible, but I'm not sure whether other operations would be feasible. Is the better fit worth the extra effort in computations? Am I doing something wrong and is there another fit which might exist that I'm not trying?

#### Calculation Continuation

Now, we can obtain:

$$
B = \rho \frac{dP}{d\rho} = 3a\rho^3 +2b\rho^2 +c \rho
$$

where $a = 102.09370879310924, b = -154.09055691274017$ and $c = 70.75750354230874$. We can now combine our equation relating the elasticity tensor components and the bulk modulus with our strain-strain equation to get:

If we then combine and solve the two equations, we find:
\begin{align}
B_{11} &= \frac{1}{3}\bigg(3B + 2\frac{\text{d}(\sigma_{11}' - P)}{\text{d}\delta}\bigg)\\
B_{12} &= \frac{1}{3}\bigg(3B - \frac{\text{d}(\sigma_{11}' - P)}{\text{d}\delta}\bigg)
\end{align}

Noting that, $B=68.85751609615611$ and $\rho = 1.0000000$:


The values shouldn't fluctuate like this. The error I think would lie in the equations of analysis in my Python script of in the deformation itself.

Noting that it would be more appropriate (to mitigate accounting for phase transitions) for our equation of state fit to be such that $\rho >1$:


The corresponding equation is: $P = a\rho^3 +b\rho^2 +c \rho +d$ where $a = 366.1320765001868, b = -959.572638315137, c = 862.3180343658629, d = -257.37654591930146$. Therefore, $B =41.56898723614981$


### Determining $B_{44}$

For an extreme shear value of $\delta = 0.5$


$$\sigma_{13}'= -B_{1313} \frac{\delta}{2} - B_{1331} \frac{\delta}{2} = -B_{44}\delta$$


**Q1:** These values seem to be negative, which doesn't agree with Belde (2023) or Frenkel & Ladd (1987).


## EDMD Simulations

The provided code is an Event-Driven Molecular Dynamics (EDMD) simulation for a system of 2D squares (more specifically, 2D polygons). EDMD is an approach where one calculates the time of the next event (e.g., particle collision) and directly jumps to that event in the simulation. The code provides an overview of the setup and the evolution of such a simulation system.

### Code Breakdown

#### Definitions:
**MAXNEIGH:** Maximum number of neighbors a particle can have.
**EXTRAEVENTS:** Number of extra events (like writing, thermostat) to allocate space for.
**M_PI:** Mathematical constant π.

#### Variables:
**System Parameters:**
**N:** Number of particles.
**packingfraction:** The fraction of the system volume filled by particles.
**maxtime:** Maximum simulation time.
**makesnapshots, writeinterval, snapshotinterval:** Parameters for data output.
**stretchy, shiftxy:** Modify particle deformations.

#### Snapshot File

N
xsize ysize zsize offset: 
x y z corners angle 
...
// Repeat of the above

### Code Visualisation

Below is a visualiser written for the current configuration. Further improvements will be needed to encorporate visualisation of transformations


### Elastic Constants

When we consider a 2D lattice with square symmetry, the behavior is somewhat analogous to the 3D cubic symmetry but with fewer dimensions. In a general 2D case, the elastic constant tensor has 16 components, represented by $C_{ijkl}$ where each index $i,j,k,l$ can take values from the set ${ {1, 2} }$. This is because the elastic tensor is a rank-4 tensor, and in two dimensions, each index can be either 1 or 2. However, due to the symmetries of a square, the number of independent components is reduced:

**Symmetries:**

Reflection about one axis (e.g., x-axis): Components like $C_{1121}$ would be equivalent to $C_{1211}$ due to this symmetry. Other similar symmetries will reduce further independent components.

90° rotation: A square lattice has rotational symmetry every 90°. This symmetry will lead to equivalences like $C_{1111} = C_{2222}$ and $C_{1122} = C_{2211}$.

**Independent Components for 2D Square Lattice:**
After considering the reflection and rotation symmetries inherent to a square lattice, we can deduce that the only independent components are:

* $C_{1111}$ (equivalent to $C_{2222}$)
* $C_{1122}$ (equivalent to $C_{2211}$)
* $C_{1212}$ (or equivalently $C_{2112}$ due to symmetry)

Similar to the cubic lattice case, these independent components have physical interpretations:
* $C_{1111}$: Resistance to stretching or compressing along one axis (e.g., x-axis).
* $C_{1122}$: Coupling between stretching in one direction and compression in the perpendicular direction.
* $C_{1212}$: Shear modulus, representing resistance to shearing deformations.

Therefore, our analysis remains the same, with the independent nontrivial components of the elasticity tensor being $B_{11}$, $B_{12}$, and $B_{44}$

### Equation of State

The code prints out the stress tensor for a given simulation in row major form:

nan  nan  nan  nan  
P_xx P_xy P_yx Pyy

Noting first that the stress files output by the code are in units $kT/\sigma^2$, and internally the length $\sigma$ is given by `\sigma = \sqrt{2} [internal length unit]`. We can convert the units of the output and compare with [Frenkel & Wojciechowski (2004)](https://https://cmst.eu/wp-content/uploads/files/10.12921_cmst.2004.10.02.235-255_Wojciechowski.pdf):



Here, $B = 137.686235594173685964$ for $\rho = 0.87$:


$a = 35002.536972903996, b = -90481.31565854346, c = 78115.48858240536$

### $B_{11}$ and $B_{12}$

\begin{align}
\mathbf{D}_\text{StrechY} &= \begin{pmatrix} \frac{1}{1+\delta} & 0 \\ 0 & 1 + \delta
\end{pmatrix}
\end{align}

Using $\epsilon = \frac{1}{2} \big(\mathbf{D}^T\mathbf{D} - 1\big)$ we can obtain:

$$
\epsilon = \frac{1}{2}\begin{pmatrix} \frac{1}{(1+\delta)^2} -1 & 0 \\ 0 & (1 + \delta)^2 -1 \end{pmatrix}
$$

As $\delta <<1$, it follows that $\frac{1}{(1+\delta)^2} \approx 1-2\delta$ and $(1+\delta)^2 \approx 1+2\delta$:

$$
\epsilon \approx \begin{pmatrix} -\delta & 0 \\ 0 & \delta \end{pmatrix}
$$

Therefore:

$$
\sigma_{ij}' - \sigma_{ij} = B_{ij22}\delta - B_{ij11}\delta
$$

Again by using the bulk modulus $B = \frac{1}{2}\left(B_{11} + B_{12} \right) = \rho \frac{dP}{d\rho}$ (see Appendix below for derivation), and by setting $i,j = 2,2$ we can obtain:

$$
\frac{\sigma_{22}' - \sigma_{22}}{\delta} = B_{11} - B_{12}
$$

and then:

$$
B_{12} =-\frac{1}{2}\left(\frac{\sigma_{22}' + P}{\delta} - 2B\right)
$$

$$
B_{12} = \frac{1}{2}\left(\frac{\sigma_{22}-\sigma_{22}'}{\delta} +2B\right) = -\frac{1}{2}\left(\frac{\sigma_{22}' + P}{\delta} - 2B\right)
$$
$$
B_{11} = \frac{1}{2}\left(\frac{\sigma_{22}' + P}{\delta} + 2B\right)
$$

Now, for a 24 hour run:


### Determining $B_{44}$

Our deformation matrix is:

$$
\mathbf{D}_\text{ShiftXY} = \begin{pmatrix} 1 & -\delta \\ 0 & 1 \end{pmatrix}
$$

Therefore:

$$
\epsilon = \begin{pmatrix} 0 & -\frac{\delta}{2} \\ -\frac{\delta}{2} & \frac{\delta^2}{2} \end{pmatrix}
$$

As $\delta<<1$ we can approximate the above as:

$$
\epsilon \approx \begin{pmatrix} 0 & -\frac{\delta}{2} \\ -\frac{\delta}{2} & 0 \end{pmatrix}
$$

Resulting in:

$$
B_{1212} = -\frac{\sigma'_{12}}{\delta}
$$

## Sketching a Derivation of Constants

### Poisson Consants
To begin, we will consider the simplest case of 2D isotropic systems and then attempt to extend our analysis to our anisotropic planar hard squares system. For $d$-dimensional systems the Poisson’s ratio is simply related to $\mu$ and $B$ ([Tretiakov & Wojciechowski, 2005](https://onlinelibrary.wiley.com/doi/epdf/10.1002/pssb.200460390))by:

$$
\nu = \frac{dB - 2\mu}{d(d - 1)B + 2\mu}
$$

For $d = 2$:

$$
\nu = \frac{B-\mu}{B+\mu}
$$

We know from previous calculations that $B =\frac{1}{2}(B_{11} + B_{12})$. At this point it is important to note at this point that it has been proved ([Wojciechowski, 2002](https://www.osti.gov/etdeweb/servlets/purl/20350535)) that for 2D isotropic systems, the limit of the above are:

$$
-1 \leq \nu \leq 1
$$

To determine $\mu$, we can follow the outline of the derivation of specific Poisson constants ([Ting & Chen, 2004,](https://t.ly/uq9d7) [Dmitriev & Shigenari, 2001](https://www.researchgate.net/publication/243566894_Poisson_Ratio_Beyond_the_Limits_of_the_Elasticity_Theory)). The Poisson's ratio is the negative ratio of transverse strain to axial strain in a material subjected to axial stress. In an anisotropic material, this ratio can vary depending on the direction of the applied stress and the direction in which the strain is measured. For a material that exhibits anisotropic behavior, the Poisson's ratio in the direction $i$ due to stress applied in the direction $j$ is denoted as $\nu_{ij}$ and can be approximated by the ratio of the stiffness tensor components in those directions, assuming the anisotropy is sufficiently small.

Specifically, $\nu_{12}$, which is the Poisson's ratio for strain in direction 1 due to stress in direction 2, can be approximated as the negative ratio of the off-diagonal stiffness tensor component $B_{12}$ to the diagonal component $B_{11}$, representing the coupling between normal strains in direction 1 due to normal stresses in direction 2:
$$
\nu_{12} \approx -\frac{B_{12}}{B_{11}}
$$

Similarly, $\nu_{21}$, which is the Poisson's ratio for strain in direction 2 due to stress in direction 1, can be approximated as the negative ratio of the same off-diagonal component $B_{12} to the diagonal component $B_{22}$, representing the coupling between normal strains in direction 2 due to normal stresses in direction 1:

$$
\nu_{21} \approx -\frac{B_{12}}{B_{22}}
$$

From the above it follows that $\mu = \frac{(B_{11} + B_{12})^2}{2(B_{11} - B_{12})}$ (see Appendix).

For anisotropic auxetic squares:

$$
\nu_{\mathbf{nm}} = \frac{\sum_{ijkl} m_i m_j n_k n_l S_{ijkl}}{\sum_{ijkl} n_i n_j n_k n_l S_{ijkl}}
$$

Where $i,j,k,l\in \{1,2\}$. For $\mathbf{n} = (0,1)$ and $\mathbf{m} = (1,0)$:

$$
\nu_{\mathbf{nm}} = \frac{-S_{12}}{S_{11}}
$$

The above can be reproduced from (2.7) in [Grima et al.(2011)](https://www.jstor.org/stable/pdf/41059587.pdf) or (4.55) in [Daniel & Ishai (1994)](https://dl.icdst.org/pdfs/files3/ff750bc24ae06fc21cf3301c2b3fcb5a.pdf)

**Q1:** Can we use the relation below to obtain each successive poisson constant?:

$$
\begin{bmatrix}
\sigma_{xx} \\
\sigma_{yy} \\
\sigma_{xy}
\end{bmatrix}
= \frac{1}{1 - \nu_{xy} \nu_{yx}}
\begin{bmatrix}
E_x & \nu_{yx} E_x & 0 \\
\nu_{xy} E_y & E_y & 0 \\
0 & 0 & G_{xy}(1 - \nu_{xy}\nu_{yx})
\end{bmatrix}
\begin{bmatrix}
\varepsilon_{xx} \\
\varepsilon_{yy} \\
2\varepsilon_{xy}
\end{bmatrix}
$$

If yes, then we again end up with $\nu_{yx} = -\frac{B_{12}}{B_{11}}$

For $\mathbf{n} = (1,1)$ and $\mathbf{m} = (0,1)$:

$$
\frac{-S_{11}}{S_{11}+S_{12}+0.5S_{44}}
$$

**Q2:** To draw up a general formula, we need to determine a $n$ and $m$ in terms of rod angle $\theta$. The derivation in [Tokmakova, 2004](https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/pssb.200460389) doesn't justify why they used those specific values (or maybe I'm missing something). I tried $n = (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})$ and $m = (\frac{1}{\sqrt{2}} \sin{\theta}, - \frac{1}{\sqrt{2}}\sin{\theta})$ but the answer is

$$
\frac{(S_{11} +S_{12} - 0.5S_{44})\sin^2{\theta}}{S_{11} +S_{12} + 0.5S_{44}}
$$

Which wouldn't make sense for $\theta = 0^{\circ}, 90^{\circ}$.

### Components of the Compliance Tensor
Given the stiffness tensor for a 2D lattice with square symmetry, the nontrivial components are $C_{11}$, $C_{12}$, $C_{21}$, $C_{22}$, and $C_{44}$ (also known as $C_{66}$). The stiffness tensor in matrix form can be represented as:

$$
\mathbf{C} = \begin{pmatrix}
C_{11} & C_{12} & 0 \\
C_{21} & C_{22} & 0 \\
0 & 0 & C_{44}
\end{pmatrix}
$$

For square symmetry, typically $C_{11} = C_{22}$ and $C_{12} = C_{21}$. We can represent this tensor in a 2D matrix form as $\mathbf{M}$:

$$
\mathbf{M} = \begin{pmatrix}
C_{11} & C_{12} \\
C_{21} & C_{22}
\end{pmatrix}
$$

The compliance tensor $\mathbf{S}$ is the inverse of the stiffness tensor $\mathbf{C}$. We first calculate the inverse of matrix $\mathbf{M}$:

$$
\mathbf{M}^{-1} = \frac{1}{\text{det}(\mathbf{M})} \begin{pmatrix}
C_{22} & -C_{12} \\
-C_{21} & C_{11}
\end{pmatrix}
$$

where $\text{det}(\mathbf{M}) = C_{11}C_{22} - C_{12}C_{21}$. Assuming $C_{11} = C_{22}$ and $C_{12} = C_{21}$, this simplifies to:

$$
\text{det}(\mathbf{M}) = C_{11}^2 - C_{12}^2
$$

and

$$
\mathbf{M}^{-1} = \frac{1}{C_{11}^2 - C_{12}^2} \begin{pmatrix}
C_{11} & -C_{12} \\
-C_{12} & C_{11}
\end{pmatrix}
$$

The compliance tensor in matrix form is then:

$$
\mathbf{S} = \begin{pmatrix}
S_{11} & S_{12} & 0 \\
S_{21} & S_{22} & 0 \\
0 & 0 & S_{44}
\end{pmatrix}
$$

where $S_{11} = S_{22} = \frac{C_{11}}{C_{11}^2 - C_{12}^2}$, $S_{12} = S_{21} = \frac{-C_{12}}{C_{11}^2 - C_{12}^2}$, and $S_{44} = \frac{1}{C_{44}}$.

### Computing $\nu_{(01)(10)}$

$$
S_{11} = S_{22} = \frac{C_{11}}{C_{11}^2 - C_{12}^2}
$$



First: Look at visualiser to check what's going on with crystal

## Meeting: Weeks of the 8th of January
**Progress Updates:**
1. Graph for three denisities:


Output

Extrapolated Poisson's Constant at Delta = 0 for ρ = 0.90: -0.12844714336404295
Extrapolated Poisson's Constant at Delta = 0 for ρ = 0.87: -0.10831773702848677
Extrapolated Poisson's Constant at Delta = 0 for ρ = 0.85: -0.0905807683108304

**Note:** 0.87 needs to be run again for the altered seed
2. Thesis feedback recieved. Will continue to edit and also begin results section
3. Cluster run for larger N sent

**Questions**
1. Thoughts on elasticity graph + worth adding error bars?
2. For $\rho = 0.85$ fluctuations more prominent. Is this because we are near the phase tranisiton density? (in Anderson, 2017 it is $0.84$ for squares), or due to increased noise at lower $\rho$. A combination of both?
3. Larger simulation will be for $N = 4096$ (box length doubles). Does simulation time vary linearly with real time at fixed density for higher $N$?
4. Poisson constants vs Effective Poisson constants
5. For graphs in result sections. Any formatting tips welcome
6. For bulk modulus derivation in main thesis, would it be advise writing it in an appendix?

## Meeting: Weeks of the 16th of January [last meeting :(]
**Progress Updates:**
1. Results section written (EDMD)
**Things to keep in mind while reading:** Some of the graphs aren't as good image resolution as I thought I will be reuploading better quality ones
Monte Carlo section is still left but most analysis done but I wasn't planning on it being extensive
**What I'd like feedback on most:** is the interpretation and phenomenological statements made
2. Analysis on $N=4096$ almost done (just got last set of results results from cluster)
**Questions**
1. Literature comparison for MC
2. Analysis for MC worth completing?
3. Poisson constants vs Effective Poisson constants
**What's left?**
1. Introduction (to be done at the end, so Wednesday)
2. Larger N and more densities (as many as possible prior to thursday)
3. Theory and Methods draft 2 (Tuesday night)
4. Abstract, layman's summary etc(done with introduction)

## Appendix

### Determining Bulk modulus for 2D square


\begin{table}[h]
\centering
\begin{tabular}{|c|c|}
\hline
\textbf{Polynomial Degree} & \textbf{Bulk Modulus} \\
\hline
2 & 26.46 \\
\hline
3 & 100.14 \\
\hline
4 & 196.07 \\
\hline
5 & 143.98 \\
\hline
6 & 150.42 \\
\hline
\end{tabular}
\caption{Bulk Moduli for Different Polynomial Degrees at \( \rho = 0.87 \)}
\label{tab:bulk_moduli}
\end{table}

### MC Simulation $B_{11}$ and $B_{12}$ with different bulk modulus

Value of $B_{11}$ and $B_{12}$ a bit different than expected, is this normal? Using a different bulk modulus $B = 30.24859883753434$:


### Calculating $\mu$ for planar hard squares:

# Define the symbols explicitly
B_11, B_12, mu = sympy.symbols('B_{11} B_{12} mu')

# Given Poisson's ratio in terms of the stiffness tensor components for an anisotropic system
nu_anisotropic = -B_12 / B_11

# Expression for mu in terms of B and nu, before substitution
mu_expression = (B * (1 - nu)) / (nu + 1)

# Substitute B with its expression in terms of B_11 and B_12
B_expression = (B_11 + B_12)/2

# Substitute nu with the anisotropic expression for Poisson's ratio
mu_expression_substituted = mu_expression.subs({B: B_expression, nu: nu_anisotropic})

# Simplify the final expression for mu
mu_simplified = sympy.simplify(mu_expression_substituted)

# Display the simplified mu in explicit form
sympy.expand(mu_simplified)

$$
\
\nu_{(01), (10)} = \frac{-S_{12}}{S_{11}} \approx -9.978
\
$$

$$
\nu_{(01), (10)} = \frac{-S_{12}}{S_{11}} = \frac{B_{12}}{B_{11}}
$$
This value is outside the typical range for isotropic materials and indicates either anisotropic behavior or a non-traditional material, as it does not align with the expected range for standard linear elastic theory. The extreme value might indicate an error with the analysis. The corresponding graph is:


Now using `np.polyfit`and rechecking the fits by using varying polynomial degrees and calculating the bulk moduli corresponding to each fit for the equation of state for $\rho = 0.87$:

| Polynomial Degree | Bulk Modulus |
| ----------------- | ------------------- |
| 2 | 26.46395221745479 |
| 3 | 100.14107312599663 |
| 4 | 196.06596906669438 |
| 5 | 143.9782977104187 |
| 6 | 150.41816979646683 |

For the case of degree three:


**Determinant of Matrix $\mathbf{M}$:**
$$
\
\text{det}(\mathbf{M}) = C_{11}^2 - C_{12}^2 = (−30.72038122130251)^2 - (306.09285240964994)^2
\]
\[
\text{det}(\mathbf{M}) = 943.8225016327 - 93703.66476238945 = -92759.84226075675
\
$$
**Calculating $S_{11}$ and $S_{12}$**

$$
\
S_{11} = \frac{C_{11}}{\text{det}(\mathbf{M})} = \frac{-30.72038122130251}{-92759.84226075675}
\]
\[
S_{12} = \frac{-C_{12}}{\text{det}(\mathbf{M})} = \frac{-306.09285240964994}{-92759.84226075675}
$$

**Calculating the Poisson's constant $\nu_{(01), (10)}$:**