# Poisson ratios polyhedra

###### tags: `PhD projects`

---

::: danger
Check whether all elastic constant runs run long enough
Compute SD for all Poisson ratios
:::
# New EDMD code for point-assymetric particles
*Found on HERA /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D_cubes/ and /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D_triprism/*
Since Franks EDMD code was designed for point-assymetric particles, we wanted to alter the code to be also able to deal with point-assymetric particles.
### Initial shape
*See notebook shapesnonpoint.nb in the folder on HERA /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D_cubes/Shapes/*

Since we cannot rely on symmetry anymore, we have to save more edges/faces/vertices for each particle:
- We save all vertices
- We save all edges that have a different direction, Per edge we then link all edges that have the same direction but a different position in the particle (this is needed if we want to find the touching point of edge-edge collisions)
- We save all different facenormals.

For a triangular prism we plot on the left thed edges, vertices and facenormals associated with that particle (where here each color edges belongs to the same group of edges). Moreover, per facenormal we save the maximum and projection of a particle on that facenormal. On the right we plot for each facenormal the respective maximum (red) and minimum (blue) projection vectors.

| Edges, vertices and facenormals |Maximum and minimum projection |
| -------- | -------- |
|  |  |

### Seperating axis theorem
For the seperating axis theorem we use the following algorithm. Given two particles A and B, rotation matrices $\mathbf{R}_A$ and $\mathbf{B}_A$ and seperation distance $\mathbf{r}_{AB}$ we loop over

**All facesnormals *L* of particle *A***
- Compute $\mathbf{L}^{rot} = \mathbf{B}^T\mathbf{L}$ to find $\mathbf{L}$ in the body frame of B.
- Project all vertices on $\mathbf{L}^{rot}$ to find maximum and mimimum projection of particle B.
- Take the earlier computed maximum and mimumum projection of the vertices of particle A on $\mathbf{L}$.
- If $d\mathbf{r}\cdot \mathbf{L} > 0$, the seperation projected on that axis is given by $\text{sep} = d\mathbf{r}\cdot \mathbf{L} - \text{maxa} + \text{minb}$;
- If $d\mathbf{r}\cdot \mathbf{L} < 0$, the seperation projected on that axis is given by $\text{sep} = - d\mathbf{r}\cdot \mathbf{L} + \text{mina} - \text{maxb}$

**All facesnormals *L* of particle *B***
- Compute $\mathbf{L}^{rot} = \mathbf{A}^T\mathbf{L}$ to find $\mathbf{L}$ in the body frame of A.
- Project all vertices on $\mathbf{L}^{rot}$ to find maximum and mimimum projection of particle A.
- Take the earlier computed maximum and mimumum projection of the vertices of particle B on $\mathbf{L}$.
- If $d\mathbf{r}\cdot \mathbf{L} > 0$, the seperation projected on that axis is given by $\text{sep} = d\mathbf{r}\cdot \mathbf{L} - \text{maxa} + \text{minb}$;
- If $d\mathbf{r}\cdot \mathbf{L} < 0$, the seperation projected on that axis is given by $\text{sep} = - d\mathbf{r}\cdot \mathbf{L} + \text{mina} - \text{maxb}$

**All cross-products between edges**
- Compute $\mathbf{L}^{rot}_A = \mathbf{A}^T\mathbf{L}$ to find $\mathbf{L}$ in the body frame of A.
- Compute $\mathbf{L}^{rot}_B = \mathbf{B}^T\mathbf{L}$ to find $\mathbf{L}$ in the body frame of B.
- Project all vertices on $\mathbf{L}^{rot}_A$ to find maximum and mimimum projection of particle A.
- Project all vertices on $\mathbf{L}^{rot}_B$ to find maximum and mimimum projection of particle B.
- Take the earlier computed maximum and mimumum projection of the vertices of particle B on $\mathbf{L}$.
- If $d\mathbf{r}\cdot \mathbf{L} > 0$, the seperation projected on that axis is given by $\text{sep} = d\mathbf{r}\cdot \mathbf{L} - \text{maxa} + \text{minb}$;
- If $d\mathbf{r}\cdot \mathbf{L} < 0$, the seperation projected on that axis is given by $\text{sep} = - d\mathbf{r}\cdot \mathbf{L} + \text{mina} - \text{maxb}$

In principle this is all we need. Throughout the entire code this new algorithm was implemented.

### Check EOS
We checked the equation of state for both hard cubes and hard triangular prisms.
#### Cubes
*See EOScubes.nb in the folder /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D/CHECKS_cubes/, the data on the cubes obtained via Franks old code can be found in the folder /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D/CHECKS_oldcubes/*

We compare the data from Franks old code, the new code and the data read from Figure 10 B in the thesis of Nena slaats.

| EOS cubes |EOS triangular prisms |
| -------- | -------- |
| |  |

#### Triangular prisms
*See EOStriangularprisms.nb in the folder /home/users/5670772/Elastic_constants/CODE_Polyhedra_EDMD_3D/CHECKS_triprisms/*

We compare the data from Marjoleins MC code, the new code and the data read from \url{https://www.nature.com/articles/nmat2959}.

Note that pressure is expressed in units of $\sigma$, where $\sigma$ is the edglength of the triangle. As mentioned before we use triangles with a ratio of $a/b = 1$, where a is the height of the prism and b is equal to two times the inradius of the polygonal face (so in this case the triangular face). Particles with an edge length $\sigma =1$ have an innerscribed circle of radius $\frac{1}{2\sqrt{3}}$. As a result, to obtain $a/b = 1$, the height of the triangular prism should be $h =\frac{1}{2\sqrt{3}}\cdot2 = \frac{1}{\sqrt{3}}$. The volume of a triangular prism in general is given by $\frac{\sqrt{3}}{4}\cdot\sigma^2\cdot h$, which is $1/4$ in the case of an prism with $\sigma=1$.

Marjolein uses particles with $\sigma=1$, but since we use particles with volume 1, the edgelength is given by $\sigma=4^{1/3} = 1.5874$ (and trivially a height of $h = 4^{1/3}\cdot\frac{1}{\sqrt{3}} = 0.9164$. Our measured pressure thus needs to be multiplied by $1.5874^3$.

The paper uses triangular prisms with $\sigma = 1$, but expresses the pressure in terms of the outerscribed circle radius, which is given by $\frac{1}{3} \sqrt{\frac{5}{3}}$. We thus need to divide their pressure by $\left[\frac{1}{3}\sqrt{\frac{5}{3}}\right]^3$.

#### Hexagonal prisms

For the triangular prisms we again use a shape that has $a/b = 1$, ie. the height is given by $h = 1.049115$ and the sides of the hexagon $\sigma = 0.605707$, therefore, the innerscribed circle has a radius of $r_\text{in} = \frac{1}{2}\sqrt{3}\sigma = 0.524558 = \frac{b}{2}$, therefore $b =1.04912 = h$
# Theory on strain and stress

Crystals undergoing a small deformation will induce an elastic response. The deformation is generally expressed in the so called strain tensor, while the response is measured in terms of the stress tensor. Both tensors will be discussed below. The proportionality constants between the two tensors are the elastic constants. Below we will discuss how these elastic constants can be measured. For that we make the assumption that, when deforming the box, all changes in free energy are due to changed distances between particles.

## Deformation

Given a solid in an undeformed box under hydrostatic (isotropic) pressure. If we deform the box using a deformation tensor $\mathbf{D}$ the distances between particles in the crystal will change. An arbitrary vector $\mathbf{r}$ will change to $\mathbf{r}'=\mathbf{D}\mathbf{r}$, with $\mathbf{r}'$ the vector in the deformed system. We can express $\mathbf{D}$ in terms of the displacement vector $\mathbf{u} = \mathbf{r}'-\mathbf{r}$. To do so we first calculate $r'^2$ in terms of $\mathbf{u}$ and $\mathbf{r}$ (using Einstein notation).
$$
\begin{align}
r'^2 &= (r_i- u_i)^2\\
&= (r_i- \frac{\partial u_i}{\partial r_k} r_k)(r_i- \frac{\partial u_i}{\partial r_l} r_l)\\
&=r^2 + 2 \frac{\partial u_i}{\partial r_k}r_k r_i + \frac{\partial u_i}{\partial r_k}\frac{\partial u_i}{\partial r_l}r_k r_l\\
&= r^2 + \left[\frac{\partial u_i}{\partial r_k} + \frac{\partial u_k}{\partial r_i} + \frac{\partial u_l}{\partial r_k} \frac{\partial u_l}{\partial r_i}\right]r_ir_k
\end{align}
$$
Where for this last term we re-assigned some indices. Given this last expression, we defined the strain tensor as
$$
\epsilon_{ij} = \frac{1}{2}\left[\frac{\partial u_i}{\partial r_j} + \frac{\partial u_j}{\partial r_i} + \frac{\partial u_k}{\partial r_i} \frac{\partial u_k}{\partial r_j}\right]
$$
Such that $r'^2 = r^2 + 2\epsilon_{ij}r_ir_j$.

::: info
Note that the strain tensor only holds deformations and not rotations. Therefore there are only 6 components (push and shift in every direction). Note that there is never twist, we only consider uniform deformations.
:::

Given this definition, we can rewrite $\mathbf{D}$ in terms of $\epsilon$. Since $r'_i = D_{ij}r_j$, $u_i = r'_i-r_i =(D_{ij}-\delta_{ij})r_j$ and thus $\frac{\partial u_i}{\partial r_j} =D_{ij}-\delta_{ij}$. Therefore
$$
\begin{align}
\epsilon_{ij} &= \frac{1}{2}\left[D_{ij} - \delta_{ij} + D_{ji}- \delta_{ji}+(D_{ki}-\delta_{ki})(D_{kj}-\delta_{kj})\right]\\
&=\frac{1}{2}\left[D_{ki}D_{kj} - \delta_{ij}\right]\\
&=\frac{1}{2}\left[\mathbf{D}^T\mathbf{D} - \mathbf{I}\right]_{ij}
\end{align}
$$

Lastly we can also express the deformed volume of the system in terms of $\mathbf{\epsilon}$, namely $V' = V\det(\mathbf{I}+\mathbf{\epsilon})$. Since for small $\epsilon$ we can use that $\det(\mathbf{I}+\mathbf{\epsilon})= 1+\text{Tr}(\epsilon)$, we can write

$$
V' = V(1+\text{Tr}(\epsilon))
$$

Defining $\mathbf{G} = \mathbf{D}^T\mathbf{D}$, we can see that we can write
$$\mathbf{G} = \mathbf{I} + 2\mathbf{\epsilon},
\tag{1}$$
which is someting that we will use later.

## Free energy
Given these tools above we can look at the effect of the deformations. To determine how much work it costs to deform a certain system, we want to know the change in free energy. Since in general the deformations that we consider are small, we can Taylor the free energy in terms of $\epsilon$

$$
\begin{align}
\frac{F(\epsilon)}{V}&= \frac{F(0)}{V} + \frac{1}{V}\frac{\partial F(\epsilon)}{\partial\epsilon_{ij}}\bigg\rvert_{\epsilon_{ij} = 0}\epsilon_{ij}+ \frac{1}{2V}\frac{\partial^2 F(\epsilon)}{\partial\epsilon_{ij}\partial\epsilon_{kl}}\bigg\rvert_{\epsilon_{ij},\epsilon_{kl} = 0}\epsilon_{ij}\epsilon_{kl}\\
&=\frac{F(0)}{V} + C_{ij}\epsilon_{ij}+ \frac{1}{2}C_{ijkl}\epsilon_{ij}\epsilon_{kl}
\end{align}
$$

With $C_{ij}$ and $C_{ijkl}$ respectively the first and second order elastic constants. In order to measure these we need to express the free energy in terms of $\epsilon$. For this we use that $F(\epsilon) = -k_BT\log Z(\epsilon)$, with
$$
Z(\epsilon) = \frac{1}{h^3N!}\int d\mathbf{r}'^Nd\mathbf{p}'^Ne^{-\beta \mathcal{H}(\mathbf{r}'^N,\mathbf{p}'^N)}
$$

Since $\frac{\partial F(\epsilon)}{\partial \epsilon} = -k_BT\frac{1}{Z(\epsilon)}\frac{\partial Z(\epsilon)}{\partial \epsilon}$, we need to take the derivative of $Z(\epsilon)$ w.r.t $\epsilon$. For that we first write $\mathbf{r}'$ and $\mathbf{p}'$ in terms of $\epsilon$. Trivially, $\mathbf{r}' = \mathbf{D}\mathbf{r}$.

To find $\mathbf{p}'$, we use that $\mathbf{p}'_i = \frac{\partial \mathcal{K}}{\partial \mathbf{r}_i'}$. Since $\mathcal{K} = \sum_i \frac{1}{2}m_i \mathbf{\dot{r'}}_i^2 = \sum_i \frac{1}{2}m_i \mathbf{\dot{r}}_i\mathbf{D}^T\mathbf{D}\mathbf{\dot{r}}_i$, we find that
$$
p^\alpha_i = m_i [\mathbf{D}^T\mathbf{D}]^{\alpha\beta}\dot{r}^{\beta} \text{ and }\dot{r}_i^{\alpha} = \frac{1}{m_i}\left([\mathbf{D}^T\mathbf{D}]^{-1}\right)^{\alpha\beta}{p}_i^{\beta}
$$
And since
$$
\mathbf{p}'_i= m_i\mathbf{\dot{r'}}_i = m_i\mathbf{D}\cdot\mathbf{\dot{r}}_i = m_i\mathbf{D}\frac{1}{m_i}\frac{1}{\mathbf{D}^T\mathbf{D}}{p}_i=\frac{1}{\mathbf{D}^T}\mathbf{p}_i
$$
This also means that we can write $\mathcal{K} = \sum_i \frac{1}{2}m_i \mathbf{\dot{r}}_i\mathbf{D}^T\mathbf{D}\mathbf{\dot{r}}_i = \sum_i \frac{1}{2m_i}\mathbf{\dot{p}}_i[\mathbf{D}^T\mathbf{D}]^{-1}\mathbf{\dot{p}}_i$
Since $\mathbf{p}'_i = \frac{1}{\mathbf{D}^T}\mathbf{p}_i$ and $\mathbf{r}'_i = \mathbf{D}\mathbf{r}_i$ the jacobian of the transformation between $d\mathbf{r}^Nd\mathbf{p}$ and $d\mathbf{r}^Nd\mathbf{p}$ is 1. Hence,
$$
\begin{align}
Z(\epsilon) &= \frac{1}{h^3N!}\int d\mathbf{r}'^Nd\mathbf{p}'^Ne^{-\beta \mathcal{H}(\mathbf{r}'^N,\mathbf{p}'^N)}\\
&=\frac{1}{h^3N!}\int d\mathbf{r}^Nd\mathbf{p}^Ne^{-\beta \left[\sum_i \frac{1}{2m_i} \mathbf{{p}}_i[\mathbf{D}^T\mathbf{D}]^{-1}\mathbf{{p}}_i + U(\mathbf{r}^N, \epsilon)\right]}
\end{align}
$$

:::info
Note that if we assume pressure is hydrostatic, the momentum transfer cannot depend on the coordinate system. In that case we can integrate out the momenta, leaving only the integral over $\mathbf{r}$. In that casethe jacobian however is no longer 1, resulting in extra factors of $\det(\mathbf{D})$, which will lead to a volume factor (see Vera's thesis). In the end this will lead to the same result, however, because of the hydrostatic assumption, the stress tensor will be less general.
:::

Now we can take derivatives
$$
\begin{align}
\frac{\partial U}{\partial \epsilon_{\alpha\beta}}& = \sum_{i<j}\frac{\partial U}{\partial r'^2_{ij}}\frac{\partial r'^2_{ij}}{\partial \epsilon_{\alpha\beta}}\\
& = \sum_{i<j}\frac{\partial U}{\partial r'_{ij}}\frac{\partial r'_{ij}}{\partial r'^2_{ij}}\frac{\partial r'^2_{ij}}{\partial \epsilon_{\alpha\beta}}
\end{align}
$$

Now we can use that $\frac{\partial U}{\partial r^2_{ij}}= f(r_{ij})$, $\frac{\partial r'^2_{ij}}{\partial \epsilon_{\alpha\beta}} = \frac{\partial (r_{ij}^2 + 2\epsilon_{\gamma\delta}r_{ij, \gamma}r_{ij, \delta})}{\epsilon_{\alpha\beta}} = 2r_{ij, \alpha}r_{ij, \beta}$, such that

$$
\begin{align}
\frac{\partial U}{\partial \epsilon_{\alpha\beta}}& = \sum_{i<j}\frac{\partial U}{\partial r'^2_{ij}}\frac{\partial r'^2_{ij}}{\partial \epsilon_{\alpha\beta}}\\
& = \sum_{i<j}f(r'_{ij})\frac{1}{ 2r'_{ij,\gamma}}2r_{ij, \alpha}r_{ij, \beta}\\
& = \sum_{i<j}f(r'_{ij})\frac{1}{ 2r'_{ij}}2D^{-1}_{\alpha,\gamma}r'_{ij, \gamma}r_{ij, \delta}[D^{-1}]^T_{\delta,\beta}\\
& =\left[\mathbf{D}^{-1} \sum_{i<j}\mathbf{f}(r'_{ij})\otimes \mathbf{r}_{ij}[\mathbf{D}^{-1}]^T\right]_{\alpha\beta}
\end{align}
$$

where we used that $f(r'_{ij})\frac{1}{ r'_{ij}}r'_{ij, \gamma} = f_\gamma(r'_{ij})$, and where $[\mathbf{f}(r'_{ij})\otimes \mathbf{r}_{ij}]_{\alpha,\beta} = f_\alpha(r'_{ij}) r_{ij, \beta}$

For the second derivative we need to look at
$$
\begin{align}
\frac{\partial (\sum_i \frac{1}{2m_i} \mathbf{{p}}_i[\mathbf{D}^T\mathbf{D}]^{-1}\mathbf{{p}}_i)}{\partial \epsilon_{\alpha\beta}}& = \sum_i \frac{1}{2m_i} \mathbf{{p}}_i\frac{\partial \mathbf{G}^{-1}}{\partial \epsilon_{\alpha\beta}}\mathbf{{p}}_i,
\end{align}
$$
where we used the definition of $\mathbf{G}$ in Eq. (1). Since $\mathbf{G}^{-1} = \frac{1}{\mathbf{I}+ 2\mathbf{\epsilon}}$, we find that
$$
\begin{align}
\sum_i \frac{1}{2m_i} \mathbf{{p}}_i\frac{\partial [\mathbf{D}^T\mathbf{D}]^{-1}}{\partial \epsilon_{\alpha\beta}}\mathbf{{p}}_i &= \sum_i \frac{1}{2m_i} {{p}}_{i,\gamma}\frac{\partial [{G}]_{\gamma,\delta}^{-1}}{\partial \epsilon_{\alpha\beta}}{{p}}_{i,\delta}\\
&= \sum_i \frac{1}{2m_i} {{p}}_{i,\gamma}\frac{-2}{G^2_{\gamma\delta}}\delta_{\alpha\gamma}\delta_{\beta\delta}{{p}}_{i,\delta}\\
&= -\sum_i \frac{1}{m_i} \left[\mathbf{p}_{i}\mathbf{G}^{-1}\right]_{\alpha}\left[\mathbf{G}^{-1}\mathbf{{p}}_i\right]_{\beta}\\
&=-\left([\mathbf{D}]^{-1} \sum_i \frac{1}{m_i} \left[\mathbf{p}'_i\otimes\mathbf{p}'_i\right][\mathbf{D}^T]^{-1}\right)_{\alpha\beta},
\end{align}
$$
where for the last part we used that $\mathbf{p}'_i= \frac{1}{\mathbf{D}^T}\mathbf{p}_i$ and that $\mathbf{G} = \mathbf{D}^T\mathbf{D}$

Both derivations above lead to the expression
$$
\begin{align}
\frac{1}{V}\frac{\partial F(\epsilon)}{\partial\epsilon_{ij}} &= -\frac{k_BT}{V}\frac{\partial\log Z(\epsilon)}{\partial\epsilon_{ij}}\\
&=\frac{k_BT}{V}\frac{1}{Z}\frac{-\beta}{h^3N!}\int d\mathbf{r}^Nd\mathbf{p}^Ne^{-\beta \left[\sum_i \frac{1}{2m_i}\mathbf{{p}}_i[\mathbf{D}^T\mathbf{D}]^{-1}\mathbf{{p}}_i + U(\mathbf{r}^N, \epsilon)\right]}\left[[\mathbf{D}]^{-1} \sum_i\left(\frac{1}{m_i} \left[\mathbf{p}'_i\otimes\mathbf{p}'_i\right]+\sum_{i<j}\mathbf{f}(r'_{ij})\otimes \mathbf{r}_{ij}\right)[\mathbf{D}^{-1}]^T\right]_{\alpha\beta}\\
&=\frac{-1}{V}\bigg\langle[\mathbf{D}]^{-1} \sum_i\left(\frac{1}{m_i} \left[\mathbf{p}'_i\otimes\mathbf{p}'_i\right]+\sum_{i<j}\mathbf{f}(r'_{ij})\otimes \mathbf{r}_{ij}\right)[\mathbf{D}^{-1}]^T\bigg\rangle_{\alpha\beta}\\
&=\frac{V'}{V}\bigg\langle[\mathbf{D}]^{-1} \sigma[\mathbf{D}^{-1}]^T\bigg\rangle_{\alpha\beta}\\
\end{align}
$$
With the stress tensor defined as
$$
\sigma_{\gamma,\delta} = -\frac{1}{V'}\sum_i\left[\frac{p_{i\gamma}p_{i\delta}}{m_i}+\sum_{j<i}r_{ij\gamma}f_{ij\delta}\right]
$$
Note that this is exactly minus the pressure tensor. For an undeformed system under hydrostratic pressure, the off diagonal terms of the pressure tensor fall away, and we end up with $\frac{1}{V'}\sum_i\frac{p_{i\gamma}p_{i\delta}}{m_i} = \rho k_BT \mathbf{I}$. In general we then find $\sigma_{\alpha\beta} = -P\delta_{\alpha\beta}$.

::: info
Note that $\sigma$ is the linear response to deformation (i.e. pressure), while we are interested in the quadratic response.l
:::

Filling our derivative in in the elastic constant equation we find that
$$
C_{\alpha\beta} = \frac{1}{V}\frac{\partial F(\epsilon)}{\partial \epsilon_{\alpha\beta}} = \frac{V'}{V}\left[[\mathbf{D}]^{-1}\sigma[\mathbf{D}^T]^{-1}\right]_{\alpha\beta}
$$

Now we can use that $V' = V\det(1+\epsilon ) = V(1+\text{Tr}(\epsilon) +\mathcal{O}(\epsilon^2))$, and that for very small deformations, we can ignore $\mathcal{O}(\mathbf{u}^2)$. As a result, the strain tensor simplifies to $\epsilon_{ij} = \frac{1}{2}\left[\frac{\partial u_i}{\partial r_j} + \frac{\partial u_j}{\partial r_i} \right]$ and thus that $\mathbf{\epsilon} = \frac{1}{2}\left[\mathbf{D}- \mathbf{I}+\mathbf{D}^T-\mathbf{I}\right]$. Assuming $\mathbf{D}$ is symmetric (one can prove that any assymatric $\mathbf{\bar{D}}$ can always be written as a rotation matrrix times a symmetric matrix $\mathbf{D}$. Since the free energy is unaffected by rotation matrices, it will only depend on the symmetric matrix $\mathbf{D}$), we can thus write

$$
\begin{align}
\mathbf{\epsilon} = \mathbf{D}-\mathbf{I}\text{ & }\mathbf{D} = \mathbf{\epsilon}+\mathbf{I}.
\end{align}
$$

As a result, up to $\mathcal{O}(\mathbf{\epsilon}^2)$, we find that $\mathbf{D}^{-1} = \mathbf{I}-\mathbf{\epsilon}$ (such that $\mathbf{D}^{-1}\mathbf{D} = \mathbf{I}+\mathcal{O}(\mathbf{\epsilon}^2)$). This means that we can rewrite the first elaxtic constant.
$$
\begin{align}
C_{\alpha\beta} &= \frac{\partial F(\epsilon)}{\partial \epsilon_{\alpha\beta}} = (1+\text{Tr}(\epsilon))\left[(\delta_{\alpha\gamma}-\epsilon_{\alpha\gamma})\left[\sigma_{\gamma\delta}(\epsilon = 0)+\frac{\partial \sigma_{\gamma\delta}}{\partial \epsilon_{\nu\rho}}\bigg\vert_{\epsilon = 0}\epsilon_{\nu\rho}\right](\delta_{\delta\beta}-\epsilon_{\delta\beta})\right]+\mathcal{O}(\epsilon^2)\\
&= \left(\sigma_{\alpha\beta}(\epsilon = 0)+\sigma_{\alpha\beta}(\epsilon = 0)\text{Tr}(\epsilon)- \epsilon_{\alpha\nu}\sigma_{\nu\beta}+ \frac{\partial \sigma_{\alpha\beta}}{\partial \epsilon_{\nu\rho}}\bigg\vert_{\epsilon = 0}\epsilon_{\nu\rho}-\sigma_{\alpha\nu}\epsilon_{\beta\nu}\right)
\end{align}
$$

By now evaluation the second derivative of this tensor w.r.t $\epsilon_{\gamma\delta}$, we find

$$
\begin{align}
C_{\alpha\beta\gamma\delta} &= \frac{1}{V}\frac{\partial^2 F(\epsilon)}{\partial \epsilon_{\alpha\beta}\partial \epsilon_{\gamma\delta}}\bigg\vert_{\epsilon=0} = \frac{\partial}{\partial \epsilon_{\gamma\delta}}\left(\sigma_{\alpha\beta}(\epsilon = 0)+\sigma_{\alpha\beta}(\epsilon = 0)\text{Tr}(\epsilon)- \epsilon_{\alpha\nu}\sigma_{\nu\beta}+ \frac{\partial \sigma_{\alpha\beta}}{\partial \epsilon_{\nu\rho}}\bigg\vert_{\epsilon = 0}\epsilon_{\nu\rho}-\sigma_{\alpha\nu}\epsilon_{\beta\nu}\right)\\
&= \left(\delta_{\gamma\delta}\sigma_{\alpha\beta}- \delta_{\alpha\gamma}\sigma_{\delta\beta}+ \frac{\partial \sigma_{\alpha\beta}}{\partial \epsilon_{\gamma\delta}}\bigg\vert_{\epsilon = 0}-\delta_{\beta\gamma}\sigma_{\alpha\delta}\right)
\end{align}
$$

If the system initially is under a hydrostatic pressure $P$, such that $\sigma_{\alpha\beta} = -P\delta_{\alpha\beta}$, we find
$$
\frac{\partial \sigma_{\alpha\beta}}{\partial \epsilon_{\gamma\delta}}\bigg\vert_{\epsilon = 0} = C_{\alpha\beta\gamma\delta}+\left(\delta_{\alpha\beta}\delta_{\gamma\beta}+ \delta_{\alpha\gamma}\delta_{\delta\beta} -\delta_{\gamma\delta}\delta_{\alpha\beta}\right)P \equiv B_{\alpha\beta\gamma\delta},
$$
where $B_{\alpha\beta\gamma\delta}$ are the **effective elastic constants**. As we can see, these effective elastic constants are what connects the strain and stress tensor.

Since the strain tensor is symmetric and the order of derivatives does not matter, we see that (we switch to $ijkl$)
$$
\frac{\partial^2 F(\epsilon)}{\partial \epsilon_{ij}\partial \epsilon_{kl}} = \frac{\partial^2 F(\epsilon)}{\partial \epsilon_{ji}\partial \epsilon_{lk}}= \frac{\partial^2 F(\epsilon)}{\partial \epsilon_{lk}\partial \epsilon_{ji}}= \frac{\partial^2 F(\epsilon)}{\partial \epsilon_{kl}\partial \epsilon_{ij}}
$$
which means that
$$
C_{ijkl}=C_{jilk}=C_{lkji}=C_{klij}
$$
Since the first two components and the second two components are interchangable we can write $C$ as a $6\times 6$ matrix $C_{ij}$. Since for this new tensor it is true that $C_{ij} = C_{ji}$, there are $21= 1+2+3+4+5+6$ independent variables.

In the mathematica notebooks on HERA */home/staff/5670772/Elastic_constants/Analysis_symmetries/ (**N.b. Note that in these notebooks we use a different Voigt notation, where $xy=4$,$xz=5$ and $yz=6$. It does not change the analysis, but it does change the naming and the place of the entries in the matrix)***, one can find the analysis for the symmetries for all the four particle types we will look at. We find
#### Squares
$$
\mathbf{B} = \begin{pmatrix}
B_{11} & B_{12} & 0 &0 &0 &0\\
B_{12} & B_{11} & 0 &0 &0 &0\\
0 & 0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$
#### Triangles
$$
\mathbf{B} =\begin{pmatrix}
B_{11} & B_{12} & 0 &0 &0 &0\\
B_{12} & B_{11} & 0 &0 &0 &0\\
0 & 0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$
where $B_{66}= \frac{B_{11}-B_{12}}{2}$
#### Hexagons
$$
\mathbf{B} =\begin{pmatrix}
B_{11} & B_{12} & 0 &0 &0 &0\\
B_{12} & B_{11} & 0 &0 &0 &0\\
0 & 0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$
where $B_{66}= \frac{B_{11}-B_{12}}{2}$
#### Cubes
$$
\mathbf{B} =\begin{pmatrix}
B_{11} & B_{12} &B_{12} &0 &0 &0\\
B_{12} & B_{11} & B_{12} &0 &0 &0\\
B_{12} & B_{12} & B_{11} &0 &0 &0\\
0 &0 & 0 &B_{66} &0 &0\\
0 &0 & 0 &0 &B_{66} &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$

#### Triangular prisms

$$
\mathbf{B} =\begin{pmatrix}
B_{11} & B_{12} &B_{23} &0 &0 &0\\
B_{12} & B_{11} & B_{23} &0 &0 &0\\
B_{23} & B_{23} & B_{33} &0 &0 &0\\
0 &0 & 0 &B_{55} &0 &0\\
0 &0 & 0 &0 &B_{55} &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$
where $B_{66}= \frac{B_{11}-B_{12}}{2}$
#### Hexagonal prisms
$$
\mathbf{B} =\begin{pmatrix}
B_{11} & B_{12} &B_{23} &0 &0 &0\\
B_{12} & B_{11} & B_{23} &0 &0 &0\\
B_{23} & B_{23} & B_{33} &0 &0 &0\\
0 &0 & 0 &B_{55} &0 &0\\
0 &0 & 0 &0 &B_{55} &0\\
0 &0 & 0 &0 &0 &B_{66}\\
\end{pmatrix}
$$
where $B_{66}= \frac{B_{11}-B_{12}}{2}$
::: danger
The tensor is the same for triangles as for hexagons right?
:::
## Bulk modulus
The bulk modulus is a measure of how much a substance will compress under a given amount of external pressure:
$$
B = -V\frac{dP}{dV}
$$
We can however also express the bulk modulus in terms of the elastic energy density

$$
F = \frac{1}{2}B_{ij}\epsilon_{i}\epsilon_{j}
$$
Given a uniform dilation $\epsilon_{ii} = \frac{1}{3}\delta$, i.e. a deformation given by
$$\mathbf{D} =
\begin{pmatrix}
1+\frac{1}{3}\delta & 0 \\
0 &1+\frac{1}{3}\delta &0\\
0 & 0 & 1+\frac{1}{3}\delta\\
\end{pmatrix}
$$ and thus a strain tensor ($\mathbf{D}^T\mathbf{D}-\mathbf{1}$) that up to first order is given by
$$\mathbf{\epsilon} =
\begin{pmatrix}
\frac{1}{3}\delta & 0 \\
0 &\frac{1}{3}\delta &0\\
0 & 0 & \frac{1}{3}\delta\\
\end{pmatrix}
$$
the Bulk modulus is given by
$$
F = \frac{1}{2}B\delta^2
$$
Of course for 2D the deformation becomes
$$\mathbf{D} =
\begin{pmatrix}
1+\frac{1}{2}\delta &0 \\
0 &1+\frac{1}{2}\delta\\
\end{pmatrix}
$$
such that
$$\mathbf{\epsilon} =
\begin{pmatrix}
\frac{1}{2}\delta & 0\\
0 &\frac{1}{2}\delta \\
\end{pmatrix}
$$
Such that the change in volume is $\delta^2$
#### Squares
Given the $\mathbf{B}$ matrix for squares we find an energy density of
$$\begin{align}
F &= \frac{1}{2}\left[2\cdot B_{11}\cdot\epsilon_{xx}\epsilon_{xx} + 2\cdot B_{12}\cdot(\epsilon_{xx}\epsilon_{yy}) + B_{66}\cdot 0 \right]\\
&= \frac{\delta^2}{2\cdot 4}\left[2\cdot B_{11} + 2\cdot B_{12}\right]\\
&= \frac{\delta^2}{4}\left[ B_{11} + B_{12}\right]
\end{align}
$$

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

#### Triangles
Given the $\mathbf{B}$ matrix for triangles we find an energy density of
$$\begin{align}
U = \frac{\delta^2}{4}\left[ B_{11} + B_{12}\right]
\end{align}
$$

and thus $B =\frac{1}{2}(B_{11}+B_{12})$
#### Cubes
Given the $\mathbf{B}$ matrix for cubes we find an energy density of
$$\begin{align}
F &= \frac{1}{2}\left[3\cdot B_{11}\cdot\frac{\delta^2}{9} + 6\cdot B_{12}\frac{\delta^2}{9} \right]\\
&= \frac{3\delta^2}{18}\left[ B_{11} + 2\cdot B_{12}\right]\\
&= \frac{\delta^2}{6}\left[ B_{11} + 2\cdot B_{12}\right]
\end{align}
$$

So $B =\frac{1}{3}(B_{11}+2B_{12})$
#### Triangular prisms
Given the $\mathbf{B}$ matrix for triangular prisms we find an energy density of
$$\begin{align}
F &= \frac{1}{2}\left[2\cdot B_{11}\cdot\frac{\delta^2}{9} + 2\cdot B_{12}\frac{\delta^2}{9}+ 4\cdot B_{23}\frac{\delta^2}{9} +B_{33}\frac{\delta^2}{9}\right]\\
&= \frac{\delta^2}{18}\left[2\cdot B_{11} + 2\cdot B_{12}+ 4\cdot B_{23} +B_{33}\right]\\
\end{align}
$$

So $B =\frac{1}{9}(2\cdot B_{11} + 2\cdot B_{12}+ 4\cdot B_{23} +B_{33})$

## Measuring strain and stress

If we have a small enough deformation, we can approximate

$$
\frac{\partial \sigma_{\alpha\beta}}{\partial \epsilon_{\gamma\delta}}\bigg\vert_{\epsilon = 0} \approx \frac{\sigma'_{\alpha\beta}(\epsilon) - \sigma_{\alpha\beta}(\epsilon = 0)}{\epsilon_{\gamma\delta}}=B_{\alpha\beta\gamma\delta},
$$

which means that

$$
\sigma'_{\alpha\beta}(\epsilon) = \sigma_{\alpha\beta} + B_{\alpha\beta\gamma\delta}\cdot\epsilon_{\gamma\delta}
$$
## Poisson ratio's
*See Poisson ratio notebooks in /home/staff/5670772/Elastic_constants/Analysis_symmetries/*

The poisson ratio is measured by compressing a system in one direction and then measuring the strain in both that direction and a perpendicular direction. The Poisson's ratio is then defined as
$$
\nu =- \frac{\partial\epsilon_\text{trans}}{\partial\epsilon_\text{axial}}
$$
where axial is the direction of the applied stress and trans is in a direction perpendicular to that direction.
This means that for 2D there are two independent poisson ratio's, namely $\nu_{xy} = - \frac{\partial\epsilon_{xx}}{\partial\epsilon_{yy}}$ and $\nu_{yx} = - \frac{\partial\epsilon_{yy}}{\partial\epsilon_{xx}}$, while in 3D (following the same reasoning) there are 6. Note however that due to the symmetry of the stress and strain tensor they are not all independent and one can connect them by the elastic moduli (which we will not do here).

Since we want to express the Poisson ratio's in terms of the elastic constants, we use that $d\sigma_{ij} = B_{ijkl}d\epsilon_{kl}$ and the fact that we can defined compliance tensor $\mathbf{S}$, which is the inverse of the elastic constants tensor, i.e. $\mathbf{S} = \mathbf{B}^{-1}$. Together this implies that $d\epsilon_{ij} = S_{ijkl}d\sigma_{kl}$.

Finally, we use that when computing the Poisson ratio, we apply stress is in one direction. This means that both the axial and trans strain can be written in terms of the elastic constants and a stress difference in the same direction (that will trivially divide away), i.e. in general

$$
\nu_{ij} = - \frac{S_{iijj}\partial\sigma_{jj}}{S_{jjjj}\partial\sigma_{jj}}= - \frac{S_{iijj}}{S_{jjjj}}= - \frac{S_{ij}}{S_{jj}}
$$

This means that to find the poisson ratio's we only have to find the compliance tensor and compute the different ratio's .

### Squares, triangles and hexagons
Given the $\mathbf{B}$ tensor, we find
$$
\nu_{12} = \nu_{21} = \frac{B_{12}}{B_{11}}
$$
### Cubes
Given the $\mathbf{B}$ tensor, we find
$$
\nu_{12} = \nu_{21} = \nu_{13} = \nu_{31} = \nu_{23}= \nu_{32} = \frac{B_{12}}{B_{11} +B_{12}}
$$
### Triprisms and hexprisms
Given the $\mathbf{B}$ tensor, we find
$$
\begin{align}
\nu_{12} &= \nu_{21} =\frac{B^2_{23} -B_{12}B_{33}}{B^2_{23} -B_{11}B_{33}}\\
\nu_{13} &= \nu_{23} =\frac{B_{23}}{B_{11} +B_{12}}\\
\nu_{31} &= \nu_{32} =\frac{(B_{11} -B_{12})B_{23}}{-B^2_{23} +B_{11}B_{33}}
\end{align}

# Measurements specific systems

## Squares
#### Measuring $B_{11}$ and $B_{12}$
To measure $B_{11}$ and $B_{12}$ we need a stretch in the y-direction
$$
\mathbf{D}_\text{StretchY} =
\begin{pmatrix}
\frac{1}{1+\delta} & 0 &\\
0 & 1+\delta \\
\end{pmatrix}
$$
This leads to a strain tensor of
$$
\mathbf{\epsilon}_\text{StretchXY} =\frac{1}{2}(\mathbf{D}^T\mathbf{D}-1)=
\begin{pmatrix}
\frac{1}{2(1-\delta)^2}-\frac{1}{2}& 0 \\
0 & \frac{1}{2}(1+\delta)^2-\frac{1}{2}\\
\end{pmatrix}
$$
Up to order $\mathcal{O}(\delta^2)$ this becomes
$$
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
-\delta& 0\\
0& \delta\\
\end{pmatrix}
$$
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\sigma'_{ij}= \sigma_{ij}- B_{ij11}\delta+B_{ij22}\delta
$$

Since we are interested in $B_{11}= B_{1111} = B_{2222}$ we have the freedom to choose $ij= 11$ or $ij= 22$, where we choose the latter
$$
\sigma'_{22}= -P - B_{12}\delta+B_{11}\delta
$$
Or rewritten:

$$
B_{11}-B_{12} = \frac{\sigma'_{22}+P}{\delta}
$$

Given that $B=\frac{1}{2}(B_{11}+ B_{12})$, we thus find

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

#### Measuring $B_{66}$

To measure $B_{66}$, we use the deformation
$$
\mathbf{D}_\text{ShiftXY} =
\begin{pmatrix}
1 & -\delta\\
0 & 1 \\
\end{pmatrix}
$$
This leads to a strain tensor of
$$
\mathbf{\epsilon}_\text{shiftXY} =
\begin{pmatrix}
0& -\frac{\delta}{2}\\
-\frac{\delta}{2} &0 \\
\end{pmatrix} + \mathcal{O}(\delta^2)
$$
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\sigma'_{ij}(\epsilon) = \sigma_{ij}- B_{ij12}\frac{\delta}{2}-B_{ij21}\frac{\delta}{2}
$$

Since we are interested in $B_{66}= B_{1212}$, we see that we can choose $ij=12$ (where $\sigma_{ij} = -P\delta_{ij}$ falls away).

$$
\sigma'_{12}(\epsilon) = -B_{66}\delta
$$
Which we can rewrite to
$$
B_{66} = -\frac{\sigma'_{12}(\epsilon)}{\delta}
$$
#### What to measure
I.e. we need to measure
- The pressure at zero strain (to find $B$ and $P$)
- A shift in $xy$-direction
- A stretch in $y$-direction
## Triangles
#### Measuring $B_{11}$ and $B_{12}$
To measure $B_{11}$ and $B_{12}$ we need a stretch in the y-direction
$$
\mathbf{D}_\text{StretchY} =
\begin{pmatrix}
\frac{1}{1+\delta} & 0 &\\
0 & 1+\delta \\
\end{pmatrix}
$$
This leads to a strain tensor of
$$
\mathbf{\epsilon}_\text{StretchXY} =\frac{1}{2}(\mathbf{D}^T\mathbf{D}-1)=
\begin{pmatrix}
\frac{1}{2(1-\delta)^2}-\frac{1}{2}& 0 \\
0 & \frac{1}{2}(1+\delta)^2-\frac{1}{2}\\
\end{pmatrix}
$$
Up to order $\mathcal{O}(\delta^2)$ this becomes
$$
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
-\delta& 0\\
0& \delta\\
\end{pmatrix}
$$
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\sigma'_{ij}= \sigma_{ij}- B_{ij11}\delta+B_{ij22}\delta
$$

Since we are interested in $B_{11}= B_{1111} = B_{2222}$ we have the freedom to choose $ij= 11$ or $ij= 22$, where we choose the latter
$$
\sigma'_{22}= -P - B_{12}\delta+B_{11}\delta
$$
Or rewritten:

$$
B_{11}-B_{12} = \frac{\sigma'_{22}+P}{\delta}
$$

Given that $B=\frac{1}{2}(B_{11}+ B_{12})$, we thus find

$$
B_{11} =B+\frac{1}{2}\left(\frac{\sigma'_{22}+P}{\delta} \right) \text{ and } B_{12} = B-\frac{1}{2}\left(\frac{\sigma'_{22}+P}{\delta} \right)
$$
#### Measuring $B_{66}$
We can find $B_{66}$ using that $B_{66}= \frac{B_{11}-B_{12}}{2}$
#### What to measure
I.e. we need to measure
- The pressure at zero strain (to find $B$ and $P$)
- A stretch in $xy$-direction
## Hexagons
Since the elastic constant tensor is the same for triangles and hexagons, we can do the same measurements.
## Cubes
#### Measuring $B_{11}$ and $B_{12}$
To measure $B_{11}$ and $B_{12}$ we need two deformations
$$
\mathbf{D}_\text{StretchXY} =
\begin{pmatrix}
1+\delta & 0 & 0\\
0 & 1-\delta & 0\\
0 & 0 & \frac{1}{1-\delta^2}
\end{pmatrix}
$$
This leads to a strain tensor of
$$
\mathbf{\epsilon}_\text{StretchXY} =\frac{1}{2}(\mathbf{D}^T\mathbf{D}-1)=
\begin{pmatrix}
\frac{\delta^2}{2} +\delta& 0 & 0\\
0 & \frac{\delta^2}{2} -\delta & 0\\
0 & 0 & \frac{1}{2(1-\delta^2)^2}-\frac{1}{2}
\end{pmatrix}
$$
Up to order $\mathcal{O}(\delta^2)$ this becomes
$$
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & -\delta & 0\\
0 & 0 & 0
\end{pmatrix}
$$
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\sigma'_{ij}(\epsilon) = \sigma_{ij}+ B_{ij11}\delta-B_{ij22}\delta
$$

Since we are interested in $B_{11}= B_{1111} = B_{2222}= B_{3333}$ and $B_{22}= B_{1122} = B_{1133}= B_{2233}$ we have the freedom to choose $ij= 11$ or $ij= 22$, where we choose the latter

$$
\sigma'_{22}(\epsilon) = -P+ B_{12}\delta-B_{11}\delta
$$
Or rewritten:

$$
B_{11}-B_{12} = - \frac{\sigma'_{22}+P}{\delta}
$$

Given that $B=\frac{1}{3}(B_{11}+ 2B_{12})$, we thus find

$$
B_{11} = B-\frac{2}{3}\left(\frac{\sigma'_{22}+P}{\delta} \right) \text{ and } B_{12} = B+\frac{1}{3}\left(\frac{\sigma'_{22}+P}{\delta} \right)
$$

#### Measuring $B_{66}$

To measure $B_{66}$, we use the deformation
$$
\mathbf{D}_\text{ShiftXZ} =
\begin{pmatrix}
1 & 0 & -\delta\\
0 & 1& 0\\
0 & 0 & 1
\end{pmatrix}
$$
This leads to a strain tensor of
$$
\mathbf{\epsilon}_\text{ShiftXZ} =
\begin{pmatrix}
0& 0 & -\frac{\delta}{2}\\
0 &0 &0 \\
-\frac{\delta}{2} & 0 & 0
\end{pmatrix} + \mathcal{O}(\delta^2)
$$
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\sigma'_{ij}(\epsilon) = \sigma_{ij}- B_{ij13}\frac{\delta}{2}-B_{ij31}\frac{\delta}{2}
$$

Since we are interested in $B_{66}= B_{1313} =B_{1331}= B_{1212}$, we see that we can choose $ij=13$ (where $\sigma_{ij} = -P\delta_{ij}$ falls away).

$$
\sigma'_{13}(\epsilon) = -B_{1313}\delta = -B_{66}\delta
$$
Which we can rewrite to
$$
B_{66} = -\frac{\sigma'_{13}(\epsilon)}{\delta}
$$
#### What to measure
I.e. we need to measure
- The pressure at zero strain (to find $B$ and $P$)
- A shift in $xz$-direction
- A stretch in $xy$-direction

## Triangular prisms
*N.B. to indicate how we measure the stress tensor we use SH for shift and ST for stretch, followed by ij that indicate the direction*

#### Measuring $B_{11}$ and $B_{12}$ and $B_{23}$
*See cubes for further derivation for the stretch matrices.* We have 4 unknown constants that we want to measure via stretches: $B_{11},B_{12},B_{23}$ and $B_{33}$.

so we need to do two more deformation, that we choose here as a stretch in the $xy$-direction and in the $xz$-direction.
$$
\mathbf{D}_\text{StretchXY} =
\begin{pmatrix}
1+\delta & 0 & 0\\
0 & 1-\delta & 0\\
0 & 0 & \frac{1}{1-\delta^2}
\end{pmatrix}
\rightarrow
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & -\delta & 0\\
0 & 0 & 0
\end{pmatrix}
$$
$$
\mathbf{D}_\text{StretchXZ} =
\begin{pmatrix}
1+\delta & 0 & 0\\
0 & \frac{1}{1-\delta^2}& 0\\
0 & 0 & 1-\delta
\end{pmatrix}
\rightarrow
\mathbf{\epsilon}_\text{StretchXZ} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & 0 & 0\\
0 & 0 & -\delta
\end{pmatrix}
$$
and a deformation only in the $x$-direction
$$
\mathbf{D}=
\begin{pmatrix}
1+\delta& 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}
$$
which leads to a stress tensor of
$$
\mathbf{\epsilon}_\text{StretchX} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & 0 & 0\\
0 & 0 & 0
\end{pmatrix}
$$
::: warning
Can we actually measure this already with the current code?
> The code can measure these without trouble, as it is just a change in initial configuration. You will have to implement the initial configuration for triangular prisms, though.
:::
Filling this in in our earlier found equation for $\sigma_{ij}'$ gives
$$
\begin{align}
\sigma'^{ST12}_{ij}(\epsilon) &= \sigma_{ij}+ B_{ij11}\delta-B_{ij22}\delta\\
\sigma'^{ST13}_{ij}(\epsilon) &= \sigma_{ij}+ B_{ij11}\delta-B_{ij33}\delta\\
\sigma'^{ST1}_{ij}(\epsilon) &= \sigma_{ij}+ B_{ij11}\delta
\end{align}
$$

We choose the $ij$ such that they are convenient for us
\begin{align}
\sigma'^{ST12}_{22}(\epsilon) &= -P + B_{2211}\delta-B_{2222}\delta =-P + B_{12}\delta-B_{11}\delta\\
\sigma'^{ST13}_{22}(\epsilon) &= -P + B_{2211}\delta-B_{2233}\delta= -P + B_{12}\delta-B_{23}\delta\\
\sigma'^{ST13}_{33}(\epsilon) &= -P + B_{3311}\delta-B_{3333}\delta= -P + B_{23}\delta-B_{33}\delta\\
\sigma'^{ST1}_{11}(\epsilon) &= -P + B_{1111}\delta= -P + B_{11}\delta
\end{align}
Solving this (using mathematica, *see notebook Symmetries triangular prism.nb*) gives us
$$
\begin{align}
B_{11} &=\frac{P+\sigma'^{ST1}_{11}}{\delta}\\
B_{12} &=\frac{2P+\sigma'^{ST1}_{11}+\sigma'^{ST12}_{22}}{\delta}\\
B_{23} &=\frac{P+\sigma'^{ST1}_{11}+\sigma'^{ST12}_{22}-\sigma'^{ST13}_{22}}{\delta}\\
B_{33} &=\frac{\sigma'^{ST1}_{11}+\sigma'^{ST12}_{22}-\sigma'^{ST13}_{22}-\sigma'^{ST13}_{33}}{\delta}
\end{align}
$$

#### Measuring $B_{55}$
*See cubes for further derivation for the shift matrices.* To measure $B_{55}$, we use a shift in the $xz$-direction that leads to
$$
\sigma'^{SH13}_{13}(\epsilon) = -B_{1313}\delta = -B_{55}\delta
$$
Which we can rewrite to
$$
B_{55} = -\frac{\sigma'^{SH13}_{13}}{\delta}
$$
#### Measuring $B_{66}$
To measure $B_{66}$ we can use that
$$B_{66}= \frac{B_{11}-B_{12}}{2}$$

#### What to measure
I.e. we need to measure
- The pressure at zero strain (to find $B$ and $P$)
- A shift in $xz$-direction
- A stretch in $xy$-direction
- A stretch in $xz$-direction

### Triangles in a hexagonal crystal

Since the EDMD code had not yet a triangular lattice I implemented that. Below I show the results for the normal system, the system with a shift and the system with a stretch. Particle colours are purely indicative.
| | $t/\tau=0$ | $t/\tau\approx20$ |
| --- | -------------------------------------------------------------------------------- | -------- |
|Normal | |  |
| Shift $4.30$ |  ||
| Stretch $0.1$ |  ||

## Hexagonal prisms
Since the elastic constant tensor is the same for triangular prisms and hexagonal prisms, we can do the same measurements.

## Error propagation
### Pressure
Since $P = \frac{1}{3}Tr[\mathbf{P}]$, for the pressure the error is given by $\sigma_P = \frac{1}{3} \sqrt{\sigma^2_{P_{xx}}+ \sigma^2_{P_{yy}}+ \sigma^2_{P_{zz}}}$.

# Results
## Squares
For a system of $N=1024$ and densities $\rho \in(0.75, 0.8, 0.84, 0.85, 0.9, 0.91, 0.925, 0.94, 0.95, 0.96)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.05 0.1)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 500.000$ timesteps.

### Evaluation data

As it turns out, sometimes the square system turns,
|$t<t_{shift}$ | $t>t_{shirt}$|
| -------- | -------- |
|  ||

To see how often and when this happens we plot the mean orientation of the particles as a function of time. We do not correct for the fact that particles have a symmetry of $\pi/2$, since we are also interested in the direction particles turn to.

|$\rho =0.75$ |$\rho =0.80$|
| -------- | -------- |
| ||
|$\rho =0.85$ |$\rho =0.90$|
|||

#### Shift turn
|$\rho =0.75$ |$\rho =0.80$|
| -------- | -------- |
|||
|$\rho =0.85$ |$\rho =0.90$|
|||

Individual particles $\rho = 0.80$ stretch

| $\delta = 0.001$ | $\delta = 0.009$ | $\delta = 0.01$ |
| -------- | -------- | -------- |
| | | |

Individual particles $\rho = 0.85$ stretch
| $\delta = 0.001$ | $\delta = 0.01$ |
| -------- | -------- |
| | |

At $\eta = 0.96$ and $\delta = 0.01$ we even get some sot of density wave through the system.
 Stretch $\eta = 0.96$, $\delta = 0.01$.

**From the analysis (see notebook Analysis_angle_nomean.nb) it becomes clear that from these densities we can use $0.86<\rho<0.95$, using the interval $0<\delta<0.008$. Shift has a slightly bigger range, but to stay at the save side we can use this.**

### Stretch

| $\beta( B_{11}-B_{12})\sigma^2$ |Poisson ratio fit $\delta < 0.01$|
| -------- | -------- |
|||
|$\beta B_{11}\sigma^2$ |$\beta B_{12}\sigma^2$ |
|||

Then we can measure the Poisson ratio as a function of $\rho$.

::: danger
What do we fit? It feels like sometimes it works better to fit with higher orders of delta too.

| Fit up to $\delta$ | Fit up to $\delta^3$ |
| -------- | -------- |
| | |

:::
### Shift
Note that we fit the shift data with the equation $y = a\delta^2 +c$. The reason is that the system should be symmetric under shift (i.e. $+\delta$ should give the same results as $-\delta$, therefore the system can not linearily depend on $\delta$).
::: info
Note that for the lower densities the influence of small deviations in the pressure measurements probably starts playing a role, hence the wiggling.
:::

## Triangles
### Adjusting the code
*Found on HERA /home/staff/5670772/Elastic_constants/CODE_polyhedra_EDMD_2D/CHECKS/Triangular_lattice/ the mov files hold the configurations associated with the snapshot below.*

For a system of $N=1008$ and densities $\rho \in(0.75, 0.8, 0.85 0.9)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 500.000$ timesteps.

### Evaluation data

We know that systems of triangles have a triatic phase between $\eta = 0.70$ and $\eta= 0.90$ (https://pubs.rsc.org/en/content/articlehtml/2015/sm/c5sm01762a). Below we have a fluid phase (continuous phase transition) and above we have a chiral phase. To check whether our system indeed has a nice triatic phase, we do an analysis of the angle of our triangles (see notebook *Analysis_angle.nb*). In the figure we plot the orientation of a random selection of particles over tie, black lines are a multiple of $i\cdot\pi/3 - \pi/6$, i.e. the symmmetries of our particle)

| Stretch $\eta = 0.75$, $\delta = 0.001$ | Stretch $\eta = 0.80$, $\delta = 0.01$ | $\eta = 0.90$, $\delta = 0.01$ |
| -------- | -------- | -------- |
| | | |

Note that these results are for the stress tensor, however, we see the same behaviour for shift. As we can see, for the intermediate packing fractions, the orientation of the particles stays constant. However, for the lower packing fractions the particles rotate a lot. Moreover, for the highest packing fraction, we see a constant deviation of the black line. This means that all our particles are rotated a little, i.e. most likely we are in a chiral phase.

**Conclusion: we can only use values between $\eta = 0.80$ and $\eta = 0.86$ for our analysis. To make sure that we don't include too large deformations, and at the same time don't have a too large influence of the pressure deviations, we only include $0.004 \leq \delta \leq 0.009$**

Note that in *Analysis_stresstensorvalues.nb* we also looked a the values of the stress tensor, these look nice (meaning that there are no weird spikes and the data seems to be well equilibrated. )
### Stretch

| $\beta( B_{11}-B_{12})\sigma^2$ |Poisson ratio|
| -------- | -------- |
|||
|$\beta B_{11}\sigma^2$ |$\beta B_{12}\sigma^2$ |
|||

We obtain the folloing Poisson ratio

### Shift
The elastic constant $\beta B_{66}\sigma^2$ we can measure directly and obtain using $B_{66} = ( B_{11} - B_{12} )/2$. On the right the open symbols represent $B_{66} = ( B_{11} - B_{12} )/2$, while the closed points are measured directly. As we can see, when including only higher $\delta$'s' they converge to the same point when $\delta\rightarrow 0$.

| $\beta B_{66}\sigma^2$ |Compare|
| -------- | -------- |
|||

## Hexagons
For a system of $N=1024$ and densities $\rho \in(0.85,0.87,0.89, 0.91)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 500.000$ timesteps.

### Evaluation data
For this simulation by accident we only printed the angle of particles every 10.000 steps, which means that a thorough analysys of the angle is hard.
::: danger
We do see that the data 'wiggles' around the equilibritum values, even for these high densities. I think we can almost not get rid of this, and the data looks fine. Therefore I think we leave it

:::

We do see however that for the higher $\delta's$ the shift simulations do not work properly (I think the shift is too strong for these high densities). Therefore we throw all data that has $\delta >0.008$ out.

At the same time, for the lower $\delta$'s, we see that the pressure influences the results, therefore we also throw out everything below $\delta < 0.003$

**Conclusion: we can use all densities that we measured and $0.002 < \delta < 0.09$**

### Stretch

| $\beta( B_{11}-B_{12})\sigma^2$ |Poisson ratio fit $\delta < 0.01$|
| -------- | -------- |
|||
|$\beta B_{11}\sigma^2$ |$\beta B_{12}\sigma^2$ |
|||

Poisson ratio
Then we can measure the Poisson ratio as a function of $\rho$.

### Shift

| $B_{66}$ function $\delta$| $B_{66}$ function $\rho$ |
| -------- | -------- |
| ||

Where again dotted lines correspond to the $B_{66}$ obtained from $B_{11}$ and $B_{12}. Note that we fit the shift data with the equation $y = a\delta^2 +c$. The reason is that the system should be symmetric under shift (i.e. $+\delta$ should give the same results as $-\delta$, therefore the system can not linearily depend on $\delta$).

Note that if we interpolate the values to $\delta = 0$ both methods give exactly the same results (see also plot on the right).

## Cubes
For a system of $N=1000$ and densities $\rho \in(0.65,0.70,0.75, 0.8, 0.85 0.9)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 50.000$ timesteps. We measure the simulations using the Andersen thermostat.

::: danger
Why sometimes runtime is longer than what I put in Franks codes?
:::
### Stretch and shift
| $\beta( B_{11}-B_{12})\sigma^2$ |$B_{66}$ |
| -------- | -------- |
|||
|Poisson ratio| Poisson ratio func density, only fit $\delta>0.005$ |
|||
|$\beta B_{11}\sigma^2$ |$\beta B_{12}\sigma^2$ |
|||

::: warning
We get other values than Laura does....
:::
## Triangular prisms

*See HERA /home/staff/5670772/Elastic_constants/DATA_triprisms/*

For a system of $N=1344$ and densities $\rho \in(0.65,0.70,0.725, 0.75,0.775 0.8,0.825, 0.85, 0.9)$ and $\delta=(0.0005, 0.001,0.0015, 0.002, 0.0025,0.003,0.0035, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 50.000$ timesteps.

### Hexagonal spacking $z$-direction

First we need to measure the hexagonal spacing for a a hexagonal crystal as a function of the density. Therefore we measure the pressure tensor for various hexagonal spacing and check for which spacing it becomes isotropic.
::: warning
Some of the hex simulations also stopped too early because of the gridsize (see the notebook Analysis_hexspacing). However, since even with lightly less datapoints, I trust the data enough to use those points (in the end one of the datapoints is used to fit a line through for which we compute the intersection with the line $y =1$. If one of the points is based on slightly less data, this is fine).
:::
| $\sigma_{11}/\sigma_{22}$ | $\sigma_{11}/\sigma_{33}$ with fit | Hexagonal spacing|
| ------------------------- | --- | -------- |
| ||  |

As we can see for some hexagonal spacings the $x$- and $y$-directions are not equal to each other anymore. This means that the crystal has probably shifted. We remove these points from the analysis.

The hexagonal spacing is given by

## Hexagonal prisms

*See HERA /home/staff/5670772/Elastic_constants/DATA_hexprisms/*

For a system of $N=1080$ and densities $\rho \in(0.70, 0.75,0.8, 0.85, 0.9)$ and $\delta=( 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 250.000$ timesteps.

### Hexagonal spacking $z$-direction

First we need to measure the hexagonal spacing for a a hexagonal crystal as a function of the density. Therefore we measure the pressure tensor for various hexagonal spacing and check for which spacing it becomes isotropic.
::: warning
Some of the hex simulations also stopped too early because of the gridsize (see the notebook Analysis_hexspacing). However, since even with lightly less datapoints, I trust the data enough to use those points (in the end one of the datapoints is used to fit a line through for which we compute the intersection with the line $y =1$. If one of the points is based on slightly less data, this is fine).
:::
| $\sigma_{11}/\sigma_{22}$ | $\sigma_{11}/\sigma_{33}$ with fit | Hexagonal spacing|
| ------------------------- | --- | -------- |
|| | |

As we can see for some hexagonal spacings the $x$- and $y$-directions are not equal to each other anymore. This means that the crystal has probably shifted. We remove these points from the analysis.

The hexagonal spacing is given by
$$
\text{hexspacing}/\sigma = 0.420269 - 0.431629 \rho\sigma^2
$$

::: danger
If we run the pressure analysis for this spacing, it doesn't look too isotropic....

Maybe look per point?? or fit $x, x^2, x^3$ through it?
If we use a higher order fit we obtain this

:::

## Understand Poisson constants in 3D
We have the following Poisson constants in the case of the triangular and hexegonal crystals
$$
\begin{align}
\nu_{12} &= \nu_{21} =\frac{B^2_{23} -B_{12}B_{33}}{B^2_{23} -B_{11}B_{33}}\\
\nu_{13} &= \nu_{23} =\frac{(B_{11} -B_{12})B_{23}}{-B^2_{23} +B_{11}B_{33}}\\
\nu_{31} &= \nu_{32} =\frac{B_{23}}{B_{11} +B_{12}}\\
\end{align}
$$

Where we know that because of stability
- $B_{11}, B_{33} > 0$.
- Most likely $|B_{11}|> |B_{12}|$.

This means that we get the following contributions that are surely positive or negative. This also imply that the sign of $B_{12}$ and $B_{23}$ determine almost everything.
$$
\begin{align}
\nu_{12} &= \nu_{21} =\frac{\color{green}{B^2_{23}} - B_{12}\color{green}{B_{33}}}{\color{green}{B^2_{23}} -\color{green}{B_{11}}\color{green}{B_{33}}}\\
\nu_{13} &= \nu_{23} =\frac{(\color{green}{B_{11}} -B_{12})B_{23}}{-\color{green}{B^2_{23}} +\color{green}{B_{11}}\color{green}{B_{33}}}\\
\nu_{31} &= \nu_{32} =\frac{B_{23}}{\color{green}{B_{11}} +B_{12}}\\
\end{align}
$$

Since for both prisms we expect and see that $B_{12} > 0$, $\nu_{12}$ will be almost always positive. The only exception is when $B_{12} \approx 0$. In that case we see that the upper part of the fraction is positive, while the lower part is negative. I think at this point, the wedging effect is almost zero, which means that the allignment that the particles undergo in the $z$-direction and as a result therefore the allignment that the particles undergo in the $y$-direction dominates. This only happens for hexagonal prisms, because the wedging effect is very small there.

For $\nu_{13}$, since in general $B_{11}B_{33}> B^3_{23}$, the lower part of the fraction will be positive. The upper part will be negative, since we found that $B_{23} < 0$ and $B_{11} - B_{12}>0$. Therefore the constant is negative.

For $\nu_{31}$ the lower part will be positive, while the upperpart is negative, i.e. the Poisson ratio is negative.

Since for the last Poisson ratio's it is mostly the sign that dominates and not the size, we can with certainty say that the Poisson ratio's are negative. For the first one, that is not the case, so there it could be it is negative instead of positive, or vice versa.
## Results


## Thorough check
- [x] Squares
- [x] Theory
- [x] Data complete
*Some pressure tensor measurements are stopped a slightly early (length list 219544 vs 222208), however this difference will have no influence. The same is true for the shift and stretch files.*
- [x] Stress tensor
- [x] Mathematica notebook
- [x] Triangles
- [x] Theory
- [x] Data complete
*Same as squares, Some pressure tensor measurements are stopped a slightly early, however this difference will have no influence. The same is true for the shift and stretch files.*
- [x] Stress tensor
- [x] Mathematica notebook
- [x] Hexagons
- [x] Theory
- [x] Data complete
*Same as squares, Some pressure tensor measurements are stopped a slightly early, however this difference will have no influence. The same is true for the shift and stretch files.*
- [x] Stress tensor
- [x] Mathematica notebook
- [ ] Cubes
- [ ] Theory
- [ ] Data complete
- Data pressure is still running for the highest packing fractions
- [ ] Stress tensor
- [ ] Mathematica notebook
- [ ] Triprisms
- [ ] Theory
- [ ] Data complete
- [ ] Stress tensor
- [ ] Mathematica notebook
- [ ] Stresstensor pressure isotropic
- [ ] Hexprisms
- [ ] Theory
- [ ] Data complete
- [ ] Stress tensor
- [ ] Mathematica notebook
- [ ] Stresstensor pressure isotropic

# To do myself
## To do new code:
- Cubes:
- [ ] Run pressure
- [ ] Analyse everything
- [ ] See if no spikes
- Triangular prisms:
- [ ] Run hexspacing
- [ ] Run pressure
- [ ] Run shift
- [ ] Run stretch
- [ ] Analyse everything
- Hexagonal prisms:
- [ ] Make hexagonal lattice
- [ ] Run hexspacing
- [ ] Run pressure
- [ ] Run shift
- [ ] Run stretch
- [ ] Analyse everything

:::info
Stress is pressure, since we are in an equilibrium the applied force must be counteracted by the pressure, so they must balance each other. Since there is no torque (we only have deformation, we loose 3 independt entries: the tensor must be symmetric. )
:::

## Runs that run too long

We stopped these simulations because then ran for too long. Probably something went wrong with them. Highlighted simulationes are started again (15-4-2024).

Find de propruns using

find /home/staff/5670772/Elastic_constants/ -type f -iname PropRun1.o3325367

==3325369: DATA_squares/shift/run_0.85_0.01==
==3325375: DATA_triangles/shift/run_0.75_0.002==
==3325427: DATA_triangles/stretch/run_0.85_0.004==
==3325481: DATA_triangles/pressure/run_0.88==
==3325484: DATA_triangles/pressure/run_0.91==
3326049: DATA_triprisms_spikecode/hexspacing/run_0.70_0.002
3326051-3326058: DATA_triprisms_spikecode/hexspacing/run_0.(74-78)_0.002
3326063-3326070: DATA_triprisms_spikecode/hexspacing/run_0.(76-90)_0.002
3326079-3326081: DATA_triprisms_spikecode/hexspacing/run_0.(8690)_0.004
3326083:DATA_triprisms_spikecode/hexspacing/run_0.72_0.005
3326093-3326094:DATA_triprisms_spikecode/hexspacing/run_0.(70-72)_0.01/
3326179: DATA_cubes_spikecode/pressure/run_0.91
3326274: DATA_cubes_spikecode/stretch/run_0.85_0.003
3326346: DATA_cubes_spikecode/shiftx/run_0.85_0.001
3326358: DATA_cubes_spikecode/shiftx/run_0.85_0.004
3326363: DATA_cubes_spikecode/shiftx/run_0.90_0.005/
3326405: DATA_triprisms_spikecode/hexspacing/run_0.74_0.02/
3326433: DATA_triprisms_spikecode/hexspacing/run_0.86_0.06/
3326446: DATA_triprisms_spikecode/hexspacing/run_0.90_0.08
3326886: DATA_triprisms_spikecode/pressure/run_0.95/
3326908: DATA_triprisms_spikecode/shiftx/run_0.65_0.005
3326997: DATA_triprisms_spikecode/stretchxz/run_0.65_0.009/
3327061: DATA_triprisms_spikecode/stretchx/run_0.80_0.001/
3327326: DATA_triprisms_spikecode/stretchxz/run_0.70_0.0005/
3327331: DATA_triprisms_spikecode/stretchxz/run_0.70_0.0035
3327611: DATA_triprisms_spikecode/stretchxy/run_0.825_0.0005
3327645: DATA_triprisms_spikecode/stretchx/run_0.725_0.006
3327662: DATA_triprisms_spikecode/stretchxy/run_0.725_0.0005/
3327700: DATA_triprisms_spikecode/stretchx/run_0.825_0.01/
3327706: DATA_triprisms_spikecode/shiftx/run_0.825_0.001
3327722: DATA_triprisms_spikecode/shiftx/run_0.725_0.004
3328130: DATA_triprisms_spikecode/stretchxy/run_0.875_0.001/
3328136: DATA_triprisms_spikecode/stretchxy/run_0.875_0.0025
3328150: DATA_triprisms_spikecode/stretchxy/run_0.875_0.008
3328199: DATA_triprisms_spikecode/shiftx/run_0.9_0.0035/
3328201: DATA_triprisms_spikecode/shiftx/run_0.9_0.004/
==3328237: DATA_squares/shift/run_0.94_0.002==
==3328271: DATA_squares/shift/run_0.84_0.008==
3328811: DATA_triprisms_spikecode/pressure/run_0.875/
==3328906: DATA_triangles/stretch/run_0.85_0.002==
==3338798: DATA_hexagons/shift/run_0.85_0.009==
==3338802: DATA_hexagons/shift/run_0.85_0.01==
==3338803: DATA_hexagons/shift/run_0.87_0.01==
==3338805: DATA_hexagons/shift/run_0.91_0.01==

# Old Poisson ratio's

## Cubes
For a system of $N=1000$ and densities $\rho \in(0.65,0.70,0.75, 0.8, 0.85 0.9)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 50.000$ timesteps. We measure the simulations using the Andersen thermostat.

### Stretch and shift
| $\beta( B_{11}-B_{12})\sigma^2$ |$B_{66}$ |
| -------- | -------- |
|||
| Poisson ratio func $\delta$ | Poisson ratio func density, only fit $\delta>0.005$ |
|||
|$\beta B_{11}\sigma^2$ |$\beta B_{12}\sigma^2$ |
|||
::: danger
Again there are spikes int he pressure... Filter them out?
 We see them mainly for $\rho = 0.80$
:::
:::info
To do:
- Run longer simulations
- Which values to include (which densities and which delta's). And how to find out, because angle analysis is less trivial.
- Look at SD's of the stress tensor
:::
::: warning
We get other values than Laura does....
:::
## Triangular prisms

*See HERA /home/staff/5670772/Elastic_constants/DATA_triprisms/*

For a system of $N=1344$ and densities $\rho \in(0.65,0.70,0.725, 0.75,0.775 0.8,0.825, 0.85, 0.9)$ and $\delta=(0.0005, 0.001,0.0015, 0.002, 0.0025,0.003,0.0035, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$ we measure the stress tensor in case of shift and in case of stretch. Note that we measure for $t/\tau = 50.000$ timesteps.

### Hexagonal spacking $z$-direction

First we need to measure the hexagonal spacing for a a hexagonal crystal as a function of the density. Therefore we measure the pressure tensor for various hexagonal spacing and check for which spacing it becomes isotropic.
| $\sigma_{11}/\sigma_{22}$ | $\sigma_{11}/\sigma_{33}$ with fit | Hexagonal spacing|
| ------------------------- | --- | -------- |
| | |  |

As we can see for some hexagonal spacings the $x$- and $y$-directions are not equal to each other anymore. This means that the crystal has probably shifted. We remove these points from the analysis.

The hexagonal spacing is given by
$$
\text{hexspacing} = 0.440489 - 0.451303 \cdot \rho
$$

::: warning
We again see peaks
 $\rho = 0.85$, stretchxy
:::
### Results
::: warning
Weird wiggles

| $B_{11}$ | $B_{12}/B_{22}$ | $B_{12}/B_{33}$ |
| -------------------------------------------------------------------------------- | -------------------------------------------------------------------------------- | --- |
|  |  |  |

:::

| $B_{11}$ | $B_{12}$ | $B_{23}$ |
| -------- | -------- | -------- |
|| | |
| | $B_{12}$ zoomed in | $B_{23}$ zoomed in|
||||
| $B_{33}$ | $B_{55}$ | $B_{66}$ |
||| |
| $B_{33}$ zoomed in | | |
||||
### Poisson ratio's
| $\nu_{12}=\nu_{21}$ |$\nu_{13}=\nu_{23}$ |$\nu_{31}=\nu_{32}$|
| -------- | -------- | -------- |
| |||

::: danger
Fit with something different than $a + bx$, since it seems that higher order effect could play a role, howere maybe we have too little data. We could also wait for more data.
:::
All Poisson ratio's

## Obtain poisson ratio's from experiments
To obtain the Poisson ratio's from experiments, we need a way to obtain the elastic constants from correlations between local deformations. For this we first need some background on DFT.

### DFT order parameter
*Obtained from https://www.educare.bz/wp-content/uploads/2018/02/Chaikin-Lubensky-Principles-of-Condensed-Matter-Physics.jb2_.pdf*

Assume we have some order parameter $\phi(\mathbf{x})$ and an external field $h(\mathbf(x))$. In general $h(\mathbf{x})$ could be a real field working on the order parameter (like an B-field in a spin system that couples to the direction of the spin), however in our case it is an auxiliry field. We can then modify our hamiltonian to
$$
\mathcal{H} = \mathcal{H}_{K} + U + U_{ext}
$$
where
$$
U_{ext} \equiv -\int d^dx \phi(\mathbf{x})h(\mathbf{x}),
$$

Given this new potential, our F is now a functional of $\phi(\mathbf{x})$
$$
Z_N[T, V, [\phi(x)]] = Tr[e^{-\beta\mathcal{H}}]
$$
Note that $\phi(\mathbf{x})$ is essentially a Lagrange multiplier, meaning that if we derive $Z_N[T, V, [h(x)]]$ w.r.t. $h(\mathbf{x})$ we obtain the average order parameter $\phi(\mathbf{x})$ (just like deriving it w.r.t $\beta$ gives us the average energy). The only main difference is that $h(\mathbf{x})$ is a scalar **field** instead of a scalar. Given this $Z_N[T, V, [h(x)]]$, we can find
$$
\frac{1}{\delta Z_N}\frac{\delta Z_N}{\beta\delta h (\mathbf{x})} = \langle \phi(\mathbf{x})\rangle \text{ and }\frac{1}{\delta Z_N}\frac{\delta^2 Z_N}{\beta\delta h (\mathbf{x})\beta\delta h (\mathbf{x'})} = \langle \phi(\mathbf{x})\phi(\mathbf{x'})\rangle
$$
such that the fluctuations in $\langle \delta \phi(\mathbf{x})\delta \phi(\mathbf{x'})\rangle$ are given by
$$
\begin{align*}
S_{\phi_\mathbf{x}\phi_\mathbf{x'}}&=\langle [\phi(\mathbf{x}) - \langle \phi(\mathbf{x})\rangle][\phi(\mathbf{x'}) - \langle \phi(\mathbf{x'})\rangle] \\
&= \langle \phi(\mathbf{x})\phi(\mathbf{x'})\rangle - \langle \phi(\mathbf{x})\rangle\langle \phi(\mathbf{x}')\rangle\\
&= \frac{1}{\delta Z_N}\frac{\delta^2 Z_N}{\beta\delta h (\mathbf{x})\beta\delta h (\mathbf{x'})} -\frac{1}{Z_N^2}\frac{\delta Z_N}{\beta\delta h (\mathbf{x})}\frac{\delta Z_N}{\beta\delta h (\mathbf{x}'')}\\
&= \frac{\delta^2 \ln Z_N}{\beta\delta h (\mathbf{x})\beta\delta h (\mathbf{x'})} \\
& =\frac{\delta \langle \phi(\mathbf{x'})\rangle}{\beta \delta h (\mathbf{x}) }\equiv \frac{\chi(\mathbf{x}, \mathbf{x}'),}{\beta},
\end{align*}
$$
with $\chi$ the susceptibility. However, often is more useful to express things in terms of the order parameter instead of the field. Therefore we can do a Legendre trasform
$$
F_\phi[T, \langle\phi(\mathbf{x})\rangle] = F_h[T, h(\mathbf{x})] +\int d^dx h(\mathbf{x})\phi(\mathbf{x})
$$
and thus $\delta F_\phi/\delta \phi(\mathbf{x}) = h(\mathbf{x})$. Now we want to write the susceptibility in terms of derivitaves to $\phi(\mathbf{x})$ instead of $h(\mathbf{x})$. Therefore we first consider
$$
\begin{align*}
\frac{\delta\phi(\mathbf{x})}{\delta\phi(\mathbf{x}'')} &=\delta(\mathbf{x} - \mathbf{x}'')\\
&= \int d^dx\frac{\delta\phi(\mathbf{x}')}{\delta h(\mathbf{x}')} \frac{\delta h(\mathbf{x}')}{\delta\phi(\mathbf{x}'')}\\
&\equiv \int d^dx\chi(\mathbf{x}, \mathbf{x}')\chi^{-1}(\mathbf{x}', \mathbf{x}''),
\end{align*}
$$
where in the last step we defined the inverse $\chi^{-1}(\mathbf{x}', \mathbf{x}'')$. As a result
$$
\begin{align*}
\chi^{-1}(\mathbf{x}, \mathbf{x}')&=\frac{\delta h(\mathbf{x}')}{\delta\phi(\mathbf{x}'')}\\
&=\frac{\delta^2 F_\phi}{\delta \phi(\mathbf{x} )\phi(\mathbf{x'})},
\end{align*}
$$
where in the last step we used that $\delta F_\phi/\delta \phi(\mathbf{x}) = h(\mathbf{x})$.

In general the $F_\phi$ will not depend trivially on $\phi$, however we can do a mean field approximation in $\phi$. Assume that in the high temperature limit the mean $\phi$ will on average be zero and that variations in $\phi$ will cost energy, such that

$$
F = \int d^d x[f(T, \langle\phi (\mathbf{x})\rangle) + \frac{1}{2}c[\nabla\langle\phi(\mathbf{x})\rangle ]^2]
$$
with $f(T, \langle\phi (\mathbf{x})\rangle)$ a local free energy density that we can expand like
$$
f(T, \langle\phi (\mathbf{x})\rangle) = a\phi + \frac{1}{2} r \phi^2 - w\phi^3 + u\phi^4 +\mathcal{O}(\phi^5)
$$
Since for high temperatures $\langle(\phi)$ should be zero when $h$ is zero (i.e. $\frac{\delta f}{\delta \phi} |_{\phi = 0} = h = 0$) the linear terms should be zero. Moreover since $\langle\phi\rangle \rightarrow \langle-\phi\rangle$ all odd terms should be zero too. Therefore
$$
f(T, \langle\phi (\mathbf{x})\rangle) = \frac{1}{2} r \phi^2 + u\phi^4 +\mathcal{O}(\phi^5)
$$

To find the susceptibility we first need to take the functional derivative of $F_\phi$ w.r.t to $\phi(\mathbf{x})$, for this we use that
$$
\begin{align*}
\frac{\delta F(\phi(\mathbf{x}),\nabla \phi(\mathbf{x}))}{\delta \phi(\mathbf{y})} &= \int d^dx \frac{\partial F}{\partial \phi(\mathbf{y})}\\
&= \int d^dx [\frac{\partial F}{\partial \phi(\mathbf{x})}\frac{\partial \phi(\mathbf{x})}{\partial \phi(\mathbf{y})} + \frac{\partial F}{\partial \nabla\phi(\mathbf{x})}\frac{\nabla \partial \phi(\mathbf{x})}{\partial \phi(\mathbf{y})}]\\
&= \int d^dx [\frac{\partial F}{\partial \phi(\mathbf{x})}\delta(\mathbf{x}-\mathbf{y}) + \frac{\partial F}{\partial \nabla\phi(\mathbf{x})}\nabla\delta(\mathbf{x}-\mathbf{y})]\\
&= \frac{\partial F}{\partial \phi(\mathbf{y})} -\nabla\frac{\partial F}{\partial \nabla\phi(\mathbf{x})},
\end{align*}
$$
where the last step was optained using partial integration.

Given our $F_\phi$ we thus find
$$
\begin{align*}
\chi^{-1}(\mathbf{x}, \mathbf{x}')&=\frac{\delta^2 F_\phi}{\delta \phi(\mathbf{x} )\phi(\mathbf{x'})},\\
&=(r + 12 u\langle\phi\rangle^2-c\nabla^2)\delta (\mathbf{x} -\mathbf{x}')\\
\end{align*}
$$

Now we can take the Fourier transform to find (see 2.3.6 from the book)
$$
\begin{align*}
\chi^{-1}(\mathbf{q})&=\frac{1}{V}\int d^dx d^dx'(r + 12 u\langle\phi\rangle^2-c\nabla^2)\delta (\mathbf{x} -\mathbf{x}')e^{i\mathbf{q}(\mathbf{x} -\mathbf{x}')}\\
&=\frac{V}{V} \int d^dx(r + 12 u\langle\phi\rangle^2-c\nabla^2)\delta (\mathbf{x} )e^{i\mathbf{q}\cdot\mathbf{x}}\\
&=r + 12 u\langle\phi\rangle^2+cq^2,
\end{align*}
$$
where in the last step we again used integration by parts. Mean that
$$
\begin{align*}
\chi(\mathbf{q})&=\frac{1}{r + 12 u\langle\phi\rangle+cq^2}\\
&=\frac{\chi}{1+ (a\xi)^2},
\end{align*}
$$
with $\xi = \ c^{1/2}[r + 12 u\langle\phi\rangle^2]^{-1/2}$ the correlation length of the fluctuations and $\chi = \chi(\mathbf{q} = 0) = \frac{1}{r + 12 u\langle\phi\rangle}$. Given this $\chi(\mathbf{q})$, we can write
$$
\begin{align}
\chi(\mathbf{x}, 0) &=\chi\int \frac{d^dq}{(2\pi)^d}
\frac{e^{i\mathbf{q}\mathbf{x}}}{1+ (q\xi)^2}\\
&\overset{2D}{=}\chi\int \frac{dq_xdq_y}{(2\pi)^2}
\frac{e^{i\mathbf{q}\mathbf{x}}}{(|\mathbf{x}|/\xi)^2+ (|\mathbf{q}||\mathbf{x}|)^2}\left(\frac{|\mathbf{x}|}{\xi}\right)^2\\
&\overset{2D}{=}\chi|x|^{-2}\int \frac{dzd\theta}{(2\pi)^2}z
\frac{e^{iz}}{\eta^2+ z^2}\eta^2 \text{ }(z \equiv |x||q|, \mathbf{q}\mathbf{x} = |x||q|\cos\theta, \eta = |\mathbf{x}|/\xi)\\
&\overset{2D}{=}\frac{\chi}{\xi^2}|x|^{-(2-2)}\int \frac{z dz}{(2\pi)^2}\int d\theta
\frac{e^{iz\cos\theta}}{\eta^2+ z^2}\\
&\overset{2D}{=}c^{-1}\int \frac{z dz}{(2\pi)^2}\int d\theta
\frac{e^{iz\cos\theta}}{\eta^2+ z^2}\\
\\
&\overset{3D}{=}\chi\int \frac{dq_xdq_ydq_z}{(2\pi)^3}
\frac{e^{i\mathbf{q}\mathbf{x}}}{(|\mathbf{x}|/\xi)^2+ (|\mathbf{q}||\mathbf{x}|)^2}\left(\frac{|\mathbf{x}|}{\xi}\right)^2\\
&\overset{3D}{=}\chi|x|^{-2}\int \frac{dzd\theta d\phi}{(2\pi)^3}z^2\cos\theta
\frac{e^{iz\cos\theta}}{\eta^2+ z^2}\eta^2 \text{ }(z \equiv |x||q|, \mathbf{q}\mathbf{x} = |x||q|\cos\theta, \eta = |\mathbf{x}|/\xi)\\
&\overset{3D}{=}\frac{\chi}{\xi^2}|x|^{-(3-2)}\int \frac{z^2 dz}{(2\pi)^3}\int d\theta\cos\theta
\frac{e^{iz\cos\theta}}{\eta^2+ z^2}\\
&\overset{3D}{=}c^{-1}|x|^{-1}\int \frac{z^2 dz}{(2\pi)^3}\int d\theta\cos\theta
\frac{e^{iz\cos\theta}}{\eta^2+ z^2}\\
&\overset{3D}{=}c^{-1}|x|^{-1}\frac{1}{\pi}e^{-|x|/\xi}
\end{align}
$$
So this tells us to that the fluctuations order paramter at $\mathbf{x}$ and $\mathbf{0}$ are correlated via the above formula. As we see that $\xi$ is indeed the lengthscale of the correlations.

### Getting the correlations given a finite box
*See https://journals.aps.org/pre/pdf/10.1103/PhysRevE.61.1072*

The $\chi$ that is used in the equations thusfar is the inifinite susceptibility $\chi$^\infty$. However, if we want to measure it from simulations we do not have a infinitely large system. Therefore we want to approximate $\chi$^\infty$ from the finite $\chi$^{L_B}$ given a box with volume $L^2_B$. Let's furthermore assume that we can ignore the fourth order $\langle \phi\rangle$. We know that
$$
\begin{align*}
\chi(\mathbf{x}, \mathbf{x}')=\beta\langle \delta\phi(\mathbf{x})\delta\phi(\mathbf{x'})\rangle
\end{align*}
$$
such that
$$
\begin{align*}
\chi(\mathbf{q} = 0)=\frac{1}{V}\int d^dxd^dx'\beta\langle \delta\phi(\mathbf{x})\delta\phi(\mathbf{x'})\rangle\equiv \chi^\infty
\end{align*}
$$
::: danger
Is there a mistake in the paper? In equation 2 we look at the greens function, however as fare as I understand it we don't consider $\langle \phi_q\phi_{-q}\rangle$, but instead $\langle \delta\phi_q\delta\phi_{-q}\rangle$. Moreover I think the dimensions are not completely right.
:::
::: warning
Question:
- Boundary conditions, and where does the bulk modulation come from?
- How does the Lagrange multplier work?
- Where do the fluctuations come back, the paper sometmies uses fluctuations and sometimes the average values. What are the fluctuations for us? Are these the strain given a certain displacement $u$?
- 2 or 3 independent parameters squares
- Do we measure C or B?
:::

Assume
Assume now we can measure the average deviation correlations in a smaller box.
$$
\begin{align*}
\chi^{L_B} \equiv \frac{1}{L_B^2}\int d^dxd^dx'\beta\langle \delta\phi(\mathbf{x})\delta\phi(\mathbf{x'})\rangle
\end{align*}
$$

### Obtainin Poisson ratio
Given our $B$ and $\mu$ we can find the Poisson ratio by using that
$B = \frac{1}{2}(B_{11} + B_{12})$ and $\mu = \frac{1}{2}(B_{11} - B_{12})$, such that

$$
\nu = \frac{B_{12}}{B_{11}} = \frac{B-\mu}{B+\mu}
$$

# To do

- [ ] Derivation Bulk constant understanding