H. Normal Modes

Normal Mode Analysis

Normal mode analysis is a technique for obtaining approximate dynamical and thermodynamical information if only single, usually stationary-point, structures of a system are available. It makes the assumption that the potential energy surface about a point is harmonic (i.e. all derivatives higher than the second are zero) which leads to a dynamical model that can be solved exactly. An essential part of a normal mode calculation is the evaluation of the matrix of second derivatives of the energy. This is expensive for large systems, and also for QC methods, which limits its application in these cases.

The program Example14.py performs a normal mode analysis for the chair conformation of cyclohexane and then creates a trajectory of structures which reproduces the motion produced by the first mode of non-zero frequency at 600 K. Such trajectories are useful for structural analyses or for inspecting mode motions with third-party visualization programs.

# . Define the molecule and its energy model.

molecule = ImportSystem ( os.path.join ( molPath, "cyclohexane_chair.mol" ) )

molecule.coordinates3 = ImportCoordinates3 ( os.path.join ( xyzPath, "cyclohexane_chair.xyz" ) )

molecule.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )

molecule.DefineNBModel ( NBModelFull.WithDefaults ( ) )

molecule.Summary ( )

# . Calculate the normal modes.

NormalModes_SystemGeometry ( molecule, modify = ModifyOption.Project )

# . Generate a trajectory for one of the modes.

trajectory = ExportTrajectory ( os.path.join ( scratchPath, "cyclohexane_chair_mode7.ptGeo" ), molecule )

NormalModesTrajectory_SystemGeometry ( molecule ,

trajectory ,

mode = 7 ,

cycles = 10 ,

frames = 21 ,

temperature = 600.0 )

The vibrational frequencies, and mode vectors, coming from a normal mode analysis can be employed in a number of applications. One of these is as input for the calculation of various thermodynamic and kinetic parameters using formulae from statistical thermodynamics. The following program, Example15.py, calculates the Gibb's free energy differences and equilibrium constants between the chair and twist-boat forms of cyclohexane within the rigid-rotor harmonic-oscillator approximation at various temperatures:

# . Methods.

def FreeEnergies ( fileName, temperatures, symmetryNumber = 1 ):

"""Calculate the potential energy for a system and its

Gibbs free energies at several temperatures."""

# . Define the molecule and its energy model.

molecule = ImportSystem ( os.path.join ( molPath, fileName + ".mol" ) )

molecule.coordinates3 = ImportCoordinates3 ( os.path.join ( xyzPath, fileName + ".xyz" ) )

molecule.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )

molecule.DefineNBModel ( NBModelFull.WithDefaults ( ) )

molecule.Summary ( )

# . Calculate the energy and normal modes.

e = molecule.Energy ( )

NormalModes_SystemGeometry ( molecule, modify = ModifyOption.Project )

# . Loop over the temperatures.

g = []

for T in temperatures:

tdics = ThermodynamicsRRHO_SystemGeometry ( molecule ,

pressure = 1.0 ,

symmetryNumber = symmetryNumber ,

temperature = T )

g.append ( tdics["Gibbs Free Energy"] )

# . Return the energies.

return ( e, g )

# . Header.

logFile.Header ( )

# . Create a sequence of temperatures.

temperatures = [ 100.0 * i for i in range ( 1, 11 ) ]

# . Get the energies for the boat and chair structures.

( eB, gBoat ) = FreeEnergies ( "cyclohexane_twistboat", temperatures, symmetryNumber = 4 )

( eC, gChair ) = FreeEnergies ( "cyclohexane_chair", temperatures, symmetryNumber = 6 )

deltaE = ( eB - eC )

# . Output the equilibrium constants.

table = logFile.GetTable ( columns = [ 25, 25 ] )

table.Start ( )

table.Title ( "Equilibrium Constants (Chair -> Twist Boat)" )

table.Heading ( "Temperature" )

table.Heading ( "Log K" )

for ( T, gC, gB ) in zip ( temperatures, gChair, gBoat ):

RT = ( Constants.Molar_Gas * T ) / 1000.0

log10K = math.log10 ( math.e ) * ( - ( gB - gC + deltaE ) / RT )

table.Entry ( "{:.4f}".format ( T ) )

table.Entry ( "{:.6g}".format ( log10K ) )

table.Stop ( )

This program is more complicated than many others in the tutorial as it contains a function, FreeEnergies, that has been written to calculate the potential energy of a species along with its Gibb's free energies at a variety of temperatures.

Exercises

• 1. Modify Example14.py to examine different modes of the chair conformation of cyclohexane and also different conformations. What does the saddle point mode of imaginary frequency look like?

  2. Calculate the rates for the chair ↔ twist-boat isomerization using a program based on Example15.py. For the rates, use a simple rigid-rotor, harmonic-oscillator transition state theory formula with a constant frequency prefactor.