Dear Dr. Case,
After some tests with my input files it seems that leap writes out incorrect ATOMS_PER_MOLECULE sections when different (i.e. TER-separated) molecules are connected by additional bonds using the "bond"-command. The hemes in my system are covalently bound to the protein via histidine ligands and cysteine linkages, and the necessary bonds I have to add myself upon loading the pdb file. Leaving out these bond commands leads to a correct ATOMS_PER_MOLECULE section. However, the bonds connect my ten hemes with two of the three protein chains (and these with each other), and the first number in the incorrect ATOMS_PER_MOLECULE section (8265) is the combined number of atoms in these twelve separate species.
Following your suggestion, I've created a simple test system that reproduces this error - a sequence of four residues incl. two cysteines (CYX) in two separate chains. Simple loading the pdb file, setting the box and writing out a topology file using saveamberparmpol leads to a correct ATOMS_PER_MOLECULE section; however, if I manually bind the two cysteines before writing out the topology file the two subsets of atoms are added up again in the ATOMS_PER_MOLECULE section. By the way, without setting the box in leap using "set [...] box { ... }", no ATOMS_PER_MOLECULE section is written at all. Is it supposed to be like that? Is a different decomposition scheme used in cluster conditions that does not require the ATOMS_PER_MOLECULE section?
I've attached the test input pdb; I used the following commands in xleap to create a topology file from it:
tmp=loadpdb [Path]/CYX_only.pdb
set tmp box { 110 100 110 }
bond tmp.2.SG tmp.4.SG
saveamberparmpol tmp [Path]/[Topologyfilename] [crdfilename]
All the best,
Marian
________________________________________
From: David A Case [case.biomaps.rutgers.edu]
Sent: 04 November 2011 17:39
To: AMBER Mailing List
Subject: Re: [AMBER] Issue with running sander with SHAKE
On Fri, Nov 04, 2011, Breuer, Marian wrote:
>
> Thank you so much for this information. I used indeed an input pdb that
> I had written out from xleap before; but I checked it and it did contain
> all necessary TER specifiers. The hemes (incl. "HER") had already been
> at the beginning above the protein atoms in the original pdb, and the
> number of protein atoms written out by xleap (8365) is not even equal to
> the actual number of protein atoms (8363; maybe xleap somehow counted
> two TER lines in between the protein atom lines that terminate two
> subchains of my protein). I thus would give parmed a try.
It would help if you could post the input pdb file plus the exact commands you
used to make the prmtop file. Or maybe the simplest system that shows the
error (say with just one HEM group?)--it doesn't have to lead to a usable
system, just an example that leads to bad pointers in the prmtop.
[My best guess: there used to be some pieces of code in LEaP that ordered
molecules by size (greatest number of atoms first, then go in descending
order). I had thought that such behavior had all been removed, but
something might be left. Since "real" PDB files almost always have heme
groups and other ligands after the protein, instead of before it, there
might be a bug lurking related to this. If so, putting the HEM groups
after the protein might make a difference.]
But if you can provide any example that doesn't work correctly, that is a
great help in fixing the problem. Using parmed is only a work-around...we
should try to track down the problem at its source.
...thx...dac
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Nov 05 2011 - 16:00:02 PDT
