This section contains two exercises that present situations in which the ubiquitin system is not at equilibrium. This section requires some knowledge of the principles of statistics 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 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. You will determine the thermal diffusivity by monitoring changes of the system's temperature and comparing them to the theoretical expression.
In order to use the temperature coupling feature of NAMD2, you need to create a file which marks the atoms subject to the temperature coupling.
The following lines
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
|mol load psf ubq_ws.psf namdbin ubq_ws_eq.restart.coor|
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
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
decay through non-linear contributions to forces between atoms. The
decay of coherence 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.13. At time the velocities of all
atoms in the system are quenched; then, at the atomic
velocities are reassigning 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:
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 (NEV) ensemble.
In this problem, you will consider the simplest temperature
echo: the set of velocities of the second reassignment is exactly the
same as the first one, i.e., without scaling.
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).
Let's take a look at some features in the configuration file:
# Continuing a job from the restart files
set inputname ../common/ubi_equil
You will now calculate the temperature autocorrelation function,
which will be used later to analyze the temperature echoes.
This will calculate the temperature autocorrelation function and put it in the file auto_tmp.dat.
A plot of the autocorrelation function will appear on your
screen (c.f. Fig. 14) in a dotted curve.
The solid curve is a fit to an exponential
function of Eq. 12. The decay (temperature
autocorrelation) time from this exponential is calculated
and displayed in the Command Window. This value will be
important for a later analysis of your simulations. The value is
stored in the file tau0.dat.
Once you know your value for , you are going to quench
the temperature (set it to zero), and monitor the recovery of the
system as it runs for time steps.
Now, you will quench the system for the second time.
You now have the temperature data for the whole temperature
quench experiment. You will use Matlab to analyze the data.
Note anything unusual? You quenched the temperature twice, and
yet you see three drops in temperature! The third one is a
Next, you will compare in the vicinity of the echo
() obtained from the MD simulation with the theoretical
prediction (14) 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.
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
At you should reassign the velocities according to a
Maxwell-Boltzmann distribution for a temperature . Make sure
that you save the reassigned velocities to a file 300.vel.
This 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: