G. 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:
# . Define the molecule and its energy model.
molecule = ImportSystem ( os.path.join ( molPath, "cyclohexane_chair.mol" ) )
molecule.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )
molecule.DefineNBModel ( NBModelFull.WithDefaults ( ) )
molecule.Summary ( )
# . Determine the starting energy.
eStart = molecule.Energy ( )
# . Optimization.
BakerSaddleOptimize_SystemGeometry ( molecule ,
logFrequency = 1 ,
maximumIterations = 250 ,
rmsGradientTolerance = 1.0e-3 )
# . Determine the final energy.
eStop = molecule.Energy ( )
# . Print the energy change.
logFile.Paragraph ( "Energy change after search = {:20.4f}".format ( eStop - eStart ) )
# . Save the coordinates.
molecule.label = "Cyclohexane saddle conformation - OPLS-AA optimized."
ExportSystem ( os.path.join ( scratchPath, "cyclohexane_saddle.xyz" ), molecule )
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:
# . Define the molecule and its energy model.
molecule = ImportSystem ( os.path.join ( molPath , "cyclohexane_chair.mol" ) )
molecule.coordinates3 = ImportCoordinates3 ( os.path.join ( scratchPath, "cyclohexane_saddle.xyz" ) )
molecule.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )
molecule.DefineNBModel ( NBModelFull.WithDefaults ( ) )
molecule.Summary ( )
# . Calculate an energy.
molecule.Energy ( )
# . Optimization.
trajectory = ExportTrajectory ( os.path.join ( scratchPath, "cyclohexane_sdpath.ptGeo" ), molecule )
SteepestDescentPath_SystemGeometry ( molecule ,
functionStep = 2.0 ,
logFrequency = 10 ,
maximumIterations = 400 ,
pathStep = 0.025 ,
saveFrequency = 10 ,
trajectory = trajectory ,
useMassWeighting = True )
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 these examples. The program above makes use of a trajectory in pDynamo's geometry trajectory format, ptGeo. However, it also has other trajectory formats that may be preferable in some circumstances, such as when wanting to visualize trajectory data with third-party programs (e.g. VMD). See the wiki for further details.
An alternative to the piecemeal approach described above is to determine a reaction path in one fell swoop. Example13.py employs a chain-of-states (CoS) algorithm. The program works by defining reactant (chair) and product (twist-boat) structures and generating an initial 11-structure guess for the reaction path using the growing string method. The trajectory containing the guess is then optimized by the CoS algorithm:
# . Define the molecule and its energy model.
molecule = ImportSystem ( os.path.join ( molPath, "cyclohexane_chair.mol" ) )
molecule.DefineMMModel ( MMModelOPLS.WithParameterSet ( "bookSmallExamples" ) )
molecule.DefineNBModel ( NBModelFull.WithDefaults ( ) )
molecule.Summary ( )
# . Assign the reactant and product coordinates.
reactants = ImportCoordinates3 ( os.path.join ( xyzPath, "cyclohexane_chair.xyz" ) )
products = ImportCoordinates3 ( os.path.join ( xyzPath, "cyclohexane_twistboat.xyz" ) )
# . Get an initial path.
trajectoryPath = os.path.join ( scratchPath, "cyclohexane_sawpath.ptGeo" )
GrowingStringInitialPath ( molecule, 11, reactants, products, trajectoryPath )
# . Optimization.
trajectory = ExportTrajectory ( trajectoryPath, molecule, append = True )
ChainOfStatesOptimizePath_SystemGeometry ( molecule ,
trajectory ,
logFrequency = 1 ,
maximumIterations = 100 ,
rmsGradientTolerance = 0.1 )
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?