Analysis of the SMD trajectory

Launch VMD by typing vmd in a terminal window.

Within VMD select the menu item Extension $\rightarrow$ Tk Console. A console window should appear with a prompt. The following Tcl/Tk commands are executed from within the TkCon window. We will use the tag » to indicate Tcl commands.

Make sure you are in the correct directory by typing:

» cd ~/Workshop/10Ala-tutorial/files/

Open the Tcl Output file generated by your SMD simulation, and store its content to a variable mytraj:

» set infile [open da_smd_tcl.out r]
» set mytraj [read $infile]
» close $infile

The content of the file da_smd_tcl.out is printed on the screen. You should see three columns:

First column: time (ps)
Second column: extension, or the end to end distance (Å)
Third column: applied force (kcal/mol/Å)

We use the following units: ps for time, Å for distance, kcal/mol for energy.

Assign columns of the mytraj to three separate arrays (lists), by running the following Tcl script:

» source read.tcl

The script defines three variables, t, z, f, corresponding to each column of mytraj respectively.

Define parameters, v (pulling velocity) and dt (time step of the data):

» set v 1
» set dt 0.1

Recall that one end of the molecule was constrained to the origin and the other end was constrained to a moving point. The distance between the two constraints was changed from 13 Å to 33 Å with the constant velocity $v$. Define an array (list) $c$ representing the distance between the two constraint points at each time step:

» set c {}
» set i 13
» while {$i < 33.1} { lappend c $i; set i [expr $i + $v * $dt] }

Plot $z-t$ and $c-t$ curves using the multiplot plugin within VMD. In the TkCon window type:

» package require multiplot
» set plot1 [multiplot -x $t -y $z -plot]
» $plot1 add $t $c -plot

Q: The extension $z$ generally lags behind $c$, the distance between the two constraints. How big is the lag in this case? What would you do if you want to reduce the lag further?

Figure: Extension versus time during the simulation, compared with the distance between the constraint point (straight line).

Now plot the force-extension curve:

» set plot2 [multiplot -x $z -y $f -plot]

Q: Is the force generally positive or negative? Why?

Figure: Force versus extension during the simulation.

The work done on the system during the pulling simulation is:

$$W(t) = v \int_0^{t'} {dt'} f(t')$$ (1)

where $v$ is the pulling velocity.

To calculate the work $w$ by numerical integration, run the following Tcl script:

» source calcwork1.tcl

The work is stored in a variable (list) called $w$.

If a pulling simulation is performed very slowly, then the process is reversible. The work done during such a reversible pulling is equal to the free energy difference of the system between initial and final states. Thus, a reversible work curve can be considered as an exact PMF. For our system, the reversible pulling speed turns out to be about 0.0001 Å/ps (10000 times slower than the one you used).

The exact PMF obtained from a reversible pulling is stored in the file Fexact.dat. Read the content of the file, and store it in a variable called Fexact, by running the follwong Tcl script:

» source load-Fexact.tcl

The array Fexact contains two columns, representing the potenial of mean force (second column) at each extension (first column).

Plot $W-c$ in blue together with $Fexact$ in red:

» set plot3 [multiplot -x $c -y $w -linecolor blue -plot]
» $plot3 add $Fexact(1) $Fexact(2) -linecolor red -plot

Figure: Work versus extension in the fast pulling simulation (blue) and reversible pulling simulation (red).

Q: How do they compare? What is the origin of the difference?

Quit VMD:

» quit

PMF calculation

The trajectory you have generated in this session is actually not practical for the PMF calculation because the pulling speed was too high. Besides, you need multiple trajectories to use Jarzynski's equality. Therefore, you will use pre-generated trajectories.

Launch VMD by typing vmd in a terminal window.

Within VMD select the menu item Extensions $\rightarrow$ Tk Console. The following commands are all executed from within the TkCon window.

Make sure you are in the correct directory by typing:

» cd ~/Workshop/10Ala-tutorial/files/

Ten trajectories obtained from SMD simulations with a pulling speed of 0.01 Å/ps are stored in da.dat. This pulling speed is 100 times slower than the one you used, but still 100 times faster than the reversible speed. Load the trajectories as well as the exact PMF:

» source load-Fexact.tcl
» source load-traj.tcl

The script load-traj.tcl defines three variable. The array $t$ contains time steps of the data. $z$ and $f$ are matrices of 10 columns, with each column representing the extension and the applied force for each trajectory, respectively.

To plot the extension-time curve for each trajectory, you need to reduce the number of data points plotted for each graph. The following script plots $z$ vs $t$ for a tenth of the data points. Each trajectory is plotted with a different color:

» source plot-extension-time.tcl

Figure: Extension versus time for 10 pulling trajectories.

And plot the force-extension curve:

» source plot-force-extension.tcl

In each figure window, you should see ten overlapping curves, with each curve corresponding to a trajectory.

Figure: Force versus extension for 10 pulling trajectories.

Define parameters: $dt$, time step; $v$, pulling speed; $T$, temperature (including Boltzmann's constant).

» set dt 0.1
» set v 0.01
» set T 0.6

Define an array $c$ representing the distance between the two constraint points at each time step:

» set c {}
» set i 13
» while {$i < 33.001} {lappend c $i; set i [expr $i + $v * $dt] }

Calculate work for each trajectory by executing the following command:

» source calcwork2.tcl

Figure: Work versus extension for 10 pulling trajectories. Fexact, obtained from reversible pulling is shown in bold.

Plot $W$ vs. $c$ together with the exact PMF:

» source plot-W.tcl

Plot the exact PMF on the same graph:

» $plot3 add $Fexact(1) $Fexact(2) -linewidth 3 -linecolor black -plot

You should see ten work curves and one curve for the exact PMF. The thick black line is the exact PMF.

Q: Work done to the system should be on average higher than the PMF. Are all the work curves higher than the PMF, or do you also see some work curves lower than the PMF?

Let's now employ Jarzynski's equality:

$$\exp{\{\frac{-[F(c(t)) - F(c(0))]}{T}\}} =~< \exp[-W(t)/T] >$$ (2)

or

$$F(c(t)) - F(c(0)) = - T \log < \exp[-W(t)/T] >$$ (3)

The right-hand side can be expanded in terms of so-called cumulants:

$$F(c(t)) - F(c(0)) = < W(t) > - \frac{1}{2T} [ < {W(t)}^2 > - {< W(t) >}^2 ] + \cdots$$ (4)

where the cumulants up to the second order are shown.

Estimate the PMF according to three different formulas: the original formula (exponential average), the first order cumulant expansion formula (average work), and the second order cumulant expansion formula. The following Tcl script will calculate and store the three different estimates in Fexp, F1, F2 variables.

» source cumulants.tcl

Plot the three estimates for the PMF together with the exact PMF. Plot the exact PMF, Fexact, in red, the average work (or the first order cumulant expansion), $F1$, in green, the second order cumulant expansion, $F2$, in blue, and the exponential average, Fexp, in black:

» set plot4 [multiplot -plot]
» $plot4 add $Fexact(1) $Fexact(2) -linecolor red -plot
» $plot4 add $c $F1 -linecolor green -plot
» $plot4 add $c $F2 -linecolor blue -plot
» $plot4 add $c $Fexp -linecolor black -plot

Q: Which estimate is the closest to the exact PMF? Does it make sense to you in view of what you learned in the lecture?

Figure: The PMF obtained from the cumulant expansion of Jarzynski's equality, compared with the PMF obtained from reversible pulling simulation (shown in red).

Quit VMD:

» quit

For more information on the PMF calculation from SMD simulations, see the paper included (paper.pdf).