In the previous lecture we've learned about the theoretical basis of chemical reaction rate theory. In particular, we looked at Transition State Theory (TST) (or the theory of the activated complex) and Kramers Theory (KT).
In this practical we are going to investigate in detail on these two theories.
To this end, we will compute molecular dynamics trajectories using the so-called Langevin equation,
which basically describes Newtonian dynamics with an additional friction term
with a time-independent friction constant plus a random force that ensures proper thermalization of our system.
The random force R(t) is Gaussian white noise with zero mean and obeys the fluctuation-dissipation relation. Note that the Langevin equation also provides the fundamental basis in Kramers Theory (see previous lecture).
All calculations performed in this practical will be done with a simple molecular dynamics
program written in Fortran 77/90 that solves the Langevin equation for single degree of freedom x (i.e. one particle) in a symmetric double minimum potential. In case you are interested in the details you can have a look at the source code langevin.f. To compute the rate after we have generated a trajectory we use a short analyzing program
called rate.f.
To start the practical download the relevant files from here followed by unpacking them with the following command
tar xzvf rate_prakt.tar.gz
Then change into the new directory and compile both programs needed with
gfortran langevin.f -o langevin gfortran rate.f -o rate
Now we are ready to run our first computation. Our first task will be to have a look
at the potential V(x) that will govern the deterministic part of our Langevin dynamics.
To plot the potential V(x) just start the program langevin by typing
./langevin md.in
and then start a Linux plot program, like for example gnuplot
gnuplot (RETURN) p "POTENTIAL" u 1:2 w l
(A symmetric double minimum potential should appear on the screen.)
The potential is a generic example for all kinds of reactive processes
that show no potential (or free energy) difference between the two
metastable regions A and B. The dynamical variable x is our order parameter
(or reaction coordinate, collective variable or meta-coordinate) that must be able
to distinguish between region A and B. Thinking in terms of macromolecules
it would be a collective coordinate like e.g. the number of native contacts,
the solvent-accesible surface or the radius of gyration. You can also view
it as a projection from 3N cartesian coordinates of your macromolecule
at hand onto a 1D (collective) coordinate. For example in the animation of
a distinct conformational transition in alanine dipeptide shown above one could
imagine to use an angle or a hydrogen-bond distance as an order parameter
to arrive at a reduced 1D description of the overall reactive process
described by the concerted movement of all 22 atoms shown.
The input file md.in of our MD program contains the necessary parameters shown here
4000 1 # nsteps,nout 0.005 # timestep 0.2 # x(t=0) 0.0 # temperature 0.0 # friction coefficient .F. # compute TST and Kramers rate
that will be read by the program langevin prior to the computation.
The parameters will be explained briefly in the following: nsteps are the total
number of molecular dynamics steps performed, nout gives the inverse
output frequency (ouptut is written only every n'th step), timestep is the discrete
time step used in the numerical integration algorithm, x(t=0) are rather self-explanatory, temperature is the reduced temperature and friction coefficient is the parameter gamma that enters the Langevin equation. Note that the program uses a reduced unit system with k_B = 1, i.e. the temperature is given
in units of energy.
Now we can start our first simulation with the input parameters and look at the resulting trajectory. Typing
./langevin md.in
there should now be a file TRAJECTORY
in your directory. We can display the content of this file by starting gnuplot and typing
p "TRAJECTORY" u 1:2 w l
This should reveal the time evolution (time on the x axis) of the coordinate x (shown on the y axis) that evolves in double minimum potential according to Newton's equation of motion (Note that the friction coefficient and temperature in the input file are both zero).
The signal we obtain with the previous input setting is representative of so-called ballistic dynamics (as opposed to diffusive dynamics we will
encouter later on). The particle oscillates between the two potential wells
with constant total energy (try in gnuplot p "TRAJECTORY" u 1:4 w l,"TRAJECTORY" u 1:5 w l to display the time evolution of the potential energy and total energy).
What is the rate k
for a process like this?
Now we can turn on the friction term by setting gamma in the md.in input file
to a finite value (e.g. 0.5) and starting the previous simulation with otherwise identical parameters. What behaviour do you see now? What is the difference to the previous example?
Try to run a trajectory with the initial coordinate x located at the barrier top and try to explain the observed behaviour.
Finally we are going to also turn on the random force (see Langevin equation) in order
to counteract the damping that is introduced by the friction term. What different physical mechanism are mimiced by the friction term in connection with the random force?
By setting temperature to 0.75 and starting the simulation again we obtain now a different time evolution of x. How many reactive event (crossing from A to B) do you observe for this short trajectory? Again try to run a trajectory with the initial coordinate x located at the barrier top and try to explain the observed behaviour. What is different compared to the previous example with R(t) = 0?
Being done with the preliminaries we can begin with computing reaction rates from our trajectories. To do so, we have to increase the number of integration steps nsteps
in md.in such that we observe a sufficiently large number of reactive events.
To compute the rate from our trajectory, we run a second program rate that reads
in TRAJECTORY file and performs a smart counting of reactive events, i.e. how often
does a the dynamical variable x move from region A to region B. We simply start
the program by typing
./rate
and it will give us the rate k which has dimension 1/(time) and units 1/(timestep), where
timestep is the specified timestep in md.in
Next we are going to analyze whether our simple 1D model shows the expected Arrhenius behaviour. In the previous lecture we learned that the reaction rate as a function of temperature follows a phenomenological law (Equation), which means that a plot ln(rate) vs. 1/T should reveal a straight line. To check this result we perform rate computations at four different temperatures (0.25 0.5 0.75 1.0) and from each extract the rate k and plot the combined data in gnuplot
p "arrhenius.dat" (1./$1):(log($2)) w lp
