Exact-two-component block-localized wave function: A simple scheme for the automatic computation of relativistic ΔSCF

Block-localized wave function (BLW) is a technique for constrained optimization of molecular orbitals in self-consistent field (SCF) calculations. In this method, molecular orbitals can be strictly localized within a molecular fragment or on the basis of an imposed symmetry,^1,2 providing a variational procedure to partition the total energy into specific components in energy-decomposition analysis (EDA).^3–5 BLW-based EDA has been widely used to understand the origin of Lewis acid–base bonding, cation–π interactions, and transition metal-ligand charge-transfer effects.⁶ Recently, we introduced a generalized approach, namely, generalized block-localized wave function (gBLW), to optimize a determinant wave function that can be specifically constrained in molecular orbital space,^7,8 which further expands the scope of its application in EDA^9–11 as well as local and charge-transfer excited states in molecular complexes.¹² 

The generalized block-localized wave function (gBLW) method can be used to optimize configuration states to be used in nonorthogonal configuration interaction, such as that used in effective valence-bond (VB) theory¹³ and multistate density functional theory (MSDFT),¹⁴ or approximations to excited states.⁸ Effective VB type of diabatic states can be useful to model chemical reactions such as the S_n2 reaction,^15,16 charge-transfer reactions,^17–21 and excited energy transfer processes in photoreceptor proteins.^22,23 In particular, Mo and co-workers investigated hyperconjugation,²⁴ aromaticity,²⁵ and the anomeric effect.²⁶ We have employed gBLW to optimize electronic excited states in the ΔSCF framework.⁸ The transformation of molecular orbitals into specific blocks in gBLW avoids the risk of variational collapse, and the computational cost is similar to that of the ground state SCF optimization by a direct inversion of the iterative subspace (DIIS) scheme that takes into account both block-localized and the full system densities; the method is complementary to other approaches such as the maximum overlap method^27–30 or time-dependent density functional theory.^30–33 

In this work, we extended the gBLW approach to the two-component relativistic electronic structure framework. In the two-component theory, the relaxation of spin axis allows for the inclusion of spin–orbit couplings in a variational approach. However, it also creates a challenging scenario for the gBLW as orbitals that belong to the same term (e.g., ²P_3/2 and ²D_5/2) can lead to variational collapse of the gBLW due to their spatial overlap. Additionally, since the two-component excited state manifold is more dense, it would be more difficult to guarantee that all of the states within a degenerate manifold be adequately represented. In this regard, it is necessary to construct all configuration states if this method is to be used to generate the determinant basis for nonorthogonal configuration interaction.

This article is organized as follows: In Sec. II, we introduce exact-two-component (X2C) wave functions along with the gBLW algorithm. We will explore strategies to allow the gBLW to obtain all excited states that are split due to spin–orbit coupling in the ΔSCF process. Then, in Sec. IV, we apply our method to atomic systems to explore their zero-field-splitting and the degree to which ΔSCF can maintain Kramers’ and total angular momentum symmetries naturally (i.e., without resorting to extraordinary algorithmic measures).

Throughout this work, we use the following notation:

• μ, ν index atomic spinor orbitals (AO’s);

• i, j index molecular orbitals (MO’s);

• A, B index blocks/subspaces.

In the exact-two-component (X2C) method,^34–53 the electronic and positronic components of the four-component Dirac equation are decoupled by a unitary transformation $U$ that block-diagonalizes the four-component Hamiltonian,

In Eq. (1), the matrix form of the Dirac four-component Hamiltonian is written in the restricted-kinetically balanced (RKB) condition.^54,55 Usually, only the two-component Hamiltonian corresponding to electronic solutions, $H+$⁠, is needed for determining the properties of chemical interest. In Eq. (1), $V$ and $T$ are the non-relativistic potential energy and kinetic energy matrices, respectively. $Wμν=14c2χμ(σ⋅p̂)V(σ⋅p̂)χν$ gives rise to relativistic corrections and spin couplings in an atomic/molecular system. Here, V is the nuclear attraction, σ is a vector of Pauli spin matrices, and $p̂$ is the momentum operator.

The solutions of the four-component RKB Dirac Hamiltonian⁵⁶ in Eq. (1) are a set of bi-spinor (four-component) molecular orbitals with the large and small components (⁠$ψiL$ and $ψiS$⁠) expanded in a set of atomic spinors centered on different nuclei (R_A),

where $χμA(r−RA)$ is a basis function centered on nucleus A. The orbital coefficients that correspond to positive and negative energy solutions can be written in a matrix form,

In X2C, the transformation matrix takes the form

In Eq. (6), the matrix $Y$ is calculated from the orbital coefficients as

In this work, we employ the one-electron X2C approach, which is a one-step procedure to construct a transformation matrix through the diagonalization of the one-electron four-component core Hamiltonian. Thus, the one-electron spin–orbit coupling is included in the Hamiltonian in the self-consistent field process. However, since the two-electron spin–orbit terms contribute with an opposite sign to the one-electron spin–orbit terms, we used the Boettger factors to scale the one-electron spin–orbit effect to semi-empirically treat the two-electron spin–orbit terms.⁵⁷ 

In the generalized block-localized wave function (X2C-gBLW) approach, we first partition the basis functions into M subspaces according to specific constraints of interest in a given application; the basis functions can be the primitive atomic spinor orbitals or a linear combination thereof. The molecular orbitals are defined as

where C_A is the coefficient matrix for the Ath subspace and T_A is a rectangular transformation matrix that enforces the subspace localization.⁸ In the X2C-gBLW method, bases {χ_μ} are atomic spinors and the resulting C_A is subject to the “picture change” or the X2C transformation from the four-component orbitals. As in the conventional BLW approach, molecular orbitals between different subspace blocks in the X2C-gBLW framework can be nonorthogonal, although the atomic spinors can be fully shared by different blocks in the latter method. Thus, the molecular orbitals may no longer be block-localized in the coefficient matrix in the atomic spinor basis. Importantly, the separation into different subspaces permits us to place different constraints on electronic configurations. For a problem with M number of subspaces, the gBLW approach leads to M eigenvalue problems that can be solved separately.

The molecular spinor coefficients are determined by solving a modified Roothaan–Hall equation for each subspace,

where $F̃A$ and $S̃A$ are the projected Fock and overlap matrices, and $ϵ̃A$ is the diagonal matrix of energy eigenvalues. The projected Fock and overlap matrices are computed using

where F is the Fock matrix in the full basis space, S is the overlap matrix, and P_A is the projection matrix of subspace A.

The projection matrix partitions the full system into subspaces so that the coefficients can be solved separately but also couples the subspaces together both in terms of energy and orthogonality. This results in a wave function that is a variational minimum subject to the constraints of the optimization. The projection matrix is defined as

where D is the full density matrix and D_A is the density matrix for subspace A that has been transformed from the subspace basis into the atomic orbital basis,

Because there is no guarantee that the molecular orbitals will be orthonormal, the density matrices are evaluated using the nonorthogonal formula,

Here, S_A is the overlap matrix of the subspace basis (⁠$SA=TA†STA$⁠). The molecular orbitals for the total system (C) can be assembled by transforming the subspace molecular orbitals (C_A) into the atomic orbital basis and concatenating them together.

The X2C-gBLW method presented in Sec. II B allows for the facile optimization of determinants corresponding to electronically excited configurations, which, in certain applications, may be considered as approximations to localized diabatic excited states. When molecular orbitals of a reference state are used as the subspace transformation matrix, it is possible to isolate a given orbital or a group of orbitals into a subspace that is not mixed with orbitals in other blocks in X2C-gBLW optimization. Consequently, this technique is particularly powerful in resolving degenerate states. For instance, spin–orbit coupling splits the six-fold degenerate ²P term into a two-fold ²P_1/2 and four-fold ²P_3/2 levels. With the X2C-gBLW approach, spin–orbit coupling is included variationally in both the molecular orbitals or the reference state and the X2C-gBLW optimized orbitals. We can exclude spinor molecular orbitals that are not the subject of interest but may cause variational collapse due to mixing in the same symmetry group. Thus, the excited or degenerate state determinant becomes a minimum in the constrained space.⁸ 

Consider the sodium atom frontier excitations (Fig. 1). The electrons in the 1s, 2s, and 2p orbitals comprise the core electrons, while the electron in the 3s orbital is the excited electron. The first excitation roughly corresponds to a spin-flip excitation of the 3s electron to yield an energy degenerate doublet state. Then, 3s → 3p excitations give rise to the ²P term. We define two subspaces for the core electrons and excited electron, respectively. The first subspace contains all ten core electrons and includes all molecular orbitals except the target (the destination of the excited electron) orbital that defines an excited configuration. The second subspace consists of only one electron, the orbital occupied by the excited electron (target orbital), and all unoccupied molecular orbitals that are higher in energy than the target orbital. To prevent the admixture of degenerate orbitals (e.g., the three 3p orbitals of the sodium atom), only the target orbital where the excited electron resides is included in the second subspace, whereas the other degenerate orbitals are excluded. Therefore, the subspace partition is dependent on the specific excited configurations. Since orbitals with lower energy than the target orbital of the excited configuration are excluded from the second subspace, optimization collapse to the ground state configuration is prevented.⁸ A schematic representation of the coefficient matrices for each of the subspaces is given in Fig. 1 for the first few excited determinants in sodium. We used an arbitrary cutoff of 0.01 hartree of the orbital energy to determine degeneracy.

The procedure described above is a form of ΔSCF algorithm that has been augmented with the X2C-gBLW Hamiltonian and automated for convenience of application in the Chronus Quantum software package.^58,59 A series of atomic systems were used to test the algorithm. Calculations on sodium were performed with the Sapporo quadruple zeta basis set with diffuse functions.⁶⁰ Along with Hartree–Fock (wave function theory),⁶¹ the BHandH⁶² and PBE0⁶³ density functionals were used. For group 13 atoms (gallium, indium, and thalium), the Sapporo triple zeta basis set was used along with HF,⁶¹ BHandH,⁶² B3LYP,⁶⁴ and PBE0.⁶³ For the scandium 2+ cation (Sc²⁺), a series of basis sets were used to test the effect that the basis set has on the optimization. The following basis sets were used: 6-31G,⁶⁵ 6-31G(D),⁶⁶ cc-pVDZ, cc-pVTZ, cc-pVQZ,⁶⁷ Sapporo-DZP, Sapporo-TZP, and Sapporo-QZP.⁶⁰ 

For all of the atomic systems, a positively charged state was used to define the transformation matrix in Eq. (8) to yield the molecular orbital basis for partition into the core and excited configuration subspaces. For sodium, gallium, indium, and thalium, the +1 cationic state was used, and for scandium, the 3+ cation was used. For these states, the core electrons are closed shells. Since our goal is to optimize a set of degenerate states, we found that the molecular orbital basis from the closed-shell cation calculations provides a good description, whereas open-shell optimization of the atoms results in orbitals biased toward the ground state.

In this work, we focus on analyzing the quality of X2C-gBLW on computing atomic zero-field splittings and Kramers’ (time-reversal) symmetry recovering, compared to experimental and relativistic MCSCF results.

One of the simplest cases is the sodium atom, which has a small but measurable relativistic effect. The molecular orbital diagram is given in Fig. 1. There is a single electron in the 3s orbital that gives rise to the ²S_1/2 doubly degenerate ground state. The ²P manifold, from the 3s → 3p excitation, is split into ²P_1/2 and ²P_3/2 states. The ²P_1/2 and ²P_3/2 states have two-fold and four-fold degeneracies, respectively. The computed energies of each of the states are given in Table I.

Overall, X2C-gBLW optimization of the block-localized determinant wave functions (denoted by HF) yields the corresponding degenerate states for each of the term symbol groups. However, we did not observe the same using density functional theory (DFT) with block-localized Kohn–Sham orbitals except for the ²S_1/2 states, which are degenerate for all methods. For the ²P_1/2 states, we observe two states that are roughly degenerate with an energy difference of 0.0007 and 0.0003 eV for BHandH and PBE0, respectively. Meanwhile, for the ²P_3/2 states we observed that the two highest energy states are degenerate, but these are not degenerate with the lowest energy states. Thus, the ²P_3/2 manifold is split into two sets of Kramers pairs. The Kramers pairs are close in energy, which suggests that the time-reversal symmetry is roughly maintained. The states within the ²P_3/2 atomic term symbol are different components of the total angular momentum operator (J) with values of −3/2, −1/2, 1/2, and 3/2 ℏ. The 3/2 and −3/2 components along with the 1/2 and −1/2 components are connected by time-reversal symmetry, respectively. Given the absence of a magnetic field, these components are completely degenerate for the exact wave function. In this article, we will refer the symmetry of the total angular momentum operator as J-symmetry. Overall, this suggests that the time-reversal symmetry is roughly maintained, but J-symmetry is not.

A simple method for quantifying the J-symmetry breaking is to compute the energy range within each atomic term symbol. In Fig. 2, we show the range (the difference between highest and lowest energies) within the ²P atomic term symbol in a logarithmic scale. In wave function theory, the range is on the order of 10⁻⁶ eV, whereas it increases to 10⁻²–10⁻³ eV using DFT. Overall, we observe a smaller deviation in energy degeneracy for the ²P_1/2 states than that for the ²P_3/2 states for all methods. Since BHandH and PBE0 include different (50% vs 25%) Hartree–Fock exchange, the observation of roughly similar performance suggests that the amount of Hartree–Fock exchange does not have a significant impact on the results.

In comparison with the experimental excitation energies, there is a smaller mean unsigned error (MUE) for all methods with DFT than that from wave function theory (WFT) calculations. Inclusion of dynamic correlation in X2C-gBLW calculations roughly reduces the MUE by two third. However, the zero-field splitting between ²P_1/2 and ²P_3/2 exhibits a different trend. The energy difference between the ²P_1/2 and ²P_3/2 states is 2.1 meV experimentally. For comparison, our computational results are 2.0, 7.6, and 16.5 meV for WFT, DFT/BHandH, and DFT/PBE0, respectively. The energies for each term symbol are obtained by averaging over energies within each term group. It is surprising that the energy splitting due to spin–orbital coupling is significantly overestimated using DFT.

We computed the ²P atomic transitions in gallium, indium, and thalium. These are all group 13 elements with a single electron in the set of p orbitals. As with the ²P states of sodium, spin–orbit coupling splits the six configurations into two-fold (²P_1/2) and four-fold (²P_3/2) degenerate states. Spin–orbit coupling increases as the atomic number increases, resulting in greater energy differences between the two groups of states. The results are given in Table II.

We observe very good degeneracy among the ²P_1/2 states for all methods, which suggests that time-reversal symmetry is maintained. Interestingly, the time-reversal symmetry is significantly better for these atoms when using DFT with an energy difference on the order of 10⁻⁵ eV. However, similar to sodium, we observe the ²P_3/2 states split roughly into Kramers pairs for all methods. The energy ranges in a logarithmic scale for the ²P_3/2 states using each method and atom are plotted in Fig. 3. Overall, we do not observe a clear trend in breaking of energy degeneracy with respect to the atom number, and the energy range goes from 10⁻¹ to 10⁻³. BHandH and B3lYP appear to give the smallest range for all of the methods, which contrasts with the sodium results, where Hartree–Fock displays the smallest range.

To understand the origin that causes the deviation from energy degeneracy, we optimized the determinants with a frozen core. Here, only the valence p orbital is allowed to relax, which is the orbital that contains the excited electron. We use the energy range of the ²P_3/2 states as a metric for breaking J-symmetry, and the results are given in Table III. We observe no breaking of J-symmetry for the frozen core. Then, we allowed core electrons with specific angular momenta to relax individually from the frozen state. Thus, each set of optimizations begins from the reference state. Here, we observe that allowing the s orbitals to relax causes the largest breaking of J-symmetry, followed by the d orbitals. Allowing the core p orbitals to relax displayed the lowest breaking of symmetry.

Kasper et al. recently showed that the state averaged complete active space self-consistent field (SA-CASSCF) is able to maintain both time-reversal symmetry and J-symmetry to at least 10⁻⁴ eV in the energy range.⁷⁰ In this case, the core orbitals share a single set of orbitals between all states. Since ΔSCF has a different set of core orbitals for each state, it is possible that the algorithm breaks J-symmetry to reach a lower variational minimum in the same way that the unrestricted Hartree–Fock breaks spin symmetry.

In order to compare with experiment, the ²P_3/2 frequencies of Table II are averaged. The results are given in Table IV. Overall, the mean unsigned error compared to experiment ranged from 0.067 to 0.040 for Hartree–Fock and PBE0, respectively. Additionally, the DFT methods show a lower mean unsigned error compared to Hartree–Fock. Hu et al. showed that dynamic correlation reduces the mean unsigned error using multireference configuration interaction singles and doubles (MRCISD).⁵¹ Overall, using DFT reduced the mean unsigned error but not as much as was observed with MRCISD. However, X2C-gBLW with DFT is significantly cheaper than MRCISD.

We tested X2C-gBLW on the scandium 2+ cation. Here, there is a single electron in the d-manifold, and spin–orbit coupling splits the d-manifold into ²D_3/2 and ²D_5/2 states. ²D_3/2 are four-fold degenerate and the ²D_5/2 are six-fold degenerate in the absence of a magnetic field. Overall, this is a more difficult system since the level of degeneracy is higher. Additionally, we computed the excitation to the 4s manifold, resulting in doubly degenerate ²S_1/2 excited states. On this system, we tested the effect due to the size of basis sets on the optimization. The X2C-gBLW energies using Hartree–Fock and the Sapporo bases are given in Table V. The remaining frequencies not displayed in Table V are provided in the supplementary material.

Overall, X2C-gBLW yielded two-fold degeneracy for all ²S_1/2 excited states regardless of the basis set following time-reversal symmetry. Meanwhile, the ²D states roughly optimized into Kramers pairs. Interestingly, for the ²D manifolds, as the size of the basis set increases, the energy difference between the sets of Kramers pairs within each term symbol also increases. As the size of the basis set increases, so does the variational freedom, and the forces that lead to J-symmetry breaking are only encouraged.

We expanded this survey to include other basis set families, such as the Pople and correlation consistent basis sets. In Fig. 4, we present the energy range in a logarithmic scale for each of the ²D term symbols and compute the same using density functional theory. The correlation consistent basis set (using Hartree–Fock) displayed similar performance to the Sapporo basis sets, and we similarly observed the same trends in the energy range within each atomic term symbol as in Table V. For both the correlation consistent and Sapporo basis set families, the ²D_3/2 states usually have smaller ranges than the ²D_5/2 states. However, for the Pople basis set, the trends are less obvious. With regards to DFT, BHandH displayed similar trends to HF in that as the basis set gets larger, so does the energy range, but there is no clear trend using the PBE0 functional. However, the energy range for PBE0 is generally larger than Hartree–Fock or BHandH.

To understand where the J-symmetry breaking originates, we froze the core electrons and only optimized the d electron. Then, we relaxed core electrons of the same angular momenta to show which set of orbitals causes a greater change in the energy range within the same J-symmetry. The results are presented in Table VI. Consistent with findings in group 13 elements, we observed exact degeneracy when the core orbitals are frozen. Additionally, the energy range increased slightly when the s orbitals are relaxed. However, relaxation of the p orbitals produces a significantly larger increase in the energy range.

To further understand this, we optimized the entire core except for the p manifold directly beneath the d orbital (i.e., the 2p manifold). Here, the energy range was only slightly perturbed relative to the frozen-core case. This small change suggests that the largest contribution to the splitting of the Kramers pairs is due to polarization of the p orbital manifold.

The effect of coupling between basis functions of different angular momenta was tested by block-localizing the basis functions with respect to angular momenta. Thus, there is a block for s, p, and d functions, respectively. The results are presented in Table VII. Using Hartree–Fock, both the J-symmetry and time-reversal symmetry are significantly improved. Overall, this suggests that the coupling between basis functions of different angular momenta is the primary cause for the breaking of J-symmetry and, consequently, time-reversal symmetry. However, this is not observed with the density functional theory results, which show little improvement. Additionally, the BHandH functional maintains symmetry significantly better than PBE0, perhaps due to the larger amount of Hartree–Fock exchange. Overall, this suggests that the nonlinear nature of the exchange correlation functional makes it more difficult to obtain correct symmetry even in the best of circumstances.

Furthermore, because these results are constrained, we are able to make a one-to-one comparison between the states using various methods. The DFT energies in Table VII are listed in the order of corresponding states to Hartree–Fock. Particularly, for PBE0, there is a ²D_5/2 Kramers’ pair that is lower in energy than the high energy ²D_3/2 Kramers’ pair. Overall, this suggests that exchange correlation functional approximations with low amounts of Hartree–Fock exchange lack the precision to adequately differentiate the ²D_3/2 and ²D_5/2 states.

Similar to spin symmetry in the non-relativistic case, time-reversal and J-symmetry of the wave function can be constrained from the beginning⁷² in an analogous way to restricted open-shell Hartree–Fock (i.e., Kramers-restricted). While there are merits to both approaches, the Kramer-unrestricted formalism (used in this article) yields the variational minimum and Kramers-restricted methods maintain the symmetry of the wave function.

In Table VIII, we compare the mean unsigned error across basis sets and methods using the unconstrained optimizations. Overall, there is no systematic dependence of the basis set on relative energies. Additionally, the BHandH and PBE0 functionals appeared to perform best when using the correlation consistent and Sapporo bases, while Hartree–Fock yielded roughly equal errors for all basis sets except for cc-pVDZ. In this data set, it was not clear that DFT performed better than Hartree–Fock. Examining the data further, we found that the largest errors in all methods were for the ²S_1/2 states. We found that this frequency was overestimated using DFT for all basis sets (these frequencies are provided in the supplementary material). Meanwhile, for HF, the Pople basis sets overestimated the frequency, while the Sapporo and correlation consistent bases tended to underestimate the frequency. The Sapporo-DZP cc-pVDZ basis sets yielded the lowest error for the ²S_1/2 states, which were also the basis sets where the mean unsigned error was the lowest.

In this article, we introduced an integration of the generalized block-localized wave function (gBLW) method within the relativistic X2C framework. We carried out X2C-gBLW calculations on a number of atomic systems consisting of a single unpaired electron. Using the gBLW algorithm, it is possible to constrain orbital mixing within a group of generalized orbital (molecular orbital) basis functions. By excluding lower energy orbitals in the excited block, we showed that variational collapse to the ground state can be voided. Furthermore, the specific partition scheme of subspace construction allows all of the complementary configurations within a degenerate manifold to be fully optimized. Overall, the computed excitation energies have mean unsigned errors of about 0.1 and 0.04 eV using Hartree–Fock and DFT, respectively. Furthermore, we observed that the optimized orbitals within a J-symmetry group can roughly be grouped into Kramers pairs that yield a range of energies instead of exactly degenerate states. This suggests that the total angular momentum symmetry is not maintained. The self-interaction error^73–76 is a common deficiency to all density functional approximations, and X2C-gBLW does not limit or remove that. Because the self-interaction error is difficult to analyze, it is unclear how the self-interaction error could be affecting the results. However, the energy spread from degeneracy for all the systems that have been examined in this work is an order of magnitude smaller than the energy splitting due to spin–orbit coupling. We found that the largest contribution to the splitting of the Kramers pairs within an atomic term symbol is due to the relaxation of the core electrons in the orbital manifold directly beneath the unpaired electron of interest. Investigations that employ the constrained excited configurations as the basis state in nonorthogonal configuration interaction calculations will be reported in a forthcoming study.

See the supplementary material for the ΔSCF frequencies of Sc²⁺ using Hartree–Fock, BHandH, and PBE0 with a variety of basis sets.

X.L. acknowledges support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, in the Heavy-Element Chemistry program (Grant No. DE-SC0021100), for the development of relativistic electronic structure methods. The development of excited state methods was supported by the Computational Chemical Sciences (CCS) Program of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division in the Center for Scalable and Predictive methods for Excitations and Correlated phenomena (SPEC) at the Pacific Northwest National Laboratory. The development of the open source software package was supported by the U.S. National Science Foundation (Grant Nos. OAC-1663636 and CHE-1856210). Work carried out at the Shenzhen Bay Laboratory was supported by a grant from the Shenzhen Municipal Science and Technology Innovation Commission (Grant No. KQTD2017-0330155106581 to J.G.).

The data that support the findings of this study are available from the corresponding author upon reasonable request.