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 runcommand inside your cluster run script looks then as follows:
$ pam noarch inp=dirac.inp mol=molecule.xyz
where the DIRAC output is redirected 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:
From group 2: Mg, Sr and Ra
From group 12: Zn, Cd and Hg
From group 18: Ne, Kr and Rn
 The exercises are:
Perform 4component HFcalculations on your series of three atoms. We focus on the outermost \(p\) orbitals: i) How does spinorbit splitting vary as a function of nuclear charge Z ? ii) How do spinorbit interaction affect orbital sizes :math:`<r^2>^{1/2}`of spinorbit partners down the series ?
Perform 4component spinfree and nonrelativistic HFcalculations on the heaviest atom of your series in order to investigate the effects of spinorbit 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.
Perform 2component relativistic HFcalculations on the heaviest atom of your series using the following Hamiltonians: i) X2C ii) secondorder DouglasKrollHess iii) ZORA and iv) scaled ZORA. Make a table of energies of all orbitals obtained with these Hamiltonians comparing with those obtained with the 4component DiracCoulomb 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 nonrelativistic case atomic orbitals have the general form
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 2component relativistic case, the orbitals have the general form
where \(\chi_{\kappa,m_{j}}\) are 2component angular functions. The quantum number \(\kappa\) requires some explanation. Spinorbit 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=\ell1/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 
4component atomic orbitals have separate radial and angular parts for both the large and small components
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 HartreeFock 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 xyzfile Ne.xyz (download
):
1
Ne 0.0 0.0 0.0
Neon has the groundstate 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 SCFcycles and the final energy
SCF  CYCLE

* Convergence on norm of error vector (gradient).
Desired convergence:1.000D07
Allowed convergence:1.000D06
* 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.91D01 0.11870000s LL SL Fri Sep 24
It. 3 128.6830824465 2.47D02 3.01D01 1.85D01 DIIS 2 0.06050000s LL SL Fri Sep 24
It. 4 128.6851704671 2.09D03 1.49D01 7.90D02 DIIS 3 0.06010000s LL SL Fri Sep 24
It. 5 128.6857206713 5.50D04 4.98D02 1.50D03 DIIS 4 0.05580000s LL SL Fri Sep 24
It. 6 128.6857208810 2.10D07 3.73D04 9.33D05 DIIS 5 0.05480000s LL SL Fri Sep 24
It. 7 128.6857208819 8.80D10 3.56D05 5.95D06 DIIS 6 0.04900000s LL SL Fri Sep 24
It. 8 128.6857208819 7.67D12 1.37D06 8.69D07 DIIS 7 0.05550000s LL SL Fri Sep 24
It. 9 128.6857208819 5.12D13 2.39D07 1.80D07 DIIS 8 0.04830000s LL SL Fri Sep 24
It. 10 128.6857208819 1.71D13 5.01D08 2.18D09 DIIS 8 0.04780000s LL SL Fri Sep 24

* Convergence after 10 iterations.
* Average elapsed time per iteration:
No 2ints : 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 \((SSSS)\) twoelectron 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 \(2p_{1/2}\) and \(2p_{3/2}\) orbitals have energies 0.851172 \(E_h\) and 0.846688 \(E_h\), respectively, corresponding to a spinorbit 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 Oneelectron 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.11124959E01 2.0000 0.22249918E01
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 ix^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. 
\(<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 \(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 spinorbit 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 nonrelativistic case since we have seen that the orbitals are complex 4component vector functions. What is possible is to plot orbital densities, see Plotting orbital densities for instructions.
So far we have done 4component relativistic calculations. It is possible to successively turn off spinorbit interaction and scalar relativistic effects using the keywords .SPINFREE and .LEVYLEBLOND. The latter calculation is based on the LévyLeblond Hamiltonian, which is a 4component nonrelativistic Hamiltonian.
It is also possible to explore various approximate 2component Hamiltonians. We shall be interested in the following ones:
The eXact 2Component Hamiltonian (X2X). Specification:
**HAMILTONIAN .X2CThe ZORA Hamiltonian. Specification:
**HAMILTONIAN .ZORA 0 0The scaled ZORA Hamiltonian. Specification:
**HAMILTONIAN .ZORA 0 1The secondorder DouglasKrollHess 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:
B3LYP
PBE
PBE0
 The exercises are:
Perform 4component DFTcalculations on the monohydrides of the noble metals Cu, Ag and Au in order to investigate the effects of spinorbit coupling and scalar relativity (again using the keywords .SPINFREE and .LEVYLEBLOND) on equilibrium distances \(r_e\), harmonic frequencies \(\omega_e\) and anharmonic constants \(\omega_e x_e\).
For AuH, explore the nonrelativistic limit by recalculating spectroscopic constants setting \(c=4000\) a.u.
For CuH, explore the ultrarelativstic 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 xyzfile 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 SCFcycles and the final energy
SCF  CYCLE

* Convergence on norm of error vector (gradient).
Desired convergence:1.000D07
Allowed convergence:1.000D06
* 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.87D02 5.13D01 1.02D+00 DIIS 2 1.83020000s LL SL Sat Sep 25
It. 4 1655.433353314 2.86D02 1.04D+00 1.59D+00 DIIS 3 1.83300000s LL SL Sat Sep 25
It. 5 1655.624584139 1.91D01 1.03D+00 6.43D01 DIIS 4 1.84550000s LL SL Sat Sep 25
It. 6 1655.669690995 4.51D02 3.81D01 1.49D01 DIIS 5 1.84250000s LL SL Sat Sep 25
It. 7 1655.671844620 2.15D03 1.06D01 4.39D02 DIIS 6 1.83240000s LL SL Sat Sep 25
It. 8 1655.672026017 1.81D04 2.30D02 2.57D03 DIIS 7 1.81630000s LL SL Sat Sep 25
It. 9 1655.672027185 1.17D06 2.69D04 1.64D03 DIIS 8 1.80900000s LL SL Sat Sep 25
It. 10 1655.672027456 2.71D07 7.47D04 3.34D04 DIIS 9 1.79780000s LL SL Sat Sep 25
It. 11 1655.672027483 2.73D08 1.05D04 4.66D05 DIIS 9 1.80330000s LL SL Sat Sep 25
It. 12 1655.672027484 4.85D10 1.24D05 1.21D05 DIIS 9 1.79280000s LL SL Sat Sep 25
It. 13 1655.672027484 2.59D11 6.16D06 1.12D06 DIIS 9 1.74590000s LL SL Sat Sep 25
It. 14 1655.672027484 5.68D12 3.45D07 2.32D07 DIIS 9 1.73510000s LL SL Sat Sep 25
It. 15 1655.672027484 7.96D12 7.70D08 5.29D08 DIIS 9 1.63110000s LL SL Sat Sep 25

* Convergence after 15 iterations.
* Average elapsed time per iteration:
No 2ints : 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 nonzero, albeit very small. In this particular case the energy correction is calculated as
where \(q^S\) is the contribution of the small components to the atomic charge.
In passing we note that the output contains an xyzfile
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 xyzinput 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 xyzfile 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
where \(v\) is the speed of the particle and \(c\) is the speed of light. We have also seen that for hydrogenlike 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 nonrelativistic 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 ultrarelativistic 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 Å).