Reaction Paths

An important problem in contemporary computational molecular
science is the characterization of pathways between the stable states of
a system. In what follows we consider two different techniques to find
pathways for the chair ↔ twist-boat isomerization of
cyclohexane using an OPLS-AA MM energy model. The stable chair
and twist-boat conformations of cyclohexane have previously been
determined by geometry optimization.

The first approach employs the two example programs, Example11.py and Example12.py.
The former searches for a first-order saddle point starting from the
chair conformation of cyclohexane, whereas the latter builds
the reaction path that passes through the saddle-point
structure. Example11.py is:

The saddle-point-location algorithm that is invoked in this program is of the mode-following
type and requires either the second (coordinate) derivatives of the
system's potential energy or an approximation to them. This limits the
application of this class of
algorithm to relatively small systems. Once the saddle point has
been located, its structure is saved in an XYZ file for use in the
subsequent program.

Example12.py takes the previously
determined saddle point structure and maps out the reaction path that
descends from it, in both directions, using a simple steepest-descent
technique. Due to the naivety of the algorithm, the path trajectory
does not terminate at the minimum-energy chair and twist-boat
conformations but only in their vicinity. In full the program is:

An important aspect of this program is the use of a trajectory
object to store the path structures that the steepest-descent technique
produces. Trajectories are ubiquitous in pDynamo (and, indeed, other
molecular modeling programs) and are employed by a great
many of its algorithms to store structural and other
information. The analysis of trajectory data is an important topic in
its own right and will be treated later in this tutorial. The program
above makes use of a trajectory of the SystemGeometryTrajectory
class. However, pDynamo has other classes of trajectory that may be
preferable in some circumstances, such as when wanting to visualize
trajectory data with third-party programs (like VMD). See the Trajectories tutorial for further details.

An alternative to the piecemeal approach described above is to determine a reaction path in one fell swoop. Example13.py employs the self-avoiding walk (SAW)
algorithm. This is a relatively old method but newer, more-efficient
ones
have been implemented in pDynamo and will be described in other
tutorials. The program works by defining reactant (chair) and product
(twist-boat) structures and generating an initial 11-structure guess for
the reaction path by linear interpolation. The trajectory
containing the guess is then optimized by the SAW algorithm:

Exercises

Employ the methods described above to map out the reaction paths
between the various conformations of the bALA molecule. How do the
results compare to those obtained by constrained φ/ψ geometry
optimization?