# Vera's project
###### tags: `bachelors projects`

## MC simulation
After looking at previous results from a year ago, we noticed that some of the results were faulty (packing fraction > 1), so I needed to check whether I had already fixed this or not.

I ran my current code with the same initialization variables as one of the faulty simulations. And the packing fraction was now more in the expected range, namely for $\beta P \sigma^3 = 10$, I ended with a packing fraction of 0.52. While this is higher than the Carnahan-Starling equation of state predicts at this pressure, it is at least a lot more realistic than the previous results.

I assume that I changed the code before I ran it again, and thus that the results I got after these faulty results are usable.

### Analyzing data
I could not find the code I used last year to analyze the data from my MC simulaions. Nor was I able to find my notes on how to analyze said data.

I'm now again studying how this works.

The idea is that we can calculate the elastic constants by looking at the fluctuations in strain.

For this, we use the Parrinello-Rahman fluctuation formula:
\begin{align}
B_{ijkl} = \frac{1}{\beta \langle V \rangle} \langle \epsilon_{ij} \epsilon_{kl} \rangle^{-1}
\end{align}
where $\mathbf{B}$ is the effective elastic constant, $\beta$ the thermodynamical beta, $\langle V \rangle$ is the average volume of the system. We also take the average over some components of the strain $\epsilon$.

We can determine these components of the strain using the following formula from Gusev and Zehnder:
\begin{align}
\epsilon_{ij} = \frac{1}{2} \big( h_{nk} \langle h \rangle^{-1}_{kj} h_{np} \langle h \rangle^{-1}_{pi} - \delta_{ij} \big)
\end{align}

here $h$ is the scaling matrix of our system. The scaling matrix consists of three vectors that form the computational box in which our particles exist. $\langle h \rangle$ is the reference shape of our box, i.e. the undeformed box (*this is in my writing from a year ago, but I didn't give a source...*)

We want to determine $B_{11}$, $B_{12}$ and $B_{44}$, or, in original notation: $B_{1111}, B_{1122}$ and $B_{1212}$.

To do this, we need to calculate $\epsilon_{11}$, $\epsilon_{22}$ and $\epsilon_{12}$.
Note that we cannot just calculate them out of context, because in the formula for $\epsilon_{ij}$, we also use the indices $k, l, n, p$. While $n, p$ indicate the sum $\Sigma_{n, p = 1}^{3}$, $k$ and $l$ come from $B_{ijkl}$.

I assume that $\langle h \rangle$ is the scaling matrix of the undeformed box, thus:
\begin{align}
h = \begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix} * C
\end{align}
where $C$ is some constant (I think in my case C = 7.5, but I wasn't able find a confirmation of that in my code),

This leads me to a problem, namely that $\epsilon_{22}$ always returns 0, because there's always $\langle h \rangle^{-1}_{kj} = \langle h \rangle^{-1}_{12} = 0$.
We know that this is not true, so I must be doing something wrong, but I don't know what.

## Calculating elastic constants (EDMD)
From the EDMD simulation, we get 2 data files as result, namely the positions of particles throughout the simulation, and the stress tensor throughout the simulation.

With this simulation are three different deformations we can apply to the box in which the particles are placed, called ShiftX, ShiftY and StrechZ. These deformations are given by the following deformation matrices:
\begin{align}
\mathbf{D}_\text{ShiftX} &= \begin{pmatrix} 1 & 0 & -\delta \\ 0 & 1 & 0 \\ 0 & 0 & 1
\end{pmatrix}, &
\mathbf{D}_\text{ShiftY} &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & -\delta \\ 0 & 0 & 1
\end{pmatrix}, \\
\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}

### Determining $B_{44}$
For determining $B_{44}$, where $B$ is the effective elastic constant, we use ShiftX. Using the function for the strain, $\epsilon = \frac{1}{2} \big(\mathbf{D}^T\mathbf{D} - \mathbf{1}\big)$, we know that the strain that corresponds with this deformation is:
\begin{align}
\epsilon = \begin{pmatrix}0 & 0 & -\frac{\delta}{2} \\ 0 & 0 & 0 \\ -\frac{\delta}{2} & 0 & \frac{\delta^2}{2}
\end{pmatrix}
\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 $\sigma_{ij}$ is an element of the stress tensor of the system before and $\sigma_{ij}'$ after the strain $\epsilon$ was applied.
We are interested in $B_{44} = B_{1313} = B_{1331}$. If we take $ij = 13$, then only for $kl = 13, 31, 33$ the second term on the left-hand side will be non-zero.
And 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}

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 $\delta$, the ShiftX we chose.

#### Trying different values for ShiftX
For 10 short simulations with different values of ShiftX in the range of $[0.0001, 0.001]$, we get the following plot:

I decided to run 5 more simulations, to see if there was a more clear choice for ShiftX going forward. This resulted in the following plot:

### Determening $B_{11}$ and $B_{12}$
We use a StrechZ deformation to determine $B_{11}$ and $B_{12}$:
\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 relation, the following is derived
\begin{align}
B_{ijkl} = \frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\epsilon_{kl}}
\end{align}

If we combine this with the equation of the strain, we get:
\begin{align}
\frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\delta} &= \frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\epsilon} \frac{\text{d}\epsilon}{\text{d}\delta}\\
&= \frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\epsilon_{11}} \frac{\text{d}\epsilon_{11}}{\text{d} \delta} + \frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\epsilon_{22}} \frac{\text{d}\epsilon_{22}}{\text{d}\delta} + \frac{\text{d}(\sigma_{ij}' - \sigma_{ij})}{\text{d}\mathcal{O}(\delta^2)} \frac{\text{d}\mathcal{O}(\delta^2)}{\text{d} \delta}\\
&= B_{ij11} - B_{ij22} + \mathcal{O}(\delta)
\end{align}

We know that before any deformation was applied, the only stress on the system is pressure, $\sigma_{ij} = -P \delta_{ij}$. As we are interested in $B_{11} = B_{1111}$ and $B_{12} = B_{1122}$ specifically, we take $ij = 11$. This gives us
\begin{align}
\frac{\text{d}(\sigma_{11}' - P)}{\text{d}\delta} = B_{11} - B_{12} + \mathcal{O}(\delta)
\end{align}

We can get $P$ from the stress using:
\begin{align}
P = \frac{\text{tr}(\sigma')}{3}
\end{align}

However, the equation above does not give us $B_{11}$ and $B_{12}$ seperately, which is what we want.
For that we need the bulk modulus $B$:
\begin{align}
B = \frac{1}{3}\big(B_{11} + 2B_{12}\big) = - V \frac{\text{d}P}{\text{d}V}
\end{align}

We can find the value of the bulk modulus by determining the equation of state of the system. For hard spheres we use the Speedy equation of state.
For hard cubes we will determine the equation of states ourselves.
But for both we will need extra simulations at each density with no deformation.

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}

#### Trying different values for StretchZ
For 10 short simulations with different values of StretchZ in the range of $[0.0001, 0.001]$, we get the following plot:

Here I also decided to run 5 more simulations, from which I got the following plot:

This evens out more nicely than the plot for ShiftX

#### Comparing results
For hard spheres, the elastic constants have already been calculated by D. Frenkel and A. Ladd ('Elastic constants of hard-sphere crystals', 1987).

We can compare the results we got to their results. However, there are a few things to keep in mind.

First of all, in some papers they talk about 'elastic constants', but are actually talking about 'effective elastic constants'. I can't remember whether this was also the case with this paper.

Second, I need to make sure that we're using the same units. I'm mostly confused by the density. I'm unsure what the unit is of the density in the simulation, and in the paper they are using 'reduced density' and I cannot figure out how it relates to 'our' density.

Despite this, I tried to just compare the results I had to the results of Frenkel and Ladd.
They have the following results:

Our results are at 'density 1.00'.

For $B_{11}$ we have:
which then for $C_{11}$ gives us:

For $B_{12}$ we have:
which then for $C_{12}$ gives us:

For $B_{44}$ we have:
which then for $C_{44}$ gives us:

\begin{align}
B_{ijkl} = C_{ijkl} - P(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl}) \label{eff.el.const.}
\end{align}

### Overview results

In the images below I have assembeled the results from the different simulations at different densities, with the results from Frenkel also plotted.

#### Results for $C_{11}$
Note that the results I get for $C_{11}$ are all much lower than the results Frenkel got.

#### Results for $C_{12}$
Note that the results I get for $C_{12}$ are all much higher than the results Frenkel got.

#### Results for $C_{44}$
The results I get for $C_{44} agree with the results Frenkel got.
At a higher density, $C_{44}$ gets bigger.
When my density is inbetween two densities of Frenkel, my value for $C_{44}$ is inbetween his values for $C_{44}$. When my density is slightly higher than his density, the value I find for $C_{44}$ is bigger than his value of $C_{44}$. So it seems that this works.

Now I need to figure out why it doens't work for $C_{11}$ and $C_{12}$

## Pictures



## To do
- [x] Check if results from MC simulation makes sense
- [x] Try to reproduce data point with EDMD
- [x] try out different shears to see which one works best
- [x] figure out again how to calculate the elastic constants
- [x] Check if results EDMD corresponds with D.Frenkel
- [ ] Figure out what goes wrong with $C_{11}$ and $C_{12}$
- [ ] Figure out what I'm doing wrong when I want to calculate the strain from a MC simulation
- [ ] Ask Frank what it means when I set the density to 1 (because spheres cannot fill an entire cube)