Molecular Dynamics Simulations

Molecular dynamics simulations provide a powerful means of investigating the properties of molecular systems. They have two principal roles: (i) to study how a system changes as a function of time; and (ii) to explore a much wider range of conformational space that is accessible to a system than the geometry optimization and reaction path methods that we have discussed up to now.

Example16.py of the pDynamo distribution performs a simple dynamics simulation of the bALA molecule in vacuum using a velocity Verlet algorithm. The simulation consists of three phases: (i) a heating phase in which the temperature is increased gradually from 10 to 300 K; (ii) an equilibration phase in which the temperature is kept at around 300 K; and (iii) a data-collection phase in which no temperature manipulation is effected and during which structural data is saved at 100 fs intervals on an external trajectory. The program is:

        # . Define the energy models.
        mmmodel = MMModelOPLS ( "bookSmallExamples" )
        nbmodel = NBModelFull ( )

        # . Generate the molecule.
        molecule = MOLFile_ToSystem ( os.path.join ( molpath, "bala_c7eq.mol" ) )
        molecule.DefineMMModel ( mmmodel )
        molecule.DefineNBModel ( nbmodel )
        molecule.Summary ( )
        molecule.Energy  ( )

        # . Optimization.
        ConjugateGradientMinimize_SystemGeometry ( molecule,                    \
                                                   maximumIterations    = 2000, \
                                                   logFrequency         =  100, \
                                                   rmsGradientTolerance =  0.1  )

        # . Define a random number generator in a given state.
        rng = Random ( )
        rng.seed ( 175189 )

        # . Heating.
        VelocityVerletDynamics_SystemGeometry ( molecule,                             \
                                                logFrequency              =      100, \
                                                rng                       =      rng, \
                                                steps                     =     1000, \
                                                timeStep                  =    0.001, \
                                                temperatureScaleFrequency =      100, \
                                                temperatureScaleOption    = "linear", \
                                                temperatureStart          =     10.0, \
                                                temperatureStop           =    300.0  )

        # . Equilibration.
        VelocityVerletDynamics_SystemGeometry ( molecule,                               \
                                                logFrequency              =        500, \
                                                steps                     =       5000, \
                                                timeStep                  =      0.001, \
                                                temperatureScaleFrequency =        100, \
                                                temperatureScaleOption    = "constant", \
                                                temperatureStart          =      300.0  )

        # . Data-collection.
        trajectory = SystemGeometryTrajectory ( os.path.join ( scratchpath, "bala_c7eq.trj" ), molecule, mode = "w" )
        VelocityVerletDynamics_SystemGeometry ( molecule,             \
                                                logFrequency =   500, \
                                                steps        = 10000, \
                                                timeStep     = 0.001, \
                                                trajectories = [ ( trajectory, 100 ) ] )

To be of use, of course, the trajectory data generated during a dynamics has to be analyzed. Example17.py shows how to do a simple analysis by reading the frames in the trajectory and, for each frame, calculating the φ and ψ dihedral angles of the molecule. These values are printed out to a table at the end of the program along with their means and standard deviations.

        # . Read the molecule definition.
        molecule = MOLFile_ToSystem ( os.path.join ( molpath, "bala_c7eq.mol" ) )
        molecule.Summary ( )

        # . Define the trajectory.
        trajectory = SystemGeometryTrajectory ( os.path.join ( scratchpath, "bala_c7eq.trj" ), molecule, mode = "r" )

        # . Loop over the frames in the trajectory.
        phi = []
        psi = []
        while trajectory.RestoreOwnerData ( ):
            phi.append ( molecule.coordinates3.Dihedral ( 4, 6,  8, 14 ) )
            psi.append ( molecule.coordinates3.Dihedral ( 6, 8, 14, 16 ) )

        # . Set up the statistics calculation.
        phistatistics = Statistics ( phi )
        psistatistics = Statistics ( psi )

        # . Output the results.
        table = logfile.GetTable ( columns = [ 20, 20, 20 ] )
        table.Start   ( )
        table.Title   ( "Phi/Psi Angles" )
        table.Heading ( "Frame" )
        table.Heading ( "Phi"   )
        table.Heading ( "Psi"   )
        for ( i, ( h, s ) ) in enumerate ( zip ( phi, psi ) ):
            table.Entry ( `i` )
            table.Entry ( "%.2f" % ( h, ) )
            table.Entry ( "%.2f" % ( s, ) )
        table.Entry ( "Mean:",               alignment = "l" )
        table.Entry ( "%.2f" % ( phistatistics.mean, ) )
        table.Entry ( "%.2f" % ( psistatistics.mean, ) )
        table.Entry ( "Standard Deviation:", alignment = "l" )
        table.Entry ( "%.2f" % ( phistatistics.standardDeviation, ) )
        table.Entry ( "%.2f" % ( psistatistics.standardDeviation, ) )
        table.Stop ( )

Exercises

pDynamo has a number of dynamics algorithms. Try some of these (e.g. LangevinDynamics_SystemGeometry) and see how their behavior compares to that of the velocity Verlet method by analyzing the simulation data in different ways.
Choosing a suitable dynamics algorithm, perform extended simulations of the bALA molecule — either a very long simulation or many small ones — and then plot the φ/ψ map of the simulation results. Is there any correspondence between the dynamic map and those determined previously using geometry optimization and reaction path calculations? If not, why not? What happens when the temperature of the simulations is changed?