next up previous
Next: Steered Molecular Dynamics Up: Analysis Previous: Equilibrium

Subsections


Non-equilibrium

In this section, you will analyze two non-equilibrium properties of proteins, the diffusion of heat and a phenomenon known as temperature echoes. This section requires some knowledge of the principles of statistical physics, but the results of the simulations can be understood intuitively even without a deep understanding of their principles.

Heat Diffusion

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...s of your simulation, the thermal diffusivity of your system.}
\end{minipage} }

For this exercise, you will use the last step of the equilibration of ubiquitin in a water sphere, file ubq_ws_eq.restart.coor. You will use the temperature coupling feature of NAMD to set the temperature of the molecules in the outer layer of the sphere to 200 K, while the rest of the bubble is set to a temperature of 300 K. In this way, you will determine the thermal diffusivity by monitoring change of the system's temperature and comparing it to the theoretical expression.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...ft[-\left(\frac{n\pi}{R}\right)^2Dt\right].
\end{displaymath}}
\end{minipage} }

* 1
Change to this problem's directory in a Terminal window:

cd ../2-6-heat-diff/
 

In order to use the temperature coupling feature of NAMD, you need to create a file which marks the atoms subject to the temperature coupling. The following lines

tCouple on  
tCoupleTemp 200  
tCoupleFile ubq_shell.pdb  
tCoupleCol B
 

in the provided configuration file ubq_cooling.conf will enable this feature, and set the atoms coupled to the thermal bath as the ones that have a value of 1.00 in the B column of the ubq_shell.pdb, which you will create below.

* 2
In the VMD TkCon window, navigate to the 2-6-heat-diff/ directory.

cd ../2-6-heat-diff/
 

3
Load the system into VMD by typing the following in the VMD TkCon window:

mol new ../common/ubq_ws.psf  
mol addfile ../1-2-sphere/ubq_ws_eq.restart.coor
 

4
Select all atoms in the system:

set selALL [atomselect top all]
 

5
Find the center of the system:

set center [measure center $selALL weight mass]
 

6
Find $x$, $y$ and $z$ coordinates of the system's center and place their values in the variables xmass, ymass, and zmass:

foreach {xmass ymass zmass} $center { break }
 

7
Select atoms in the outer layer:

set shellSel [atomselect top "not ( sqr(x-$xmass)  
+ sqr(y-$ymass) + sqr(z-$zmass) <= sqr(22) ) "]
 

8
Set beta parameters of the atoms in this selection to 1.00:

$shellSel set beta 1.00
 

9
Create the pdb file that marks the atoms in the outer layer by ``1.00'' in the beta column:

$selALL writepdb ubq_shell.pdb
 

10
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...:
\\ \\
{\tt namd2 ubq\_cooling.conf > ubq\_cooling.log \&}
}
\end{minipage} }

* 11
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/ubq_cooling.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 12
Use namdstats.tcl to find out how the system's temperature changes with time. In the VMD TkCon window, type:

source ../2-3-energies/namdstats.tcl  
data_time TEMP ubq_cooling.log
 

The file TEMP.dat will be written with the temperature of the system over time.

* 13
You can now quit VMD.

* 14
In a Terminal window, type xmgrace TEMP.dat

* 15
In the xmgrace program, choose the Data  $\rightarrow$ Transformations  $\rightarrow$ Geometrical transform $\ldots$ menu item. In the Scale X: window, type 2.0, which corresponds to the time step used in your simulation. Press Accept, which will rescale the plot such that the $x$-axis is now in units of fs. Now, press the AS button in the upper left hand corner of the xmgrace window to auto-scale your plot.

* 16
Double-click on the y-axis of the plot and enter a value of 320 in the Stop field. Click Accept.

* 17
You will now fit the simulated temperature dependence to the theoretical expression. For this, you will use the Non-linear curve fitting feature of xmgrace.

* 18
Choose the Data  $\rightarrow$ Transformations  $\rightarrow$ Non-linear curve fitting... menu item. In the Source  $\rightarrow$ Set window, click on the last line, which corresponds to the simulated time dependence.

* 19
In the Main tab, type or copy and paste over the Formula window (make sure it is on one line):

y = 200 + 66.87*(exp(-0.0146*a0*x) +0.25*exp(-0.25*0.0146*a0*x)  
+1/9*exp(-1/9*0.0146*a0*x) +1/16*exp(-1/16*0.0146*a0*x)  
+1/25*exp(-1/25*0.0146*a0*x) +1/36*exp(-1/36*0.0146*a0*x)  
+1/49*exp(-1/49*0.0146*a0*x) +1/64*exp(-1/64*0.0146*a0*x)  
+1/81*exp(-1/81*0.0146*a0*x) +1/100*exp(-1/100*0.0146*a0*x))
 

* 20
In the Parameters drop-down menu, choose 1. Windows will appear below. Click several times on Apply, to get a better fit. This will fit the curve and get a value for the parameter a0, which is the thermal diffusivity.

Figure 13: Cooling of ubiquitin in a water sphere. The red line shows a non-linear fit of the data to the theoretical expression.
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_unit02_coolunix}
}
\end{center}
\end{figure}

* 21
Multiply the value of the a0 parameter by 0.1 to get the thermal diffusivity in cm$^2$s$^{-1}$. You should get a value of around 0.45$\times10^{-2}$cm$^2$s$^{-1}$. How does your result compare to the thermal diffusivity of water $D=1.4\times10^{-3}$cm$^2$s$^{-1}$?

* 22
Compute the thermal conductivity of the system, $K$, by the following formula: $K=\rho\: c_V D$, where $c_V$ is the specific heat (that was calculated in Section 2.1.5) and $D$ is the thermal diffusivity (assume that the density of the system, $\rho$, is 1 g/mL).

* 23
When you are finished investigating the plot, close xmgrace, saving your plot if you wish.

Temperature echoes

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...teins by generating
temperature echoes using MD simulations.}
\end{minipage} }

The motions of atoms in globular proteins (e.g., ubiquitin), referred to as internal dynamics, comprise a wide range of time scales, from high frequency vibrations about their equilibrium positions with periods of several femtoseconds to slow collective motions which require seconds or more, leading to deformations of the entire protein.

The internal dynamics of these proteins on a picosecond time scale (high frequency) can be described as a collection of weakly interacting harmonic oscillators referred to as normal modes. Since normal modes are formed by linear superposition of a large number of individual atomic oscillations, it is not surprising that the internal dynamics of proteins on this time scale has a delocalized character throughout the protein. The situation is similar to the lattice vibrations (phonons) in a crystalline solid. Experimentally, there exist ways to synchronize, through a suitable signal or perturbation, these normal modes, forcing the system in a so-called (phase) coherent state, in which normal modes oscillate in phase. The degree of coherence of the system can be probed with a second signal which through interference with the coherent normal modes may lead to resonances, referred to as echoes, which can be detected experimentally. However, the coherence of atomic motions in proteins decays through non-linear contributions to forces between atoms. This decay develops on a time scale $\tau_d$ which can be probed, e.g., by means of temperature echoes, and can be described by employing MD simulations.

In a temperature echo the coherence of the system is probed by reassigning the same atomic velocities the system had at an earlier time and then looking to an echo in the temperature at time $\tau_e$, as a result of such reassignment. An example is shown in Fig. 14. At time $t_1=0$ the velocities of all atoms in the system are quenched; then, at $t_2=\tau$ the atomic velocities are reassigned again (quenched) to their values at time $t_1=0$. As a result, a temperature echo, i.e., a sharp dip in $T(t)$, is detected at a subsequent time $\tau_e = \tau$ after $t_2$ (i.e., at time $t=2\tau$).

In this section, you will employ MD simulations to generate temperature echoes in ubiquitin by applying the velocity reassignment just described. By modeling ubiquitin as a large collection of weakly interacting harmonic oscillators (normal modes), you will find that:

Figure 14: Temperature quench echo.
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_quench-schem}
}
\end{center}
\end{figure}

First you will generate temperature echoes in MD simulations, and then you will analyze the echoes in the framework of the normal mode approximation.

% latex2html id marker 11786
\fbox{
\begin{minipage}{.2\textwidth}
\includegr...
...o as temperature echo(es)
(Fig.~\ref{fig:tut_quench-schem}).
}
\end{minipage} }

The phenomenon of temperature echoes has a simple explanation within the normal modes description of a protein. The first velocity reassignment enforces phase coherence for the oscillators. Because of the frequency dispersion of the normal modes (as well as the deviation from the harmonic approximation) the phase coherence of the protein will decay in time with a characteristic dephasing time $\tau_{d} \sim 1$ps. The second velocity reassignment after a delay time $\tau$ is a ``probing signal'' which will test the degree of coherence of the system at the instant of time it was applied. The depth of the echo and the instant of time at which it occurs are quantitative characteristics of the coherence of the internal dynamics of proteins.

In order to generate temperature echoes, one needs to equilibrate ubiquitin at the desired initial temperature $T_0=300$K, e.g., by using the velocity rescaling method. Once the system is equilibrated, the thermostat is removed and all the following simulations are carried out in the microcanonical (NVE) ensemble.

In this section, you will consider two types of temperature echoes. In the first, you will reassign all atomic velocities to zero, and in the second, you will reassign the atomic velocities to a 300 K distribution (the same as the initial temperature of the system) and the second reassignment will be exactly the same as the first.

Temperature quench echo

The temperature quench echo is obtained by resetting to zero the velocity of all atoms of the protein at both times $t_1=0$ and $t_2=\tau$. You should start from the pre-equilibrated protein in vacuum at $T_0=300K$; the required files are located in the common directory.

You need to run the simulations in the microcanonical ensemble (NVE) by using the configuration file equil.conf in the 2-7-echoes/01_equil_NVE directory. The simulation will run for 500 time steps (fs).

* 1
Change to the 2-7-echoes/01_equil_NVE directory in a Terminal window.

Let's take a look at some features in the equil.conf configuration file. You can read the file using nedit.

* 2
Close the configuration file and nedit.

3
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...erminal window:
\\ \\
{\tt namd2 equil.conf > equil.log \&}
}
\end{minipage} }

* 4
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/equil.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 5
Open VMD by typing vmd in the Terminal window.

\fbox{
\begin{minipage}{.17\textwidth}
\includegraphics[width=2.0 cm, height=2...
...on the VMD icon under {\tt Applications} in the {\tt Finder}.}
\end{minipage} }

* 6
Once your job is done, you need to get the temperature data from the output file. In the VMD TkCon window (Extensions $\rightarrow$ Tk Console), navigate to the 2-7-echoes/01_equil_NVE directory if you are not already there.

* 7
Get the temperature data at every timestep using our namdstats.tcl script:

source ../../2-3-energies/namdstats.tcl  
data_time TEMP equil.log
 

This will write the data into the file TEMP.dat.

* 8
You will now calculate the temperature autocorrelation function, which will be used later to analyze the temperature echoes.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...0)]>^2}
\nonumber \\ &\approx &exp(-t/\tau_0)
\end{eqnarray*}}
\end{minipage} }

* 9
A tcl script has been provided for you to calculate the temperature autocorrelation function. Source the script by typing in the VMD TkCon window:

source auto-corr.tcl
 

This will put the temperature autocorrelation function in the file auto-corr.dat.

* 10
In a Terminal window, type

xmgrace auto-corr.dat
 

A plot of the temperature autocorrelation function will appear.

You will now fit the simulated temperature autocorrelation function to a decaying exponential approximation (see the Science Box). For this, you will use the Non-linear curve fitting feature of xmgrace.

* 11
Choose the Data  $\rightarrow$ Transformations  $\rightarrow$ Non-linear curve fitting... menu item. In the Source  $\rightarrow$ Set window, click on the only data set shown, which corresponds to the autocorrelation function.

* 12
In the Main tab, type the following into the Formula window:

y = exp(-x/a0)
 

* 13
In the Parameters drop-down menu, choose 1. Windows will appear below. Click several times on Apply, to get a better fit. This will fit the curve and get a value for the parameter a0. This parameter is the decay (temperature autocorrelation) time, $\tau_0$, and it will be important for a later analysis of your simulations.

14
If you prefer, double-click on the title area (above the graph), the x-axis, and the y-axis to enter suitable labels for your graph.

* 15
When you are finished investigating the plot, close xmgrace, saving your plot if you wish.

Figure 15: Temperature autocorrelation function from simulation (black) and a decaying exponential approximation (red).
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_autotempunix}
}
\end{center}
\end{figure}

* 16
Go to the next directory, 02_quencha, by typing in a Terminal window:

cd ../02_quencha
 

* 17
You will perform the temperature quenching experiment for the value $\tau = 200 fs$. You are going to quench the temperature (set it to zero), and monitor the recovery of the system as it runs for $\tau$ time steps. If you would like to probe other values of $\tau$, you should edit your configuration file echo.conf using nedit, changing the value of tau:

set tau 200
 

18
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... Terminal window:
\\ \\
{\tt namd2 echo.conf > echo.log \&}
}
\end{minipage} }

* 19
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/echo.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 20
Once your job is done, you need to again get the temperature data from the output file. In the VMD TkCon window, navigate to the 2-7-echoes/02_quencha directory.

* 21
Get the temperature data at every timestep using our namdstats.tcl script, which we have already sourced:

data_time TEMP echo.log
 

This will write the data into the file TEMP.dat.

Now, you will quench the system for the second time.

* 22
Change to the directory 03_quenchb in a Terminal window.

Note that the configuration file echo.conf indicates that the simulation will run for $3\tau$ steps: run [expr 3*$tau]

23
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... Terminal window:
\\ \\
{\tt namd2 echo.conf > echo.log \&}
}
\end{minipage} }

* 24
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/echo.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 25
Again, get the temperature data from the output file. In the VMD TkCon window, navigate to the 2-7-echoes/03_quenchb directory and get the temperature data at every timestep using our namdstats.tcl script:

data_time TEMP echo.log
 

This will write the data into the file TEMP.dat.

* 26
In a Terminal window, change the name of this output file so as not to confuse it with the others:

mv TEMP.dat temp3.dat
 

* 27
Copy the previous temperature output to this directory, changing their names as you do so:

cp ../01_equil_NVE/TEMP.dat ./temp1.dat  
cp ../02_quencha/TEMP.dat ./temp2.dat
 

* 28
Merge the three temperature data files into temp.dat:

cat temp[1-3].dat > temp.dat
 

You now have the temperature data for the whole temperature quench experiment. You will use xmgrace to analyze the data.

* 29
Examine the temperature profile across all three simulations by typing at a Terminal window:

xmgrace temp.dat
 

Note anything unusual? You quenched the temperature twice, and yet you see three drops in temperature! The third one is a temperature echo. See Fig. 16.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...rac{1}{2}C_{T,T}(\vert t-2\tau\vert)\right]\;.
\end{equation}}
\end{minipage} }

Next, you will compare $T(t)$ in the vicinity of the echo ($t=2\tau$) obtained from the MD simulation with the theoretical prediction (2) involving the temperature autocorrelation function obtained also from the MD simulations. The degree of agreement between these two results is a measure of the accuracy of the harmonic approximation.

* 30
You will use a tcl script to calculate the harmonic approximation of the temperature echo, given by Eq. 1. Source the script by typing in the VMD TkCon window:

source harmonic.tcl
 

This will write the data for the harmonic approximation to the file harmonic.dat.

* 31
In xmgrace, view the data by clicking Data $\rightarrow$ Import $\rightarrow$ ASCII$\ldots$ and choosing the file harmonic.dat. Click OK. Then, double click on the x-axis and change the axis Stop value to 1400 for better visualization.

You need to find the depth of the echo. From the graph, you can see that the echo occurs around 900 fs.

* 32
Open the file temp.dat by typing less temp.dat.

* 33
Scroll down to the temperature data near 900 fs, and find the minimum local value of temperature around that time. Record the value and calculate the depth of the echo by subtracting the value from 75:

Echo depth = 75 - Minimum Temperature
 

75 is the average temperature after the second quench, and it corresponds to 1/4 of the initial temperature $T_0$=300K.

Figure 16: A comparison of the simulation and the expression from harmonic approximation at the vicinity of the echo.
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_quenchechounix}
}
\end{center}
\end{figure}

You have calculated the echo depth for a particular value of the delay time $\tau = 200 fs$. You may wish to repeat the experiment using other values of $\tau$ (100-800 fs). (You must alter the harmonic.tcl script.) If this is done and the data is plotted and fitted to a single exponential $\exp(-\tau/\tau_{d})$, one obtains the so-called dephasing time $\tau_{d}$, which is an inherent property of the system.

* 34
When you are finished investigating the plot, close xmgrace, saving your plot if you wish.

Velocity replacement echo

Here you will generate a temperature echo by replacing the velocities of the atoms at time $t_2$ with the ones you assign to them at time $t_1$. Moreover, by choosing these velocities according to the Maxwell-Boltzmann distribution corresponding to $T_0=300$K, i.e., the temperature of the equilibrated system ($t<0$), there will be no discontinuities (jumps) in $T(t)$ at $t_1$ and $t_2$, yet you will observe a temperature echo at a later time, i.e., a sharp dip in the temperature at time $t=\tau+\tau_e$. Your goal is to determine $\tau_e$ and the depth $\Delta T$ of the echo. To this end, you will need to follow a similar procedure to the one in the previous exercise.

* 35
In a Terminal window, navigate to the directory 04_consta.

* 36
Reassign the velocities at $t_1$ according to a Maxwell-Boltzmann distribution for a temperature $T_1=300K$, saving the reassigned velocities to a file 300.vel. In order to do this, the following is added to the end of your configuration file echo.conf:

run 0  
output 300
 

These settings will not really run a simulation, but will assign initial velocities with the desired distribution and save them to 300.vel. Following these, the normal run command is used to perform a simulation for tau time steps:

run $tau
 

37
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... Terminal window:
\\ \\
{\tt namd2 echo.conf > echo.log \&}
}
\end{minipage} }

* 38
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/echo.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 39
In the VMD TkCon window, navigate to the 04_consta directory. Use the script namdstats.tcl to get the temperature data as before and put it in the file TEMP.dat. Note that at the beginning there is a duplicate entry for the first time step. You should delete one of these using your text editor:

500 300.5656  
500 300.5656  
501 301.0253  
502 302.5395  
...
 

* 40
In a Terminal window, go to the directory 05_constb.

41
For this simulation, you should use the file 300.vel you generated before as the velocity restart file. This is included in the configuration file echo.conf as:

velocities ../04_consta/300.vel
 

This will reassign the velocities to the exact same distribution they had at the beginning of the previous simulation. This simulation will run for 3$\tau$ time steps.

42
Run the Molecular Dynamics simulation

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... Terminal window:
\\ \\
{\tt namd2 echo.conf > echo.log \&}
}
\end{minipage} }

* 43
If you do not run the simulation yourself, you will need to obtain the data for the temperature from the provided log file. Copy this log file from the directory example-output/.

cp example-output/echo.log .
 

If you run the simulation yourself, your log file will be generated by NAMD, and you do not need to copy it.

* 44
Again, in the VMD TkCon window, navigate to the 05_constb directory and use the script namdstats.tcl to output the temperature data to the file TEMP.dat.

* 45
In a Terminal window, change the name of this output file and copy the other temperature file to this directory:

mv TEMP.dat temp2.dat  
cp ../04_consta/TEMP.dat ./temp1.dat
 

* 46
Merge the two temperature data files into temp.dat:

cat temp[1-2].dat > temp.dat
 

* 47
Note that there are now two temperature values at 700 fs. You should open the file temp.dat with your text editor and delete one.

* 48
View the temperature data and observe the temperature echo using
xmgrace:

xmgrace temp.dat
 

* 49
Try to locate the time of the echo and find its depth by looking at the temperature data as before:

less temp.dat
 

This time, the depth of the echo is measured by:

Echo depth = 300 - Minimum Temperature
 

* 50
Compare your result with the pre-made results stored in the directory example_output for $\tau = 200$ fs. How would you explain your findings?

* 51
You can now close VMD and xmgrace.


next up previous
Next: Steered Molecular Dynamics Up: Analysis Previous: Equilibrium
namd@ks.uiuc.edu