Practical 1: Introduction to molecular dynamics simulations


A. Introduction

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.

Go back to Contents

B. 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:

Go back to Contents

C. 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, kate).

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

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: The first analysis step during or after a simulation is usually a visual inspection of the trajectory. For this we will use VMD (Visual Molecular Dynamics). Other possibilities would be Pymol, or gOpenmol.

vmd argon_start.pdb traj.xtc
In the main VMD window, select the menu Graphics -> Representations. In the new "Graphical Representations" window, choose "VDW" in the drop-down menu "Drawing Method". Back in the main VMD window, click the Play button, and play with the "speed" bar. 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 -maxwarn 2
Now start the annealing simulation:
mdrun -v -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:
vmd argon_start.pdb traj.xtc
Do you see the formation of liquid argon? In the following, let us also visualize the simulation with PyMol. To do so, we first need to convert the trajectory to PDB format:
trjconv -s -f traj.xtc -o traj.pdb
and select "0". Then view the trajectory with 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
Question: Why is there a drop in the potential energy?

Question: 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 -maxwarn 2
mdrun -v -c argon_anneal.gro -nice 0
and analyse like above with g_energy and xmgrace.

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

Question: 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 simulation:

grompp -f heatup.mdp -c argon_anneal.gro -p -maxwarn 2
Note that we use the final structure of the previous simulation as a starting configuration.
mdrun -v -c argon_heatup.gro -nice 0
Analyse the results as above with g_energy and xmgrace.

Question: 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.

Go back to Contents

D. 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 -maxwarn 2
mdrun -v -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.

Question: 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

(Note that if the last command didn't work (didn't produce output for more than 30 seconds), please download a working version of g_rdf from here, and make it executable with:

chmod +x g_rdf
and repeat the g_rdf command from above).

Now view the radial distribution function with:

xmgrace rdf_liquid.xvg
Question: 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 -maxwarn 2
mdrun -v -c anneal2.pdb -nice 0
Check the animation (with pymol) and the potential energy.

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

Question: Can you explain the observed differences?

Go back to Contents

E. Optional

Using a simulation strategy similar to above, try to answer the following questions:
For questions or feedback please contact Bert de Groot /