Combining Molecular Dynamics with Electronic Structure Calculations methods

Dr. Udo W. Schmitt, Max-Planck Institute for Biophysical Chemistry, Goettingen

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-Juergen 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


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


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


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:



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





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


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


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


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.

Go to Editor View