  Using density-functional perturbation theory and the Gr\"uneisen formalism,
we directly calculate the linear thermal expansion coefficients (TECs) of a
hexagonal bulk system MoS$_2$ in the crystallographic $a$ and $c$ directions.
The TEC calculation depends critically on the evaluation of a
temperature-dependent quantity $I_i(T)$, which is the integral of the product
of heat capacity and $\Gamma_i(\nu)$, of frequency $\nu$ and strain type $i$,
where $\Gamma_i(\nu)$ is the phonon density of states weighted by the
Gr\"uneisen 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 $\Gamma(\nu)$. We highlight that the Gr\"uneisen parameter associated with
a biaxial strain may not be the same as the average of Gr\"uneisen 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.
