When you minimized systems in simple quadratic or Morse potential, you could analytically differentiate the energy expression twice to obtain the Hessian. If you were really attentive, you might have figured out that you even could differentiate the variational energy expression for the hydrogen atom's Schrodinger equation twice and minimize it using the Newton-Rhapson method. But even in case of diatomic molecules, we do not have a simple analytical formula that relates the quantum mechanical energy to the distance between nuclei. However, the potential can be calculated point-by-point by solving the electronic Schrodinger equation for several internuclear distances. In some cases, the complex relationship between energy and molecular geometry can be analytically differentiated; computational chemists then say that analytic gradients and the Hessian are available. For example, programs suchas Gaussian and PC GAMESS / Firefly can calculate analytic first and second derivatives at the Hartree-Fock level. Such analytical derivatives are not always possible, for example PC GAMESS cannot calculate analytic derivatives at the MP2 level; Gaussian03 can calculate analytical MP2 second derivatives but not analytical CCSD(T) second derivatives. When analytic derivatives are not avalable, numeric differentiation of energies usng finite difference approximatios allows to calculate derivatives numerically. This is trivial for diatomics and practical with small molecules but becomes inefficient as the number of geometric variables increases.
We know from elementary geometry how to calculate the slope of line based on the coordinates of the two points that define the line: rise over run. The slope is the first derivative. It turns out that the logic that leads to the "rise over run" formula can be generalized to calculate higher derivatives. The method of obtaining derivatives of the function from the function values is called the finite difference method. To calculate the second derivative, minimally three points are needed, to calculate the third derivative minimally four points are neeted and so on. In practice, better accuracy is obtained by using more points; for example accurate first and second derivative at the minimum energy geometry can be calculated by evaluating the energy at five-point stencil about the energy minimum.
The computation of derivatives via simple finite difference methods requires a set of evenly spaced energy values. The first derivative should be zero at the minimum; the second derivative determines the infrared wavenumber. The calculation of infrared frequency is straighforward after some combersome unit conversion (from Hartrees and Angstroms to cm-1. Here is an example code:
#!/usr/bin/env python
# Program to perform numeric vibrational analysis on F2:
from numpy import *
# Universal Constants and Conversion Factors
Hartree_to_J = 4.35974394e-18 # 2*R_inf * h * c; CODATA 2006
Ang_per_Meter = 1e+10
hp = 6.62606896e-34 # Planck constant in J*s, CODATA 2006
c = 29979245800 # Speed of light in cm / s
amu = 1.660538782e-27 # Atomic mass unit, CODATA 2006
# User should enter the minimum energy distance here
r_min = 1.40
# Atomic masses
m1 = 19.992440
m2 = 19.992440
mu = m1 * m2 / (m1 + m2)
print "Reduced mass in a.m.u. is:", mu
red_mass = mu * amu # in kg
# Return the first and second derivatives via O(h4) finite central difference
def Der12(r):
h = 0.01
r_p1 = r + 1*h
r_p2 = r + 2*h
r_m1 = r - 1*h
r_m2 = r - 2*h
# User could substitute energy values below from their QM calculations:
H_oo = Hamiltonian(r)
H_p1 = Hamiltonian(r_p1)
H_p2 = Hamiltonian(r_p2)
H_m1 = Hamiltonian(r_m1)
H_m2 = Hamiltonian(r_m2)
force_r = ( H_m2 - 8*H_m1 + 8*H_p1 - H_p2) / (12*h)
hessian_r = ( -H_m2 + 16*H_m1 - 30*H_oo + 16*H_p1 - H_p2) / (12*h*h)
return H_oo, force_r, hessian_r
def VibAnal(r):
d = Der12(r) # array of (energy,force,Hessian)
E = d[0] # Energy
F = d[1] # Force (first derivative)
H = d[2] # Hessian (second derivative)
print "Step Distance Energy Force Hessian "
print "0", r, E, F, H # Values
if (abs(F)) < 0.001:
freq = sqrt(H * Hartree_to_J * Ang_per_Meter * Ang_per_Meter / red_mass ) / (2 * pi * c)
print "**************************************************************"
print "* NUMERIC FREQUENCY ANALYSIS "
print "* STATIONARY POINT DISTANCE R =", r," "
print "* HARMONIC FREQUENCY =", freq, "cm-1"
print "**************************************************************"
else:
print "Force too large:", F
VibAnal(r_min)
You can use this code when answering your homework questions about fluorine.