# Overview of project

The primary goals of this Bachelor project will be the analysis of both hard-like cubic systems with the Weeks-Chandler-Andersen potential, as well as the analysis of hard hexagonal prisms. Particularly, the focus will be on discovering the effective elastic constant components $C_{IJ}$.

# Preliminary coding

## Estimating pi

DO LATER

### Using direct sampling

Pi is approximated as: 3.820000, 3.394000, 3.300800, 3.308680

### Using a Markov Chain

#### Multiple runs inside square

Datapoint MSD_hit Rejection Rate
1 0.000000 0.000000
2 0.000005 0.004064
3 0.000004 0.004012
4 0.000006 0.004556
5 0.000006 0.004464
6 0.000007 0.004281
7 0.000004 0.004740
8 0.000007 0.004519
9 0.000004 0.004310
10 0.000005 0.004472

Overall, it seems like line 3, 7 and 9 have the lowest MSD.


#### Singular run inside discrete lattice

This simulation was for a lattice of 10x10 with 1000 distinct simulations and 100000 per simulation.

971763 971236 972329 970530 970191 970234 969653 970084 971313 972123
969277 970043 970517 969369 970009 969319 968566 969861 970223 970314
971914 970309 969578 970083 969942 969300 968401 969468 968626 968454
967811 968525 967739 969531 968201 968289 968640 968748 968527 965075
963913 964969 968085 967736 966806 967220 965867 964880 964987 963223
964831 965456 966778 965749 965558 965583 963189 961789 964682 962986
965690 964926 963743 965141 965730 964892 963805 961540 962938 960842
962752 962553 962111 963290 964124 964584 963642 961070 962820 962992
962936 963700 962228 961092 962041 962525 962077 960541 963499 963959
963026 961975 962588 960282 961580 961842 961440 958407 959791 959952

Highest: 972329, Lowest: 958407
The biggest difference in percentage visited is: 0.013922

### Determining correlation time of discrete lattice as a function of lattice size with transfer matrix

## Random Walk

x y z <r_squared>
74.672821 46.253738 92.260437 110635.812500
74.908638 45.985584 92.682388 110597.289062
75.144928 46.846046 92.721680 110486.859375
74.749260 46.427227 92.129257 110493.859375
75.072433 46.682816 93.091599 110494.437500
74.906555 47.393089 93.537117 110528.875000
74.094864 48.207237 94.469543 110553.765625
73.981979 48.194424 94.355881 110495.609375
73.950813 47.785164 94.170998 110505.148438
74.831345 48.445995 93.283752 110306.187500
74.174370 47.664593 93.405930 110376.773438
73.814537 47.195496 92.683167 110371.328125
73.314171 47.630558 92.246407 110379.671875
73.904236 47.182426 92.483162 110424.101562
72.970703 47.737896 93.258705 110508.515625
73.611420 48.001904 93.792519 110531.148438
73.280006 47.550419 92.969246 110497.546875
72.888283 46.951351 93.514786 110470.109375
72.192734 46.302738 92.912979 110396.625000

This dataset is much larger, but I have not plotted it yet, so just to show that it works I have placed a fraction of the output file here. Additionally, periodic boundary conditions are also implemented. The amount of simulations was 20, amount of steps 10000, box_size 100 and dimensionality 3.

I plotted it :)



# Actual coding

## Monte Carlo simulation of hard-like spheres in NVT ensemble

### Preliminary coding (fillCL.c)

After creating two different programs that place the hard-like spheres in both a random structure (I was really just doodling with code, but I think to be honest it wouldn't be too different from a cubic structure for a huge lattice) and in an ordered structure, I tried to figure out on my own where to find the file named in the exercise. The original link to it was dead, so I used another file linked indirectly by my supervisor, Dr. Filion, named "emptyMC.c". This file was very similar to what I needed to create, but it had a few key differences it seemed. Eventually, I decided to cannibalize this C-program with implementations of my own for a few functions.

The following is the output of the emptyMC.c file:


### Energy

The equation for the interparticle Weeks-Chandler-Andersen (WCA) potential is as follows:

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

The energy for each interparticle potential is equal to one another in the starting configuration (FCC lattice). At a first glance however, there appeared to be very tiny differences between a few energies of these pairs.

Therefore, I have increased the number of decimals of the coordinates to be 18, so that any rounding errors on the double parameters are minimalized. This seemed to yield an energy of 12.000000 for every particle in the box. The unequal energy was thus the result of a rounding error.

After a simulation with parameters cycles=100000, steps=108, Nparticles=108, the following is the output of particle[0] per cycle:


However, for the average energy per particle, calculated as the total energy of the system divided by the amount of particles, we have:


This has the burn-in at the start and the stabilization or equilibrated afterwards, as is usually seen for Monte Carlo simulations.

### Virial pressure

Next, we want to calculate the force between particles to eventually be able to measure the pressure. This force is quite simply just the derivative of the WCA potential with respect to the distance $r$:

$$\begin{equation}
\frac{d\, V_{WCA}}{d\, r} = 4 \varepsilon \left[-12\sigma^{12}r^{-13} + 6 \sigma^{6}r^{-7}\right]
\end{equation}$$

With this expression of the force, we can now calculate the virial pressure of the simulation:

$$
P = \frac{Nk_BT}{V} - \frac{1}{3V} \left< \sum_{j<i}^{N} r_{ij} \cdot \frac{d\, V_{WCA}}{d\, r_{ij}} \right>
$$

Here, the first term is a constant in our simulation. The next term should be calculated per particle pair.

The simulation produced the following plot for the pressure:


This is in line with what I expect. Starting from a perfect FCC crystal, shifting particle positions will increase the pressure from its absolute minimum. However, the pressure seems significantly lower.

After some bugfixing, both the energy and pressure have been adjusted to correspond to the expected values.

#TODO: add images

Next, the volume of the box will be slowly increased to see if the pressure decreases along $1/V^3$. From what I gather at a quick inspection of the curvature of the plot, this seems to indeed be the case. Below is a plot over ten increases of the volume.


The corresponding pressure tensor is

\begin{equation}
P = \frac{Nk_BT}{V}\mathbf{I} - \frac{1}{V} \left< \sum_{j<i}^{N} r_{ij} \otimes \frac{d\, V_{WCA}}{d\, r_{ij}} \right>
\end{equation}

After struggling for (what feels like) weeks, Frank's suggestions have been able to fix, I think, my code regarding calculating the pressure tensor. The following data is with 10000 initial cycles and 100000 cycles.

Box: (4.24, 4.24, 4.24), Npart: 108
Starting total energy of system = 648.000000
78.520984       -0.007989       0.003654
-0.007989       78.514349       0.005759
0.003654        0.005759        78.518754

Compared to Frank's data for an equal density of 1.414:

density                pressure               energy
1.4142099999999975	78.47994952164261	7.44478030521642

The pressure seems to roughly align with his data, differing only by 0.04-0.05.

The next steps are calculating the equations of state and implementing deformation functions for the simulation box.

### Elastic constant components

The elastic constant components are reduced from 81 components to a mere 3 components in Voigt notation (TODO: add text).

For these components we also need to find the bulk modulus $B$:

\begin{equation}
B = \rho \dfrac{\partial P}{\partial \rho}
\end{equation}

This bulk modulus can be obtained from the equation of state. Also, let's look at the appropriate box deformations to calculate the elastic constant components.

#### Box deformations

Below you can find the implementation of ```deformbox()```, which changes particle positions (and box sizes) accordingly for the selected deformation type. The deformation selected for the simulation is a hard-coded variable at the top of the code called ```StretchType``, which is coded to have an integer value of 1 to 3. This integer value corresponds to a type of deformation; in order, these are `ShiftXZ`, `ShiftYZ`, and 'StretchXY'.

void deformbox() 
{
	int particle;

	deltavcur = deltav + (deltavmax-deltav)*scaling/maxsimulations;

	switch (StretchType)
	{
        case 1: // 1 -> ShiftXZ
            for (particle = 0; particle < N; ++particle) {
                r[particle][0] += -deltavcur*box[2]; 
            }
            break;
        case 2: // 2 -> ShiftYZ
            for (particle = 0; particle < N; ++particle) {
                r[particle][1] += -deltavcur*box[2];
            }
            break;
        case 3: // 3 -> StretchXY
            box[0] *= 1+deltavcur;
            box[1] *= 1-deltavcur;
            box[2] *= 1/(1-(deltavcur*deltavcur));

            for (particle = 0; particle < N; ++particle) {
                r[particle][0] *= 1+deltavcur;
                r[particle][1] *= 1-deltavcur;
                r[particle][2] *= 1/(1-(deltavcur*deltavcur));
            }
            break;
        default:
            printf("Invalid StretchType! No deformation selected!");
            break;
	}
}

Initially, only ```StretchType = 3```, or StretchXY, is tested, as the majority of the data we want is stored within this type of deformation.

#### Bulk modulus

As mentioned before, we need to find the equation of state if we want to find the bulk modulus $B$. We then fit the equation of state with $\texttt{SciPy}$ to find the value of $\dfrac{\partial P}{\partial \rho}$.

##### Equation of state

The equation of state seems to show relatively lower values for the pressure than Frank's data suggests. To me this suggests that either our difference in methods to calculate differs, or my code contains yet another bug.

The bug was in double counting measurements in the C-code, but that is now fixed. The equation of state seems to almost perfectly align with Frank's data, as shown below:


### Calculating the deformed stress tensor

Ultimately, the simulation should be able to yield clear values of the elastic constant of the system. The path to these values is through calculating the effective elastic constant $B_{ijkl}$, which relate to the elastic constant $\mathbf{C}$ as follows:

\begin{equation}
B_{ijkl} = C_{ijkl} - P \left( \delta_{ik}\delta_{jl} + \delta_{il} \delta_{jk} - \delta_{ij} \delta_{kl} \right)
\end{equation}

Due to symmetries and rotational symmetries, of the 3x3x3x3 = 81 components fo the effective elastic constant $B_{ijkl}$, only 3 are non-zero. Written in Voigt notation, these components are for ${ijkl} = {11}, {12}, {44}$.

\begin{equation}
\sigma_{ij}^\prime = \sigma_{ij} + \sum_{k,l=1,2,3} B_{ijkl}\epsilon_{kl}
\end{equation}

We want to obtain the stress tensor of the system after the deformations, $\sigma_{ij}^\prime$, next, and for that we need to work with both our obtained stress tensor without these applied deformations and with something not yet before named here: a sum over new indices $k,l$ containing the product of the effective elastic constant $B_{ijkl}$ and the strain tensor $\epsilon_{kl}, which will be discussed later.

For the next step, the box deformations will be tested. As mentioned, we start with the `StretchXY` deformation at `delta = 0.1` to qualify our results for a Taylor expansion.

While the code is running, let's restate Vera's theoretical backgrounds on how to obtain the strain tensor $\epsilon$:

#### Theory behind the calculations

##### Effective elastic constant

##### Strain tensor

The strain tensor $\epsilon$ is commonly writen as

\begin{equation}
\epsilon = \dfrac{1}{2}\left(\dfrac{\partial u_i}{\partial r_j} + \dfrac{\partial u_j}{\partial r_i} + \dfrac{\partial u_k}{\partial r_i}\dfrac{\partial u_k}{\partial r_j}\right)
\end{equation}

and

\begin{equation}
\epsilon = \dfrac{1}{2}\left(\mathbf{D}^T\mathbf{D}-1\right)
\end{equation}

Where $\mathbf{D}$ is the deformation matrix, which will be discussed extensively in a moment. Vera's thesis provided the proof that these two equations of the strain tensor $\epsilon$ are equivalent.

For a `StretchXY` deformation...:

\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\frac{\delta ^2}{2}+\delta & 0 & 0 \\
0 & \frac{\delta ^2}{2}-\delta & 0 \\
0 & 0 & \frac{1}{2 \left(1-\delta ^2\right)^2}-\frac{1}{2} \\
\end{array}
\right)
\end{align}

Because $\delta$ is chosen to be lower than 1, we can throw out the squared deltas and just use

\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\delta & 0 & 0 \\
0 & -\delta & 0 \\
0 & 0 & 0 \\
\end{array}
\right) + \mathcal{O}(\delta^2)
\end{align}

After filling in ..., we obtain

\begin{equation}
\sigma_{ij}^\prime - \sigma_{ij} = B_{ij11}\delta - B_{ij22}\delta + \mathcal{O}(\delta^2)
\end{equation}

Without any deformations, there are no external forces, and therefore $\sigma_{ij} = -P\delta_{ij}$.

For this specific deformation, our only interest will be in $B_{11} = B_{1111}$ and $B_{12} = B_{1122}$, so it follows that the deformed stress tensor will be:

\begin{equation}
\sigma_{ij}^\prime + P\delta_{ij} = B_{11}\delta - B_{22}\delta + \mathcal{O}(\delta^2)
\end{equation}

Now we enter for $ij = 11$

\begin{equation}
\begin{split}
\sigma_{11}^\prime + P &= B_{1111}\delta - B_{1122}\delta + \mathcal{O}(\delta^2)\\
B_{11} + B_{22} &= \dfrac{\sigma_{11}^\prime + P}{\delta} - \mathcal{O}(\delta^2)
\end{split}
\end{equation}

This gives us both $B_{11}$ and $B_{22}$, but not separated. For this we need the bulk modulus $B$

\begin{equation}
B = \dfrac{1}{3}\left(B_{11}+2B_{12}\right) = \dfrac{1}{3}\left(C_{11} + 2C_{12} + P \right)
\end{equation}

This bulk modulus is also defined as the modulus of hydrostatic compression

\begin{equation}
B = -V \dfrac{\mathrm{d} P}{\mathrm{d}V}
\end{equation}

This expression shows that the bulk modulus can be obtained by looking at the rate of change of the equation of state.

Fitting the equation of state yields

a = 2.955952649515459711e+02
b = -6.980109880032289311e+02
c = 5.408608377893353918e+02
d = -1.264562692458837461e+02

Using:

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

This means that at density $\rho = 0.80000$, the corresponding bulk modulus is

\begin{equation}
\begin{split}
B = 3 \cdot 2.955952649515459711 \cdot e^{2} \cdot 1.00000^3 + 2 \cdot -6.980109880032289311 \cdot e^{2} \cdot 1.00000^2 \\
+ 5.408608377893353918 \cdot e^{2} \cdot 1.00000 = 31.6246566375
\end{split}
\end{equation}

## Event-driven molecular simulation

### Hexagonal prisms

Frank and Laura have provided me with the event-driven simulation code needed for calculating the elastic constants of a lattice of hard hexagonal prisms. This code takes into account collisions and proper directions of movement and velocity after collisions.

#### Calculating hexspacingz

For a specific shape of the hexagonal prisms, we quickly encounter a problem: the 2D-shape of the hexagonal is known, but *a priori* the spacing in the z-direction between two instances of such hexagonal prisms is not. Thus, we need to find the correct value for the `hexspacingz` for which the simulation is approximately isotropic, i.e. at its lowest free energy. In other words: at this density the 3D conformation of the hexagonal prism will result in the system that is energetically most favored (most stable). This in turn means that this simulation correlates best to a real-life version of the same unit cell.

We try for different values of `hexspacingz` and run these simulations for 10 minutes each. This gives each simulation time to calibrate to their lowest energy state. We plot the average saved stress tensor components $\sigma_{11} - \sigma_{33}$ and $\sigma_{22} - \sigma_{33}$ as a function of `hexspacingz` and then fit a second-order polynomial function using `scipy.optimize` to find the intersect with the point $\dfrac{\Delta_{1133} + \Delta_{2233}}{2} =0$. This value is the value at which the crystal is isotropic.

Running this on my laptop, for a density of 0.7, the following plot yields that for $y = 0$, the intersect is at `hexspacingz = [0.10823, 0.1828]`. Each point was a simulation that ran for 10 minutes.


Intercept 11-33 is:  0.10827743999887512
Intercept 22-33 is:  0.10822535972878912

There is still a lot of noise on this plot. Therefore, I will switch to using the `Thor` cluster from now on to easily run it for longer amounts of time (8 hours). This way, the datapoints should become less noisy. Also, by simulating over a smaller range of values for `hexspacingz`, the accuracy increases significantly.

By using Brent's method, we can find an estimate of the best value of `hexspacingz` for an isotropic crystal over which to pick a proper range of 20 different values for `hexspacingz`.

After using the cluster for the aforementioned 20 values of `hexspacingz`, we obtained the following plot of `hexspacingz` as a function of density:


For specifically any density


Close to the optimal value of `hexspacingz`, we can see that the stress tensor components are nearly equal


##### Bulk modulus $B$

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

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

Expressing the bulk modulus in terms of equation of state:

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

So, our next step is to calculate the equation of state to find the bulk modulus $B$.

HERE

#### Remaining components of the effective elastic constant $B_{ijkl}$

First off, I found that Vera's thesis may have numbered the elements of the effective elastic constant $B_{ijkl}$ in Voigt notation incorrectly, switching up element 4 (32) with element 6 (21).

For a hexagonal prism there is the same reflection symmetry in the x-axis, as well as the same reflection symmetry in the y-axis. After these, of the 36 components, of which there are 21 independent components (symmetrical), only 9 independent elements will remain. However, beyond that, the hexagonal prisms only have a rotational symmetry around the z-axis of $\dfrac{\pi}{3}$. These will give the following equalities:

#### Shell scripts for automation

The optimal value of `hexspacingz` is dependent on the density. Normally, this would mean that someone would have to eyeball and interpolate the correct value of `hexspacingz` per density. However, this seems very much like a waste of my time, so I have opted to create a shell script that applies Brent's method to find the value where the difference between strain11 and strain 22 with strain33 is closest to 0 (for 5 decimals). It will also loop over a range of densities. Additionally, to optimize the speed, it will run 7 processes concurrently. I believe it may not yet be perfected, but here it is:

#!/bin/bash

# Set the LC_NUMERIC locale to use a period as the decimal separator
LC_NUMERIC="en_US.UTF-8"
command 2>/dev/null

#set -x

# Define the maximum runtime in seconds (10 minutes = 600 seconds)
max_runtime_hexspacingz=60
max_runtime_pressure=600

# Define the start, end, and step values for the sequence
start=0.830000
end=0.850000
step=0.000100

start_hsz=0.02000
end_hsz=0.11000

# Initialize the current value with the start value
current=$start

# Initialize the running process counter
process=0
running_processes=0

tolerance=0.00001

# Set initial values for hexspacingz and avg_strain_difference
hexspacingz=$(bc -l <<< "scale=5; ($start_hsz + $end_hsz) / 2" | awk '{gsub(/,/, "."); printf "%.5f", $0}')
avg_strain_difference=100000  # Some very high number

echo "scale=10" | bc -l

run_for_best_hexspacingz() {
    local_hexspacingz=$1
    local_density=$2
    local_unique_identifier=$3

    local_density=$(bc -l <<< "scale=5; $local_density" | awk '{gsub(/,/, "."); printf "%.5f", $0}')
    local_hexspacingz=$(bc -l <<< "scale=5; $local_hexspacingz" | awk '{gsub(/,/, "."); printf "%.5f", $0}')

    filename="stress.init3.d$local_density.s$local_hexspacingz.x0.0000.y0.0000.z0.0000"

    make -B OUTPUT_FILE="$local_unique_identifier" > /dev/null 2>&1
    timeout "$max_runtime_pressure" ./"$local_unique_identifier" "$local_density" "$local_hexspacingz" > /dev/null 2>&1 &
    pid=&!
    wait $pid
    rm -f "$unique_identifier"
}

# Define the function to calculate avg_strain_difference
calculate_avg_strain_difference() {
    local_hexspacingz=$1
    local_density=$2
    local_unique_identifier=$3

    local_density=$(bc -l <<< "scale=5; $local_density" | awk '{gsub(/,/, "."); printf "%.5f", $0}')
    local_hexspacingz=$(bc -l <<< "scale=5; $local_hexspacingz" | awk '{gsub(/,/, "."); printf "%.5f", $0}')

    filename="stress.init3.d$local_density.s$local_hexspacingz.x0.0000.y0.0000.z0.0000"
    make -B OUTPUT_FILE="$local_unique_identifier" > /dev/null 2>&1
    timeout "$max_runtime_hexspacingz" ./"$local_unique_identifier" "$local_density" "$local_hexspacingz" > /dev/null 2>&1 &
    pid=&!
    wait $pid

    # Update hexspacingz and avg_strain_difference
    avg_strain_difference=$(python latticespacingzverifier.py "$filename" | bc)

    rm -f "$local_unique_identifier"
    echo "$avg_strain_difference"
}

# Function to handle each value of current
run_for_current() {
    local local_current=$1
    local local_hexspacingz=$2
    local process=$3

    local local_start_hsz=$start_hsz
    local local_end_hsz=$end_hsz

    unique_identifier="md_$process"
    echo "Running process for hexspacingz value for density $local_current"

    # Perform Brent's method to find the root
    while (( $(bc <<< "$local_end_hsz - $local_start_hsz > $tolerance") )); do
        # Check if the maximum number of concurrent processes (7 in this case) has been reached

        file_local_density=$(bc -l <<< "scale=5; $local_density" | awk '{gsub(/,/, "."); printf "%.5f", $0}')
        file_local_hexspacingz=$(bc -l <<< "scale=5; $local_hexspacingz" | awk '{gsub(/,/, "."); printf "%.5f", $0}')

        file_filename="stress.init3.d$file_local_density.s$file_local_hexspacingz.x0.0000.y0.0000.z0.0000"

        midpoint=$(bc -l <<< "scale=5; ($local_start_hsz + $local_end_hsz) / 2" | awk '{gsub(/,/, "."); printf "%.5f", $0}')
        avg_strain_difference=$(calculate_avg_strain_difference "$midpoint" "$local_current" "$unique_identifier")
        
        rm -f "$file_filename"

        if (( $(awk -v strain="$avg_strain_difference" -v tol="$tolerance" 'BEGIN { if (+strain <= +tol && +strain >= -tol) print 1; else print 0; }') )); then
            break
        elif (( $(awk -v strain="$avg_strain_difference" -v tol="$tolerance" 'BEGIN { if (+strain > +tol) print 1; else print 0; }') )); then
            echo "Too high value for $local_current: current hexspacingz: $midpoint, current strain difference: $avg_strain_difference"
            local_end_hsz=$midpoint
        else
            echo "Too low value for $local_current: current hexspacingz: $midpoint, current strain difference: $avg_strain_difference"
            local_start_hsz=$midpoint
        fi
        # For bug fixing
    done
    
    echo "Found the root for $local_current: hexspacingz is $midpoint"
    echo "Best found strain difference: $avg_strain_difference"
    echo "Rerunning simulation for longer time"
    $(run_for_best_hexspacingz "$midpoint" "$local_current" "$unique_identifier")
}

# while-loop for the density
while (( $(bc <<< "$current <= $end") )); do
    # Check if the maximum number of concurrent processes (7 in this case) has been reached
    if [ $running_processes -lt 7 ]; then
    (
        # Run the function in a subshell
        run_for_current "$current" "$hexspacingz" "$process"
    )&

    # Increment the current value by the step
    current=$(bc <<< "$current + $step")

    # Increment the running process counter
    ((running_processes++))
    ((process++))

    else
        # Wait for any background process to finish before starting a new one
        wait -n

        # Decrement the running process counter
        ((running_processes--))
    fi
done

# Wait for any remaining background processes to finish
wait

echo "All processes have completed."

In its current form, Brent's method will take at its longest 17 minutes, and with the final simulation taking longer, each datapoint will take at worst 27 minutes. Which means that for 200 datapoints, it will take 5400 minutes or 90 hours. Luckily, as it runs 7 processes simultaneously, this computation time is reduced to a meager 12.86 hours or 12 hours and 52 minutes.

The result is the following:


##### Without deformations

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

$B = \frac{2}{9}\left(C_{11} + C_{12} + 2C_{13} + \frac{1}{2}C_{33}\right)$

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

Expressing the bulk modulus in terms of equation of state:

$B = \frac{2}{9}\left(C_{11} + C_{12} + 2C_{13} + \frac{1}{2}C_{33}\right) = \rho \frac{dP}{d\rho}$

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

#### Reduction of the effective elastic components

When applying any symmetry operations on the hexagonal prism, you would expect that the strain components also change accordingly. More specifically, by applying symmetry operations, one can argue that some components are equal to other components and even that some components are, in fact, zero. This means we can reduce the amount of components from 6x6 to potentially just a few. First, let's state the equation for the transformed effective elastic components $\bar{C}_{ijkl}$ after any appropriate symmetry operation:

$\bar{C}_{ijkl} = \dfrac{\partial \bar{x}_i}{\partial x_m} \dfrac{\partial \bar{x}_j}{\partial x_n} \dfrac{\partial \bar{x}_k}{\partial x_p} \dfrac{\partial \bar{x}_l}{\partial x_q} C_{mnpq}$

Note the use of Einstein summation used in this equation.
For the most part, we should be able to follow Vera's initials symmetry operations as they are also present for the hexagonal prisms. This means that, after the x-reflection and y-reflection, we are left with only those components that arent 1jkl, 111l, 2jkl, 222l.

\begin{pmatrix}
C_{1111} & C_{1122} & C_{1133} & 0 & 0 & 0 \\
C_{1122} & C_{2222} & C_{2233} & 0 & 0 & 0 \\
C_{1133} & C_{2233} & C_{3333} & 0 & 0 & 0 \\
0 & 0 & 0 & 2C_{2323} & 0 & 0 \\
0 & 0 & 0 & 0 & 2C_{1313} & 0 \\
0 & 0 & 0 & 0 & 0 & 2C_{1212}
\end{pmatrix}

or

\begin{pmatrix}
C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\
C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\
C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\
0 & 0 & 0 & 2C_{44} & 0 & 0 \\
0 & 0 & 0 & 0 & 2C_{55} & 0 \\
0 & 0 & 0 & 0 & 0 & 2C_{66}
\end{pmatrix}

The hexagonal prisms also have a reflection in the z-axis. Thusly, we know that $C_{3jkl} = 0$ when $j,k,l \neq 3$, but also $C_{333l} = 0$ when $l \neq 3$. This does not lead to any further reductions to the amount of individual nonzero components, but it is named at any rate to show all possible symmetry operations are applied.

The last remaining symmetry operation is a rotation around the z-axis by $\dfrac{\pi}{3}$. This rotation has the following rotation operator

\begin{equation}
\mathbf{R} = \begin{bmatrix}
\cos(\theta) & -\sin(\theta) & 0 \\
\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\end{equation}

where $\theta = \dfrac{\pi}{3}$, of course. the elastic constant tensor in Voigt notation is now:

\begin{pmatrix}
C_{1111} & C_{1122} & C_{1133} & 0 & 0 & 0 \\
C_{1122} & C_{2222} & C_{2233} & 0 & 0 & 0 \\
C_{1133} & C_{2233} & C_{3333} & 0 & 0 & 0 \\
0 & 0 & 0 & C_{2323} & 0 & 0 \\
0 & 0 & 0 & 0 & C_{1313} & 0 \\
0 & 0 & 0 & 0 & 0 & C_{1212}
\end{pmatrix}

Let's work through all 9 of them, so I don't have to touch Mathematica.

##### 1. $\bar{C}_{1111}$

From a rotation around the z-axis by $\dfrac{\pi}{6}$, we find the following equalities (using Einstein summation):

\begin{multline}
\bar{C}_{1111} = \cos^{4}(\theta) \cdot C_{1111} - \cos^{3}(\theta)\sin(\theta) \cdot C_{1112} - \cos^{3}(\theta)\sin(\theta) \cdot C_{1121} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{1122} - \cos^{3}(\theta)\sin(\theta) \cdot C_{1211} + \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{1212} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{1221} - \cos(\theta)\sin^{3}(\theta) \cdot C_{1222} - \cos^{3}(\theta)\sin(\theta) \cdot C_{2111} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{2112} + \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{2121} - \cos(\theta)\sin^{3}(\theta) \cdot C_{2122} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{2211} - \cos(\theta)\sin^{3}(\theta) \cdot C_{2212} - \cos(\theta)\sin^{3}(\theta) \cdot C_{2221} \\
+ \sin^{4}(\theta) \cdot C_{2222}
\end{multline}

In voigt notation, this is equal to

\begin{multline}
\bar{C}_{11} = \cos^{4}(\theta) \cdot C_{11} - \cos^{3}(\theta)\sin(\theta) \cdot C_{16} - \cos^{3}(\theta)\sin(\theta) \cdot C_{16} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{12} - \cos^{3}(\theta)\sin(\theta) \cdot C_{61} + \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{66} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{66} - \cos(\theta)\sin^{3}(\theta) \cdot C_{62} - \cos^{3}(\theta)\sin(\theta) \cdot C_{61} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{66} + \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{66} - \cos(\theta)\sin^{3}(\theta) \cdot C_{62} \\
+ \cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{21} - \cos(\theta)\sin^{3}(\theta) \cdot C_{26} - \cos(\theta)\sin^{3}(\theta) \cdot C_{26} \\
+ \sin^{4}(\theta) \cdot C_{22}
\end{multline}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{multline}
\bar{C}_{11} = \cos^{4}(\theta) \cdot C_{11} - 4\cos^{3}(\theta)\sin(\theta) \cdot C_{16} + 2\cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{12} \\
+ 4\cos^{2}(\theta)\sin^{2}(\theta) \cdot C_{66} - 4\cos(\theta)\sin^{3}(\theta) \cdot C_{62} + \sin^{4}(\theta) \cdot C_{22}
\end{multline}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{11} = \dfrac{1}{16} C_{11} - \dfrac{4\sqrt{3}}{16}C_{16} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} - \dfrac{12\sqrt{3}}{16} C_{62} + \dfrac{9}{16} C_{22}
\end{equation}

then:

\begin{equation}
\bar{C}_{11} = \dfrac{1}{16} C_{11} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} + \dfrac{9}{16} C_{22}
\end{equation}

> [name=Frank] Agree

##### 2. $\bar{C}_{2222}$

Likewise, for $\bar{C}_{2222}$:

\begin{multline}
\bar{C}_{2222} = \sin^{4}(\theta) \cdot C_{1111} + \sin^{3}(\theta)\cos(\theta) \cdot C_{1112} + \sin^{3}(\theta)\cos(\theta) \cdot C_{1121} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1122} + \sin^{3}(\theta)\cos(\theta) \cdot C_{1211} + \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1212} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1221} + \sin(\theta)\cos^{3}(\theta) \cdot C_{1222} + \sin^{3}(\theta)\cos(\theta) \cdot C_{2111} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2112} + \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2121} + \sin(\theta)\cos^{3}(\theta) \cdot C_{2122} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2211} + \sin(\theta)\cos^{3}(\theta) \cdot C_{2212} + \sin(\theta)\cos^{3}(\theta) \cdot C_{2221} \\
+ \cos^{4}(\theta) \cdot C_{2222}
\end{multline}

In voigt notation, this is equal to

\begin{multline}
\bar{C}_{22} = \sin^{4}(\theta) \cdot C_{11} + \sin^{3}(\theta)\cos(\theta) \cdot C_{16} + \sin^{3}(\theta)\cos(\theta) \cdot C_{16} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{12} + \sin^{3}(\theta)\cos(\theta) \cdot C_{61} + \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin(\theta)\cos^{3}(\theta) \cdot C_{62} + \sin^{3}(\theta)\cos(\theta) \cdot C_{61} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin(\theta)\cos^{3}(\theta) \cdot C_{62} \\
+ \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{21} + \sin(\theta)\cos^{3}(\theta) \cdot C_{26} + \sin(\theta)\cos^{3}(\theta) \cdot C_{26} \\
+ \cos^{4}(\theta) \cdot C_{22}
\end{multline}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{multline}
\bar{C}_{22} = \sin^{4}(\theta) \cdot C_{11} + 4\sin^{3}(\theta)\cos(\theta) \cdot C_{16} + 2\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{12} \\
+ 4\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + 4\sin(\theta)\cos^{3}(\theta) \cdot C_{62} + \cos^{4}(\theta) \cdot C_{22}
\end{multline}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{22} = \dfrac{9}{16} C_{11} + \dfrac{12\sqrt{3}}{16}C_{16} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} + \dfrac{4\sqrt{3}}{16} C_{62} + \dfrac{1}{16} C_{22}
\end{equation}

then:

\begin{equation}
\bar{C}_{22} = \dfrac{9}{16} C_{11} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} + \dfrac{1}{16} C_{22}
\end{equation}
> [name=Frank] Agree
##### 3. $\bar{C}_{3333}$

Likewise, for $\bar{C}_{3333}$:

\begin{equation}
\bar{C}_{3333} = \left(1\right) \left(1\right) \left(1\right) \left(1\right) C_{3333}
\end{equation}
In voigt notation, this is equal to

\begin{equation}
\bar{C}_{33} = C_{33}
\end{equation}
> [name=Frank] Agree
##### 4. $\bar{C}_{1122}$

Likewise, for $\bar{C}_{1122}$:

\begin{multline}
\bar{C}_{1122} = \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1111} + - \sin(\theta)\cos^{3}(\theta) \cdot C_{1112} + -\sin(\theta)\cos^{3}(\theta) \cdot C_{1121} + \\
\cos^{4}(\theta) \cdot C_{1122} + \sin^{3}(\theta)\cos(\theta) \cdot C_{1211} + - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1212} + \\
-\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1221} + \sin(\theta)\cos^{3}(\theta) \cdot C_{1222} + \sin^{3}(\theta)\cos(\theta) \cdot C_{2111} + \\
-\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2112} + - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2121} + \sin(\theta)\cos^{3}(\theta) \cdot C_{2122} + \\
\sin^{4}(\theta) \cdot C_{2211} + -\sin^{3}(\theta)\cos(\theta) \cdot C_{2212} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{2221} + \\
\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2222}
\end{multline}

In voigt notation, this is equal to

\begin{multline}
\bar{C}_{12} = \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{11} + - \sin(\theta)\cos^{3}(\theta) \cdot C_{16} + -\sin(\theta)\cos^{3}(\theta) \cdot C_{16} + \\
\cos^{4}(\theta) \cdot C_{12} + \sin^{3}(\theta)\cos(\theta) \cdot C_{61} + - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \\
-\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin(\theta)\cos^{3}(\theta) \cdot C_{62} + \sin^{3}(\theta)\cos(\theta) \cdot C_{61} + \\
-\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin(\theta)\cos^{3}(\theta) \cdot C_{62} + \\
\sin^{4}(\theta) \cdot C_{21} + -\sin^{3}(\theta)\cos(\theta) \cdot C_{26} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{26} + \\
\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{22}
\end{multline}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{multline}
\bar{C}_{12} = \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{11} + 2\left( \sin^{3}(\theta)\cos(\theta) - \sin(\theta)\cos^{3}(\theta)\right) \cdot C_{16} +\\
\left(\sin^{4}(\theta) + \cos^{4}(\theta)\right) \cdot C_{12} - 4\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \\
2\left(\sin(\theta)\cos^{3}(\theta) - \sin^{3}(\theta)\cos(\theta)\right) \cdot C_{62} +\\
\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{22}
\end{multline}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{12} = \dfrac{3}{16} C_{11} + \dfrac{4\sqrt{3}}{16} C_{16} + \dfrac{10}{16} C_{12} - \dfrac{12}{16} C_{66} - \dfrac{4\sqrt{3}}{16} C_{62} + \dfrac{3}{16} C_{22}
\end{equation}

then:

\begin{equation}
\bar{C}_{12} = \dfrac{3}{16} C_{11} + \dfrac{10}{16} C_{12} - \dfrac{12}{16} C_{66} + \dfrac{3}{16} C_{22}
\end{equation}
> [name=Frank] Agree
##### 5. $\bar{C}_{1133}$

Likewise, for $\bar{C}_{1133}$:

\begin{equation}
\bar{C}_{1133} = \cos^{2}(\theta) \cdot C_{1133} + \sin(\theta)\cos(\theta) \cdot C_{1233} + \sin(\theta)\cos(\theta) \cdot C_{2133} + \sin^{2}(\theta) \cdot C_{2233}
\end{equation}

In voigt notation, this is equal to

\begin{equation}
\bar{C}_{13} = \cos^{2}(\theta) \cdot C_{13} + \sin(\theta)\cos(\theta) \cdot C_{63} + \sin(\theta)\cos(\theta) \cdot C_{63} + \sin^{2}(\theta) \cdot C_{23}
\end{equation}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{equation}
\bar{C}_{13} = \cos^{2}(\theta) \cdot C_{13} + 2 \sin(\theta)\cos(\theta) \cdot C_{63} + \sin^{2}(\theta) \cdot C_{23}
\end{equation}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{13} = \dfrac{1}{4} C_{13} + \dfrac{2\sqrt{3}}{4} C_{63} + \dfrac{3}{4} C_{23}
\end{equation}

then:

\begin{equation}
\bar{C}_{13} = \dfrac{1}{4} C_{13} + \dfrac{3}{4} C_{23}
\end{equation}
> [name=Frank] Agree
##### 6. $\bar{C}_{2233}$

Likewise, for $\bar{C}_{2233}$:

\begin{equation}
\bar{C}_{2233} = \sin^{2}(\theta) \cdot C_{1133} + - \sin(\theta)\cos(\theta) \cdot C_{1233} + - \sin(\theta)\cos(\theta) \cdot C_{2133} + \cos^{2}(\theta) \cdot C_{2233}
\end{equation}

In voigt notation, this is equal to

\begin{equation}
\bar{C}_{23} = \sin^{2}(\theta) \cdot C_{13} - \sin(\theta)\cos(\theta) \cdot C_{63} - \sin(\theta)\cos(\theta) \cdot C_{63} + \cos^{2}(\theta) \cdot C_{23}
\end{equation}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{equation}
\bar{C}_{23} = \sin^{2}(\theta) \cdot C_{13} - 2\sin(\theta)\cos(\theta) \cdot C_{63} + \cos^{2}(\theta) \cdot C_{23}
\end{equation}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{23} = \dfrac{3}{4} C_{13} - \dfrac{2\sqrt{3}}{4} C_{63} + \dfrac{1}{4} C_{23}
\end{equation}

then:

\begin{equation}
\bar{C}_{23} = \dfrac{3}{4} C_{13} + \dfrac{1}{4} C_{23}
\end{equation}
> [name=Frank] Agree
##### 7. $\bar{C}_{1212}$

Likewise, for $\bar{C}_{1212}$:

\begin{multline}
\bar{C}_{1212} = \sin^{2}\cos^{2}(\theta) \cdot C_{1111} + -\sin(\theta)\cos^{3}(\theta) \cdot C_{1112} + \sin^{3}(\theta)\cos(\theta) \cdot C_{1121} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1122} + - \sin(\theta)\cos^{3}(\theta) \cdot C_{1211} + \cos^{4}(\theta) \cdot C_{1212} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{1221} + \sin(\theta)\cos^{3}(\theta) \cdot C_{1222} + \sin^{3}(\theta)\cos(\theta) \cdot C_{2111} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2112} + \sin^{4}(\theta) \cdot C_{2121} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{2122} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{2211} + \sin(\theta)\cos^{3}(\theta) \cdot C_{2212} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{2221} \\
+ \sin^{2}\cos^{2}(\theta) \cdot C_{2222}
\end{multline}

In voigt notation, this is equal to

\begin{multline}
\bar{C}_{66} = \sin^{2}\cos^{2}(\theta) \cdot C_{11} + -\sin(\theta)\cos^{3}(\theta) \cdot C_{16} + \sin^{3}(\theta)\cos(\theta) \cdot C_{16} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{12} + - \sin(\theta)\cos^{3}(\theta) \cdot C_{61} + \cos^{4}(\theta) \cdot C_{66} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin(\theta)\cos^{3}(\theta) \cdot C_{62} + \sin^{3}(\theta)\cos(\theta) \cdot C_{61} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{66} + \sin^{4}(\theta) \cdot C_{66} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{62} \\
+ - \sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{21} + \sin(\theta)\cos^{3}(\theta) \cdot C_{26} + - \sin^{3}(\theta)\cos(\theta) \cdot C_{26} \\
+ \sin^{2}\cos^{2}(\theta) \cdot C_{22}
\end{multline}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{multline}
\bar{C}_{66} = \sin^{2}\cos^{2}(\theta) \cdot C_{11} + 2 \left(\sin^{3}(\theta)\cos(\theta)-\sin(\theta)\cos^{3}(\theta)\right) \cdot C_{16} \\
- 2\sin^{2}(\theta)\cos^{2}(\theta) \cdot C_{12} + \left(\sin^{2}(\theta) - \cos^{2}(\theta)\right)^{2} \cdot C_{66} \\
+ 2\left(\sin(\theta)\cos^{3}(\theta) - \sin^{3}(\theta)\cos(\theta)\right) \cdot C_{62} + \sin^{2}\cos^{2}(\theta) \cdot C_{22}
\end{multline}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{66} = \dfrac{3}{16} C_{11} + \dfrac{4\sqrt{3}}{16}C_{16} - \dfrac{6}{16} C_{12} + \dfrac{4}{16} C_{66} - \dfrac{4\sqrt{3}}{16} C_{62} + \dfrac{3}{16} C_{22}
\end{equation}

then:

\begin{equation}
\bar{C}_{66} = \dfrac{3}{16} C_{11} - \dfrac{6}{16} C_{12} + \dfrac{4}{16} C_{66} + \dfrac{3}{16} C_{22}
\end{equation}
> [name=Frank] I think the $C_{66}$ prefactor should be 4/16. This also seems consistent with your equation before filling in $\pi/3$.

##### 8. $\bar{C}_{1313}$

Likewise, for $\bar{C}_{1313}$:

\begin{equation}
\bar{C}_{1313} = \cos^{2}(\theta) \cdot C_{1313} + \sin(\theta)\cos(\theta) \cdot C_{1323} + \sin(\theta)\cos(\theta) \cdot C_{2313} + \sin^{2}(\theta) \cdot C_{2323}
\end{equation}

In voigt notation, this is equal to

\begin{equation}
\bar{C}_{55} = \cos^{2}(\theta) \cdot C_{55} + \sin(\theta)\cos(\theta) \cdot C_{54} + \sin(\theta)\cos(\theta) \cdot C_{45} + \sin^{2}(\theta) \cdot C_{44}
\end{equation}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{equation}
\bar{C}_{55} = \cos^{2}(\theta) \cdot C_{55} + 2\sin(\theta)\cos(\theta) \cdot C_{54} + \sin^{2}(\theta) \cdot C_{44}
\end{equation}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{55} = \dfrac{1}{4} C_{55} + \dfrac{2\sqrt{3}}{4} C_{54} + \dfrac{3}{4} C_{44}
\end{equation}

then:

\begin{equation}
\bar{C}_{55} = \dfrac{1}{4} C_{55} + \dfrac{3}{4} C_{44}
\end{equation}

> [name=Frank] Agree
##### 9. $\bar{C}_{2323}$

Likewise, for $\bar{C}_{2323}$:

\begin{equation}
\bar{C}_{2323} = \sin^{2}(\theta) \cdot C_{1313} - \sin(\theta)\cos(\theta) \cdot C_{1323} - \sin(\theta)\cos(\theta) \cdot C_{2313} + \cos^{2}(\theta) \cdot C_{2323}
\end{equation}

In voigt notation, this is equal to

\begin{equation}
\bar{C}_{44} = \sin^{2}(\theta) \cdot C_{55} - \sin(\theta)\cos(\theta) \cdot C_{54} - \sin(\theta)\cos(\theta) \cdot C_{45} + \cos^{2}(\theta) \cdot C_{44}
\end{equation}

Ordering this, and taking that the matrix is symmetrical, yields

\begin{equation}
\bar{C}_{44} = \sin^{2}(\theta) \cdot C_{55} - 2 \sin(\theta)\cos(\theta) \cdot C_{54} + \cos^{2}(\theta) \cdot C_{44}
\end{equation}

Filling in $\theta=\dfrac{\pi}{3}$ yields

\begin{equation}
\bar{C}_{44} = \dfrac{3}{4} C_{55} - \dfrac{2\sqrt{3}}{4} C_{54} + \dfrac{1}{4} C_{44}
\end{equation}

then:

\begin{equation}
\bar{C}_{44} = \dfrac{3}{4} C_{55} + \dfrac{1}{4} C_{44}
\end{equation}
> [name=Frank] Agree

#### All equations

In conclusion, we find the following set of equations:
\begin{equation}
\begin{cases}
\bar{C}_{11} = \dfrac{1}{16} C_{11} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} + \dfrac{9}{16} C_{22} \\
\bar{C}_{22} = \dfrac{9}{16} C_{11} + \dfrac{6}{16} C_{12} + \dfrac{12}{16} C_{66} + \dfrac{1}{16} C_{22} \\
\bar{C}_{33} = C_{33} \\
\bar{C}_{12} = \dfrac{3}{16} C_{11} + \dfrac{10}{16} C_{12} - \dfrac{12}{16} C_{66} + \dfrac{3}{16} C_{22} \\
\bar{C}_{13} = \dfrac{1}{4} C_{13} + \dfrac{3}{4} C_{23} \\
\bar{C}_{23} = \dfrac{3}{4} C_{13} + \dfrac{1}{4} C_{23} \\
\bar{C}_{66} = \dfrac{3}{16} C_{11} - \dfrac{6}{16} C_{12} + \dfrac{4}{16} C_{66} + \dfrac{3}{16} C_{22} \\
\bar{C}_{55} = \dfrac{1}{4} C_{55} + \dfrac{3}{4} C_{44} \\
\bar{C}_{44} = \dfrac{3}{4} C_{55} + \dfrac{1}{4} C_{44}
\end{cases}
\end{equation}

>[name=Frank] Okay, so the above agrees with my math. From here, and imposing $\bar{C} = C$, we can trivially obtain the following equations:
\begin{equation}
\begin{split}
C_{13} = C_{23}\\
C_{44} = C_{55}
\end{split}
\end{equation}
by looking at the equations with only two terms on the right-hand side.

>[name=Frank]Additionally, subtracting the first two equations, we obtain:
\begin{equation}
C_{11}-C_{22} = -\dfrac{1}{2}C_{11} + \dfrac{1}{2}C_{22},
\end{equation}
which gives us $C_{11} = C_{22}$ as you also had before.

>[name=Frank] Then, using that last equality, the equation for $C_{66}$ can be rewritten as:
\begin{equation}
\dfrac{3}{4}C_{66} = \dfrac{3}{8}C_{11} - \dfrac{3}{8}C_{12},
\end{equation}
which can be simplified to $C_{66} = (C_{11} - C_{12})/2$.
(Note that all other equations are redundant: if you substitute in the results above, they should just reduce to being always true.)

This leads to a few distinct linear system, but let's look at the following:

\begin{bmatrix}
\dfrac{1}{16} & \dfrac{6}{16} & \dfrac{9}{16} & \dfrac{12}{16} \\
\dfrac{3}{16} & \dfrac{10}{16} & \dfrac{3}{16} & -\dfrac{12}{16} \\
\dfrac{9}{16} & \dfrac{6}{16} & \dfrac{1}{16} & \dfrac{12}{16} \\
\dfrac{3}{16} & -\dfrac{6}{16} & \dfrac{3}{16} & \dfrac{4}{16}
\end{bmatrix}
\begin{bmatrix}
C_{11} \\
C_{12} \\
C_{22} \\
C_{66}
\end{bmatrix}
=
\begin{bmatrix}
\bar{C}_{11} \\
\bar{C}_{12} \\
\bar{C}_{22} \\
\bar{C}_{66}
\end{bmatrix}

Don't look at what is below here

\begin{equation}
\begin{split}
C_{1111} = C_{2222}\\
C_{1122} = C_{2211}\\
C_{1133} = C_{2233}\\
C_{2323} = C_{1313}\\
\end{split}
\end{equation}

This further reduces our elastic strain tensor to the following:

\begin{pmatrix}
C_{1111} & C_{1122} & C_{1133} & 0 & 0 & 0 \\
C_{1122} & C_{1111} & C_{1133} & 0 & 0 & 0 \\
C_{1133} & C_{1133} & C_{3333} & 0 & 0 & 0 \\
0 & 0 & 0 & C_{2323} & 0 & 0 \\
0 & 0 & 0 & 0 & C_{2323} & 0 \\
0 & 0 & 0 & 0 & 0 & C_{1212}
\end{pmatrix}

Thus, of the 9 remaining independent components, after these equalities, only 6 will remain. These are $C_{11}$, $C_{12}$, $C_{13}$, $C_{33}$, $C_{44}$ and $C_{66}$.

### Bulk modulus

Equation 43 of Kittel gives us an equation for the energy density of the crystal:

\begin{multline}
U = \dfrac{1}{2} \left( C_{11} \left( e_{xx}^{2} + e_{yy}^{2} \right) + C_{12} \left( e_{xx}e_{yy} + e_{yy}e_{xx} \right) + C_{13} \left( e_{xx}e_{zz} + e_{zz}e_{xx} + e_{yy}e_{zz} + e_{zz}e_{yy} \right) + \right.\\
\left. C_{33} \left( e_{zz}^{2} \right) + C_{44} \left( e_{yz}^{2} + e_{xz}^{2} \right) + \left( \dfrac{C_{11} - C_{12}}{2} \right) \left( e_{xy}^{2} \right) \right)
\end{multline}

Which can be rewritten in a simpler form

\begin{multline}
U = \dfrac{1}{2} \left( C_{11} \left( e_{xx}^{2} + e_{yy}^{2} + \dfrac{1}{2} e_{xy}^{2} \right) + C_{12} \left( 2e_{xx}e_{yy} - \dfrac{1}{2} e_{xy}^{2} \right) + C_{13} \left( 2e_{xx}e_{zz} + 2e_{yy}e_{zz} \right) + \right.\\
\left. C_{33} \left( e_{zz}^{2} \right) + C_{44} \left( e_{yz}^{2} + e_{xz}^{2} \right) \right)
\end{multline}

We consider now the uniform dilation $e_{xx} = e_{yy} = e_{zz} = \frac{1}{3}\delta$. For this deformation, the accompanying energy density of the hexagonal prisms crystal will be

\begin{equation}
U = \dfrac{1}{18} \left( 2 C_{11} + 2 C_{12} + 4 C_{13} + C_{33} \right) \delta^{2}
\end{equation}

From the definition of the bulk modulus $U = \frac{1}{2} B \delta^{2}$, we can quickly find that

\begin{equation}
B = \dfrac{1}{9} \left( 2 C_{11} + 2 C_{12} + 4 C_{13} + C_{33} \right)
\end{equation}

#### Something

I have found the equation of state. I still have to actually remove the burn-in (i'm sticking with that word). Also, I notice that some hexspacingz values are still a bit off:



#### Elastic constant components derivations

The goal is to derive proper deformations for the hexagonal prisms, that are able to account for all remaining independent components of $B_{ijkl}$. Let us start with the components as defined by Belde.

Equation 2.53:

\begin{equation}
\dfrac{\partial \sigma_{ij}}{\partial \epsilon_{kl}}{\Big |}_{\epsilon = 0} = C_{ijkl} + \left(\delta_{ij}\delta_{kl} - \delta_{jl}\delta_{ik} - \delta_{il}\delta_{jk}\right)P = B_{ijkl}
\end{equation}

This leads to equation 2.61

\begin{equation}
\sigma^\prime_{ij} = \sigma_{ij} + B_{ijkl}\epsilon_{kl}
\end{equation}

This is the backbone to the derivation of every elastic constant component.

##### StretchXY

Equation 2.64

\begin{equation}
\sigma^\prime_{ij} - \sigma_{ij} = B_{ij11}\delta - B_{ij22}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

Entering $ij = 11, 22, 33$ gives

\begin{equation}
\sigma^\prime_{11} - \sigma_{11} = B_{1111}\delta - B_{1122}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

\begin{equation}
\sigma^\prime_{22} - \sigma_{22} = B_{2211}\delta - B_{2222}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

\begin{equation}
\sigma^\prime_{33} - \sigma_{33} = B_{3311}\delta - B_{3322}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

This gives us expressions for $B_{11}$, $B_{12}$, $B_{22}$, $B_{31}$ and $B_{32}$ all with one deformation.

> [name=Frank] Kind of... but note that according to your derivation above, B11=B22 anyway, meaning the first two lines essentially give you the same equality. Also, B31=B32 from the same derivation, which means that in the last line the linear term vanishes.
> It may be necessary to also to a stretchXZ to get a handle on B33.

##### StretchXZ

Similarly to StretchXY, we define a new type of deformation not named in Vera's thesis: StretchXZ. Like StretchXY, this deformation changes positively in the x-direction, while now also shrinking equally in the z-direction. This deformation is defined as:

\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\frac{\delta ^2}{2}+\delta & 0 & 0 \\
0 & \frac{1}{2 \left(1-\delta ^2\right)^2}-\frac{1}{2} & 0 \\
0 & 0 & \frac{\delta ^2}{2}-\delta \\
\end{array}
\right)
\end{align}

Because $\delta$ is chosen to be lower than 1, we can throw out the squared deltas and just use

\begin{align}
\mathbf{\epsilon} = \left(
\begin{array}{ccc}
\delta & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & -\delta \\
\end{array}
\right) + \mathcal{O}(\delta^2)
\end{align}

Using equation 2.64 yet again yields

\begin{equation}
\sigma^\prime_{ij} - \sigma_{ij} = B_{ij11}\delta - B_{ij33}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

Entering $ij = 33$ gives

\begin{equation}
\sigma^\prime_{33} - \sigma_{33} = B_{3311}\delta - B_{3333}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

This gives us expressions for $B_{33}$.

Also, entering $ij = 11$ gives

\begin{equation}
\sigma^\prime_{11} - \sigma_{11} = B_{1111}\delta - B_{1133}\delta + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

##### ShiftXZ

\begin{equation}
\sigma^\prime_{13} = -B_{1313}\dfrac{\delta}{2} -B_{1331}\dfrac{\delta}{2} = - B_{55}\delta
\end{equation}

##### ShiftYZ

\begin{equation}
\sigma^\prime_{23} = -B_{2332}\dfrac{\delta}{2} -B_{2323}\dfrac{\delta}{2} = - B_{44}\delta
\end{equation}

##### ShiftXY

We only need to find $B_{66}$, as the last remaining component. From `ShiftXZ` and `ShiftYZ` it becomes clear a deformation `ShiftXY` fits the bill:

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

According to $\epsilon = \dfrac{1}{2}\left(\mathbf{D}^{T}\mathbf{D} - \mathbf{I}\right)$, the corresponding strain will be

\begin{equation}
\epsilon =
\begin{pmatrix}
0 & -\frac{\delta}{2} & 0 \\
-\frac{\delta}{2} & 0 & 0 \\
0 & 0 & \frac{\delta^{2}}{2}
\end{pmatrix} =
\begin{pmatrix}
0 & -\frac{\delta}{2} & 0 \\
-\frac{\delta}{2} & 0 & 0 \\
0 & 0 & 0
\end{pmatrix} + \mathcal{O}\left(\delta^{2}\right)
\end{equation}

\begin{equation}
\sigma^\prime_{12} = -B_{1221}\dfrac{\delta}{2} -B_{1212}\dfrac{\delta}{2} = - B_{66}\delta
\end{equation}

We may not even need it, but it is good to have defined anyways.
>[name=Frank] Yes, this would be an option, but is a bit more tedious to implement than another stretch. Since we can simply get $B_{66}$ from its relation to $B_{12}$ and $B_{11}$, we do indeed not need this.

#### Expressions of the elastic constant components

From the defined deformations and the expression of the bulk modulus, we should now be capable of deriving the expressions leading to the values of the elastic constant components.

The bulk modulus must be used for the expression of each component present in its expression, so let us start with these:

$$
B_{11} = \dfrac{1}{9} \left( 9B + 2 \dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXY}} + 5\dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}} + \dfrac{\mathrm{d} \left( {\sigma_{33}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}}\right)
$$

$$
B_{12} = \dfrac{1}{9} \left( 9B - 7 \dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXY}} + 5\dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}} + \dfrac{\mathrm{d} \left( {\sigma_{33}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}}\right)
$$

$$
B_{13} = \dfrac{1}{9} \left( 9B + 2 \dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXY}} - 4\dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}} + \dfrac{\mathrm{d} \left( {\sigma_{33}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}}\right)
$$

$$
B_{33} = \dfrac{1}{9} \left( 9B + 2 \dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXY}} - 4\dfrac{\mathrm{d} \left( {\sigma_{11}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}} - 8\dfrac{\mathrm{d} \left( {\sigma_{33}}'-P \right)}{\mathrm{d}\delta}_{\mathrm{StretchXZ}}\right)
$$

This leaves us with only $B_{44}$ to figure out. Luckily, this is simply

$$
B_{44} = -\dfrac{\mathrm{d} {\sigma_{23}}'}{\mathrm{d}\delta}_{\mathrm{ShiftXZ}}
$$

### Elastic constants

The value of $B_{44}$ for a density of $\rho = 0.72000$ is as follows:


with an average value of $B_{44} = -6.05$.

##### Meeting questions/remarks

I need help with figuring out the reduction of independent elastic constant components. Specifically, how to find the solution from the linear system of equations for $C_{11} = C_{22}$ and $C_{66} = \frac{1}{2} C_{11} - C_{12}$.

Is the value of $B_{44}$ looking good. Did I correctly do what Laura asked me to do?