Calculation of electric-field gradients based on higher-order generalized Douglas-Kroll transformations

In this paper, the calculation of electric-field-like properties based on higher-order Douglas-Kroll-Hess (DKH) transformations is discussed. The electric-field gradient calculated within the Hartree-Fock self-consistent field framework is used as a representative property. The properties are expressed as an analytic first derivative of the four-component Dirac energy and the nth-order DKH energy, respectively. The differences between a "forward" transformation of the relativistic energy or the "back transformation" of the wave function is discussed in some detail. Detailed test calculations were carried out on the electric-field gradient at the halogen nucleus in the series HX (X=F,Cl,Br,I,At) for which extensive reference data are available. The DKH method is shown to reproduce (spin-free) four-component Dirac-Fock results to an accuracy of better than 99% which is significantly closer than previous DKH studies. The calculations of both the Hamiltonian and the property operator are shown to be essentially converged after the second-order transformation, even for elements as heavy as At. In addition, we have obtained results within the density-functional framework using the DKHZ and zeroth-order regular approximation (ZORA) methods. The latter results included picture-change effects at the scalar relativistic variant of the ZORA-4 level and were shown to be in quantitative agreement with earlier results obtained by van Lenthe and Baerends. The picture-change effects are somewhat smaller for the ZORA method compared to DKH. For heavier elements significant differences in the field gradients predicted by the two methods were found. Based on comparison with four-component Dirac-Kohn-Sham calculations, the DKH results are more accurate. Compared to the spin-free Dirac-Kohn-Sham reference values, the ZORA-4 formalism did not improve the results of the ZORA calculations.     

Exact decoupling of the Dirac Hamiltonian. II. The generalized Douglas-Kroll-Hess transformation up to arbitrary order

In order to achieve exact decoupling of the Dirac Hamiltonian within a unitary transformation scheme, we have discussed in part I of this series that either a purely numerical iterative technique (the Barysz-Sadlej-Snijders method) or a stepwise analytic approach (the Douglas-Kroll-Hess method) are possible. For the evaluation of Douglas-Kroll-Hess Hamiltonians up to a pre-defined order it was shown that a symbolic scheme has to be employed. In this work, an algorithm for this analytic derivation of Douglas-Kroll-Hess Hamiltonians up to any arbitrary order in the external potential is presented. We discuss how an estimate for the necessary order for exact decoupling (within machine precision) for a given system can be determined from the convergence behavior of the Douglas-Kroll-Hess expansion prior to a quantum chemical calculation. Once this maximum order has been accomplished, the spectrum of the positive-energy part of the decoupled Hamiltonian, e.g., for electronic bound states, cannot be distinguished from the corresponding part of the spectrum of the Dirac operator. An efficient scalar-relativistic implementation of the symbolic operations for the evaluation of the positive-energy part of the block-diagonal Hamiltonian is presented, and its accuracy is tested for ground-state energies of one-electron ions over the whole periodic table. Furthermore, the first many-electron calculations employing sixth up to fourteenth order DKH Hamiltonians are presented.     

One-electron contributions to the g-tensor for second-order Douglas-Kroll-Hess theory

The electric g-tensor is a central quantity for the interpretation of electron paramagnetic resonance spectra. In this paper, a detailed derivation of the 1-electron contributions to the g-tensor is presented in the framework of linear response theory and the second-order Douglas-Kroll-Hess (DKH) transformation. Importantly, the DKH transformation in the presence of a magnetic field is not unique. Whether or not the magnetic field is included in the required Foldy-Wouthuysen transformation, different transformation matrices and, consequently, Hamiltonians result. In this paper, a detailed comparison of both approaches is presented, paying particular attention to the mathematical properties of the resulting Hamiltonians. In contrast to previous studies that address the g-tensor in the framework of DKH theory, the resulting terms are compared to those of the conventional Pauli theory and are given a physical interpretation. Based on these mathematical and physical arguments, we establish that the proper DKH transformation for systems with constant magnetic fields is based on a gauge-invariant Foldy-Wouthuysen transformation, i.e., a Foldy-Wouthuysen transformation including the magnetic field. Calculations using density functional theory (DFT) are carried out on a set of heavy, diatomic molecules, and a set of transition-metal complexes. Based on these calculations, the performance of the relativistic calculation with and without inclusion of picture-change effects is compared. Additionally, the g-tensor is calculated for the Lanthanide dihydrides. Together with the results from the other two molecular test sets, these calculations serve to quantify the magnitude of picture-change effects and elucidate trends across the periodic table.     

Density functional based structure optimization for molecules containing heavy elements: analytical energy gradients for the Douglas-Kroll-Hess scalar relativistic approach to the LCGTO-DF method

The self-consistent scalar-relativistic linear combination of Gaussian-type orbitals density functional (LCGTO-DF) method has been extended to calculate analytical energy gradients. The method is based on the use of a unitary second order Douglas-Kroll-Hess (DKH) transformation for decoupling large and small components of the full four-component Dirac-Kohn-Sham equation. The approximate DKH transformation most common in molecular calculations has been implemented; this variant employs nuclear potential based projectors and it leaves the electron-electron interaction untransformed. Examples are provided for the geometry optimization of a series of heavy metal systems which feature a variety of metal-ligand bonds, like Au₂, AuCl, AuH, Mo(CO)₆ and W(CO)₆ as well as the d¹⁰ complexes [Pd(PH₃)₂O₂] and [Pt(PH₃)₂O₂]. The calculated results, obtained with several gradient-corrected exchange-correlation potentials, compare very well with experimental data and they are of similar or even better accuracy than those of other high quality relativistic calculations reported so far.     

Exact decoupling of the Dirac Hamiltonian. I. General theory

Exact decoupling of positive- and negative-energy states in relativistic quantum chemistry is discussed in the framework of unitary transformation techniques. The obscure situation that each scheme of decoupling transformations relies on different, but very special parametrizations of the employed unitary matrices is critically analyzed. By applying the most general power series ansatz for the parametrization of the unitary matrices it is shown that all transformation protocols for decoupling the Dirac Hamiltonian have necessarily to start with an initial free-particle Foldy-Wouthuysen step. The purely numerical iteration scheme applying X-operator techniques to the Barysz-Sadlej-Snijders (BSS) Hamiltonian is compared to the analytical schemes of the Foldy-Wouthuysen (FW) and Douglas-Kroll-Hess (DKH) approaches. Relying on an illegal 1/c expansion of the Dirac Hamiltonian around the nonrelativistic limit, any higher-order FW transformation is in principle ill defined and doomed to fail, irrespective of the specific features of the external potential. It is shown that the DKH method is the only valid analytic unitary transformation scheme for the Dirac Hamiltonian. Its exact infinite-order version can be realized purely numerically by the BSS scheme, which is only able to yield matrix representations of the decoupled Hamiltonian but no analytic expressions for this operator. It is explained why a straightforward numerical iterative extension of the DKH procedure to arbitrary order employing matrix representations is not feasible within standard one-component electronic structure programs. A more sophisticated ansatz based on a symbolical evaluation of the DKH operators via a suitable parser routine is needed instead and introduced in Part II of this work.     

Exact decoupling of the Dirac Hamiltonian. IV. Automated evaluation of molecular properties within the Douglas-Kroll-Hess theory up to arbitrary order

In Part III [J. Chem. Phys. 124, 064102 (2005)] of this series of papers on exact decoupling of the Dirac Hamiltonian within transformation theory, we developed the most general account on how to treat magnetic and electric properties in a unitary transformation theory on the same footing. In this paper we present an implementation of a general algorithm for the calculation of magnetic as well as electric properties within the framework of Douglas-Kroll-Hess theory. The formal and practical principles of this algorithm are described. We present the first high-order Douglas-Kroll-Hess results for property operators. As for model properties we propose to use the well-defined radial moments, i.e., expectation values of r(k), which can be understood as terms of the Taylor-series expansion of any property operator. Such moments facilitate a rigorous comparison of methods free of uncertainties which may arise in a direct comparison with experiment. This is important in view of the fact that various approaches to two-component molecular properties may yield numerically very small terms whose approximate or inaccurate treatment would not be visible in a direct comparison to experimental data or to another approximate computational reference. Results are presented for various degrees of decoupling of the model properties within the Douglas-Kroll-Hess scheme.     

Douglas–Kroll–Hess Theory: a relativistic electrons-only theory for chemistry

A unitary transformation allows to separate (block-diagonalize) the Dirac Hamiltonian into two parts one part: solely describes electrons, while the other gives rise to negative-energy states, which are the so-called positronic states. The block-diagonal form of the Hamiltonian no longer accounts for the coupling of both kinds of states. The positive-energy (‘electrons-only’) part can serve as a ‘fully’ relativistic electrons-only theory, which can be understood as a rigorous basis for chemistry. Recent developments of the Douglas–Kroll–Hess (DKH) method allowed to derive a sequence of expressions, which approximate this electrons-only Hamiltonian up to arbitrary-order. While all previous work focused on the numerical stability and accuracy of these arbitrary-order DKH Hamiltonians, conceptual issues and paradoxa of the method were mostly left aside. In this work, the conceptual side of DKH theory is revisited in order to identify essential aspects of the theory to be distinguished from purely computational consideration.     

Local relativistic exact decoupling

We present a systematic hierarchy of approximations for local exact decoupling of four-component quantum chemical Hamiltonians based on the Dirac equation. Our ansatz reaches beyond the trivial local approximation that is based on a unitary transformation of only the atomic block-diagonal part of the Hamiltonian. Systematically, off-diagonal Hamiltonian matrix blocks can be subjected to a unitary transformation to yield relativistically corrected matrix elements. The full hierarchy is investigated with respect to the accuracy reached for the electronic energy and for selected molecular properties on a balanced test molecule set that comprises molecules with heavy elements in different bonding situations. Our atomic (local) assembly of the unitary exact-decoupling transformation--called local approximation to the unitary decoupling transformation (DLU)--provides an excellent local approximation for any relativistic exact-decoupling approach. Its order-N(2) scaling can be further reduced to linear scaling by employing a neighboring-atomic-blocks approximation. Therefore, DLU is an efficient relativistic method well suited for relativistic calculations on large molecules. If a large molecule contains many light atoms (typically hydrogen atoms), the computational costs can be further reduced by employing a well-defined nonrelativistic approximation for these light atoms without significant loss of accuracy. We also demonstrate that the standard and straightforward transformation of only the atomic block-diagonal entries in the Hamiltonian--denoted diagonal local approximation to the Hamiltonian (DLH) in this paper--introduces an error that is on the order of the error of second-order Douglas-Kroll-Hess (i.e., DKH2) when compared with exact-decoupling results. Hence, the local DLH approximation would be pointless in an exact-decoupling framework, but can be efficiently employed in combination with the fast to evaluate DKH2 Hamiltonian in order to speed up calculations for which ultimate accuracy is not the major concern.     

Analytic energy gradients for the spin-free exact two-component theory using an exact block diagonalization for the one-electron Dirac Hamiltonian

We report the implementation of analytic energy gradients for the evaluation of first-order electrical properties and nuclear forces within the framework of the spin-free (SF) exact two-component (X2c) theory. In the scheme presented here, referred to in the following as SFX2c-1e, the decoupling of electronic and positronic solutions is performed for the one-electron Dirac Hamiltonian in its matrix representation using a single unitary transformation. The resulting two-component one-electron matrix Hamiltonian is combined with untransformed two-electron interactions for subsequent self-consistent-field and electron-correlated calculations. The "picture-change" effect in the calculation of properties is taken into account by considering the full derivative of the two-component Hamiltonian matrix with respect to the external perturbation. The applicability of the analytic-gradient scheme presented here is demonstrated in benchmark calculations. SFX2c-1e results for the dipole moments and electric-field gradients of the hydrogen halides are compared with those obtained from nonrelativistic, SF high-order Douglas-Kroll-Hess, and SF Dirac-Coulomb calculations. It is shown that the use of untransformed two-electron interactions introduces rather small errors for these properties. As a first application of the analytic geometrical gradient, we report the equilibrium geometry of methylcopper (CuCH(3)) determined at various levels of theory.     

The employment of relativistic adapted Gaussian basis sets in Douglas–Kroll–Hess scalar calculations with diatomic molecules

The prolapse-free relativistic adapted Gaussian basis sets (RAGBSs), developed by our research group on the basis of the four-component approach, are used for the first time in Douglas–Kroll–Hess 2nd order scalar relativistic calculations (DKH2) of simple diatomic molecules containing Hydrogen and the halogens from Fluorine up to Iodine: HX and X₂, where X = F, Cl, Br, and I. To this end, the RAGBSs were contracted with the general contraction scheme to triple-, quadruple-, and quintuple-zeta sets. Polarization functions were also added to the basis sets by optimization with the configuration interaction method including single and double excitations into the DKH2 environment, DKH2-CISD. The molecular properties were then calculated with the coupled cluster electronic correlation treatment and the DKH2 scalar relativistic method, DKH2-CCSD(T), and indicated that our RAGBSs should be contracted as quadruple-zeta basis sets. The results achieved with the DKH2-CCSD(T) calculations and the selected quadruple-zeta RAGBSs are able to reproduce the experimental data of equilibrium distances, dissociation energies, and harmonic vibrational frequencies with root-mean-square (rms) errors of 0.015 Å, 3.6 kcal mol⁻¹, and 21.7 cm⁻¹, respectively.     

Picture change error correction in the radial distributions of canonical orbital densities and total electron density of radon atom: the effect of the size of nucleus and the basis set limit

The 2nd order Douglas-Kroll-Hess (DKH2) and the Infinite Order Two Component (IOTC) radial distributions of electron density of canonical Hartree-Fock (HF) orbitals of radon atom are presented. Furthermore, the total electron density is revisited. The picture change error (PCE) correction is investigated by analytical means. The point charge model of nucleus and the Gaussian nucleus model are employed. The basis set is extrapolated by means of including tight s and also p Gaussians within the original triple zeta basis set. It is found that the DKH1 PCE corrected DKH2 total electron and s orbital contact densities are negative for the point charge model of nucleus if tight enough s Gaussians are included in the basis set. It is shown that this failure is caused due to the missing terms of the second order Douglas-Kroll transformation for the DKH2 electron density. PCE is found the most striking in the DKH2/IOTC electron density of s orbitals close to the nucleus. The radial distributions of the 2-component p _1/2 orbital densities are considerably affected by PCE at the nucleus as well. Furthermore, the PCE corrected DKH2/IOTC scalar p orbital densities have a non-zero value of electron density at nucleus and can be considered as an spin-orbit (SO) average of the p _1/2 and p _3/2 orbitals. The d and f orbitals are affected by PCE in the vicinity of the nucleus only little. The PCE corrected DKH2 and IOTC radial distributions of orbital densities are nodeless, which is completely in agreement with the radial distribution of the analytic or numeric DCH orbital densities.     

A fully relativistic method for calculation of nuclear magnetic shielding tensors with a restricted magnetically balanced basis in the framework of the matrix Dirac-Kohn-Sham equation

A new relativistic four-component density functional approach for calculations of NMR shielding tensors has been developed and implemented. It is founded on the matrix formulation of the Dirac-Kohn-Sham (DKS) method. Initially, unperturbed equations are solved with the use of a restricted kinetically balanced basis set for the small component. The second-order coupled perturbed DKS method is then based on the use of restricted magnetically balanced basis sets for the small component. Benchmark relativistic calculations have been carried out for the (1)H and heavy-atom nuclear shielding tensors of the HX series (X=F,Cl,Br,I), where spin-orbit effects are known to be very pronounced. The restricted magnetically balanced basis set allows us to avoid additional approximations and/or strong basis set dependence which arises in some related approaches. The method provides an attractive alternative to existing approximate two-component methods with transformed Hamiltonians for relativistic calculations of chemical shifts and spin-spin coupling constants of heavy-atom systems. In particular, no picture-change effects arise in property calculations.     

Orbitally invariant internally contracted multireference unitary coupled cluster theory and its perturbative approximation: theory and test calculations of second order approximation

A unitary wave operator, exp (G), G(+) = -G, is considered to transform a multiconfigurational reference wave function Φ to the potentially exact, within basis set limit, wave function Ψ = exp (G)Φ. To obtain a useful approximation, the Hausdorff expansion of the similarity transformed effective Hamiltonian, exp (-G)Hexp (G), is truncated at second order and the excitation manifold is limited; an additional separate perturbation approximation can also be made. In the perturbation approximation, which we refer to as multireference unitary second-order perturbation theory (MRUPT2), the Hamiltonian operator in the highest order commutator is approximated by a Mo?ller-Plesset-type one-body zero-order Hamiltonian. If a complete active space self-consistent field wave function is used as reference, then the energy is invariant under orbital rotations within the inactive, active, and virtual orbital subspaces for both the second-order unitary coupled cluster method and its perturbative approximation. Furthermore, the redundancies of the excitation operators are addressed in a novel way, which is potentially more efficient compared to the usual full diagonalization of the metric of the excited configurations. Despite the loss of rigorous size-extensivity possibly due to the use of a variational approach rather than a projective one in the solution of the amplitudes, test calculations show that the size-extensivity errors are very small. Compared to other internally contracted multireference perturbation theories, MRUPT2 only needs reduced density matrices up to three-body even with a non-complete active space reference wave function when two-body excitations within the active orbital subspace are involved in the wave operator, exp (G). Both the coupled cluster and perturbation theory variants are amenable to large, incomplete model spaces. Applications to some widely studied model systems that can be problematic because of geometry dependent quasidegeneracy, H4, P4, and BeH(2), are performed in order to test the new methods on problems where full configuration interaction results are available.     

Resolution of identity Dirac-Kohn-Sham method using the large component only: Calculations of g-tensor and hyperfine tensor

A new relativistic two-component density functional approach, based on the Dirac-Kohn-Sham method and an extensive use of the technique of resolution of identity (RI), has been developed and is termed the DKS2-RI method. It has been applied to relativistic calculations of g and hyperfine tensors of coinage-metal atoms and some mercury complexes. The DKS2-RI method solves the Dirac-Kohn-Sham equations in a two-component framework using explicitly a basis for the large component only, but it retains all contributions coming from the small component. The DKS2-RI results converge to those of the four-component Dirac-Kohn-Sham with an increasing basis set since the error associated with the use of RI will approach zero. The RI approximation provides a basis for a very efficient implementation by avoiding problems associated with complicated integrals otherwise arising from the elimination of the small component. The approach has been implemented in an unrestricted noncollinear two-component density functional framework. DKS2-RI is related to Dyall's [J. Chem. Phys. 106, 9618 (1997)] unnormalized elimination of the small component method (which was formulated at the Hartree-Fock level and applied to one-electron systems only), but it takes advantage of the local Kohn-Sham exchange-correlation operators (as, e.g., arising from local or gradient-corrected functionals). The DKS2-RI method provides an attractive alternative to existing approximate two-component methods with transformed Hamiltonians (such as Douglas-Kroll-Hess [Ann. Phys. 82, 89 (1974); Phys. Rev. A 33, 3742 (1986)] method, zero-order regular approximation, or related approaches) for relativistic calculations of the structure and properties of heavy-atom systems. In particular, no picture-change effects arise in the property calculations.     

A general nonlinear expansion form for electronic wave functions

A new expansion form is presented for electronic wave functions. The wave function is a linear combination of product basis functions, and each product basis function in turn is formally equivalent to a linear combination of configuration state functions that comprise an underlying linear expansion space. The expansion coefficients that define the basis functions are nonlinear functions of a smaller number of variables. The expansion form is appropriate for both ground and excited states and to both closed and open shell molecules. The method is formulated in terms of spin-eigenfunctions using the graphical unitary group approach (GUGA), and consequently it does not suffer from spin contamination.     

Benchmarking the multipole shielding polarizability/reaction field approach to solvation against QM/MM: applications to the shielding constants of N-methylacetamide

We present a benchmark study of a combined multipole shielding polarizability/reaction field (MSP/RF) approach to the calculation of both specific and bulk solvation effects on nuclear magnetic shielding constants of solvated molecules. The MSP/RF scheme is defined by an expansion of the shielding constants of the solvated molecule in terms of electric field and field gradient property derivatives derived from single molecule ab initio calculations. The solvent electric field and electric field gradient are calculated based on data derived from molecular dynamics simulations, thereby accounting for solute-solvent dynamical effects. The MSP/RF method is benchmarked against polarizable quantum mechanics/molecular mechanics (QM/MM) calculations. The best agreement between the MSP/RF and QM/MM approaches is found by truncating the electric field expansion in the MSP/RF approach at the linear electric field level which is due to the cancelation of errors. In addition, we investigate the sensitivity of the results due to the choice of one-electron basis set in the ab initio calculations of the property derivatives and find that these derivatives are affected by the basis set in a way similar to the shielding constants themselves.     

Calculations of electronic structure of the UF<Subscript>6</Subscript> molecule and the UO<Subscript>2</Subscript> crystal with a relativistic pseudopotential

The influence of relativistic effects on the properties of uranium hexafluoride was considered. Detailed comparison of the spectrum of one-electron energies obtained in the nonrelativistic (by the Hartree-Fock method), relativistic (by the Dirac-Fock method), and scalar-relativistic (using a relativistic potential of the uranium atom core) calculations was carried out. The methods of optimization of atomic basis in the LCAO calculations of molecules and crystals are discussed which make it possible to consider distortion of atomic orbitals upon the formation chemical bonds. The influence of the atomic basis optimization on the results of scalar-relativistic calculations of the molecule UF₆ properties is analyzed. Calculations of the electronic structure and properties of UO₂ crystals with relativistic and nonrelativistic pseudopotentials are fulfilled.