Resources‎ > ‎Tutorials‎ > ‎

Transition Path Sampling


Overview

Transition path sampling (TPS) is a method used to study rare events in computer simulations.  This tutorial demonstrates the application of TPS to study nucleation in the two-dimensional Ising model.  It is based around a simple Ising model simulation code that I cooked up for the purpose of demonstration (i.e., it is not optimized and I will not promise that it is correct).  

Background

TPS is a method to sample rare events in a way that does not require the presumption of reaction coordinates.  It is therefore distinct from rare-events techniques that require the relevant order parameters be identified a priori, such as metadynamics and umbrella sampling (see Tutorial 4).  Some relevant references are outlined below:
Prerequisites

This tutorial demonstrates the sampling of nucleation pathways in the context of the Ising model. I assume that the reader is already at least somewhat familiar with the methods described in the references linked above. The tutorial makes use of the following software:

Step 0: Ising model "Hello World"

I assume that the reader is more than familiar with the Ising model model, but if not, see, for example, http://en.wikipedia.org/wiki/Ising_model.  

Change directories into the step0 subfolder of tutorial4: 

tutorial start

$ git clone https://github.com/askeys/tutorials.git
$ cd tutorials/tutorial3

The enclosed code generates an Ising model trajectory over which the system nucleates a droplet and changes phase:

step0.cpp

#include <lmc/lmc.h>
using namespace lmc;

//set the side length of the square lattice
int L = 50;

int main(int argc, char** argv)
{
    //Initialize the lattice for the system to run on:
    LatticeSquare lattice(L*L);

    //Set the random number generator:
    RandomNumberGenerator48 rng;
    
    //Initialize the simulation:
    SimulationIsingModel ising(lattice, rng);
    
    //Set up a list of initial spins (all -1):
    std::vector<double> s0(L*L, -1.0);
    ising.setSpins(s0);

    //Set the inverse temperature less than Tc ~ 2.3
    ising.setBeta(1.0/2.0);

    //Set a positive external field that favors +1 spins
    ising.setField(1.5);
    
    //Write the trajectory to an xyz file
    VisualizerXYZ xyz("traj.xyz");
    
    //Write 500 frames spaced by 10 MC sweeps
    for (int i=0; i<500; i++) {
        xyz.visualize(ising);
        ising.run(10);
    }
}

To compile the code, link with the lmc library (installed during the prerequisites). Example:

Compiling and running step 0


$ c++ -O3 step0.cpp -llmc -o step0
$ ./step0 
$ vmd traj.xyz 


The terminal should launch a VMD trajectory that can be played as a movie:

Notice that the waiting time for the rare nucleation event greatly exceeds the time scale of the event itself. 

Step 1: Plugging the Ising Model Engine Into the TPS Library

This demo uses a transition path sampling library that is designed to use any arbitrary simulation code as an engine to run trajectories.  To use a simulation code as an engine, the library employs polymorphism and virtual functions (see, for example, this page for a quick introduction to these constructs).  In short, this means that the library supplies a class who's member functions can be overwritten to provide functionality for a particular type of simulation.  The prototype (.h file) for the plugin looks like this:

TpsSimulationPlugin.cpp

class TpsSimulationPlugin
{
    public:

        TpsSimulationPlugin();
        virtual ~TpsSimulationPlugin();

//Required:
        virtual void run(int nsteps=1);
        virtual void writeRestart(const char* filename);
        virtual void readRestart(const char* filename);
        virtual void freeRestart(const char* filename);
        virtual void copyRestart(const char* src_file, const char* dest_file);
        virtual double computeHamiltonian();
        virtual double getExpectedTemperature();
        
//Required for deterministic dynamics:
        virtual void invertVelocities();

//... (Optional functions, etc)

};

Functions with the "virtual" keyword can be overridden in a derived class by defining a derived member function with the same name and arguments.  The functions labeled "required" must be overwritten in a derived class and re-implemented using the desired simulation engine.  

Change directories to step1.  The prototype for our Ising model simulation plugin (TpsSimulationPluginLibLMC.h) looks like this:

TpsSimulationPluginLibLMC.h


#include <tps/tps.h>
#include <lmc/lmc.h>
#include <string>
#include <vector>
#include <map>

class TpsSimulationPluginLibLMC : public TpsSimulationPlugin
{
    public:    

        //re-implemented from TpsSimulationPlugin:
        TpsSimulationPluginLibLMC(lmc::Simulation&);
        void run(int);
        void writeRestart(const char*);
        void readRestart(const char*);
        void freeRestart(const char*);        
        void copyRestart(const char*, const char*);        
        double computeHamiltonian();
                
        //extra helper functions specific to this plugin:
        int getDimensions();
        lmc::Simulation& getSimulation();
    
    protected:
        lmc::Simulation& _sim;
        std::map< std::string, std::vector<double> > _saved_timeslices;                        
};

TpsSimulationPluginLibLMC& safeDowncastToLibLMCPlugin(TpsSimulationPlugin&);

Notice that all of the required virtual functions are implemented in the derived class.  The implementation for each function is given in the TpsSimulationPluginLibLMC.cpp file.  The basic functionalities of the simulation required by the TPS library can now be controlled through the plugin.  An example is given in step1.cpp:

step1.cpp

#include "TpsSimulationPluginLibLMC.h"
#include <lmc/lmc.h>
#include <cassert>

using namespace lmc;
using namespace std;

int L = 50;

int main(int argc, char** argv)
{
    //Initialize simulation engine:
    LatticeSquare lattice(L*L);
    RandomNumberGenerator48 rng;    
    std::vector<double> s0(L*L, -1.0);
    SimulationIsingModel ising(lattice, rng);
    ising.setSpins(s0);
    ising.setBeta(1.0/2.0);
    ising.setField(1.5);
    
    //Plug engine into TPS framework:
    TpsSimulationPluginLibLMC sim(ising);

    //Make sure the simulation is reading and writing restart data correctly:
    sim.writeRestart("start.dat");
    double ei = sim.computeHamiltonian();
    sim.run(1000);
    double em = sim.computeHamiltonian();
    sim.readRestart("start.dat");
    double ee = sim.computeHamiltonian();    
    //(Energies should be the same before and after restarting):
    assert(em != ei);
    assert(ei == -4*L*L + 1.5*L*L);
    assert(ee == ei);
    
    //Make sure simulation is working properly through the TPS interface:
    sim.run(10000);
    double e = 0.0;
    for (int i=0; i<100; i++) {
        sim.run(100);
        e += sim.computeHamiltonian();
    }
    e/=100;
    //The energy should belower now:
    assert(e < -4*L*L);
}

The demo is compiled by linking with the LMC library as well as the TPS library: 

step1

$ c++ -O3 TpsSimulationPluginLibLMC.cpp step1.cpp -llmc -ltps
$ ./a.out

The demo code simply checks for errors in restarting from a file and computing energies.  If everything is configured correctly, the test program should exit with no errors.

Step2: Sample the Transition Path Ensemble (TPE)

The first step in the TPS formalism is to sample the ensemble of reactive trajectories that connect the two stable basins A and B (see the references listed at the top of this page for additional details).  This is accomplished by creating an initial reactive trajectory and then performing shooting and shifting moves to reach an equilibrium distribution of trajectories.  An example path sampling setup for the Ising model is given below:

step2.cpp (main)

int main(int argc, char** argv)
{
    //Initialize simulation engine:
    LatticeSquare lattice(L*L);
    RandomNumberGenerator48 rng48;    
    std::vector<double> s0(L*L, -1.0);
    SimulationIsingModel ising(lattice, rng48);
    ising.setSpins(s0);
    ising.setBeta(1.0/2.0);
    ising.setField(1.5);
    
    //Plug engine into TPS framework:
    TpsSimulationPluginLibLMC sim(ising);
    
    //Declare the type of trajectory: (100 steps total separated by 1 MC cycle) 
    int trajectory_length = 100;
    int step_size = 1;
    TpsTrajectoryUniformStep traj(sim, trajectory_length, step_size);

    //A class to get the first trajectory (just wait until an event happens) 
    TpsInitializerBruteForce init;

    //A class to hold the ensemble of reactive trajectories (TPE) 
    TpsTrajectoryEnsemble tpe(traj, 0, false);

    //A class to compute the order parameter and define different basins
    TpsOrderParameterNucleusSize op;
    
    //Random number generator for TPS
    TpsRNG48 rng(1, 2, 3);    
    
    //Initialize the first trajectory:
    cerr << "Initializing...\n";
    MyTpsAlgorithm tps(tpe, rng, op, init);
    cerr << "Sampling...\n";
    
    //Do shooting moves and shifting moves:
    TpsTrialMoveShotForward fshot;
    TpsTrialMoveShotBackward bshot;    
    TpsTrialMoveShift shift(0, LTRAJ, STEP);
    
    tps.addTrialMove(fshot, 0.5);
    tps.addTrialMove(bshot, 0.5);
    tps.addTrialMove(shift, 1.0);
    
    //Do path samping:
    int number_of_trajectories = 100;
    for (int i=0; i<NTRAJ; i++) {
        cout << i << "\n";
        tps.doStep();
        //trajectories are created in order, starting with index 0;
        //we can erase the old ones to save memory
        tpe.eraseTrajectories(0, tpe.getLastTrajectory().getID()-2);
    }
}

A few aspects of this code snippet require additional discussion.  First, the order parameter that defines the stable basins depends heavily on the system and therefore must be defined.  This is done using virtual inheritance (as described previously in the context of the simulation engine).  The order parameter re-implements the hA and hB indicator functions from the TpsOrderParameter class.  In the case of the Ising model, state A is defined as having a maximum nucleus size 5 or smaller and state B is defined as having a maximum nucleus size greater than 200:

step2.cpp (order parameter)

//Define a class to compute the nucleus size:
class TpsOrderParameterNucleusSize : public TpsOrderParameter
{
    public:
   
        bool hA(TpsTrajectory& traj)
        {
            //nucleus size at the beginning of the trajectory:
            int n = evaluate(traj, 0);
            
            //nucleus must be smaller than N=5 to be in basin A
            if (n < 5) {
                return true;
            }
            else {
                return false;
            }
        }

        bool hB(TpsTrajectory& traj)
        {   
            //nucleus size at the end of the trajectory:
            int n = evaluate(traj, traj.getNumberOfTimeslices()-1);
                        
            //nucleus must be larger than N=200 to be in basin B
            if (n > 200) {
                return true;
            }
            else {
                return false;
            }
        }
        // ... (functions to actually do the calculation)
};

Another point that requires clarification is the class called MyTpsAlgorithm.  The class overrides the normal TpsAlgorithmTPS class by adding a function that prints a trajectory for ever 10 trajectories accepted.  This is accomplished by overriding the virtual function called callbackOnTrialMoveAccepted(), which is called by the library every time a trial move (i.e., a shooting or shifting move) is accepted:

step2.cpp TPS algorithm

//Define a class for visualizing trajectories:
class MyTpsAlgorithm : public TpsAlgorithmTPS
{
    public: 
        MyTpsAlgorithm(TpsTrajectoryEnsemble& t, 
            TpsRNG& r, TpsOrderParameter& f, TpsInitializer& i)
        :    TpsAlgorithmTPS(t, r, f, i)
        {
            _count = 0;
        }
        
        //overwrite this function to visualize the path after every 10 steps:
        void callbackOnTrialMoveAccepted()
        {
            if (_count++ % 10 == 0) {
                //print the trajectory to a file
            }        
        }
        
    private:
        int _count;

};

The code is run by issuing the following shell commands:

step2

$ c++ -O3 step2.cpp TpsSimulationPluginLibLMC.cpp -llmc -ltps
$ ./a.out
$ vmd -m *.xyz

The resulting trajectories look similar to the initial trajectory, except that they evolve from their initial path resulting in a slightly different realization of the transition at each step.

Step 3: Collect Members of the Transition State Ensemble (TSE)

Each reactive trajectory must pass through a transition state on its way to transforming between stable states.  The transition state sits atop a free energy barrier and is therefore is equally likely to evolve towards either stable basin, A or B. Members of the transition state ensemble are therefore defined as microstates for which a randomly generated trajectory has a 50% chance to commit to either basin.  In the context of our Ising model example, the transitino state represents the critical nucleus.  The code to compute the transition states using the TPS library is shown below:

step3.cpp

    // ... same as in step 2 up to this point 
    // Compute the committer.  The settings are explained in libtps documentation
    TpsAnalysisCommitter committer(op, trajectory_length/5);
    committer.setUseShortcut(true);
    TpsAnalysisTransitionStates tse(committer, TpsAnalysisTransitionStates::STEPWISE);
    tse.setAlpha(2.0);
    tse.setNmax(50);

    // ... do path sampling    
    tse.analyze(tpe.getLastTrajectory());
    // print the transition state configuration to a file
    // ...

This code is the only difference between step 2 and step 3.  Compile and run step 3 as follows:

step3

c++ -O3 step3.cpp TpsSimulationPluginLibLMC.cpp -llmc -ltps
$ ./a.out
$ vmd transition_states.xyz

The panel below illustrates transition states for the Ising model under the conditions we consider.  The critical nucleus size is not particularly large, consisting of between 5 and 20 lattice sites.  

Notice that the nucleus evolves and moves around the box.  These are indications that our sampling of trajectory space is ergodic.

Step 4: Compute Committer and Determine Relevant Reaction Coordinates

The ultimate goal in applying TPS is often to determine the relevant reaction coordinates that guide the evolution of a system from state A to state B.  This topic is complex and its proper treatment is beyond the scope of this tutorial.  Instead, the reader is directed to the following references.
In a rough sense, the viability of a reaction coordinate is captured by crossover behavior in the committer function for a given order parameter.  The panel below illustrates the committer Pas a function of the nucleus size N.  Clearly, if N is < 10, the system almost never commits to B and if the nucleus size is > 30 the system almost always commits to B.  The critical nucleus size is approximately N=12.
Significantly more work must be done to determine if N is the only relevant reaction coordinate.  For example, probability distributions of the committer should be considered for different realizations of microstates with the critical value of the reaction coordinate.  If the distribution is sharply peaked at 50%, the order parameter is a good reaction coordinate.  As it turns out, at least one more reaction coordinate is necessary to properly describe nucleation in the Ising model.  This problem was studied in detail by Pan and Chandler who's paper provides a concrete calculation to accompany the proof of concept examples described here.
    Comments