5.4. Influence of correlation on geometry: HNO\(_3\) molecule#
The influence of electron correlation on structural parameters such as bond lengths and angles is not negligible. Hartree-Fock calculations may give bond lengths that deviate considerably from the experimental value. Comparably cheap methods such as MPn may often remedy sufficiently for this problem, but some smaller deviations may remain. The most rigorous approach in such cases would be e.g. a full configuration interaction (full CI/FCI), but this is computationally untractable for any larger molecule using a reasonable basis. One may therefore restrict the method to singly and doubly excited Slater determinants, assuming that the main effects of correlation are captured well enough in such a truncation. However, such truncations introduce size-consistency problems.
In this exercise, we will assess and comment on the performance of HF, MP2 and MP3 with respect to the prediction of bond lengths and bond angles in HNO\(_3\) in reference to the experimental values (Phys. Chem. Chem. Phys., 2014 16, 19437).
Exercise 8
Assess and comment on the performance of the different methods with respect to experimental values. Can you think of possible reasons for differences between the experimental and computed values?
import psi4
import py3Dmol
import pandas as pd
import numpy as np
import sys
sys.path.append("..")
from helpers import *
psi4.core.clean_options()
hno3 = psi4.geometry("""
0 1
symmetry c1
O 1.326 -0.373 -0.127
N 0.015 -0.015 0.415
O -0.329 1.292 -0.145
O -0.944 -0.973 -0.135
H -1.811 -0.738 0.223""")
drawXYZ(hno3)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
psi4.set_options({'reference':'uhf','guess':'gwh'})
hno3_hf = hno3.clone()
psi4.core.set_output_file(f'hno3_hf.log', False)
E_hf, history_hf = psi4.optimize('scf/6-31+G*', molecule=hno3_hf, return_history=True)
Change in internal coordinate of 5.72e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 2.
Energy ratio indicates iffy step.
Intrafrag trust radius decreased to 0.25.
Change in internal coordinate of 6.13e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 1.
Change in internal coordinate of 5.12e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 2.
Change in internal coordinate of 7.09e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 1.
Change in internal coordinate of 5.98e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 2.
Change in internal coordinate of 7.95e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 1.
Change in internal coordinate of 5.85e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 6.
Change in internal coordinate of 6.91e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 5.
Change in internal coordinate of 8.18e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 4.
Change in internal coordinate of 8.34e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 3.
Change in internal coordinate of 1.18e+00 exceeds limit of 5.00e-01.
Skipping Hessian update for step 2.
Previous geometry is closer to target in internal coordinates, so using that one.
Best geometry has RMS(Delta(q)) = 4.41e-04
Change in internal coordinate of 5.04e-01 exceeds limit of 5.00e-01.
Skipping Hessian update for step 7.
Energy has increased in a minimization.
Energy ratio indicates iffy step.
Intrafrag trust radius decreased to 0.188.
Previous geometry is closer to target in internal coordinates, so using that one.
Best geometry has RMS(Delta(q)) = 2.61e-07
Unable to completely converge to displaced geometry.
RMS(dx): 3.156e-06 Max(dx): 8.093e-06 RMS(dq): 4.063e-06
Previous geometry is closer to target in internal coordinates, so using that one.
Best geometry has RMS(Delta(q)) = 4.06e-06
Unable to completely converge to displaced geometry.
RMS(dx): 3.007e-06 Max(dx): 6.698e-06 RMS(dq): 4.336e-06
Previous geometry is closer to target in internal coordinates, so using that one.
Best geometry has RMS(Delta(q)) = 4.34e-06
Optimizer: Optimization complete!
hno3_mp2 = hno3_hf.clone()
psi4.core.set_output_file(f'hno3_mp2.log', False)
E_mp2, history_mp2 = None, None # run MP2 calculation here, need to use return_history=True to return the history as output (used later)
hno3_mp3 = hno3_mp2.clone()
psi4.core.set_output_file(f'hno3_mp3.log', False)
E_mp3, history_mp3 = None, None # run MP3 calculation here, need to use return_history=True to return the history as output (used later)
#save last optimized coordinates in a dictionary
opt_coord = {
'hf':np.array(history_hf['coordinates'][-1]),
'mp2':np.array(history_mp2['coordinates'][-1]),
'mp3':np.array(history_mp3['coordinates'][-1])
}
You can calculate the distance between two points is the L2 norm (Euclidean distance) as
a = np.array([0,0,0])
b = np.array([1,1,1])
np.linalg.norm(a-b)
this has been implemented in the calculate_bond function contained in helpers.py imported above.
The calculation for the angle is also possible using numpy. You can use the following function:
def calculate_angle(a,b,c):
ba = a - b
bc = c - b
cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
angle = np.arccos(cosine_angle)
return np.degrees(angle)
drawXYZ_labeled(hno3)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
# we use the opt_coord dictionary to calculate the properties using python
# the central N has index 1, the oxygens have index 0,2,3
# the units are in bohr, we convert to Å
results = {'hf':{},'mp2':{},'mp3':{},'exp':{'r(O0-N) Å':1.198, 'r(O2-N) Å':1.213, 'r(O3-N) Å':1.410, 'φ(O-N-O) °':130.2}, }
for method in ['hf', 'mp2', 'mp3']:
xyz = opt_coord[method] # store the geometry array optimized with the current method in a temp variable
results[method]['r(O0-N) Å'] = calculate_bond(xyz[0], xyz[1]) * psi4.constants.bohr2angstroms
results[method]['r(O2-N) Å'] = calculate_bond(xyz[2], xyz[1]) * psi4.constants.bohr2angstroms
results[method]['r(O3-N) Å'] = calculate_bond(xyz[3], xyz[1]) * psi4.constants.bohr2angstroms
results[method]['φ(O-N-O) °'] = calculate_angle(xyz[0], xyz[1], xyz[2])
pd.DataFrame.from_dict(results)