Assessing exchange-correlation functional performance in the chalcogenide lacunar spinels GaM$_4$Q$_8$ (M = Mo, V, Nb, Ta; Q = S, Se)

TL;DR

This study systematically evaluates various density functional theory exchange-correlation functionals to determine their effectiveness in modeling the structural, electronic, magnetic, and optical properties of chalcogenide lacunar spinels, highlighting the strengths and limitations of each.

Contribution

It provides a comprehensive comparison of LDA, GGA, meta-GGA, and hybrid functionals in accurately predicting properties of GaM$_4$Q$_8$ lacunar spinels, guiding future computational studies.

Findings

LDA underperforms in predicting lattice constants and magnetism.

GGA functionals like PBE and PBEsol perform reasonably well.

Hybrid HSE06 accurately predicts lattice constants but overestimates band gaps.

Abstract

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.

Tables

Table 1. Table 1: Experimental Jahn-Teller ( T JT subscript 𝑇 JT T_{\mathrm{JT}} ) and Curie ( T C subscript 𝑇 C T_{\mathrm{C}} ) transition temperatures and unit cell volumes (V) for the distorted R ​ 3 ​ m 𝑅 3 𝑚 R3m vanadium and molybdenum lacunar spinels. The vanadium (molybdenum) chalcogenides exhibit acute (obtuse) angluar distortions away from the ideal cubic 60 ∘ . θ m subscript 𝜃 𝑚 \theta_{m} and α r ​ h subscript 𝛼 𝑟 ℎ \alpha_{rh} are obtained at temperatures below T JT subscript 𝑇 JT T_{\mathrm{JT}} , and no significant structural changes have been observed at around T C subscript 𝑇 C T_{\mathrm{C}} .

Compound	$T_{\mathrm{JT}}$ (K)	$T_C$ (K)	${\alpha }_{rh}$ (∘)	${\theta }_m$ (∘)	V${}_{}{}^{F\overline{4}3m}{}_{\text{cell}}{}^{}$ (Å3)	V${}_{}{}^{R3m}{}_{\text{cell}}{}^{}$ (Å3)	Ref.
GaV4S8	44	12.7	59.6	58.4	225.6	224.3	Widmann et al., 2017, Powell et al.,2007
GaV4Se8	41	17.5	59.6	57.7	260.7	259.6	Bichler, 2010,Fujima et al.,2017
GaMo4S8	45	19.5	60.5	61.6	230.1	230.0	Francois et al., 1991,Powell et al.,2007
GaMo4Se8	45	23	60.6	61.4	263.3	262.2	François et al.,1992

Table 2. Table 2: Energy differences (in eV/f.u.) of different magnetic configurations compared with non-magnetic calculations for the cubic phase. E σ subscript 𝐸 𝜎 E_{\sigma} denotes the energy of the highly-polarized state with 5  μ B subscript 𝜇 𝐵 \mu_{B} or 7  μ B subscript 𝜇 𝐵 \mu_{B} per formula unit. An ‘–’ indicates that the state was not stable. PS is an abbreviation for the PBEsol functional and the number in parenthesis is the value of the on-site Coulomb interaction used in the GGA + U 𝑈 +U method.

	GaV4S8		GaV4Se8
	$E_{\sigma }$	$E_{\mu =1{\mu }_B}$	$E_{\sigma }$	$E_{\mu =1{\mu }_B}$
LDA	–	–	–	–
LDA$+U\left(1.0\right)$	–	0.005	-0.172	-0.007
LDA$+U\left(2.0\right)$	-0.404	-0.008	-0.68	0.06
LDA$+U\left(3.0\right)$	-0.946	-0.03	-1.264	-0.051
PBEsol	–	-0.001	-0.072	-0.008
PS$+U\left(1.0\right)$	-0.305	-0.014	-0.562	-0.028
PS$+U\left(2.0\right)$	-0.823	-0.035	-1.123	–
PS$+U\left(3.0\right)$	-1.641	–	-2.174	–
PBE	-0.095	-0.010	-0.327	-0.022
SCAN	-0.875	-0.049	-1.196	-0.064
HSE06	-1.355	0.005	-1.742	-0.092

Table 3. Table 3: Electronic band splitting of the triply degenerate bands in the cubic lacunar spinels at different levels of DFT. γ 𝛾 \gamma quantifies the splitting between the two spin channels. Δ Δ \Delta is the value of the band splitting among the triply degenerate bands at the X point in momentum space near E F subscript 𝐸 𝐹 E_{F} . For the vanadium compounds, these values are tabulated for different spin-magnetic moment states separately. An ‘–’ indicates that the state was not stable.

compound		LDA	PBEsol	PBE	PBE$+U\left(1.0\right)$	PBE$+U\left(2.0\right)$	SCAN	HSE06
GaMo4S8	$\gamma$	0.09	0.15	0.16	0.22	0.28	0.22	0.56
	$\Delta$	0.69	0.65	0.55	0.55	0.55	0.59	0.67
GaMo4Se8	$\gamma$	0.11	0.14	0.16	0.22	0.29	0.22	0.52
	$\Delta$	0.50	0.48	0.38	0.39	0.40	0.44	0.52
GaV4S8	$\gamma$	–	–	1.13	1.58	2.01	1.76	3.1
(5 ${\mu }_B$)	$\Delta$	–	–	0.48	0.44	0.33	0.43	0.46
GaV4S8	$\gamma$	–	0.12	0.23	0.32	–	0.37	0.95
(1 ${\mu }_B$)	$\Delta$	–	0.52	0.41	0.40	–	0.43	0.55
GaV4Se8	$\gamma$	–	1.07	1.17	1.65	2.60	1.82	3.2
(5 ${\mu }_B$)	$\Delta$	–	0.42	0.36	0.34	0.20	0.35	0.51
GaV4Se8	$\gamma$	–	0.20	0.23	0.34	–	0.40	0.81
(1 ${\mu }_B$)	$\Delta$	–	0.37	0.30	0.28	–	0.31	0.41
GaNb4S8	$\Delta$	0.88	0.82	0.69	0.71	0.71	0.73	0.86
GaNb4Se8	$\Delta$	0.63	0.59	0.50	0.52	0.53	0.55	0.66
GaTa4Se8	$\Delta$	0.74	0.71	0.60	0.63	0.63	0.65	0.69

Table 4. Table 4: Wyckoff positions of the transition metals in rhombohedral GaMo 4 S 8   and GaV 4 S 8   after structural relaxation with different DFT functionals. The z 𝑧 z value of the 3 ​ a 3 𝑎 3a and 9 ​ b 9 𝑏 9b sites in space group no. 160 are labeled z 1 subscript 𝑧 1 z_{1} , z 2 subscript 𝑧 2 z_{2} , respectively. PS is an abbreviation for the PBEsol functional and the number in parenthesis is the value of the on-site Coulomb interaction used in the GGA + U 𝑈 +U method.

	GaMo4S8
	$3a$ ($z_1$)	$9b$ ($x$)	$9b$ ($z_2$)
experimentalPowell et al. (2007)	0.4014	0.1956	0.2023
LDA	0.3982	0.1960	0.2022
PBEsol	0.4012	0.1951	0.2012
PS$+U\left(1.0\right)$	0.4014	0.1950	0.2011
PS$+U\left(2.0\right)$	0.4015	0.1949	0.2010
PS$+U\left(3.0\right)$	0.4016	0.1948	0.2010
PBE	0.4020	0.1956	0.2009
SCAN	0.4029	0.1962	0.2006
HSE06	0.4025	0.1959	0.2007
	GaV4S8
	$3a$ ($z_1$)	$9b$ ($x$)	$9b$ ($z_2$)
experimentalPowell et al. (2007)	0.3910	0.1937	0.2013
LDA	0.3944	0.1946	0.1998
PBEsol	0.3913	0.1969	0.2005
PS$+U\left(1.0\right)$	0.3877	0.1966	0.2005
PS$+U\left(2.0\right)$	0.3856	0.1958	0.2005
PS$+U\left(3.0\right)$	0.3834	0.1951	0.2008
PBE	0.3888	0.1972	0.2005
SCAN	0.3852	0.1963	0.2004
HSE06	0.3839	0.1956	0.2005

Table 5. Table 5: The unit cell rhombohedral angle α r ​ h subscript 𝛼 𝑟 ℎ \alpha_{rh} (in degrees) for the R ​ 3 ​ m 𝑅 3 𝑚 R3m phases and the corresponding apical bond angle θ m subscript 𝜃 𝑚 \theta_{m} (in degrees) for the transition-metal cluster at different levels of DFT functional. PS is an abbreviation for the PBEsol functional and the number in parenthesis is the value of the on-site Coulomb interaction used in the GGA + U 𝑈 +U method.

	GaMo4S8		GaMo4Se8		GaV4S8		GaV4Se8
	${\alpha }_{rh}$	${\theta }_m$	${\alpha }_{rh}$	${\theta }_m$	${\alpha }_{rh}$	${\theta }_m$	${\alpha }_{rh}$	${\theta }_m$
experimental	60.47	61.60	60.57	61.43	59.62	58.38	59.56	57.72
LDA	60.00	60.00	60.78	62.76	60.00	60.00	59.55	57.71
PBEsol	60.70	62.29	60.80	62.91	59.56	57.67	59.33	56.59
PS$+U\left(1.0\right)$	60.75	62.50	60.81	63.00	59.28	56.39	59.20	56.00
PS$+U\left(2.0\right)$	60.76	62.59	60.81	63.06	59.16	55.91	59.09	55.60
PS$+U\left(3.0\right)$	60.77	62.66	60.81	63.12	58.99	55.26	58.87	54.75
PBE	60.73	62.53	60.79	63.07	59.36	56.62	59.26	56.09
SCAN	60.76	62.80	60.84	63.27	59.21	55.72	59.12	55.47
HSE06	60.77	62.74	60.79	63.15	59.05	55.38	58.94	55.04

Equations

$$E_{xc}^{HSE06}=\frac{1}{4}{E}_x^{\mathrm{HF},\mathrm{SR}}+\frac{3}{4}{E}_x^{\mathrm{PBE},\mathrm{SR}}+E_x^{\mathrm{PBE},\mathrm{LR}}+E_c^{\mathrm{PBE}}$$

Full text

I Introduction

The lacunar spinel family GaM4Q8 (M = Mo, V, Nb, Ta; Q = S,Se) have garnered attention for decades owing to their fascinating properties, which include metal-insulator transitionsCamjayi et al. (2014), the capability to host skyrmion latticesKézsmárki et al. (2015), and multiferroismWidmann et al. (2017). GaV4S8 and GaMo4S8 are the most well studied materials in this family; they exhibit Jahn-Teller-type structural phase transitions at $\approx$40 K upon cooling, followed by spontaneous magnetic ordering below their Curie temperatures $T_{\mathrm{C}}$ Pocha et al. (2000). The multiple phase transitions – metallic-to-insulating and paramagnetic-to-ferromagnetic – connecting distinct physical states make these transition-metal compounds ideal candidate materials for novel electronic platformsCario et al. (2010).

After decades of continuous studies on various properties of the lacunar spinel compounds, the mechanism of these phase transitions as well as the proper theoretical approaches to describe various electronic states are still unclear. For instance, while the vanadium and molybdenum compounds can undergo symmetry-lowering structural phase transitions at low temperatureSahoo and Rastogi (1993); Francois et al. (1991); François et al. (1992), the niobium and tantalum lacunar spinels remain in the high-symmetry cubic phase over a broad temperature rangeAbd-Elmeguid et al. (2004). One of the possible reasons for this behavior in the family may be attributed to variations in the strength of electron-electron interactions,Sahoo and Rastogi (1993) since electron-correlation effects are expected to be stronger in 3$d$ rather than 5$d$ transition metals. However, there is also evidence that local structural distortions in lacunar spinel compounds could lead to insulating states even in the absence of strong correlationSieberer et al. (2007). In addition, a number of members within the lacunar spinel family exhibit interesting resistive-switching behavior Corraze et al. (2013), making them potential materials for resistive random-access memory (RRAM) materials. Much of the literature attributes the aforementioned features to the special tetrahedral transition-metal clusters within the unit cellLe Beuze et al. (1995); yet, how and why it supports all of these properties remains to be agreed upon. Rastogi et al. (1984); Rastogi and Wohlfarth (1987); Pocha et al. (2000) In order to have a better understanding of the structure-property relationships among the lacunar spinels, a qualitative and possibly quantitative investigation of electron-correlation effects and structural distortions within these materials is needed.

Density functional theory (DFT) simulations are widely used in solid-state materials research owing to the efficiency and accuracy they achieve by replacing the original many-electron interaction problem with an auxiliary independent-particle problem through a suitably constructed exchange-correlation potential ($V_{xc}$). Because DFT simulations can capture the interplay of structural effects on electron-electron interactions and its dependence on determining the ground state, it is an ideal method to study the lacunar spinels with many internal atomic, spin, and orbital degrees-of-freedoms. However, no available $V_{xc}$ can provide the exact description of exchange and correlation, which necessitates benchmarking both common and state-of-the-art density functionals against available experimental data. To that end, it becomes possible to identify the optimal functional for describing and predicting properties in the lacunar spinel family.

In this work, we systematically benchmark the performance of DFT $V_{xc}$ functionals in describing the lacunar spinel family at four rungs of “Jacob’s ladder”, specifically the local density approximation (LDA), the generalized gradient approximation (GGA) as implemented by Perdew-Burke-Ernzerhof (PBE) and PBE revised for solids (PBEsol), the meta-GGA functional SCAN, and the hybrid functional HSE06. Our aim is to identify the best description of the lacunar spinel family from first-principles DFT simulations and where compromises on performance are made so as to to facilitate future studies and predictions. To that end, we investigate the functional dependency of lattice parameters, magnetism, electronic structures, optical properties, and lattice dynamics in both the cubic and Jahn-Teller distorted rhombohedral phases. Our main conclusion is that GGA and higher level $V_{xc}$ functionals are more reasonable than LDA in predicting almost all properties assessed. The GGA functionals with an on-site Coulomb interaction (GGA$+U$) values of $U\approx 2$ eV quantitatively improves functional performance for the electronic structures of the rhombohedral phases. Spin-orbit coupling (SOC) interactions lift orbital degeneracies in the electronic structures of the cubic phases and enable a semiconducting phase to emerge with Hubbard $U$ values ranging from 2 eV to 3 eV. However, SOC does not significantly impact the electronic structures of the rhombohedral phase, where orbital symmetry is already broken by lattice distortions. SCAN and HSE06 are able to predict accurate lattice parameters, but HSE06 leads to band gaps significantly larger than experimental estimations. Our findings suggest that the predicted physical properties of the lacunar spinel family are highly $V_{xc}$ functional dependent. Therefore, it is important to benchmark different $V_{xc}$ performance on properties of interest before further studies. It is likely that the coupling of internal degrees of freedom in lacunar spinels, e.g., local cluster distortion, intra- and inter-cluster electronic and magnetic interactions, underlie the observed fascinating behavior as well as our reported high sensitivity to $V_{xc}$ in this materials family.

II Materials and Methods

II.1 Crystal structure, electrical, and magnetic properties

The crystal structure of the lacunar spinel, also referred to as an A-site deficient spinel (AM4Q8), is derived from the regular spinel (AM2Q4 composition) by removing the interpenetrating FCC A-site sublattice as depicted in Figure 1(a). Upon removing half of the A-site cations occupying the tetrahedral holes in the regular spinel, the space group loses inversion symmetry, reducing from $Fd\bar{3}m$ to $F\bar{4}3m$ (space group no. 216). The structure then undergoes additional internal displacements and spontaneous strains: the previously equidistant M-M network breaks into isolated tetrahedral transition-metal clusters with chalcogenide ligands [M4Q4]5+. In order to quantitatively describe the internal degrees of freedom in the crystal structure, we define $\alpha_{rh}$ as the interaxial angle of the rhombohedral unit cell, and $\theta_{m}$ as the M2-M1-M3 angle centering the apical metal atom along the $C_{3v}$ axis of the M4 cluster, as shown in Figure 1(b).

At room temperature, all lacunar spinels studied here exhibit cubic $F\bar{4}3m$ symmetry with $\alpha_{rh}=\theta_{m}=60^{\circ}$. GaV4S8, GaV4Se8, GaMo4S8, and GaMo4Se8, however, undergo symmetry-lowering structural Jahn-Teller (JT) transitions at $\approx$40 K from $F\bar{4}3m$ to $R3m$ (space group no. 160) Sahoo and Rastogi (1993); Francois et al. (1991); François et al. (1992), followed by spontaneous magnetic ordering at a lower Curie temperature $T_{\mathrm{C}}$ [Figure 1(d)]. These displacive distortions lead to a unit cell of slightly different volume, lattice parameters, and rhombohedral angle $\alpha_{rh}$. The geometry of the metal cluster within the unit cell is also distorted away from its cubic structure and the $\theta_{m}$ angle deviates from the ideal cubic value (60∘). The experimental crystallographic data for the distorted vanadium and molybdenum lacunar spinels are tabulated in Table 1, where we also provide the Jahn-Teller and Curie transition temperatures.

The lacunar spinels are reported to be narrow-bandwidth semiconductors with $\approx$0.2 eV band gaps that vary with temperaturePocha et al. (2000); Cario et al. (2010); Widmann et al. (2017). Early work showed that the valence bands mainly consist of transition-metal $d$ orbitalsShanthi and Sarma (1999). Since the transition-metal M4 clusters are relatively distant from each other with about 4 Å inter-cluster separation [Figure 1(c)], the low-energy valence electronic structure can be described using a molecular orbital (MO) diagram for the cluster (Figure 2). In the cubic phase, the valance bands are triply degenerate with $t_{2}$ symmetry. V, Nb, and Ta chalcogenide lacunar spinels all exhibit $t_{2}^{1}$ occupancies whereas Mo exhibits $t_{2}^{5}$ filling, indicating susceptibility to a first-order Jahn-Teller distortion. After the Jahn-Teller structural distortion, the triply degenerate $t_{2}$ orbital splits into two sets of orbitals, $a_{1}$ and $e$. The relative energy of the two sets of orbitals is occupancy-dependent; the $a_{1}$ orbital is more stable in the vanadium compounds whereas the $e$ orbitals are preferentially stabilized in the molybdenum compounds.

In addition to these structural transitions, the vanadium and molybdenum compounds show spontaneous magnetic ordering at $T_{\mathrm{C}}$ when in the rhombohedral phase. GaV4S8 is also reported to have a complex magnetic phase diagram at low temperatureKézsmárki et al. (2015). The effective local magnetic moment in both the paramagnetic and ferromagnetic phases corresponds to approximately 1 unpaired electron per unit cell, and is mostly localized about the transition-metal clusterPocha et al. (2000) rather than on the individual atomic sites comprising the cluster. GaNb4S8, GaNb4Se8, GaTa4Se8 are paramagnetic at ambient conditions with effective magnetic moments of 1.76 $\mu_{B}$, 1.6 $\mu_{B}$, and 0.7 $\mu_{B}$ per clusterPocha et al. (2005); Abd-Elmeguid et al. (2004). No structural phase transition or spontaneous magnetic ordering are reported in these compounds down to 1.6 KAbd-Elmeguid et al. (2004). Last, we note that the family of materials is also often referred to as Mott insulators owing to the large distance between transition-metal clustersPocha et al. (2000) and not typically because of strong electron-electron interactionsSieberer et al. (2007) although they likely play some role. The semiconducting behavior is typically attributed to variable range hopping (VRH) conductionSahoo and Rastogi (1993) among these separated metal clusters. Nonetheless, the microscopic mechanisms behind the semiconducting nature, as well as the multiple phase transitions, are still under active investigationWidmann et al. (2017); Reschke et al. (2017).

II.2 Exchange-correlation functionals

We use exchange-correlation potentials ($V_{xc}$) at four different levels of approximation to assess the structure and properties of the chalcognide lacunar spinels. The functionals examined include LDA, GGA as implemented by Perdew-Burke-Ernzerhof (PBE)Perdew (1996), and PBE revised for solids (PBEsol)Csonka et al. (2009), meta-GGA functional SCAN as implemented by Sun et al.,Sun et al. (2015) and Heyd-Scuseria-Ernzerhof hybrid functional HSE06Heyd et al. (2005). The $V_{xc}$ in LDA is not derived from first principles, but from Monte Carlo simulations of the uniform electron gas. The functional solely depends upon the local electron density in space and usually provides a good approximation for simple materials (including metals) with electronic states that vary slowly in space. However, the LDA potentials decay rapidly for finite systems while the true exchange-correlation potential has significant non-local contributions; this behavior often leads to overestimation of the binding energyVan de Walle and Ceder (1999) and underestimation of lattice constants in solidsHaas et al. (2009).

To improve on the LDA, GGA functionals that take the gradient of electron density $\nabla n(\mathbf{r})$ into consideration have been developed. The PBE and PBEsol functionals improve the binding energy by roughly an order of magnitude, but have a general tendency to overestimate lattice constantsHaas et al. (2009). Since LDA and GGA functionals are well-known to be unable to predict the insulating state of Mott insulatorsPaier et al. (2006) with strong correlations and nonlocal exchange, the beyond DFT method, DFT$+U$, is typically used to account for such interactions among the localized $d$ electrons. The on-site Coulomb interaction term $U$ favors the on-site occupancy matrix towards fillings that are fully occupied or fully unoccupied and hence a more localized electronic structure within the correlated manifold. Here, we use the GGA functionals PBE and PBEsol with on-site Coulomb interaction (GGA$+U$) with $U$ values of 1.0, 2.0 and 3.0 eV on the the M-metal sites using the formalism introduced by Dudarev et al. Dudarev et al. (1998) to assess the effect of electron correlation in the M4 clusters. The range of $U$ values is based on results from previous computational studiesZhang et al. (2017); Müller et al. (2006); Pocha et al. (2005) and our own preliminary assessments, where reasonable we focused on band gap and magnetic moment predictions.

It comes naturally from the previous two rungs of Jacob’s ladder that the second-order derivative of the electron density should be considered. Meta-GGA functionals are essentially an extension to GGAs whereby the Laplacian of the electron density $\nabla^{2}n(\mathbf{r})$ is also considered. In practice, the kinetic energy density $\tau(\mathbf{r})=\sum_{i=1}^{N_{occ}}\frac{1}{2}\left|\nabla\psi_{i}(\mathbf{r})\right|^{2}$ is used, where the summation runs over the occupied Kohn-Sham orbitals $\psi_{i}(\mathbf{r})$. The recently developed meta-GGA functional SCAN (strongly constrained and appropriately normed semi-local density function) fulfills all known constraints required by the exact density functional, and is reported to have achieved remarkable accuracy for systems where the exact exchange-correlation hole is localized around its electronSun et al. (2015).

Hybrid DFT functionals incorporate a portion of exact exchange interaction from Hartree-Fock (HF) theory with that of a local or semi-local density functional. The semi-empirical hybrid functional B3LYP has been widely used for finite chemical systems and shown more accurate results in thermochemical and electronic propertiesTirado-Rives and Jorgensen (2008); Di Valentin et al. (2006). In periodic solid state systems, one route to incorporate an exact exchange interaction is by means of range separation. In the range separated HSE06 hybrid functional, the short-range (SR) exchange interaction consists of partial contributions from exact exchange and the PBE functional. The long-range (LR) part of the Fock exchange term is replaced by that from the semi-local PBE functional. The correlation term from PBE is used in the HSE06 hybrid functional. The resulting exchange-correlation energy expression is

[TABLE]

The inclusion of exact-exchange interactions in hybrid functionals also partly fixes the self-interaction problem in pure DFT functionals, and can provide accurate descriptions of lattice parameters, bulk moduli and band gaps in periodic systemsHe and Franchini (2012); Hummer et al. (2009); Ramzan et al. (2013).

II.3 Computational details

We perform DFT simulations as implemented in the Vienna Ab initio Simulation Package (VASP)Kresse and Furthmüller (1996); Kresse and Joubert (1999). The projector augmented-wave (PAW) potentialsBlöchl (1994) are used for all elements in our calculations with the following valence electron configurations: Ga ($3d^{10}4s^{2}4p^{1}$), Mo ($4s^{2}4p^{6}4d^{5}5s^{1}$), V ($3s^{2}3p^{6}3d^{4}4s^{1}$), Nb ($4s^{2}4p^{6}4d^{4}5s^{1}$), Ta ($5p^{6}5d^{4}6s^{1}$), S ($3s^{2}3p^{4}$), and Se ($4s^{2}4p^{4}$). Based on convergence test with respect to $k$-point meshes in reciprocal space and plane wave basis set cutoff energies, we use a $\Gamma$-centered $6\times 6\times 6$ mesh with a 500 eV kinetic energy cutoff. For HSE06 calculations, we use a $4\times 4\times 4$ $k$-point mesh and a 400 eV kinetic energy cutoff due to the high computational cost and convergence difficulties for the spin-polarized calculations. Since the lacunar spinels are small-gap semiconductors, we employ Gaussian smearing with a small 0.05 eV width. For density-of-state calculations, we use the tetrahedron method with Blöchl corrections.Blöchl et al. (1994)

We perform full lattice relaxations with different DFT functionals until the residual forces on an individual atom are less than 1.0 meVÅ*-1*. The experimental crystal structures of the lacunar spinels GaV4S8, GaV4Se8, GaMo4S8, GaMo4Se8, GaNb4S8, GaNb4Se8, and GaTa4Se8 are obtained from the Inorganic Crystal Structure Database (ICSD) Bergerhoff et al. (1987) and used as initial inputs for these geometry relaxations. GaTa4S8 is not included here since the structure is not experimentally reported. Both high-temperature cubic and low-temperature rhombohedral phases are investigated for all target compounds. Crystal structures of the rhombohedral phase Nb and Ta compounds are obtained by making a small displacement to their cubic atomic positions along the symmetry-lowering pathway (i.e., from $F\bar{4}3m$ to $R3m$), followed by DFT structural relaxations. All experimental and DFT-relaxed crystal structures are available electronically at Ref.  git, .

The effect of on-site Coulomb interactions on the crystal structures is also investigated at the LDA and GGA functional level. Since the lacunar spinels exhibit various magnetic properties, we also initialize the calculations with multiple possible magnetic configurations for the lattice relaxations. This is a necessary process owing to the multiple metastable spin configurations accessible. The magnetic configuration with the lowest energy is reported as the DFT ground state and used to compare with other functional results. Spin-orbit interactions are also considered in our electronic structure simulations owing to their potentially significant impact on the orbital structure of 4$d$ and 5$d$ transition metals. For spin-orbit coupling (SOC) calculations, we use the fully-relaxed crystal structures from the aforementioned non-SOC simulations. The magnetic moment is set to be 1$\mu_{B}$ per formula unit along the (111) direction for both the cubic and rhombohedral phases.

Zone center ($\mathbf{k}=\mathbf{0}$) phonon frequencies and eigendisplacements for both the cubic and rhombohedral phases of GaV4S8 (within primitive cells) are obtained using the frozen-phonon method with pre- and post-processing performed with the Phonopy packageTogo and Tanaka (2015).

III Results and discussions

III.1 $F\bar{4}3m$ cubic phase

III.1.1 Lattice parameters

The crystal structures of the lacunar spinels with cubic symmetry are fully relaxed with DFT using the different $V_{xc}$ potentials. Figure 3 shows the volumetric error for the DFT ground state unit cell volume relative to the experimental room temperature data. For molybdenum and vanadium compounds, we report the cell volumes of ferromagnetic spin structures with magnetic moments of $1\,\mu_{B}$ and $5\,\mu_{B}$ per formula unit, respectively. The niobium and tantalum compounds are non-magnetic at all DFT functionals levels. See subsubsection III.1.2 for the detailed descriptions of the magnetic moment configurations.

In general, the LDA and PBEsol functionals underestimate the lattice parameters, while PBE predicts larger lattice constants compared with experimental data. LDA has relatively larger deviations (4 % or higher) compared with the GGA results, it is a well-known problem that LDA tends to underestimate the lattice constants. We also check the effect of on-site Coulomb interactions (LDA and GGA$+U$) on lattice parameters with $U$ values up to 3 eV. With increasing on-site Coulomb interaction strength, we find the lattice parameters follow a monotonic increasing trend for both the LDA and GGA functionals. We only show the trend for PBEsol in Figure 3 owing to its similarity with the others. Therefore, a reasonable Hubbard $U$ value could quantitatively improve the lattice parameter predictions in the LDA and PBEsol functionals.

Interestingly, the vanadium compounds exhibit cell volumes that are the most sensitive to the choice of the $U$ value among the lacunar spinels. For instance, the difference in volumetric error induced by $U=3.0$ eV is less than 2 % in GaMo4S8, but is almost 8 % in GaV4S8. The highly spin-polarized electronic state used for the vanadium compounds may be a possible cause of the different sensitivity on the on-site Coulomb interactions.

GGA functionals with $U<3.0$ eV generally predict reasonable cubic lattice constants with less than 4 % error in the cell volumes. The meta-GGA functional SCAN and hybrid functional HSE06 have smaller errors in predicting lattice constants, which give less than 2% error for all 7 compounds studied here. Considering the high computational cost of structural relaxations with HSE06, SCAN should be preferred over HSE06 for lattice parameter estimation unless one requires a specific accuracy requirement or improved forces.

These results suggest that most of the DFT functionals are able to predict reasonable cubic phase crystal structures in the lacunar spinel family with less than 4 % error in the volumes. Generally, we recommend using GGA functionals with a tunable Hubbard $U$ value of 1 eV to 3 eV for lattice parameter predictions. SCAN and HSE06 give more accurate lattice constants compared with lower-level functionals, while SCAN is preferable based on a compromise between accuracy and efficiency.

III.1.2 Magnetism

Experimentally, the vanadium and molybdenum compounds exhibit paramagnetism above their Curie temperatures and exhibit spontaneous magnetic ordering at low temperaturePocha et al. (2000). The magnetically ordered phases can host multiple fascinating magnetic states, including ferromagnetism and complex spin textures (skyrmion lattices) Kézsmárki et al. (2015). Those complex magnetic structures are not considered here. The niobium and tantalum compounds show very weak magnetism and do not exhibit spontaneous magnetic ordering down to 1.6KAbd-Elmeguid et al. (2004). Since the transition-metal clusters are relatively far from each other with a distance of around 4 Å, the inter-cluster magnetic interactions are expected to be quite small. Here we use a ferromagnetic spin configuration on all metal sites within the cluster to model the magnetically ordered phases.

From our DFT simulations, different transition metal clusters are able to hold various magnetic configurations. For the molybdenum compounds, we are only able to stabilize one ferromagnetic configuration in the cubic phase which corresponds to 1 $\mu_{B}$ per primitive cell. The magnetic moments are evenly distributed about the four molybdenum atoms in the Mo4 cluster with negligible contributions from other atomic species. In contrast, the vanadium compounds show numerous stable magnetic configurations (Table 2). Apart from the same ferromagnetic configuration as in the molybdenum compounds, we also find a highly spin-polarized state in the cubic phase. To the best of our knowledge, the electronic structures of this state has not yet been reported before. Recent neutron diffraction studies show that there is one single spin distributed across the V4 cluster instead of residing on a single vanadium ionJef .

In the highly spin-polarized state, the magnetic moment could be $\,5\mu_{B}$ or $7\,\mu_{B}$ per formula unit (f.u.), depending on the DFT functional used. The spin-moments are evenly distributed about the transition-metal cluster, with approximately 1.25 $\mu_{B}$ magnetic moment localized on each vanadium atom. This state is significantly lower in energy than the ferromagnetic configuration with 1 $\mu_{B}$ per formula unit in our DFT simulations.

In several cases, we are not able to stabilize some of the magnetic configurations for vanadium compounds (indicated by ‘–’ in Table 2). For example, LDA only converges to non-magnetic configurations, and PBE$+U=2.0$ eV cannot stabilize the state with 1 $\mu_{B}$ per cluster. In cases where both states can be stabilized, however, the more strongly spin-polarized state is always significantly more stable than the other two configurations (Table 2). We also observe a trend that the highly polarized state is more favored with larger on-site Coulomb interactions or with higher level DFT functionals. In addition, the $\mu=1\,\mu_{B}$ state is usually energetically closer to the non-magnetic state than the highly spin-polarized state. These ground state magnetic configurations are also sensitive to the $V_{4}$ cluster volume, which we show varies with different levels of DFT functional (Figure 4). A larger $V_{4}$ cluster usually supports a higher magnetic moment, while a smaller volume leads to reduced or quenched moments. Our findings show that local structure and magnetic moments are correlated with each other and should be assessed carefully because both depend on the choice of exchange-correlation functional. A recent study utilizing dynamical mean-field theory simulations showed similar results, where the significance of electron correlations in describing the MO Mott physics and structural properties of GaV4S8  is also reportedKim et al. (2018).

There is also evidence that local cluster distortions still exist above the Jahn-Teller temperatureWang et al. (2015), and that the symmetry-broken V4 cluster could lead to different magnetic configurations that are in better agreement with experimental resultsZhang et al. (2017). Why this occurs is attributed to the physics of the distorted phase described in subsection III.2. To that end, we suggest high-resolution detection methods (e.g. pair distribution function) be used to probe the local structures of cubic phase lacunar spinels.

Last, the niobium and tantalum compounds are always non-magnetic in our calculations, regardless of the initial magnetic configuration or choice of DFT functional. This may be a consequence of strong but geometrically frustrated antiferromagnetic interactions in the cubic Nb and Ta clustersPocha et al. (2005) or due to a reduction in the on-site Hund’s interactions, which drives moment formation, from the greater hybridization from the extended $4d$ and $5d$ orbitals.

III.1.3 Electronic structures

We next use the relaxed cubic crystal structure and ground state magnetic configuration of each compound and examine the electronic structures (Figure 5). According to the idealized charge distribution in Ga3+[M4X4]*5+*X${}^{2-}_{4}$, the number of electrons per V4, Nb4, and Ta4 cluster is 7 (since they are in the same column of the periodic table) while there are 11 electrons for a Mo4 cluster; these electrons fill the cluster orbitals depicted in Figure 2.

We selectively show the electronic structures of GaV4S8, GaMo4S8, GaNb4Se8, and GaTa4Se8  in Figure 5 because their S/Se counterpart compounds with the same transition-metal cluster exhibit similar band properties. From our PBE-DFT band structures and projected DOSs for the molybdenum, niobium, and tantalum compounds (Figure 5 a[1,3,4]), we find six valence bands mainly consisting of transition-metal $d$-orbital character with relatively small contribution from the anion $p$ orbitals. The triply degenerate band ($t_{2}$ MO symmetry) is higher in energy than the doubly ($e$ MO) and singly ($a_{1}$ MO) degenerate bands. The band degeneracy and ordering agree well with the cluster MO descriptions of the low-energy electronic structure in these compounds.

The DFT ground state electronic structures of the cubic phase vanadium compounds, however, are significantly different from the other chalcogenides in the lacunar spinel family (Figure 5 a2). All six valence bands in the spin-up channel (green bands) are fully occupied, while only the lower part of the spin-down channel is partially occupied. The triply degenerate spin-down bands are shifted $\approx$1 eV above the Fermi level. Interestingly, the metastable magnetic state with 1 $\mu_{B}$ per cluster exhibits band dispersions that are more similar to the rest of the family (Figure 6) and the magnetic moment of 1 $\mu_{B}$ agrees better with experimental results. It remains unknown whether this DFT ground state in the cubic phase is stable and experimentally accessible; further low-temperature neutron-based scattering measurements, for example, could be used to probe the existence of this spin configuration.

We next quantitatively assess the impact of different $V_{xc}$ as well as on-site Coulomb interactions on the electronic structures by defining two parameters, $\gamma$ and $\Delta$, as shown in Figure 6, which describe the key features in the band structure. $\gamma$ corresponds to the energy difference between different spin-channels of the triply degenerate valence band at the $\Gamma$ point. $\Delta$ quantifies the magnitude of the splitting among the triply-degenerate minority-spin bands at the X point, $k=(1/2,0,1/2)$, near $E_{F}$. The values of $\gamma$ and $\Delta$ for the chalcogenide lacunar spinels at different levels of DFT theory are tabulated in Table 3. For the non-magnetic Nb and Ta compounds, we only report $\Delta$.

All of the cubic phase lacunar spinels are metallic from band theory without considering spin-orbit interactions. Figure 5 shows that the Fermi level, $E_{F}$, is always located within the valence bands, regardless of the magnetic configuration or DFT functional. Specifically, the cubic phase V and Mo compounds are predicted to be half-metals as only one spin channels cross the Fermi level whereas the other spin-channel is fully gapped. In either group VB or VIB transition metal compounds, there is an odd number of electrons in three degenerate bands (Figure 2). For the low spin-polarized states with 1 $\mu_{B}$ magnetic moment per formula unit, we then find that the Fermi level crosses this set of triply degenerate bands and metallicity is protected by the $F\bar{4}3m$ crystal symmetry. Here the splitting of these triply degenerate valence bands throughout the Brillouin zone is quite small; although the $\Delta$ value is functional dependent, it does not exceed 0.7 eV. There is also a small trend of increasing splitting between different spin channels ($\gamma$) with higher levels DFT functionals. We attribute this to the more accurate exchange interactions captured with the more advanced functionals.

The flat valence bands derived from these cluster orbitals lead to large effective masses, and these electrons should be highly localized in real space. This is in agreement with the fact that the transition-metal clusters are far from each other within the unit cell, and the electrons are highly localized within the cluster. One of the possible conduction mechanisms for the lacunar spinels is through variable-range hopping (VRH) Sahoo and Rastogi (1993). It is for the same reason that these compounds have been called “Mott insulators”Pocha et al. (2000).

For the highly-polarized magnetic state in the vanadium compounds, $\gamma$ is much larger than $\Delta$, which makes it different from the rest of the family. In this case, the $\Delta$ term may not be that important since the triply degenerate band is no longer the highest occupied band. The two bands crossing the Fermi level are the $a_{1}$ and $e$ orbitals in the spin-down channel. It is therefore possible to obtain a semiconducting state by shifting the $e$-symmetry orbitals to higher energy and fully occupying the $a_{1}$ orbital. Indeed, we find such a state in GaV4Se8 using the SCAN functional, where the band gap is approximately 60 meV. Whether this highly-polarized state is experimentally accessible, however, remains unclear.

We next report results with SOC included in our simulations. The band structures and DOSs with the PBE functional are shown in Figure 5 b[1-4]. Orbital degeneracy is partly broken compared with the non-SOC band structures. The broken symmetry here is vital for reproducing a semiconducting state since it enables further orbital splitting by increasing electron-electron interactions. Figure 5 c[1-4] show the electronic structures with PBE$+$SOC and a $U$ value of 3.0 eV, where all four compounds exhibit a small but finite band gap. It is interesting to note that both SOC and on-site Coulomb interactions are necessary in order to produce a semiconducting cubic phase for all compounds studied. Intuitively, SOC serves the purpose of symmetry-breaking in the highly-symmetric cubic phase while on-site Coulomb interactions localize electrons and increase repulsion between bands, which eventually lead to a semiconducting state in the cubic lacunar spinels. Although the electron-correlation effect (modeled by the Hubbard $U$) is typically considered more important in 3$d$ transition metals, spin-orbit interactions are more significant in 5$d$ transition metals. Indeed, the lacunar spinel compounds investigated, which include transition metals from the 3$d$, 4$d$, and 5$d$ rows, exhibit similar yet non-identical behaviors. This behavior could be the outcome of competing SOC and on-site Coulomb interactions within these transition-metal cluster systems. It has been shown that spin-orbit coupling effect within the lacunar spinel system could lead to exciting physics (e.g. spin-orbital entangled molecular $j_{\mathrm{eff}}$ states).Kim et al. (2014); Jeong et al. (2017)

Our findings in the cubic phase lacunar spinels indicate that different DFT functionals, as well as various internal electron-electron, spin-orbital interactions, can lead to qualitatively different interpretations of their electronic and magnetic properties. Therefore, extra care in the exchange-correlational functional selection should be taken before pursuing extensive DFT simulations on this family.

III.2 $R3m$ distorted phase

III.2.1 Lattice parameters

In this section, we investigate the DFT functional dependency of properties in the distorted rhombohedral phase. Since only molybdenum and vanadium compounds are reported to exhibit Jahn-Teller-type structural distortions, we benchmark the $V_{xc}$ performance in predicting lattice parameters against available experimental data of GaV4S8, GaV4Se8, GaMo4S8, and GaMo4Se8. In all cases, we use a ferromagnetic spin configuration with 1$\mu_{B}$ magnetic moment per unit cell in our structural relaxations; see subsubsection III.2.3 for a detailed discussion of the magnetic moment configurations.

The Jahn-Teller structural phase transition reduces the crystal symmetry from space group $F\bar{4}3m$ to $R3m$ and occurs with a change in unit cell volume. The relative error of the fully relaxed unit cell volumes for the molybdenum and vanadium compounds are shown in Figure 7. Here, we observe a similar trend as found in the cubic phase. The LDA and PBEsol functionals underestimate the ground state lattice volume, while LDA shows larger deviations from the experimental data. Moreover, structural relaxations of the rhombohedral phases of GaV4S8  and GaMo4S8  with LDA converge to non-magnetic cubic structures, regardless of the initial magnetic moment configurations. LDA is able to stabilize a ferromagnetic configuration in the rhombohedral phase only with on-site Coulomb interactions (LDA+$U$).

PBE overestimates the lattice constants of all four compounds. With increasing value of the on-site Coulomb interactions, the lattice parameters also increase slightly. In the rhombohedral phase, the cell volume of vanadium compounds is not as sensitive to the Hubbard-$U$ value as in the cubic phase, presumably because the electronic structure is semiconducting in the $R3m$ symmetry. SCAN and HSE06 functional again perform quite well with regards to the lattice parameters with less than 2% error.

III.2.2 Internal degrees of freedom

The occupied Wyckoff sites of the transition metals also split upon the transition into the rhombohedral phase, leading to one apical site [M1 in Figure 1(b)] along the $C_{3v}$ distortion axis and three basal atoms [M2, M3, M4 in Figure 1(b)] forming a plane perpendicular to the $C_{3v}$ axis. The Wyckoff positions of the transition metals in GaMo4S8 and GaV4S8 with $R3m$ symmetry (space group no. 160) after structural relaxation with different exchange-correlation functionals are tabulated in Table 4. The selenide compounds show similar functional dependencies and are not shown here. In general, the changes in Wyckoff positions with respect to functional are quite small. However, we find that the $z_{1}$ value in GaV4S8  has a significantly higher functional dependency over that in GaMo4S8 (Figure 8). Both increasing the value $U$ as well as going to higher levels of exchange-correlation functionals favor larger structural distortions in GaV4S8, i.e., keeping the apical V atom far away from the center of the tetrahedral transition-metal cluster. The $z_{1}$ Wyckoff position of the Mo atoms is also largely insensitive to the choice of the DFT functional, possibly owing to the reversed distortion in GaMo4S8, where steric effects might prohibit further distortion.

After the structural phase transition, both the rhombohedral angle $\alpha_{rh}$ and bond angle $\theta_{m}$ in the transition-metal cluster diverge from the cubic 60*∘*, leading to a greater number of internal degrees of freedom in the distorted phase. The latter is correlated with the change in occupied Wyckoff sites of the transition metals . We record these internal bond angles of the four lacunar spinels after structural relaxation using different DFT functionals, the results are shown in Table 5.

Almost all DFT functionals (except for LDA) predict similar results for $\alpha_{rh}$ compared with the experimental data, but in general they give larger local metal-cluster distortions. $\theta_{m}$ values are 1$\sim$2 degrees larger than experimentally reported in the molybdenum compounds, whereas the vanadium compounds show a similar but reversed trend in $\theta_{m}$, i.e., 1$\sim$2 degrees smaller. This agreement is reasonable, and the difference in internal coordinates compared with experiment might come from low-resolution experimental characterization Sieberer et al. (2007). It is also possible that the structural phase transition is incomplete at low temperatureFrancois et al. (1991). We also find a minor trend that higher-level functionals, as well as larger Hubbard-$U$ values, favor larger structural distortions. Since our DFT simulations are performed at 0 K, while lab characterizations are performed at finite temperature, our results are more likely to capture the correct ground state structure where thermal expansion effects are small.

It is interesting to note that the vanadium and molybdenum compounds show reversed structural distortions across the phase transition. This can be explained from the valence MO diagram in Figure 2. In the cubic phase, there is either 1 electron (GaV4S8) or 5 electrons (GaMo4S8) in the valence $t_{2}$ orbitals. Such electronic configurations are Jahn-Teller active, whereby a structural distortion accompanied by orbital-degeneracy lifting could further stabilize the system. Owing to the different electron occupations in the vanadium and molybdenum compounds, their favored electronic configurations require a reversed ordering of the $a_{1}$ and $e$ orbitals. Therefore, the spontaneous structural distortion permits each compound to lift its orbital degeneracy and achieve its favored electronic configuration.

We now summarize the structural benchmark assessment of the cubic and rhombohedral phases. We recommend using the GGA$+U$ method for lattice structure relaxations with a Hubbard-$U$ value of approximately 2 eV to 3 eV. LDA functional should be used with on-site Coulomb interactions for both cubic and rhombohedral phases. The SCAN functional is another reasonable choice that predicts accurate lattice structures. Structural relaxations with HSE06 give very accurate lattice constants, but its high computational costs may be prohibitive if only trying to obtain reasonable crystal structures.

III.2.3 Magnetism

We find that the DFT ground state of both the molybdenum and vanadium compounds in the distorted $R3m$ structure are ferromagnetic with 1 $\mu_{B}$ per formula unit. (Simulation of the complex magnetic phase diagram of lacunar spinels is out of the scope of this work, readers with interest should refer to Refs. Franke et al., 2018; Widmann et al., 2017; Butykai et al., 2017; Ruff et al., 2015; Zhang et al., 2019; Kézsmárki et al., 2015; Ehlers et al., 2016; Leonov and Kézsmárki, 2017.) In the molybdenum compounds, the magnetic moment is evenly distributed about all four Mo atoms, which is the same as what we found in the cubic phase. In the vanadium compounds, however, the apical V atom along the $C_{3v}$ symmetry axis has a large local magnetic moment. The other three basal V atoms have relatively smaller moments that are anti-aligned to the apical spin. This results in a ferrimagnetic configuration in the V4 cluster, giving a net-magnetic moment of 1 $\mu_{B}$ per formula unit. For instance, in the rhombohedral phase of GaV4S8 with the PBE functional, the magnetic moment on the apical V atom is 1.3 $\mu_{B}$ while the three basal V atoms contribute each $-$0.1 $\mu_{B}$. Thus, the net-magnetic moment in one formula unit of GaV4S8  is 1 $\mu_{B}$. Some DFT functionals (e.g. PBE$+U(1.0)$) are able to stabilize a ferromagnetic configuration in GaV4S8  similar to that of GaMo4S8, i.e. evenly distributed, but this magnetic configuration is less energetically favorable compared with the ferrimagnetic configuration. We report properties of the rhombohedral phase vanadium compounds using the ferrimagnetic configuration in the remainder of the manuscript. A recent work that used random-phase approximation correctly reproduced the ground state of GaV4Se8  and explored the coupling between magnetism and structure.Schueller et al. (2019)

III.2.4 Electronic structures

We examine the electronic structures of rhombohedral GaV4S8, GaMo4S8, GaNb4Se8, and GaTa4Se8  in this section, because there is evidence that symmetry-breaking in the transition-metal clusters without distortion of the lattice parameters could lead to different magnetic configurations.Zhang et al. (2017) Such small local distortions may also be challenging to detect with low-resolution characterization techniques; for that reason, we hypothesize that the niobium and tantalum-based lacunar spinels could also exhibit a distorted rhombohedral phase. Therefore, we slightly distort the cubic niobium and tantalum lacunar spinel structures along the symmetry-breaking pathway, and use these geometries as the initial structure for structural relaxations. The structural relaxation settings using different DFT functionals are similar to those used for the cubic phase.

Figure 9 a[1-4] presents the electronic band structures and projected DOSs of these four compounds. The triply degenerate valence bands in the cubic phase split into two sets of orbitals with $a_{1}$ and $e$ symmetry. In the molybdenum compounds, the minority spin $a_{1}$ orbital shifts to higher energy, above the Fermi level, such that five valence electrons occupy the three majority spin orbitals and the minority $e$ orbitals. The vanadium, niobium, and tantalum compounds exhibit different orbital occupations and structural distortions; the $a_{1}$ orbital is further stabilized to lower energy relative to the other five orbitals, and the only one valence electron occupies the $a_{1}$ orbital. Remarkably, we find that the PBE functional is able to open up a small band gap without any on-site Coulomb interactions in the distorted phase GaV4S8  and GaMo4S8. For GaNb4Se8  and GaTa4Se8, the lowest conduction band (minority $a_{1}$) barely touches the Fermi level. However, LDA predicts an unreasonable metallic ground state for these compounds. This finding indicates that the structural distortion alone is sufficient to lift the orbital degeneracy and open a semiconducting gap without strong electron-correlation effect—apparently the additional electron density gradient in $V_{xc}$ through the enhancement factor provides an improved description. In other words, the rhombohedral lacunar spinels may not be strictly described as “Mott” insulators. Our findings confirm the importance of local structural distortions on electronic structures in lacunar spinelsSieberer et al. (2007),

Meanwhile, different GGA functionals qualitatively describe the rhombohedral electronic structures differently. For instance in GaV4S8, PBEsol predicts a metallic state, while PBE opens a small gap of 0.09 eV. In GaMo4S8, PBEsol gives a very small band gap of 0.02 eV while PBE predicts a band gap of 0.13 eV, which is much closer to the experimentally estimated value of 0.2 eV. Therefore, DFT-GGA simulations on this family of compounds should be performed with extra caution with attention focused on the role of the enhancement factor in reducing the self-interaction error.Perdew et al. (1999) We recommended that when using the PBEsol functional to simulate the electronic structures of the lacunar spinels, a slightly larger ($\sim$1 eV) Hubbard $U$ value is used than that for PBE.

Figure 9 b[1-4] shows the electronic structures with the PBE functional and SOC. The effect of including SOC is similar to that in the cubic phase; orbital degeneracy is lifted, but the overall band structures remain similar. Figure 9 c[1-4] shows the effect of now adding a Hubbard on-site Coulomb interaction of 2.0 eV. All four compounds now exhibit a clear band gap, where the conduction band is pushed to higher energy owing to stronger electron-electron interactions. Interestingly, we find that SOC does not seem to play a decisive role in predicting reasonable electronic structures in the rhombohedral phase, whereas the GGA functional alone could predict qualitatively correct behavior. It is possible that SOC plays a less significant role here since crystal symmetry is already broken in the rhombohedral phase, unlike in the highly symmetric cubic phase. Our findings here support our previous hypothesis about the roles of SOC and electron-correlation effect in predicting semiconducting phases.

Figure 10 shows the different band gaps predicted using different DFT functionals. The decreasing band gap in the V-Nb-Ta series from 3$d$ to 5$d$ agrees well with our physical intuition, where electron-correlation effects are expected to decrease. The reason why molybdenum compounds show large band gaps might be caused by different orbital occupations – more valence electrons lead to larger orbital repulsion, which pushes the conduction band to a higher energy level, leading to a larger band gap. We also observe a higher functional rather than compositional dependency on the band gap. With an increasing Hubbard-$U$ value, the band gap increases monotonically. Since a larger on-site Coulomb interaction effectively increases the repulsion between bands with the same spin, it results in a larger gap between the highest occupied and lowest unoccupied bands.

We also observe an interesting similarity in Figure 8 and Figure 10, where the functional dependency of the Wyckoff position $3a$ ($z_{1}$) in GaV4S8  is similar to the trend in the band gap. A larger structural distortion in GaV4S8  also leads to a higher electronic band gap. GaMo4S8  however, does not exhibit such correlated properties. The distinct behaviors of the Mo and V compounds indicate rather different relationships between the structural distortion and the ground state electronic structures. Niobium and tantalum compounds have relatively smaller band gaps compared with vanadium and molybdenum ones, this is consistent with experimental estimation of band gaps. It is also clear that including SOC has a negligible effect on ground state band gap in all four compounds studied here.

SCAN predicts band gaps close to the experimentally suggested 0.2$\pm$0.1 eV valueCario et al. (2010), whereas HSE06 finds approximately a 1.0 eV band gap for the vanadium and molybdenum chalcogenides, and around 0.7 eV for the niobium and tantalum compounds. Since the hybrid functionals partially correct the self-interaction problem in DFT, it is expected to predict more accurate band gaps than lower rung functionals. The larger portion of non-local and range-separated exact exchange interactions included in HSE06, however, might also destroy the balance within the transition-metal cluster, causing the large deviations in the band gaps of the lacunar spinels. It has also been reported that the ferromagnetic ground state is determined by the symmetric exchange interactions Zhang et al. (2017), which could possibly explain the different behaviors of HSE06 from lower-level functionals. More careful experimental characterization of the distorted phase band gaps is required to have a better understanding of which $V_{xc}$ performs the best.

We next examine the effect of the on-site Coulomb interactions and exact exchange interactions on the electronic structures of rhombohedral GaMo4Se8 (Figure 11). The band structure and DOS of GaMo4Se8 using PBE functional with Hubbard-$U$ values of 0.0, 1.0, 2.0, and 3.0 eV are shown in the first five panels. The four panels starting from the right of Figure 11 correspond to the DOS obtained using the HSE06 functional with different portions of exact exchange included as indicated in parenthesis. In general, we observe very similar DOS for the occupied bands. The three valence bands in the spin-up channel are slightly shifted to lower energy relative to the Fermi level, $E_{F}$ with either larger $U$ values or larger amounts of exact exchange. The orbitals beneath these valence orbitals, approximately located at $-1$ eV, are always lower in energy in our HSE06 calculations. Because the HSE06 functional treats all orbitals on the same footing, these lower energy orbitals are also ‘corrected’ in a self-consistent manner, whereas the on-site Coulomb interaction through the $+U$ correction basically forces integer occupancy among only the correlated orbitals.

In addition, we find an increasing trend in the band gap with larger $U$ values or greater contributions of exact exchange to $V_{xc}$. We find that $U=1.0$-2.0 eV leads to very similar electronic structures obtained with HSE06 with 5-10 % exact exchange. Our findings here suggest that a GGA$+U$ functional could be used as an alternative method to study electronic structures in lacunar spinels by reproducing the low-energy electronic structure obtained from a hybrid functional but at lower computational cost. The limitation is that lower lying orbitals that may be of interest are not corrected and therefore cannot exactly reproduce the results of the hybrid functional. Based on our simulation results, we do not suggest using HSE06 functional for electronic structure simulations in the lacunar spinel family.

III.2.5 Optical conductivity

We compute the optical conductivity of the ferrimagnetic rhombohedral phase GaV4S8 and compare our DFT results with the experimental data Reschke et al. (2017) in Figure 12. The experimental data shows the first optical transition occurs at $\approx$2,700 cm*-1* (black symbols), corresponding to an approximate 0.33 eV optical band gap. The optical conductivity then plateaus at approximately 800 $\Omega^{-1}\,\mathrm{cm}^{-1}$ for higher frequencies. Our DFT simulations are able to semi-qualitatively capture the plateau structure, but do not quantitatively reproduce the optical conductivity. With increasing values of $U$, the plateau shifts to higher frequency and this behavior coincides with a larger optical gap as expected from the aforementioned band gap dependencies on the exchange-correlation functional. SCAN functional performs similar to PBE with $U=1.0$-$2.0$ eV. PBE with $U=2.0$ eV gives an optical gap closest to the experimental value. We note that because DFT is a single-particle ground state theory, it may not be the optimal tool to study excited state properties, such as optical conductivity. More accurate simulations, for example, could be pursued by solving the Bethe-Salpeter equation using the GW quasiparticle energiesLiu et al. (2018).

III.2.6 Lattice dynamics

Last, we investigate the exchange-correlation functional dependency on the phonon frequencies in GaV4S8. We present our computed normal mode frequencies for both the cubic and rhombohedral phase with different functionals in Figure 13: LDA, PBE, PBE$+U=1.0$ eV, and SCAN. These calculated values are compared with the experimental Raman/IR frequencies at 80 K reported in Ref. Hlinka et al., 2016, which appear in the first column of Figure 13. Note that the LDA results for the rhombohedral phase are not shown, because the $R3m$ structure is unstable at the LDA level.

We find that the cubic phase phonon frequencies generally agree well with experimental IR/Raman characterization data. LDA and PBE perform reasonably well in reproducing the phonon frequencies. However, with PBE$+U=1.0$ eV, we find a significant decrease in the frequency of the lowest $T_{2}$ phonon mode. The same behavior is also obtained with SCAN functional. Interestingly, this $T_{2}$ mode mainly corresponds to the distortion of the transition-metal cluster along the symmetry-breaking pathway. This could be an evidence of electron-phonon coupling induced structure instabilityRastogi et al. (1984). Although no Raman/IR data is available for the rhombohedral phase, we still see the same phonon mode softening with functional choice upon going from PBE to SCAN. The major difference between different functionals is at the low-frequency domain, where the vibrational modes are mainly related to the transition-metal (V4) clusters. Our finding here suggest that lattice dynamics in the lacunar spinels also have non-negligible functional dependency. More experimental (temperature-dependent) data, however, is required to ascertain the functional that best reproduces the lattice dynamical properties.

IV conclusions

In conclusion, LDA underperforms the other functionals and we recommended to use it only with on-site Coulomb interactions added. The GGA functionals (PBE and PBEsol) perform reasonably well, and the results can be quantitatively improved with on-site Coulomb interactions explicitly added. The meta-GGA functional SCAN is another alternative choice that works well and does not require extra parameterization. Last, the hybrid functional HSE06 predicts accurate lattice structures, but leads to a large electronic band gap in the low-temperature rhombohedral phase. Owing to its high computational cost as well as large deviation in electronic structure predictions, we do not recommend using this hybrid functional for the lacunar spinel family.

All exchange-correlation functionals predict reasonable lattice constants in both the cubic and rhombohedral polymorphs of the chalcogenide lacunar spinels. For electronic structure simulations, the cubic phase is always metallic from band theory and exhibits a narrow transition-metal-derived bandwidth at the Fermi level. Spin-orbit interactions are necessary to predict a semiconducting state in the cubic phase, but not in the rhombohedral phase, at the DFT level. At the LDA and GGA level, on-site Coulomb interactions of 2 eV to 3 eV are recommended to obtain quantitatively improved results. We also found that the PBE functional without on-site Coulomb interactions could predict stable semiconducting states for the rhombohedral phase. Our results obtained with SCAN are similar to PBE$+U(2.0)$ and thus can be safely used in simulations. We also find a highly spin-polarized DFT ground state in GaV4S8, which differs from available experimental data, motivating additional investigations of the magnetic order. We found that the single-particle DFT simulations of the optical conductivity do not give a quantitatively satisfying description of GaV4S8; more sophisticated methods such as with the $GW$ method may be necessary to treat the excited state properties in the lacunar spinels. The LDA and PBE functional, however, perform well in predicting cubic phase phonon frequencies in GaV4S8.

Acknowledgements.

This work in supported by the National Science Foundation (NSF) under award number DMR-1729303. Ab initio DFT simulations are performed on Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant number ACI-1548562 and the DoD-HPCMP (Copper cluster). The authors thank E. Schueller, J. Zuo, R. Seshadri, and S. Wilson for useful discussions.