next up previous
Next: Non-equilibrium Up: Analysis Previous: Analysis

Subsections


Equilibrium

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.


RMSD for individual residues

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...ith each residue colored according to its average RMSD value.}
\end{minipage} }

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 User field for each residue.

\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
... contain pertinent PDB information!) like the User field can.}
\end{minipage} }

In the previous unit, you had VMD open with the equilibration trajectory of the ubiquitin periodic system loaded. You will repeat this process:

* 1
Load the structure file. First, open the Molecule File Browser from the File $\rightarrow$ New Molecule... menu item. Browse for and Load the file ubq_wb.psf in the directory common.

* 2
The menu at the top of the Molecule File Browser window should now show 2: ubq_wb.psf. This ensures that the next file that you load will be added to that molecule (molecule ID 2). Now, browse for the file ubq_wb_eq.dcd located in the directory 1-3-box/ and click on Load again. You can close the Molecule File Browser window.

* 3
Go to the TkCon window in VMD, and navigate to the 2-1-rmsd/ directory by typing cd ../2-1-rmsd/

* 4
The script you will use is called residue_rmsd.tcl. In the TkCon window, type:

source residue_rmsd.tcl
 

This command itself does not actually perform any calculations. It reads 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 <mol> <sel_resid>
 


where <mol> is the molecule in VMD that you select (if not specified, defaults to the top molecule), and <sel_resid> is a list of the residue numbers in that selection. The procedure employs the equations for RMSD shown in Unit 1.

* 5
For this example, you will select all the residues in the protein. The list of residue numbers can be obtained by typing in the TkCon window:

set sel_resid [[atomselect top "protein and alpha"] get resid]
 

The above command will get the residue numbers of all the $\alpha $-carbons in the protein. Since there is just one and only one $\alpha $-carbon per residue, it is a good option. The command will create a list of residue numbers in the variable $sel_resid.

* 6
Call the procedure to calculate the RMSD values of all atoms in the newly created selection:

rmsd_residue_over_time top $sel_resid
 

You should see the molecule wiggle as each frame gets aligned to the initial structure (This may happen very fast in some machines). At the end of the calculation, you will get a list of the average RMSD per residue. This data is also printed to the file residue_rmsd.dat.

The procedure also sets the value of the User field 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 User value.

* 7
Open the Representations window using the Graphics $\rightarrow$ Representations... menu item.

* 8
In the Atom Selection window, type protein. Choose Tube as Drawing Method, and in the Coloring Method drop-down menu, choose submenu Trajectory, then User, and again User. Now, click on the Trajectory tab, and in the Color Scale Data Range, type 0.40 and 1.00. Click the Set button.

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 somewhat counterintuitive, so you will change the color scale.

* 9
Choose the Graphics $\rightarrow$ Colors$\ldots$ menu item. Choose the Color Scale tab. In the Method pull-down menu, choose BGR. Your figure should now look similar to Fig. 9.

Figure 9: Ubiquitin colored by the average RMSD per residue. Red denotes more mobile residues, and blue denotes residues which moved less during equilibration.
Image tut_residue_rmsd

* 10
Open Excel by clicking Start $\rightarrow$ Programs and finding your version of Excel.

* 11
In Excel, click File $\rightarrow$ Open$\ldots$, navigate to the directory 2-1-rmsd/, and choose the file residue_rmsd.dat. Make sure the dropdown menu Files of type is set to All Files.

* 12
When the Text Import Wizard appears, choose Delimited for your file type and click Finish. Your data is tab delimited, so it should appear in Columns A (residue number) and B (RMSD).

* 13
Make a plot of the data by clicking Insert $\rightarrow$ Chart$\ldots$, choosing XY (Scatter) as your Chart Type, selecting the last Chart sub-type in the lower-right, and clicking Finish. Excel should plot Column A on the x-axis and B on the y-axis, and the plot of RMSD data for each residue will appear.

Take a look at the RMSD distribution. You can see regions where a set of residues shows less mobility. Compare with the structural feature of your protein. You will find a correlation between the location of these regions and secondary structure features like $\alpha $ helices or $\beta $ sheets.

* 14
Save the files you created in Excel if you wish, and close Excel.

* 15
Delete the current displayed molecule in VMD using the VMD Main window by clicking on ubq_wb.psf in the window, and selecting Molecule $\rightarrow$ Delete Molecule.

Maxwell-Boltzmann Energy Distribution

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...esponds to the Maxwell distribution for a given temperature. }
\end{minipage} }

For this exercise, you need a file that contains the velocities of every atom for one instant of the simulation. You will use a file corresponding to the last frame of an equilibration of ubiquitin. The simulation that generated this file was run exactly as the one you performed in directory 1-3-box/, except that the time step used here was 1fs. You can compare the configuration file used for that simulation to the one in the example-output/ of this directory. As discussed in the first unit, one should use rigidBonds for water in any production simulation since water molecules have been parametrized as rigid molecules. But for the illustration of Maxwell-Boltzmann energy distribution, the effect of force field can be negligible. To be simple, we are not using the option rigidBonds. Here, we provide you with the file ubq_wb_eq_1fs.restart.vel. Alternatively, you can use any velocity file, as long as constraints were not used in the simulation.

1
You will now load the structure file. First, open the Molecule File Browser from the File $\rightarrow$ New Molecule... menu item. Browse for and Load the file ubq_wb.psf in the directory common.

2
The menu at the top of the window should now show 3: ubq_wb.psf. This ensures that the next file that you load will be added to that molecule (molecule ID 3). Now, browse for the file ubq_wb_eq_1fs.restart.vel in the directory 2-2-maxwell/. This time you need to select the type of file. In the Determine file type: pull-down menu, choose NAMD Binary Coordinates, and click on Load again. Close the Molecule File Browser.

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.

3
In the VMD TkCon window, navigate to the 2-2-maxwell/ directory by typing cd ../2-2-maxwell/

4
First, you will create a selection with all the atoms in the system. Do this by typing in the TkCon window:

set all [atomselect top all]
 

5
Now, you will open a file energy.dat for writing:

set fil [open energy.dat w]
 

6
For each atom in the system, calculate the kinetic energy $\epsilon_k=\frac{1}{2}m v^2$, and write it to a file

foreach m [$all get mass] v [$all get {x y z}] {  
   puts $fil [expr 0.5 * $m * [vecdot $v $v] ]  
}
 

7
Close the file:

close $fil
 

8
Outside VMD, use WordPad to take a look at your new energy.dat file. When you open the file, make sure the dropdown menu Files of type is set to All Files. It should look something like the list below (although the values will not necessarily be the same):

1.26764386127  
0.189726622306  
0.268275447408  
...
 

9
Quit WordPad.

* 10
If you chose not to go through the steps above, there is a VMD script that reproduces them and will create the file energy.dat. In the TkCon Window, type:

source get_energy.tcl
 

This will source the script get_energy.tcl, and produce the energy.dat file.

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 Excel. If you are more familiar with another tool, you may use it instead.

* 11
Quit VMD.

* 12
In Excel, click File $\rightarrow$ Open$\ldots$, navigate to the directory 2-2-maxwell/, and choose the file energy.dat. Make sure the dropdown menu Files of type is set to All Files.

* 13
When the Text Import Wizard appears, choose Delimited for your file type and click Finish. The energy data will appear in Column A.

We will construct a histogram plot of the energy values and fit a Maxwell-Boltzmann distribution to it to find the average temperature of your simulations.

* 14
Construct a column of bin labels for the histogram. In cell B1, enter the number 0, and in cell B2, enter the formula =B1+0.1 and hit Enter. Click on cell B2, the click and hold its lower right corner and drag down to row 71 and release. You should have a list of values from 0 to 7 filling Column B. These are our histogram bins.

* 15
We will now use the Data Analysis tool to create the histogram data. Click Tools $\rightarrow$ Data Analysis$\ldots$. If that option is not present, you will need to add in the Data Analysis tool by clicking Tools $\rightarrow$ Add Ins$\ldots$ and selecting Analysis ToolPak from the list.

* 16
Select Histogram from the Data Analysis window and click OK. Type in the following cell ranges and click OK:

Input Range: $A$1:$A$7051  
Bin Range: $B$1:$B$71
 

A new sheet will be created with your energy bins and the number of energy measurements which occur within the given bin. We could plot these values, but we would prefer to plot the normalized frequency instead.

To obtain the normalized frequencies, we will obtain the area under our histogram and divide each frequency by it.

* 17
In the new sheet, in cell C2, type the formula =B2*0.1. Then click and hold the lower-right corner of cell C2 and drag down to row 72 to fill in the area values for each bin.

* 18
In cell D2, type the formula =sum(C2:C72). This is the area under your histogram.

* 19
Create the normalized frequencies by typing the following formula into cell E2: =B2/D$2. Extend this input down to row 72 as you did for Column C before.

* 20
Make a histogram plot of the data by clicking Insert $\rightarrow$ Chart$\ldots$, choosing Column as your Chart Type, selecting the first Chart sub-type in the upper-left, and clicking Next.

* 21
In Step 2, click the Series tab and remove all existing series. Then add your own by clicking the Add button. Click the red button in the Values: field, select every data value in Column E (Rows 2-72), and click the red button to complete. Do the same for the Category (X) axis labels: using the bins in Column A. Click Finish.

* 22
Modify the chart in anyway you like. We suggest right-clicking on the histogram bars, selecting Format Data Series$\ldots$, clicking the Options tab, and setting the Gap width: to 0.

Figure 10: Maxwell-Boltzmann distribution for kinetic energies.
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_unit02_mboltzwin}
}
\end{center}
\end{figure}

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!)

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...\exp\left(-\frac{\epsilon_k}{k_{B}T}\right)}
\end{equation} }
\end{minipage} }

* 23
In cell F2, type the formula::

= (2 / SQRT(pi() * G$2$\wedge$3)) * SQRT(A2) * EXP(-A2 / G$2)
 

Enter the value 0.5 in cell G2, then extend the formula in F2 down to cell F72. This will create the Maxwell-Boltzmann distribution which we will fit to our raw data. The parameter in cell G2 corresponds to $k_{B}T$ in units of kcal/mol. (These are the energy units that NAMD uses.)

* 24
In the histogram chart, right-click on the histogram bars, and select Source Data$\ldots$. Click the Series tab and add a series with Values: in Column F, Rows 2-72. Click OK.

* 25
Right-click the red histogram bars you just created, select Chart Type$\ldots$, choose Line as your Chart Type, and choose the first Chart sub-type in the upper-left. Click OK.

You may want to change the line weight of your Maxwell-Boltzmann fit by right-clicking the line in your histogram and selecting Format Data Series$\ldots$.

We have the Maxwell-Boltzmann expression plotted, but the fit is not optimal. Indeed, if we change the value for $k_{B}T$ in cell G2, we will see the fit change. We will now determine the value of $k_{B}T$ that yields the optimal fit.

* 26
In cell H2, type the formula =(F2-E2)$\wedge$2. Then click and hold the lower-right corner of cell H2 and drag down to row 72 to fill in the error values for each bin.

* 27
Type =sum(H2:H72) into cell I2. This is the total square error for our fit, and it is what we seek to minimize. We will do so using the Solver tool in Excel.

* 28
Click Tools $\rightarrow$ Solver$\ldots$. If that option is not present, you will need to add in the Solver tool by clicking Tools $\rightarrow$ Add Ins$\ldots$ and selecting Solver Add-in from the list.

* 29
In the Solver Parameters window, set the following cell ranges and click Solve:

Set Target Cell: $I$2  
Equal To: Min  
By Changing Cells: $G$2
 

Solver will find the value for $k_{B}T$ (cell G2) which optimally fits the Maxwell-Boltzmann distribution to our raw energy data. Your plot should look like Fig. 10. Obtain the temperature $T$ for this distribution with $k_{B} = 0.00198657 \frac{kcal}{mol K}$. Your result is hopefully close to $310 $K!

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 $\pm10K$ even in case of perfect sampling.

* 30
Save the files you created in Excel if you wish, and close Excel.

Energies

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...different internal energies) over the course of a simulation.}
\end{minipage} }

You will look at the kinetic energy and all internal energies: bonded (bonds, angles and dihedrals) and non-bonded (electrostatic, van der Waals) over the course of the water box simulation which you ran in the previous unit. This section will introduce you to a tcl script called namdstats.tcl which is very useful for calculating average energy data for a simulation or extracting energy data over time.

* 1
In a Terminal, go to the directory 2-3-energies/

We will investigate the simulation you ran in directory ../1-3-box/. Namdstats.tcl will search all the information in the NAMD output file and extract the energy data for you. You should investigate it with your text editor to familiarize yourself with it.

* 2
Open VMD by clicking Start $\rightarrow$ Programs $\rightarrow$ VMD.

* 3
In the VMD TkCon window, navigate to the 2-3-energies/ using the cd command if you are not already there.

* 4
To use the script, type

source namdstats.tcl  
data_avg ../1-3-box/ubq_wb_eq.log 101 last
 

The first line will source the script itself, and the second will call a procedure which will calculate the average for all output variables in the log file from the first logged timestep after 101 to the end of the simulation.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...data stream>} is one of the energy values that NAMD outputs.
}
\end{minipage} }

* 5
The output of namdstats.tcl should appear similar to that in Table 1. Note that the standard units used by NAMD are Å for length, kcal/mol for energy, Kelvin for temperature, and bar for pressure.


BOND: 232.561128
ANGLE: 690.805988
DIHED: 397.649676
IMPRP: 41.250136
ELECT: -23674.863868
VDW: 1541.108344
BOUNDARY: 0.0
MISC: 0.0
KINETIC: 4579.147764
TOTAL: -16192.340832
TEMP: 313.429264
TOTAL2: -16173.821268
TOTAL3: -16174.663916
TEMPAVG: 313.140436
PRESSURE: -47.513352
GPRESSURE: -19.105588
VOLUME: 71087.006196
PRESSAVG: 3.340104
GPRESSAVG: 3.4207
% latex2html id marker 20192
$\textstyle \parbox{0.75\textwidth}{\caption{Typic...
...r length, kcal/mol for energy, Kelvin for temperature, and bar for pressure.
}}$

The temperature of your system should be close to 310 K, the temperature we specified for the water box simulation. It will not be exact because temperature fluctuations occur due to the finite size of the system. We will investigate these fluctuations in the next section.

Temperature distribution

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Ob...
...in an NVE ensemble, and analyze the temperature distribution.}
\end{minipage} }

* 1
In a Terminal window, go to the directory 2-4-temp/.

For this exercise, you need a simulation that is long enough to sample the quantity of interest very well. Here, we need to sample 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 example-output/ubq-nve.log. However, if you choose to run your own simulation, you can find all the necessary files and instructions below. We recommend that you read the whole section even if you don't perform the simulations.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...displaymath}where $N$\ is the number of atoms in the system.
}
\end{minipage} }

The simulation is performed in the microcanonical ensemble (NVE, i.e., constant number, volume, and energy). The configuration file provided to start the NAMD run is called ubq-nve.conf. The main features of this configuration file are:

2
Run the Molecular Dynamics simulation

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

* 3
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/.

copy example-output\ubq-nve.log .
 

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

* 4
In the VMD TkCon window, navigate to the 2-4-temp/ directory by typing cd ../2-4-temp/

* 5
In order to obtain the data for the temperature from the log file we will again use the script namdstats.tcl, which we have already sourced. Do this by typing:

data_time TEMP ubq-nve.log
 

Each timestep and the corresponding temperature of the system is now in the file TEMP.dat. The script will include all timesteps and it may take some time to finish. Note there is no minimization of the system since this was already performed in the initial simulation.

In order to plot the data, you can use Excel.

* 6
Open Excel.

7
Click File $\rightarrow$ Open$\ldots$, navigate to the directory 2-4-temp/, and choose the file TEMP.dat. Make sure the dropdown menu Files of type is set to All Files.

8
When the Text Import Wizard appears, choose Delimited for your file type and click Next. Under the category Delimiters, click Tab. The Data Preview window should show your timestep and temperature data separated into two columns. Click Finish. The timestep and temperature values will appear in Columns A and B, respectively.

We will construct a histogram plot of the temperature values and fit a Gaussian distribution to it.

* 9
Construct a column of bin labels for the histogram. In cell C1, enter the number 300, and in cell C2, enter the formula =C1+0.25 and hit Enter. Click on cell C2, then click and hold its lower right corner and drag down to row 121 and release. You should have a list of temperature values from 300 to 330 filling Column C. These are our histogram bins.

* 10
We will now use the Data Analysis tool to create the histogram data. Click Tools $\rightarrow$ Data Analysis$\ldots$. If that option is not present, you will need to add in the Data Analysis tool by clicking Tools $\rightarrow$ Add Ins$\ldots$ and selecting Analysis ToolPak from the list.

* 11
Select Histogram from the Data Analysis window and click OK. Type in the following cell ranges and click OK:

Input Range: $B$1:$B$10000  
Bin Range: $C$1:$C$121
 

A new sheet will be created with your temperature bins and the number of temperature measurements which occur within the given bin. We will plot these values.

* 12
Make a histogram plot of the data by clicking Insert $\rightarrow$ Chart$\ldots$, choosing Column as your Chart Type, selecting the first Chart sub-type in the upper-left, and clicking Next.

* 13
In Step 2, click the Series tab. Remove the first series, and then click the Add button to add a new series. Click the red button in the Values: field, select every data value in Column B (Rows 2-122), and click the red button to complete. Do the same for the Category (X) axis labels: using the bins in Column A. Click Finish.

* 14
Modify the chart in anyway you like. We suggest right-clicking on the histogram bars, selecting Format Data Series$\ldots$, clicking the Options tab, and setting the Gap width: to 0.

Does this distribution look any familiar? Your distribution should look like a Gaussian distribution.

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... broad. For $T=300K$\ and $N=1000$\ one has $\sigma{\sim}8K$.}
\end{minipage} }

* 15
You will now fit your distribution to a normal distribution. We will do this in a similar manner as for the Maxwell-Boltzmann distribution of energies.

* 16
In cell C2, type the formula::

= D$2 * EXP(-((A2-D$3)$\wedge$2)/D$4)
 

The values in Column D are parameters for the Gaussian distribution. Enter initial guesses of 400, 315, and 100 in cells D2, D3, and D4, respectively. Then extend the formula in C2 down to cell C122. This will create the Gaussian distribution which we will fit to our raw data. The parameters correspond to the normalization constant, the average temperature, and $\sigma^2$.

* 17
In the histogram chart, right-click on the histogram bars, and select Source Data$\ldots$. Click the Series tab and add a new series with Values: in Column C, Rows 2-122. Click OK.

* 18
Right-click the red histogram bars you just created, select Chart Type$\ldots$, choose Line as your Chart Type, and choose the first Chart sub-type in the upper-left. Click OK.

You may want to change the line weight of your Gaussian fit by right-clicking the line in your histogram and selecting Format Data Series$\ldots$.

We have the Gaussian plotted, but the fit is not optimal. If we change any of the values for parameters Column D, we will see the fit change. We will now determine the optimal parameters.

* 19
In cell E2, type the formula =(C2-B2)$\wedge$2. Then click and hold the lower-right corner of cell E2 and drag down to row 122 to fill in the error values for each bin.

* 20
Type =sum(E2:E122) into cell F2. This is the total square error for our fit, and it is what we seek to minimize using the Solver tool.

* 21
Click Tools $\rightarrow$ Solver$\ldots$. If that option is not present, you will need to add in the Solver tool by clicking Tools $\rightarrow$ Add Ins$\ldots$ and selecting Solver Add-in from the list.

* 22
In the Solver Parameters window, set the following cell ranges and click Solve:

Set Target Cell: $F$2  
Equal To: Min  
By Changing Cells: $D$2:$D$4
 

Solver will find the parameter values which optimally fit the Gaussian distribution to our raw temperature data.

Figure 11: Fluctuation of Temperature in a Finite Sample
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_unit02_tempwin}
}
\end{center}
\end{figure}

* 23
Compare the average value of the temperature obtained by this method with the one you would obtain using namdstats by typing in the VMD TkCon window:

data_avg ubq-nve.log
 

* 24
Now, look at the deviation $\sigma^2=2T^2/3N$ which can be calculated from your data. Note how it is dependent on the size of the system!

Take home message: You may note from this example that single proteins, on one hand, behave like infinite thermodynamic ensembles, reproducing the respective mean temperature. On the other hand, they show signs of their finiteness in the form of 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!

* 25
Save the files you created in Excel if you wish, and close Excel.

Specific Heat

\framebox[\textwidth]{
\begin{minipage}[r]{0.9\textwidth}
\noindent{\textbf{Objective:} Find the specific heat of ubiquitin.}
\end{minipage} }

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.)

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
...^2}\rangle - \langle{E}\rangle^2}{k_{B}T^2}
\end{displaymath}}
\end{minipage} }

For this exercise, you need to carry out a long simulation in the NVT (canonical) ensemble that samples sufficiently the averages $\langle{E^2}\rangle$ and $\langle{E}\rangle$ that arise in the definition of $c_V$. 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. You may wish to work only through the steps marked with *, that correspond to the analysis of the data. Nevertheless, we encourage you to read through the whole section since it will provide you with tools for future analysis.

* 1
In a Terminal windows, change to the directory 2-5-spec-heat/.

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.

2
Run the Molecular Dynamics simulation

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

\fbox{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2.3 cm, height=2....
... use, as normally you won't need this massive kind of output.}
\end{minipage} }

* 3
If you do not run the simulation yourself, you will need to obtain the trajectory that has been provided. Copy the simulation output from the directory example-output/ by typing in a Terminal window:

copy example-output\ubq-nvt* .
 

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

You now have an equilibration run of ubiquitin in the NVT ensemble, and you want to look at the fluctuation in the energy of ubiquitin. 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, however 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, rather than for the entire system, as is reported in the log file. One way to do this is to run a separate simulation with only the part of the system you are interested in, but this may be unrealistic from a biological standpoint. Fortunately, VMD has a plugin which will call NAMD and allow you to calculate energies for any part of your system that you specify.

* 4
In the VMD TkCon window, navigate to the 2-5-spec-heat/ directory:

cd ../2-5-spec-heat/
 

* 5
If you did not run the simulation, and copied output from the example-output directory, please skip this step. If you ran the simulation yourself, perform this step and skip the next one.
Load the psf file and trajectory for your system.

mol new ../common/ubq_ws.psf  
mol addfile ubq-nvt.dcd
 

* 6
Load the psf file and trajectory for the example-output system.

mol new ../common/example-output/ubq_ws.psf  
mol addfile ubq-nvt.dcd
 

* 7
In the VMD Main window, open the NAMD Energy plugin by clicking Extensions $\rightarrow$ Analysis $\rightarrow$ NAMD Energy.

Figure 12: NAMD Energy GUI (Graphical User Interface)
\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.5]{pictures/tut_unit02_namdenergy}
}
\end{center}
\end{figure}

* 8
In the NAMDEnergy window, select the following inputs: All other parameters are chosen by default and are consistent with the specifications of your simulation. For more information on NAMD Energy's options, click Help $\rightarrow$ Help...

* 9
Click Run NAMDEnergy to perform the analysis. Output will be written in the vmd console window. When you are finished, close the NAMDEnergy window. (Note that NAMD Energy may give you an error message if it cannot find NAMD installed. You may need to point it to the location of NAMD on your computer.) This analysis may take several minutes to be completed.

Now you have a file called myenergy.dat which contains the energy of only the ubiquitin protein atoms in your system. In order to calculate the specific heat of the molecule, we must derive both the average and squared-average total energy of the molecule from this data. We will use a tcl script to do so.

* 10
In the VMD TkCon window, source the script to calculate the average total energy of ubiquitin for your simulation:

source average.tcl
 

The script will search your output file for the total energy data for the system, calculate both the average and squared-average total energy over all timesteps, and output the values for you. You should become familiar with the script to fully utilize the power of Tcl scripting in analyzing NAMD output. This script will teach you how to read input files, extract data from them, and perform calculations on that data.

* 11
Calculate the specific heat for ubiquitin using the equation provided in the Energy fluctuations and specific heat box at the beginning of this section. Consider $k_{B} = 0.00198657 \frac{kcal}{mol K}$ and $T=310$K. What is your result? Try to convert your results to units of specific heat used everyday, i.e., $J/(kg \symbol{23}C)$ using the conversion factor $1~J= 1.43846\times10^{20}~kcal/mol$. Note that this units correspond to the specific heat per unit mass. You calculated the fluctuation in the energy of the whole protein, so you must take into account the mass of your protein by dividing your obtained value by $m=1.4219\times10^{-23}kg$. Compare with some specific heats that are given to you in Table 2.


Substance Specific Heat
  $J/(kg \symbol{23}C)$
Human Body(average) 3470
Wood 1700
Water 4180
Alcohol 2400
Ice 2100
Gold 130
Strawberries 3890
% latex2html id marker 20328
$\textstyle \parbox{0.75\textwidth}{\caption{ Spec...
...home/sh\_table.html and http://www.eng.auburn.edu/users/wfgale/usda\_course
}}$

* 12
Delete the current displayed molecule in VMD using the VMD Main window by clicking on ubq_ws.psf in the window, and selecting Molecule $\rightarrow$ Delete Molecule.


next up previous
Next: Non-equilibrium Up: Analysis Previous: Analysis
namd@ks.uiuc.edu