Geometric properties: NO_3\cdot radical

6.2. Geometric properties: NO\(_3\cdot\) radical#

In this set of exercises, you will assess the performance of various state-of-the-art density functionals in the prediction of geometric properties, and you will again compare your results to both wavefunction theory and experimental data.

6.2.1. Structural Parameters of NO\(_3\cdot\): Performance of DFT, HF and MP2#

Before diving into the theoretical aspects of exchange-correlation functionals, you will put a representative selection of functionals to good use in a notoriously tricky system, the NO\(_3\cdot\) radical. Nitrite radicals are highly reactive species that are rapidly destroyed by sunlight, but as the sun sets, they start to play an important role in chemical transformations (in what is called the night-time chemistry of the atmosphere). [Ref: Phys. Chem. Chem. Phys., 2014, 16, 19437-19445]

There exist various experimental and theoretical studies of NO\(_3\cdot\) , with the experiment indicating a fully symmetric D\(_3^h\) structure with equal N-O bond lengths of 1.24 Å and O-N-O bond angles of 120\(^\circ\).

In this exercise, you will be comparing the performance of various DFT exchange-correlation functionals, Hartree-Fock theory and MP2 in predicting these structural parameters.

First, we’ll import the required modules and set Psi4 parameters.

import psi4
import pandas as pd

import sys
sys.path.append("..")

from helpers import *

psi4.set_num_threads(2)
psi4.set_memory('2 GB')

  Threads set to 2 by Python driver.

  Memory set to   1.863 GiB by Python driver.

2000000000

This is our starting geometry.

radical = psi4.geometry("""
symmetry c1
0 2
 N     0.000000     0.000000     0.000000
 O     1.400000     0.000000     0.000000
 O    -1.000000     1.000000     0.000000
 O    -1.000000    -1.000000     0.000000
""")

drawXYZ(radical)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

We will perform calculations at HF, MP2 and DFT level. We will also visualize MP2 and DFT orbitals and compare them (see Visualizing the orbitals subsections of MP2 and DFT sections).

6.2.1.1. HF#

First, we’ll do a Hartree Fock calculation. As this is a difficult system we make use of a second order SCF (SOSCF) convergence method implemented in Psi4.

SOSCF makes use of the fact that the energy that we want to minimize in the SCF calculation has a gradient (and a hessian) both with respect to geometry changes and orbital changes. In a standard SCF calculation as you get close to convergence the orbital changes will be very small. When using SOSCF, Psi4 uses the orbital hessian to make the next orbital step when performing the SCF calculation. In the code, the energy is written as a taylor expansion of orbitals and truncated for all terms higher than second order. This helps in converging tricky systems but requires that one is already close to a minimum. So by setting soscf to true by default the soscf procedure is switched on when the gradient RMS is below \(1.0 \times 10^{-2}\). By specifying soscf_max_iter we increase the number of SOSCF steps from the default of 5 to 20.

method = "hf"

psi4.core.set_output_file(f'radical_opt_{method}.log', False)
psi4.set_options({'reference':'uhf','guess':'gwh', 'soscf':True, 'soscf_max_iter':20,'maxiter':200})
E_hf = psi4.optimize('hf/6-31+G*', molecule=radical)
E_hf

6.2.1.2. MP2#

method = "mp2"


psi4.core.clean_options()
psi4.core.set_output_file(f'radical_opt_{method}.log', False)
psi4.set_options({'reference':'uhf','guess':'read', 'mp2_type':'conv', 'maxiter':200, 'soscf':True, 'soscf_max_iter':20})
E_mp2, wfn = psi4.optimize(f'{method}/6-31+G*', molecule=radical, return_wfn=True)

E_mp2

6.2.1.3. Geometrical parameters#

Let’s have a look at the geometrical parameters predicted by MP2 calculation. The following cell uses the calculate_bond and calculate_angle helper functions for that. The atom numbering can be checked with another helper function, drawXYZ_labeled.

# Convert molecule coordinates to numpy array
r_numpy = radical.to_arrays()[0]*hartree2A 

print('MP2 values')
print(f'r1: {calculate_bond(r_numpy[0], r_numpy[1]):.3f}')
print(f'r2: {calculate_bond(r_numpy[0], r_numpy[3]):.3f}')
print(f'r3: {calculate_bond(r_numpy[0], r_numpy[2]):.3f}')
print(f'O-N-O angle: {calculate_angle(r_numpy[1], r_numpy[0], r_numpy[2]):.3f}')
print(f'O-N-O angle: {calculate_angle(r_numpy[3], r_numpy[0], r_numpy[2]):.3f}')

drawXYZ_labeled(radical)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

6.2.1.4. Visualizing the orbitals#

In order to assess the reliability of an electronic structure calculation, it is usually not sufficient to simply glance at the optimised structure and energies. One should as well examine the orbitals. This is also valid for DFT: Although the Kohn-Sham orbitals lack, in a strict sense, any physical interpretation, experience shows that they are very close to wavefunction-based single particle orbitals. Thus, the DFT Kohn-Sham orbitals provide a useful interpretative tool to visualise changes in the electronic structure between different species (or in a chemical transformation), but one should bear in mind that they have to be interpreted with care.

In Psi4 we will use the frontier_orbitals setting, which produces cube representations of the frontier molecular orbitals. For closed shell species, the highest occupied (HOMO) and the lowest unoccupied (LUMO) alpha orbitals (ie. \(\psi_\alpha(r)\)) are printed, while for open shell species a total of (\(4+M_s\)) orbitals are printed (\(\alpha\) and \(\beta\) spin for both lowest virtual (LVMO) and highest doubly occupied orbitals (DOMO), along with all singly occupied (SOMO) orbitals). The plot below will print \(\alpha\) spin orbitals in red and \(\beta\) spin orbitals in blue.

You can visualize the cube files using the commands findCubeFiles(pathToDirectory) and showOrbitals(pathToDirectory)

Using cubeprop you can also access the electron density (density), the basis functions (basis_functions), the electrostatic potential (esp) and a descriptor calculated from the frontier orbitals which is \(f^2(r) = \rho_{\text{LUMO}}(r) - \rho_{\text{HOMO}}(r)\) (dual_descriptor). This descriptor is a good measure of nucleophilicity and electrophilicity and is like both Fukui functions combined.

psi4.set_options({'cubeprop_tasks':['frontier_orbitals']})
psi4.cubeprop(wfn)

orbitals = findCubeFiles('')
show_orbitals(mol=radical)

6.2.1.5. DFT#

Run the same geometry optimization for two of the following density functionals: BLYP, BP86, PBE, B3LYP, B97-2, M06-L, M06-2X, TPSS, mPW1PW91, wb97x-d.

For this small system, Hartree-Fock and MP2 were surprisingly fast; however, you may rest assured that for larger molecules and basis sets, all the DFT methods will outperform MP2 in computational efficiency, with some of them even beating Hartree-Fock

psi4.core.clean_options()
psi4.core.set_output_file(f'radical_opt_{method}.log', False)
psi4.set_options({'reference':''}) # Choose a reference here, for DFT it can be Unrestricted KS (uks), restricted (rks) or restricted open shell (roks)
E_dft, wfn = psi4.optimize(f'method/6-31+G*', molecule=radical, return_wfn=True)

# Convert molecule coordinates to numpy array
r_numpy = radical.to_arrays()[0]*hartree2A 

print('DFT: Method') #insert your method here (used just for printing)
print(f'r1: {calculate_bond(r_numpy[0], r_numpy[1]):.3f}')
print(f'r2: {calculate_bond(r_numpy[0], r_numpy[3]):.3f}')
print(f'r3: {calculate_bond(r_numpy[0], r_numpy[2]):.3f}')
print(f'O-N-O angle: {calculate_angle(r_numpy[1], r_numpy[0], r_numpy[2]):.3f}')
print(f'O-N-O angle: {calculate_angle(r_numpy[3], r_numpy[0], r_numpy[2]):.3f}')

psi4.set_options({'cubeprop_tasks':['frontier_orbitals']})
psi4.cubeprop(wfn)

orbitals = findCubeFiles()
show_orbitals(mol=radical)

6.2.2. Results comparison#

Now that you have performed the required calculations for HF, MP2 and (some) DFT case, let’s compare the performances of the different methods You can find a table with the experimental values for the different geometrical parameters below

Method:	\(\phi_1\) [\(^\circ\)]	\(\phi_2\) [\(^\circ\)]	\(\mathbf{r}_1 (O-N)\) [Å]	\(\mathbf{r}_2 (O-N)\) [Å]	\(\mathbf{r}_3 (O-N)\) [Å]	Symmetry
Exp.	120	120	1.24	1.24	1.24	D\(_3^h\)

Exercise 5

Fill the above table with the results of your calculations. Compare the performance of the methods in predicting accurate
structures.     Is there a trend relating the complexity of the exchange-correlation functional (LDA, GGA, hybrid, ...) to the quality of your results?

**Bonus**: For the DFT methods, specify what approximations are
used in the exchange functional (you may refer to {ref}`dfttheory`). 

Exercise 6

Take a screenshot of the MP2 and DFT SOMOs. How do the orbitals differ from each other? How does the symmetry of the SOMO relate to the predicted structure?

Exercise 7

Explain the difference between static and dynamic correlation. Relate this to the different results you obtained for the nitrate radical. What kind of correlation do HF, MP2 and DFT take into account?

Exercise 8

If you needed a highly accurate structure and energy for NO\(_3\cdot\), which method would you use? Why?