# Final overview: Poisson ratios polyhedra

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

In this project we measure the Poisson ratio of different hard polyhedra particles. To do so, we measure the elastic constant tensor, which in turn grands us with the Poisson ratio.

## Theory
Expanding the free energy of a crystal in terms of small deformations $\epsilon_{ij}$, we can write

$$
\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}
$$

By first expressing the free energy in terms of $\epsilon$, we can compute these first and second order derivatives. The first order linear response will be associated with the pressure(/stress) tensor. As a result, when computing the second order derivative, we can link the pressure to the strain. If the system initially is under a hydrostatic pressure $P$, such that $\sigma_{\alpha\beta} = -P\delta_{\alpha\beta}$, one can show that
$$
\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**. Because every crystal is associated with many symmetries. First, since the strain tensor is symmetric (strain does not capture rotation, which is the assymetric part of deformation) 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
$$
B_{ijkl}=B_{jilk}=B_{lkji}=B_{klij}
$$
So we can write $B$ as a two-rank tensor $B_{ij} = B_{ji}$ (this is intuitively true because otherwise the potential energy would not be conservative, i.e. applying strain in a different order would change the potential energy). This also follows from the fact that both strain and stress need to be symmetric. This latter follows from the fact that otherwise an infinitely small volume element would experience a torque.

:::info
Strain by definition is symmetric: Given a certain displacement field
the strain is defined as the symmetric part of the tensor, i.e. up to first order in deformation is $\frac{1}{2}(\mathbf{D}^T+\mathbf{D} -I)$, i.e
while rotation is captured in the odd part of the tensor is $\frac{1}{2}(\mathbf{D}^T-\mathbf{D})$, i.e.

The difference between e.g. the first and second shear (where volumechange due to shear is equal in both cases)
is thef fact that the third case is pure shear, while, the second case needs shear and rotation. I.e. the shear matrix in both cases is respectively
while the rotation in both cases is respectively
For shear there is a symmetry line, while for rotation is there is a symmetry axis perpendicular.
:::

### Elastic tensor matrices
*See /home/rinskealkemade/Documents/Elastic constants/Analysis_symmetries/Summetries<>.nb*

Using the crystal symmetries, for squares, triangles and hexagons, this $\mathbf{B}$ matrix is given by
$$
\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 in the case of triangles and hexagons $B_{66}= \frac{B_{11}-B_{12}}{2}$.
For cubes, we find
$$
\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}
$$
and for triangular and hexagonal prims
$$
\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}$

### 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}=V\frac{d^2F}{dV^2}
$$
From our $F$ expansion we however also know that
$$
\begin{align}
\frac{F(\epsilon)}{V}&=\frac{F(0)}{V} + C_{ij}\epsilon_{ij}+ \frac{1}{2}C_{ijkl}\epsilon_{ij}\epsilon_{kl}
\end{align}
$$
in equilibrium we are in a potential minimum, meaning that this first term will fall away, i.e. the free energy density $f = F/V$ is

$$
\begin{align}
f&=f_0 +\frac{1}{2}C_{ijkl}\epsilon_{ij}\epsilon_{kl}
\end{align}
$$
for a uniform expansion in all directions
$$\mathbf{\epsilon} =
\begin{pmatrix}
\frac{1}{3}\delta & 0 \\
0 &\frac{1}{3}\delta &0\\
0 & 0 & \frac{1}{3}\delta\\
\end{pmatrix}
$$
we thus find for cubes $f = \frac{\delta^2}{6}\left[ B_{11} + 2\cdot B_{12}\right]$. At the same time, we can write $V_n = V_o\det(\mathbf{I}+\mathbf{\epsilon} )\approx V_o(1+ \text{Tr}(\mathbf{\epsilon}))+\mathcal{O}(\epsilon^2)$. Therefore $dV= V_0d\text{Tr}(\mathbf{\epsilon})$, which in our case is equal to $dV =V_od\delta$. As a result,
$$
B = V_o\frac{d^2F}{dV^2}=\frac{V_o}{V_o^2}\frac{d^2F}{d\delta^2}=\frac{d^2F/V_o}{d\delta^2} = \frac{d^2\frac{\delta^2}{6}\left[ B_{11} + 2\cdot B_{12}\right]}{d\delta^2} = \frac{1}{3}\left[ B_{11} + 2\cdot B_{12}\right]
$$
Similarily we find for squares, triangles and hexagons $B =\frac{1}{2}(B_{11}+B_{12})$
For triangular and hexagonal prisms, we cannot use the bulk modulus, because in that case the strain tensor under uniform pressure will not have diagonal elements that are identical.

### Measuring the elastic constants
To measure the $\mathbf{B}$'s, we expand the earlier found relation between stress and strain to
$$
\sigma'_{\alpha\beta}(\epsilon) = \sigma_{\alpha\beta} + B_{\alpha\beta\gamma\delta}\cdot\epsilon_{\gamma\delta}
$$
This means that we measure the strain in the system before and after a small deformation. We can obtain the stresses since we know that in equilibrium external stress that is applied f.e. by doing a deformation must be equal to the internal pressure exerted by the system. I.e. we can measure the pressure in the system to obtain minus the stress.

Expressing the strain in terms of some parameter $\delta$, we can up to first order in $\delta$ write
$$
\sigma'_{\alpha\beta}(\delta) = \sigma_{\alpha\beta} + \text{linear combination}_{\gamma\delta}(B_{\alpha\beta\gamma\delta})\cdot\delta + \mathcal{O}(\delta^2)
$$
Assuming that we can expand the measure stress $\sigma'_{\alpha\beta}(\delta)$ also in terms of $\delta$, we find
$$
\sigma'_{\alpha\beta} = \sigma_{\alpha\beta} + A_1\delta + A_2\delta^2 +\mathcal{O}(\delta^3)
$$
Since we are interested in the linear combination of the $B_{\gamma\delta}$'s, we find
$$
\text{linear combination}_{\gamma\delta}(B_{\alpha\beta\gamma\delta}) = A_1
$$
In other words; to find this linear combination of $B_{\gamma\delta}$'s, we have to expand the measured stress of the system in terms of $\delta$ and then obtain the first order.

By doing multiple deformations, we can find different linear relation of $B_{\gamma\delta}$'s, which together with the Bulk modulus give us a system of equations to find all te different $B$ components individually.
### Poisson ratio

*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{\partial\epsilon_{ii}}{\partial\epsilon_{jj}}=- \frac{S_{iijj}\partial\sigma_{jj}}{S_{jjjj}\partial\sigma_{jj}}= - \frac{S_{iijj}}{S_{jjjj}}= - \frac{S_{ij}}{S_{jj}}
$$
where the uniaxial stress is applied in the the $jj$ direction (note that Poisson’s ratio is defined as the ratio of strains that result from applying a uniaxial stress in one direction. In other words, to measure Poisson's ratio, you need to apply stress in only one direction (uniaxial stress) and then observe how the material strains in both the direction of the applied stress and perpendicular to it.

### 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}
$$

## Specific deformations
*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*

### Squares, triangles and hexagons
#### 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},
\text{ such that }
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
-\delta& 0\\
0& \delta\\
\end{pmatrix} + \mathcal{O}(\delta^2)
$$
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
$$
Expanding our measured $\sigma'_{22}$ in terms of $\delta$ grants us with the relation

$$
B_{11}-B_{12} = A_1^\text{ST}
$$
where $A_1$ is the first order in the expansion. Given that $B=\frac{1}{2}(B_{11}+ B_{12})$, we thus find

$$
B_{11} =B+\frac{1}{2}A_1^\text{ST} \text{ and } B_{12} = B-\frac{1}{2}A_1^\text{ST}
$$

#### Measuring $B_{66}$

To measure $B_{66}$, we use the deformation
$$
\mathbf{D}_\text{ShiftXY} =
\begin{pmatrix}
1 & -\delta\\
0 & 1 \\
\end{pmatrix}, \text{ such that }
\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, \text{ i.e. }
B_{66} = -A_1^\text{SH}
$$
Note that for triangles and hexagons we can also use We can find $B_{66}$ using that $B_{66}= \frac{B_{11}-B_{12}}{2}$.

### 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},
\text{ such that }
\mathbf{\epsilon}_\text{StretchXY} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & -\delta & 0\\
0 & 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_{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, \text{ such that }
B_{11}-B_{12} = - A_1
$$
Given that $B=\frac{1}{3}(B_{11}+ 2B_{12})$, we thus find
$$
B_{11} = B-\frac{2}{3}A_1^\text{ST} \text{ and } B_{12} = B+\frac{1}{3}A_1^\text{ST}
$$

#### 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},
\text{ such that }
\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
\text{ i.e. }
B_{66} = -A_1^\text{SH}
$$
### Triangular and hexagonal prisms

#### 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}_\text{StretchX}=
\begin{pmatrix}
1+\delta& 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}
\rightarrow
\mathbf{\epsilon}_\text{StretchX} =
\begin{pmatrix}
\delta& 0 & 0\\
0 & 0 & 0\\
0 & 0 & 0
\end{pmatrix}
$$
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} = A_1^{ST(11)1}\\
B_{12} &=\frac{2P+\sigma'^{ST1}_{11}+\sigma'^{ST12}_{22}}{\delta} = A_1^{ST(11)1} + A_1^{ST(22)12}\\
B_{23} &=\frac{P+\sigma'^{ST(11)1}_{11}+\sigma'^{ST12}_{22}-\sigma'^{ST13}_{22}}{\delta} = A_1^{ST(11)1} + A_1^{ST(22)12}-A_1^{ST(22)13} \\
B_{33} &=\frac{\sigma'^{ST1}_{11}+\sigma'^{ST12}_{22}-\sigma'^{ST13}_{22}-\sigma'^{ST13}_{33}}{\delta} = A_1^{ST(11)1} + A_1^{ST(22)12}-A_1^{ST(22)13}-A_1^{ST(33)13}
\end{align}
$$
where $A^{ST(ii)jk}$ implies that we measure the $ii$ stress component under a stretch deformation of $jk$.

#### 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} \text{ i.e. } B_{55} = -A_1^{SH}
$$
#### Measuring $B_{66}$
To measure $B_{66}$ we can use that
$$B_{66}= \frac{B_{11}-B_{12}}{2}$$

## Measure $A_1$
In order to measure the first order expansion term, we use the following procedure: We measure the stress as a function of $\delta$, through which we then fit a polynomial. Note that in this dataset we also include the stress measured at $\delta = 0$, which we obtain from the pressure measurements. The polynomial we fit has the form
$$
\text{Fit}^{\mathcal{O(\delta^n)}}(\delta) = \sum_{i = 0}^n A_i\delta^i,
$$
where consider $n\in(2,3,4,5)$. We see orders fit the data very well, so we choose to go with the lowest order $n = 2$.

## Details on measurements
For all the 2D systems we measure the stress in the deformed and non-deformed systems in an EDMD simulation for $t/\tau = 5\cdot10^5$. For all the 3D systems we measure the stress in the deformed and non-deformed systems in an EDMD simulation for $t/\tau = 2.5\cdot10^5$. In both cases we measure every $t/\tau =1$. The first quarter of the data is thrown away. Below we discuss what packing fractions and values of $\delta$ we take into account
### Squares
*See notebook Analysis_angle.nb in the square DATA*
For a system of $N=1024$ and packing fractions $\eta \in(0.875, 0.9, 0.91, 0.925, 0.94, 0.95)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008)$
From the analysis (see notebook Analysis_angle_nomean.nb) it becomes clear that we lower packing fractions the particles in the system starts turning because that lowers the free energy. At higher packing fractions or higher $\delta$, the system starts forming wave like configurations.

### Triangles

For a system of $N=1008$ and packing fractions $\rho \in(0.80, 0.81, 0.84, 0.85, 0.86)$ and $\delta=(0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01)$

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. Therefore, we can only use values between $\eta = 0.80$ and $\eta = 0.86$ for our analysis (although we did not check for values between $0.86$ and $0.90$).

### Hexagons
For a system of $N=1024$ and packing fractions $\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)$. The angle data shows that particles wiggle quit a bit around their equilibrium value (which makes sense, since the particles are closer to being round than the cubes and triangles). In the simulation data we see that the shift simulation does not work anymore for $\delta > 0.008$ (I think the shift is too strong here). Therefore we do not include those in the shift analysis.
### Cubes
We measure for a system of $N=1000$ and packing fractions $\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)$.

### Triangular prisms

We measure for a system of $N=1344$ and packing fractions $\rho \in(0.65,0.70,0.725, 0.75,0.775, 0.8,0.825, 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)$.

#### Hexagonal spacking $z$-direction

*See /home/rinskealkemade/Documents/Elastic constants/DATA_triprisms/Analysis_hexspacing.nb*

To measure measure the hexagonal spacing for a a hexagonal crystal as a function of the density, we measure the pressure tensor for various hexagonal spacing and check for which spacing it becomes isotropic.

**Fitting the ratio worked better than fitting the spacing (the spacing also has a denisty component which has to be fit). Howver, the initial code just looks at lattice spacing. Therefore, we ratiolattice.cc code reads in the hexspacing and then outputs the ratio of the lattice constants for that ratio in ratiolatticeconstants.dat.**

| $\sigma_{11}/\sigma_{22}$ | $(\sigma_{11} +\sigma_{22})/2/\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
We also checked that our data matches that of Marjolein

### Hexagonal prisms
For a system of $N=1080$ 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 = 250.000$ timesteps.

#### Hexagonal spacking $z$-direction

*See /home/rinskealkemade/Documents/Elastic constants/DATA_hexprisms/Analysis_hexspacing_ratio.nb*

To measure measure the hexagonal spacing for a a hexagonal crystal as a function of the density, we measure the pressure tensor for various ratio's between the $a$ and $c$ lattic constantand check for which spacing it becomes isotropic.

**Fitting the ratio worked better than fitting the spacing (the spacing also has a denisty component which has to be fit). Howver, the initial code just looks at lattice spacing. Therefore, we ratiolattice.cc code reads in the hexspacing and then outputs the ratio of the lattice constants for that ratio in ratiolatticeconstants.dat.**

This is read in by the notebook, as well as the presure tensor. We found the system is very sensitive to the ratio, so we use a very high fit:
$$ratio = -19727.09908965223 + 206587.652315288 \eta - 945479.602501588 \eta ^2+ 2.4700956319096233\cdot 10^6 \eta ^3 - 4.0290989752363632\cdot 10^6 \eta ^4+ 4.201788014308065\cdot 10^6 \eta ^5 - 2.7358396829915103\cdot 10^6 \eta ^6+ 1.0168488910612391\cdot 10^6 \eta ^7 - 165174.70619932798` \eta ^8$$

| Pressure | Fit | Difference with data |
| -------- | --- | -------- |
||||

Moreover, we found that at low densities the system is quite unstable. Therefore, we only include packing fractions from $\eta=0.70$. Note that the fit picture includes the fit only on the interval we consider.

We also found that the pressure is not isotropic in the $xy$-plane:

This can most likely be attribruted to finite size effects (the fact that the box sizes have different lengths). To test this, we performed simulationes with various boxratios with correct unit cell ratio's at $\eta=0.75$. *See /home/rinskealkemade/Documents/Elastic constants/DATA_hexprisms/pressure/test*

From here it becomes clear that the ratio between pressures in the $x$- and $y$- direction indeed depends on the size of the axis. Below I show the ratio in the pressure in the $x$ and $y$ direction as a function of the ratio of the axis for the correct hexspacing.
In this figure all circular plot markers correspond to systems of around 1000 particles, while the square corresponds to a system of ~8000 particles. Clearly we see that the effect of the axis ratio disappears if the system becomes large enough.