A Visual Representation for Accurate Local Basis Set Construction and Optimization: A Case Study of SrTiO₃ with Hybrid DFT Functionals

Abstract

The linear combination of atomic orbitals (LCAO) method is advantageous for calculating important bulk and surface properties of crystals and defects in/on them. Compared to plane wave calculations and contrary to common assumptions, hybrid density functional theory (DFT) functionals are actually less costly and easier to implement in LCAO codes. However, choosing the proper basis set (BS) for the LCAO calculations representing Guassian-type functions is crucial, as the results depend heavily on its quality. In this study, we introduce a new basis set (BS) visual representation, which helps us (1) analyze the collective behavior of individual atoms’ shell exponents (s, p, and d), (2) better compare different BSs, (3) identify atom-type invariant relationships, and (4) suggest a robust method for building a local all-electron BS (denoted as BS1) from scratch for each atom type. To compare our BS1 with the others existing in the literature, we calculate the basic bulk properties of ${\mathrm{SrTiO}}_3$ (STO) in cubic and tetragonal phases using several hybrid DFT functionals (B3LYP, PBE0, and HSE06). After adjusting the exact Hartree–Fock (HF) exchange of PBEx, HSEx, and the state-of-the-art meta-GGA hybrid ${\mathrm{r}}^2$SCANx functionals, we find the ${\mathrm{r}}^2$SCAN15 and HSE27 for BS1, with the amount of exact HF exchange of 0.15 and 0.27, respectively, perform equally well for reproducing several most relevant STO properties. The proposed robust BS construction scheme has the advantage that all parameters of the obtained BS can be reoptimized for each new material, thus increasing the quality of DFT calculation predictions.

1. Introduction

Strontium titanate (STO) is a prototypical material in the perovskite family and is a principal focus for technology and research. Its unique electronic and lattice dynamics properties [1] have established it as a model material in various applications for many years (e.g., [2,3,4,5]). Understanding these properties on an atomistic level requires the best available methods. Virtually every up-to-date density functional theory (DFT) method is used to study STO’s challenging and fascinating properties, a trend likely to continue into the foreseeable future.

We choose STO as a model material for basis set (BS) testing for two main reasons. Firstly, it is a $d^0$ transition metal oxide, making the bandgap problem relevant here. Secondly, it exhibits a specific phonon frequency behavior at the $\Gamma$-point and a high dielectric constant. Below about 105 K, STO undergoes a structural phase transition from cubic $Pm\overline{3}m$ (space group (SG) 221) to tetragonal $I4/mcm$ (SG 140) phase [6]. The phase transformation is associated with antiferrodistortive (AFD) displacements when oxygen octahedra rotate in antiphase around the [001] direction. According to a group-theoretical analysis, the transition from cubic to tetragonal phase occurs due to the softening of the $R_4^{+}$ phonon mode [7]. Nevertheless, STO is also known as incipient ferroelectric, in which quantum fluctuations suppress the ferroelectric (FE) mode [8] related to the phonon frequency at the $\Gamma$-point.

Spectroscopic ellipsometry experiments suggest that the direct gap of the STO increases from 3.7 to 3.8 eV with temperature without any noticeable effect at the structural phase transition temperature [9,10,11]. The indirect gap in the cubic phase, determined from the comparatively weaker signal, decreases from 3.23 eV with temperature [9,11].

Thus, the theoretical STO model has to account for many geometric and electronic properties in the unified approach. There exist numerous DFT studies of perfect bulk STO [7,12,13,14,15,16,17], employing different DFT methods and, thus, different degrees of accuracy. We distinguish three different paradigms in applications of DFT to the studies on STO: Hubbard-model-based approaches, $GW$ calculations, and hybrid density functionals.

The ONETEP code with a linear scaling approach to the DFT calculations, atom-centered BS (Wannier functions), and DFT + $U^{d,p}$ + $J^{d,p}$ method (which means the application of U- and J-parameters to both Ti-$3d$ and O-$2p$ electrons) as used by Winszewski et al. [15] suggest significant improvements in the description of cubic STO electronic and structural properties and of defects behavior therein. Even though the method implies a ’reasonable’ band gap value for the oxygen vacancy and hydrogen interstitial behavior, the direct/indirect band gap is still smaller by 0.53/0.33 eV compared with the experimental value. A better agreement with the experimental indirect band gap was achieved by applying three different U-values to Sr-$4d$, O-$2p$, and Ti-$3d$ electrons in a very recent study of Derkaoui et al. [17]. Alternatively, in the DFT + U + V calculations of Ricca et al. [16] extended for the intersite interactions, the Hubbard U- and V-parameters are obtained self-consistently within the density functional perturbation theory, resulting in a band gap close, albeit with slightly smaller lattice parameters, to the experimental value (see our analysis in Section 5).

Contrary to standard GGA and DFT + U (with the application of only one U-parameter for Ti 3d electrons) calculations [15,16], earlier $GW$ calculations tend to significantly overestimate the band gap [18,19,20] with some exceptions [21,22,23]. For example, Sponza et al. [22] showed the direct band gap of 3.76 eV calculated with the $G_0{W}_0$ approximation, whereas Hamann and Vanderbilt [23] obtained the indirect band gap of 3.32 eV from the QS$GW$ calculations with the maximally localized Wannier functions for quasiparticles of the cubic STO. However, we would also like to mention a more recent study by Bhandari et al. [24] based on the extended quasiparticle self-consistent $GW$ (i.e., QS$GW$ + LPC + 0.8$\Sigma$), which showed the calculated band gap values for cubic and tetragonal STO consistent with the experiments. They also emphasized including Sr-$4p$ and Ti-$3p$ states as the valence electrons.

Even though all the calculation methods discussed above suggest essential contributions to the studies of STO, the hybrid DFT functionals within the framework of LCAO still represent the most effective approach in terms of calculation time and accuracy. First, such calculations allow for the simultaneous treatment of several material properties. Second, the point defects require large unit cells and low symmetries, which makes their calculations demanding. However, the hybrid functionals can be efficiently implemented using LCAO, which significantly reduces the calculation time compared to the plane wave calculations without imposing additional parameters. We usually do not experience problems with the LCAO calculations with the hybrid DFT functionals and large supercells. Therefore, the LCAO approach is worthy of study and development as an important alternative and/or complement to more popular plane-wave calculations. Our team at the Institute of Solid State Physics of the University of Latvia [25] has significantly contributed to applying the LCAO approach to solids since the 1980s. Nevertheless, it needs to be recognized that finding an optimal BS for the LCAO calculations for crystals is a challenge.

In the present study, we discuss the basic BS properties overlooked so far in the literature and, based on that, suggest a practical method for all-electron BS construction from scratch. We will show that our approach to the local BS construction and optimization is robust and easy to follow. In addition, we will make proper comparisons with the other BSs for STO existing in the literature. We will apply the state-of-the-art meta-GGA hybrid ${\mathrm{r}}^2$SCANx [26,27] functional available in the CRYSTAT23 code for the first time to the calculation of STO properties in the cubic and tetragonal phases. We will compare four different hybrid density functionals with and without non-local terms and adjusted amounts of Hartree–Fock exact exchange, giving the best agreement with the experimental values of several STO properties.

2. Local Basis Set Selection

In the LCAO approach, a set of Gaussian-type functions is selected from a database (e.g., CRYSTAL website [28] or BSE library [29]). Notably, all-electron BS of triple-zeta valence with polarization quality, known as pob-TZVP [30,31] and pob-TZVP-rev2 (an improved version with reduced BSSE error) [32,33], have recently been developed. In our experience, the pob-TZVP BSs have been effective for NASICON-structured compounds ${\mathrm{Na}}_{1+2\mathrm{x}}{\mathrm{Mn}}_{\mathrm{x}}{\mathrm{Ti}}_{2-\mathrm{x}}{\left({\mathrm{PO}}_4\right)}_3$ [34] and for H impurities in (La,Sr)Fe${\mathrm{O}}_{3-\delta }$ [35] without needing a preliminary optimization of their Gaussian-type functions. In this study, the pob-TZVP-rev2 BS is also used for comparison.

There is, however, a clear distinction between the local BSs for atoms in crystals and molecules. Atomic or molecular BSs usually contain diffuse exponents with values smaller than $0.1$ ${\mathrm{Bohr}}^{-2}$. As a rule of thumb, in the case of crystals, they are often removed [30] without a complete BS re-optimization. In particular, the contracted exponents, which describe the inner-core electrons, are left unchanged or are replaced by pseudopotentials. However, many DFT studies show that even optimization of the valence exponents considerably improves agreement with the experimental results, especially for bulk materials. The optimization of the BS is carried out in the following papers [32,36,37,38,39,40,41]. Furthermore, the discussion of the role of local BS optimization and comparison between the BSs is discussed in [7,42,43]. The main criterion for BS optimization in these studies is the system’s total energy. BS optimizer recently introduced in CRYSTAL23 additionally uses the penalty function to avoid linear dependency problem [44].

Using identic BS for the same atom in different compounds is convenient. Nevertheless, such BSs, particularly those made for one system (such as an atom, molecule, or crystal), must be used cautiously for a different system. It needs to be said that we stand for the material-specific approach for the local BS optimization and not for the local BS transferability between the materials. These considerations have been discussed previously [39,45,46], while these works significantly contributed to the field, we aim to complement them with a more detailed analysis of such local BS properties as the ratios of Gaussian-type functions exponents and contraction coefficients. Based on this, we suggest a robust scheme for obtaining new BS or reoptimizing existing ones for each material-specific case.

3. Methodology

3.1. DFT Functionals

The first-principles calculations are performed using the CRYSTAL23 [44,47] computer code within the DFT formalism. The single-particle wave functions are expanded as a linear combination of Bloch functions, which, in turn, are a linear combination of atomic orbitals (Gaussian-type functions). A comparison between the exchange-correlation functionals and BSs lies at the heart of the present study. As is well known, the hybrid density functionals remain an essential tool to reproduce the bandgap values in perovskite oxides [48,49,50]. In its turn, the hybrid density functionals can be expressed in the general form suggested by Becke in [51,52], and, thus, re-written as [44]

$$\begin{array}{c}\hfill E_{xc}=E_x^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}}+A(E_x^{\mathrm{HF}}-E_x^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}})+(1-A)B(E_x^{\mathrm{DFA}}-E_x^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}})+E_c^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}}+C(E_c^{\mathrm{DFA}}-E_c^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}}),\end{array}$$

(1)

where $E_{x,c}^{\mathrm{L}\left(\mathrm{S}\right)\mathrm{DA}}$ the local density functional exchange and correlation contributions and $E_{x,c}^{\mathrm{DFA}}$ the semi-local density functional exchange and correlation contributions (such as DFA = GGA). The non-local density functional exchange and correlation contributions are given by B and C parameters, respectively, whereas A finds the amount of exact Hartree–Fock exchange $E_x^{\mathrm{HF}}$. So, the one parameter (A) global density functional PBEx is defined from Equation (1) using $B=1$ and $C=1$ and, thus, reads for the exchange-correlation part

$$E_{xc}^{\mathrm{PBEx},{\mathrm{r}}^2\mathrm{SCANx}}=A{E}_x^{\mathrm{HF}}+(1-A)E_x^{\mathrm{DFA}}+E_c^{\mathrm{DFA}},$$

(2)

where DFA = PBE or ${\mathrm{r}}^2$SCAN [26] in the present study. A successor of PBEx, range-separated hybrid DFT functional, has demonstrated very good results in many cases even in its reduced form HSEx [53]:

$$E_{xc}^{\mathrm{HSEx}}=A{E}_x^{\mathrm{HF},\mathrm{SR}}\left(\omega \right)+(1-A)E_x^{\mathrm{PBE},\mathrm{SR}}\left(\omega \right)+E_x^{\mathrm{PBE},\mathrm{LR}}\left(\omega \right)+E_c^{\mathrm{PBE}},$$

(3)

where SR and LR mean short and long-range components, respectively. As the Coloumb operator is attenuated by the complementary error function for the short-range component and by the error function for the long-range component [54,55], there is also an additional range separation parameter $\omega$ due to erf($\omega r$). A very popular HSE06 functional is often defined when A = 0.25 and $\omega$ = 0.11 ${\mathrm{Bohr}}^{-1}$ [44,53].

We begin our analysis of STO properties with the more standard and widely used hybrid functionals to effectively compare with the literature data in Section 5.1. For this, we have chosen the three-parameter global B3LYP (B = 0.90, C = 0.81, and A = 0.20 in Equation (1)) functional, one-parameter global PBE0 functional (A = 0.25 in Equation (2)) and range-separated HSE06 functional based on the global PBE0 functional (A = 0.25 in Equation (3)). It turned out that the hybrid B3LYP, HSE06, and PBE0 functionals are among the most popular hybrid functionals for calculating perovskite systems [56,57]. In more recent DFT studies, it has been demonstrated that the next category (rung) of functionals on the Jacob’s ladder [55,58], i.e., strongly constrained and normed meta-GGA SCAN functional [59], combined with the hybrid functional idealogy provides us with certain improvements in oxides properties and calculation efficiency. Therefore, the present study intends to compare the three hybrid density functionals in Section 5.2 and Section 5.3, namely meta-GGA ${\mathrm{r}}^2$SCANx functional [27], range-separated HSEx functional, and PBEx functional. In the present study, all three functionals (i.e., ${\mathrm{r}}^2$SCANx, PBEx, and HSEx) suppose adjusting the exact exchange A amount to the indirect band gap value STO of 3.40 eV in the cubic phase. Notice also that variations in the exact HF exchange were considered in previous studies in the literature for perovskites [48,49,56] and other materials [27]. In the present study, however, such an approach gave us consistent lattice parameters and band gaps with the experiments even though the value of 3.40 eV for the indirect band gap is slightly increased compared to the experiments (see discussion in Section 5.2).

3.2. DFT Calculation Input Parameters

We consider STO in the cubic $Pm\overline{3}m$ (SG 221) and tetragonal $I4/mcm$ (SG 140) phases that have been experimentally and theoretically well studied. In the $Pm\overline{3}m$ model, the atoms occupy the following Wyckoff position (WP)

$$\begin{array}{ccc}\hfill \mathrm{Sr}:1\mathrm{b}& & (1/2,1/2,1/2),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \mathrm{Ti}:1\mathrm{a}& & (0,0,0),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \mathrm{O}:3\mathrm{d}& & (1/2,0,0).\hfill \end{array}$$

(6)

The tetragonal $I4/mcm$ model is obtained by the transformation $\left(\begin{array}{ccc}1& 1& 0\\ \overline{1}& 1& 0\\ 0& 0& 2\end{array}\right)$ of cubic lattice parameters without coordinate translation. In the $I4/mcm$ the WP of atoms split into [60]

$$\begin{array}{ccc}\hfill \mathrm{Sr}:4\mathrm{b}& & (0,1/2,1/4),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \mathrm{Ti}:4\mathrm{c}& & (0,0,0),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill \mathrm{O}:4\mathrm{a}& & (0,0,1/4),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \mathrm{O}:8\mathrm{h}& & (0.25+\delta ,0.25-\delta ,0).\hfill \end{array}$$

(10)

The oxygen 8 h WP is described by a single free parameter $\delta$, which is connected with the antiphase rotation of oxygen octahedra, ${\theta }_z$, i.e., AFD mode,

$${\theta }_z=arctan\left(4\delta \right).$$

(11)

In the absence of rotation, $\delta =0$, the tetragonal model coincides with the cubic one.

For the calculations with SGs $Pm\overline{3}m$ and $I4/mcm$, we use SH12 = 12 × 12 × 12 and SH8 = 8 × 8 × 8 Monkhorst–Pack k-point mesh [61], respectively. The chosen k-point mesh for the tetragonal phase is related to the fact that its primitive unit cell is given by the equal lattice parameters in all three directions. On the contrary, the primitive and conventional unit cells coincide in the cubic STO. The SCF convergence threshold for the total energy is set to ${10}^{-10}$ Hartree [44]. Integration is performed on a predefined pruned grid consisting of 99 radial and a maximum of 1454 angular points (XXLGRID) (to achieve accurate convergence of geometry); DFT density and DFT grid weight tolerances are kept as 8 and 16, respectively, [44]. Such highly accurate calculations are important for the calculations of forces and consequent phonon calculations, as STO obeys a non-trivial competition between the AFD and ferroelectric instabilities. Even though we do not aim at explaining the ferroelectric instability often predicted by DFT calculations (e.g., [14,43] and references therein), the comparison of BSs enforces us to obtain very accurate phonon frequencies.

Complete geometry optimization is performed until the energy difference between the two steps is less than the threshold (TOLDEE) ${10}^{-10}$ Hartree, root-mean-square of the gradient (TOLDEG) and displacement (TOLDEX) are $0.00003$ Hartree/Bohr and $0.00012$ Bohr, respectively, using no trust radius to limit displacement (NOTRUSTR). For an accurate comparison of system energy, we use the FIXINDEX option [44].

Tight tolerance set T16 (16 16 16 16 32) for Coulomb and exchange sums (five TOLINTEG parameters) is selected as a default in our calculations. This allows us to achieve the $R_4^{+}$ mode frequency saturation, see

Table 1. ${\mathrm{SrTiO}}_3$ soft mode $R_4^{+}$ (calculated phonon frequency in ${\mathrm{cm}}^{-1}$) dependence on TOLINTEG and SHRINK calculation parameters for the ${\mathrm{r}}^2$SCAN15 (15%HF) functional and BS1 in the $Pm\overline{3}m$ model. Due to the used supercell, the Monkhorst–Pack k-point mesh setting corresponds to $I4/mcm$ model.

TOLLINTEG Set	SH08	SH10	SH12	SH14
T12 (12:12:12:12:24)	1	$2i$	$1i$	0
T14 (14:14:14:14:28)	$5i$	$3i$	$5i$	$5i$
T16 (16:16:16:16:32)	$11i$	$11i$	$11i$	$11i$
T18 (18:18:18:18:36)	$10i$	$9i$	$9i$	$9i$

. The phonon calculations are performed using the supercell (direct) method implemented in CRYSTAL23 [44]. To perform the phonon calculations at the $\Gamma$ and R points of the cubic STO, we use the supercell due to $\left(\begin{array}{ccc}1& 1& 0\\ 1& 0& 1\\ 0& 1& 1\end{array}\right)$ transformation [7].

Notice that the bulk modulus was calculated according to a Voigt–Reuss–Hill scheme as implemented in CRYSTAL23 [44].

4. Visual Representation and the Local Basis Set Construction Method

The atomic orbitals, $\phi$, in CRYSTAL23, consist of shells, $\lambda$ (not to be confused with the principal quantum number), that are defined as normalized linear combinations of normalized real solid harmonic, $X_l^m\left(\mathbf{r}\right)$, Gaussian-type functions, $G({\alpha }_i^{\lambda };\mathbf{r})=exp(-{\alpha }_i^{\lambda }{\mathbf{r}}^2)$,

$$\begin{array}{ccc}\hfill {\phi }_{\lambda }^{lm}& =& N_{\lambda }\sum _i{d}_i^{\lambda }{c}_i^{lm}{X}_l^m\left(\mathbf{r}\right)G({\alpha }_i^{\lambda };\mathbf{r}),\mathrm{where}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill l& =& 0\left(s\right),1\left(p\right),2\left(d\right),3\left(f\right),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill m& =& -l\dots l,\hfill \end{array}$$

(14)

and $c_i^{lm}$ and $N_{\lambda }$ are normalisation coefficients, $d_i^{\lambda }$ are contraction coefficients, and index i runs over primitive functions of contractions [44].

We will construct the all-electron BSs for the Sr, Ti, and O atoms (denoted as BS1) for STO in the cubic $Pm\overline{3}m$ model from scratch. The steps are outlined in Figure 1. To start, we select a trial number of shells, $\lambda$, for each orbital type l (refer to Equation (13)) for each atom type. We construct an initial set of uncontracted nonredundant exponents, ${\alpha }_1^{\lambda }$, with a consecutive exponent ratio of approximately 2.5 for a clean start. (A more accurate starting guess for ${\alpha }_1^{\lambda }$ can be obtained by noting similar patterns for each orbital type independently on atom type in ratio plots displayed in Figure 2 (middle) or Supplementary Figure S1 (middle)). These exponent values are then adjusted to minimize the total system energy using the ${\mathrm{r}}^2$SCANx functional with the amount of exact exchange A = 0.15. For uncontracted exponents ($i=1$), the contraction coefficients $d_1^{\lambda }$ are inconsequential due to the normalization coefficient $N_{\lambda }$ (see Equation (12)). This straightforward approach allows us to obtain optimized exponent values that depend solely on the number of selected exponents and are independent of their initial values.

During the optimization process, we reduce the number of exponents if linear dependence is observed (i.e., the exponent ratio becomes smaller than 2). We also check if adding new exponents could reduce the total system energy (Figure 1). Optimally, all types of uncontracted exponents for all atoms should be determined simultaneously to prevent potential instabilities. For instance, an overestimated number of O s-type exponents might destabilize the Sr and Ti exponents. However, it is important to note the CRYSTAL23 limitations, allowing for a maximum of 31 shells (uncontracted exponents) per atom in bulk calculations (this limit is absent in molecule calculations). If the number of exponents exceeds this threshold, self-consistent procedures may need to be implemented.

Once uncontracted nonredundant exponent number and their optimized values ${\alpha }_1^{\lambda }$ are obtained, we keep them unchanged. We also do not reoptimise them after creating contractions, as discussed below. Visual representation comes into play during and after the optimization process to analyze the obtained exponents conveniently and simultaneously for all three atoms. However, let us consider all nonredundant (due to contractions) Sr s-type exponents as an example, see Supplementary Table S1, and abbreviate them as ${\alpha }_1^{\lambda }={\alpha }_i$, where here index $i\equiv \lambda$ and it runs over all uncontracted shells $\lambda$, Equation (12). Then the consecutive exponent ratio plot ${\alpha }_1/{\alpha }_2,{\alpha }_2/{\alpha }_3,\dots$ demonstrates a smooth decrease up to $i=11$, see middle graph in Figure 2a. Further, at $i=12,15,17$, indexes exponent ratio values have small maximums. These changes for $i>11$ are better seen when plotting the angle ${\varphi }_i$. Here ${\varphi }_i$ is defined as the angle among each of three consecutive points (${\alpha }_{i-1}/{\alpha }_i$, ${\alpha }_i/{\alpha }_{i+1}$,${\alpha }_{i+1}/{\alpha }_{i+2}$) to identify maximums in the middle graph of Figure 2a if any. So, the angle behavior changes smoothly in the absence of maximums (bottom graph in Figure 2a) and shows pronounced jumps at the maximum points. Therefore, the angle calculation appeared to bring in helpful information for further BS optimization. In a first approximation, it could be assumed that these maximums separate the exponent groups describing $1s$, $2s$, $3s$, and $4s$ electrons. However, the actual situation is more complicated since all uncontracted exponents could contribute to all electron descriptions.

There are numerous ways to create contractions of exponents ${\alpha }_i$. We follow Occam’s razor principle and retain only the most relevant exponents in contractions. Thus, groups of exponents with contraction coefficients $d_i^{\lambda }$ are formed to describe, e.g., $1s$, $2s$ and partly $3s$ Sr electrons and only the contraction coefficients are optimized (Figure 1). Their ratios are shown in the top Figure 2a. As discussed earlier, the same exponent could contribute to the description of several electrons and, therefore, appear in several contractions with different contraction coefficients, e.g., Sr ${\alpha }_i$ exponents with the indexes $i=7\dots 12$ appear in both shells $\lambda =1$ and $\lambda =2$. The Mulliken charge of each shell $q_{\lambda =1}=2.0$, $q_{\lambda =2}=2.0$, and $q_{\lambda =2}=1.1$$e^{-}$ allows us to assess the contribution of each shell in the description of the corresponding electrons. Shells $\lambda =1$ and 2 carry full charges of $1s$ and $2s$ orbitals, while $\lambda =3$ shell contains only a fraction of the expected $3s$ shell charge. Similarly, we create and analyze the Sr p- (Figure 2b), d- (Figure 2c), f-type functions. Exponents and contraction coefficients for Ti and O ions are given in Supplementary Figure S1. Furthermore, we give full BS1 in Supplementary Table S1 for ease of reference.

Although the total number of exponents and contractions could be selected slightly differently, we find that the

$$\begin{array}{ccc}\hfill \mathrm{Sr}:& & 29s23p10d1f/7s5p4d1f\left[20s14p10d1f\right],\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \mathrm{Ti}:& & 25s14p8d1f/5s4p4d1f\left[17s10p8d1f\right],\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \mathrm{O}:& & 16s7p1d/4s4p1d\left[12s7p1d\right]\hfill \end{array}$$

(17)

combinations accurately describe STO properties, where on the left side of the slash, we give the total number of Gaussian-type functions, while on the right side, we give the number of shells, $\lambda$, for each shell type, l. The nonredundant exponent number for each shell type is given in square brackets. For example, Sr basis contains 29 s-type basis functions (from which 20 are nonredundant) that are arranged in three contractions, see Figure 2a, and four uncontracted exponents leading to a total of 7 s-shells, Equation (15).

For a different model system, the BS reoptimization is straightforward. First, the contractions are removed, and then all nonredundant exponents are reoptimized. Second, the contractions are created with new exponents. Finally, only contraction coefficients are reoptimized.

The basis set BS1 is compared with three more BSs taken from the literature. Firstly, for the Ti and Sr atoms, we use the BS built on CRENBL small-core relativistic effective core potentials [62,63] and the all-electron O BS that is abbreviated as PBE0 optimized BS + ${\mathrm{Sr}}_d$ in Ref. [7] (denoted here as BS2, we give BS2 in Supplementary Table S1 for ease of reference).

Secondly, we consider the BS proposed for STO by Piskunov et al. [12] with the Hay-Wadt effective small core [64] potential-based BSs for Sr and Ti, and all electron BS for O (denoted as BS3, additional information on it can be found in Supplementary Table S1). The extension of oxygen (from BS3) set added for two Gaussian d-exponents could be found in [65,66,67] where the same oxygen BS is considered for the two limiting cases of the ${\mathrm{O}}_2$ molecule and ${\mathrm{SrTiO}}_3$ bulk crystal.

Thirdly, we use the pob-TZVP-rev2 BS for Sr that is built using the full-relativistic effective core potential taken from the ECP pseudopotential database [33,68], while for Ti and O, we use all-electron pob-TZVP-rev2 BSs [32] (denoted as BS4, additional information on it can be found in Supplementary Table S1).

It is instructive to apply the visual representation to analyze different BSs. For a comparison with BS1, we select oxygen p- and titanium d-type exponents from BS2-4 since they contain only a single contraction and, thus, closely resemble the BS1 construction method (Figure 3). Despite being of different lengths and obtained with different optimization procedures, the BSs exponents’ ratios demonstrate remarkable similarity. At the same time, their absolute values depend on the number of exponents and might differ considerably. Thus, the ratio of O p-type exponents at $i=1$ is ${\alpha }_1/{\alpha }_2\approx 4.6$ for BS1, where ${\alpha }_1$ and ${\alpha }_2$ values are 165.29 and 35.729 ${\mathrm{Bohr}}^{-2}$, respectively, while for BS3 they are 49.43 and 10.47 ${\mathrm{Bohr}}^{-2}$, respectively. Moreso, we find that this p-type exponents ratio is almost independent of the type of BS optimization (BS1-4 in Figure 3a) and weakly dependent on atom type, since for BS1 Sr p-type exponents, we obtain $\approx 4.2$ (Figure 2b) and for Ti p-type exponents we obtain $\approx 4.3$ (Supplementary Figure S1b). We also note that the smallest p-type value saturates at $\approx 0.17$ ${\mathrm{Bohr}}^{-2}$ for BS1 and BS3, which indicates that the optimal number of O p-type exponents is reached.

Titanium d-type exponents at $i=1$ cluster around ${\alpha }_1/{\alpha }_2\approx 3.6$ for BS1 with exponent values of 249.06 and 68.579 ${\mathrm{Bohr}}^{-2}$, while for BS3, those are 21.429541 and 6.08722431 ${\mathrm{Bohr}}^{-2}$, respectively, (Figure 3b). Similarly to the p-type exponents case, the ratio of d-type at $i=1$ for titanium is close to that of strontium ${\alpha }_1/{\alpha }_2\approx 3.3$ (Figure 2c). The tail part of the exponents’ ratios curve of BS1 and BS2 have a similar shape (Figure 3b); they have local maximums at $i=6$ and 3, respectively. For BS4, we find that at $i=4$, the ratio is ${\alpha }_4/{\alpha }_5\approx 1.6$, formed by the last exponent in the contraction ${\alpha }_4=1.35885903030$ and the first uncontracted exponent ${\alpha }_5=0.86367514000$ ${\mathrm{Bohr}}^{-2}$. In case of problems, these exponents should be crosschecked since their ratio is less than the recommended safe threshold of 2 to avoid linear dependence.

The visual representation thus has several advantages. It allows us to (i) detect the possible linear dependency problems if the exponent ratio, ${\alpha }_i/{\alpha }_{i+1}$, becomes smaller than 2, (ii) inspect and compare BS with a different number of exponents, and (iii) identify invariant relations between BS exponents. Additionally, it could facilitate the development of robust methods for constructing self-consistent BS, allowing complete BS optimization for each particular material, functional, or calculation accuracy.

5. Simulation Results and Discussion

Our analysis of the main results is divided into three Subsections. Before discussing the exact exchange adjustment A for each BS1-4 and different chosen hybrid functional combinations, we calculate the basic STO properties with these BSs and more popular functionals in Section 5.1. We use three different well-known hybrid functionals, namely B3LYP, PBE0, and HES06, to evaluate the structural and electronic properties of STO in two experimentally established cubic, $Pm\overline{3}m$, and tetragonal, $I4/mcm$ (AFD), phases (

Table 2. ${\mathrm{SrTiO}}_3$ pseudo-cubic lattice constants, $a_0$, $b_0$, $c_0$, tetragonality ratio, $c_0/a_0-1$, rotation angle of the oxygen octahedra around the z axis, ${\theta }_z$, and in-/direct band gaps, $E_g^i$/$E_g^d$. Energy difference, $\Delta E^{STO}$, relative to the cubic $Pm\overline{3}m$ model per formula unit. All the properties are calculated for BS1-4 and the standard hybrid DFT functionals.

Parameters	B3LYP (20%HF)				PBE0 (25%HF)				HSE06 (25%HF)
	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4
$Pm\overline{3}m$(SG 221)
$a_0$, Å	$3.935$	$3.937$	$3.938$	$3.925$	$3.899$	$3.902$	$3.901$	$3.891$	$3.900$	$3.903$	$3.903$	$3.892$
$E_g^i$, eV	$3.48$	$3.45$	$3.72$	$3.53$	$3.97$	$3.92$	$4.22$	$4.03$	$3.27$	$3.22$	$3.50$	$3.32$
$E_g^d$, eV	$3.84$	$3.77$	$4.05$	$3.90$	$4.32$	$4.24$	$4.54$	$4.39$	$3.66$	$3.57$	$3.86$	$3.72$
$I4/mcm$(SG 140)
$a_0$, Å	$3.931$	$3.931$	$3.936$	$3.925$	$3.898$	$3.898$	$3.901$	$3.891$	$3.899$	$3.900$	$3.903$	$3.892$
$(c_0/a_0-1){10}^3$	$2.6$	$3.2$	$1.1$	$0.3$	$0.6$	$2.2$	$0.2$	$0.1$	$0.8$	$2.2$	$0.3$	$0.1$
${\theta }_z$, °	$3.6$	$4.2$	$2.4$	$0.9$	$1.7$	$3.3$	$0.9$	$0.0$	$2.0$	$3.3$	$1.2$	$0.0$
$E_g^d$, eV	$3.53$	$3.51$	$3.74$	$3.54$	$3.98$	$3.96$	$4.22$	$4.03$	$3.28$	$3.26$	$3.50$	$3.32$
$\Delta E^{STO}$,  meV	$-0.90$	$-2.98$	$-0.43$	$-13.90$	$-0.13$	$-1.22$	$-0.04$	$-0.00$	$-0.15$	$-1.26$	$-0.05$	$-0.00$

). From numerous STO properties, we have chosen the most relevant and important ones for our present “DFT method” + “BS” analysis and comparison with the experiments. We aimed to present an analysis of several STO properties simultaneously, including the pseudocubic lattice constant, $a_0$, in-/direct band gaps, $E_g^i/E_g^d$, tetragonality ratio, $c_0/a_0-1$, and, correlated with it, angle of rotation of the O octahedra around the z-axis, ${\theta }_z$. In Section 5.2, the same properties are then re-evaluated due to the adjusted A-values for PBEx (Equation (2)), HSEx (Equation (3)), and ${\mathrm{r}}^2$SCANx (Equation (2)). Section 5.3 also analyzes the most relevant phonon frequencies at the cubic phase’s $\Gamma$- and R-points and $\Gamma$-point only of the tetragonal phase.

5.1. Basic STO Properties Calculated with B3LYP, PBE0, HSE06

As was expected, the obtained results are scattered, as shown in Table 2, and general trends are hardly observable. In other words, each combination of functional and BS provides different values for the structural properties chosen. Notice that the experimental values for the comparison are given in Table 3. It may be said that if one parameter agrees well with the experiment, the others deviate. This contradictory behavior may be due to the poor quality of the BS. Therefore, the BS optimization needed to be reconsidered and done in a more sophisticated manner, as is already pointed out above, and is the main subject of the present study.

Our initial analysis already demonstrates a significant trend. Specifically, the B3LYP functional noticeably overestimates the lattice constant and band gap values for both phases. This finding aligns with previous LCAO calculations on STO using the B3LYP functional [12]. In contrast, the PBE0 and HSE06 functionals provide more accurate lattice parameters. Additionally, the HSE06 functional yields good band gap values, while the PBE0 functional significantly overestimates them in both phases. The oxygen octahedra rotation angle (${\theta }_z$ in Table 2) in the SG $I4/mcm$ also demonstrates differences between the HSE06 and PBE0 functionals. Notably, the lattice parameters calculated using the HSE06 functional closely match those obtained from plane wave calculations using the HSE functional with the screening parameter $\omega$ = 0.15 ${\mathrm{Bohr}}^{-1}$ and A = 0.25 [14]. Furthermore, the hybrid functionals better describe the rotation angle ${\theta }_z$ equal to 1.2–3.3° for in the present LCAO calculations depending on the BS and 2.63° in the plane wave calculations [14] with the HSE-like functionals compared to 6.47° in the DFT + U + V calculations Ref. [16].

Various effects can be observed due to differences in the BSs. The “tetragonality” quantities characterize the degree of tetragonality and stabilization of the AFD phase for different calculation methods. It would be particularly beneficial to estimate them to analyze the effect of BS and functional. So, the tetragonality ratio $c_0/a_0-1$ and rotation angle ${\theta }_z$ appeared to be sensitive to the choice of BS. Moreover, the stabilization energy of the AFD phase ($\Delta E^{STO}$ in Table 2) depends on the BS and functional used, resulting in variations from almost 0.00 to 1.26 meV/f.u. in the LCAO calculations. Comparatively, our calculated values for $\Delta E^{STO}$ are slightly smaller than in the plane wave calculations with the HSE functional (2 meV in [14]; however, it is not indicated directly in [14] whether the energy difference is given per f.u.). Notably, BS4 does not correctly reproduce the AFD phase with PBE0 and HSE06 functionals. Another significant finding is that the band gap of the AFD phase, calculated with the HSE06 functional, is underestimated by 0.08–0.14 eV for BS1-2 and BS4. Only calculations using BS3 with the HSE06 functional resulted in an overestimation of the band gap of the AFD phase by 0.10 eV compared to the experimental data. Moreover, it is seen from the results in Table 2 that the calculated standard STO properties do not allow the selection of an optimal BS and functional combination since each functional gives slightly different predictions and transferability (invariance of the results for different functionals) is absent. Therefore, to highlight the differences, we refine our calculations by adjusting the exact HF exchange amount in Section 5.2 and considering vibrational properties in Section 5.3.

Overall, we exclude the B3LYP functional from further considerations in Section 5.2 and utilize the PBEx and HSEx functionals in subsequent calculations together with the meta-GGA hybrid ${\mathrm{r}}^2$SCAN-type functional and adjusted A-values.

5.2. Basic STO Properties Calculated with Adjusted Exact HF Exchange and $r^2$SCANx Functionals

As mentioned in Section 3.1, we adjust the exact Hartree-Fock exchange contribution, A, in the selected three functionals (${\mathrm{r}}^2$SCANx, PBEx, and HSEx) so that the indirect gap in the $Pm\overline{3}m$ model matches the $E_g^i=3.40$ eV for better functional comparison. This adjustment is mainly influenced by the fact that the band gap is more responsive to the changes of A-value than the lattice parameters. The adjusted A-value is similar for BS1-2 and BS4 within the same functional. However, for BS3 we find systematically lower A values compared to the other three functionals: $A\approx 0.11$ (BS3) vs. 0.14–0.15 (BS1-2, BS4) for ${\mathrm{r}}^2$SCANx, ≈0.16 (BS3) vs. 0.18–0.19 (BS1-2, BS4) for PBEx, and ≈0.24 (BS3) vs. 0.26–0.27 (BS1-2, BS4) for HSEx (Table 3). This difference can be attributed to the fact that BS3 was optimized in pure HF calculations with $A=1$ (Equation (2)) [12].

The lattice parameters and other properties for the two STO phases, as shown in Table 3, are calculated using the adjusted A for each functional and BS. Indeed, we did not observe a significant impact of A on the lattice parameters for either the standard hybrid functionals (PBE0 and HSE06) in Table 2 or for the hybrid functionals with the adjusted A (PBEx and HSEx) in Table 3. Typically, the DFT + U-like calculations tend to overestimate the lattice constant (e.g., 3.96 Åin the DFT + $U^{d,p}$ + $J^{d,p}$ [15] for the cubic phase) when the Hubbard parameters are adjusted to match the band gap. However, the current hybrid DFT calculations suggest a much better agreement between the lattice constant and the experiments. Nonetheless, the bulk modulus B for the cubic phase is slightly overestimated with the hybrid functionals (see also Table 1 in [15]) regardless of LCAO or plane wave basis. Additionally, it is evident that PBEx results in a smaller B (and somewhat larger $a_0$) regardless of BS, though it still overestimates the experimental value.

The “tetragonality” quantities, i.e., $\Delta E^{STO}$, $c_0/a_0-1$ and ${\theta }_z$, increase for the PBEx functional after the adjusted-A value is applied compared with PBE0. Conversely, these quantities decrease for the HSEx functional except ${\theta }_z$-value for BS2. It is likely due to a reduced exact exchange for PBEx ($A\approx 0.19$ (BS1)) and the opposite effect for HSEx ($A\approx 0.27$ (BS1)). Notice that HSE06 and PBE0 are characterized by A = 0.25. The basis set BS4 did not show significant improvements with the adjusted A-value regarding the stabilization of the AFD phase; all the corresponding “tetragonality” quantities are underestimated. The shortcomings of BS4 are well understood and explained in Section 4 and Figure 3 with the help of the visualization technique as proposed in the present work.

The most striking result is the almost insensitivity of STO properties to the functional change from ${\mathrm{r}}^2$SCANx to HSEx for BS1-3 despite significantly different amounts of exact exchange. Our approach with the adjusted A-value improves the transferability of BS1-3 between the two functionals Table 3, which was not the case for the standard functionals in Table 2. However, the obtained picture could be completed only with the inclusion of calculated phonon properties. Below, we present the analysis of the phonon frequencies calculated at the $\Gamma$- and R-points for BS1-4 with PBEx, HSEx, and ${\mathrm{r}}^2$SCANx and adjusted A-values.

5.3. Vibrational Properties of STO

We focus on the phonon modes and frequencies below approximately 100 ${\mathrm{cm}}^{-1}$, which are the most challenging to calculate. In Table 3, we provide different notations of irreducible representations for the phonon frequencies to make it clearer for the reader and simplify the comparison with the literature. Thus, we use the form Bethe (Mulliken) [69] for the notations for the tetragonal $I4/mcm$ model. For the symmorphic SG $Pm\overline{3}m$, we also include the corresponding Bouckaert-Smoluchowski-Wigner (BSW) notations [69,70] in addition in the form Bethe (Mulliken, BSW).

According to the group-theoretical analysis, the cubic-tetragonal phase transition occurs due to the unstable mode $R_4^{+}$ ($T_{1g}$, $R_{15}^{\prime }$) [71] in the cubic phase [7]. It is a three-fold degenerate mode (due to the x, y, and z directions) associated with the oxygen octahedra rotations around those axes. The unstable (imaginary frequency) phonon mode $R_4^{+}$ is obtained consistently with BS1 in combination with all three considered functionals (${\mathrm{r}}^2$SCANx, PBEx, and HSEx) (Table 3), and that agrees well with plane wave calculation results [7]. However, for other BSs, the mode stability depends on the selected functional. For example, the unstable mode $R_4^{+}$ is obtained with BS2 only with the ${\mathrm{r}}^2$SCANx and PBEx functionals. Albeit the phonon calculations with HSEx predict marginal stability of $R_4^{+}$ for BS2-4, the relaxation in the tetragonal $I4/mcm$ model leads to the appearance of tetragonality and rotated oxygen octahedra structure with a non-zero ${\theta }_z$ for BS2-BS3 and a zero ${\theta }_z$ for BS4 (Table 3). Similarly, the phonon calculations with BS3 predict unstable mode $R_4^{+}$ only with ${\mathrm{r}}^2$SCANx, but a geometry relaxation with HSEx and PBEx in $I4/mcm$ leads to oxygen octahedra rotations. All the functionals predict stable mode $R_4^{+}$ with BS4, but only ${\mathrm{r}}^2$SCANx gives rotated oxygen octahedra geometry.

Relaxing the unstable mode $R_4^{+}$ in the order parameter direction (a,0,0) leads to the tetragonal $I4/mcm$ phase, where $R_4^{+}$ splits into two Raman active modes, $R_4^{+}\to A_{1g}+E_g$ [72]. The non-degenerate $A_{1g}$ mode (associated with O ion displacements) and the double-degenerate $E_g$ mode (connected with O-Sr displacements) are stable in all our calculations. Their frequencies relation due to $f_{A_{1g}}$ > $f_{E_g}$ agrees with the experiments for BS1-3 with all considered functionals (Table 3). Experimental results demonstrate that below the phase transition temperature, both $A_{1g}$ and $E_g$ modes gain intensity on cooling and finally harden in the low-temperature limit [73,74], while quantitative comparison with DFT calculations is challenging due to scattered experimental results, we note that BS1 and BS3 (with ${\mathrm{r}}^2$SCANx and HSEx) align well with the experiments [73].

Let us consider the triple degenerate mode ${\Gamma }_4^{-}$ ($T_{1u}$, ${\Gamma }_{15}$) in the cubic phase that corresponds to ferrodistortive (FD) displacements of both Sr and Ti and all O ions in opposite directions along x, y, and z directions. According to the group-theoretical analysis [7,72] the infrared active ${\Gamma }_4^{-}$ mode splits into two active infrared modes in the tetragonal $I4/mcm$ phase ${\Gamma }_4^{-}\to A_{2u}+E_u$. The double degenerate $E_u$ mode corresponds to FE displacements when all Sr and Ti and all O ions move in opposite directions along x and y axes. However, the non-degenerate $A_{2u}$ mode corresponds to FE displacements along the z axis. Experiments show that temperature lowering decreases intensities of $A_{2u}$ and $E_u$ modes and increases deviation from the Curie-Weiss law below the phase transition, which predicts the condensation out of those modes and thus the appearance of a new FE phase. Instead, the low-temperature experiments demonstrate the saturation of these two modes due to quantum fluctuations [74,75].

Contrary to the unstable mode $R_4^{+}$ that leads to the $I4/mcm$ model, the expected nature of mode ${\Gamma }_4^{-}$ in the DFT calculations still needs to be clarified since quantum fluctuations of ions are not included in the DFT calculations. As we focus on the BS construction and optimization in the present study, we leave this question out for future studies. However, we do analyze a consistency between the calculated phonon frequencies for ${\Gamma }_4^{-}$ in the cubic phase and $A_{2u}$, $E_u$ in the tetragonal phase to compare the calculation methods. We expect the sign conservation for accurate BS for these modes since FE mode is suppressed in both $Pm\overline{3}m$ and $I4/mcm$ models. Indeed, we have observed this behavior in the present study for all three functionals only for BS1 predicting unstable ${\Gamma }_4^{-}$. Such behavior coincides with our earlier study using the standard PBE functional in the plane wave calculations [7]. BS2 demonstrated consistency with ${\mathrm{r}}^2$SCANx and HSEx (except PBEx) functionals predicting stable ${\Gamma }_4^{-}$ mode. BS3 showed consistency for all functionals; however, it predicted different natures of ${\Gamma }_4^{-}$—stable with ${\mathrm{r}}^2$SCANx, but unstable with HSEx and PBEx.

Overall, from the analysis of calculated STO properties, which included the basic structure properties, electronic band gaps, and phonon frequencies, we propose using BS1 combined with the functionals ${\mathrm{r}}^2$SCAN15 or HSE27 for accurate STO calculations within the LCAO formalism. Notice that the reduced amount of HF exchange of 15% is very much in line with the ${\mathrm{r}}^2$SCANh (amount of the HF exact exchange 10%) functional suggested earlier in [27].

Table 3. Comparison of ${\mathrm{SrTiO}}_3$ characteristic parameters obtained with the parametrized ${\mathrm{r}}^2$SCANx, HSEx, and PBEx functionals with other DFT calculations and experiments. Adjusted HF exact exchange, A, for each BS and functional, providing an indirect gap of $\sim 3.40$ eV for the T16 set in the SG $Pm\overline{3}m$. ${\mathrm{SrTiO}}_3$ pseudo-cubic lattice constants, $a_0$, $b_0$, $c_0$, tetragonality ratio, $c_0/a_0-1$, rotation angle of the oxygen octahedra around the z axis, ${\theta }_z$, bulk modulus, B, and in-/direct band gaps, $E_g^i$/$E_g^d$. In the cubic $Pm\overline{3}m$ model we consider the ${\Gamma }_4^{-}$ ($T_{1u}$) and $R_4^{+}$ ($T_{1g}$) phonon modes at $\Gamma$ and R points, respectively. In the tetragonal $I4/mcm$ model we consider two Raman modes ${\Gamma }_5^{+}$ ($E_g$), ${\Gamma }_1^{+}$ ($A_{1g}$), and two infrared modes ${\Gamma }_5^{-}$ ($E_u$), ${\Gamma }_2^{-}$ ($A_{2u}$). We provide Bethe (Mulliken and for $Pm\overline{3}m$ also BSW) notations of irreducible representations [69]. Energy difference, $\Delta E^{STO}$, relative to the cubic $Pm\overline{3}m$ model per formula unit.

Parameters	${\mathrm{r}}^2$SCANx (A)				HSEx (A)				PBEx (A)				Exp
	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4
A,  %	$14.7$	$13.9$	$10.8$	$14.1$	$27.1$	$27.8$	$23.5$	$26.2$	$18.9$	$19.3$	$16.2$	$18.3$
$Pm\overline{3}m$ (SG 221)
$a_0$,  Å	$3.895$	$3.901$	$3.901$	$3.888$	$3.897$	$3.899$	$3.905$	$3.890$	$3.908$	$3.910$	$3.915$	$3.901$	$3.895$ a
B,  GPa	202	203	200	207	198	202	196	204	191	196	191	198	$179$ b
${\Gamma }_4^{-}$ ($T_{1u}$, ${\Gamma }_{15}$),  ${\mathrm{cm}}^{-1}$	$101i$	112	64	44	$121i$	88	$48i$	48	$142i$	32	$83i$	$56i$	soft c
$R_4^{+}$ ($T_{1g}$, $R_{15}^{\prime }$),  ${\mathrm{cm}}^{-1}$	$12i$	$50i$	$35i$	24	$22i$	10	61	79	$41i$	$29i$	53	72	soft c
$E_g^i$,  eV	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.2\dots 3.3$ d
$E_g^d$,  eV	$3.77$	$3.75$	$3.76$	$3.78$	$3.79$	$3.75$	$3.76$	$3.80$	$3.76$	$3.71$	$3.74$	$3.77$	$3.7\dots 3.8$ d
$I4/mcm$ (SG 140)
$a_0$,  Å	$3.894$	$3.897$	$3.900$	$3.888$	$3.896$	$3.896$	$3.904$	$3.890$	$3.906$	$3.905$	$3.913$	$3.901$	$3.895$ a
$c_0$,  Å	$3.896$	$3.904$	$3.902$	$3.889$	$3.898$	$3.904$	$3.906$	$3.891$	$3.911$	$3.916$	$3.917$	$3.901$	$3.899$ a
$(c_0/a_0-1){10}^3$	$0.5$	$1.8$	$0.4$	$0.1$	$0.6$	$1.9$	$0.5$	$0.1$	$1.5$	$2.7$	$1.1$	$0.1$	$0.6$, $0.3$ e
${\theta }_z$, °	$1.6$	$3.2$	$1.7$	$0.9$	$1.6$	$3.1$	$1.7$	$0.0$	$2.7$	$3.7$	$2.4$	$0.0$	$2.1$, $1.4$ e
${\Gamma }_5^{+}$ ($E_g$),  ${\mathrm{cm}}^{-1}$	16	17	33	44	18	10	22	49	29	7	14	35	$40$c, $15$ f
${\Gamma }_1^{+}$ ($A_{1g}$),  ${\mathrm{cm}}^{-1}$	40	88	55	33	39	78	42	47	64	90	56	31	$52$c, $48$ f
${\Gamma }_5^{-}$ ($E_u$),  ${\mathrm{cm}}^{-1}$	$89i$	114	81	59	$124i$	74	$76i$	$31i$	$143i$	$32i$	$101i$	$80i$	$8$ g
${\Gamma }_2^{-}$ ($A_{2u}$),  ${\mathrm{cm}}^{-1}$	$100i$	121	85	68	$121i$	76	$78i$	$33i$	$141i$	$25i$	$102i$	$80i$	$17$ g
$E_g^d$,  eV	$3.41$	$3.43$	$3.41$	$3.41$	$3.41$	$3.43$	$3.41$	$3.40$	$3.43$	$3.45$	$3.42$	$3.40$	$3.4$ d
$\Delta E^{STO}$,  meV	$-0.10$	$-1.40$	$-0.24$	$-0.01$	$-0.07$	$-0.98$	$-0.07$	$-0.00$	$-0.54$	$-1.90$	$-0.38$	$-0.00$

^a Data from Ref. [76] extrapolated to 0 K. ^b Ref. [77]. ^c Two soft, temperature-dependent modes split below phase transition as ${\Gamma }_4^{-}\left(T_{1u}\right)\to A_{2u}+E_u$ and $R_4^{+}\left(T_{1g}\right)\to A_{1g}+E_g$. They could be approximated as $f\left(T_{1u}\right)=5.6{(T-31)}^{1/2}$ and $f\left(A_{1g}\right)=11.2{(132-T)}^{1/3}$, respectively, see Ref. [74]. ^d Refs. [6,11,78,79,80]. ^e Ref. [81] at $4.2$ and 77 K, respectively. ^f Ref. [73]. ^g Ref. [75].

Table 3. Comparison of ${\mathrm{SrTiO}}_3$ characteristic parameters obtained with the parametrized ${\mathrm{r}}^2$SCANx, HSEx, and PBEx functionals with other DFT calculations and experiments. Adjusted HF exact exchange, A, for each BS and functional, providing an indirect gap of $\sim 3.40$ eV for the T16 set in the SG $Pm\overline{3}m$. ${\mathrm{SrTiO}}_3$ pseudo-cubic lattice constants, $a_0$, $b_0$, $c_0$, tetragonality ratio, $c_0/a_0-1$, rotation angle of the oxygen octahedra around the z axis, ${\theta }_z$, bulk modulus, B, and in-/direct band gaps, $E_g^i$/$E_g^d$. In the cubic $Pm\overline{3}m$ model we consider the ${\Gamma }_4^{-}$ ($T_{1u}$) and $R_4^{+}$ ($T_{1g}$) phonon modes at $\Gamma$ and R points, respectively. In the tetragonal $I4/mcm$ model we consider two Raman modes ${\Gamma }_5^{+}$ ($E_g$), ${\Gamma }_1^{+}$ ($A_{1g}$), and two infrared modes ${\Gamma }_5^{-}$ ($E_u$), ${\Gamma }_2^{-}$ ($A_{2u}$). We provide Bethe (Mulliken and for $Pm\overline{3}m$ also BSW) notations of irreducible representations [69]. Energy difference, $\Delta E^{STO}$, relative to the cubic $Pm\overline{3}m$ model per formula unit.

Parameters	${\mathrm{r}}^2$SCANx (A)				HSEx (A)				PBEx (A)				Exp
	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4	BS1	BS2	BS3	BS4
A,  %	$14.7$	$13.9$	$10.8$	$14.1$	$27.1$	$27.8$	$23.5$	$26.2$	$18.9$	$19.3$	$16.2$	$18.3$
$Pm\overline{3}m$ (SG 221)
$a_0$,  Å	$3.895$	$3.901$	$3.901$	$3.888$	$3.897$	$3.899$	$3.905$	$3.890$	$3.908$	$3.910$	$3.915$	$3.901$	$3.895$ a
B,  GPa	202	203	200	207	198	202	196	204	191	196	191	198	$179$ b
${\Gamma }_4^{-}$ ($T_{1u}$, ${\Gamma }_{15}$),  ${\mathrm{cm}}^{-1}$	$101i$	112	64	44	$121i$	88	$48i$	48	$142i$	32	$83i$	$56i$	soft c
$R_4^{+}$ ($T_{1g}$, $R_{15}^{\prime }$),  ${\mathrm{cm}}^{-1}$	$12i$	$50i$	$35i$	24	$22i$	10	61	79	$41i$	$29i$	53	72	soft c
$E_g^i$,  eV	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.40$	$3.2\dots 3.3$ d
$E_g^d$,  eV	$3.77$	$3.75$	$3.76$	$3.78$	$3.79$	$3.75$	$3.76$	$3.80$	$3.76$	$3.71$	$3.74$	$3.77$	$3.7\dots 3.8$ d
$I4/mcm$ (SG 140)
$a_0$,  Å	$3.894$	$3.897$	$3.900$	$3.888$	$3.896$	$3.896$	$3.904$	$3.890$	$3.906$	$3.905$	$3.913$	$3.901$	$3.895$ a
$c_0$,  Å	$3.896$	$3.904$	$3.902$	$3.889$	$3.898$	$3.904$	$3.906$	$3.891$	$3.911$	$3.916$	$3.917$	$3.901$	$3.899$ a
$(c_0/a_0-1){10}^3$	$0.5$	$1.8$	$0.4$	$0.1$	$0.6$	$1.9$	$0.5$	$0.1$	$1.5$	$2.7$	$1.1$	$0.1$	$0.6$, $0.3$ e
${\theta }_z$, °	$1.6$	$3.2$	$1.7$	$0.9$	$1.6$	$3.1$	$1.7$	$0.0$	$2.7$	$3.7$	$2.4$	$0.0$	$2.1$, $1.4$ e
${\Gamma }_5^{+}$ ($E_g$),  ${\mathrm{cm}}^{-1}$	16	17	33	44	18	10	22	49	29	7	14	35	$40$c, $15$ f
${\Gamma }_1^{+}$ ($A_{1g}$),  ${\mathrm{cm}}^{-1}$	40	88	55	33	39	78	42	47	64	90	56	31	$52$c, $48$ f
${\Gamma }_5^{-}$ ($E_u$),  ${\mathrm{cm}}^{-1}$	$89i$	114	81	59	$124i$	74	$76i$	$31i$	$143i$	$32i$	$101i$	$80i$	$8$ g
${\Gamma }_2^{-}$ ($A_{2u}$),  ${\mathrm{cm}}^{-1}$	$100i$	121	85	68	$121i$	76	$78i$	$33i$	$141i$	$25i$	$102i$	$80i$	$17$ g
$E_g^d$,  eV	$3.41$	$3.43$	$3.41$	$3.41$	$3.41$	$3.43$	$3.41$	$3.40$	$3.43$	$3.45$	$3.42$	$3.40$	$3.4$ d
$\Delta E^{STO}$,  meV	$-0.10$	$-1.40$	$-0.24$	$-0.01$	$-0.07$	$-0.98$	$-0.07$	$-0.00$	$-0.54$	$-1.90$	$-0.38$	$-0.00$

^a Data from Ref. [76] extrapolated to 0 K. ^b Ref. [77]. ^c Two soft, temperature-dependent modes split below phase transition as ${\Gamma }_4^{-}\left(T_{1u}\right)\to A_{2u}+E_u$ and $R_4^{+}\left(T_{1g}\right)\to A_{1g}+E_g$. They could be approximated as $f\left(T_{1u}\right)=5.6{(T-31)}^{1/2}$ and $f\left(A_{1g}\right)=11.2{(132-T)}^{1/3}$, respectively, see Ref. [74]. ^d Refs. [6,11,78,79,80]. ^e Ref. [81] at $4.2$ and 77 K, respectively. ^f Ref. [73]. ^g Ref. [75].

6. Conclusions

Optimizing or at least double-checking all the exponents and contraction coefficients of Gaussian-type functions for each new material is essential. However, this is rarely done because of the lengthy and cumbersome BS optimization process and the resultant value dependence on many subjective choices made during optimization. Based on the invariant exponent ratio relationships established through BS visual representation, the proposed robust local BS construction method in LCAO calculations has significant practical implications. This method is demonstrated for STO due to its challenging properties for DFT calculations, making it an exciting subject of study. Thus, the obtained local BS1 and three other well-known local BSs are used to calculate the properties of STO in its cubic and tetragonal phases.

The calculated STO properties include basic structure properties: the so-called “tetragonality” properties, i.e., the tetragonality ratio and rotation angle, electronic band gaps and phonon frequencies at $\Gamma$- and R-points of cubic STO phase, and only $\Gamma$-point of tetragonal STO phase. The present study also thoroughly addresses the band gap problem using several hybrid DFT functionals (PBEx, HSEx, and ${\mathrm{r}}^2$SCANx) with the adjusted amount of exact exchange. By comparing all these properties for four BSs with the existing theoretical and experimental literature, we show that all electron BS1 - made with the suggested BS construction method - demonstrates results comparable with the experiments for several properties simultaneously as chosen and is transferable between the functionals HSEx and ${\mathrm{r}}^2$SCANx. The suggested amount of HF exact exchange for BS1 is 15% for ${\mathrm{r}}^2$SCANx and 27% for HSEx.

Notice that we have recently used a more elaborated BS construction and optimization method that included simultaneous all BS exponent and contraction coefficient optimization for ${\mathrm{Sr}}_2{\mathrm{FeO}}_4$ [40] and ${\mathrm{LiYF}}_4$ [41]. However, here we demonstrate that a more robust scheme provides accurate and consistent results for STO that are comparable with other state-of-the-art BSs. The advantage of the proposed robust BS construction scheme is that obtained BS can be easily decontracted (retaining only nonredundant exponents), nonredundant exponents optimized, and then back-contracted (with only contraction coefficient optimization) for each new material. Such a method allows users to verify all BS parameters for transferability to new materials (and also different functionals, calculation parameters) and thus increase the accuracy of DFT calculations.