Introduction¶

In order to run a DIRAC calculation you need an input file of arbitrary name, e.g. dirac.inp which will be setup in the following and which contains the computational directives. In addition, a second file is needed providing the molecular structure. DIRAC can handle two formats but for convenience we here focus on the standard XYZ molecular input file, named for example molecule.xyz. The run-command inside your cluster run script looks then as follows:

$ pam --noarch --inp=dirac.inp --mol=molecule.xyz

where the DIRAC output is re-directed to the file dirac_molecule.out. Make sure that you rename the output file with a proper name if you want to take a look at what you did back home.

More information on running DIRAC calculations can be found in the First steps in DIRAC territory section of this documentation.

If not stated otherwise we will always use a Gaussian nuclear model in DIRAC.

Note

It is important that you understand the inputs. Before running calculations, carefully go through the input with reference to the manual that you find on the DIRAC pages

Note

You are expected to make a report presenting your results, also including some theoretical background as well as computational details.

Atomic calculations¶

Note

Exercise A In these exercises we look at relativity in atomic systems. Three columns of the periodic table will be explored: 1. From group 2: Mg, Sr and Ra 2. From group 12: Zn, Cd and Hg 3. From group 18: Ne, Kr and Rn

The exercises are: 1. Perform 4-component HF-calculations on your series of three atoms. We focus on the outermost $p$ orbitals: i) How does spin-orbit splitting vary as a function of nuclear charge Z ? ii) How do spin-orbit interaction affect orbital sizes $< r^{2} >^{1 / 2} ‘ o f s p i n - o r b i t p a r t n e r s d o w n t h e s e r i e s ? 2. P e r f o r m 4 - c o m p o n e n t * s p i n - f r e e * a n d * n o n - r e l a t i v i s t i c * H F - c a l c u l a t i o n s o n t h e h e a v i e s t a t o m o f y o u r s e r i e s i n o r d e r t o i n v e s t i g a t e t h e e f f e c t s o f s p i n - o r b i t c o u p l i n g a n d s c a l a r r e l a t i v i s t i c e f f e c t s . i) M a k e a t a b l e s h o w i n g o r b i t a l e n e r g i e s a n d o r b i t a l s i z e s o f a l l o r b i t a l s f o r e a c h o f t h e t h r e e H a m i l t o n i a n s . i i) W h a t o r b i t a l s c o n t r a c t (s t a b i l i z e) a n d e x p a n d (d e s t a b i l i z e) ? i i i) P l o t t h e d e n s i t i e s o f t h e * o u t e r m o s t * : m a t h : ‘ p$ orbitals with the three Hamiltonians. 3. Perform 2-component relativistic HF-calculations on the heaviest atom of your series using the following Hamiltonians: i) X2C ii) second-order Douglas-Kroll-Hess iii) ZORA and iv) scaled ZORA. Make a table of energies of all orbitals obtained with these Hamiltonians comparing with those obtained with the 4-component Dirac-Coulomb Hamiltonian. Helpful information for these exercises is given below !

We shall start gently by looking at atoms. It will be useful to recall the form and notation for atomic orbitals. In the non-relativistic case atomic orbitals have the general form

ψ_{n ℓ m} (r) = R_{n ℓ} (r) Y_{ℓ m} (θ, ϕ)

where appears the principal quantum number $n$ , the azimuthal quantum number $ℓ$ (orbital angular momentum) and the magnetic quantum number $m$ . Orbitals are indicated using the letter code $s, p, d, f, \dots$ for orbital angular momentum $ℓ = 0, 1, 2, 3, \dots$ . Complete specification requires multiplying these spatial functions with spin $α$ or $β$ functions

In the 2-component relativistic case, the orbitals have the general form

ψ_{n κ m_{j}} (r) = R_{n κ} (r) χ_{κ, m_{j}} (θ, ϕ)

where $χ_{κ, m_{j}}$ are 2-component angular functions. The quantum number $κ$ requires some explanation. Spin-orbit interaction leads to a coupling of orbital and spin angular momentum to form total angular momentum $j = ℓ + s$ . The associated quantum number $j$ can take the values $j = ℓ + 1 / 2$ and $j = | ℓ - 1 / 2 |$ . However, just giving $j$ does not specify an orbital uniquely, since it can for instance not distinguish $s_{1 / 2}$ and $p_{1 / 2}$ orbitals (both have $j = 1 / 2$ ). This is possible with the quantum number $κ$ according to the following table:

	$s_{1 / 2}$	$p_{1 / 2}$	$p_{3 / 2}$	$d_{3 / 2}$	$d_{5 / 2}$	$f_{5 / 2}$	$f_{7 / 2}$
$κ$	-1	+1	-2	+2	-3	+3	-4

4-component atomic orbitals have separate radial and angular parts for both the large and small components

\begin{array}{r} ψ_{n κ m_{j}} (r) = [\begin{array}{c} ψ^{L} \\ ψ^{S} \end{array}] = [\begin{array}{c} R^{L} (r) χ_{κ, m_{j}} (θ, ϕ) \\ i R^{S} (r) χ_{- κ, m_{j}} (θ, ϕ) \end{array}] \end{array}

with the notation for orbitals dictated by the quantum number of the large component.

With these introductory remarks in mind, we are ready to do our first calculation.

We will run a Hartree-Fock calculation on the neon atom. We shall use the input file Ne.inp (download):

**DIRAC
.WAVE FUNCTION
.ANALYZE
.PROPERTIES
**INTEGRAL
*READIN
.UNCONT
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
4 6
**ANALYZE
.MULPOP
**PROPERTIES
*EXPECTATION VALUE
.ORBANA
.OPERATOR
 'XXSECMOM'
.OPERATOR
 'YYSECMOM'
.OPERATOR
 'ZZSECMOM'
**MOLECULE
*BASIS
.DEFAULT
 dyall.2zp
**END OF

and the xyz-file Ne.xyz (download):

1

Ne 0.0 0.0 0.0

Neon has the ground-state configuration $1 s^{2} s^{2} 2 p^{6}$ . Note that under the keyword (.CLOSED SHELL) we specify the occupation of gerade and ungerade orbitals separately, that is:

.CLOSED SHELL
4 6

For this particular calculation we have chosen to use the dyall.2zp, which is a double zeta basis set for SCF calculations:

*BASIS
.DEFAULT
dyall.2zp

The minimal commmand for running this calculations is:

pam --inp=Ne.inp --mol=Ne.xyz

Looking in the output, we find information about the SCF-cycles and the final energy

                                   SCF - CYCLE
                                   -----------

* Convergence on norm of error vector (gradient).
  Desired convergence:1.000D-07
  Allowed convergence:1.000D-06

* ERGVAL - convergence in total energy
* FCKVAL - convergence in maximum change in total Fock matrix
* EVCVAL - convergence in error vector (gradient)
--------------------------------------------------------------------------------------------------------------------------------
           Energy               ERGVAL    FCKVAL    EVCVAL      Conv.acc    CPU          Integrals   Time stamp
--------------------------------------------------------------------------------------------------------------------------------
It.    1    -65.37954964130      0.00D+00  0.00D+00  0.00D+00               0.00641500s   Atom. scrpot   Fri Sep 24
It.    2    -128.6584266675      6.33D+01  7.20D+00  8.91D-01               0.11870000s   LL SL          Fri Sep 24
It.    3    -128.6830824465      2.47D-02  3.01D-01  1.85D-01   DIIS   2    0.06050000s   LL SL          Fri Sep 24
It.    4    -128.6851704671      2.09D-03 -1.49D-01  7.90D-02   DIIS   3    0.06010000s   LL SL          Fri Sep 24
It.    5    -128.6857206713      5.50D-04  4.98D-02  1.50D-03   DIIS   4    0.05580000s   LL SL          Fri Sep 24
It.    6    -128.6857208810      2.10D-07 -3.73D-04  9.33D-05   DIIS   5    0.05480000s   LL SL          Fri Sep 24
It.    7    -128.6857208819      8.80D-10 -3.56D-05  5.95D-06   DIIS   6    0.04900000s   LL SL          Fri Sep 24
It.    8    -128.6857208819      7.67D-12  1.37D-06  8.69D-07   DIIS   7    0.05550000s   LL SL          Fri Sep 24
It.    9    -128.6857208819      5.12D-13  2.39D-07  1.80D-07   DIIS   8    0.04830000s   LL SL          Fri Sep 24
It.   10    -128.6857208819     -1.71D-13 -5.01D-08  2.18D-09   DIIS   8    0.04780000s   LL SL          Fri Sep 24
--------------------------------------------------------------------------------------------------------------------------------
* Convergence after   10 iterations.
* Average elapsed time per iteration:
      No 2-ints    :    0.01110000s
      LL SL        :    0.06116667s


                                   TOTAL ENERGY
                                   ------------

   Electronic energy                        :    -128.68572088190120

   Other contributions to the total energy
   Nuclear repulsion energy                 :       0.00000000000000
   SS Coulombic correction                  :       0.00000000000000

   Sum of all contributions to the energy
   Total energy                             :    -128.68572088190120

The entry SS Coulombic correction refers to a default approximation in which the calculations of the expensive $(S S | S S)$ two-electron integrals are replaced by a simple energy correction see [Visscher1997a] for further details. It is strictly zero for atoms.

The final energy is followed by a list of orbital eigenvalues. In the present calculation we have activated Mulliken population analysis and find further information about the occupied orbitals in the output

                                Fermion ircop E1g
                                -----------------


* Electronic eigenvalue no.   1: -32.815238678865   (Occupation : f = 1.0000)  s 1/2;  1/2
============================================================================================

* Gross populations greater than 0.00010

Gross     Total   |    L Ag Ne s      B3uNe _small   B2uNe _small   B1uNe _small
-----------------------------------------------------------------------------------
 alpha    0.9992  |      0.9988         0.0000         0.0000         0.0004
 beta     0.0008  |      0.0000         0.0004         0.0004         0.0000

* Electronic eigenvalue no.   2: -1.9345569982051   (Occupation : f = 1.0000)  s 1/2;  1/2
============================================================================================

* Gross populations greater than 0.00010

Gross     Total   |    L Ag Ne s
--------------------------------------
 alpha    0.9999  |      0.9999
 beta     0.0001  |      0.0000


** Total gross population of fermion ircop E1g **

Gross      Total    |    L Ag Ne s      B3uNe _small   B2uNe _small   B1uNe _small
-----------------------------------------------------------------------------------
 total     4.00000  |      3.99725        0.00092        0.00092        0.00092


				Fermion ircop E1u
                                -----------------

* Electronic eigenvalue no.   1: -0.8511723193420   (Occupation : f = 1.0000)  p 1/2;  1/2
============================================================================================

* Gross populations greater than 0.00010

Gross     Total   |    L B3uNe px     L B2uNe py     L B1uNe pz     Ag Ne _small
-----------------------------------------------------------------------------------
 alpha    0.3334  |      0.0000         0.0000         0.3333         0.0001
 beta     0.6666  |      0.3333         0.3333         0.0000         0.0000

* Electronic eigenvalue no.   2: -0.8466879891865   (Occupation : f = 1.0000)  p 3/2; -3/2
============================================================================================

* Gross populations greater than 0.00010

Gross     Total   |    L B3uNe px     L B2uNe py
-----------------------------------------------------
 alpha    0.0001  |      0.0000         0.0000
 beta     0.9999  |      0.4999         0.4999

* Electronic eigenvalue no.   3: -0.8466879891700   (Occupation : f = 1.0000)  p 3/2;  1/2
============================================================================================

* Gross populations greater than 0.00010

Gross     Total   |    L B3uNe px     L B2uNe py     L B1uNe pz
--------------------------------------------------------------------
 alpha    0.6666  |      0.0000         0.0000         0.6666
 beta     0.3334  |      0.1666         0.1666         0.0000


** Total gross population of fermion ircop E1u **

Gross      Total    |    L B3uNe px     L B2uNe py     L B1uNe pz     Ag Ne _small
-----------------------------------------------------------------------------------
 total     6.00000  |      1.99977        1.99977        1.99977        0.00041


*** Total gross population ***

Gross      Total    |    L Ag Ne s      L B3uNe px     L B2uNe py     L B1uNe pz     Ag Ne _small   B3uNe _small   B2uNe _small
--------------------------------------------------------------------------------------------------------------------------------
 total    10.00000  |      3.99725        1.99977        1.99977        1.99977        0.00041        0.00092        0.00092

Note again the separation on gerade and ungerade orbitals. We shall focus on the valence $p$ orbitals of neon and the heavier homologues. In the present calculation we note that the $2 p_{1 / 2}$ and $2 p_{3 / 2}$ orbitals have energies -0.851172 $E_{h}$ and -0.846688 $E_{h}$ , respectively, corresponding to a spin-orbit splitting of 0.004484 $E_{h}$ , or 984.2 ${cm}^{- 1}$ (the converstion factor is 2.19474625E+5).

We have also chosen to have get some information about orbital sizes. We do not have direct access to radial expectation values $< r >$ , but can calculate the rms $< r^{2} >^{1 / 2}$ by calculating $< r^{2} >=< x^{2} > + < y^{2} > + < z^{2} >$ . This is specified in the input as:

*EXPECTATION VALUE
.ORBANA
.OPERATOR
'XXSECMOM'
.OPERATOR
'YYSECMOM'
.OPERATOR
'ZZSECMOM'

(see One-electron operators for syntax). We have also used the keyword .ORBANA to ask for contributions from individual orbitals. In the output we therefore for instance find

Operator XXSECMOM        :
 --------------------------



                    Expectation value for individual orbitals
                    -----------------------------------------


                Matrix element     Occ.
    E1g    1    0.11124959E-01     2.0000      0.22249918E-01
    E1g    2    0.320721858726     2.0000      0.641443717451
    E1u    1    0.405856904615     2.0000      0.811713809231
    E1u    2    0.489542686491     2.0000      0.979085372983
    E1u    3    0.326361790990     2.0000      0.652723581981
    -------------------------------------------------------------------------------------
    Total :                                    3.107216399820 a.u.      s0 = F   t0 = F
    -------------------------------------------------------------------------------------

One column gives the matrix elements $⟨ i | x^{|} i ⟩ > f o r o c c u p i e d o r b i t a l s : m a t h :$ varphi_i`, whereas the final column gives the matrix element multiplied with the occupation number of the orbital.

We are interested in the rms of the valence p orbitals. We make the following table

Orb.	$< x^{2} >$	$< y^{2} >$	$< z^{2} >$	$< r^{2} >$	$< r^{2} >^{1 / 2}$
E1u 1	0.405856904615	0.405856904615	0.405856904615	1.218	1.103
E1u 2	0.489542686491	0.489542686491	0.244771343246	1.223	1.106
E1u 3	0.326361790990	0.326361790990	0.571133134233	1.223	1.106

From the Mulliken analysis output we see that the three listed orbitals correspond to $2 p_{1 / 2, 1 / 2}$ , $2 p_{3 / 2, 1 / 2}$ , and $2 p_{3 / 2, 1 / 2}$ , respectively, where the second subscript corresponds to the value of $| m_{j} |$ . We see that the $2 p_{3 / 2}$ orbitals all have the same size ( $< r^{2} >^{1 / 2} = 1.106 a_{0}$ ) and are slightly bigger than the spin-orbit partner $2 p_{1 / 2}$ ( $< r^{2} >^{1 / 2} = 1.103 a_{0}$ ).

It may also be of interest to have a look at the shape of the atomic orbitals. Here things are somewhat more complicated than in the non-relativistic case since we have seen that the orbitals are complex 4-component vector functions. What is possible is to plot orbital densities, see Plotting orbital densities for instructions.

So far we have done 4-component relativistic calculations. It is possible to successively turn off spin-orbit interaction and scalar relativistic effects using the keywords .SPINFREE and .LEVY-LEBLOND. The latter calculation is based on the Lévy-Leblond Hamiltonian, which is a 4-component non-relativistic Hamiltonian.

It is also possible to explore various approximate 2-component Hamiltonians. We shall be interested in the following ones:

The eXact 2-Component Hamiltonian (X2X). Specification:
**HAMILTONIAN
.X2C
The ZORA Hamiltonian. Specification:
**HAMILTONIAN
.ZORA
0  0
The scaled ZORA Hamiltonian. Specification:
**HAMILTONIAN
.ZORA
0  1
The second-order Douglas-Kroll-Hess Hamiltonian Specification:
**HAMILTONIAN
.DKH2

Molecular calculations¶

Note

Exercise M : In these exercises we look at relativistic effects on the spectroscopic constants of noble metal monohydrides. Experimental data can be found at these NIST pages by giving the chemical formula and selecting Constants of diatomic molecules. Note that we are looking at spectroscopic constants for the ground state, normally denoted X and found at the bottom of the table. Three different DFT functionals will be explored: 1. B3LYP 2. PBE 3. PBE0

The exercises are: 1. Perform 4-component DFT-calculations on the monohydrides of the noble metals Cu, Ag and Au in order to investigate the effects of spin-orbit coupling and scalar relativity (again using the keywords .SPINFREE and .LEVY-LEBLOND) on equilibrium distances $r_{e}$ , harmonic frequencies $ω_{e}$ and anharmonic constants $ω_{e} x_{e}$ . 2. For AuH, explore the non-relativistic limit by recalculating spectroscopic constants setting $c = 4000$ a.u. 3. For CuH, explore the ultra-relativistic limit by recalculating spectroscopic constants setting $c = 20$ a.u. Helpful information for these exercises is given below !

As a sample calculation, let us look at CuH. These are simple diatomic species, so we shall simply obtain spectroscopic constants by calculating the energy at a selection of interatomic distances about the expected minimum. We can then calculate the equilibrium distance $r_{e}$ , harmonic frequency $ω_{e}$ and anharmonic constant $ω_{e} x_{e}$ using the utility program TWOFIT. The manual pages provide theoretical background. The executable twofit.x is in the same directory as the DIRAC executable.

We will run a DFT calculation on CuH, here using the B3LYP functional. We shall use the input file CuH.inp (download):

**DIRAC
.WAVE FUNCTION
**INTEGRALS
*READIN
.UNCONT
**HAMILTONIAN
.DFT
B3LYP
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
30
**MOLECULE
*BASIS
.DEFAULT
dyall.3zp
**END OF

and the xyz-file CuH.xyz (download):

2

Cu 0.0 0.0 0.0
H  0.0 0.0 1.46

The minimal commmand for running this calculations is:

pam --inp=CuH.inp --mol=CuH.xyz

Looking in the output, we find information about the SCF-cycles and the final energy

                                   SCF - CYCLE
                                   -----------

* Convergence on norm of error vector (gradient).
  Desired convergence:1.000D-07
  Allowed convergence:1.000D-06

* ERGVAL - convergence in total energy
* FCKVAL - convergence in maximum change in total Fock matrix
* EVCVAL - convergence in error vector (gradient)
--------------------------------------------------------------------------------------------------------------------------------
           Energy               ERGVAL    FCKVAL    EVCVAL      Conv.acc    CPU          Integrals   Time stamp
--------------------------------------------------------------------------------------------------------------------------------
It.    1    -956.5681555153      0.00D+00  0.00D+00  0.00D+00               0.17635000s   Atom. scrpot   Sat Sep 25
It.    2    -1655.540673820      6.99D+02 -1.36D+02  6.00D+00               2.29200000s   LL SL          Sat Sep 25
It.    3    -1655.461930487     -7.87D-02 -5.13D-01  1.02D+00   DIIS   2    1.83020000s   LL SL          Sat Sep 25
It.    4    -1655.433353314     -2.86D-02  1.04D+00  1.59D+00   DIIS   3    1.83300000s   LL SL          Sat Sep 25
It.    5    -1655.624584139      1.91D-01 -1.03D+00  6.43D-01   DIIS   4    1.84550000s   LL SL          Sat Sep 25
It.    6    -1655.669690995      4.51D-02  3.81D-01  1.49D-01   DIIS   5    1.84250000s   LL SL          Sat Sep 25
It.    7    -1655.671844620      2.15D-03 -1.06D-01  4.39D-02   DIIS   6    1.83240000s   LL SL          Sat Sep 25
It.    8    -1655.672026017      1.81D-04  2.30D-02  2.57D-03   DIIS   7    1.81630000s   LL SL          Sat Sep 25
It.    9    -1655.672027185      1.17D-06  2.69D-04  1.64D-03   DIIS   8    1.80900000s   LL SL          Sat Sep 25
It.   10    -1655.672027456      2.71D-07  7.47D-04  3.34D-04   DIIS   9    1.79780000s   LL SL          Sat Sep 25
It.   11    -1655.672027483      2.73D-08  1.05D-04  4.66D-05   DIIS   9    1.80330000s   LL SL          Sat Sep 25
It.   12    -1655.672027484      4.85D-10  1.24D-05  1.21D-05   DIIS   9    1.79280000s   LL SL          Sat Sep 25
It.   13    -1655.672027484      2.59D-11  6.16D-06  1.12D-06   DIIS   9    1.74590000s   LL SL          Sat Sep 25
It.   14    -1655.672027484     -5.68D-12 -3.45D-07  2.32D-07   DIIS   9    1.73510000s   LL SL          Sat Sep 25
It.   15    -1655.672027484      7.96D-12 -7.70D-08  5.29D-08   DIIS   9    1.63110000s   LL SL          Sat Sep 25
--------------------------------------------------------------------------------------------------------------------------------
* Convergence after   15 iterations.
* Average elapsed time per iteration:
      No 2-ints    :    0.19430000s
      LL SL        :    1.82906429s


                                   TOTAL ENERGY
                                   ------------

   Electronic energy                        :    -1666.1830818801755

   Other contributions to the total energy
   Nuclear repulsion energy                 :       10.5110541891692
   SS Coulombic correction                  :        0.0000002073215

   Sum of all contributions to the energy
   Total energy                             :    -1655.6720274836848

We may now see that the entry SS Coulombic correction is non-zero, albeit very small. In this particular case the energy correction is calculated as

E_{i n t e r}^{S S} = \frac{q_{C u}^{S} q_{H}^{S}}{R_{C u - H}}

where $q^{S}$ is the contribution of the small components to the atomic charge.

In passing we note that the output contains an xyz-file

  Cartesian coordinates in XYZ format (Angstrom)
  ----------------------------------------------

    2

Cu     0.0000000000   0.0000000000  -0.0230135093
H      0.0000000000   0.0000000000   1.4369864907

which is modified compared to the input file. This is because DIRAC checked the original xyz-input and found linear symmetry

                          SYMGRP:Point group information
                          ------------------------------

Full group is:  C(oo,v)
Represented as: C2v

DIRAC uses linear supersymmetry in these calculations; this means that the underlying integrals are adapted to the $C_{2 v}$ subgroup, but the Fock matrix has the full block structure due to linear symmetry. In the same manner, the atomic calculations above uses atomic supersymmetry.

We are now ready to scan the potential surface of CuH in the vicinity of the expected minimum. A pedestrian way is to run the calculation point by point, modfifying the input xyz-file for each new geometry. However, it is also possible to automatize the scan, e.g. (using bash)

for dist in 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.49
do
cat <<EOF > ${dist}.xyz
2

Cu 0.0 0.0 0.0
H  0.0 0.0 ${dist}
EOF
pam --inp=CuH --mol=${dist}.xyz
done

Total energies can then be extracted for instance using:

grep 'Total energy                             :' CuH_*out > CuH.pot

After editing the file CuH.pot (download) reads:

40    -1655.6709986236738
41    -1655.6713227352741
42    -1655.6715820283903
43    -1655.6717794281951
44    -1655.6719177206721
45    -1655.6719995609726
46    -1655.6720274836848
47    -1655.6720039141965
48    -1655.6719311794373
49    -1655.6718115179062

A polynomial fit of order 6, using TWOFIT (output here), gives the spectroscopic constants $r_{e} = 1.4604$ Å, $ω_{e} = 1958.2 {cm}^{- 1}$ and $ω_{e} x_{e} = 34.6 {cm}^{- 1}$ . From the NIST page we find $r_{e} = 1.46263$ Å, $ω_{e} = 1941.26 {cm}^{- 1}$ and $ω_{e} x_{e} = 37.51 {cm}^{- 1}$ , so our results are not bad !

We have seen that relativistic effects are dictated by the Lorentz factor

γ = \frac{1}{\sqrt{1 - \frac{v^{2}}{c^{2}}}}

where $v$ is the speed of the particle and $c$ is the speed of light. We have also seen that for hydrogen-like atoms the speed of the $1 s$ orbital is $Z$ (nuclear charge) in atomic units, to be compared with the speed of light which is around 137 in the same units. The non-relativistic limit can be attained by setting the speed of light to infinity. This is not possible on a computer, so we settle for a large, but finite value. It is also possible to investigate the ultra-relativistic case by lowering the speed of light. In DIRAC the speed of light can be adjusted using the .CVALUE keyword.

Note

When extracting spectroscopic constants using TWOFIT it is important that the minimum $r_{e}$ is in the interval of selected points, ideally somewhere in the middle, so you may have to add further points to your scan to achieve this. Relativistic effects are important, so when searching for an unknown minimum, you may first make a rather coarse grid of points (spacing 0.1 Å) and use TWOFIT to guide you in the right direction, where you can put a finer grid (spacing 0.01 Å).