Force Field Analysis

PGOPHER can do a traditional force field analysis, calculating vibrational frequencies, the l matrix, isotope shifts, rotational constants, centrifugal distortion constants and Coriolis coupling constants from force constants and geometry. The force constants and geometry can be derived by fitting to any combination of the calculated quantities. The calculations are similar to those done by the ASYM20/ASYM40 program of Hedberg and Mills (1993, 2000), though implemented in a more flexible way, and these papers give a good outline of the methodology and how the calculations  are performed. The VIBCA program (http://www.ifpan.edu.pl/~kisiel/vibr/vibr.htm) also performs similar calculations.

The force field can be expressed in terms of valence or symmetry coordinates. A typical analysis will include several isotopologues and might have a data file structured something like this:

Note the various single and multiple objects:

• A single Global Variables object which contains geometric variables and force constants common to all isotopologues.
  
• A z Matrix for each isotopologue, which specifies the geometry. The use of variables means this will be identical for all isotopologues,  except for the masses.
• An (optional) Internal Coordinates object for each isotopologue, which specifies the valence coordinates and force constants. The use of variables means this will be identical for all isotopologues.
• An (optional) Symmetry Coordinates object for each isotopologue, which specifies the combinations of valence coordinates to use. The use of variables means  this will only differ for isotopologues of different symmetry.

A Cartesian Coordinates object may also be present; the exact combination of objects used depends on the required calculation. Nuclear Co-ordinates and Vibrational Modes objects will also be present, automatically created and updated from the force field analysis objects listed above.

The current version of the program does not explicitly make use of symmetry for the force field analysis, but it is implicit in the way the force constants and geometry are specified. The current version of PGOPHER is restricted to harmonic terms, but the program is structured so that anhamonic terms can be added in future. The units are currently less flexible than for rotational calculations; the overall mixture units must be set to cm-1, the geometry specified in Angstroms and degrees and force constants (for stretches) in millidyne/Å = aJ/Ų (1 millidyne/Å = 1 aJ/Ų = 100 J/m² = 100 N/m). For a full discussion of the units see Hedberg and Mills (1993).

where λ is a diagonal matrix, and the diagonal elements are the force constants in the normal coordinates, Q. Note that the 5 or 6 lowest frequency modes will correspond to translations and rotations, and will be excluded from most subsequent calculations. The test for this is on λ; values smaller than SmallCoefficient of the largest are discarded. The potential energy thus has the required diagonal form:

The dimensionless normal coordinates are related to the normal coordinates by a scaling factor, γ_i^1/2:

q_i = γ_i^1/2 Q_i

If the scaling factor is chosen to be γ_i = 2πλ_i^1/2/h the total energy has a particularly simple form:

T + V = ½ Σ_i (P_i² + λ_i Q_i²) = ½ Σ_i (2π)^-1 h λ_i^½ (p_i² + q_i²) = ½ Σ_i hν_i(p_i² + q_i²)

which also makes q_i dimensionless as required and identifies λ_i^½ as 2πν_i. The transformation matrix between these dimensionless normal coordinates and Cartesian displacements, which we define as d, is then given by:

Δx = M^-½ l Q = M^-½ l γ^-½ q = d q

For a harmonic oscillator this the range of the zero point motion is ±1, so the elements of this matrix give the range of the zero point motion. The transformation matrix between dimensionless normal coordinates and internal coordinates, d^int, follows from this:

R = B Δx = B M^-½ l γ^-½ q = B d q = d^int q

The matrix U is specified in the Symmetry Coordinates object. The potential energy can also be expressed in terms of symmetry coordinates:

Note the use of F for the matrix of force constants expressed in symmetry coordinates; it is related to the force constant matrix, f, used above by:

Either f or F can be given but the calculation of normal modes from symmetry coordinates requires F. If non redundant coordinates are used, U will be square so the above equation can be inverted by working out U^-1. If redundant coordinates are used, they must be chosen such that UU^T = 1, in which case F = U f U^T.

When using symmetry coordinates a slightly different approach is used to find the normal modes. The (standard) G matrix is calculated from:

G = U B M^-1B^T U^T

and Choleski decomposition is used to factorise the G matrix such that:

G = C C^T

where C is a lower triangular matrix. The matrix C^T F C has the same eigenvalues,  λ, as the GF matrix, and the eigenvectors, V, can be used to calculate the transformation matrix, L, between the symmetry and normal coordinates:

S = C V Q = L Q

and the transformation matrix, d^sym, between the symmetry and dimensionless normal coordinates:

S = L Q = L γ^-½ q = d^sym q

Given L the following can be calculated:

(L^-1)^T = F L λ^-1
G^-1 = (L^-1)^T L^-1
A = M^-1B^TG^-1 = M^-1B^T(L^-1)^T L^-1

where A is the matrix that transforms symmetry coordinates to Cartesian coordinates:

Δx = A S

In practice the matrix A L = M^-1B^T(L^-1)^T is calculated as the l matrix can be calculated from:

l = M^½ A L