On Thu, 2015-07-02 at 14:44 +0100, Hannes Loeffler wrote:
> Hi,
>
> I have found an inconsistency with leap when creating the prmtop for the
> same molecule (neopentane in this case) but in different atom order.
> The difference in the two prmtops is in the list of excluded atoms. See
> attachment.
There is nothing inherently wrong with this. I grant that it's a little
surprising, but everything is indeed consistent. If you compute
energies and forces for these two systems, you will find that they are
identical.
The reason you seemingly have "extra" exclusions is due to the way that
the Amber prmtop file stores exclusions. Each exclusion is stored once,
and only once in the exclusion list. And that exclusion is always added
to the exclusion list belonging to the atom with a *lower* index. When
you use this convention, there are many atoms that have *no* assigned
exclusions since all of their bonded partners have lower indexes.
The convention that Amber uses for atoms that have no exclusions is to
put a "1" in the NUMBER_EXCLUDED_ATOMS section for that atom and put a
"0" for that atom's "segment" of the EXCLUDED_ATOMS_LIST section (which
points to no atom by virtue of Amber beginning its indexing from 1 in
the standard Fortran convention). To show how atom ordering can affect
the apparent number of exclusions (and the NNB pointer), consider
methane, and the following two orders for the atoms: C, H1, H2, H3, H4
and H1, C, H2, H3, H4. (Note in this example I am only considering 1-2
exclusions for simplicity and clarity. Hopefully it is clear that this
demo readily extends to the real situation where 1-3 and 1-4 pairs are
also considered.)
In the first case, C is bonded to each H and "owns" those exclusions
(since it has the lowest index). So it will seemingly have 4
exclusions. The remaining hydrogens have no exclusions, so they will
get a placeholder "1" and be excluded from atom 0. So
NUMBER_EXCLUDED_ATOMS will look like [4, 1, 1, 1, 1] and
EXCLUDED_ATOMS_LIST will look like [2, 3, 4, 5, 0, 0, 0, 0]. And NNB
will be 8.
Now consider the second case, H1 is bonded to C and owns that exclusion
(since it comes first). C owns the remaining exclusions. And the last
3 hydrogens get 0 place-holders. So the NUMBER_EXCLUDED_ATOMS will look
like [1, 3, 1, 1, 1] and EXCLUDED_ATOMS_LIST will look like [2, 3, 4, 5,
0, 0, 0]. And NNB will be 7. Basically, what happened was a "dummy"
exclusion was effectively replaced by a real one, which is where the
"extra" exclusions came from in the first ordering.
In your examples, NNB differed by two between the two systems, so this
must have happened twice. Rest assured that when I used ParmEd to make
sure every atom's position was exactly the same, then computed the
energies, they were identical to well past the full single precision
that was used to compute the energies and forces. The relevant section
from my IPython session is shown below as a demonstration:
In [42]: for atom in y.atoms:
...: for atom2 in x.atoms:
...: if atom.name == atom2.name:
...: atom.xx, atom.xy, atom.xz = atom2.xx, atom2.xy, atom2.xz
...:
In [43]: pmd.tools.energy(y, 'omm', 'decompose').execute()
Bond = 59.5205786 Angle = 0.5010947
Dihedral = 0.0039229 Nonbond = 10.2030082
TOTAL = 70.2286044
In [44]: pmd.tools.energy(x, 'omm', 'decompose').execute()
Bond = 59.5205786 Angle = 0.5010947
Dihedral = 0.0039229 Nonbond = 10.2030080
TOTAL = 70.2286042
This is one of the very many quirks of the topology file that I learned
in the process of writing ParmEd. You can see where I duplicated this
logic here:
https://github.com/ParmEd/ParmEd/blob/master/parmed/amber/_amberparm.py#L1546-L1554
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jul 02 2015 - 08:00:02 PDT
