Re: [AMBER] mmpbsa number of atoms issues

From: Jason Swails <>
Date: Sun, 15 Sep 2013 00:05:31 -0400

On Sat, Sep 14, 2013 at 6:44 AM, Ayesha Fatima <>wrote:

> Dear All,
> I am trying to calculate the binding energy for my complex. I run into a
> familiar issue of incorrect atom numbers between prmtop and crd files. I
> have searched the tutorials and also the mailing list but i cannot find a
> specific solution. I am sure it is easy but i am unable to put my finger on
> it.
> The complex i simulated has water molecules in the binding site which i
> want to retain. when i want to generate the prmtop files, i am confused.
> the order of the atoms in the complex is protein-ligand-water in binding
> site. I prepared the prmtop files for the complex for mmpbsa by stripping
> the added water besides the water in the binding site. The prmtop file
> prepared for the prtotein alone also has the water in the binding site. The
> prmtop file of the ligand is of the ligand only, nothing extra.
> Q1 Where should i get the value for the RSTART AND RSTOP AND the LSTART and
> LSTOP? Which pdb file to consider. is it the one which is used in the
> simulation? This is by far the most confusing part since this the file to
> generate the snapshots because when i want to write the RSTART and RSTOP, i
> am not sure what should be written.

I'm not very familiar with the Perl script, but you can check to see if you
can specify multiple values of RSTART and RSTOP. If not, I don't see any
easy way of using to compute your binding free energy (it may be
worth studying the manual and examples in this case).

> Q2 How do i account for the water molecules in the binding site? Should
> they be removed also when preparing the prmtop files for the mmpbsa?

They are typically removed, but if the waters play a crucial, structural
role then implicit solvent models are unlikely to treat them properly
(except perhaps 3D-RISM). Continuum models treat all solvent-occupied
regions as bulk solvents, in general (the GB and PB approaches, that is).
 If a water molecule behaves very differently in your protein than it does
in bulk, your answers may be very dependent on that water molecule being

You can always test your ideas by computing binding free energies of
different systems (with and without explicit waters).

> Q3 Could there be an issue with the order of the complex atoms? Should it
> be protein-ligand-water in binding site or protein-water in binding site-
> ligand?

It's best not to think of your system as a 3-part system. You have 2 parts
-- a receptor and a ligand. Whether you consider the water as part of the
receptor or ligand depends on what you want to calculate and what you are
trying to learn. In either case, tleap will likely create the topology
file with the water at the end (but otherwise the residue sequence will be
the same as it is in the original PDB file).

> Q4 if i change the order of atoms, should i run the md again? by right i
> should, but is there a way around the problem?

Have you tried I know that it can handle such systems without
having to jump through many hoops. (And changing the order of atoms in a
topology file is no easy task -- especially when tleap tends to put waters
at the very end).

Finally, you can also compute MM/PBSA free energies 'by hand' -- that is,
don't use either or and run the necessary energy
calculations manually. If you have never done this before, it may be a
useful exercise -- it should yield insight into the theoretical basis of
the method as well as what the scripts are automating at each step (which
should in turn help you learn how to use the scripts more easily and


Jason M. Swails
Rutgers University
Postdoctoral Researcher
AMBER mailing list
Received on Sat Sep 14 2013 - 21:30:02 PDT
Custom Search