Comment on 'The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals'

In [1, figure 1] the workflow used to calculate the structure and properties of the materials in C2DB the authors stated that the dynamical stability condition for a structure is not satisfied when elastic constants are negative. Unfortunately, it is an incorrect condition. Moreover, in [1, equation (3)] the authors, for reasons difficult to understand, disregarded shear deformations and calculated only the planar elastic stiffness coefficients C₁₁, C₂₂, and C₁₂, what makes the aforementioned analysis incomplete and insufficient. In addition, even these calculated coefficients in C2DB are erroneous, i.e. the stiffness tensor does not have a proper symmetry resulting from the symmetry of the crystal.

In order to explain what the problem is, some facts from the theory of 2D linear elasticity and elastic stability analysis should be recalled.

The generalised Hooke's law is the linear strain-stress tensor relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{hl} \sigma_{ij}=C_{ijkl}\varepsilon_{kl}\rightarrow \boldsymbol{\sigma}=\mathbb{C}\boldsymbol{\varepsilon}, \nonumber \end{align} \tag{ 1 }$$

where $\boldsymbol{\sigma}$ is the second-rank Cauchy stress tensor, $\mathbb{C}$ is the fourth-rank anisotropic elastic stiffness tensor and $\boldsymbol{\varepsilon}$ is the second-rank small strain tensor ($i,j,k=1,2,3$ for 3D and $i,j,k=1,2$ for 2D problems), from Einstein summation convention repeated indices are implicitly summed over.

From the symmetry of $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ it follows that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{cs1} C_{ijkl}=C_{jikl}=C_{ijlk}, \nonumber \end{align} \tag{ 2 }$$

and from the thermodynamic requirement of existing of a strain energy density function $U(\boldsymbol{\varepsilon})$ (hyperelastic material) [2] such that

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{se} U=\frac{1}{2}\frac{\partial^2 U}{\partial \varepsilon_{ij} \partial \varepsilon_{kl}}\bigg|_{\boldsymbol{\varepsilon}=0}\varepsilon_{ij}\varepsilon_{kl} = \frac{1}{2}C_{ijkl}\varepsilon_{ij}\varepsilon_{kl}\rightarrow U=\frac{1}{2}\boldsymbol{\varepsilon}\mathbb{C}\boldsymbol{\varepsilon}, \nonumber \end{align} \tag{ 3 }$$

additionally

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{cs2} C_{ijkl}=C_{klij}, \nonumber \end{align} \tag{ 4 }$$

and hence number of independent components of four-rank C_ijkl reduces to 21 in 3D [3] and to 6 in 2D [4, 5]. In relations (1) and (3) the fourth-rank tensor notation, employing fourth-rank Cartesian tensor in three or two dimensions, is used. Also different notations for the generalized Hooke's law, relation (1), are in use. The most popular is a non-tensorial Voigt notation that employs in 2D 3x3 matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \begin{array}{@{}c@{}} {\sigma_{11}} \nonumber \\ {\sigma_{22}} \nonumber \\ {\sigma_{12}} \nonumber \\ \end{array} \right]=\left[ \begin{array}{@{}ccc@{}} {C_{1111}} & {C_{1122}} & {C_{1112}} \nonumber \\ {C_{1122}} & {C_{2222}} & {C_{2212}} \nonumber \\ {C_{1112}} & {C_{2212}} & {C_{1212}} \nonumber \\ \end{array} \right] \left[ \begin{array}{@{}c@{}} {\varepsilon_{11}} \nonumber \\ {\varepsilon_{22}} \nonumber \\ 2{\varepsilon_{12}} \nonumber \\ \end{array} \right],\ \label{eqn:HookVoight} \nonumber \end{align} \tag{ 5 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \begin{array}{@{}c@{}} {\hat{\sigma}_{1}} \nonumber \\ {\hat{\sigma}_{2}} \nonumber \\ {\hat{\sigma}_{3}} \nonumber \\ \end{array} \right]=\left[ \begin{array}{@{}ccc@{}} {\hat{C}_{11}} & {\hat{C}_{12}} & {\hat{C}_{13}} \nonumber \\ {\hat{C}_{12}} & {\hat{C}_{22}} & {\hat{C}_{23}} \nonumber \\ {\hat{C}_{13}} & {\hat{C}_{23}} & {\hat{C}_{33}} \nonumber \\ \end{array} \right] \left[ \begin{array}{@{}c@{}} {\hat{\varepsilon}_{1}} \nonumber \\ {\hat{\varepsilon}_{2}} \nonumber \\ {\hat{\gamma}_{3}} \nonumber \\ \end{array} \right]\rightarrow \boldsymbol{\hat{\sigma}}=\boldsymbol{\hat{c}}\boldsymbol{\hat{\varepsilon}}.\ \label{eqn:HookVoight2} \nonumber \end{align} \tag{ 6 }$$

The less popular is a second-rank tensor, called also orthonormal or Mandel, notation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \begin{array}{@{}c@{}} {\sigma_{11}} \nonumber \\ {\sigma_{22}} \nonumber \\ \sqrt{2} {\sigma_{12}} \nonumber \\ \end{array} \right]=\left[ \begin{array}{@{}ccc@{}} {C_{1111}} & {C_{1122}} & \sqrt{2}{C_{1112}} \nonumber \\ {C_{1122}} & {C_{2222}} & \sqrt{2}{C_{2212}} \nonumber \\ \sqrt{2}{C_{1112}} & \sqrt{2}{C_{2212}} & 2{C_{1212}} \nonumber \\ \end{array} \right] \left[ \begin{array}{@{}c@{}} {\varepsilon_{11}} \nonumber \\ {\varepsilon_{22}} \nonumber \\ \sqrt{2} {\varepsilon_{12}} \nonumber \\ \end{array} \right],\ \label{eqn:HookOrthonormal} \nonumber \end{align} \tag{ 7 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[ \begin{array}{@{}c@{}} {{\sigma}_{1}} \nonumber \\ {{\sigma}_{2}} \nonumber \\ {{\sigma}_{3}} \nonumber \\ \end{array} \right]=\left[ \begin{array}{@{}ccc@{}} {{C}_{11}} & {{C}_{12}} & {{C}_{13}} \nonumber \\ {{C}_{12}} & {{C}_{22}} & {{C}_{23}} \nonumber \\ {{C}_{13}} & {{C}_{23}} & {{C}_{33}} \nonumber \\ \end{array} \right] \left[ \begin{array}{@{}c@{}} {{\varepsilon}_{1}} \nonumber \\ {{\varepsilon}_{2}} \nonumber \\ {{\varepsilon}_{3}} \nonumber \\ \end{array} \right]\rightarrow \boldsymbol{\sigma}=\boldsymbol{{c}}\boldsymbol{\varepsilon}.\ \label{eqn:HookOrthonormal2} \nonumber \end{align} \tag{ 8 }$$

The difference between Voigt and second-rank tensor notation is not only by factors of 2 and its square root but is more fundamental. In the Voigt notation, the elements of matrix $\boldsymbol{\hat{c}}$ in equation (6) are not the elements of a second-rank tensor, whereas in the second-rank tensor notation the elements of $\boldsymbol{{c}}$ in equation (8) are the elements of a second-rank tensor in six dimensions for 3D and three dimensions for 2D problems. The fourth-rank tensor notation (1) and second-rank tensor notation (8) are tensorially equivalent [4, 6].

In two-dimensional space, there are five different cell lattice types:

• I.  
  Oblique (parallelogram) (a$\neq$b, $\measuredangle\neq90^{\circ}$),
• II.  
  Rectangular (a$\neq$b, $\measuredangle=90^{\circ}$),
• III.  
  Centered rectangular or diamond (a$\neq$b, $\measuredangle=90^{\circ}$),
• IV.  
  Square (a  =  b, $\measuredangle=90^{\circ}$),
• V.  
  Rhombic or hexagonal (a  =  b, $\measuredangle=120^{\circ}$).

It is clear that symmetry aspects are important in the study of physical phenomena. From symmetry principle: if a crystal is invariant with respect to certain symmetry elements, any of its physical properties must also be invariant with respect to the same symmetry elements and Curie laws, it results that the symmetries of the physical properties of the material may not be lower than the symmetry of the crystal, but may be higher [7, 8].

The symmetry classification of linear elastic materials is not related to crystallography. This is due to the properties of fourth-rank Euclidean symmetric tensors (from the linearity of phenomenological Hooke's law and the properties of two, three-dimensional Euclidean space) [9]. For 3D linear hyperelastic materials, there are eight classes of symmetry and four classes of symmetry for 2D [4, 9].

Necessary and sufficient elastic stability conditions, also called Born stability conditions, in various 3D crystal systems are gathered in [10], but from my best knowledge, there is no such work for 2D crystal systems.

In general, the unstressed crystalline structure is stable with no external loads and in the harmonic approximation, if and only if two independent conditions are fulfilled:

• 1.  
  All its phonon modes have positive frequencies $\boldsymbol{\omega}$ for all wave vectors $\boldsymbol{q}$ (dynamical stability):

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ds} \boldsymbol{\omega^2(q)}>0. \nonumber \end{align} \tag{ 9 }$$
• 2.  
  The strain energy density function, given by the quadratic form (3), is always positive (elastic stability):

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{es} U(\boldsymbol{\varepsilon})>0, \forall{\boldsymbol{\varepsilon}}\neq 0. \nonumber \end{align} \tag{ 10 }$$

It is worth pointing out that some authors incorrectly identify elastic stability (10) with dynamic stability (9) for the long wave limit (i.e. for vanishing wavevectors $\boldsymbol{q}\rightarrow$0) [11, 12]. In the mathematical elasticity this phonon condition is called strong ellipticity and does not imply positive definiteness of the strain energy density function (3), but the opposite implication occurs [3].

It would be quite difficult to check the positive definiteness of the quadratic form (10) directly and it can, therefore, be replaced by equivalent easier conditions [10]:

• 1.  
  All eigenvalues of tensor $\boldsymbol{{c}}$ in second-rank tensor notation (8) are positive,or
• 2.  
  All the leading principal minors of tensor $\boldsymbol{{c}}$ in (8) (determinants of its upper-left k by k submatrix) are positive (Sylvester's criterion).

After this theoretical introduction we can give the form of elastic stiffness tensor $\boldsymbol{{c}}$ in the second-rank tensor notation (8) and the necessary and sufficient elastic stability conditions (10) for all four classes of symmetry for 2D hyperelastic materials.

• 1.  
  Full symmetry (isotropy) $\rightarrow$ Hexagonal lattice (V)(2 elastic constants)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {C}_{IJ} \rightarrow\left[ \begin{array}{@{}ccc@{}} {{C}_{11}} & {{C}_{12}} & 0 \nonumber \\ {{C}_{12}} & {{C}_{11}} & 0 \nonumber \\ 0 & 0 & C_{11}-C_{12} \nonumber \\ \end{array} \right], \label{eqn:Hisotropy} \nonumber \end{align} \tag{ 11 }$$

  C₁₁  >  0 & $C_{11} > \left|C_{12}\right|$ or $\lambda_I=(C_{11} + C_{12})>0$ & $\lambda_{II}=(C_{11} - C_{12})>0$.
• 2.  
  Symmetry of a square, (tetragonal) $\rightarrow$ Square lattice (IV)(3 elastic constants)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {C}_{IJ} \rightarrow\left[ \begin{array}{@{}ccc@{}} {{C}_{11}} & {{C}_{12}} & 0 \nonumber \\ {{C}_{12}} & {{C}_{11}} & 0 \nonumber \\ 0 & 0 & C_{33} \nonumber \\ \end{array} \right], \label{eqn:Htetragonal} \nonumber \end{align} \tag{ 12 }$$

  C₁₁  >  0 & C₃₃  >  0 & $C_{11} > \left|C_{12}\right|$ or $\lambda_I\,=$ $(C_{11} + C_{12})>0$ & $\lambda_{II}=(C_{11} - C_{12})>0$ & $\lambda_{III}=C_{33}>0$.
• 3.  
  Symmetry of a rectangle, (orthotropy)$\rightarrow$ Rectangular (II) & Centered rectangular lattice (III)(4 elastic constants)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {C}_{IJ} \rightarrow\left[ \begin{array}{@{}ccc@{}} {{C}_{11}} & {{C}_{12}} & 0 \nonumber \\ {{C}_{12}} & {{C}_{22}} & 0 \nonumber \\ 0 & 0 & C_{33} \nonumber \\ \end{array} \right], \label{eqn:Horthotropy} \nonumber \end{align} \tag{ 13 }$$

  C₁₁  >  0 & C₃₃  >  0 & $C_{11}C_{22} > {C^2_{12}}$ or $\lambda_I=\frac{1}{2} \left(C_{11} + C_{22}+\sqrt{4C^2_{12}-(C_{11}-C_{22}){}^2}\right)>0$ & $\lambda_{II}=\frac{1}{2} \left(C_{11} + C_{22}-\sqrt{4C^2_{12}-(C_{11}-C_{22}){}^2}\right)>0$ & $\lambda_{III}=C_{33}>0$.
• 4.  
  No symmetry (anisotropy) $\rightarrow$ Oblique lattice (I)(6 elastic constants)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {C}_{IJ} \rightarrow\left[ \begin{array}{@{}ccc@{}} {{C}_{11}} & {{C}_{12}} & {{C}_{13}} \nonumber \\ {{C}_{12}} & {{C}_{22}} & {{C}_{23}} \nonumber \\ {{C}_{13}} & {{C}_{23}} & {{C}_{33}} \nonumber \\ \end{array} \right], \label{eqn:Hanisotropy} \nonumber \end{align} \tag{ 14 }$$

  C₁₁  >  0 & $C_{11}C_{22} > {C^2_{12}}$ & det(C_IJ)  >  0 or $\lambda_I>0$ & $\lambda_{II}>0$ & $\lambda_{III}>0$ (e.g. from the Cardano formula [13]).

The problem can arise if we find C₁₃ and/or C₂₃ other than zero: it is hard to say, in this case, if there is no symmetry at all or, possibly, we have chosen a wrong axis [4]. To avoid this it is recommended to check for all crystals the most general stability condition for anisotropy (14).

The above considerations are not only of a general nature, selected examples of erroneous stiffness tensors and incorrectly verified elastic stability can be found in the Computational 2D Materials Database (C2DB).

As it was written earlier, crystal symmetry implies symmetries of its physical properties, and hence the symmetries of tensors, e.g. the stiffness tensor. The conditions for elastic stability were given in equations (11)–(14).

For example, we can find in the C2DB database:

• Au₂O₂: https://cmrdb.fysik.dtu.dk/c2db/row/Au2O2-GaS-NMSpace group:P-6m2, C₁₁  =  86.93 N m⁻¹, C₂₂  =  87.90 N m⁻¹ and C₁₂  =  103.62 N m⁻¹.Because it is a Hexagonal lattice (V) the stiffness tensor $\boldsymbol{{c}}$ must be isotropic (1) and C₁₁ must be equal to C₂₂. Although all calculated elastic constants are positive, the crystal is not elastically stable because not all required stability conditions (equation (11)) are fulfilled.
• Ta₂Se₂: https://cmrdb.fysik.dtu.dk/c2db/row/Ta2Se2-GaS-FMSpace group:P-6m2, C₁₁  =  75.15 N m⁻¹, C₂₂  =  75.81 N m⁻¹ and C₁₂  =  −45.67 N m⁻¹.Because it is a Hexagonal lattice (V) the stiffness tensor $\boldsymbol{{c}}$ must be isotropic (1) and C₁₁ must be equal to C₂₂. Although calculated elastic constant C₁₂ is negative, the crystal is elastically stable because all mandatory stability conditions (equation (11)) are satisfied.
• Re₂O₂: https://cmrdb.fysik.dtu.dk/c2db/row/Re2O2-FeSe-NMSpace group:P4/nmm, C₁₁  =  17.70 N m⁻¹, C₂₂  =  16.18 N m⁻¹ and C₁₂  =  239.42 N m⁻¹.Because it is a Square lattice (IV) the stiffness tensor $\boldsymbol{{c}}$ must have symmetry of a square (2) and C₁₁ must be equal to C₂₂ (the difference here is more than 9%). Although all calculated elastic constants are positive, the crystal is not elastically stable because not all stability requirements (equation (12)) are met.

This work was partially supported by the National Science Centre (NCN – Poland) Research Project: UMO-2016/21/B/ST8/02450.