Current density functional framework for spin–orbit coupling: Extension to periodic systems

Spin–orbit coupling induces a current density in the ground state, which consequently requires a generalization for meta -generalized gradient approximations. That is, the exchange-correlation energy has to be constructed as an explicit functional of the current density and a generalized kinetic energy density has to be formed to satisfy theoretical constraints. Herein, we generalize our previously presented formalism of spin–orbit current density functional theory [Holzer et al. , J. Chem. Phys. 157 , 204102 (2022)] to non-magnetic and magnetic periodic systems of arbitrary di-mension. Besides the ground-state exchange-correlation potential, analytical derivatives such as geometry gradients and stress tensors are implemented. The importance of the current density is assessed for band gaps, lattice constants, magnetic transitions, and Rashba splittings. For the latter, the impact of the current density may be larger than the deviation between different density functional approximations.

To account for the change in Fermi hole curvature, the kinetic energy density τ is generalized with the current density ⃗ j.In a two-component (2c) non-collinear formalism, [41][42][43][44][45][46][47][48][49][50][51] this results in the generalized current density based on the spin-up and down quantities. 36These are formed with the particle and spin-magnetization contributions, e.g. the spin-up and down electron density n ↑,↓ follows as with the particle density n, the spin-magnetization vector ⃗ m, and the spin density s.Therefore, the exchange-correlation (XC) energy of a "pure" functional 2 explicitly depends on the current density where f XC describes the density functional approximation and γ ↑↓ (⃗ r) = 1 4 ⃗ ∇n ↑ (⃗ r) • ⃗ ∇n ↓ (⃗ r) , hence γ ↑↓ = γ ↓↑ .For a Kramers-restricted system, the spin-magnetization vector and the particle current density vanish.However, the spin-current density is generally non-zero and thus it is still necessary to form the generalized kinetic energy density.
In this work, we will extend our previous CDFT formulation 36 to non-magnetic and magnetic periodic systems.We assess the importance of the current density for band gaps, cell structures, magnetic transitions, and Rashba splittings with common meta-GGAs.Together with previous studies on the impact of spin-orbit-induced current densities for meta-GGAs 36,38,[52][53][54] this will further help to set guidelines and recommendations for the application of CDFT to the different properties and functionals for discrete and periodic systems.

II. THEORY
The one-particle density matrix associated with the twocomponent spinor functions ⃗ where V FBZ is the volume of the first Brillouin zone (FBZ), ε i the energy eigenvalue, and ε F the Fermi level.The spinor functions are expanded with Bloch functions, φ µ , based on atomic orbital (AO) basis functions, χ µ , as N UC denotes the number of electrons in the unit cell (UC) and ⃗ L the lattice vector.Thus, all density variables are available from the AO density matrix in position space given by with the expansion coefficients c µi and σ , σ ′ ∈ {α, β }.The complete two-component AO density matrix reads In the spirit of Bulik et al., 47 the real symmetric (RS), real antisymmetric (RA), imaginary symmetric (IS), and imaginary antisymmetric (OA) linear combinations are formed.Of course, the same-spin antisymmetric contributions are zero.The electron density and its derivatives are available from the symmetric linear combinations The particle current density ⃗ j p and the spin-current densities ⃗ j u , with u ∈ {x, y, z}, are obtained from the antisymmetric linear combinations Following our molecular ansatz, 36 the scalar part of the XC potential is obtained as and the spin-magnetization part with u ∈ {x, y, z} reads The spin blocks of the Kohn-Sham equations are formed by combining the scalar and spin-magnetization contributions with the respective Pauli matrices, c.f. Refs.36 and 56.After transformation to the k space, the Kohn-Sham equations can be solved as usual.
For non-magnetic or closed-shell system, ⃗ m and ⃗ j p vanish so that a Kramers-restricted framework can be applied and time-reversal symmetry 57 may be exploited.However, the spin current densities are still non-zero and hence the spin-up and spin-down quantities ⃗ j ↑,↓ contribute to the XC potential through the diamagnetic or quadratic terms. 366][67] The numerical integration of the exchange-correlation potential is carried with the algorithm of Ref. 59.Note that the current density also leads to antisymmetric contributions.Integration weights are constructed according to Stratmann et al. 68 Geometry gradients and stress tensors are implemented based on previous work, as it essentially involves further derivatives of the basis functions. 56,61,63We note that all integrals and the exchange potential for the Kohn-Sham equations are evaluated in the position space, which allows to exploit sparsity.The increase in the computational cost by calculating the current density contributions on a grid is modest-especially compared to the inclusion of the current density through Hartree-Fock exchange with hybrid functions.

III. COMPUTATIONAL METHODS
0][71][72][73] Two Pt atoms are placed in the unit cell with the cell parameter d ranging from 4.0 Å to 6.0 Å. Calculations are performed at the TPSS/dhf-SVP-2c 74,75 level employing Dirac-Fock effective core potentials (DF-ECPs), 76 replacing 60 electrons (ECP-60).A Gaussian smearing of 0.01 Hartree 77 is used and a k-mesh with 32 points is applied.Integration grids, convergence settings, etc. are given in the Supplementary Material.We note in passing that the Karlsruhe dhf-type basis sets were optimized for discrete systems, 75 however, the corresponding pob-type basis sets for periodic calculations [78][79][80][81][82][83] are not yet available with tailored extensions for spin-orbit two-component calculations. 84Therefore and for consistency with previous studies, we apply the dhf-type basis sets.

A. Linear Pt Chain
The one-dimensional linear Pt chain is a common reference system for a transition to a magnetic system.][71][72][73] This is also confirmed at the scalar and spin-orbit TPSS level in Fig. 1.For small cell parameters from d = 4.0 to d = 4.9 Å, the open-shell initial guess converges to a closedshell non-magnetic solution in the self-consistent field (SCF) procedure.The most favorable total energy is found for d = 4.7 Å in agreement with previous studies based on GGA functionals. 56,72Further, CDFT and DFT lead to a very similar potential energy surface or a similar behavior of the relative energy with respect to the cell parameter d.
At the scalar level, the transition to a magnetic material occurs between d ≈ 5.2 Å and d ≈ 5.3 Å. Inclusion of spin-orbit coupling shifts this transition to a smaller d parameter of about 5.0 Å.Here, the current-dependent variant of TPSS (cTPSS) leads to a lower total energy for both closed-shell and openshell solutions, i.e. a more negative energy, and a larger magnetic moment.The impact of the current density is generally larger for the magnetic solution than for the closed-shell state, see also the Supplementary Material for detailed results.For the closed-shell solution, the difference in the total energy by inclusion of the current density is too small to be visible in panel (b) of Fig. 1.This finding is qualitatively in line with our previous studies on molecular systems. 36nclusion of the current density also consistently leads to a larger magnetic moment.For instance, TPSS predicts a magnetic moment of 1.04 µ B /atom at d = 6.0 Å, whereas a value of 1.13 µ B /atom is found with cTPSS.Generally, an increase in the range of 5-8 % is observed after the transition to a magnetic solution.

B. Transition-Metal Dichalcogenide Monolayers
Transition-metal dichalcogenide monolayers have many interesting physical properties such as the quantum spin Hall, 98 non-linear anomalous Hall, 99 and Rashba effects. 100In the H2 phase, time-reversal symmetry holds for the non-magnetic systems but space-inversion symmetry is lost.Thus, spinorbit coupling lifts the degeneracy of the valence band at the K point in the Brillouin zone.For the Mo and W systems, this Rashba splitting is very pronounced and values between 0.1 and 0.5 eV are obtained with relativistic all-electron methods. 93,101and gaps and Rashba splittings at the K point obtained with DFT and CDFT approaches are listed in Tab.I.One the one hand, the impact on the band gaps is rather small, with TASK and r 2 SCAN showing the largest changes but not exceeding 0.1 eV.Here, the deviation between the different DFAs is far larger.One the other hand, the Rashba splitting is very sensitive to the inclusion of the current density.For instance, the results change from 0.436 eV to 0.483 eV and 0.465 eV to 0.566 eV for WTe 2 with r 2 SCAN and TASK, respectively.
The most pronounced current-density effects are found for TASK throughout all systems, which is in line with molecular studies. 36,53,102Here, the current density increase the Rashba splitting by 25 % on average.r 2 SCAN ranks second in this regard with 13 % followed by M06-L with 6 %.For PKZB, Tao-Mo and TPSS, changes of only 1-3 % are observed.Among the different monolayers, MoTe 2 leads to the largest impact of the current density for all DFAs, with a relative deviation of 36 % between DFT and CDFT for TASK.The impact of the current density for the DFAs can be rationalized by the enhancement factor. 102,103 Silver Halide Crystals The band gaps and optimized lattice constant of AgI with various meta-GGAs are listed in Tab.II.Results for AgCl and AgBr are presented in the Supplementary Material.Overall, the current density is of minor importance for the band gaps.Changes are in the order of meV.These results are in qualitative agreement with the two-dimensional MoTe 2 monolayer, which consists of atoms from the same row of the periodic table of elements.
Likewise the current density does not lead to major changes for the cell structure and the lattice constant.The small impact on the lattice constant can be rationalized by the comparably small impact of spin-orbit coupling on the cell structure. 56D3 dispersion correction with Becke-Johnson damping 8,[104][105][106][107] (D3-BJ) leads to much larger changes than the application of CDFT.Therefore, other computational parameters than the inclusion of the current density for meta-GGAs are more important for the cell structures of the silver halide crystals.

V. CONCLUSION
We have extended our previous formulation of spin-orbit current density functional theory to periodic systems of arbitrary dimension.The impact of the current density was assessed for various properties of non-magnetic and magnetic systems.Here, the band gaps and lattice constants are not notably affected.In contrast, the Rashba splitting which is only due to spin-orbit coupling is substantially affected.The inclusion of the current density for a given functional may lead to larger changes than the deviation of the results among different functionals.
With the present work, CDFT is now applicable to a wide range of chemical and physical properties of discrete and periodic systems, including analytic first-order property calculations.We generally recommend to include the current density for r2 SCAN and the Minnesota functionals, if available in the used electronic structure code.For TASK, inclusion of the current density is clearly mandatory as it leads to substantial changes of the results.

SUPPLEMENTARY MATERIAL
Supporting Information is available with all computational details and data.

II. Two-Dimensional Transition-Metal Dichalcogenide Monolayers
A. Computational Details B. Results

III. Three-Dimensional Silver Halide Crystals
A. Computational Details

References
TABLE I. Total SCF energies in Hartree obtained with the 2c KR DFT, 2c KU DFT, 2c KR CDFT, and 2c KU CDFT approach, and energies for the limit of a vanishing Gaussian smearing (limit).The expectation values of spin S x are listed for the open-shell calculations.SK, 21 TPSS, 6 Tao-Mo, 22 and PKZB 23 density functional approximations (DFAs).Note that we use Libxc 24-26 for M06-L, r 2 SCAN, TASK, Tao-Mo, and PKZB.The dhf-TZVP-2c orbital and auxiliary basis sets are applied. 2Originally, the dhf-type basis sets only employ ECPs for elements beyond Kr.That is, Dirac-Fock ECPs are employed for Mo (ECP-28), W (ECP-60), and Te (ECP-28). 1,27,28However, this leads to rather large deviations from the relativistic all-electron calculations at the PBE level for MoSe 2 and WSe 2 . 29,30See Tabs.II and III.Therefore, we applied a small-core Dirac-Fock ECP for Se (ECP-10) as well. 27,31The corresponding results are shown in Tabs.IV and V.
Large integration grids 3,4 (grid size 3) are employed for the numerical integration of the exchange-correlation potential. 5Again, default settings are applied for RI-J, CFMM, and the transformation to the orthogonal basis set for the diagonalization of the Kohn-Sham or Fock matrix. 8,9The SCF procedure is converged with a threshold of 10 −7 Hartree and a k-mesh of 33 × 33 points is used.Structures are taken from Ref. 29.

FIG. 1 .
FIG. 1.(a) Dependence of the energy on the cell parameter in units of Hartree per atom for one-component (1c) restricted Kohn-Sham (RKS) and 1c unrestriced Kohn-Sham (UKS) at the TPSS/dhf-SVP-2c level.(b) Dependence of the energy on the cell parameter in units of Hartree per atom for 2c KR and 2c KU with the spin aligned along x (S x ).(c) Magnetic moment in units of Bohr's magneton µ B per atom for the spin contribution of 1c UKS, 2c KU S x DFT and CDFT approach.Periodicity is along the x direction for one-dimensional systems. 55The open-shell solutions are energetically favored compared to the respective closed-shell solutions and the S x alignment is preferred over S y,z .Results without the current-independent TPSS functional are taken from Ref. 56.Detailed results are listed in the Supplementary Material.

TABLE I .
Band gaps and Rashba splittings of the valence band at the K point for transition-metal dichalcogenide monolayers at the 1c DFT, 2c DFT, and 2c CDFT level with the dhf-TZVP-2c basis set and DF-ECPs for all atoms.All values in eV.Results for MoS 2 and WS 2 are listed in the Supplementary Material.

TABLE II .
Optimized lattice constant a (in Å, rocksalt structure) of three-dimensional AgI and band gaps (in eV) at high symmetry points of the first Brillouin zone with various density functional approximations and the dhf-SVP basis sets.A "c" indicates the currentdependent variant of the DFA.Results for DFAs without inclusion of the current density are taken from Ref. 56, except for M06-L.Calculations are performed without dispersion correction (no D3) and with the D3 correction using Becke-Johnson damping (D3-BJ).

TABLE II .
Band gaps and Rashba splittings of the valence band at the K point for transition-metal dichalcogenide monolayers at the 1c DFT, 2c DFT, and 2c CDFT level with the dhf-TZVP-2c basis set and DF-ECPs for Mo and Te.All values in eV.

TABLE III .
Band gaps and Rashba splittings of the valence band at the K point for transition-metal dichalcogenide monolayers at the 1c DFT, 2c DFT, and 2c CDFT level with the dhf-TZVP-2c basis set and DF-ECPs for W and Te.All values in eV.

TABLE IV .
Band gaps and Rashba splittings of the valence band at the K point for transition-metal dichalcogenide monolayers at the 1c DFT, 2c DFT, and 2c CDFT level with the dhf-TZVP-2c basis set and DF-ECPs for Mo, Se, and Te.All values in eV.

TABLE V .
Band gaps and Rashba splittings of the valence band at the K point for transition-metal dichalcogenide monolayers at the 1c DFT, 2c DFT, and 2c CDFT level with the dhf-TZVP-2c basis set and DF-ECPs for W, Se, and Te.All values in eV.

TABLE VI .
2ptimized lattice constant a (in Å, rocksalt structure) of three-dimensional AgCl and band gaps (in eV) at high symmetry points of the first Brillouin zone with various density functional approximations and the dhf-SVP basis sets.2ADirac-FockECP is applied for Ag (ECP-28).35A"c" indicates the currentdependent variant of the DFA.Results for DFAs without inclusion of the current density are taken from

TABLE VII .
2ptimized lattice constant a (in Å, rocksalt structure) of three-dimensional AgBr and band gaps (in eV) at high symmetry points of the first Brillouin zone with various density functional approximations and the dhf-SVP basis sets.2ADirac-FockECP is applied for Ag (ECP-28).35A"c"indicates the currentdependent variant of the DFA.Results for DFAs without inclusion of the current density are taken from Ref. 16, except for M06-L.Calculations are performed without dispersion correction (no D3) and with the D3 correction using Becke-Johnson damping (D3-BJ).39,40

TABLE VIII .
2ptimized lattice constant a (in Å, rocksalt structure) of three-dimensional AgI and band gaps (in eV) at high symmetry points of the first Brillouin zone with various density functional approximations and the dhf-SVP basis sets.2Dirac-FockECPs are applied for Ag (ECP-28) and I (ECP-28).35,36A"c"indicates the current-dependent variant of the DFA.Results for DFAs without inclusion of the current density are taken from Ref. 16, except for M06-L.Calculations are performed without dispersion correction (no D3) and with the D3 correction using Becke-Johnson damping (D3-BJ).39,40