Re: [AMBER-Developers] FW: [AMBER] Issue with running sander with SHAKE

From: Duke, Robert E Jr <>
Date: Wed, 7 Dec 2011 02:36:00 +0000

If memory serves, I did hit some problems with incorrect molecule definition in the prmtop in pmemd development, and I responded by basically creating a default mode in pmemd whereby pmemd checks over the bonding patterns and redefines molecules. This definitely created differences in the virial, but struck me as much more correct (I believe Darden concurred; we discussed it). I did also leave the old mode intact, but you have to request the old mode by setting a new namelist variable, no_intermolecular_bonds, to 0. There were also implications for frames definitions in extra points, so when using extra points ff's, setting no_intermolecular_bonds to 0 was not permitted. Anyway - just a head's up - I think pmemd by default already handles this issue.
Regards - Bob

From: case []
Sent: Tuesday, December 06, 2011 5:43 PM
To: Breuer, Marian
Subject: Re: [AMBER-Developers] FW: [AMBER] Issue with running sander with SHAKE

On Tue, Dec 06, 2011, Breuer, Marian wrote:
> Of course, this does
> unfortunately not provide any explanation why sander seemingly bothers
> about molecules (rather than residues) in the first place.

Sander handles constant pressure simulations via a "molecular virial" model,
in which internal forces inside a "molecule" can be ignored. For constant
pressure runs, it also distributes the atoms among processors, dividing things
along molecule boundaries. But this code (and similar code in LEaP) assumes
that the atoms in a molecule are contiguous in the atom list.

> Andrew Voronkov (who described a similar problem on the mailing list
> recently) told me that in his system (protein + ligand + zinc atom) the
> zinc atom was covalently bound to the protein but the cofactor was not,
> and as he described on the mailing list, the solution to his problem
> was to place the zinc atom between protein and cofactor rather than
> after the cofactor. I'd then assume that this was an analogous situation
> (avoiding bound atoms to be split in the atom list).

I agree that this is the same problem.

I'm cc-ing this to the Amber developers' list; it is rather odd that this
problem has not shown up earlier, since the code involved is quite old.

On the to-do list (increasing level of difficulty):

1. LEaP definitely should check for this, and certainly should not create
an ATOMS_PER_MOLECULE data structure that is incorrect. At an even more
basic minimum, Dave Cerutti's [or Jason's] prmtop validation/editing program
could also flag this problem.]

2. The Users' manual should document this requirement, and error messages
about its violation should point to such documentation.

3. The pmemd sanity check on ATOMS_PER_MOLECULE should be backported to

4. (Longer term) better barostats use an "atomic" model for volume
fluctuations, avoiding the rather arbitrary division into molecules.

> But an additional reason for my slow reply: my memory of the code is that
> the splitting of atoms should depend on residues, not molecules. So I want
> to see if I can understand how your original problem (with > 14 nodes) arose.

As per usual, my memory was wrong. For constant pressure runs, splitting
of atoms among processors depends on molecules, not residues [as comments
in the set.f file (that I think I wrote!) clearly indicate.] If the
information in the prmtop file is wrong, things can be badly messed up.

...thanks for your input....dac

AMBER-Developers mailing list
AMBER-Developers mailing list
Received on Tue Dec 06 2011 - 19:00:03 PST
Custom Search