Vibrational Quantities - Example of Nitrogen molecule

Back to Vibrational Relationships - Derivation

A nitrogen molecule, N₂, provides a good, simple, quantitative introduction to molecular vibration quantities.  The MNDO method will be used mainly because it is very old and stable, and all the steps described here can be reproduced using any copy of MOPAC. Several different quantities will be used in this example.  To assist in conversions, a set of conversion factors are useful.

The starting point will be the N₂ molecule at its equilibrium distance, 1.1038 Ångstroms.

Data-set

The following dataset was used as the source for many of the numbers mentioned here.  Keywords DFORCE, FMAT, LARGE, NOSYM, and NOREOR were used in order to provide expanded output, normally only FORCE would be used.

DFORCE MNDO FMAT LARGE NOSYM NOREOR
  Calculate the vibrational frequency of N2
  
 N     0.00000000 +0   0.00000000 +0   0.00000000 +0
 N     1.10380157 +1   0.00000000 +0   0.00000000 +0

Build the Hessian matrix

Hessian matrices for vibrational calculations are constructed using the gradients of the displaced Cartesian coordinates for each atom, and the following expression, for a description of how this is constructed, see Hessian matrix.

In this equation, Δx_i = Δx_j = 0.00833333 Å. This step-size was selected automatically, and similar step-size would work just as well. Only one matrix element, H(1,4), is important, this represents the N-N stretch along the "X" axis.

Four gradients are present in this equation, their values for Hessian element H(1,4) are:

Ideally, gradients 1 and 4 should be the same, as should gradients 2 and 3, but because of numerical imprecision there are small differences.  A 30% reduction in numerical imprecision can be achieved by averaging gradients 1 and 4, and 2 and 3.

Matrix element H(1,4) can then be calculated using these gradients, giving H(1,4) = -4455.175 kcal x mol^-1 x Å^-2

Convert to SI

Converting the Hessian from (kcal mol^-1 Å^-2) to (millidynes Å^-1) involves multiplying the Hessian by a factor.  This factor can be expressed variously as:

When this is done, the matrix element becomes H(1,4) = -4455.175 x 6.9477 x 10^-3 = -30.9520 millidynes x Å^-1

Mass-weighting the Hessian

The final step in preparing the matrix for diagonalization involves mass-weighting.  This operation is straightforward:

For all matrix elements:

Nitrogen atoms have a mass of M_i = 14.0067 amu, so matrix element H(1,4) is:

Convert the mass-weighted Hessian into normal modes of vibration

When the mass-weighted Hessian is diagonalized, the result is a set of eigenvalues, ε_i, and associate eigenvectors, ψ_i.  The eigenvectors are the coefficients of the wave-function that represents the normal modes of vibration, and form an orthonormal set.  The eigenvalues represent the force constants for the modes divided by the effective mass.  Nitrogen has only one non-trivial vibration, so:

and the eigenvector: ψ₁ is (0.7071) x (N₁(x) - N₂(x)).  This represents the N - N stretch mode.

Recap of force constant and reduced mass

Eigenvalues, ε_i, of the mass-weighted Hessian are related to the force-constant, k_i, and reduced mass, μ_i, for a normal mode of vibration by:

ε_i = k_i/μ_i 

Calculation of the Reduced Mass for the N-N stretch

Given the normal mode eigenvector, the reduced mass for the vibration can be calculated from:

μ_i = Σ_A((ψ_i(A)(x))² + (ψ_i(A)(y))² + (ψ_i(A)(z))²)² M_A amu

where "A" refers to atom A, and "x", "y", and "z" refer to the respective Cartesian coordinate of atom "A".

For the vibration in N₂, this sum is:

Calculation of the Force Constants for the N-N stretch

Given ε₁ and μ₁, the force constant, k₁ is:

Convert eigenvalues to wavenumbers (frequencies) in reciprocal centimeters

The eigenvalues are then converted into vibrational frequencies, ν_i, in cm^-1 via:

Here, "c" is in centimeters per second, N is in moles, and ε is in millidynes x Å^-1 x amu^-1

Convert wavenumbers to period of oscillation (T) in femtoseconds

The period for a molecular vibration is normally given in femtoseconds, fs; one femtosecond is 10^-15 of a second. The conversion from wavenumber to period is given by:

Convert wavenumbers to Angular Frequency, ω

Angular frequency is the angular displacement per second, or 2π/T:

Convert frequency into kcal/mol

Zero point energy

The zero-point energy of a polyatomic molecule is the sum of the zero-point energies of all its normal modes.  The zero-point energy of a normal mode is half of its vibrational energy.  For N₂, the vibration of 2738.8 cm^-1, corresponds to an energy of 7.83 kcal mol^-1, so its zero-point energy is 3.915 kcal mol^-1.

Convert energy represented by frequency into velocity of the atoms

The starting point for this operation is the familiar E = ½mv². For this operation, the cgs system will be used.  Rearranging E = ½mv² gives the velocity, v = (2E/m)^½. At this point, E is the energy in ergs per mole, and "m" is the reduced mass of the vibration in amu. These quantities can be re-worded as ergs per molecule and reduced mass in grams. Expanding and solving for the velocity gives:

This velocity can be split into velocity per atom, v_A, using the vibrational wavefunction, ψ_i.  First, the wavefunction must be re-normalized so that the total motion is 1.0 Ångstroms.  This is done by dividing the coefficients for each atom by the square-root of the atomic mass:

 ψ_iA^' = ψ_iA x m_A^-½

Then re-normalizing using a normalization factor N,

N = (Σ_A (ψ_iA^'_(x)²  + ψ_iA^'_(y)²  + ψ_iA^'_(z)²   )^½)^-1

to give:

ψ_i^" = N x ψ_i^'

Finally, the starting velocity of each atom can be calculated from:

 v_A = v x ψ_i^"(A)

For nitrogen, the wavefunction from the mass-weighted Hessian was ψ₁ = 0.7071 x N₁(x) - 0.7071 x N₂(x).  After re-normalizing, this becomes ψ₁^" = 0.5 x N₁(x) - 0.5 x N₂(x), and the velocity of each atom becomes:

v₁ = 1.52935 x 10⁵ cm x sec^-1

and

v₂ = -1.52935 x 10⁵ cm x sec^-1

in the "x" direction.  This calculation is a bit tedious, but an easier way to obtain the same result is to examine the output file from a DRC calculation, see below.

 Travel

In a molecular vibration the system starts at the energy minimum, then each atom moves along a path defined by the normal mode eigenvector.  At the turning point - the point where atomic motion stops - the sum of the distances that the atoms are from their original position is the travel distance.  In the case of the N₂ molecule, this is the distance from N₁ to its equilibrium position plus the distance of N₂ from its equilibrium position.  This quantity can also be calculated using the energy of the vibration and the force constant for the mode.  The cgi system is used in this calculation:

The travel can be partitioned into contributions from individual atoms using the same strategy as that used in the previous section on velocity, but a faster way to obtain the same result is to examine the output file from a DRC calculation, see below..

Maximum acceleration

Acceleration is at a maximum at the turning point of the energy in a molecular vibration.  At that point, the geometry is at its maximum distortion from the equilibrium.  Assuming simple harmonic motion, then the acceleration, a, is given by:

a = -(2πf)²x

where f = frequency in hertz. = 1/(12.17 *10^-15) = 8.22 x 10¹³ hertz, and x, the travel or distance from the geometry at the energy minimum, is a = 0.05929 Å = 5.929 x 10^-10 cm.  Solving this gives:

a = 2.667 x 10²³ x 5.929 x 10^-10 cm x s^-2  = 1.581 x 10¹⁴ cm x s^-2 

This acceleration can be compared with the acceleration due to gravity.  1G = 9.807 x 10² cm x s^-2, therefore the maximum acceleration is 1.612 x 10¹¹ G

Animation of a single vibrational cycle for N₂

An animation of a single complete vibration of the N-N stretch normal mode and the output can be used in verifying many of the calculations shown here.  

To reproduce this cycle, use the following data-set:

IRC=1 DRC T-priority large bigcycles=1  mndo HTML 
N2 one complete vibrational cycle
(To see N-N distance, click on "Script" at the bottom.)

N
N 1.10380157 1 

The step-size used in this animation is 0.1 femtosecond or 1 x10^-16 s.
Small differences can be seen when the calculated results and the equivalent results obtained from the animation or output are compared.  Some of these are due to round-off, and some are due to the difference between simple harmonic motion and the calculated trajectory, but in general the agreement is very good, indicating that the trajectory was almost simple harmonic.