Next: Steered Molecular Dynamics
Up: Analysis
Previous: Equilibrium
Subsections
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.
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.
- * 1
- Change to this problem's directory in a Terminal window:
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.
- 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 , and 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:
- 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
- * 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 Transformations Geometrical transform 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 -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 Transformations Non-linear curve fitting... menu item. In the Source 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.
|
- * 21
- Multiply the value of the a0 parameter by 0.1 to get the
thermal diffusivity in cms. You should get a value of around
0.45cms. How does your result compare
to the thermal diffusivity of water
cms?
- * 22
- Compute the thermal conductivity of the system, , by the following
formula:
, where is the specific heat (that was
calculated in Section 2.1.5) and is the thermal diffusivity
(assume that the density of the system, , is 1 g/mL).
- * 23
- When you are finished investigating the plot, close xmgrace, saving your plot if you wish.
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 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 ,
as a result of such reassignment. An example is shown in
Fig. 14. At time the velocities of all
atoms in the system are quenched; then, at the atomic
velocities are reassigned again (quenched) to their values at time
. As a result, a temperature echo, i.e., a sharp dip in
, is detected at a subsequent time after
(i.e., at time ).
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:
- the dephasing time of the oscillators is about one picosecond;
- the temperature echo can be expressed in terms of the temperature
autocorrelation function; and
- the depth
of the echoes decay exponentially
with the delay time , i.e.,
, the exponential decay being
determined by the dephasing time .
Figure 14:
Temperature quench echo.
|
First you will generate temperature echoes in MD simulations, and then you will analyze the echoes in the framework of the normal mode approximation.
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
ps.
The second velocity reassignment after a delay time 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 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.
The temperature quench echo is obtained by resetting to zero
the velocity of all atoms of the protein at both times and
. You should start from the pre-equilibrated protein in
vacuum at ; 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.
- The simulation takes restart files, i.e., the coordinate file and the velocity file, from an equilibration simulation at .
# Continuing a job from the restart files |
|
set inputname ../common/ubi_equil |
|
binCoordinates $inputname.restart.coor |
|
extendedSystem $inputname.xsc
|
|
- Notice that the temperature line is commented out. The
initial temperature is pre-defined by the initial velocities.
#temperature $temperature |
|
binvelocities $inputname.restart.vel
|
|
- Since you are running an NVE simulation, there is no temperature
coupling or pressure control.
- As discussed previously, one should use rigidBonds for water in any production simulation since water molecules have been parametrized as rigid molecules. But for the illustration of temperature quench echo, the effect of force field can be negligible. To be simple, the rigidBonds option is turned off.
- * 2
- Close the configuration file and nedit.
- 3
- Run the Molecular Dynamics simulation
- * 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.
- * 6
- Once your job is done, you need to get the temperature data from the output file. In the VMD TkCon window (Extensions 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.
- * 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:
This will put the temperature autocorrelation function in the file auto-corr.dat.
- * 10
- In a Terminal window, type
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 Transformations Non-linear curve fitting... menu item. In the Source 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:
- * 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, , 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).
|
- * 16
- Go to the next directory, 02_quencha, by typing in a Terminal window:
- * 17
- You will perform the temperature quenching experiment for the value . You are going to quench the temperature (set it to zero), and monitor the recovery of the system as it runs for time steps. If you would like to probe other values of , you should edit your configuration file echo.conf using nedit, changing the value of tau:
- 18
- Run the Molecular Dynamics simulation
- * 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:
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 steps: run [expr 3*$tau]
- 23
- Run the Molecular Dynamics simulation
- * 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:
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:
- * 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:
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.
Next, you will compare in the vicinity of the echo () 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:
This will write the data for the harmonic approximation to the file harmonic.dat.
- * 31
- In xmgrace, view the data by clicking Data Import ASCII 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 =300K.
Figure 16:
A comparison of the simulation and the expression from harmonic approximation at the vicinity of the echo.
|
You have calculated the echo depth for a particular value of the delay time . You may wish to repeat the experiment using other values of (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
, one obtains the so-called dephasing time , which is an inherent property of the system.
- * 34
- When you are finished investigating the plot, close xmgrace, saving your plot if you wish.
Here you will generate a temperature echo by replacing the
velocities of the atoms at time with the ones you assign to
them at time . Moreover, by choosing these velocities according
to the Maxwell-Boltzmann distribution corresponding to K,
i.e., the temperature of the equilibrated system (), there will
be no discontinuities (jumps) in at and , yet you
will observe a temperature echo at a later time, i.e., a sharp dip in
the temperature at time . Your goal is to determine
and the depth 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 according to a
Maxwell-Boltzmann distribution for a temperature , 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:
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:
- 37
- Run the Molecular Dynamics simulation
- * 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 time steps.
- 42
- Run the Molecular Dynamics simulation
- * 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:
- * 49
- Try to locate the time of the echo and find its depth by looking at the temperature data as before:
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 fs. How would you explain your findings?
- * 51
- You can now close VMD and xmgrace.
Next: Steered Molecular Dynamics
Up: Analysis
Previous: Equilibrium
namd@ks.uiuc.edu