Using density-functional perturbation theory and the Grüneisen formalism, we directly calculate the linear thermal expansion coefficients (TECs) of a hexagonal bulk system ${\mathrm{MoS}}_2$ in the crystallographic $a$ and $c$ directions. The TEC calculation depends critically on the evaluation of a temperature-dependent quantity $I_i\left(T\right)$, which is the integral of the product of heat capacity and ${\mathrm{\Gamma }}_i\left(\nu \right)$, of frequency $\nu$ and strain type $i$, where ${\mathrm{\Gamma }}_i\left(\nu \right)$ is the phonon density of states weighted by the Grüneisen parameters. We show that to determine the linear TECs we may use minimally two uniaxial strains in the $z$ direction and either the $x$ or $y$ direction. However, a uniaxial strain in either the $x$ or $y$ direction drastically reduces the symmetry of the crystal from a hexagonal one to a base-centered orthorhombic one. We propose to use an efficient and accurate symmetry-preserving biaxial strain in the $xy$ plane to derive the same result for $\mathrm{\Gamma }\left(\nu \right)$. We highlight that the Grüneisen parameter associated with a biaxial strain may not be the same as the average of Grüneisen parameters associated with two separate uniaxial strains in the $x$ and $y$ directions due to possible preservation of degeneracies of the phonon modes under a biaxial deformation. Large anisotropy of TECs is observed where the linear TEC in the $c$ direction is about 1.8 times larger than that in the $a$ or $b$ direction at high temperatures. Our theoretical TEC results are compared with experiment. The symmetry-preserving approach adopted here may be applied to a broad class of two lattice-parameter systems such as hexagonal, trigonal, and tetragonal systems, which allows many complicated systems to be treated on a first-principles level.

• Received 4 July 2016
• Revised 13 September 2016

DOI:https://doi.org/10.1103/PhysRevB.94.134303