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:
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 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 ?
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.
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
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
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
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 [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 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. |
\(<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 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:
The eXact 2-Component Hamiltonian (X2X). Specification:
**HAMILTONIAN .X2CThe ZORA Hamiltonian. Specification:
**HAMILTONIAN .ZORA 0 0The scaled ZORA Hamiltonian. Specification:
**HAMILTONIAN .ZORA 0 1The 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:
B3LYP
PBE
PBE0
- The exercises are:
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\).
For AuH, explore the non-relativistic limit by recalculating spectroscopic constants setting \(c=4000\) a.u.
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
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
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 Å).