# How to calculate 2D elastic constants of boron nitride monolayers from uniaxial strains?

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


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. • Maybe it would be helpful for diagnosing if we could see explicitly how the 3 by 3 lattice vector matrix in the POSCAR file is changed for the deformed structures. Perhaps the plots of the energy against strain with your obtained data points and the fitted quadratic curve might be helpful too. Apr 30 at 13:46
• @CWTan Thank you for advice. I add some details related to the fitting process. Apr 30 at 14:21
• Thanks, I didn't spot anything immediately off with the way you're deforming the structures, but I'm suspecting a sensitivity issue. Since you're using energies and the strain is small, and your trying to extract a second derivative, even small numerical errors in the energy might lead to a bad fit. I gather from your tag that you're using VASP - what EDIFF are you using? And how tightly did you do the convergence tests for ENCUT and KPOINTS? You could compare that to the energy differences of your deformed structures. Apr 30 at 14:38
• @CWTan I used EDIFF = 1E-8 and ENCUT = 600 eV. The kpoint sampling is a 60 * 60 * 1 Gamma-centered Monkhorst-Pack k-mesh. I'll check the calculation condition again. Thanks. Apr 30 at 15:07
• If the conversation gets too long here, I recommend the VASP chat room: chat.stackexchange.com/rooms/109983/vasp Apr 30 at 19:00