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.
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). GROMACS is
already installed on the machines of the CIP pool, so no need to download or
install any software.
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.
Go back to Contents
Condensation and boiling of Argon
For the Argon simulations, three files have already been
prepared. Download argon_start.pdb, argon.top 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).
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/argon_start.pdb
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/argon.top
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/md.mdp
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:
gmx grompp -f md.mdp -c argon_start.pdb -p argon.top
Please note that the commands in the gray boxes 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:
gmx 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
to see which files were written by mdrun. We see the following files:
- traj_comp.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 VMD
as a trajectory viewer. Other possibilities would be Pymol, or gOpenmol.
vmd argon_start.pdb traj_comp.xtc
Select Graphics -> Representations and then as "Drawing Method": VDW.
As standard the view on the system is shown as perspective. To change
this click Display -> Orthographic.
To display the box click Extensions -> Tk Console and type in
You can display the trajectory as a movie by pressing the arrow symbol at the bottom.
To change the speed of the movie you find a controller next to the arrow symbol.
We have simulated Argon at 100K, which is above
the boiling point, so we have simulated argon gas. Close vmd when you are done.
Before we continue, rename traj_comp.xtc so that we can easily find it back later:
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:
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/anneal.mdp
gmx grompp -f anneal.mdp -c argon_start.pdb -p argon.top -maxwarn 1
Now start the annealing simulation:
gmx 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_comp.xtc
Now analyse the energies and temperature from the simulation:
gmx energy -o temperature.xvg
and select "8" for the temperature (and "0" to stop the selection).
Watch the result with:
It can be seen that the actual temperature during the simulation
indeed monotonously decreases as the simulation progresses. Now
analyse the potential energy:
gmx energy -o potential.xvg
and select "4" (and "0" to stop the selection).
Watch the result with:
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 gedit, or any of your favourite editors,
like xedit, nedit, kate or emacs):
nsteps = 2500000
annealing_time = 0 5000
and save "anneal.mdp". And run the simulation:
gmx grompp -f anneal.mdp -c argon_start.pdb -p argon.top -maxwarn 1
gmx mdrun -v -c argon_anneal.gro -nice 0
and analyse like above with gmx 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:
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/heatup.mdp
gmx grompp -f heatup.mdp -c argon_anneal.gro -p argon.top -maxwarn 1
Note that we use the final structure of the previous simulation as a
starting configuration.
gmx mdrun -v -c argon_heatup.gro -nice 0
Analyse the results as above with gmx 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
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.
Edit your topology file (argon.top) and change the number of atoms (last line) from
100 to 216.
Start the simulation of liquid argon with:
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/liquid.gro
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/94k.mdp
gmx grompp -f 94k.mdp -c liquid.gro -p argon.top
gmx 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 msd:
gmx msd -s -o msd.xvg -f traj_comp.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.
The diffusion constant can be read in the msd.xvg file (look for the line with D[System] = ...)
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.
gmx rdf -s -f traj_comp.xtc -b 100 -o rdf_liquid.xvg
Select "0" , "0" followed by Ctrl-D to start the calculation.
Now view the radial distribution function with:
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:
wget http://www3.mpibpc.mpg.de/groups/de_groot/compbio1/p1/anneal2.mdp
gmx grompp -f anneal2.mdp -c liquid.gro -p argon.top -maxwarn 1
gmx mdrun -v -c anneal2.pdb -nice 0
Check the animation (with vmd) and the potential energy.
vmd liquid.gro traj_comp.xtc
exit
gmx energy -o potential.xvg
xmgrace potential.xvg
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:
gmx rdf -s -f traj_comp.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
Optional
Using a simulation strategy similar to above, try to answer the following
questions:
-
Pressure dependence
From te phase diagram, it is clear that in addition to the temperature,
the phase behaviour of argon is also strongly dependent on the pressure. Is
this properly accounted for in the simulations? Try to see how the boiling
point and the melting point at 1000 bar differ from the phase transitions at an
ambient pressure of 1 bar. Is the observation consistent with experiment?
(Hint: look for the word "pressure" in the gromacs input .mdp file. If
necessary, also check here).
-
Sublimation
At low pressures, there is no liquid state of argon, but instead there is a direct
transition from the solid state to the gas phase, a process called
sublimation. Is this also observed in the simulation?
-
Simulation of Neon
How are the results different for the much smaller neon? (Hint: use this
set of simulation parameters - neon.top).
-
Collision of two droplets
Investigate if you can use the kinetic energy of a moving droplet to
induce a phase transition. In order to investigate this, shoot two droplets
towards each other with different opposite velocities. With which velocities do
you get a larged, merged droplet and from which velocities do you get a
transition into the gas phase? Hint: use this structure
file containing two droplets to start from. You'll also need
this index file, this topology
file, and this parameter file
Go back to Contents
Further references
- M. Karplus and A. McCammon. Molecular Dynamics simulations of
biomolecules Nature structural biology 9: 646-652 (2002).[link]
- D.C. Rapaport The Art of Molecular Dynamics Simulations - 2nd edn
Cambridge University Press (2004).
- H. Scheraga, M. Khalili and A. Liwo Protein-Folding Dynamics:
Overview of Molecular Simulation Techniques Annual Review of Physical Chemistry 58: 57-83 (2007).[link]
For questions or feedback please contact Bert de Groot / bgroot@gwdg.de
