  We perform systematic density functional theory (DFT) calculations to assess
the performance of various exchange-correlation potentials $V_{xc}$ in
describing the chalcogenide GaM$_4$Q$_8$ lacunar spinels (M=Mo, V, Nb, Ta; Q=S,
Se). We examine the dependency of crystal structure (in cubic and rhombohedral
symmetries), electronic structure, magnetism, optical conductivity, and lattice
dynamics in lacunar spinels at four different levels of $V_{xc}$: the local
density approximation (LDA), generalized gradient approximation (GGA),
meta-GGA, and hybrid with fractional Fock exchange. We find that LDA
significantly underperforms GGA and higher level functionals in predicting
lattice constants. LDA also fails to properly describe the magnetism in the
family, and for some compositions it does not find the ground state distorted
crystal structure. The Perdew-Burke-Ernzerhof (PBE) and PBE revised for solids
(PBEsol) GGA functionals perform reasonably well in predicting lattice
constants as well as the electronic structures. We find that the GGA with an
on-site Coulomb interaction (GGA$+U$) is unnecessary to produce a
semiconducting state in the distorted polar $R3m$ phase. Plus Hubbard $U$
values ranging from 1$\sim$2 eV, however, improve the quantitative performance
of the GGA functionals. The meta-GGA functional SCAN predicts reasonable
lattice constants and electronic structures; it exhibits behavior similar to
the GGA$+U$ functionals for small $U$ values. The hybrid functional HSE06 is
accurate in predicting the lattice constants, but leads to a band gap of
$\approx$1 eV in the rhombohedral phase, which is higher than the experimental
estimation of 0.2 eV. Our findings suggest that accurate qualitative and
quantitative simulations of the lacunar spinel family with DFT requires careful
attention to the nuances of the exchange-correlation functional and considered
spin structures.
