Introduction to molecular dynamics simulations

Note: if you encounter an english word that you don't know, you can

visit the LEO online dictionary to

translate it into german. Otherwise, just as with any other question,


If you encounter a problem accessing the internet, please configure your

browser to use the following proxyserver:

www-cache.cip.physik.local with port 3128

Often it is necessary to understand the dynamics of a biomolecule in

order to understand its function. For proteins, much information can

usually be derived from the structure, or even based solely on the

sequence, but the detailed functional mechanism often includes a

structural transition or atomic fluctuations of the protein or

substrate that can only be understood if the underlying atomistic

dynamics are known in detail. The functional dynamics

of proteins typically take place on a picosecond to millisecond

timescale. Unfortunately, there are hardly experimental techniques

available with a high enough time and spatial resolution to follow the

atomistic motions step by step, in order to be able to view proteins

"in action", like under a microscope. Therefore, computer simulations

of the atomic motions are the methods of choice to study biomolecular

dynamics. Such simulations allow insights into the dynamics involved

in e.g. enzyme catalysis, molecular motors, molecular transport, or

molecular recognition, in atomic detail. Today we will learn how to

set up such a so-called molecular dynamics simulation, first on a

simple system (argon gas), after which we will turn to a real protein.

MD simulation of Argon

As a first case study, we consider the dynamics of Argon gas. Argon

is a noble gas, characterised by a high degree of inertness (i.e. it

is chemically relatively unreactive). Physically, it is a relatively

simple compound, as it is mono-atomic both as gas and as fluid (no

covalent bonds that form molecules) and because it is apolar.

Therefore, it can be modeled as a neutral compound. Concerning the

"force-field", as it has been introduced during the lectures, it means

that in argon the only relevant force-field terms are the

Pauli-repulsion (that prevent atoms to overlap) and London-dispersion

(induced dipole-dipole interactions), which together are described by

a Lennard-Jones potential. Hence, in simulation jargon, Argon and

other noble gases are also refered to as Lennard-Jones fluids.

The simulations will be carried out with the GROMACS simulation package. On the

GROMACS homepage you can find both

the software and documentation (online reference and paper manual).

To run a simulation, three input files are usually required:

  • a structure file, containing the atomic coordinates of the system

    to be simulated

  • a "molecular topology", containing the atomic simulation

    parameters (force-field) and a description of the bonds, bond angles

    etc (if any) of the simulation system

  • the simulation parameters: type of simulation, number of steps,

    simulation temperature, etc.

Condensation and boiling of Argon

For the Argon simulations, three files have already been

prepared. Download argon_start.pdb, and md.mdp to your local working

directory without changing the filenames. Please feel free to

have a look at these files (on the command prompt, you can open these

files with any of the following programs: more, less, xedit, emacs, vi,


Gromacs uses the preprocessor grompp to collect the data from

the three input files, check the internal consistency, and write one

unified input file, topol.tpr On the command prompt in a linux

terminal or shell, type:

grompp -f md.mdp -c argon_start.pdb -p

Please note that the commands in bold print can be easily transferred to the

command prompt with copy-and-paste (select text by dragging the mouse over it

with the left mouse button pressed, and paste by pressing the middle mouse


Carefully check the output to see if any errors or warnings ocurred.

Note that, sice we didn't explicitly specify a name for the output

file, the default name of topol.tpr was used.

If everything was OK, then start the simulation:

mdrun -s topol.tpr -v -c argon_1ns.gro -nice 0

As mdrun says, 500,000 MD steps are carried out, each of 2 femtoseconds, making

a total of 1 nanosecond. Do a

ls -lrt

to see which files were written by mdrun. We see the following files:

  • traj.xtc - the trajectory to be used for analyses

  • traj.trr - the trajectory to be used for a restart in case of a crash

  • ener.edr - energies

  • md.log - a LOG file of mdrun

  • argon_1ns.gro - the final coordinates of the simulation

The first analysis step during or after a simulation is usually a

visual inspection of the trajectory. For this we will use the gromacs

trajectory viewer ngmx. Other possibilities would be VMD , Pymol, or gOpenmol.

ngmx -s topol.tpr -f traj.xtc

Select a group for visualisation ("system" might be a good choice),

and then, under "Display", click "Animate". Now a animation player

button appears near the bottom of the screen, which can be used to

view the trajectory. If the movie displays too fast, change the refresh rate

under "Display -> Options", and specify a Wait time of e.g. 3 ms.

We have simulated Argon at 100K, which is above

the boiling point, so we have simulated argon gas. Before we continue,

rename traj.xtc so that we can easily find it back later:

mv traj.xtc gas.xtc

In the next step,

we are slowly going to cool down the system (simulated annealing), to

see if we can simulate the condensation to liquid Argon. For this,

download anneal.mdp and start grompp with:

grompp -f anneal.mdp -c argon_start.pdb -p

Now start the annealing simulation:

mdrun -c argon_anneal.gro -nice 0

(note that we do not need to specify "-s topol.tpr" - as mdrun needs a

tpr file to run, it will automatically search for the default

filename: topol.tpr. Also note, that, although we now write to the

same trajectory files as before, the old files are not overwritten,

but are backed up automatically). Check the result with:

ngmx -s topol.tpr -f traj.xtc
Do you see the formation of

liquid argon? To see better what's going on, let's view the trajectory

with pymol. To do so, we first need to convert the trajectory to PDB


trjconv -s -f traj.xtc -o traj.pdb

and select "0". Then view the trajectory with pymol:

/usr/global/whatif/pymol/pymol traj.pdb

On the white command prompt in pymol, type:

show spheres

show cell

and press the "play" button on the bottom-right of the main pymol

window. Now analyse the energies and temperature from the simulation:

g_energy -o temperature.xvg

and select "7" for the temperature (and "0" to stop the selection).

Watch the result with:

xmgrace temperature.xvg

It can be seen that the actual temperature during the simulation

indeed monotonously decreases as the simulation progresses. Now

analyse the potential energy:

g_energy -o potential.xvg

and select "4" (and "0" to stop the selection).

Watch the result with:

xmgrace potential.xvg


Why is there a drop in the potential energy?


What is the boiling temperature of Argon?

Now repeat the annealing simulation, but extend the length of the

simulation by a factor of five. Do this by modifying the following

lines in "anneal.mdp" (open "anneal.mdp" in your favourite editor,

like xedit, nedit, kate or emacs):

nsteps = 2500000

annealing_time = 0 5000

and save "anneal.mdp". And run the simulation:

grompp -f anneal.mdp -c argon_start.pdb -p

mdrun -c argon_anneal.gro -nice 0

and analyse like above with g_energy and xmgrace.


Why is the final potential energy lower than in the previous simulation?


What is the boiling temperature this time? Why is the apparent boiling point higher than in the previous simulation?

Now, let us simulate the boiling instead of condensation of argon. For

this, download heatup.mdp and start the


grompp -f heatup.mdp -c argon_anneal.gro -p

Note that we use the final structure of the previous simulation as a

starting configuration.

mdrun -c argon_heatup.gro -nice 0

Analyse the results as above with g_energy and xmgrace.


What is now the apparent boiling temperature? Can you explain the

difference to the cooling simulations? Which of the three values would you

expect to be most realistic? Check your results by searching google

for the true boiling temperature of Argon.

Freezing and melting of Argon

In the annealing simulations, we have not only observed the

condensation to the liquid phase, but at lower temperatures also

the freezing towards the solid phase. To concentrate on this phase

transition more explicitly, we are now going to simulate the condensed

phase. A simulation box has already been prepared.

Download liquid.gro, and 94k.mdp to your local working directory. Update your

topology ( and change the number of atoms (last line) from

100 to 216. Start the simulation of liquid argon with:

grompp -f 94k.mdp -c liquid.gro -p

mdrun -c liquid.pdb -nice 0

To analyse the dynamic behaviour of liquid argon, we'll calculate the

diffusion constant from the simulation. The Einstein formula can be

used, which relates the mean atomic displacement to the diffusion

constant. For this, we'll use the program g_msd:

g_msd -s -f traj.xtc -trestart 50 -b 100

The option '-trestart 50' means that time windows of 50 ps are used

for the analysis (to enhance the statistics), and '-b 100' means that

the first 100 ps are excluded from analysis. Select '0' or '1' when

asked for a group.


What is the calculated diffusion constant? How does it compare to the

experimental value of 2.43x10-5 cm2/s?

One way to analyse the structural properties of a liquid is the

so-called radial distribution function, which shows the probability of

finding an atom center at a certain distance around an atom.

g_rdf -s -f traj.xtc -b 100 -o rdf_liquid.xvg

and view the radial distribution function with:

xmgrace rdf_liquid.xvg


What is the nearest-neighbor distance for two argon molecules? Why is

the second-neighbor peak not exactly at twice the nearest-neighbor distance?

Now let us slowly cool the system. Download anneal2.mdp and start the simulation:

grompp -f anneal2.mdp -c liquid.gro -p

mdrun -c anneal2.pdb -nice 0

Check the animation (with pymol) and the potential energy.


When does the freezing begin? When would you have expected the

freezing to start? (check the phase diagram above).

Compute the radial distribution function for the solid state,

by taking into account only the second half of the simulation:

g_rdf -s -f traj.xtc -b 500 -o rdf_solid.xvg

and view the result together with the RDF of the liquid state:

xmgrace rdf_solid.xvg rdf_liquid.xvg

(the solid state is in black, the liquid state in red).


Can you explain the observed differences?

Go to Editor View