# 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 :math:<r^2>^{1/2}of spin-orbit partners down the series ? 2. Perform 4-component spin-free and non-relativistic HF-calculations on the heaviest atom of your series in order to investigate the effects of spin-orbit coupling and scalar relativistic effects. i) Make a table showing orbital energies and orbital sizes of all orbitals for each of the three Hamiltonians. ii) What orbitals contract (stabilize) and expand (destabilize) ? iii) Plot the densities of the outermost $$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 $\psi_{n\ell m}\left(\boldsymbol{r}\right)=R_{n \ell}\left(r\right)Y_{\ell m}\left(\theta,\phi\right)$ where appears the principal quantum number $$n$$, the azimuthal quantum number $$\ell$$ (orbital angular momentum) and the magnetic quantum number $$m$$. Orbitals are indicated using the letter code $$s,p,d,f,\ldots$$ for orbital angular momentum $$\ell=0,1,2,3,\ldots$$ . Complete specification requires multiplying these spatial functions with spin $$\alpha$$ or $$\beta$$ functions In the 2-component relativistic case, the orbitals have the general form $\psi_{n\kappa m_{j}}\left(\boldsymbol{r}\right)=R_{n\kappa}\left(r\right)\chi_{\kappa,m_{j}}\left(\theta,\phi\right)$ where $$\chi_{\kappa,m_{j}}$$ are 2-component angular functions. The quantum number $$\kappa$$ requires some explanation. Spin-orbit interaction leads to a coupling of orbital and spin angular momentum to form total angular momentum $$\mathbf{j} = \boldsymbol{\ell} + \mathbf{s}$$. The associated quantum number $$j$$ can take the values $$j=\ell+1/2$$ and $$j=|\ell-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 $$\kappa$$ 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}$$ $$\kappa$$ -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{split}\psi_{n\kappa m_{j}}\left(\boldsymbol{r}\right)=\left[\begin{array}{c}\psi^{L}\\\psi^{S}\end{array}\right] =\left[\begin{array}{c}R^{L}\left(r\right)\chi_{\kappa,m_{j}}\left(\theta,\phi\right)\\{\color{red}i}R^{S}\left(r\right)\chi_{-\kappa,m_{j}}\left(\theta,\phi\right)\end{array}\right]\end{split}$ 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 $$1s^2s^22p^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 $$(SS|SS)$$ two-electron integrals are replaced by a simple energy correction see 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 $$2p_{1/2}$$ and $$2p_{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 $$\mbox{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 $$\langle i|x^|i\rangle> for occupied orbitals :math:$$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.     $$^{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 $$2p_{1/2,1/2}$$, $$2p_{3/2,1/2}$$, and $$2p_{3/2,1/2}$$, respectively, where the second subscript corresponds to the value of $$|m_j|$$. We see that the $$2p_{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 $$2p_{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: 1. The eXact 2-Component Hamiltonian (X2X). Specification: **HAMILTONIAN .X2C  2. The ZORA Hamiltonian. Specification: **HAMILTONIAN .ZORA 0 0  3. The scaled ZORA Hamiltonian. Specification: **HAMILTONIAN .ZORA 0 1  4. 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 $$\omega_e$$ and anharmonic constants $$\omega_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-relativstic 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 $$\omega_e$$ and anharmonic constant $$\omega_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^{SS}_{inter}=\frac{q^S_{Cu}q^S_H}{R_{Cu-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_{2v}$$ 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:

1.40    -1655.6709986236738
1.41    -1655.6713227352741
1.42    -1655.6715820283903
1.43    -1655.6717794281951
1.44    -1655.6719177206721
1.45    -1655.6719995609726
1.46    -1655.6720274836848
1.47    -1655.6720039141965
1.48    -1655.6719311794373
1.49    -1655.6718115179062


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

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

$\gamma = \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 $$1s$$ 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 Å).