This section contains exercises that will let you analyze properties of a system in equilibrium. The exercises mostly use pre-equilibrated trajectories from different ensembles to calculate properties of the ubiquitin system.
You will use a script within VMD that will allow you to compute the average RMSD of each residue in your protein, and assign this value to the B column of the pdb file, which is typically reserved for the temperature factor of each residue.
In the previous unit, you had VMD open with the equilibration trajectory of the ubiquitin cubic system loaded. You will repeat this process:
source residue_rmsd.tcl |
This command itself does not actually perform any calculations. Instead, it executes the script residue_rmsd.tcl which contains a procedure called rmsd_residue_over_time. Calling this procedure will calculate the average RMSD for each residue you select over all the frames in a trajectory. The procedure is called as:
rmsd_residue_over_time molsel_resid |
set sel_resid [[atomselect top "protein and alpha"] get resid] |
rmsd_residue_over_time top $sel_resid |
The procedure also sets the value of the B column of all the atoms of the residues in the selection to the computed RMSD value. You will now color the protein according to this value. You will be able to recognize which residues are free to move more and which ones move less during the equilibration.
In order to use these new values to learn about the mobility of the residues, you will create a representation of the protein colored by B value.
You should now see the protein colored according to average RMSD values. The residues displayed in blue are more mobile while the ones in red move less. This is counter intuitive, so you will change the color scale.
xmgrace residue_rmsd.dat |
For this exercise, you will use the last step of the equilibration of ubiquitin. Specifically, you need the file that contains the velocities of every atom in the last frame of this equilibration. The name of the file is ubq_wb_eq.restart.vel, you created it in directory 1-3-box/. We will get the necessary information from this file and the structure file ubq_wb.psf with the help of VMD.
The molecule you loaded looks terrible! That is because VMD is reading the velocities as if they were coordinates. It does not matter how it looks, because what we want is to make VMD write a useful file: one that contains the masses and velocities for every atom.
set all [atomselect top all] |
set fil [open energy.dat w] |
foreach m [$all get mass] v [$all get {x y z}] { | |
puts $fil [expr 0.5 * $m * [vecdot $v $v] ] | |
} |
close $fil |
nedit energy.dat |
0.39707419313 | |
0.391385065921 | |
... |
vmd -dispdev text -e get_energy.tcl |
You now have a file of raw data that you can use to fit the obtained energy distribution to the Maxwell-Boltzmann distribution. You can obtain the temperature of the distribution from that fit. In the rest of this section, we will show how to do this with the 2-D graphics package xmgrace. If you are more familiar with another tool, you may use it instead.
Now, you will fit a Maxwell-Boltzmann distribution to this plot and find out the temperature that the distribution corresponds to. (Hopefully, the temperature you wanted to have in your simulation!)
y = (2/ sqrt(Pi * a03 )) * sqrt(x) * exp (-x / a0) |
The deviation from your value and the expected one is due to the poor sampling as well as to the finite size of the protein; the latter implies that temperature can be ``measured'' only with an accuracy of about even in case of perfect sampling.
You will look at the kinetic energy and all internal energies: bonded (bonds, angles and dihedrals) and non-bonded (electrostatic, van der Waals). You should remember that molecular dynamics treats bonds, angles, dihedrals and impropers as harmonic oscillators. You will determine the values of these energies at different temperatures, and find their dependence on this parameter.
To get values over a range of temperatures, you need to look at the output of several simulations, with the same conditions but different temperature. You were provided with a sample configuration file energies.conf for temperature K. The content of the file is the same as the equilibration configuration file you ran before, except that this time it will have a different value for the temperature. You were assigned a temperature according to your last name.
You need to change the temperature in the script for the temperature you were asked to simulate.
nedit energies.conf |
set temperature your_temperature |
When your job simulation is completed, you will need to get the average value of the energies at
. As you saw, the output file from NAMD has a lot of information. A script called namdstats does all the work for you! The usage of namdstats is:
namdstats [from first][to last] file |
which calculates the average for all output variables in file from frame first to frame last.
namdstats from 101 energies.log |
For this exercise, you need a simulation that is long enough to samle well the quantity of interest, here the fluctuations in kinetic energy, and hence in the temperature. Instead of doing such a simulation yourself, we recommend that you use the provided file ubq-nve.log. However, if you choose to run your own simulation, you can find all the necessary files and instructions below. We recommend you to read the whole section even if you don't perform the simulations.
The simulation is performed in the microcanonical ensemble (NVE, i.e., constant N, V and E). The configuration file provided to start the NAMD run is called ubq-nve.conf. The main features of this configuration file are:
namddat TEMP ubq-nve.log |
In order to plot the data, you can use xmgrace. You need to remove the line with the title as well as the first line of data, that corresponds to the initial temperature T=0.
nedit data.dat |
Does this distribution look any familiar? Your distribution should look like a Gaussian distribution.
This will fit the curve and get values for the parameters a0, a1 and a2, where a0 is a normalization constant, a1 the average temperature, and a2 is .
The window that appeared contains the new value of the parameters, as well as some statistical measures, including the Correlation Coefficient, which is a measure of the fit.
namdstats ubq-nve.log |
Take home message: You may note from this example that single proteins, on the one hand, behave like infinite thermodynamic ensembles, reproducing the respective mean temperature; and, on the other hand, they show signs of their finiteness, i.e. temperature fluctuations. In general, it is interesting to note that the velocities of the atoms and, therefore, the kinetic energy of the protein serves as its thermometer!
Specific heat is an important property of thermodynamic systems. It is the amount of heat needed to raise the temperature of an object per degree of temperature increase per unit mass, and a sensitive measure of the degree of order in a system. (Near a bistable point, the specific heat becomes large.)
For this exercise, you need to carry out a long simulation in the NVT (canonical) ensemble that samples sufficiently the averages
and
that arise in the definition of . To produce the data for the total energy of the molecule over time, you need to work through the steps presented below. This exercise is long and computationally intensive. We strongly recommend that you work only through the steps marked with *, that correspond to the analysis of the data. Nevertheless, we encourage you to read through the whole method since it will provide you with tools for future analysis.
In order to produce an NVT run, you will restart the equilibration you performed in Unit 1 for the water sphere. All necessary files are present in your directory.
You have now an equilibration run of ubiquitin in the NVT ensemble, and you want to look at the fluctuation in the energy of ubiquitin. The problem is that in order to calculate the specific heat for the protein only (as opposed to the whole system), you need the total energies for the atoms in the protein only, and the NAMD output for the run you just performed contains the energies for all atoms in the system.
You would like NAMD to write the energies for only a part of the system. This feature is currently available in NAMD only for interaction energies between two parts of a system. A way to accomplish this would be to have a system with the protein only, but it does not make sense to equilibrate the protein in this way. However, using a trick, you can make NAMD write the energies for the protein only given the equilibration you already performed.
You will first create a dcd with the protein only. For this, you will use a VMD script called mk_protein_dcd.tcl. This script writes a pdb for each frame of the trajectory, then loads all the pdb files, and writes a dcd file.
vmd -dispdev text -e mk_protein_dcd.tcl |
Now that you have a dcd file for the protein only, you will use the following trick: for each frame of your pre-equilibrated trajectory, start a NAMD simulation that will run for 0 time steps taking the coordinates from the dcd. This will make NAMD write the energies for every frame on the dcd in the output file.
The configuration file for running NAMD in this setup is ubq-get-energy.conf. The way NAMD loops over all frames in the trajectory, taking the coordinates from the current frame and running for 0 time steps is included in this file:
# runs 0 time steps for each frame in the dcd | |
# opens the dcd file to read the coordinates | |
coorfile open dcd ubq-nvt.dcd | |
set i 0 | |
while { ![coorfile read] } { | |
incr i | |
firsttimestep $i | |
run 0 | |
} | |
coorfile close |
Now you have a log file with the energy of the protein. You need to retrieve the values for the total energy and put them into a file:
namddat TOTAL ubq-get-energy.log |
Now that you have a data file data.dat with the energy of the protein over time, you can calculate the energy fluctuation. For this you need:
Substance | Specific Heat |
Human Body(average) | 3470 |
Wood | 1700 |
Water | 4180 |
Alcohol | 2400 |
Ice | 2100 |
Gold | 130 |
Strawberries | 3890 |