On Thu, 2015-01-15 at 09:18 +0100, Francesco Pietra wrote:
> I am trying to verify chemical bonds commanded from leap with
>
> ---VMD
> ---parmed.py
> ---parm7 file
>
> I assume that in parm7 file:
> BONDS INC HYDROGEN
> This section contains a list of every bond in the system in which at least
> one atom is Hydrogen. Each bond is identified by 3 integers—the two
> atoms involved in the bond and the index into the BOND FORCE CONSTANT
> and BOND EQUIL VALUE. For run-time efficiency, the atom indexes are actu-
> ally indexes into a coordinate array, so the actual atom index A is
> calculated
> from the coordinate array index N by A = N/3 + 1. (N is the value in the
> topology file)
> %FORMAT(10I8)
> There are 3 × NBONH integers in this section.
> BONDS WITHOUT HYDROGEN
> This section contains a list of every bond in the system in which neither
> atom
> is Hydrogen. It has the same structure as BONDS INC HYDROGEN described
> above.
> %FORMAT(10I8)
> There are 3 × NBONA integers in this section.
>
>
> Following is the initial portion of my parm7 file (ambertools14):
> FLAG
> BONDS_WITHOUT_HYDROGEN
> %FORMAT(10I8)
>
> 30 33 1 30 36 2
Let's start with these. 30 33 1 means that atoms 30/3+1 and 33/3+1 are
defined by the first bond. These are atoms 11 and 12 -- the carbonyl
bond. These are parametrized by the first entry in BOND_EQUIL_VALUE and
BOND_FORCE_CONSTANT.
The next bond is 30 36 2, meaning that atoms 30/3+1 and 36/3+1 are
bonded to each other. These are atoms 11 and 13 -- the amide bond
between the first two residues. They are parametrized by the second
entry in BOND_EQUIL_VALUE and BOND_FORCE_CONSTANT.
> Could you please instruct me how to verify bonds in the above list for the
> two initial residues in the corresponding PDB file?
Note that tleap typically re-orders atoms in the original PDB file so
that you can't use the original PDB file to compare atom numbers like
this. If you save a PDB from tleap or use ambpdb (or cpptraj) to create
a PDB file (which it appears you have done), then that PDB file works
for this.
> ATOM 1 N ALA 1 39.886 -67.858 8.967 1.00 0.00
> ATOM 2 H1 ALA 1 40.139 -66.882 8.910 1.00 0.00
> ATOM 3 H2 ALA 1 39.435 -68.128 8.104 1.00 0.00
> ATOM 4 H3 ALA 1 40.681 -68.381 9.306 1.00 0.00
> ATOM 5 CA ALA 1 38.841 -68.004 10.020 1.00 0.00
> ATOM 6 HA ALA 1 37.934 -68.417 9.578 1.00 0.00
> ATOM 7 CB ALA 1 39.276 -69.022 11.092 1.00 0.00
> ATOM 8 HB1 ALA 1 39.319 -68.531 12.064 1.00 0.00
> ATOM 9 HB2 ALA 1 38.557 -69.840 11.130 1.00 0.00
> ATOM 10 HB3 ALA 1 40.261 -69.416 10.841 1.00 0.00
> ATOM 11 C ALA 1 38.486 -66.675 10.668 1.00 0.00
> ATOM 12 O ALA 1 39.265 -65.721 10.620 1.00 0.00
> ATOM 13 N MET 2 37.288 -66.634 11.251 1.00 0.00
> ATOM 14 H MET 2 36.670 -67.430 11.184 1.00 0.00
> ATOM 15 CA MET 2 36.788 -65.488 12.003 1.00 0.00
> ATOM 16 HA MET 2 37.498 -64.664 11.926 1.00 0.00
> ATOM 17 CB MET 2 35.425 -65.022 11.462 1.00 0.00
> ATOM 18 HB2 MET 2 34.698 -65.828 11.565 1.00 0.00
> ATOM 19 HB3 MET 2 35.087 -64.154 12.027 1.00 0.00
> ATOM 20 CG MET 2 35.418 -64.613 9.987 1.00 0.00
> ATOM 21 HG2 MET 2 35.721 -65.461 9.373 1.00 0.00
> ATOM 22 HG3 MET 2 34.414 -64.298 9.701 1.00 0.00
> ATOM 23 SD MET 2 36.551 -63.245 9.652 1.00 0.00
> ATOM 24 CE MET 2 35.781 -61.890 10.516 1.00 0.00
> ATOM 25 HE1 MET 2 36.464 -61.506 11.274 1.00 0.00
> ATOM 26 HE2 MET 2 35.541 -61.096 9.808 1.00 0.00
> ATOM 27 HE3 MET 2 34.866 -62.237 10.995 1.00 0.00
> ATOM 28 C MET 2 36.619 -65.876 13.464 1.00 0.00
> ATOM 29 O MET 2 36.237 -67.031 13.793 1.00 0.00
>
>
> Embarrassed enough to pose such a naive question, however I am in trouble
> in verifying the active center with transition metals (where not all bonds
> indicated by both leap log and parmed.py are displayed on VMD; I know that
> in such a case VMD only relies on bonds registered in the parm7 file,
> abandoning any guess on interatomic distances)
Going through each bond entry like this will be much too tedious. That
was part of the motivation behind writing ParmEd (and, originally,
rdparm). ParmEd does exactly what I described above when listing the
bonds. So if ParmEd reports that a bond exists, then the prmtop defines
that bond.
VMD *should* report all bonds when loading a prmtop (but if you load a
PDB file -- even generated from a prmtop file -- then VMD will do a
simple bond-by-distance criteria when displaying them). There are some
representations (e.g., "dynamic bonds") that ignores this "bond"
information and simply uses a distance-based cutoff. But most
representations should respect the bond info in the topology file.
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 Jan 15 2015 - 05:30:03 PST
