Tutorial on Electron Transfer

In this part of the tutorial you will construct an energy gap function through the molecular dynamics study of an electron transport protein.
The computational demands of this section of the tutorial are rather large compared to what you have studied in previous tutorials. You may want to immediately get started on your first molecular dynamics run of the day and then read on further while your simulation is running.

Starting your simulation

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...forget to keep an eye on your simulation while you read on. } \end{minipage} }$

The system you are currently simulating is the electron transport protein cytochrome c

from the purple bacterium Rhodobacter sphaeroides. You have three main steps to perform in this section:

After starting your simulation, let us first take a closer look at cytochrome c

Interlude: Structure of cytochrome c

$% latex2html id marker 3811 \fbox{ \begin{minipage}{.2\textwidth} \includegrap... ...cytochrome c$_2$\ is shown in Fig.~\ref{fig:cytc2withheme}. } \end{minipage} }$

To explore the structure of cytochrome c

we will use VMD. Start VMD and load the required coordinate and parameter files by typing

in your working directory. This loads VMD with a state file that allows you to view the equilibrated cytochrome c

structure in the beginning of the simulation. (You do not need to type these lines as they are taken from the state file you have just loaded.)
mol new cyt_reduced.psf type psf mol addfile SPH_1CXCequi.pdb type pdb mol representation Cartoon 2.100000 12.000000 5.000000 mol color ColorID 7 mol selection protein mol representation Bonds 0.300000 6.000000 mol color Name mol selection resname HEMP

This should give you a picture that resembles Fig. 7. You may also turn on the water molecules to note that this system is solvated in a water box. Indeed your current simulation has started from an already minimized and equilibrated state with periodic boundary conditions. Let us now take a closer look at the simulation configuration you are currently running.

**Figure 7:** Cytochrome c in the cartoon representation with its heme group highlighted.
$\begin{figure}\begin{center} \par\par\latex{ \includegraphics[scale=0.5]{pictures/cytc2withheme} } \end{center} \end{figure}$

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ... correspond to the {\em reduced} state of cytochrome c$_2$. } \end{minipage} }$

Interlude: a closer look at your configuration files

The NAMD configuration file cytreduced.namd you executed earlier continues your simulation from a restart point.
set inputname cyt_red_init bincoordinates $inputname.rst.coor ExtendedSystem $inputname.rst.xsc binvelocities $inputname.rst.velWithout a restart point you would have had to worry about velocity relaxations and discard the initial part of your already short simulation.

Closer examination of the configuration file will reveal use of periodic boundary conditions. Also of interest is the unusual and ordinarily costly fact that you are writing a DCD frame at each and every time step. outputEnergies 1 dcdfreq 1 This is needed to be able to probe the time dependence of the energy gap function you will compute below.

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...he system for the oxidized and reduced states respectively. } \end{minipage} }$

The output of the first namd run will be written to a DCD file called cytreduced.dcd, which in turn will be read by a second NAMD run described below.

Re-running NAMD to read your previous trajectory

Now you will start a new simulation to read the trajectory created by the first one to compute the interation between two sets of atoms in your simulation. In section 2.1.5 of the NAMD tutorial of past week, you have seen how to compute the specific heat of a protein by a similar technique.

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...member to keep an eye on your simulation while you read on. } \end{minipage} }$

The configuration file you just executed, cytrun2.namd, refers to a new structure file, cyt_deltaQ.psf, which contains the information about the charge differences between the reduced and oxidized forms of the HEME group. Hence, the total electrostatic energy computed by this run, appearing in the ELECT column of the output file, will correspond to the energy gap function defined above.

$\fbox{ \begin{minipage}{.2\textwidth} \includegraphics[width=2.3 cm, height=2.... ...hout modifying the partial charges in the rest of the system.} \end{minipage} }$

$\framebox[\textwidth]{ \begin{minipage}{.2\textwidth} \includegraphics[width=2... ... similar method was employed for specific heat computations. } \end{minipage} }$

From the trajectory to the energy gap function

In order to extract the gap function information, after your second NAMD run is over, type:

This writes the gap function to the file data.dat. Now cut the header (i.e. the first line) of the output data.dat with an editor, so we can read it easily with a script. For example, you may use pico

Delete the first line TS ELECT so that the file consists only of two columns of numbers.

As in the previous part of the tutorial you will use Mathematica to display your results. In the limited time frame of a tutorial session, it is not possible to produce and analyze a long molecular dynamics trajectory. For your convenience, a longer 10 ps version of the energy gap function data file, dataLONG.dat, is provided for you. In order to use Mathematica to extract the energy gap function from the data file. First type

and then select Kernel $\rightarrow$ Evaluation $\rightarrow$ Evaluate Notebook (See Fig. 8.). This will read the data file containing the energy gap function you have created earlier, as well as a longer data file containing the data from a 10 ps simulation. You can now examine your data more closely. In particular, note that the histogram of the energy gap function over time gives a Gaussian distribution. This concludes the last section of this tutorial.

**Figure 8:** Running a *Mathematica* script to study the energy gap function.
$\begin{figure}\begin{center} \par\par\latex{ \includegraphics[scale=0.5]{pictures/mma} } \end{center} \end{figure}$

$% latex2html id marker 3943 \fbox{ \begin{minipage}{.2\textwidth} \includegrap... ...ugh which the trajectory is divided for sampling. \end{list}} \end{minipage} }$

**Figure 9:** The energy gap function for cytochrome c obtained through a MD simulation. (Left: for a 500 fs trajectory; Right: for a 10 ps trajectory)
$\begin{figure}\begin{center} \par\par\latex{ \includegraphics[scale=0.7]{pictures/gapplot} } \end{center} \end{figure}$

**Figure 10:** The histogram of the the energy gap function, $\epsilon (t_j)$ , sampled over a time interval of 10 ps.
$\begin{figure}\begin{center} \par\par\latex{ \includegraphics[scale=0.7]{pictures/gaphistog} } \end{center} \end{figure}$