Practical: Combining Molecular Dynamics with Electronic Structure Calculations

This practical was written by Udo W. Schmitt.


In the previous lecture we've learned about Density Functional Theory-based electronic structure methods as a viable alternative to the Hartree-Fock approach.

Today we are going to combine the already introduced molecular dynamics (MD) approach with two electronic structure methods, which will lead us to the field of so-called ab initio molecular dynamics or "on the fly" (direct) MD. This is opposed to the common force field techniques that are used to represent the potential energy surface, where the potential energy (and thus the forces on the nuclei) are determined before we start our simulation. In ab initio dynamics the forces are computed "on the fly", i.e., at every new set of coordinates generated during our molecular dynamics run. To this end, for every new set of coordinates an independent electronic structure calculation needs to be performed in order to compute the interaction energy and the forces on the atoms (nuclei).

The package we going to use for the electronic structure computation is again MOLPRO (see previous practical), which has been developed by Peter Knowles and Hans-Jürgen Werner over the last 15 years. MOLPRO has been used in the previous practical as well.

In order to be able to perform "on the fly" dynamics, we also need a MD program that acts as a driver for the electronic structure calculation. Here we going to use speedEVB, which is a molecular dynamics packages that provides a numerically efficient implementation of the Multistate Empirical Valence Bond (MSEVB) model together with a combined direct dynamics/QM/MM interface. The term QM/MM stands for Quantum Mechanics/Molecular Mechanics and was together with the MSEVB model introduced in the previous lecture .

The files needed for this practical can be downloaded here. You will need to unpack the archive by typing:

tar xzvf prakt2.tar.gz
ulimit -s unlimited

Direct Molecular Dynamics using the Multistate Empirical Valence Bond (MSEVB) model

In this ection we are going to perform MD calculations on proton transfer through a so-called water wire. To be able to do so, we need a model of our potential energy surface that allows for explicitly break and form chemical bond. A model that allows for that in an elegant way is the Multistate Empirical Valence Bond (MSEVB) to describe the reactive potential energy surface (PES). The MSEVB model (also see previous Lecture) provides a simple and numerical efficient way of modeling chemical reaction dynamics in aqueous systems. The MSEVB approach is nowadays used extensively in biomolecular simulations to study proton transport phenomena in molecular detail.


A water wire is a onedimensional alignment of water molecules that are held together via hydrogen bonds. In case the water molecules are aligned properly, a proton can migrate through the water wire by the so-called Grotthus mechanism. Grotthus mechanism means that though a proton is transfered from one end to the other of the wire, the movment of the individual atoms is quite small, i.e., a few tenth of an Angstroem.

Aligned water wires are found in a variety of proteins, like e.g. in Bacteriorhodopsin or Cytochrome c Oxidase, where they constitute a crucial part of the proteins functionality. That means that lack of the water wire will lead to misfunction of the respective protein. In this course we are going to focus on the water wire plus one excess proton only and are mimicking the protein via a cylindrical trapping potential that will keep the water molecules aligned in single file.

First have a look at the starting structure located in the directory sim/ by typing

vmd start.pdb 

and locate the water molecule that carries the excess proton. Now we are going to perform a MSEVB molecular dynamics simulation with this starting structure

speedEVB wire.in

The program will report all its events to the screen output and will finish after 1000 MD steps. You can have a look at the generated trajectory by starting the VMD script:

vmd -e wire.vmd

Locate the excess proton and following its dynamical path through the water wire. How many bonds are formed and broken along the pathway? Besides the trajectory additional information has also been written to the file system. Have a look at the energies, i.e., kinetic energy, potential energy and total energy, by typing

gnuplot kin_energy.gplot
gnuplot pot_energy.gplot
gnuplot energy.gplot

What do you expect for the total energy of a system that follows Newton's (or Hamilton's) equation of motion?

In the directory transfer/ you can start a shorter simulation run by typing

speedEVB wire.in

Now you can have a look at the amplitudes of the individual ground state electronic expansion coefficients for the localized valence bond states as a function of simulation time. The amplitudes give you the fraction of the total charge of +1 (the proton charge) that resides on each water molecule at a given time. Typing will display the results:

gnuplot evb_amplitudes.gplot

How does the time evolution of the electronic expansion coefficients reflect the proton transfer events seen in the previous VMD movie?

VMD also allows you to visualizes the electronic expansion coefficients as changing color/size on the individual water molecules. Simply type

vmd -e amplitudes.vmd

Identify on how many water molecules the excess charge is distributed at a given time. Can the transfer process be described in a picture where the excess charge of +1 is localized on a single water molecule?

Direct Molecular Dynamics using Density Functional Theory in Combination with a QM/MM ansatz

With the same system as in the previous example - a protonated water wire - we are now going to perform mixed quantum-classical (QM/MM) molecular dynamics simulation. The first thing we have to define is the extend of the QM region, i.e., the part of the total system that is treated by means of an electronic structure method. The rest of our system will be handled by standard force field. The QM region will be the core region of our protonated water cluster that is embbeded in the overall water wire and is simply defined in the input files (fort.10) as:

QM_REGION
  10
  1 2 3 4 5 6 7 8 9 10



That just means that the first ten atoms will be treated as quantum-mechanically (QM_REGION), and the rest will be considered as being part of the MM treatment. See the picture to the right that shows the QM (H7O3+) cluster in thick ball-and-stick representation.

For the QM computation, which will be performed with MOLPRO, we will have to specify the electronic structure method that will be used (Hartree-Fock, Density Functional Theory, Moeller-Plesset etc.) and an appropriate basis set. This is done in the input file "qm.com":

basis=midi

geomtyp=xyz
GEOMETRY=qm.xyz

set,charge=+1

ks,b3lyp

Compare to last weeks practical only the command "ks, b3lyp" (which stands for Kohn-Sham) appears as the new DFT key word, and "b3lyp" specifies the density functional that is used during the calculation.

To start our simulation, which will perform 40 molecular dynamics steps, just type

speedEVB_qmmm wire.in

You can follow the screen output: this calculation will take a couple of minutes.

After successful finish you can look at the generated trajectory by typing

vmd -e wire_qm.vmd

QM/MM Molecular Dynamics of Protonated Water Clusters in Bacteriorhodopsin

Our last task for this semester will be to compute a short trajectory using the technqiues introduced above on a model of the light-driven proton pump Bacteriorhodopsin. Bacteriorhodopsin absorbes light in the visible range and converts this energy via a cascade of proton transfer steps into chemical energy. It can be regarded as a simple model system for larger protein complexes that perform light-driven energy conversion, like for example the Photosynthetic Reaction Center that is found in plants.

All necessary input files are given in the working directory "bR_QMMM", so you only need to type

speedEVB_qmmm bR.in

to start the simulation. In total 20 steps will be performed and can be displayed using VMD

vmd -e bR.vmd

Shown are the protein in a so-called "Ribbon" representation that highlights secondary structure elements and the QM region consisting of 3 water molecules and one excess proton at the extracellular entrance area of Bacteriorhodopsin.

tail -100 fort.10

will show you the simple definition of QM atoms. Note that we could select any region within the molecule to be treated at the high-level of theory - Hartree-Fock or Density Functional Theory - in order to compute observables with high accuracy.