# Band structures and $\mathbb{Z}_{2}$ invariants of two-dimensional
transition metal dichalcogenide monolayers from fully-relativistic
Dirac–Kohn–Sham theory using Gaussian-type orbitals

Marius Kadek<EMAIL_ADDRESS>Department of Physics, Northeastern
University, Boston, Massachusetts 02115, USA  Hylleraas Centre for Quantum
Molecular Sciences, Department of Chemistry, UiT The Arctic University of
Norway, N-9037 Tromsø, Norway  Baokai Wang  Department of Physics,
Northeastern University, Boston, Massachusetts 02115, USA  Marc Joosten
Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, UiT
The Arctic University of Norway, N-9037 Tromsø, Norway  Wei-Chi Chiu
Department of Physics, Northeastern University, Boston, Massachusetts 02115,
USA  Francois Mairesse  Laboratory of Theoretical Chemistry, University of
Namur, B-5000 Namur, Belgium  Michal Repisky  Hylleraas Centre for Quantum
Molecular Sciences, Department of Chemistry, UiT The Arctic University of
Norway, N-9037 Tromsø, Norway  Department of Physical and Theoretical
Chemistry, Faculty of Natural Sciences, Comenius University, Bratislava,
Slovakia  Kenneth Ruud  Hylleraas Centre for Quantum Molecular Sciences,
Department of Chemistry, UiT The Arctic University of Norway, N-9037 Tromsø,
Norway  Norwegian Defence Research Establishment, P.O. Box 25, 2027 Kjeller,
Norway  Arun Bansil  Department of Physics, Northeastern University, Boston,
Massachusetts 02115, USA  et al.

###### Abstract

Two-dimensional (2D) materials exhibit a wide range of remarkable phenomena,
many of which owe their existence to the relativistic spin–orbit coupling
(SOC) effects. To understand and predict properties of materials containing
heavy elements, such as the transition-metal dichalcogenides (TMDs),
relativistic effects must be taken into account in first-principles
calculations. We present an all-electron method based on the four-component
Dirac Hamiltonian and Gaussian-type orbitals (GTOs) that overcomes
complications associated with linear dependencies and ill-conditioned matrices
that arise when diffuse functions are included in the basis. Until now, there
has been no systematic study of the convergence of GTO basis sets for periodic
solids either at the nonrelativistic or the relativistic level. Here we
provide such a study of relativistic band structures of the 2D TMDs in the
hexagonal (2H), tetragonal (1T), and distorted tetragonal (1T’) structures,
along with a discussion of their SOC-driven properties (Rashba splitting and
$\mathbb{Z}_{2}$ topological invariants). We demonstrate the viability of our
approach even when large basis sets with multiple basis functions involving
various valence orbitals (denoted triple- and quadruple-$\zeta$) are used in
the relativistic regime. Our method does not require the use of
pseudopotentials and provides access to all electronic states within the same
framework. Our study paves the way for direct studies of material properties,
such as the parameters in spin Hamiltonians, that depend heavily on the
electron density near atomic nuclei where relativistic and SOC effects are the
strongest.

##  I Introduction

Two-dimensional (2D) materials [1, 2] are solids in which atoms or compounds
are bound together by strong bonds (_e.g._ covalent or ionic) along two
spatial dimensions, confining electron transport to a plane. In the out-of-
plane dimension, weaker (van der Waals) forces enable the synthesis of
materials with thicknesses of only a few atomic layers [3, 4, 5, 6]. 2D
materials have recently become very attractive as they exhibit remarkable
transport [7, 8], topological [9, 10, 11], thermoelectric [12], and
optoelectronic [13] properties that can be exploited to develop novel devices
for quantum computing [14, 15], field-effect transistors [16, 17], low-power
logic devices [18], and strain-controllable light-emitting devices [19].
Advantages of 2D materials over conventional solids stem from their thin
surfaces which enable considerable manipulation of their properties that can
be controlled using defects, adatoms, and electric-field gating [20, 21, 22].
In addition, atomically thin materials can be manually assembled to form
multilayered heterostructures [23, 24, 17, 11] to combine functionalities of
individual layers [25].

Monolayers of transition-metal dichalcogenides (TMDs) [26] are 2D materials of
type $MX_{2}$, where $M$ is a transition-metal atom (Mo, W, …), and $X$ is a
chalcogen atom (S, Se, Te), and constitute basic building blocks for many
heterostructures. TMDs host many exotic physical phenomena, such as the
quantum spin Hall (QSH) [11] and nonlinear anomalous Hall [27] effects,
higher-order topology [28], giant Rashba spin-splittings of valence bands
[29], as well as various correlated phases, _e.g._ charge density waves [30,
31], superconductivity and ferromagnetism [32]. TMDs typically contain
elements from the lower part of the periodic table, where the relativistic
theory of electrons is unavoidable. Many interesting properties of the TMDs
owe their existence to the relativistic link between the spin and orbital
degrees of freedom, _i.e._ the spin–orbit coupling (SOC) effects. The ability
to control and tune SOC in TMDs opens new possibilities for spintronics and
valleytronics devices based on non-magnetic materials, where spin is
manipulated by electrical means only [33, 34, 35, 36]. Furthermore, SOC
generates opposite effective Zeeman fields at the $K$ and $K$’ valleys of the
TMDs, which enables formation of Cooper pairs and opens the possibility of
observing topological superconductivity and Majorana fermions [37, 38].
Recently, antisite defects in TMDs were proposed to be suitable for hosting
solid-state spin qubits, where SOC enables transitions between different spin
configurations required for qubit operations [39].

First-principles computational approaches enable systematic, parameter-free,
material-specific predictions in novel materials. However, due to the presence
of heavy elements in the TMDs, it is important that relativistic and SOC
effects are accounted for in the theoretical framework. In this connection,
the four-component (4c) Dirac Hamiltonian [40, 41, 42]

$\hat{H}=\begin{pmatrix}V(\bm{r})&c\bm{\sigma}\cdot\bm{\pi}\\\
c\bm{\sigma}\cdot\bm{\pi}&V(\bm{r})-2c^{2}\end{pmatrix},$ (1)

is the commonly accepted gold standard, which is exact in the limit of
noninteracting relativistic electrons. Here, $\bm{\sigma}$ are the Pauli
matrices, $\bm{\pi}\equiv-i\bm{\nabla}+\bm{A}$ is the canonical momentum
operator of the electron, $c$ is the speed of light, and atomic units
($e=\hbar=m_{e}=1$) are used. The scalar and vector potentials $V$ and
$\bm{A}$, respectively, describe the interaction of the electron with external
electromagnetic fields, where $V$ typically contains interactions with the
nuclear charges. The lower components of the Dirac bispinor wave functions are
associated with the negative-energy states, and various strategies for their
elimination lead to a number of approximate two-component (2c) schemes [43,
44].

The multiple wave-function components needed in relativistic simulations
significantly increase the computational cost as larger matrix and vector
dimensions are involved. Kohn–Sham (KS) density functional theory (DFT) [45,
46] brings relativistic calculations to an affordable level and remains the
method of choice for large solid-state systems containing heavy elements.
However, we note the work of Yeh _et al._ [47] on the self-consistent $GW$
method at the relativistic 2c level of theory. In KS theory, the potential in
Eq. (1) also contains terms that depend on the electron density and its
gradient: the mean-field Coulomb potential generated by other electrons, and
the exchange–correlation potential [48, 41, 49, 50, 51]. The relativistic
corrections to the instantaneous electron–electron Coulomb interactions are
here neglected. If SOC is treated self-consistently, expensive KS Hamiltonian
matrix constructions must be carried out in the multicomponent regime. As a
consequence, most internal parts of quantum mechanical codes must be
reconsidered and adapted. A transparent theoretical formalism that isolates
the SOC terms, reduces computational cost, and simplifies the implementation
of fully relativistic theories can be achieved using quaternion algebra [52,
53, 54, 55].

The most commonly applied strategies incorporate SOC in relativistic
pseudopotentials and the related projector-augmented wave methods [56, 57, 58,
59, 60], where the oscillating wave function with complicated nodal structure
in the core region close to the nuclei is replaced by a smooth pseudo-
wavefunction. However, all-electron approaches are necessary to describe
Rashba-like spin-splitting induced by distortions of the wave function close
to the nucleus [61], where nuclear spins interact with the electrons [62], or
in situations where high accuracy is desired [63, 64]. Linearized augmented
plane-wave (LAPW) [65, 66, 67] and full-potential linear muffin-tin orbital
(LMTO) [68, 69, 70] methods enable all-electron calculations by constructing
muffin-tin spheres around all atoms and expressing the all-electron wave
functions using orbitals inside these spheres and plane waves in the
interstitial region. In such cases, treatment of SOC must be handled
separately in the two distinct regions and is sometimes neglected outside of
the muffin-tin spheres. An alternative all-electron approach to LAPW and LMTO
can be formulated using a finite-element basis set [71] or by expanding the
Bloch wave functions $\psi_{n}(\bm{k};\bm{r})$ using linear combinations of
atom-centered real-space basis functions $\chi_{\mu}(\bm{r})$. Techniques that
employ various choices for $\chi_{\mu}$, such as the numerical [72, 73, 74],
Slater-type [75, 76], or Gaussian-type basis set [77, 54], have been developed
for periodic systems at the 4c [72, 73, 54] as well as the approximate
relativistic [77, 75, 76, 74] levels of the theory.

In the 4c theory, the basis functions $\chi_{\mu}(\bm{r})$ need to represent
Dirac bispinors, and the basis expansion of the 4c Bloch wave functions
$\psi_{n}(\bm{k};\bm{r})$ takes the form:

$\psi_{n}(\bm{k};\bm{r})=\frac{1}{\sqrt{|\mathcal{K}|}}\sum_{\bm{R},\mu}e^{i\bm{k}\cdot{\bm{R}}}\chi_{\mu}(\bm{r}-\bm{R})c^{\mu}_{n}(\bm{k}),$
(2)

where $n$ is the band index, $\bm{k}$ is the reciprocal-space vector (quasi-
momentum), $\bm{R}$ denotes the Bravais lattice vector of the respective unit
cell in the translationally invariant system ($\bm{R}=\bm{0}$ is the reference
unit cell), index $\mu$ runs over the scalar basis functions in the unit cell,
and the normalization constant is the inverse square root of the volume of the
primitive reciprocal unit cell $|\mathcal{K}|$. The expansion coefficients
$c^{\mu}_{n}(\bm{k})$ and the ground-state electronic structure are obtained
by solving the matrix form of the reciprocal-space Dirac–Kohn–Sham (DKS)
generalized eigenvalue equation

$H(\bm{k})c(\bm{k})=S(\bm{k})c(\bm{k})\varepsilon(\bm{k}),$ (3)

where $\varepsilon(\bm{k})$ is the diagonal matrix of eigenvalues (band
energies), and $H(\bm{k})$ and $S(\bm{k})$ are the reciprocal-space DKS
Hamiltonian and overlap matrices, respectively, with the elements:

$\displaystyle H_{\mu\mu^{\prime}}(\bm{k})$
$\displaystyle=\sum_{\bm{R}}e^{i\bm{k}\cdot{\bm{R}}}\int_{\mathbb{R}^{3}}\chi^{\dagger}_{\mu}(\bm{r})\hat{H}\chi_{\mu^{\prime}}(\bm{r}-\bm{R})d^{3}\bm{r},$
(4a) $\displaystyle S_{\mu\mu^{\prime}}(\bm{k})$
$\displaystyle=\sum_{\bm{R}}e^{i\bm{k}\cdot{\bm{R}}}\int_{\mathbb{R}^{3}}\chi^{\dagger}_{\mu}(\bm{r})\chi_{\mu^{\prime}}(\bm{r}-\bm{R})d^{3}\bm{r}.$
(4b)

Due to the dependence of $H(\bm{k})$ on the electron density and its gradient
that are determined from $c(\bm{k})$, Eq. (3) must be solved self-
consistently, though the SOC terms are often neglected in the self-consistency
procedure and only included _a posteriori_ as a perturbation to save
computational cost [78].

Our goal here is to provide an all-electron approach for modeling heavy-
element-containing 2D materials based on the self-consistent treatment of SOC
and Gaussian-type orbitals (GTOs) commonly used in quantum chemistry
calculations of properties, electronic structures, and electron correlations
in molecules. The use of GTOs for solids offers several advantages. For
instance, a unified representation of wave functions in a variety of systems
(molecules, polymers, 2D materials, crystals) allows for building on the
existing quantum chemistry approaches, such as the accurate treatment of
correlation [79, 80, 81, 82], reduced cost of evaluating the exact exchange
[83], and efficient algorithms that scale linearly with system size [84, 85].
In 2D materials that are only a few atomic layers thick, Bloch functions are
straightforwardly constructed to satisfy the Bloch theorem only across the two
periodic dimensions by explicitly restricting $\bm{R}$ and $\bm{k}$ in Eq. (2)
to 2D lattice vectors. As a consequence, the system is not artificially
replicated in the nonperiodic dimension. For comparison with calculations
employing plane waves, this situation would correspond to the limit of
infinitely large vacuum layers (or zero hopping) between the replicated
images. Finally, explicit self-consistent treatment of SOC for all electronic
states on an equal footing without adopting the muffin-tin or pseudopotential
approximations provides a pathway for modeling x-ray absorption spectroscopy
[86] and magnetic response properties associated with nuclear spins [87, 88].

Modelling solid-state systems with GTOs in the nonrelativistic framework was
pioneered by Pisani and Dovesi [89, 90, 91], later joined by several other
groups [84, 85, 92]. However, standard GTO basis sets are constructed by
optimizing calculations on atoms and thus contain basis functions with small
exponents to describe the asymptotic behavior of the atomic wave functions.
These diffuse functions severely hamper the application of GTOs to solids and
cause numerical instabilities [93, 84, 94, 95]. The numerical issues can be
circumvented by removing the most diffuse basis functions at the start of the
calculation [96, 97, 98] or by constructing system-specific basis sets using
reoptimized Gaussian exponents and contraction coefficients [94, 99]. The
former strategy risks producing low quality results [54, 100] whereas the
latter strategy sacrifices the transferability of the so-constructed basis
sets and requires availability of advanced code features (basis set
optimization). Pre-optimized all-electron basis sets are also not available
for heavy elements [94]. The importance of numerical stability is even more
pronounced in relativistic theories where the additional wave-function
components increase the variational freedom, and the restricted kinetic
balance (RKB) condition [101] must be satisfied to prevent the collapse of the
spectral gap between the positive- and negative-energy states.

In this paper, we demonstrate that accurate and converged results can be
obtained for 2D TMD monolayers in various structural phases using commonly
available all-electron valence triple-$\zeta$ basis sets in the fully
relativistic framework without the need for modifying or reoptimizing the
basis functions. Our approach builds on the quaternion algebra-based theory
[54] implemented here with careful numerical considerations to ensure
robustness of the implementation and the resolution-of-the-identity (RI)
approximation for the Coulomb four-center integrals [102, 103, 104] that
reduces the computational cost of otherwise time-consuming calculations by
more than three orders of magnitude. We show how the canonical
orthogonalization can be used to construct numerically well-behaved
orthonormal bases in momentum space in the 4c setting and derive a criterion
for assessing the quality of the RKB condition in such bases. We discuss a
parallel algorithm for the electronic structure solver exhibiting minimal
input/output (IO) disk communication and memory overhead that is applicable to
unit cells with several thousand basis functions and thousands of $\bm{k}$
points. The method presented here is implemented in the ReSpect code [55] and
used to calculate the band structures of the hexagonal (2H), tetragonal (1T),
and distorted tetragonal (1T’) phases of selected TMDs. For the 1T’ phase, we
evaluate the $\mathbb{Z}_{2}$ invariant within our real-space GTO-based scheme
and confirm the findings of Qian _et al._ [11] that the 1T’ phase of $MX_{2}$
is topologically nontrivial.

The data presented in this study as well as the scripts used for pre- and
post-processing of the input and output files are available in the ZENODO
public repository [105]. The ReSpect code used in this study is available upon
a reasonable request free of charge, see Ref. [106].

##  II Methods

###  II.1 Construction of orthonormal momentum-space basis

Throughout this section, we assume Einstein’s implicit summation over doubly
repeating indices.

In order to avoid variational collapse associated with an incomplete basis
representation of the lower components of the wave function [107] and to
obtain the correct nonrelativistic limit of the kinetic energy operator in a
finite basis, we impose the RKB condition [101] for the basis bispinors
$\chi_{\mu}(\bm{r})$ in real space. Hence, we require that [108, 49, 109, 53,
55, 43]

$\chi_{\mu}(\bm{r})\equiv\begin{pmatrix}\mathbb{I}_{2}&0_{2}\\\
0_{2}&\frac{1}{2c}\bm{\sigma}\cdot\bm{p}\end{pmatrix}g_{\mu}(\bm{r}-\bm{A}_{\mu}),$
(5)

where $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix,
$\bm{p}\equiv-i\bm{\nabla}$ is the electron momentum operator, and $g_{\mu}$
are scalar basis functions centered on atomic positions $\bm{A}_{\mu}$. Due to
the multicomponent structure of $\hat{H}$ and $\chi_{\mu}$,
$H_{\mu\mu^{\prime}}$ and $S_{\mu\mu^{\prime}}$ in Eq. (4) are $4\times 4$
complex matrices for each pair $\mu,\mu^{\prime}$ and for each $\bm{k}$. For
the functions $g_{\mu}$, we choose the primitive spherical GTOs

$g_{\mu}(\bm{r})\equiv\mathcal{N}Y_{lm}(\vartheta,\varphi)e^{-\alpha\bm{r}^{2}},$
(6)

where $\mathcal{N}$ is the normalization constant, $\alpha$ is the Gaussian
exponent, and $Y_{lm}(\vartheta,\varphi)$ are the spherical harmonics. The GTO
basis is commonly implemented in many quantum chemistry codes [91, 110, 111,
112, 113, 114, 92]. For computational reasons, integrals in Eqs. (4) are
evaluated analytically by solving recurrence relations [115] formulated in
terms of the Cartesian GTOs [116]. The resulting integrals are then
transformed to spherical GTOs.

To construct an orthonormal basis in momentum space, let us first define a
nonorthogonal basis from
$\chi_{\mu\bm{R}}(\bm{r})\equiv\chi_{\mu}(\bm{r}-\bm{R})$ as

$\chi_{\mu}(\bm{k})=\frac{1}{\sqrt{|\mathcal{K}|}}\sum_{\bm{R}}e^{i\bm{k}\cdot{\bm{R}}}\chi_{\mu\bm{R}}.$
(7)

A new basis is obtained using the transformation:

$\varphi_{p}(\bm{k})=\chi_{\mu}(\bm{k})B^{\mu}_{p}(\bm{k}).$ (8)

From here, for clarity, we drop the dependence of all matrices on $\bm{k}$.
The matrix $B$ can be chosen so that $\varphi_{p}$ are orthonormal, _i.e._
$B^{\dagger}SB=\mathbb{I}$. We perform the canonical orthogonalization,

$B=Us^{-1/2}\Lambda,$ (9)

where $U$ is the unitary matrix of eigenvectors and $s$ is the diagonal matrix
of eigenvalues, both obtained by diagonalizing the overlap matrix $S$.
$\Lambda$ is a rectangular matrix consisting of the square identity matrix and
zero rows that correspond to the basis functions that are projected out. In
case no projections are required, $\Lambda$ is simply the identity matrix.

In the relativistic 4c theory, $\chi_{\mu}$ take the bispinor form of Eq. (5),
and the $B^{\mu}_{p}$ become $4\times 4$ matrices for each $\mu,p$. However,
let us now consider a general basis set of bispinors $\chi_{\mu}$ with
$\chi^{L}_{\mu}$ and $\chi^{S}_{\mu}$ denoting the $2\times 2$ large (upper)
and small (lower) components of the basis, respectively. Similarly, let
$\varphi^{L}_{p}\equiv\chi^{L}_{\mu}X^{\mu}_{p}$ and
$\varphi^{S}_{p}\equiv\chi^{S}_{\mu}Y^{\mu}_{p}$ denote the transformed
(_e.g._ orthonormalized) large and small components of the basis,
respectively, $X$ and $Y$ forming a block-diagonal matrix $B$ in Eq. (8). The
DKS Hamiltonian takes the matrix form:

$H=\begin{pmatrix}V_{LL}&c\Pi_{LS}\\\
c\Pi_{SL}&V_{SS}-2c^{2}S_{SS}\end{pmatrix},$ (10)

where $V_{LL},V_{SS},\Pi_{LS},\Pi_{SL}$ are the matrix representations of
their respective operators in the $\chi_{\mu}$ basis, and $S_{SS}$ is the
overlap matrix for the small-component basis. Ensuring numerical stability
requires projecting out basis functions $\varphi^{L}_{p}$ and
$\varphi^{S}_{p}$ corresponding to the smallest eigenvalues of the 4c overlap
matrix. At the same time, removing $\varphi^{S}_{p}$ functions degrades the
RKB condition that is no longer satisfied exactly in the new basis, which can
lead to the emergence of artificial in-gap states. Here, we proceed by showing
how we numerically track this basis truncation error. First, we express the
RKB condition in the matrix form using the arbitrary basis $\chi_{\mu}$. If we
write the 4c eigenvalue problem with the Hamiltonian matrix in Eq. (10) and
eliminate the small-component wave-function coefficients, we obtain

$\left[V_{LL}+\frac{1}{2}\Pi_{LS}\mathcal{B}^{-1}_{SS}(\varepsilon)\Pi_{SL}\right]c_{L}=\varepsilon
S_{LL}c_{L},$ (11)

where

$\mathcal{B}_{SS}(\varepsilon)=S_{SS}+\frac{\varepsilon
S_{SS}-V_{SS}}{2c^{2}},$ (12)

$c_{L}$ are the coefficients of the large component of the wave function, and
$\varepsilon$ is the one-electron energy. Imposing the requirement that the
nonrelativistic Hamiltonian is recovered in the limit $c\rightarrow\infty$
gives the following condition:

$T_{LL}=\frac{1}{2}\Pi_{LS}S^{-1}_{SS}\Pi_{SL},$ (13)

where $T_{LL}$ is the nonrelativistic kinetic energy matrix in the basis of
the large component. We note that the choice of bispinors $\chi_{\mu}$, as in
Eq. (5), satisfies Eq. (13) analytically, since
$\Pi_{LS}=\Pi_{SL}=\frac{1}{c}T_{LL}$ and $S_{SS}=\frac{1}{2c^{2}}T_{LL}$. In
the transformed basis $\varphi_{p}$, Eq. (13) reads:

$X^{\dagger}T_{LL}X=\frac{1}{2}X^{\dagger}\Pi_{LS}Y\left(Y^{\dagger}S_{SS}Y\right)^{-1}Y^{\dagger}\Pi_{SL}X.$
(14)

This RKB condition remains valid as long as two conditions are met: The former
basis $\chi_{\mu}$ satisfies Eq. (13), and the inverse $Y^{-1}$ exists, _i.e._
the new basis spans the same space as the original one. By requiring that the
new basis is orthonormal ($Y^{\dagger}S_{SS}Y=\mathbb{I}$), this condition
simplifies to:

$X^{\dagger}T_{LL}X=\frac{1}{2}X^{\dagger}\Pi_{LS}YY^{\dagger}\Pi_{SL}X,$ (15)

and the existence of $Y^{-1}$ implies that $S^{-1}_{SS}=YY^{\dagger}$.
However, orthonormalization procedures that involve removing eigenvectors in
the orthonormal basis of $\varphi^{S}_{p}$ by using a rectangular $\Lambda$ in
Eq. (9) do not generally preserve the RKB condition. Hence, in the 4c
framework, improving the conditioning of matrices by projecting out redundant
basis functions can sacrifice the RKB between the $L$ and $S$ components.

To assess the quality of the orthonormal basis for the lower components of the
wave function, we calculate how much of the RKB condition is lost by the
transformation to the reduced orthonormal basis. To this end, we evaluate

$e=\max\left[X^{\dagger}\left(\frac{1}{2c^{2}}T_{LL}YY^{\dagger}T_{LL}-T_{LL}\right)X\right],$
(16)

for each $\bm{k}$. Here, $\max$ indicates the largest matrix element. We found
that large values of $e$ ($>10^{-5}$) could be associated with the emergence
of artificial states in the band structure. To keep the error in the RKB
small, for every $\varphi^{S}_{p}$ function that is projected out, we also
remove the corresponding $\varphi^{L}_{p}$ function regardless of how large
the eigenvalue of the overlap matrix $S_{LL}$ is for this function. This means
that the $\varphi^{L}_{p}$ must be removed even if it is not causing
overcompleteness of the large-component subspace. We note that these measures
are necessary for moderately sized (triple-$\zeta$ and quadruple-$\zeta$)
basis sets, but the numerical issues are not observed in case the smaller
double-$\zeta$ basis are used. Finally, the process of orthonormalization and
removal of redundant functions heavily depends on $\bm{k}$, so the
implementation of the electronic structure solver must have the flexibility to
account for matrices with different sizes for each $\bm{k}$.

###  II.2 Parallel electronic structure solver

Construction of the ground-state one-electron wave functions within the self-
consistent field (SCF) theory consists of two distinct steps that are iterated
repeatedly until self-consistency is reached. The DKS Hamiltonian and overlap
matrices in Eqs. (4) are first assembled in real space using the GTO basis.
These are then transformed to reciprocal space, and Eq. (3) is solved in the
orthonormal basis using matrix algebra (multiplications and diagonalizations).
Large-scale or high-throughput calculations of band structures of solids
require efficient algorithms that offer good (ideally linear) scaling with
respect to the number of computer cores used as well as the size of the unit
cell. This must be accomplished for both of the SCF solver steps. The demand
for efficiency is especially important for all-electron 4c calculations, where
the matrix sizes are much larger and even otherwise negligible matrix
operations become time-consuming.

Evaluation of the DKS Hamiltonian matrix elements is dominated by the
electron–electron Coulomb contributions for DFT simulations with pure
exchange–correlation functionals. We adopt the resolution-of-the-identity (RI)
approach [102, 103, 104] (also known as the density fitting procedure)
combined with multipole expansions [117, 85] to significantly reduce the
overall cost of computing the Coulomb terms. The RI method approximates the
electron density by a linear combination of auxiliary basis functions centered
on atoms, as opposed to the exact treatment of the four-center Coulomb
integrals where the orbital products that constitute the density are expanded
individually. Our RI implementation differs from earlier work [85, 118]. The
theoretical challenges associated with the divergent terms in the Coulomb
metric matrix of periodic systems together with details of our RI method will
be presented elsewhere.

Once the real-space matrix elements are calculated, the remaining steps –
transformation to reciprocal space, orthonormalization procedure, matrix
diagonalization – is performed independently for each $\bm{k}$. We exploit
factorization of the reciprocal-space tasks by employing message-passing-
interface (MPI) directives which allow high level of parallelization to be
achieved due to almost no communication needed among various MPI processes.
Each process utilizes multiple OpenMP threads when calling internally
parallelized matrix libraries. However, in order to construct the electron
density,

$\rho(\bm{r})=\sum_{n}\int_{\mathcal{K}}f_{n}(\bm{k})\psi^{\dagger}_{n}(\bm{k};\bm{r})\psi_{n}(\bm{k};\bm{r})d^{3}\bm{k},$
(17)

where $f_{n}(\bm{k})$ is the occupation number of the $n$-th band and
$\psi_{n}(\bm{k};\bm{r})$ are the Bloch functions defined in Eq. (2). Band
energies $\varepsilon_{n}(\bm{k})$, from which the occupations are determined,
must be known for all $\bm{k}$. Hence, the evaluation of the electron density
(needed for subsequent SCF iterations) occurs after all the wave-function
coefficients $c^{\mu}_{n}(\bm{k})$ and energies are found. In the all-electron
4c framework, the coefficient matrices are too large to be stored in memory
for all $\bm{k}$ and storing the matrices on disk hampers the MPI
parallelization due to significantly increased I/O communication. To this end,
we partition the electron density into two terms:

$\rho(\bm{r})=\rho_{\text{c}}(\bm{r})+\rho_{\text{v}}(\bm{r}),$ (18)

where, for $\rho_{\text{c}}$, we restrict the sum $\sum_{n}$ in Eq. (17) to
contain only those bands that are assumed to be fully filled
($f_{n}(\bm{k})=1$) during the entire SCF procedure. $\rho_{\text{c}}$
contains contributions to the electron density from the majority of bands –
these terms are evaluated on-the-fly after diagonalization independently for
each $\bm{k}$ point without the need to store the full coefficient matrices.
The remaining terms in the density are incorporated in $\rho_{\text{v}}$ and
calculated in a separate loop over $\bm{k}$ after all
$\varepsilon_{n}(\bm{k})$ and $f_{n}(\bm{k})$ are known, which allows for a
small number of partially filled bands [$f_{n}(\bm{k})$ to be varied across
$\bm{k}$]. The coefficients $c^{\mu}_{n}(\bm{k})$ needed for $\rho_{\text{v}}$
can be stored in memory as they form narrow rectangular matrices with only one
dimension increasing with the system size. In addition, each MPI process only
needs to keep the coefficients for the subset of $\bm{k}$ points that the
process handles, which is advantageous for multinode calculations within the
distributed-memory architectures since the coefficients do not need to be
shared or communicated among nodes. The parallel scheme described here enables
seamless all-electron fully relativistic 4c calculations of band structures
and density-of-states with a large number of $\bm{k}$ points using thousands
of cores.

###  II.3 Evaluation of $\mathbb{Z}_{2}$ invariant

The QSH phase of 2D materials with time-reversal symmetry is characterized by
a nontrivial topological order associated with a nonzero $\mathbb{Z}_{2}$
index [119]. Identification of topological materials from numerical
simulations is not straightforward [120]. However, for systems with inversion
symmetry, the calculation of the $\mathbb{Z}_{2}$ invariant can be simplified
using Fu and Kane’s method [121] based on the knowledge of the parity of the
Bloch wave functions of the bulk crystal at the time-reversal invariant
momenta (TRIM) of the Brillouin zone. In this section, we describe how the
$\mathbb{Z}_{2}$ invariant is calculated in the 4c real-space GTO basis.

First, we discuss the construction of the momentum-space parity operator
$P_{\mu\mu^{\prime}}(\bm{k})$. Let $\mathcal{I}$ denote the space-inversion
operator with the inversion center at $\bm{G}$, defined by its action on a
scalar function $f$ as $\mathcal{I}f(\bm{r})\equiv f(2\bm{G}-\bm{r})$. From
the requirement that the parity-transformed bispinor wave function
$\psi(\bm{r},t)$ satisfies the same Dirac equation, it is possible to identify
the form of the 4c parity operator $\hat{P}$ as:

$\hat{P}\psi(\bm{r},t)=\eta\beta\psi(2\bm{G}-\bm{r},t)\equiv\eta\beta\mathcal{I}\psi(\bm{r},t),$
(19)

where $\beta=\operatorname{diag}(\mathbb{I}_{2},-\mathbb{I}_{2})$ and $\eta$
is an arbitrary phase factor ($|\eta|=1$) that can be introduced by the parity
transformation. Neglecting the phase gives $\hat{P}=\beta\mathcal{I}$. The 4c
matrix form of the parity operator is obtained by realizing that
$\mathcal{I}\bm{p}\mathcal{I}^{\dagger}=-\bm{p}$ and letting $\hat{P}$ act on
the RKB basis in Eq. (5), _i.e._

$\hat{P}\chi_{\mu}\equiv\hat{P}\begin{pmatrix}\mathbb{I}_{2}&0_{2}\\\
0_{2}&\frac{1}{2c}\bm{\sigma}\cdot\bm{p}\end{pmatrix}g_{\mu}=\begin{pmatrix}\mathbb{I}_{2}&0_{2}\\\
0_{2}&\frac{1}{2c}\bm{\sigma}\cdot\bm{p}\end{pmatrix}\mathcal{I}g_{\mu}.$ (20)

It follows that the parity matrix $P_{\mu\mu^{\prime}}(\bm{k})$ can be
obtained using Eq. (4b) for the overlap matrix with an additional application
of the inversion $\mathcal{I}$ on the scalar basis function
$g_{\mu^{\prime}}$. For an orbital from Eq. (6) with angular momentum number
$l$ that is centered at $\bm{A}$ (which here denotes the position of an atom
$\bm{A}_{\mu}$ in a unit cell $\bm{R}$), this inversion gives:

$\mathcal{I}g_{\mu}(\bm{r}-\bm{A})=(-1)^{l}g_{\mu}\left(\bm{r}+\bm{A}-2\bm{G}\right).$
(21)

Finally, the $\mathbb{Z}_{2}$ invariant is calculated using the eigenvalues
$\xi_{2m}$ of the parity $P_{\mu\mu^{\prime}}(\bm{k})$ as [121]:

$(-1)^{\nu}=\prod_{i=1}^{4}\prod_{m=1}^{N}\xi_{2m}(\Gamma_{i}),$ (22)

where $N$ is the number of occupied Kramers pairs, and $\Gamma_{i}$ labels the
four TRIM points in 2D.

We conclude this section by noting a few considerations required for a
numerically stable and reliable calculation of the parity matrix and its
eigenvalues when using the method described above. First, the sparse matrix
storage scheme outlined in Sec. II.1 based on the locality of basis function
products $\chi^{\dagger}_{\mu}(\bm{r})\chi_{\mu^{\prime}}(\bm{r}-\bm{R})$
needs to be modified to respect the nonlocal nature of the inversion operator.
Specifically, the list of significant elements of the inversion matrix is not
the same as for the overlap matrix but it is related via a reflection.
Likewise, it must be ensured that the atomic basis set chosen for the
inversion-equivalent atoms is the same. Finally, it is not automatically
guaranteed that the Bloch functions in Eq. (2) that are solutions of the DKS
equation are also eigenfunctions of the parity operator for Bloch states with
degenerate energies. Parity eigenvalues can be obtained by diagonalizing the
sub-blocks of the parity matrix that correspond to these degenerate energy
levels. However, in numerical simulations, energy values obtained from
diagonalization routines are not necessarily exactly degenerate, but rather
appear as near-degenerate. At the same time, eigenvectors corresponding to
these near-degenerate levels are not uniquely defined, so an energy threshold
must be introduced to identify the degeneracies.

##  III Results and Discussion

We calculated the electronic band structures of six TMDs (MoS2, MoSe2, MoTe2,
WS2, WSe2, WTe2) in the 2H, 1T, and 1T’ structural phases using our 4c all-
electron GTO method. We compared our 2H and 1T results with the corresponding
band structures presented in Ref. [6] obtained from the method of Te Velde and
Baerends [122]. For the 1T’ phase, we used the Vasp program package [123] to
investigate how our all-electron method performs compared to the
pseudopotential approximation. Due to the presence of heavy atoms, all TMDs
considered exhibit significant scalar relativistic and SOC effects. We
investigated three different levels of theory: nonrelativistic (nr) with
infinite speed of light, Dirac-based 4c scalar relativistic (sr) without SOC,
and the Dirac-based 4c fully relativistic (fr) with SOC included. In all our
calculations, we employed the nonrelativistic GGA-type XC functional PBE
[124]. For the 4c calculations, the Gaussian finite nucleus model of Visscher
and Dyall [125] was used instead of the point nuclei in order to regularize
the singularity of the lower components of the wave function evaluated at the
atomic centers. All electronic structure optimizations were performed using
the relaxed unit cell geometries taken from Refs. [6] (2H and 1T phases) and
[11] (1T’ phase). For the momentum-space Brillouin zone integrations, we used
the $\Gamma$-centered mesh of $\bm{k}$ points of $11\times 11$ for the 2H and
1T’ phases and $33\times 33$ for the metallic 1T phase. When comparing the
band structures $\varepsilon^{(1)}_{n}(\bm{k})$ and
$\varepsilon^{(2)}_{n}(\bm{k})$ obtained by using two different methods or
computational settings, we evaluate the maximum difference between the energy
eigenvalues as well as the root-mean-square deviation (“band delta”) [78],

$\Delta_{\text{b}}(\mathcal{W})=\sqrt{\frac{1}{N_{E}}\sum_{\begin{subarray}{c}\bm{k},n\\\
\varepsilon^{(1)}_{n}(\bm{k})\in\mathcal{W}\\\
\varepsilon^{(2)}_{n}(\bm{k})\in\mathcal{W}\end{subarray}}\left(\varepsilon^{(1)}_{n}(\bm{k})-\varepsilon^{(2)}_{n}(\bm{k})\right)^{2}},$
(23)

where $\mathcal{W}=[\varepsilon_{\text{l}},\varepsilon_{\text{u}}]$ is an
energy window chosen for comparison, the summations run over all energy states
along a path of $\bm{k}$ points for which both $\varepsilon^{(1)}_{n}(\bm{k})$
and $\varepsilon^{(2)}_{n}(\bm{k})$ lie inside the window $\mathcal{W}$, and
$N_{E}$ is the total number of such states.

###  III.1 Numerical and basis-set accuracy

Numerical instabilities in simulations involving periodic systems with GTOs
originate from two distinct, albeit related, reasons [93, 84]. First, the
numerical accuracy of the KS Hamiltonian matrix can be insufficient when
various approximations necessary for solid-state simulations are introduced
for basis sets with diffuse functions. Diffuse basis functions can extend over
a very large number of unit cells, and the overlaps between such functions
exhibit a much slower decay than their more localized counterparts. If these
overlaps are neglected too early in calculations, the summations in Eqs. (4)
are not converged, causing the overlap matrix $S_{\mu\mu^{\prime}}(\bm{k})$ to
be indefinite, _i.e._ having negative eigenvalues. This problem can be
mitigated by including a large number of basis function products
$\chi^{\dagger}_{\mu}(\bm{r})\chi_{\mu^{\prime}}(\bm{r}-\bm{R})$, which in
turn increases the memory requirements for storage of matrices. To this end,
in the periodic module of ReSpect, we implemented sparse data structures in
real space in combination with quaternion algebra to efficiently handle
spin–orbit-coupled wave-function components. This versatile approach allows us
to retain the accuracy necessary for numerical robustness with a very small
additional memory footprint.

The second problem is associated with overcompleteness (linear dependencies)
of the basis that causes matrices in Eq. (3) to be ill-conditioned and
introduces errors and instabilities in matrix operations. This occurs when the
eigenvalues of the momentum-space overlap matrix drop below $10^{-7}$, which
is often the case for larger basis sets. For instance, the smallest eigenvalue
of the overlap matrix of WTe2 in the 1T’ phase is $10^{-11}$, while the
remaining eigenvalues span over 12 orders of magnitude. As a result, very
large elements $\sim 10^{8}$ accumulate in the density matrix, which prohibits
the successful completion of the ground-state convergence procedure. To avoid
such errors, we construct an orthonormal basis so that the redundant basis
functions are projected out. However, we found that truncating the space
spanned by orthonormal basis functions used for describing the lower wave-
function components leads to the formation of artificial states inside the
band gap for some $\bm{k}$ points in the band structure. This can be
understood by realizing that the RKB condition is violated by the truncation
process. To numerically control how much of the RKB is lost, we derived the
matrix form of the RKB condition in a general basis and the expression for the
truncation error. More details on the process of constructing the orthonormal
basis in the 4c framework is described in Sec. II.1.

For all atoms, we employed the relativistic triple-$\zeta$ (TZ) basis sets
developed by Dyall for $d$ [126, 127, 128] and $p$ elements [129, 130]. The
basis sets used were uncontracted and contained high-angular momentum
correlating functions for several outer shells and functions for dipole
polarization of the valence shells. Our initial tests using the smaller
double-$\zeta$ (DZ) basis from the same family yielded results with
inconsistent quality across $\bm{k}$ points – for instance, the band energy at
the $M$ point of 2H MoS2 was improved by 20 meV when the TZ basis was used,
while for other $\bm{k}$ points, this difference was only about 6 meV.
Similarly, the differences in the band structures of 2H WS2 (including the
band gap) between the DZ and TZ basis sets were 20-30 meV. Hence, to ensure
that our results are well-converged for all $\bm{k}$ points and that the
target accuracy is sufficient for reliable comparisons with other methods, we
chose the TZ basis in all calculations reported here.

The basis sets employed contained very diffuse functions with several
exponents of $<10^{-2}$ a.u. The extent of the most diffuse function on Mo and
W was 18.2 and 17.2 Å, respectively, causing the basis function to span more
than a hundred unit cells of the 2D film. This in turn results in several
hundred thousand charge distributions interacting with themselves and their
periodic replicas when the four-center Coulomb integrals are evaluated. Such
calculations are currently unfeasible on common supercomputers as they can
consume a notable portion of allocated resources. However, within the RI
approximation, the charge distributions only interact with the electron
density that is expanded using an atom-centered auxiliary basis. Since the
number of the three-center terms needed is smaller than the number of four-
center terms typically by three orders of magnitude, the evaluation of the
Coulomb operator becomes affordable. Tests that we performed on various
systems including selected 2H and 1T’ structures studied here in both the nr
and fr settings indicate that the maximum error introduced by the RI
approximation is very small—only a few tens of $\mu$eV.

Our implementation allowed us to explore the validity of removing diffuse
functions from the GTO basis sets before starting the electronic-structure
optimization — an approach that is commonly employed for accelerating and
numerically stabilizing calculations of periodic systems using GTOs [96, 97,
98]. From the TZ basis sets, we deleted all exponents $<0.1$ which for most
atoms meant removing one to two functions in each of the $s$, $p$, $d$, and
$f$ shells. Even though this modification reduced the computational cost by a
factor between five and six in case of the nr theory and a factor of four in
case of the 4c theory, we observed lowered accuracy of the eigenenergies of
some bands and $\bm{k}$ points. For instance, the difference between the
energy of the $\Gamma$ point of the valence band of 2H WTe2 obtained from the
calculations using the truncated and full basis sets was 0.1 eV. The virtual
bands located 2.5 eV above the Fermi level were, in general, poorly described
with the truncated basis set. Hence, our results indicate that removing the
diffuse functions in the start of the calculation is neither justified nor
necessary.

###  III.2 Stability of phases

Figure 1:  Relative total energy of the 1T’ phase of TMD monolayers per
$MX_{2}$ calculated with respect to the 2H phase. The results were obtained
for nr, sr, and fr Hamiltonians and show the relativistic effects on the
structural stability.

Before inspecting individual structural phases and validating our approach, we
compare the total energies of each phase. Our results confirm the findings of
Qian _et al._ [11] — for all systems except WTe2, the 2H phase was found to be
the ground-state structure with the lowest energy, followed by the higher-
energy 1T’ phase. The metallic 1T phase is thermodynamically unstable and the
systems undergo Jahn–Teller distortion, relaxing to the 1T’ phase [11]. In the
case of WTe2, the energy of the 1T’ phase is below the 2H phase, which makes
1T’ the most stable phase.

Our study allows for a uniform treatment of relativistic effects for all atoms
and all electrons, as opposed to the calculations employing pseudopotentials,
which offer less control over how relativity is handled for light and heavy
elements. We found that it is the inclusion of the scalar relativistic effects
that decreases the energy of the 1T’ phase with respect to the 2H phase.
Including the SOC further lowers the energy of the 1T’ phase. For the heaviest
system, WTe2, this energy lowering becomes sufficient to change the ground-
state phase from 2H to 1T’. Figure 1 summarizes our results for all systems at
the nr, sr, and fr levels of theory and shows the relative energy differences
between the 1T’ and 2H phases. Hence, without relativistic effects, the
ground-state structure of all six systems would be 2H.

###  III.3 2H and 1T phases

Table 1:  Band gaps and SOC-induced (Rashba) splitting of the valence band at the $K$ point of TMD monolayers in the 2H phase calculated for nr, sr, and fr Hamiltonians. The values in the round and square brackets are taken from Refs. [6] and [131], respectively. | Band gap (Eg, eV) | Rashba splitting (eV)
---|---|---
| nr | sr | fr | fr
MoS2 | $1.820$ | $1.812$ | $1.734$ | $0.146$
|  | $(1.82)$ | $(1.74)$ | $[0.147]$
MoSe2 | $1.560$ | $1.550$ | $1.442$ | $0.185$
|  | $(1.56)$ | $(1.45)$ | $[0.186]$
MoTe2 | $1.171$ | $1.154$ | $1.018$ | $0.218$
|  | $(1.15)$ | $(1.01)$ |
WS2 | $2.060$ | $1.984$ | $1.628$ | $0.416$
|  | $(1.98)$ | $(1.64)$ | $[0.433]$
WSe2 | $1.741$ | $1.630$ | $1.300$ | $0.460$
|  | $(1.63)$ | $(1.33)$ | $[0.463]$
WTe2 | $1.332$ | $1.194$ | $0.863$ | $0.479$
|  | $(1.18)$ | $(0.87)$ |

The 2H and 1T phases of the TMDs both contain three atoms (one metal and two
chalcogens) in the primitive unit cell. However, the structural differences
between the two phases lead to their distinct electronic properties. Whereas
the 1T phase is space-inversion-symmetric with a metallic electronic
structure, the 2H systems are semiconductors with broken space inversion. In
combination with the strong SOC, this broken inversion symmetry leads to
“giant” spin–orbit-induced Rashba splittings in 2H $MX_{2}$ [29].

To validate the approach presented in this work, we calculated the band
structures at the nr, sr, and fr levels of theory using all-electron bases.
The resulting band diagrams traversing high symmetry $\bm{k}$ points obtained
with and without SOC can be found in Figs. S1 (2H) and S2 (1T) in Ref. [132].
The nonrelativistic results differ significantly from those obtained from the
relativistic theories and are not shown. For the 2H and 1T phases, sr band
structures are available in the 2D materials Atlas of Miró, Audiffred, and
Heine [6] that we used for comparison. We conclude that our results agree very
well at the sr level with the Atlas results.

In the MoTe2 and WTe2 band structures of the 1T phase shown in Fig. S2 [132],
it is possible to see very small SOC-induced “Rashba-like” spin splittings of
the order of tens of meV. However, such splittings should not be present for
the inversion-symmetric 1T phase, and all the bands should be strictly doubly
degenerate. We found that this observation can be attributed to small
deviations (about 0.04 Å) in the unit cell geometry of MoTe2 and WTe2 from the
tetragonal lattice, that is, the 1T structure used does not exhibit exact
inversion symmetry. The geometries of the other systems do not show any
deviations from 1T, and we observed exact double degeneracy of bands. Finally,
we note that including the SOC in the computational framework lifted the
degeneracy at the intersection points of the overlapping valence and
conduction bands of the metallic 1T phase. The values of the SOC-induced
vertical gaps at the $K$ point are reported in Table SI of Ref. [132].

Table 1 contains the values of the nr, sr, and fr band gaps of the 2H phase,
as well as the Rashba splittings of the valence band at the $K$ point. We
attribute the small discrepancies between the results obtained with the method
presented here and the Atlas approach to differences in the methodologies. In
the Atlas work, the relativistic effects are treated with the zeroth-order
regular approximation (ZORA) [133], which is a 2c technique that requires
numerical integration schemes due to the appearance of the potential in the
denominator of the Hamiltonian. In contrast, the Hamiltonian used here is 4c
with most of the integrals evaluated analytically. The basis sets employed
should not cause any significant deviations, since both sets of calculations
were carried out with TZ quality bases with several polarization functions.
Finally, we note that the Atlas results included empirical D3 treatment of
London dispersion interactions [134] that we did not consider in this work.
The fact that no larger differences between the two approaches were observed
even when the dispersion corrections were omitted indicates that the band gaps
in thin monolayers are dominated by covalent bonds, which are well described
without these corrections. However, the dispersion forces are important when
multiple layers of 2D materials are held together by van der Waals forces as,
for example, in bilayers and heterostructures [135]. In the case of WSe2, for
which we observed the largest difference of 30 meV in the fr setting, we
further benchmarked against a calculation with a larger quadruple-$\zeta$ (QZ)
basis set. The band energies improved only marginally, _i.e._ by units of meV.
Hence, we believe that our band structures are well converged with respect to
the basis set as well as the relativistic Hamiltonian and should serve as the
reference results for these systems for the PBE functional.

###  III.4 1T’ phase

Figure 2:  Structure of the 1T’ phase of monolayer TMDs.

Figure 3:  Comparison of band structures of $MX_{2}$ monolayers in the 1T’
phase obtained from the ReSpect (dashed lines) and Vasp (full lines) codes at
the relativistic level of theory including SOC. From the top left to bottom
right: MoS2, MoSe2, MoTe2, WS2, WSe2, and WTe2. The horizontal dashed black
line marks the Fermi level. The path traversing high-symmetry $\bm{k}$-points
in the reciprocal-space unit cell was chosen according to Ref. [136]. Table 2:
Root-mean-square deviation between the band structures of TMD monolayers
obtained from the ReSpect and Vasp codes at the fr level of theory (including
SOC) using Eq. (23). The valence region $\mathcal{W}_{\text{val}}$ ranges from
$-10$ to 0 eV (shifted to the Fermi level) and the conduction region
$\mathcal{W}_{\text{con}}$ ranges from 0 to 5 eV. Band delta (meV)
---
| MoS2 | MoSe2 | MoTe2 | WS2 | WSe2 | WTe2
$\Delta_{\text{b}}(\mathcal{W}_{\text{val}})$ | $12.0$ | $8.4$ | $11.7$ | $34.6$ | $25.6$ | $20.2$
$\Delta_{\text{b}}(\mathcal{W}_{\text{con}})$ | $40.5$ | $28.4$ | $32.8$ | $26.2$ | $23.4$ | $42.0$

The 1T’ phase of $MX_{2}$ is formed from the 1T structure by a spontaneous
lattice distortion in which the unit cell period is doubled along one in-plane
direction, which creates zigzag chains along the perpendicular in-plane
direction [137] (see Fig. 2). The TMDs were shown to host the QSH effect in
the 1T’ structure, _i.e._ the existence of time-reversal-symmetry protected
edge states that enable conduction of electrical currents on the surface of
the material [11]. From a computational viewpoint, the 1T’ phase has not been
studied in the literature to the same extent as the 2H and 1T counterparts.
Here, we calculated the electronic band structures of the 1T’ $MX_{2}$ using
the all-electron 4c method presented in this work as well as the
pseudopotential method as implemented in the Vasp package [123]. For the Vasp
calculations, we used an energy cutoff of 500 eV and the Brillouin zone
$\bm{k}$-point sampling of $11\times 15\times 1$. We also tested a denser mesh
of $\bm{k}$ points without observing significant differences. To avoid
artificial interactions between the periodic images of the $MX_{2}$ layers
that appear when 2D systems are studied using plane waves with periodic
boundary conditions, a vacuum region of 20 Å was applied in the direction
perpendicular to the 2D film.

The band structures obtained from the ReSpect and Vasp codes with SOC are
presented in Fig. 3. The values of the band delta differences between the two
codes are shown in Table 2, evaluated separately for the occupied (valence)
and virtual (conduction) energy regions. In general, the differences are of
the order of tens of meV and are larger in the virtual region. The biggest
discrepancies close to the Fermi level are seen in WS2, while the overall
largest band delta values are found in MoS2 and WTe2. For comparison, we
performed tests of the numerical convergence of the GTO method by tuning
various parameters of the ReSpect code. In particular, we modified the number
of removed momentum-space basis functions in the orthonormal basis (see Sec.
II.1), tightened the lattice-sum convergence thresholds, as well as turned off
the RI approximation of the Coulomb integrals. Such modifications had a
negligible effect on our GTO results; the band delta values between the higher
accuracy calculations and the results shown here reached tens of $\mu$eV.
Bigger differences were observed only when we benchmarked the WS2 TZ
calculation against the QZ basis set. Still, the energies of the bands close
to the Fermi level changed by 1–6 meV, and the values of the band deltas for
the intervals $-10$ to $0$ eV and $0$ to $5$ eV were
$\Delta_{\text{b}}(\mathcal{W}_{\text{val}})\approx 4$ meV and
$\Delta_{\text{b}}(\mathcal{W}_{\text{con}})\approx 10$ meV, respectively. To
conclude, despite the minor improvements to the band structure that are
possible when larger basis sets are employed, our TZ results can be considered
as the reference for benchmarking performance of other approximate methods,
for instance, for studying the accuracy of various pseudopotentials.

We explored the topological nature of the 1T’ phase of the $MX_{2}$ systems by
first verifying that the band-gap opening is due to the SOC. Fig. S3 in Ref.
[132] shows that energy degeneracies occurring at points where the valence
band crosses the conduction band are lifted when SOC is included. A band gap
opens in all systems except MoTe2 and WTe2 that remain metallic despite the
band inversion. The lattice distortion in the 1T’ phase of the 2D TMDs does
not break the inversion symmetry. Hence, we used the parity of the Bloch
states evaluated at the TRIM points to calculate the $\mathbb{Z}_{2}$
invariant [121]. Our results for the parities of the individual Kramers pairs
agree with those obtained by Qian [11], which yields $(-1)^{\nu}=1$ for all
six TMDs and confirms the existence of the surface edge states hallmarking the
quantum spin Hall phase.

##  IV Conclusion and outlook

2D TMDs are attracting high current interest for hosting exotic phenomena that
promise applications in spintronic devices and quantum computing. Theoretical
predictions and understanding of such applications relies on accurate
description of the relativistic effects driven by the presence of heavy atoms
in these materials. Accordingly, we have presented here an all-electron
approach that builds on the most widely used basis sets in quantum chemistry
to obtain a scheme that can accurately capture relativistic effects as well as
the properties of electron wave functions close to the nuclei within a uniform
theoretical framework applicable to both solids and molecules.

GTOs provide a convenient link between the finite and periodic or extended
systems, yet their use for solids has been met with many technical
difficulties that sparked skepticism in the community; it has even been
suggested that GTOs should be avoided altogether for condensed matter
calculations [100]. However, in this work, we show how the common limitations
associated with GTOs can be overcome, and that the convergence with respect to
the basis-set limit is possible even in the fully relativistic 4c setting. To
mitigate numerical and performance issues, our method is based on the
quaternion algebra, linear-scaling data structures, RI approximation of the
Coulomb integrals, orthonormal momentum-space bases that sufficiently preserve
the RKB condition, and a parallel electronic structure solver with minimal
required IO operations. Our in-depth analysis of the 2H, 1T, and 1T’ phases of
2D TMDs indicates that our all-electron 4c results can serve as the reference
for developing pseudopotentials or approximate 2c relativistic techniques. We
confirm the existence of the quantum spin Hall effect in the 1T’ phase of
$MX_{2}$ by evaluating the $\mathbb{Z}_{2}$ invariant within our real-space
formulation. Finally, we attribute the origin of the lower ground-state energy
of WTe2 in the 1T’ compared to the 2H phase to (scalar) relativistic effects.

Despite the good agreement between our band structures for the 2H and 1T
phases and those obtained by using the 2c ZORA method [6], our 4c approach
offers several advantages. The exact Dirac Hamiltonian contains the mass-
velocity term that is missing in the ZORA approximation [43], but it is
important for an accurate description of relativistic effects in heavy
elements and for evaluating properties in the presence of external magnetic
fields [138]. Also, the 4c DKS equation does not contain the potential in the
denominator, so that complicated numerical integration of most integrals
involved can be avoided.

Our study is important for opening a pathway for further developments that
leverage quantum chemistry methods for exploring properties that are strongly
affected by the relativistic effects without sacrificing access to the core
regions of the electron charge and spin densities. For instance, accurate
relativistic first-principles calculations of the spin Hamiltonian parameters
(_e.g._ hyperfine coupling constants, Zeeman interactions, or zero-field
splittings) are needed for determining spin dynamical properties of spin
defects in semiconductors [62] to help find viable qubit hosting materials.
Our all-electron method can also be used for generating accurate
pseudopotentials, particularly for systems that contain heavy elements in high
oxidation states.

###### Acknowledgements.

This work was supported by the Research Council of Norway through its Centres
of Excellence scheme (Grant No. 262695), a Research Grant (Grant No. 315822),
and its Mobility Grant scheme (Grant No. 301864), as well as the use of
computational resources provided by UNINETT Sigma2 – The National
Infrastructure for High Performance Computing and Data Storage in Norway
(Grant No. NN4654K). In addition, this project received funding from the
European Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie Grant Agreement No. 945478 (SASPRO2), and the Slovak Research
and Development Agency (Grant No. APVV-21-0497). The work at Northeastern
University was supported by the US Department of Energy (DOE), Office of
Science, Basic Energy Sciences Grant No. DE-SC0022216 and benefited from
Northeastern University’s Advanced Scientific Computation Center and the
Discovery Cluster and the National Energy Research Scientific Computing Center
through DOE Grant No. DE-AC02-05CH11231.

## References

  * Butler  _et al._ [2013] S. Z. Butler, S. M. Hollen, L. Cao, Y. Cui, J. A. Gupta, H. R. Gutiérrez, T. F. Heinz, S. S. Hong, J. Huang, A. F. Ismach, _et al._ , Progress, challenges, and opportunities in two-dimensional materials beyond graphene, ACS nano 7, 2898 (2013).
  * Han [2016] W. Han, Perspectives for spintronics in 2D materials, APL Mater. 4, 032401 (2016).
  * Frindt [1966] R. Frindt, Single crystals of MoS2 several molecular layers thick, J. Appl. Phys. 37, 1928 (1966).
  * Geim and Novoselov [2007] A. K. Geim and K. S. Novoselov, The rise of graphene, Nat. Mater. 6, 183 (2007).
  * Sarma  _et al._ [2011] S. D. Sarma, S. Adam, E. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Rev. Mod. Phys. 83, 407 (2011).
  * Miró  _et al._ [2014] P. Miró, M. Audiffred, and T. Heine, An atlas of two-dimensional materials, Chem. Soc. Rev. 43, 6537 (2014).
  * Tombros  _et al._ [2007] N. Tombros, C. Jozsa, M. Popinciuc, H. T. Jonkman, and B. J. Van Wees, Electronic spin transport and spin precession in single graphene layers at room temperature, Nature 448, 571 (2007).
  * Li  _et al._ [2013] X. Li, J. T. Mullen, Z. Jin, K. M. Borysenko, M. B. Nardelli, and K. W. Kim, Intrinsic electrical transport properties of monolayer silicene and MoS2 from first principles, Phys. Rev. B 87, 115418 (2013).
  * Kane and Mele [2005a] C. L. Kane and E. J. Mele, Quantum spin Hall effect in graphene, Phys. Rev. Lett. 95, 226801 (2005a).
  * Liu  _et al._ [2011] C.-C. Liu, W. Feng, and Y. Yao, Quantum spin Hall effect in silicene and two-dimensional germanium, Phys. Rev. Lett. 107, 076802 (2011).
  * Qian  _et al._ [2014] X. Qian, J. Liu, L. Fu, and J. Li, Quantum spin Hall effect in two-dimensional transition metal dichalcogenides, Science 346, 1344 (2014).
  * Hor  _et al._ [2009] Y. Hor, A. Richardella, P. Roushan, Y. Xia, J. Checkelsky, A. Yazdani, M. Hasan, N. P. Ong, and R. J. Cava, $p$-type Bi2Se3 for topological insulator and low-temperature thermoelectric applications, Phys. Rev. B 79, 195208 (2009).
  * Wang  _et al._ [2012] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotechnol. 7, 699 (2012).
  * Nilsson  _et al._ [2008] J. Nilsson, A. Akhmerov, and C. Beenakker, Splitting of a Cooper pair by a pair of Majorana bound states, Phys. Rev. Lett. 101, 120403 (2008).
  * Fu and Kane [2009] L. Fu and C. L. Kane, Josephson current and noise at a superconductor/quantum-spin-Hall-insulator/superconductor junction, Phys. Rev. B 79, 161408 (2009).
  * Radisavljevic  _et al._ [2011] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Single-layer MoS2 transistors, Nat. Nanotechnol. 6, 147 (2011).
  * Britnell  _et al._ [2012] L. Britnell, R. Gorbachev, R. Jalil, B. Belle, F. Schedin, A. Mishchenko, T. Georgiou, M. Katsnelson, L. Eaves, S. Morozov, _et al._ , Field-effect tunneling transistor based on vertical graphene heterostructures, Science 335, 947 (2012).
  * Solomon  _et al._ [2015] P. Solomon, B. Bryce, M. Kuroda, R. Keech, S. Shetty, T. Shaw, M. Copel, L.-W. Hung, A. Schrott, C. Armstrong, _et al._ , Pathway to the piezoelectronic transduction logic device, Nano Lett. 15, 2391 (2015).
  * Li  _et al._ [2015] H. Li, A. W. Contryman, X. Qian, S. M. Ardakani, Y. Gong, X. Wang, J. M. Weisse, C. H. Lee, J. Zhao, P. M. Ajayan, _et al._ , Optoelectronic crystal of artificial atoms in strain-textured molybdenum disulphide, Nat. Commun. 6, 1 (2015).
  * Brar  _et al._ [2011] V. W. Brar, R. Decker, H.-M. Solowan, Y. Wang, L. Maserati, K. T. Chan, H. Lee, Ç. O. Girit, A. Zettl, S. G. Louie, _et al._ , Gate-controlled ionization and screening of cobalt adatoms on a graphene surface, Nat. Phys. 7, 43 (2011).
  * Natterer  _et al._ [2012] F. D. Natterer, F. Patthey, and H. Brune, Ring state for single transition metal atoms on boron nitride on Rh (111), Phys. Rev. Lett. 109, 066101 (2012).
  * Tsai  _et al._ [2013] W.-F. Tsai, C.-Y. Huang, T.-R. Chang, H. Lin, H.-T. Jeng, and A. Bansil, Gated silicene as a tunable source of nearly 100% spin-polarized electrons, Nat. Commun. 4, 1 (2013).
  * Ponomarenko  _et al._ [2011] L. Ponomarenko, A. Geim, A. Zhukov, R. Jalil, S. Morozov, K. Novoselov, I. Grigorieva, E. Hill, V. Cheianov, V. Fal’Ko, _et al._ , Tunable metal–insulator transition in double-layer graphene heterostructures, Nat. Phys. 7, 958 (2011).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Van der Waals heterostructures, Nature 499, 419 (2013).
  * Haigh  _et al._ [2012] S. Haigh, A. Gholinia, R. Jalil, S. Romani, L. Britnell, D. Elias, K. Novoselov, L. Ponomarenko, A. Geim, and R. Gorbachev, Cross-sectional imaging of individual layers and buried interfaces of graphene-based heterostructures and superlattices, Nat. Mater. 11, 764 (2012).
  * Manzeli  _et al._ [2017] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2D transition metal dichalcogenides, Nat. Rev. Mater. 2, 1 (2017).
  * Kang  _et al._ [2019] K. Kang, T. Li, E. Sohn, J. Shan, and K. F. Mak, Nonlinear anomalous Hall effect in few-layer WTe2, Nat. Mater. 18, 324 (2019).
  * Wang  _et al._ [2019] Z. Wang, B. J. Wieder, J. Li, B. Yan, and B. A. Bernevig, Higher-Order Topology, Monopole Nodal Lines, and the Origin of Large Fermi Arcs in Transition Metal Dichalcogenides XTe2 (X $=$ Mo, W), Phys. Rev. Lett. 123, 186401 (2019).
  * Zhu  _et al._ [2011] Z. Zhu, Y. Cheng, and U. Schwingenschlögl, Giant spin-orbit-induced spin splitting in two-dimensional transition-metal dichalcogenide semiconductors, Phys. Rev. B 84, 153402 (2011).
  * Wilson  _et al._ [1974] J. Wilson, F. Di Salvo, and S. Mahajan, Charge-density waves in metallic, layered, transition-metal dichalcogenides, Phys. Rev. Lett. 32, 882 (1974).
  * Flicker and Van Wezel [2015] F. Flicker and J. Van Wezel, Charge order from orbital-dependent coupling evidenced by NbSe2, Nat. Commun. 6, 1 (2015).
  * Zhu  _et al._ [2016] X. Zhu, Y. Guo, H. Cheng, J. Dai, X. An, J. Zhao, K. Tian, S. Wei, X. C. Zeng, C. Wu, _et al._ , Signature of coexistence of superconductivity and ferromagnetism in two-dimensional NbSe2 triggered by surface molecular adsorption, Nat. Commun. 7, 1 (2016).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of MoS2 and other group-VI dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012).
  * Mak  _et al._ [2014] K. F. Mak, K. L. McGill, J. Park, and P. L. McEuen, The valley Hall effect in MoS2 transistors, Science 344, 1489 (2014).
  * Pulkin and Yazyev [2016] A. Pulkin and O. V. Yazyev, Spin-and valley-polarized transport across line defects in monolayer MoS2, Phys. Rev. B 93, 041419 (2016).
  * Xu  _et al._ [2014] X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Spin and pseudospins in layered transition metal dichalcogenides, Nat. Phys. 10, 343 (2014).
  * Zhou  _et al._ [2016] B. T. Zhou, N. F. Yuan, H.-L. Jiang, and K. T. Law, Ising superconductivity and Majorana fermions in transition-metal dichalcogenides, Phys. Rev. B 93, 180501 (2016).
  * He  _et al._ [2018] W.-Y. He, B. T. Zhou, J. J. He, N. F. Yuan, T. Zhang, and K. T. Law, Magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides, Commun. Phys. 1, 1 (2018).
  * Tsai  _et al._ [2022] J.-Y. Tsai, J. Pan, H. Lin, A. Bansil, and Q. Yan, Antisite defect qubits in monolayer transition metal dichalcogenides, Nat. Commun. 13, 1 (2022).
  * Dirac [1928] P. A. Dirac, The quantum theory of the electron, Proc. R. Soc. A 117, 610 (1928).
  * Saue and Helgaker [2002] T. Saue and T. Helgaker, Four-component relativistic Kohn–Sham theory, J. Comput. Chem. 23, 814 (2002).
  * Belpassi  _et al._ [2011] L. Belpassi, L. Storchi, H. M. Quiney, and F. Tarantelli, Recent advances and perspectives in four-component Dirac–Kohn–Sham calculations, Phys. Chem. Chem. Phys. 13, 12368 (2011).
  * Saue [2011] T. Saue, Relativistic Hamiltonians for chemistry: a primer, ChemPhysChem 12, 3077 (2011).
  * Liu [2014] W. Liu, Advances in relativistic molecular quantum mechanics, Phys. Rep. 537, 59 (2014).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Yeh  _et al._ [2022] C.-N. Yeh, A. Shee, Q. Sun, E. Gull, and D. Zgid, Relativistic self-consistent GW: Exact two-component formalism with one-electron approximation for solids, Phys. Rev. B 106, 085121 (2022).
  * Rajagopal and Callaway [1973] A. Rajagopal and J. Callaway, Inhomogeneous electron gas, Phys. Rev. B 7, 1912 (1973).
  * Komorovsky  _et al._ [2008] S. Komorovsky, M. Repisky, O. L. Malkina, V. G. Malkin, I. Malkin Ondik, and M. Kaupp, A fully relativistic method for calculation of nuclear magnetic shielding tensors with a restricted magnetically balanced basis in the framework of the matrix Dirac–Kohn–Sham equation, J. Chem. Phys. 128, 104101 (2008).
  * Engel [2002] E. Engel, Relativistic Density Functional Theory: Foundations and Basic Formalism, in  _Theoretical and computational chemistry_, Vol. 11, edited by P. Schwerdtfeger (Elsevier, 2002) Chap. 10, pp. 523–621.
  * Dyall and Faegri [2007] K. G. Dyall and K. Faegri,  _Introduction to relativistic quantum chemistry_ (Oxford University Press, New York, 2007) pp. 268–274.
  * Saue  _et al._ [1997] T. Saue, K. Fægri, T. Helgaker, and O. Gropen, Principles of direct 4-component relativistic SCF: application to caesium auride, Mol. Phys. 91, 937 (1997).
  * Konecny  _et al._ [2018] L. Konecny, M. Kadek, S. Komorovsky, K. Ruud, and M. Repisky, Resolution-of-identity accelerated relativistic two-and four-component electron dynamics approach to chiroptical spectroscopies, J. Chem. Phys. 149, 204104 (2018).
  * Kadek  _et al._ [2019] M. Kadek, M. Repisky, and K. Ruud, All-electron fully relativistic Kohn-Sham theory for solids based on the Dirac-Coulomb Hamiltonian and Gaussian-type functions, Phys. Rev. B 99, 205103 (2019).
  * Repisky  _et al._ [2020] M. Repisky, S. Komorovsky, M. Kadek, L. Konecny, U. Ekström, E. Malkin, M. Kaupp, K. Ruud, O. L. Malkina, and V. G. Malkin, ReSpect: Relativistic spectroscopy DFT program package, J. Chem. Phys. 152, 184101 (2020).
  * Kleinman [1980] L. Kleinman, Relativistic norm-conserving pseudopotential, Phys. Rev. B 21, 2630 (1980).
  * Bachelet and Schlüter [1982] G. B. Bachelet and M. Schlüter, Relativistic norm-conserving pseudopotentials, Phys. Rev. B 25, 2103 (1982).
  * Dal Corso and Conte [2005] A. Dal Corso and A. M. Conte, Spin-orbit coupling with ultrasoft pseudopotentials: Application to Au and Pt, Phys. Rev. B 71, 115106 (2005).
  * Nichols  _et al._ [2009] P. Nichols, N. Govind, E. J. Bylaska, and W. A. De Jong, Gaussian Basis Set and Planewave Relativistic Spin-Orbit Methods in NWChem, J. Chem. Theory Comput. 5, 491 (2009).
  * Dal Corso [2010] A. Dal Corso, Projector augmented-wave method: Application to relativistic spin-density functional theory, Phys. Rev. B 82, 075116 (2010).
  * Bihlmayer  _et al._ [2006] G. Bihlmayer, Y. M. Koroteev, P. M. Echenique, E. V. Chulkov, and S. Blügel, The Rashba-effect at metallic surfaces, Surf. Sci. 600, 3888 (2006).
  * Ghosh  _et al._ [2021] K. Ghosh, H. Ma, M. Onizhuk, V. Gavini, and G. Galli, Spin–spin interactions in defects in solids from mixed all-electron and pseudopotential first-principles calculations, Npj Comput. Mater. 7, 1 (2021).
  * Gmitra  _et al._ [2009] M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and J. Fabian, Band-structure topologies of graphene: Spin-orbit coupling effects from first principles, Phys. Rev. B 80, 235431 (2009).
  * Sichau  _et al._ [2019] J. Sichau, M. Prada, T. Anlauf, T. Lyon, B. Bosnjak, L. Tiemann, and R. Blick, Resonance microwave measurements of an intrinsic spin-orbit coupling gap in graphene: A possible indication of a topological state, Phys. Rev. Lett. 122, 046403 (2019).
  * Andersen [1975] O. K. Andersen, Linear methods in band theory, Phys. Rev. B 12, 3060 (1975).
  * MacDonald  _et al._ [1980] A. MacDonald, W. Picket, and D. Koelling, A linearised relativistic augmented-plane-wave method utilising approximate pure spin basis functions, J. Phys. C 13, 2675 (1980).
  * Wimmer  _et al._ [1981] E. Wimmer, H. Krakauer, M. Weinert, and A. Freeman, Full-potential self-consistent linearized-augmented-plane-wave method for calculating the electronic structure of molecules and surfaces: O 2 molecule, Phys. Rev. B 24, 864 (1981).
  * Christensen [1978] N. Christensen, Spin-orbit projected d densities-of-states of Pd, Ag, Pt, and Au, J. Phys. F 8, L51 (1978).
  * Ebert  _et al._ [1988] H. Ebert, P. Strange, and B. L. Gyorffy, Spin-polarized relativistic LMTO method, J. Appl. Phys. 63, 3052 (1988).
  * Rehn  _et al._ [2020] D. A. Rehn, J. M. Wills, T. E. Battelle, and A. E. Mattsson, Dirac’s equation and its implications for density functional theory based calculations of materials containing heavy elements, Phys. Rev. B 101, 085114 (2020).
  * Motamarri  _et al._ [2020] P. Motamarri, S. Das, S. Rudraraju, K. Ghosh, D. Davydov, and V. Gavini, DFT-FE–A massively parallel adaptive finite-element code for large-scale density functional theory calculations, Comput. Phys. Commun. 246, 106853 (2020).
  * Suzuki and Nakao [1999] S. Suzuki and K. Nakao, A fully relativistic full-potential LCAO method for solids, J. Phys. Soc. Jpn. 68, 1982 (1999).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).
  * Zhao  _et al._ [2021] R. Zhao, V. W.-z. Yu, K. Zhang, Y. Xiao, Y. Zhang, and V. Blum, Quasi-four-component method with numeric atom-centered orbitals for relativistic density functional simulations of molecules and solids, Phys. Rev. B 103, 245144 (2021).
  * Philipsen  _et al._ [1997] P. Philipsen, E. Van Lenthe, J. Snijders, and E. Baerends, Relativistic calculations on the adsorption of CO on the (111) surfaces of Ni, Pd, and Pt within the zeroth-order regular approximation, Phys. Rev. B 56, 13556 (1997).
  * Zhao  _et al._ [2016] R. Zhao, Y. Zhang, Y. Xiao, and W. Liu, Exact two-component relativistic energy band theory and application, J. Chem. Phys. 144, 044105 (2016).
  * Boettger and Trickey [2007] J. Boettger and S. Trickey, First-principles calculation of the spin-orbit splitting in graphene, Phys. Rev. B 75, 121402 (2007).
  * Huhn and Blum [2017] W. P. Huhn and V. Blum, One-hundred-three compound band-structure benchmark of post-self-consistent spin-orbit coupling treatments in density functional theory, Phys. Rev. Mater. 1, 033803 (2017).
  * Pisani  _et al._ [2012] C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz, and A. Erba, Cryscor: a program for the post-Hartree–Fock treatment of periodic systems, Phys. Chem. Chem. Phys. 14, 7615 (2012).
  * Booth  _et al._ [2016] G. H. Booth, T. Tsatsoulis, G. K.-L. Chan, and A. Grüneis, From plane waves to local Gaussians for the simulation of correlated periodic systems, J. Chem. Phys. 145, 084111 (2016).
  * McClain  _et al._ [2017] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, Gaussian-based coupled-cluster theory for the ground-state and band structure of solids, J. Chem. Theory Comput. 13, 1209 (2017).
  * Rebolini  _et al._ [2018] E. Rebolini, G. Baardsen, A. S. Hansen, K. R. Leikanger, and T. B. Pedersen, Divide–Expand–Consolidate Second-Order Møller–Plesset Theory with Periodic Boundary Conditions, J. Chem. Theory Comput. 14, 2427 (2018).
  * Crowley  _et al._ [2016] J. M. Crowley, J. Tahir-Kheli, and W. A. Goddard III, Resolution of the band gap prediction problem for materials design, J. Phys. Chem. Lett. 7, 1198 (2016).
  * Kudin and Scuseria [2000] K. N. Kudin and G. E. Scuseria, Linear-scaling density-functional theory with Gaussian orbitals and periodic boundary conditions: Efficient evaluation of energy and forces via the fast multipole method, Phys. Rev. B 61, 16440 (2000).
  * Łazarski  _et al._ [2015] R. Łazarski, A. M. Burow, and M. Sierka, Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole methods, J. Chem. Theory Comput. 11, 3029 (2015).
  * Kadek  _et al._ [2015] M. Kadek, L. Konecny, B. Gao, M. Repisky, and K. Ruud, X-ray absorption resonances near L2, 3-edges from real-time propagation of the Dirac–Kohn–Sham density matrix, Phys. Chem. Chem. Phys. 17, 22566 (2015).
  * Komorovskỳ  _et al._ [2006] S. Komorovskỳ, M. Repiskỳ, O. L. Malkina, V. G. Malkin, I. Malkin, and M. Kaupp, Resolution of identity Dirac-Kohn-Sham method using the large component only: Calculations of g-tensor and hyperfine tensor, J. Chem. Phys. 124, 084108 (2006).
  * Repiskỳ  _et al._ [2009] M. Repiskỳ, S. Komorovskỳ, O. L. Malkina, and V. G. Malkin, Restricted magnetically balanced basis applied for relativistic calculations of indirect nuclear spin–spin coupling tensors in the matrix Dirac–Kohn–Sham framework, Chem. Phys. 356, 236 (2009).
  * Pisani and Dovesi [1980] C. Pisani and R. Dovesi, Exact-exchange Hartree–Fock calculations for periodic systems. I. Illustration of the method, Int. J. Quantum Chem. 17, 501 (1980).
  * Dovesi  _et al._ [1983] R. Dovesi, C. Pisani, C. Roetti, and V. Saunders, Treatment of Coulomb interactions in Hartree-Fock calculations of periodic systems, Phys. Rev. B 28, 5781 (1983).
  * Dovesi  _et al._ [2020] R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D’arco, Y. Noël, M. Rérat, P. Carbonniere, _et al._ , The CRYSTAL code, 1976–2020 and beyond, a long story, J. Chem. Phys. 152, 204111 (2020).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , Recent developments in the PySCF program package, J. Chem. Phys. 153, 024109 (2020).
  * Suhai  _et al._ [1982] S. Suhai, P. Bagus, and J. Ladik, An error analysis for Hartree-Fock crystal orbital calculations, Chem. Phys. 68, 467 (1982).
  * Peintinger  _et al._ [2013] M. F. Peintinger, D. V. Oliveira, and T. Bredow, Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations, J. Comput. Chem. 34, 451 (2013).
  * Zhu and Chan [2021] T. Zhu and G. K.-L. Chan, All-electron Gaussian-based G 0 W 0 for valence and core excitation energies of periodic systems, J. Chem. Theory Comput. 17, 727 (2021).
  * Peralta  _et al._ [2005] J. E. Peralta, J. Uddin, and G. E. Scuseria, Scalar relativistic all-electron density functional calculations on periodic systems, J. Chem. Phys. 122, 084108 (2005).
  * Usvyat  _et al._ [2011] D. Usvyat, B. Civalleri, L. Maschio, R. Dovesi, C. Pisani, and M. Schütz, Approaching the theoretical limit in periodic local MP2 calculations with atomic-orbital basis sets: the case of LiH, J. Chem. Phys. 134, 214105 (2011).
  * Perry  _et al._ [2001] J. K. Perry, J. Tahir-Kheli, and W. A. Goddard III, Antiferromagnetic band structure of La 2 CuO 4: Becke-3–Lee-Yang-Parr calculations, Phys. Rev. B 63, 144510 (2001).
  * Daga  _et al._ [2020] L. E. Daga, B. Civalleri, and L. Maschio, Gaussian basis sets for crystalline solids: All-purpose basis set libraries vs system-specific optimizations, J. Chem. Theory Comput. 16, 2192 (2020).
  * Jensen [2013] F. Jensen, Analysis of energy-optimized gaussian basis sets for condensed phase density functional calculations, Theor. Chem. Acc. 132, 1 (2013).
  * Stanton and Havriliak [1984] R. E. Stanton and S. Havriliak, Kinetic balance: A partial solution to the problem of variational safety in Dirac calculations, J. Chem. Phys. 81, 1910 (1984).
  * Baerends  _et al._ [1973] E. J. Baerends, D. Ellis, and P. Ros, Self-consistent molecular Hartree—Fock—Slater calculations I. The computational procedure, Chem. Phys. 2, 41 (1973).
  * Jung  _et al._ [2005] Y. Jung, A. Sodt, P. M. Gill, and M. Head-Gordon, Auxiliary basis expansions for large-scale electronic structure calculations, Proc. Natl. Acad. Sci. U.S.A. 102, 6692 (2005).
  * Varga [2008] Š. Varga, Long-range analysis of density fitting in extended systems, Int. J. Quantum Chem. 108, 1518 (2008).
  * Kadek  _et al._ [2022] M. Kadek, B. Wang, M. Joosten, W.-C. Chiu, F. Mairesse, M. Repisky, K. Ruud, and A. Bansil, Dataset for ”Band structures and Z2 invariant of 2D transition metal dichalcogenides from fully relativistic Dirac–Kohn–Sham theory using Gaussian-type orbitals”, ZENODO (2022).
  * [106] The ReSpect code is available at www.respectprogram.org.
  * Schwarz and Wallmeier [1982] W. Schwarz and H. Wallmeier, Basis set expansions of relativistic molecular wave equations, Mol. Phys. 46, 1045 (1982).
  * Dyall [1994] K. G. Dyall, An exact separation of the spin-free and spin-dependent terms of the Dirac–Coulomb–Breit Hamiltonian, J. Chem. Phys. 100, 2118 (1994).
  * Repisky  _et al._ [2015] M. Repisky, L. Konecny, M. Kadek, S. Komorovsky, O. L. Malkin, V. G. Malkin, and K. Ruud, Excitation energies from real-time propagation of the four-component Dirac–Kohn–Sham equation, J. Chem. Theory Comput. 11, 980 (2015).
  * Balasubramani  _et al._ [2020] S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, _et al._ , TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations, J. Chem. Phys. 152, 184107 (2020).
  * Apra  _et al._ [2020] E. Apra, E. J. Bylaska, W. A. De Jong, N. Govind, K. Kowalski, T. P. Straatsma, M. Valiev, H. J. van Dam, Y. Alexeev, J. Anchell, _et al._ , NWChem: Past, present, and future, J. Chem. Phys. 152, 184102 (2020).
  * Saue  _et al._ [2020] T. Saue, R. Bast, A. S. P. Gomes, H. J. A. Jensen, L. Visscher, I. A. Aucar, R. Di Remigio, K. G. Dyall, E. Eliav, E. Fasshauer, _et al._ , The DIRAC code for relativistic molecular calculations, J. Chem. Phys. 152, 204104 (2020).
  * Olsen  _et al._ [2020] J. M. H. Olsen, S. Reine, O. Vahtras, E. Kjellgren, P. Reinholdt, K. O. Hjorth Dundas, X. Li, J. Cukras, M. Ringholm, E. D. Hedegård, _et al._ , Dalton Project: A Python platform for molecular-and electronic-structure simulations of complex systems, J. Chem. Phys. 152, 214115 (2020).
  * Belpassi  _et al._ [2020] L. Belpassi, M. De Santis, H. M. Quiney, F. Tarantelli, and L. Storchi, BERTHA: Implementation of a four-component Dirac–Kohn–Sham relativistic framework, J. Chem. Phys. 152, 164118 (2020).
  * Obara and Saika [1986] S. Obara and A. Saika, Efficient recursive computation of molecular integrals over Cartesian Gaussian functions, J. Chem. Phys. 84, 3963 (1986).
  * Boys [1950] S. F. Boys, Electronic wave functions-I. A general method of calculation for the stationary states of any molecular system, Proc. R. Soc. A 200, 542 (1950).
  * Sierka  _et al._ [2003] M. Sierka, A. Hogekamp, and R. Ahlrichs, Fast evaluation of the Coulomb potential for electron densities using multipole accelerated resolution of identity approximation, J. Chem. Phys. 118, 9136 (2003).
  * Sun  _et al._ [2017] Q. Sun, T. C. Berkelbach, J. D. McClain, and G. K.-L. Chan, Gaussian and plane-wave mixed density fitting for periodic systems, J. Chem. Phys. 147, 164119 (2017).
  * Kane and Mele [2005b] C. L. Kane and E. J. Mele, Z2 topological order and the quantum spin hall effect, Phys. Rev. Lett. 95, 146802 (2005b).
  * Gresch  _et al._ [2017] D. Gresch, G. Autes, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, Z2Pack: Numerical implementation of hybrid Wannier centers for identifying topological materials, Phys. Rev. B 95, 075146 (2017).
  * Fu and Kane [2007] L. Fu and C. L. Kane, Topological insulators with inversion symmetry, Phys. Rev. B 76, 045302 (2007).
  * Te Velde and Baerends [1991] G. Te Velde and E. Baerends, Precise density-functional method for periodic structures, Phys. Rev. B 44, 7888 (1991).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Visscher and Dyall [1997] L. Visscher and K. G. Dyall, Dirac–Fock atomic electronic structure calculations using different nuclear charge distributions, At. Data Nucl. Data Tables 67, 207 (1997).
  * Dyall [2007] K. G. Dyall, Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the 4d elements Y–Cd, Theor. Chem. Acc. 117, 483 (2007).
  * Dyall [2004] K. G. Dyall, Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the 5d elements Hf–Hg, Theor. Chem. Acc. 112, 403 (2004).
  * Dyall and Gomes [2010] K. G. Dyall and A. S. Gomes, Revised relativistic basis sets for the 5d elements Hf–Hg, Theor. Chem. Acc. 125, 97 (2010).
  * Dyall [2002] K. G. Dyall, Relativistic and nonrelativistic finite nucleus optimized triple-zeta basis sets for the 4p, 5p and 6p elements, Theor. Chem. Acc. 108, 335 (2002).
  * Dyall [2006] K. G. Dyall, Relativistic quadruple-zeta and revised triple-zeta and double-zeta basis sets for the 4p, 5p, and 6p elements, Theor. Chem. Acc. 115, 441 (2006).
  * Kośmider  _et al._ [2013] K. Kośmider, J. W. González, and J. Fernández-Rossier, Large spin splitting in the conduction band of transition metal dichalcogenide monolayers, Phys. Rev. B 88, 245436 (2013).
  * [132] See Supplemental Material at [URL will be inserted by publisher] for the direct SOC-induced band gaps in the 1T phase of MX2, and band structures in the 2H, 1T, and 1T’ phases plotted with and without the SOC.
  * Van Lenthe  _et al._ [1996] E. Van Lenthe, R. Van Leeuwen, E. Baerends, and J. Snijders, Relativistic regular two-component Hamiltonians, Int. J. Quantum Chem. 57, 281 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010).
  * Liu  _et al._ [2023] J. Liu, T. Shen, J.-C. Ren, S. Li, and W. Liu, Role of van der Waals interactions on the binding energies of 2D transition-metal dichalcogenides, Appl. Surf. Sci. 608, 155163 (2023).
  * Setyawan and Curtarolo [2010] W. Setyawan and S. Curtarolo, High-throughput electronic band structure calculations: Challenges and tools, Comput. Mater. Sci. 49, 299 (2010).
  * Eda  _et al._ [2012] G. Eda, T. Fujita, H. Yamaguchi, D. Voiry, M. Chen, and M. Chhowalla, Coherent atomic and electronic heterostructures of single-layer MoS2, ACS nano 6, 7311 (2012).
  * Visscher  _et al._ [1999] L. Visscher, T. Enevoldsen, T. Saue, H. J. A. Jensen, and J. Oddershede, Full four-component relativistic calculations of NMR shielding and indirect spin–spin coupling tensors in hydrogen halides, J. Comput. Chem. 20, 1262 (1999).