This notebook contains a tutorial to perform MD simulations of an Octamethyl Silsesquioxane molecule with NAMD and the AMBER force field

me8t8

Table of Contents

Octa-methyl Silsesqioxane
Carry out the simulation
Plot thermodynamic and energetic quantities
Visualize the trajectory
Carry out the simulation in the NVE ensemble
Remove global translations and rotations with AmberTools
Additional tutorials from the NAMD/VMD webpage
View/Edit the contents on an ASCII file

Section

• subsection

Examples of HTML and Markdown syntax</br>

(Top)

Octa-methyl Silsesqioxane (mPOSS)

The mPOSS molecule is composed of a cubic cage where silicon atoms occupy the cube vertices, and oxygen atoms are located in the cube edges (see below figure). Thus, each Si atom has a tetrahedral coordination to three O atoms and one methyl group. Methyl substitution by other chemical species makes POSS molecules highly versatile, with applications as organic solvents, polymer dispersants, catalysts, nanocomposites, diodes, and many other uses. In particular, mPOSS has found application as a coating for carbon fibers and low-dielectric films.

mPOSS molecule composed of Si (yellow), O (red), C (cyan), and H (grey) atoms. Nine different chains of consecutive O-Si-C-H covalent bonds can be constructed for each methyl group, due to the three different oxygen and hydrogen atoms that can be selected at the extremes of the chain.

mPOSS molecule non-vibrational degrees of freedom are restricted to dicrete rotational diffusion of the methyl groups (-CH3). In the AMBER force field the barrier to rotation is described by a dihedral 4-body term.

$V(\phi)=K[1+cos(3\phi)]$
Where $\phi$ is the dihedral angle defined by one of the nine combinations that can be formed with the four linked atoms O, Si, C, and H.

(Top)

Carry out the simulation

We will simulate a single molecule. The necessary files to carry out the simulation are:

• me8t8.pdb (list of atoms and molecule name)
• me8t8.prmtop (force field and topology)
• run0.restart.coor (initial coordinates)
• run1.conf (NAMD run options)

Inspect file run1.conf (there are many ways in Linux to view/edit the contents on an ASCII file) How long is the simulation? How many snapshots, or frames, will be saved? Which ensemble are we running?

1 In the terminal, go to the directory containing the files, then create a subdirectory in your scratch area and copy the files there.

2 Load the NAMD module. This will make available the namd command to actually run the simulation

Now run the simulation. We ask that you limit your simulation to four cores. It should take about 6 minutes

After the simulation is done, it will print out in the terminal the time it took to complete. New files are generated:

• run1.log #thermodynamic quantities are saved here. (ASCII file).
• run1.dcd #snapshots are saved here (binary file).
• run1.restart.coor #coordinates of the system when the simulation finished, in high precision
• run1.restart.vel #velocities of the system when the simulation finished, in high precision
• run1.restart.xsc #geometric boundaries of the system when the simulation finished

One can elongate the simulation by running a second simulation using as starting conformation files files run1.restart.coor, run1.restart.vel,run1.restart.xsc

(Top)

Plot thermodynamic and energetic quantities

We will use the VMD program to visualize results. Thermodynamic quantities are stored in file run1.log and can be plotted with the namdplot plugin. In the terminal, start vmd:

In case the display window shows up off-screen, issue the following command in the terminal to bring the top of the display window into view:

The NAMD Plot plugin is accessible from the Extensions tab:

The interface is pretty self explanatory. Loading of the file and plotting the desired quantities are accessible from the File menu:

You can inspect the different energy terms (BOND, ANGLE, DIHEDRAL, IMPROP, ELEC, VDW) as well as total energy (TOTAL) and temperature (TEMP). Since this is a simulation in "infinite" vacuum, there is no relevant information on volume or pressure. Plotting quantities like total energy and temperature is a miniminal and easy test to check that the simulation did not go wrong.

(Top)

Visualize the trajectory

Next, we will visualize the simulation as a movie. We will also use VMD for this purpose and you can use the current VMD session.
First we need to load the PDB file me8t8.pdb. From the File menu select New Molecule and load the molecule by browsing to the file and then click the Load button.

Although the PDB contains only the names of the atoms and the coordinates of a conformation, VMD is smart enough to figure out which atoms are bonded to each other, thus creating a topology for the molecule and will show the bonds. Once we have topology, we can load the time evolution of coordinates with the trajectory file run1.dcd. In the main window, select molecule "me8t8.pdb" with a left click. The selection will be highlighted. Then, right-click and select Load Data Into Molecule" from the pop up menu.

The Molecule File Browser will appear. Browse to file run1.dcd and then click in Load.

There are 4000 frames to be loaded so it will take about a minute to load them all. Once finished, the Main window will show that there are 4001 frames loaded: 1 for the initial PDB, and 4000 for the trajectory.

The lower portion of the windows feature the progress bar. By sliding the bar you can select any frame. Right now it shows the last frame, which is frame with index 4000. Notice that indexes begin at zero, thus frame number 4001 has index 4000.
You probably have noticed that the molecule has vanished from the OpenGL Display window!

The reason is the following: the coordinates stored in the PDB file run1.pdb are numerically very different than the tipical coordinates in trajectory file run1.dcd. The molecule "makes a big jump" in transitioning from frame with index=0 (the PDB file) to frame with index=1 (the first frame of the trajectory). It's like if the molecule suddenly jumped behind our back and we could not see it any more. To fix this, let's rewind the movie by pressing the "Jump to the beginning" button:

We are back to the conformation of the PDB file (index=0). Now let's advance to the first frame of the trajectory by pressing the "Step forward" button:

As expected, the molecule has vanished from the OpenGL Display Window. We can change our perspective by resetting the view from the Display menu:

Now we are looking at the first frame of the trajectory and are ready to watch the movie. Click in the "play forward" button (red circle in the picture). Also, you can modify how fast the movie plays and whether you want to skip frames with the tools highlighted in the green ellipse in the picture below:

The molecule does not stay put, it rotates as a whole and its center of mass translates. This is typical of molecules subject to Brownian motion, but since the molecule is surrounded only by vacuum, why is undergoing Brownian motion?
The answer is the molecule is surrounded by "ghost particles" due to the heat bath we imposed to maintain a constant temperature. The section in configuration file run1.conf:

Constant Temperature Control
temperature   $temperature
langevin            on    (do langevin dynamics)
langevinDamping     5.0   (damping coefficient (gamma) of 1/ps (5.0 for equilibrium runs))
langevinHydrogen    off   (don't couple langevin bath to hydrogens)
langevinTemp  $temperature (target temperature)

prescribes Langevin dynamics where random forces prevent conservation of total angular momentum (hence the rotational diffusion) and total momentum (hence the diffusion of the center of mass).

(Top)

Carry out the simulation in the NVE ensemble

As a comparison, let's run the simulation at constant energy (NVE ensemble). Think for a minute how would you modify configuration file run1.conf for such simulation.

The appropriate run1.conf file is located in your $HOME/SHUG2015/molecular_dynamics/examples/me8t8/NVE/ directory. The only difference with the run at constant temperature is the commenting out of the lines corresponding to the Constant Temperature Control section. Notice we have left out a line uncommented:

temperature   $temperature

This line will assign to each atom a velocity vector with random orientation and modulus sample from a Maxwell distribution of velocities at temperature $temperature. In order to start the simulation, we need initial velocities in addition to initial coordinates.
Create another directory in your scratch area, for instance:

Run the simulation and use the NamdPlot plugin in VMD to check the evolution of the total energy (quantity named "TOTAL") versus time.

• Is it constant? If not, why?
• How do the fluctuations of the total energy in the NVE simulation compare to the fluctuations of the total energy in the previous NVT simulation? What about the fluctuations of the temperature?

(Top)

Remove global translations and rotations with AmberTools

As we saw in the simulation at constant temperature, the molecule undergoes Brownian motion. It's diffusion coefficient depends on parameter LangevingDamping we used in the Constant Temperature Control section of file run1.conf.
If we want to look at the internal dynamics of the molecule, where it is expected that the particular heat bath has much less influence, we have to remove these global rotations and translations. We will use AmberTools, the analysis tools provided by the AMBER developers. In particular, we will use the terminal command cpptraj that allows transformations of a molecular dynamics trajectory as well as computation of many quantities. For details, read the cpptraj section in the Amber14 guide.
First, let's open the file containing the commands that will be passed to cpptraj, $HOME/SHUG2015/molecular_dynamics/examples/me8t8/rms2first.cpptraj.

trajin run1.dcd  
rms first time 1.0 out rms2first.dat  
trajout run1_rms2first.dcd dcd  
quit

The first line commands cpptraj to load the trajectory onto memory. The second line will perform an RMS superposition of each frame in the trajectory to the first frame. We also indicate we want to output the RMS values onto file rms2first.dat, and indicate that consecutive frames in the trajectory are separated by 1ps. The third lane will save the trajectory in file run1_rms2first.dcd. This is the trajectory without global translations and rotations.
To run cpptraj, we have to load the amber module:

cpptraj needs also the topology file me8t8.prmtop, although we could have used the PDB file me8t8.pdb. After we have entered these lines in the terminal, we can plot the RMS values:

The RMS values increase rapidly from zero (the RMS of the first frame to itself must be zero) to a value $<RMS>=$1.0A.

Next, use VMD to visualize the new trajectory. First load the trajectory:

Rewind the trajectory to frame with index=1 (review the instructions in section Visualize the trajectory ) and then play the movie. What are the motions of the molecule now?

(Top)

Additional tutorials from the NAMD/VMD webpage

There are many tutorials to master both NAMD and VMD from the training webpage.

(Top)

View/Edit the contents on an ASCII file

For those not familiar with the Linux operating system, we enumerate here a few ways to view and/or edit an ASCII file:

1 File is shown in the terminal:

2 File is shown in a separate window

(Top)

Section

(Top)

Subsection

(Top)

Markdown Syntax Examples

local link: link</br> remote link: http://ambermd.org font face="courier new"
$$S_{model}(Q,E)=A(Q)\cdot S_{elastic}(E) + B(Q)\cdot S_{simulation}(Q,E)\otimes S_{elastic}(E) + C(Q)+D(Q)\cdot E$$

 Quoted text 

image caption

some text