4.2 Heavy elements ($Z\gtrsim$30): Modifications for scalar relativity

Actually, this section could easily apply to all elements. As outlined in Section 3.8, it seems that the relativistic atomic_zora scalar keyword and the underlying “atomic ZORA” approach as implemented in FHI-aims (Equations (55) and (56) of Ref. [36]) perform on par with the best available scalar-relativistic benchmark methods (see Refs. [195, 148]). Just use this approach in production, unless there is a specific reason not do do so (a few methods in FHI-aims that are still being implemented may initially only support a non-relativistic treatment).

Note that total energy differences will be very large between different relativistic treatments, so it is best to stick to a single level of scalar relativity for all calculations.

With H and O, the simple H₂O testcase of the preceding Section 4.1 involves only elements so light that relativistic effects on their total energies can still be ignored in production calculations.

The rule of thumb that relativistic effects do not matter much for quantities derived from total energies (binding energies, geometries) holds up to elements number $Z$=20-30 (Ca or Zn), depending on the reqired accuracy. For example, the DFT-LDA lattice parameter for fcc Cu is 3.54 Å in the non-relativistic case, but 3.52 Å if relativistic effects are included. In any case, relativistic effects should be accounted for in accurate calculations beyond these elements. As described in more detail in Ref. [36], this process is handled by FHI-aims at different levels of approximation, with some overhead resulting compared to the simple non-relativistic case.

In a nutshell, the physical impact of relativity for heavy elements is noticeable not just as a total-energy offset, but importantly in total energy differences, and in particular impacts also geometries. The underlying reason is that core and valence orbitals “see” the near-nuclear region (where relativistic effects are most important) differently, depending on their angular momentum $l$. Somewhat simplistically put, the increased relativistic mass of the electron near the nucleus results in a tendency for all orbitals to “shrink” compared to the non-relativistic case, but to a different degree for different orbitals. The shrinkage thus changes not just atom sizes, but also the nature of bonding itself. A detailed discussion of relativistic effects is beyond the scope of this manual (and can be found in many excellent reviews, e.g., Ref. [256]) – but bear in mind that relativistic effects should not simply be shrugged off!

As noted in Section 3.8, we recommend “atomic ZORA” as the blanket default where possible, invoked by the keyword

relativistic atomic_zora scalar

Additionally, the effects of spin-orbit coupling may be included in energy levels / band structures, densities of states, etc. post-selfconsistently, i.e., after the self-consistent calculation is complete. The keyword to do this is:

include_spin_orbit

That’s it.

We illustrate the practical use of atomic ZORA and spin-orbit coupling for the Au dimer in DFT-LDA, found in the testcases/Au_dimer directory. Importantly, this follows the always recommended two-step approach for relaxations: First, a light prerelaxation (saves much time for steps of the relaxation algorithm that will be completely irrelevant for the final result) and second, a tight post-relaxation for final results.

The testcase here uses semilocal DFT, which is relatively cheap. For the much more costly hybrid density functionals, tight settings can be prohibitively expensive, and intermediate settings (where available) are often completely sufficient.

In the subdirectory relax_light, a quick but sufficiently accurate prerelaxation is set up, at the atomic_zora level. The control.in file is set up to use relativistic atomic_zora scalar and relax_geometry bfgs 1.e-2 to converge forces down to 10^-2 eV/Å. Since this is intended to be a quick prerelaxation run (starting with an arbitrary separation of 3 Å of both atoms), the “light” species default settings for Au are used for the species subsettings. Compared to the much more accurate “tight” settings,

• •

  the radial and angular integration grids are significantly reduced,

• •

  the Hartree potential expansion is capped at l_hartree=4,

• •

  the cutoff potential onset and width in cut_pot are reduced to a minimum that we still consider safe (3.5 Å / 1.5 Å, respectively),

• •

  the basis set is the spdf section of tier 1 only.

. Among these settings, the reduction of the basis set to spdf has the biggest impact on the resulting equilibrium geometry, $d$=2.464 Å. (The difference to the tight result below is quite minor – 0.012 Å.)

After this relaxation run is complete, a file geometry.in.next_step is written out by FHI-aims. This file contains the final, relaxed geometry from the “light” prerelaxation run just performed, as well as the Hessian matrix estimate created by the relax_geometry bfgs relaxation algorithm. This file can serve as an improved guess for geometry.in in the next, typically more expensive “tight” post-relaxation step.

In the subdirectory postrelax_tight, the final geometry from the previous step is used as the starting point for an accurate post-relaxation using the “tight” species_defaults settings. Note that, while all other settings are tightly converged at this level, the basis set convergence should normally still be tested explicitly. The tier 1 level (without the $h$ function) used here for Au is quite sufficient for an accurate result. The binding distance converges after two additional relaxation steps, at $d$=2.452 Å.

We also added in spin-orbit coupling for completeness using the include_spin_orbit keyword.

In conclusion, this section illustrates how a quick pre-relaxation followed by a safe calculation of the relevant energy differences can be combined. The key point in this endeavour is to ensure that the final relaxed geometry is accurate with as little computational overhead as possible.