Predictions of Electronic, Transport, and Structural Properties of Magnesium Sulfide (MgS) in the Rocksalt Structure ()
1. Background and the Motivation for This Work
Magnesium sulfide (MgS) belongs to the group of alkaline earth sulfides. The interests of scientists in alkaline earth sulfides have been increasing due to their potential applications in multicolor thin-films electroluminescent and magneto-optical devices [1]. Due to its large band gap, MgS has wide applications in blue-light emitting diodes and optical storage devices [2][3]. It is a favorable material for active solar blind UV detection [4]. It is also useful in lithium-ion batteries as electrode material [5]. MgS exists in three structures (rock salt, zinc blende, and wurtzite), with the rock salt structure as the most stable one. Many theoretical [1][6]- [20]studies have been carried out on MgS. Table 1 below lists both calculated and experimental band gap values for MgS. The Hartree-Fock method has been applied by Pandey et al. [1]to compute the electronic band structure. Their study found an indirect band gap value of 6.48 eV. Drief et al. [6]reported first principle calculations for rock salt MgS, employing a local density functional approximation (LDA) potential and the full potential linearized augmented plane wave (FP-LAPW) method. They predicted an indirect band gap of 2.208 eV. There are several other calculated, [7][8][9][10][11], indirect band gap values of rock salt MgS, obtained with ab-initio LDA potentials. Five (5) of these values were within the range of 2.5 to 2.7 eV. Recently, Tairi et al. [12]studied the structural and electronic properties of rock salt MgS by using
![]()
Table 1. Theoretical band gaps of rock salt magnesium sulfide (MgS) in the literature.
aReference [6]; bReference [7]; cReference [8]; dReference [9]; eReference [10]; fReference [11]; gReference [12]; hReference [13]; iReference [14]; jReference [15]; kReference [16]; lReference [17]; mReference [18]; nReference [19]; oReference [20]; pReference [1].
the full potential linearized augmented plane wave method. They obtained an indirect band gap of 2.745 eV and 4.248 eV with a generalized gradient approximation (GGA) and mBJ potentials, respectively. Several other groups [13][14][15][16]have studied the properties of rock salt MgS by utilizing ab-initio GGA potentials. Their calculated band gap values range from 2.6 - 2.79 eV. An indirect band gap of 2.7 eV was found by Kalpana et al. [17]utilizing tight binding linear muffin tin orbitals (LMTO), with GGA potentials. Ching et al. [18]calculated the band structure of MgS by applying the self-consistent orthogonalized linear combination of atomic orbital (OLCAO) method. Their calculated, indirect band gap was 4.59 eV. The indirect band gap range of MgS was found to be 4.0 - 4.80 eV, [11][12][19][20]using hybrid potentials.
Experimental investigations performed on this material are very limited. Taleatu et al. [21]experimentally determined a direct band gap for MgS thin films (3.54 eV, 3.73 eV, and 3.14 eV) using X-ray diffraction technique. They have suggested that there is difficulty in the determination of the nature and value of the band gap of MgS. They noted difficulties stemming from the spontaneous oxidation of MgS in air and its reactivity with water.
The disagreement among theoretical results indicates the need for further investigations of the band gap and related properties of MgS. The above disagreement and current and potential applications of MgS motivated our ab-initio, self-consistent calculations of electronic, transport, and bulk properties of MgS. Our computational method has been successful in correctly describing and predicting properties of semiconductors in the past [22]. Hence, this work is expected to provide an accurate DFT description of MgS, as most previous, ab-initio DFT calculations are presumed to have underestimated the band gap of MgS. Further, in the case of MgS, a reliable experimental value of the band gap of bulk MgS is not available. Indeed, the experimental values noted above pertain to thin films and are subject to significant uncertainties, according to the authors. We describe below our method that led to highly accurate results in dozens of previous ab-initio DFT calculations.
2. Computational Method
We recall the fundamental theorems of DFT, along with the correct understanding of them, in order to facilitate the discussions on our method. For a system of N electrons in a box, subject to an external potential v(R), the first theorem states that v(R) is a unique functional of the charge density n(R), except for an additive constant. A corollary of this theorem is that the energy content of the Hamiltonian E(Ψ) = <Ψ|H|Ψ> is a unique functional of the charge density n(R). This energy content of the Hamiltonian is simply the sum of the occupied energies. Hence, a second corollary of the first theorem is that the spectrum of the Hamiltonian is a unique functional of the charge density [22].
The second DFT theorem states that E(Ψ’) = <Ψ'|H|Ψ'> reaches its minimum if Ψ'
Ψ of the correct ground state, for arbitrary variations of Ψ' that keep the number of particles constant [22]. Specifically, the second theorem means that E(Ψ') > E(Ψ), where Ψ is for the ground state, for all Ψ' ≠ Ψ. For over 50 years, the community at large understood that the minimization of E (Ψ’), using a single basis set Ψ', provides the DFT solution for the system under study. It does not. The minimization resulting from self-consistency with a single basis set only leads to a stationary solution among an infinite number of stationary solutions. Indeed, every reasonable basis set (without linear dependency) leads to a stationary solution upon the attainment of self-consistency. The generalized minimization unveiled by Bagayoko [22]is required to attain, verifiably, the ground state of the system. This generalized minimization entails several self-consistent calculations with successively augmented basis sets until the occupied energies totally cease to decrease. Then, and only then, do we have the ground state of the system, and hence, the DFT description of the material. The first basis set to produce this description, indubitably, generates the ground state charge density upon reaching self-consistency.
Our calculations employed a method that ensures the attainment of the ground state of a material without the use of very large basis sets that are over-complete for the description of the ground state [22]- [31]. The Ceperley and Alder [32]LDA potential was parameterized by Vosko et al. [33]. We used the linear combination of atomic orbitals (LCAO) and the Bagayoko, Zhao, and Williams (BZW) method, [23][24][25][29]as enhanced by Ekuma and Franklin (BZW-EF) [28]. We employed the electronic structure package developed at the US Department of Energy Laboratory in Ames, Iowa [34][35].
Our calculations started with a basis set slightly larger, by one orbital, than the minimum basis set. We performed self-consistent Calculation I with this small basis set. The second self-consistent calculation was done with a basis set that consists of the first basis set, from Calculation I, augmented by one orbital corresponding to an excited state in the ionic species of the system. A comparison of the occupied energies from the two calculations shows that the occupied energies from Calculation II were lower than or equal to their corresponding values from Calculation I. This comparison follows the setting of the Fermi level to zero in both sets of eigenvalues. This lowering of occupied energies, from their values from Calculation I, indicates that the basis set used in Calculation I does not lead to the ground state. After augmenting the basis set of Calculation II with one orbital, we performed Calculation III and compared the resulting occupied energies with those from Calculation II. The successive, self-consistent calculations, with augmented basis sets, were stopped only after three consecutive ones led to the same occupied energies.
This situation signifies that we reached the lowest possible values of the occupied energies, i.e., the ground state of the material, as required by the second DFT theorem.
The above generalized minimization of the energy, one that is far more than the minimization resulting from self-consistent iterations with a single basis set, is required by the second DFT theorem [22]. Calculations with augmented basis sets larger than that of the third of the referenced three (3) calculations, also led to the same occupied energies. While the above three (3) successive calculations, and others with larger, augmented basis sets, led to the same occupied energies, they generally produced some different (i.e., lower), unoccupied energies, including some of the lowest ones. Generally, the larger the basis set, the lower were some unoccupied energies as compared to their counterparts obtained with the first of the calculations leading to the ground state. Clearly, the question arises as to which of the calculations producing the ground state provides the DFT description of the material. Bagayoko [22]showed, using a corollary of the first DFT theorem, that the first of the three (3) consecutive calculations leading to the ground state, with the smallest basis set, is the one producing the DFT description of the material. The basis set of this calculation is referred to as the optimal basis set. Indeed, the corollary simply states that the spectrum of the Hamiltonian is, like the external potential, a unique functional of the charge density. Thus, the occupied energy content of the Hamiltonian is a unique functional of the charge density [22]. The sum of the occupied energies in this spectrum is the energy content of the Hamiltonian. In light of this corollary, unoccupied energies different (generally lower) from their corresponding values obtained with the optimal basis set cannot belong to the spectrum of the Hamiltonian, a unique functional of the charge density.
There is another way to prove that the first of the calculations producing the ground state is the calculation providing the DFT description of the material. This proof rests on a theorem for eigenvalue due to Rayleigh. The theorem applies to eigenvalues of the same equation obtained with two different basis sets where the larger one includes entirely the small one. In this condition, the theorem states that the ordered eigenvalues obtained with the larger basis set are lower than or equal to their corresponding values obtained with the smaller basis set. All the calculations producing the ground state lead to the same self-consistent charge density and the same Hamiltonian, within computational uncertainties. We note that the same Hamiltonian does not mean the same Hamiltonian matrix, as the dimension of the matrix varies with the size of the basis set. Hence, after the first calculation producing the ground state, other calculations leading to the ground state and some different (lowered) unoccupied energies do not provide the DFT description of the material. Indeed, the unoccupied energies lowered from their corresponding values obtained with the optimal basis set result from a mathematical artifact stemming from the Rayleigh theorem; they are not physically correct eigenvalues of the Hamiltonian which did not change from its value obtained with the optimal basis set.
Computational details in this work follow. Magnesium sulfide has a cubic structure. It is in the space group
. The positions of the Mg and S ions in the irreducible zone are (0, 0, 0) and (1/2, 1/2, 1/2), respectively. The lattice constant we utilized in our self-consistent calculations is 5.200Å [36]. Our first self-consistent calculations were for Mg2+ and S2−; they provided the input orbitals for the study of solid MgS. We expanded the radial parts of the wave function in terms of Gaussian functions with even-tempered Gaussian exponents. The numbers of Gaussian orbitals utilized for the s, p, and d orbitals for Mg2+ were 18, 16, and 16, respectively, and 18, 18, and 16, respectively, for S2−. The smallest and largest exponents were respectively 0.1822 and 0.11 × 106, for the Mg2+, and 0.1489 and 0.44 × 105 for the S2−. Our criterion for self-consistency is to have no more than a difference of 10−5 between potentials from two consecutive iterations. It typically took 60 iterations to satisfy this condition. We utilized 81 k points in the irreducible Brillouin zone for the production of the final, self-con- sistent bands. Extensive testing, over the years, has established the 20 divisions for each of the intervals L-Γ, Γ-X, X-K, and K-Γ that are amply adequate for the description of this structure.
3. Results
The valence orbitals for the successive, self-consistent calculations for MgS are shown in Table 2, along with the resulting band gaps. As per the above explanation of our method, Calculation III is the one that provides the DFT description of MgS. While Calculations IV-VI also led to the absolute minima of the occupied energies (i.e., the ground state), they produced some low-laying unoccupied energies lower than their corresponding values obtained in Calculation III, with the optimal basis set. The electronic energies from Calculation III are shown in Table 3. The electronic energy bands, densities of states, effective masses, and bulk modulus discussed below are from Calculation III.
Figure 1 compares the results from Calculations III (solid lines) and IV (dashed lines). The basis set of Calculation IV consists of that of Calculation III plus the 4s0 on S2−. The occupied bands from the two calculations are the same. Calculations V and VI also produced the same occupied energies as Calculation III.
Figure 2 and Figure 3 respectively show the calculated, total density of states (DOS) and the partial densities of states (pDOS), derived from the bands from Calculation III. We found a total valence bandwidth of 12.597 eV. From the
![]()
Table 2. The successive calculations of the BZW-EF method, with the valence orbitals of Mg2+ and S2− in Columns 2 and 3, respectively. Column 4 and 5 show the total number of valence functions and the band gap, respectively. The employed, room temperature experimental lattice constant was 5.200Å. The superscript zero indicates an empty orbital. The optimal basis set is that of Calculation III. The corresponding, indirect band gap for bulk MgS, from Γ to X, is 3.278 eV, at room temperature.
![]()
Figure 1. Calculated, electronic energy bands of MgS, from Calculations III (solid lines) and IV (dashed lines). The perfect superposition of the occupied bands from Calculations III, IV, and IV signifies the attainment of the absolute minima of the occupied energies, i.e., the ground state, by Calculation III.
![]()
Figure 2. Total density of states (DOS) for rock salt MgS, derived from the bands from Calculation III. Zero, on the horizontal axis, indicates the position of the Fermi energy.
![]()
Table 3. Electronic energies at high symmetry points in the Brillouin zone for bulk rock salt MgS, as obtained from Calculation III, with the optimal basis of the LDA-BZW-EF method. The lattice constant is 5.200 Å, at room temperature.
![]()
Figure 3. Results for partial densities of states for rock salt MgS, derived from the bands from Calculation III. The zero on the horizontal axis indicates the position of the Fermi energy.
DOS figure, the width of the lowest laying valence band is 1.85 eV, while that of the group of upper valence bands is about 5.0 eV. As per the partial densities of states in Figure 3, Mg-s and Mg-p dominate in the lowest, conduction band. There is a small contribution from S-p. In contrast, S-p strongly dominates upper valence bands, with minimal contributions from Mg-p and Mg-s. Sulfur s (S-s) clearly dominates in the lowest laying valence band. These results are in qualitative agreement with findings of Drief et al. [6]for the valence states.
Effective masses are important quantities in the theoretical calculation of transport properties; such properties include electrical conductivities and the Seebeck coefficient [37]. As per the content of Table 4, we have performed calculations of electron effective masses around the minimum of conduction band, at X, and around the lowest conduction bands at Γ. The effective masses of the light and heavy holes were derived from the uppermost valence bands at Γ. The calculated effective masses are shown in Table 4, for the various directions, in units of the mass of the electron (mo). The effective masses of heavy Hole 1 and heavy Hole 2 are equal, except in the (Γ-K)110 direction. Their difference in that direction is due to the splitting of the bands in the (Γ-K)110 direction by the Coulomb crystal field. The quasi-isotropic nature of the electron effective masses at Γ is in contrast to that of the clearly anisotropic ones at X. The calculation performed by Ghebouli et al. [14]reported the electron effective mass for rock salt MgS, at the X point, to be 0.2179 mo. This value is smaller than our predicted ones of 0.324 and 0.319 mo at X, in the transverse direction. Our result is 0.503 for the longitudinal direction, at X. We found no experimental values for these effective masses. We expect future measurements to confirm our predictions in Table 4.
Figure 4 shows the total energy versus the lattice constant. The range of the lattice constant over which we calculated the total energy was from 5.10 to 5.30 Å. Our predicted equilibrium lattice constant, at the minimum of the total energy, is 5.183Å. With this lattice constant, we predicted a zero temperature band gap of 3.512 eV, larger than the room temperature value by 0.234 eV. For our purposes, no general fitting of the total energy to an equation of state is necessary; we fitted the curve in the immediate vicinity of the minimum in order to calculate the bulk modulus. Our calculated bulk modulus of 79.76 GPa is in excellent agreement with the experimental value of 78.9 ± 3.7 GPa [36].
![]()
Figure 4. The graph of the total energy versus the lattice constant, for rock salt MgS, as obtained from calculations using the optimal basis set.
![]()
Table 4. Calculated, effective masses for rock salt MgS, in units of free electron mass (mo): Me indicates an electron effective mass at the X-point and Γ-point, in the conduction band; Mhh and Mlh represent the heavy and light holes effective masses, respectively.
4. Discussion
Our calculated, indirect band gap is 3.278 eV. This value is higher than the predicted value from calculations with ab-initio LDA and GGA potentials, with ranges of 2.208 - 2.69 eV, for LDA, and 2.64 - 2.79 eV, for GGA. One cannot draw conclusions from results of calculations with ad hoc potentials, given their variability with the adjustable parameters. The Green function and dressed Coulomb approximation (GWA) is not DFT theory. For this reason, we do not make comparisons between our LDA findings and those from GWA (or GW) calculations. We suggest the following explanation for the differences between our results and the previous, ab-initio DFT ones in Table 1.
To our knowledge, previous DFT calculations mostly employed single basis sets. The chance for a self-consistent calculation with a single basis to lead to the ground state of the material is practically zero. Rather, such a calculation leads to one of the potentially infinite number of stationary solutions. Our calculations in Table 2 produced such stationary solutions. The progressive increase of the basis set allowed us to identify the first one among them that 1) yields the ground state energy and 2) is not over-complete for the description of the ground state. Basis sets that contain the optimal one and are over-complete produce the same occupied energies obtained with the optimal basis set, but they generally produce some unoccupied energies, including lowest laying ones, that are lower than their counterparts obtained with the optimal basis set. Consequently, these basis sets produce band gaps that are smaller than corresponding, experimental ones. This phenomenon, we suggest, explains the differences between our calculated band gap and those from previous ab-initio DFT (i.e., LDA or GGA) calculations as shown in Table 1.
Similarly, our calculated electron effective masses are larger than the previous results discussed above. In general, a lowering of unoccupied bands results in a decrease of their degree of flatness. A decrease of the degree of flatness around the conduction band minimum clearly leads to a decrease of the electron effective mass. Hence, the above unphysical lowering of unoccupied energies is responsible for the underestimation of the electron effective masses in a way that is linked to the underestimation of the band gap.
5. Conclusion
We have employed the LCAO formalism and an LDA potential to perform ab-initio, self-consistent calculations to predict electronic and related properties of rock salt MgS. By searching for and by verifiably attaining the ground state of the material, as required by the second DFT theorem, our results have the full, physical content of DFT and agree with available experimental results. Our calculated, room temperature indirect band gap is 3.278 eV. Our predicted equilibrium lattice constant and zero temperature band gap are 5.183Å and 3.512 eV, respectively. The calculated bulk modulus (79.76 GPa) agrees very well with the measured value of 78.9 ± 3.7 GPa. It is hoped that our predictions will spur new experimental investigations of properties of MgS. As noted in the background and motivation section, we only found bang gap measurements for thin films of MgS. Possible methods for measuring band gaps include optical absorption, optical transmittance, spectroscopic ellipsometry, and UV to visible and near infrared absorption, with diffuse reflectance.
Acknowledgements
The funding sources for this work are the US Department of Energy (DOE), National Nuclear Security Administration (NNSA) (Award No. DE-NA0003679), the National Science Foundation (NSF, Award No. 1503226DE), LaSPACE, and LONI-SUBR.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.