Scalar Relativistic Hamiltonians

Zeroth-order regular approximation (ZORA)

Code author: Nathan Gillispie and Daniel R. Nascimento

Section author: Nathan Gillispie and Daniel R. Nascimento

Source: zora.cc zora.h

ZORA is a perturbative approximation of the full relativistic Hamiltonian for DFT and wavefunction-based methods. It is commonly used to provide a more accurate total energy in DFT calculations.

When ZORA is used in PSI4, it creates a scalar relativistic kinetic energy integral to be used in the SCF procedure. It has been tested for HF and DFT references with the driver methods energy and optimize. Our implementation employs a grid-based scheme, therefore analytic energy gradients are not supported.

Equations and implementation details are based on a paper by Pak, Dada and Nascimento [Pak:2025:094110].

Usage

To use the ZORA Hamiltonian with default grid settings, use option RELATIVISTIC with ZORA. To change the number of grid points, use the ZORA_RADIAL_POINTS and ZORA_SPHERICAL_POINTS options.

set {
    reference rhf
    scf_type pk
    relativistic zora
    zora_radial_points 160
    zora_spherical_points 1202
}

Note

The number of spherical points must be a Lebedev number. See Grid Selection for a list of all options.

It may be useful to compute the non-relativistic kinetic integral using the grid points to compare against the analytic kinetic integrals. To do this, use the ZORA_NR_DEBUG option.

set {
    relativistic zora
    zora_nr_debug true
}

See Keywords for more options.

Theory

In short, the FW-transformed Dirac Hamiltonian is perturbatively expanded with respect to an expression ($E/(2m{c}^2-V)$) involving the potential $V$. To the zeroth-order, we get the ZORA Hamiltonian. Below is the ZORA Kohn-Sham (ZKS) Hamiltonian, however, results are analagous for other methods.

(1)$${\hat{h}}^{\text{ZKS}}=\frac{{\mathbf{p}}^2}{2}+\mathbf{p}\left(\frac{\kappa -1}{2}\right)\mathbf{p}+\frac{{\kappa }^2}{4c^2}\sigma \cdot (\mathrm{\nabla }{v}^{\text{KS}}\times \mathbf{p})+v^{\text{KS}}$$

$\mathbf{p}$ is the momentum operator, $\sigma$ are the Pauli matrices, $\kappa =[1-v^{\text{KS}}/(2c^2){]}^{-1}$, and $v^{\text{KS}}$ is the usual Kohn-Sham potential. In this form, the terms of the ZORA Hamiltonian are separated into the classical kinetic energy, scalar relativistic correction to the kinetic energy, spin-orbit, and KS potential terms, respectively. Were this implemented as is, the relativistic terms would depend on the KS potential, which is impracticable due to gauge-invariance and convergence issues.

The procedure given by Van Wüllen [vanWullen:1998:392] avoids these problems by replacing $v^{\text{KS}}$ in the relativistic terms of (1) with an effective potential $v_{\text{eff}}(\mathbf{r})$. This is the potential experienced at a point due to nuclear attraction and electron repulsion from a model potential that reproduces the nuclear cusp on each atom. The model potential comes from a model basis of Gaussian-type orbitals, whose values were obtained from NWChem [NWChem:2020]. This means that geometry is the only thing affecting $v_{\text{eff}}(\mathbf{r})$.

The ZORA scalar relativistic kinetic integral $T^{\text{SR}}$ is the first two terms of (1) in the atomic orbital basis. Once $v_{\text{eff}}(\mathbf{r})$ is computed on a grid, $T^{\text{SR}}$ can be calculated as

$$T_{\mu \nu }^{\text{SR}}=\int d{\mathbf{r}}^3\frac{c^2}{2c^2-v_{\text{eff}}(\mathbf{r})}\mathrm{\nabla }{\chi }_{\mu }^{†}(\mathbf{r})\cdot \mathrm{\nabla }{\chi }_{\nu }(\mathbf{r}),$$

given atomic orbitals $\chi (\mathrm{r})$. Notice that as $v_{\text{eff}}\to 0$, $T^{\text{SR}}$ simply becomes the non-relativistic kinetic energy. This is the basis behind the ZORA_NR_DEBUG option. Because $T^{\text{SR}}$ is only evaluated once before the SCF procedure, don’t cheap out on the grid! In practice, the time spent computing the ZORA Hamiltonian is relatively negligible.

Limitations

• Spin-orbit coupling effects are not available because they require a complex generalized SCF procedure.

• ZORA theory allows for adjustment of the molecular orbital energies. This provides an important correction for linear response calculations. Currently, this is not implemented in PSI4.

• Analytic energy gradients are not available.

Keywords

RELATIVISTIC

Relativistic Hamiltonian type

• Type: string

• Possible Values: NO, X2C, ZORA

• Default: NO

ZORA_RADIAL_POINTS

Number of radial points for the ZORA effective potential grid. The ZORA calculation is relatively fast, and only happens once, so don’t cheap out on the radial points!

• Type: integer

• Default: 140

ZORA_SPHERICAL_POINTS

Number of spherical points for the ZORA effective potential grid

• Type: integer

• Default: 2030

ZORA_PRUNING_SCHEME

Pruning scheme for the ZORA effective potential grid. P_slater is the best option if you must prune, but none is recommended. Robust and Treutler are not recommended for the ZORA grid as they cut too many points near the nuclear cusp.

• Type: string

• Possible Values: NONE, P_SLATER, ROBUST, LOG_SLATER, TREUTLER

• Default: NONE

ZORA_BASIS_TOLERANCE

Basis tolerance for the ZORA effective potential grid

• Type: conv double

• Default: 1e-12

ZORA_NR_DEBUG

Compute the non-relativistic kinetic energy with the ZORA code. Useful when comparing analytic and grid-based methods.

• Type: boolean

• Default: false

Exact two-component (X2C)

Code author: Prakash Verma and Francesco A. Evangelista

Section author: Prakash Verma, Wallace D. Derricotte, and Francesco A. Evangelista

The X2C approach is a convenient way to introduce scalar relativistic effects in DFT and wave function-based methods. PSI4 implements the spin-free one-electron version of X2C, which produces a modified one-electron Hamiltonian $H_{\mathrm{X}2\mathrm{C}}$:

$$H_{\mathrm{X}2\mathrm{C}}=T_{\mathrm{X}2\mathrm{C}}+V_{\mathrm{X}2\mathrm{C}}$$

that is a sum of a kinetic energy ($T_{\mathrm{X}2\mathrm{C}}$) and potential energy ($V_{\mathrm{X}2\mathrm{C}}$) operator. Our implementation is equivalent to the one reported by Cheng and Gauss [Cheng:084114]. X2C calculations require the use of special (alternatively fully uncontracted) basis sets designed for relativistic calculations. Common choices include the Dunning Douglass–Kroll basis sets (cc-pVXZ-DK, cc-pCVXZ-DK, cc-pwCVXZ-DK) and Roos’ ANO basis sets.

Note

See also Interface to DKH by A. Wolf, M. Reiher, and B. A. Hess for another relativistic Hamiltonian.

A First Example

The following is a simple input that will perform a Hartree–Fock calculation using the X2C Hamiltonian.

molecule {
  H
  F 1 0.92
}

set {
    scf_type pk
    basis cc-pvdz
    relativistic x2c
}

energy('hf')

This computation yields the following result:

@RHF Final Energy:  -100.10007984692388

 => Energetics <=

  Nuclear Repulsion Energy =              5.1767335622934780
  One-Electron Energy =                -150.7611816259664579
  Two-Electron Energy =                  45.4843682167491039
  Total Energy =                       -100.1000798469238902

while a non-relativistic calculation yields the following energy:

@RHF Final Energy:  -100.01928891411315

 => Energetics <=

  Nuclear Repulsion Energy =              5.1767335622934780
  One-Electron Energy =                -150.6645256529074572
  Two-Electron Energy =                  45.4685031765008461
  Total Energy =                       -100.0192889141131474

Basis sets options

The X2C module in PSI4 supports different combinations of basis set. By default, if the input file specifies only BASIS, then the X2C module will solve the modified Dirac equation in an uncontracted basis and then recontract the X2C Hamiltonian in the original basis. Alternatively, the user can use BASIS_RELATIVISTIC to specify a different basis set to solve the modified Dirac equation.

set {
    basis cc-pvdz-dk
    basis_relativistic cc-pvtz-dk
    relativistic x2c
}

It is recommended that when employing the X2C relativistic Hamiltonian, that you use a fully decontracted basis set. This can be done simply in the input by adding “-decon” to the name of the primary basis you want to use for the calculation as detailed in Decontracted Basis Sets. Publications resulting from the use of X2C should cite the following publication: [Verma:2015]

Theory

X2C is based on exact decoupling of positive-energy ( $h_{++}^{FW}$ ) and negative-energy ($h_{--}^{FW}$ ) blocks of the Dirac Hamiltonian ($h^D$).

$$\begin{array}{r}{U}^{†}{h}^{\mathrm{D}}U=U^{†}\left(\begin{array}{cc}{h}_{LL}& h_{LS}\\ h_{SL}& h_{SS}\end{array}\right)U=\left(\begin{array}{cc}{h}_{++}^{\mathrm{F}\mathrm{W}}& 0\\ 0& h_{--}^{\mathrm{F}\mathrm{W}}\end{array}\right)\end{array}$$

The transformation ( $U$ ) is obtained from the solutions of the Dirac equation in kinetically balanced basis [Kutzelnigg:1984] treatment. In the X2C treatment, the positive-energy block of the Hamiltonian ( $h_{++}^{FW}$ ) is given by the sum of a transformed kinetic ($T_{\mathrm{X}2\mathrm{C}}$) and potential energy ( $V_{\mathrm{X}2\mathrm{C}}$ ) contribution. Relativistic kinetic energy ( $T_{\mathrm{X}2\mathrm{C}}$ ) and nuclear-electron interaction potential ( $V_{\mathrm{X}2\mathrm{C}}$ ) is given in terms of non-relativisitc kinetic ($T={\hat{p}}^2/2$) energy and nuclear-electron interaction potential ($V$), coupling matrix ( $X$) and renormalization matrix ( $R$).

$$T_{\mathrm{X}2\mathrm{C}}=R^{†}(TX+X^{†}T-X^{†}TX)R$$

$$V_{\mathrm{X}2\mathrm{C}}=R^{†}(V+\frac{1}{4c^2}{X}^{†}{W}^{\text{SF}}X)R$$

The coupling matrix ( $X=C^S(C^L{)}^{-1}$ ) is obtained from the large ($C^{\mathrm{L}}$) and small ($C^{\mathrm{S}}$) components of the $N$ positive energy solutions of the Dirac equation. The renormalization matrix $R=S^{-1/2}(S^{-1/2}\tilde{S}{S}^{-1/2}{)}^{-1/2}{S}^{1/2}$, depends on the modified overlap matrix $\tilde{S}=S+\frac{1}{2c^2}{X}^{†}TX$. The integrals $W_{\mu \nu }^{\mathrm{S}\mathrm{F}}=⟨{\chi }_{\mu }|\hat{p}\cdot (\hat{V}\hat{p})|{\chi }_{\nu }⟩$ can be easily computed as derivatives of the nuclear-electron attraction integrals with respect to nuclear coordinates. Existing nonrelativistic electronic structure code can be extended to include scalar relativistic effects treated with the X2C method by replacing nonrelativistic kinetic and potential energy with the corresponding X2C operators $T_{X2C}$ and $V_{X2C}$. It is important to note that fully uncontracted basis in needed for the construction of X2C Hamiltonian as Foldy-Wouthuysen (FW [FW:1950]) transformation is obtained in kinetically balance basis.

Keywords

RELATIVISTIC

Relativistic Hamiltonian type

• Type: string

• Possible Values: NO, X2C, ZORA

• Default: NO

BASIS_RELATIVISTIC

Auxiliary basis set for solving Dirac equation in X2C and DKH calculations. Defaults to decontracted orbital basis.

• Type: string

• Default: No Default