I am trying to calculate 2D elastic constants of a hexagonal boron nitride (hBN) monolayer. The unit cell is the primitive cell and the Cartesian coordinate is set as follows.
The elastic constants $C_{j} \:\: (j=x, y)$ have relations with the corresponding strain, $$ \frac{(E-E_{0})}{A_{0}} = \frac{1}{2}C_{j}\varepsilon_{j}^{2}, $$ where $E_{0}$ and $A_{0}$ is the energy and area at the equilibrium.
To find $C_{j}$, I made four deformed structures of the strain $\varepsilon_{j} = \Delta a_{j}/a_{0,j} = \pm 0.005,\; \pm 0.01$. The schematic diagrams for the deformation is the following.
POSCARs of deformed structures have the same fractional coordinates written in the original POSCAR of the optimized geometry, but the lattice vectors are changed according to the above strain. I conducted single-point calculations with deformed structures using VASP and fitted the above equation with a quadratic function.
However, the problem is that $C_{x}$ and $C_{y}$ is not the same (I obtained $C_{x}=326\; \mathrm{Nm^{-1}}$ and $C_{y}=300\; \mathrm{Nm^{-1}}$). Those two values are to be the same due to the hexagonal symmetry [Mouhat, F.; Coudert, F.-X. Necessary and Sufficient Elastic Stability Conditions in Various Crystal Systems. Phys. Rev. B 2014, 90, 224104.].
I suspect the deformed structures I made are wrong somehow, but I am not sure how to define the deformation correctly.
EDIT: I added more data related to the fitting process.
A. The POSCAR lattice vector matrices in the x-direction
$\varepsilon = 0$ (the optimized structure):
2.5122293853799835 -0.0000000000000046 -0.0000000000000000
-1.2561146926399998 2.1756544679020906 0.0000000000000000
0.0000000000000000 -0.0000000000000000 15.0000000000000000
$\varepsilon = 0.005$:
2.52479053 -0.00000000 -0.00000000
-1.25611469 2.17565447 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = 0.01$:
2.53735168 -0.00000000 -0.00000000
-1.25611469 2.17565447 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = -0.005$:
2.49966824 -0.00000000 0.00000000
-1.25611469 2.17565447 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = -0.01$:
2.48710709 -0.00000000 0.00000000
-1.25611469 2.17565447 0.00000000
0.00000000 -0.00000000 15.00000000
B. The POSCAR lattice vector matrices in the y-direction
$\varepsilon = 0.005$:
2.51222939 -0.00000000 -0.00000000
-1.25611469 2.18653274 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = 0.01$:
2.51222939 -0.00000000 -0.00000000
-1.25611469 2.19741101 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = -0.005$:
2.51222939 -0.00000000 -0.00000000
-1.25611469 2.16477620 0.00000000
0.00000000 -0.00000000 15.00000000
$\varepsilon = -0.01$:
2.51222939 -0.00000000 -0.00000000
-1.25611469 2.15389792 0.00000000
0.00000000 -0.00000000 15.00000000
C. The quadratic fittings
The x axis is dimension-less and the unit of y axis is $\mathrm{eV}/\text{Å}^{2}$. Thus, the elastic constants in $\mathrm{Nm}^{-1}$ is given by factoring $\sim 16 \times 2 = 32$ to the 2nd order coefficient.
