<center>https://www.youtube.com/watch?v=qrvJkm_ucP4</center><center><sub>A Simulation Which I did in past: Nogo66 protein interating with DMPC membrane in water(not shown for clarity)</sub></center>



Hi guys! After a long gap, @dexterdev is back with the fourth part of **Classical Molecular Dynamics (CMD) series**. Today let us set up a simulation system in a software called NAMD. For those who missed my initial posts related to CMD, please see it in links below:
[Part 1 (The Fundamentals)](https://steemit.com/steemstem/@dexterdev/classical-molecular-dynamics-series-part-1-the-fundamentals)
[Part 2 (The Force Field!)](https://steemit.com/steemstem/@dexterdev/classical-molecular-dynamics-series-part-2-the-force-field)
[Part 3 (Solving the molecular dynamics equation)](https://steemit.com/steemstem/@dexterdev/classical-molecular-dynamics-series-part-3-solving-the-molecular-dynamics-equation)

One thing to be careful about this article is that this is a ***HOW TO tutorial***, rather than explaining concepts. So if you guys are confused about certain technical terms etc, don't worry. I will be explaining them in the future parts. Also, you can ask me doubts in comments, and I will try answering those. :)

## What is NAMD?

[NAMD](http://www.ks.uiuc.edu/Research/namd/) is a Molecular Dynamics software developed from ***[Theoretical and Computational Biophysics](http://www.ks.uiuc.edu/)*** group from University of Illinois at Urbana-Champaign. 
Quoting from their website:
> NAMD, recipient of a 2002 Gordon Bell Award, is a parallel molecular dynamics code designed for high-performance simulation of large biomolecular systems.

This group has also designed a software called VMD to visualize the simulations etc.

Firstly, you may want to download NAMD and VMD from this webpage: http://www.ks.uiuc.edu/Development/. Mine is a 64-bit Linux machine. Yours may be a different one. So be careful in choosing your selection. You may also need to register with the website to download these files. But they are all free software. 

*[ONE THING TO BE CAREFUL: Depending on the size of the molecule and length of the simulation, most of the simulation cases are only possible on high-performance computers like GPU machines/parallel processor machines with multiple cores. Like I use 128 or 256 cores for a job in our Institute cluster facility. So in practice, you can only check if your simulation is running or not by using your desktop machine/laptop.]*

## Let us do a simulation first!

I am a fan of learning by doing. So let us begin straight away. The example system which we will set up today is that of a phospholipid membrane (specifically DMPC one) in water. To know what DMPC membrane is, please see one of my previous articles:
[Phospholipids: The unsung heroes in biomolecules!](https://steemit.com/steemstem/@dexterdev/phospholipids-the-unsung-heroes-in-biomoleucles)

So to start the simulation of molecules, we need initial structures, right? We are going to initiate animation in the dead structures based on the laws of physics(basically Newton's second law). The [force-field](https://steemit.com/steemstem/@dexterdev/classical-molecular-dynamics-series-part-2-the-force-field) which we are going to use are [CHARMM](http://mackerell.umaryland.edu/charmm_ff.shtml)36 version.

### Intial structure

For this example tutorial, we are considering a phospholipid membrane (DMPC) in water. So the initial structures are in the forms of PDB and PSF formats. I will explain what those are. PDB structure basically contains the coordinate information of atoms in the molecule. 

So we need structure. In the case of bilayer lipid membranes, we can go to [CHARMM-gui membrane builder](http://www.charmm-gui.org/?doc=input/membrane) website. Let us select our options there:

- Check the **Membrane only system** option, first.
- And click on it. You will be directed [here](http://www.charmm-gui.org/?doc=input/membrane_only&step=1). Please see the below options for screenshots. You should choose **Heterogeneous Lipid** option first, then **Hydration number** as [26](http://www1.lsbu.ac.uk/php-cgiwrap/water/pfp.php3?page=http://www1.lsbu.ac.uk/water/membrane_hydration.html), **Length of XY based on**: `Ratios of lipid components`, **Length of X and Y**: `90`(angstrom units=10<sup>-10</sup> metres, you can choose your length instead of 90) and finally select the DMPC lipid from PC option with 1:1 ratio(which will incorporate 127 upper and 127 lower lipids).
<center>![](https://i.imgur.com/XKNcs6T.png)
</center><center>A screenshot from CHARMM-gui website showing the options to choose</center>

<center>![](https://i.imgur.com/dszqn4Z.png)
</center><center>Screenshot from CHARMM-gui website showing the next set of options to choose</center>

- Click the next step. I am changing the IONS option to 0.15 Molar/liter NaCl from KCl. But it doesn't matter much. You can see the xyz dimensions are approximately 90 x 90 x 58 Angstrom<sup>3</sup>. Click next again. (CHARMM-gui will add water and ions on top and bottom of lipid bilayer. )

- And download the output files in .tgz format. You can see a bunch of files in it. There is also PDB file in it. Let us open it using a text editor and VMD software. Visualize the step4_lipid.pdb (which you got from .tgz file after the last step) output in VMD. You can see something similar to below(DMPC lipid bilayer):

<center>![](https://i.imgur.com/d9liCh8.png)</center><center>VMD visualization of PDB file.Type `vmd step4_lipid.pdb` in terminal</center>

<center>![](https://i.imgur.com/cvOqdlW.png)
</center><center>For more information regarding PDB format see [here](https://en.wikipedia.org/wiki/Protein_Data_Bank_(file_format))</center> 

- Choose the options as in below snapshot. We will go through the concept of ensemble etc later. For now, we are going to choose NPT ensemble (Constant number of particles, pressure, and temperature). Forcefield is selected as CHARMM36m(m for modified). Temperature T=303.15.

<center>![](https://i.imgur.com/NUKuUmc.png)</center>

- ***A note about other file types:*** PDB file has the initial coordinate information. PSF file will have the atom's charge, mass, and structure related information etc. Besides PSF and PDB files, we have following files:

    - topology file(.rtf file) in order to create PSF file from PDB file. (this step is not needed for this particular example since CHARMM-gui creates PSF file for you.)
    - parameter file(.prm file) in order to specify parameter values related to bonds, angles, dihedrals, improper dihedrals and Lenard-Jones non-bonding interactions. (You may see my previous MD articles to know what I meant by dihedrals, Lenard-Jones potential etc)
    - configuration file where we will input all necessary arguments for the simulation including above-mentioned files too.

So I want to tell you guys one more thing. If the above file types and other details seem daunting for you, let me assure you that we will come to each of these things again in future. Let me remind you that the focus of this article is just to get start doing simulations.

<center>![](https://i.imgur.com/xkWGrZx.png)</center>

### A typical CONFIGURATION file

```
##############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# DMPC patch equilibration (intitial structure from CHARMM-gui)


#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

# structure and coordinate files are PDB and PSF file generated by CHARMM-gui
structure          ./step5_assembly.psf
coordinates        ./step5_assembly.pdb

set temperature    303.15
set outputname    DMPC_patch_NPT
#DCD and other files output name


firsttimestep       0


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input (CHARMM force-field ON)
paraTypeCharmm      on
parameters           ./par_water_ions.prm
parameters           ./par_all36_lipid.prm

temperature         $temperature

# Force-Field Parameters (will discuss in future)
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.
switching           on
switchdist          10.
pairlistdist        15.0
margin               8.0


# Integrator Parameters
# fullElectFrequency*timestep <=4.0
# stable time steps:
#       with rigidBonds=all: electro:6fs;short-rangeNB:2fs;bonded:2fs
#       otherwise          :         2fs               2fs        1fs
timestep            2.0  ;# 2fs/step (our time step is 2 femto second)
rigidBonds         all  ;# Needed for 2fs steps
useSettle           on
nonbondedFreq       1
fullElectFrequency  2
stepspercycle       10
#########################
zeroMomentum        yes
#########################


# Constant Temperature Control (Thermostat settings to keep fixed temperature)
langevin            off    ;# do langevin dynamics
langevinDamping     5     ;# damping coefficient (gamma) of 5/ps
langevinTemp        $temperature
langevinHydrogen    off    ;# don't couple langevin bath to hydrogens

#temperature coupling for temperature coupling
tCouple             on
tCoupleTemp         $temperature

if {1} {
# Periodic Boundary Conditions (specify the dimensions)
# center of the periodic box; NOT ORIGIN!!!
cellBasisVector1 101.5 0 0
cellBasisVector2 0 101.5 0
cellBasisVector3 0 0 71
cellOrigin 0 0 0
}
# The above numbers are approximated from the data we analysed through VMD.
wrapAll             on


# PME (Particle mesh ewald for effiecient calculation of full-system electrostatic calculations in periodic boundary conditions)
# multiples of 2,3,5 & >=dimensions above
PME                 yes
PMEGridSizeX        120
PMEGridSizeY        120
PMEGridSizeZ        90


# Constant Pressure Control (variable volume)
useGroupPressure      yes ;# needed for rigidBonds
useFlexibleCell       no
useConstantArea       no

langevinPiston       on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.
langevinPistonDecay   50.
langevinPistonTemp    $temperature


# Output
outputName          $outputname

restartfreq         500000 ;# every 1 ns(capturing restart points if at all the simualtion breaks, so that you can restart from this point)
dcdfreq             10000  #(capture MD trajectory data with frequency 10000 steps)
xstFreq             10000
outputEnergies      10000
outputPressure      10000


#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

# Minimization (to take the structure to a local minima using conjugate gradient method)
minimize            5000
reinitvels          $temperature ;#Reinitialize velocities(temperature is dependent on velocities of atoms)

run 10000000 ;# run the simulation for 20 ns
```
Now the task of NAMD software is to run this config file. For this purpose create a folder with necessary files in it. Like PDB, PSF, PRM and restart-files(necessary in case the simulation is restarted from a previous run input).

I know that there are a lot of missing points here. There are unclear things in the config file. But my soul aim in this article is to get started with NAMD.

The files `par_all36_lipid.prm` and `par_water_ions.prm` are available from [here](http://mackerell.umaryland.edu/charmm_ff.shtml) in compressed formats. But for convenience you can get it from [here](https://raw.githubusercontent.com/dexterdev/Molecular_Dynamics/master/par_water_ions.prm) and [here](https://raw.githubusercontent.com/dexterdev/Molecular_Dynamics/master/par_all36_lipid.prm).

### Running the config file
`devt@BUSYBOX ~/home/devt/SIM $ namd2 DMPC_patch_NPT_eq.config > log`

Below is a line from the log file generated:
`Info: Initial time: 1 CPUs 0.437884 s/step 2.53405 days/ns 211.227 MB memory`
From the log file, we can read that it takes 2.5 days for 1ns with 1CPU approximately! This is the case with 49805 atoms in the system.(We can get this number of total atoms either by querying using VMD software or opening the PDB file in a text editor and inspecting it.) We can improve this by using HPC(high-performance computer) machines.

## What this article will enable you to do?

This article simply allows you to see how a simulation run can be realized using NAMD. In the coming parts, we will explore some simulations of real systems like proteins etc by reproducing results from some research papers etc.

## Coming up in the next part...

- A detailed description of each filetype including generated output files
- More VMD and NAMD techniques/details
- How can we run the simulations faster?(if there is an access to HPC)
- Visualizing the simulation results and doing some basic analysis etc
- Periodic Boundary Conditions, Neutralizing total charge in the system etc
- A small system, which anyone can try out in their machine
- And more!



## References [Books for further reading]

###  Textbook references for learning theory of Molecular Dynamics:

- "Statistical Mechanics: Theory and Molecular Simulations" by Mark E. Tuckerman
- "Molecular Modelling: Principles and Applications" by Andrew R. Leach
- "Computer Simulation of Liquids" by D. J. Tildesley and M.P. Allen

###  References specific to NAMD and VMD:

- [NAMD userguide](http://www.ks.uiuc.edu/Research/namd/cvs/ug.pdf)
- [NAMD tutorial](http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix.pdf)
- [VMD userguide](http://www.ks.uiuc.edu/Research/vmd/vmd-1.9.1/ug.pdf)
- [VMD tutorial](http://www.ks.uiuc.edu/Training/TutorialsOverview/vmd/vmd-tutorial.pdf)

### CHARMM-gui related papers:

- http://www.charmm-gui.org/?doc=citations


---
### <center>Join #steemSTEM</center>
<center>Join the active science community #steemSTEM at discord: https://discord.gg/BZXkmWw</center>
<center>![](https://i.imgur.com/Ya06oXn.jpg)</center><center>Image courtesy: @elvisxx71</center>

#### <center>And to steemSTEM beginners:</center>
<center>You can ask for help in our discord page. There are people ready to help you there.</center>

<center>https://steemitimages.com/DQmPtnKZW6cQWBAzpA8yiqTDXSsVH9zEUyxFiek2kGEoA8g/bear.gif</center><center>gif courtesy: @rocking-dave</center>


---

**<center>All images without image sources are my creations :)</center>**

## <center>Follow me @dexterdev</center>

---

     ____ _______  ______ _________ ____ ______    
    /  _ /  __\  \//__ __/  __/  __/  _ /  __/ \ |\
    | | \|  \  \  /  / \ |  \ |  \/| | \|  \ | | //
    | |_/|  /_ /  \  | | |  /_|    | |_/|  /_| \// 
    \____\____/__/\\ \_/ \____\_/\_\____\____\__/