One-Step Relativistic Driven Similarity Renormalization Group Multireference Perturbation Theory

Abstract: We present an efficient implementation of a one-step relativistic second-order multireference perturbation theory based on the multireference driven similarity renormalization group (MR-DSRG) using the exact two-component (X2C) Hamiltonian, which we denote X2C-DSRG-MRPT2. We show that the X2C-DSRG-MRPT2 method can accurately capture spin--orbit coupling (SOC) effects in the electronic structure of strongly correlated systems containing elements across the periodic table. We further demonstrate that the X2C-DSRG-MRPT2 method, through its variational treatment of SOC effects, can yield spin--orbit splittings with mean absolute percentage errors consistently below 7% with respect to experimental values for systems containing up to sixth row elements. With its modest computational scaling (fifth power in system size) and high accuracy, X2C-DSRG-MRPT2 provides a promising avenue for the routine treatment of relativistic effects in strongly correlated molecular systems.

• The paper introduces X2C-DSRG-MRPT2, which combines an SNSO-X2C Hamiltonian, SOC-inclusive CASSCF orbital optimization, and DSRG-MRPT2 dynamic correlation in one relativistic multireference framework.
• Benchmarks show mean absolute percentage errors below 7% across systems up to the sixth row, including 5.4% for p-block atoms, 4.9% for transition-metal atoms, 3.2% for diatomics, and 3.74% for sixth-row elements.
• The results demonstrate that variational SOC treatment is essential for heavy elements, while density fitting and efficient complex-arithmetic implementation make correlation inexpensive enough for routine calculations on modest hardware.

Overview

This paper by Zhao and Evangelista presents X2C-DSRG-MRPT2, a one-step relativistic second-order multireference perturbation theory that combines the exact two-component (X2C) Hamiltonian with the multireference driven similarity renormalization group (MR-DSRG) formalism. The central claim is that a fully variational treatment of spin–orbit coupling (SOC) at the CASSCF level, followed by DSRG-MRPT2 dynamic correlation, yields spin–orbit splittings with mean absolute percentage errors (MAPE) consistently below 7% for systems containing elements up to the sixth row of the periodic table. The implementation, carried out in complex arithmetic with density fitting and state averaging, is available in the open-source Forte2 package. The work addresses a recognized gap: while two-component multireference methods exist (SHCI, uncontracted X2C-MRCI/MRPT2, X2C-MC-PDFT, X2C-DMRG), one-step approaches in which SOC is treated variationally throughout orbital optimization and correlation remain scarce.

The SNSO-X2C Hamiltonian

The relativistic framework is built on the matrix modified Dirac equation under restricted kinetic balance, block-diagonalized via the X2C decoupling transformation to yield the positive-energy two-component Hamiltonian $h^+$. Because transforming the two-electron integrals exactly would incur four-component cost, the authors retain non-relativistic two-electron integrals and compensate for the resulting two-electron picture-change (2ePC) errors using the screened-nuclear spin–orbit (SNSO) approximation: the spin-dependent part of $h^+$ is scaled by empirical parameters $Q(n,l)$ parameterized against four-component Dirac–Hartree–Fock calculations with the Coulomb–Breit operator. This choice is pragmatic but carries an acknowledged limitation: the regression analysis against experiment indicates an intrinsic error of a few wavenumbers in the SNSO-X2C Hamiltonian itself, so spectroscopic accuracy ($\sim$1 cm$^{-1}$) may not be systematically attainable even with an exact correlation treatment. Atomic mean-field and molecular mean-field alternatives are noted as more rigorous but substantially more complex or expensive.

X2C-CASSCF

The reference wave function is obtained from X2C-CASSCF, which variationally optimizes both CI coefficients and molecular spinors in the presence of the full (spin-free plus spin-dependent) SNSO-X2C Hamiltonian. The implementation uses a density-fitted quasi-second-order two-step algorithm: a Davidson–Liu eigensolver with a complex-arithmetic Harrison–Zarrabian direct-CI $\sigma = Hc$ evaluation for the CI step, and a complex L-BFGS optimizer over skew-Hermitian orbital rotations with gradients and diagonal Hessians adapted to the spinor case. Convergence to $10^{-8}$–$10^{-10}$ $E_h$ is achieved in fewer than 20 macro-iterations for most systems tested. State-averaging over multiple roots uses ensemble-averaged density matrices; Kramers symmetry is not explicitly exploited, though prior work suggests SA-CASSCF recovers Kramers degeneracy essentially fully through orbital optimization.

Approximate schemes for SOC inclusion

A key methodological contribution is a systematic comparison of four approximations that defer activation of the SNSO-X2C Hamiltonian to later stages of the calculation, relative to the full X2C-DSRG-MRPT2 scheme:

Scheme	CASSCF	DSRG-MRPT2
sf-X2C-CASSCF-SO	sf-X2C	— (SO only at final CI iteration)
sf-X2C-DSRG-MRPT2-SO (A)	sf-X2C	sf-X2C (SO at effective-Hamiltonian diagonalization)
sf-X2C-DSRG-MRPT2-SO (B)	sf-X2C	SNSO-X2C
sf-X2C-DSRG-MRPT2-SO (C)	sf-X2C + final SO CI	SNSO-X2C
X2C-DSRG-MRPT2	SNSO-X2C throughout	SNSO-X2C

Benchmarking on second- to sixth-row $p$-block zero-field splittings (ZFS) produces a result that contradicts the natural expectation that accuracy improves monotonically as SOC is introduced earlier. All approximate schemes perform acceptably through the fourth row but degrade sharply for fifth- and sixth-row elements, where they predict qualitatively wrong splittings. Moreover, schemes (B) and (C)—which compute DSRG amplitudes in the presence of the SNSO-X2C Hamiltonian—perform worse than sf-X2C-CASSCF-SO and scheme (A), which do not. The authors attribute this to inconsistent perturbative corrections when amplitudes are formed with a different Hamiltonian than the orbitals were optimized with. The practical implication is that orbital optimization in the presence of SOC is essential for systematic accuracy across the periodic table, and the cheaper spin-free schemes cannot be relied upon for heavy elements.

Benchmark results

p-block atoms. With the decontracted cc-pVTZ basis (matching published four-component data), X2C-DSRG-MRPT2 achieves an MAE of 19.3 cm$h^+$0 and MAPE of 5.4% over 15 splittings, compared to 15.5 cm$h^+$1/4.9% for the much more expensive four-component DSRG-MRPT2 using Dirac–Coulomb–Breit. The error correlation coefficient between the two methods is 0.9497 with an inter-method MAE of only 6.73 cm$h^+$2, validating SNSO-X2C as an efficient surrogate for the four-component Hamiltonian. X2C-DSRG-MRPT2 also matches the near-FCI-quality 4C-iCIPT2 (MAPE 6.5%), indicating that MRPT2 captures the essential dynamic correlation for these systems.

With the near-complete-basis ANO-RCC set extended to 21 splittings including thallium, X2C-DSRG-MRPT2 attains the lowest MAE among all compared methods (56.0 cm$h^+$3, MAPE 6.2%), outperforming even 4C-iCIPT2 (62.6 cm$h^+$4). For sixth-row elements (Tl, Pb, Bi, Po), the MAPE drops further to 3.74%—lower than for lighter elements—and exceeds DKH1-SO-L-PDFT in both MAE and MAPE. Notably, unlike the Breit–Pauli Hamiltonian, whose perturbative expansion is poorly convergent at high $h^+$5, the SNSO-X2C Hamiltonian maintains its accuracy for heavy elements.

Transition metal atoms. For Cu, Ag, Au ($h^+$6) and Sc, Y, La ($h^+$7), X2C-DSRG-MRPT2 gives an MAE of 96.1 cm$h^+$8 and MAPE of 4.9%, comparable to DKH2-NEVPT2 and uncontracted X2C-MRCISD. A methodological advantage emphasized here is that the one-step approach requires only minimal active spaces and state averaging over states of interest, whereas state-interaction methods needed far larger active spaces and dozens to nearly a hundred state-averaged roots to reach comparable accuracy—an obstacle to black-box application.

Diatomics. Across twelve open-shell diatomics up to IO, X2C-DSRG-MRPT2 achieves an MAE of 31.7 cm$h^+$9 (MAPE 3.2%), the best among the multireference methods compared and approaching single-reference EOM-CCSD-based results (MAE 17.7–23.0 cm$Q(n,l)$0).

TlH potential energy curves. Potential energy curves for six low-lying states of TlH agree closely with ic-MRCI+Q-SO reference curves across 1.4–3.8 Å, correctly reproducing the weakly bound $Q(n,l)$1 state, a shoulder in the $Q(n,l)$2 curve, and a crossing near 2.0 Å. Spectroscopic constants for the $Q(n,l)$3 ground state ($Q(n,l)$4 = 1.8706 Å, $Q(n,l)$5 = 1379.1 cm$Q(n,l)$6, $Q(n,l)$7 = 2.04 eV) are comparable to uc-X2C-MRPT2 and superior to 4C-CASPT2, against experimental values of 1.872 Å, 1390.7 cm$Q(n,l)$8, and 2.06 eV.

Efficiency

The method scales as $Q(n,l)$9 in core/virtual spinors (with an $\sim$0 active-space term), with a prefactor roughly 64 times larger than the non-relativistic counterpart due to two-component spinors and complex arithmetic. Density fitting avoids storage of four-index intermediates for the dominant contractions. A representative TlH computation completes in about 630 seconds on a laptop (Apple M2, 16 GB RAM) using ~1.2 GB of memory, with the DSRG-MRPT2 step costing only ~30 s versus ~533 s for X2C-CASSCF. Dynamic correlation therefore adds negligible overhead relative to the reference calculation—a strong argument for routine applicability.

Role of dynamical correlation

An appendix-level analysis comparing X2C-CASSCF and X2C-DSRG-MRPT2 errors across all systems shows that both methods give approximately unbiased error distributions, but the correlated results are markedly more tightly clustered around zero error, with the largest improvements occurring for heavy-element systems such as IO, Tl, and Po. Orbital relaxation alone (as in sf-X2C-CASSCF-SO) does not suffice; consistent dynamical correlation on top of SOC-inclusive orbital optimization is required. Conversely, adding dynamical correlation inconsistently—as in scheme (C)—can degrade accuracy relative to no correlation treatment at all.

Limitations and open questions

Several caveats bear directly on the reported accuracies. First, the SNSO screening parameters are empirical, and the intrinsic few-wavenumber error of the SNSO-X2C Hamiltonian limits the path to spectroscopic accuracy regardless of the correlation treatment; atomic or molecular mean-field treatments of 2ePC effects remain untested within this framework. Second, the comparison to four-component results relies on matching basis sets and frozen-core choices dictated by the four-component calculations, and the cc-pVTZ comparisons are explicitly acknowledged as not approaching the complete-basis limit. Third, the implementation does not exploit Kramers symmetry, doubling the effective cost relative to a quaternion-symmetric formulation. Finally, extension to third order (X2C-DSRG-MRPT3), which showed consistent improvements in the four-component case, remains an open question addressed in ongoing work.

Conclusion

This paper establishes X2C-DSRG-MRPT2 as an accurate and computationally modest one-step approach to SOC in strongly correlated systems. Its principal findings are that variational SOC treatment at the CASSCF level is necessary for reliable splittings beyond the fourth row, that the SNSO-X2C Hamiltonian reproduces four-component Dirac–Coulomb–Breit spectra with high fidelity, and that DSRG-MRPT2 dynamic correlation can be added at near-zero additional cost. With MAPEs below 7% across the periodic table—including 3.74% for sixth-row atoms—the method is competitive with the most accurate multireference alternatives while remaining tractable on commodity hardware.