Re: [AMBER] Curious behavior of cpptraj in ambertools 13.1

From: Jason Swails <jason.swails.gmail.com>
Date: Wed, 19 Jun 2013 13:39:55 -0400

On Wed, Jun 19, 2013 at 1:10 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:

> Hi,
>
> On Wed, Jun 19, 2013 at 10:49 AM, Fabrício Bracht <fabracht1.gmail.com>
> wrote:
> > Hello. I am experiencing some issues with ambertools 13 while trying to
> > analyse different systems. The first system is a "normal" pmemd.cuda
> > simulation of a protein in water. When I type the command cpptraj -i
> > rmsd.ptraj -p HOH_wat.prmtop, I get the following error:
>
> CPPTRAJ currently processes flags on the command line in order. This
> means that what is happening is your input file (rmsd.ptraj) is being
> read in before your topology file (HOH_wat.prmtop). However, since
> some of your input depends on the topology file (namely the 'trajin'
> command), it fails.
>
> > Error: Atom 5311 was assigned a lower molecule # than previous atom.
> > Error: This can happen if bond information is incorrect or missing.
> > Error: Detected # of molecules is 7615. If this is incorrect and your
> > Error: topology does not have bond information (e.g. PDB file), try
> > Error: increasing the bond search cutoff offset (currently 0.200000).
> > Error: e.g. 'parm HOH_wat.prmtop bondsearch <new offset>'
> > Error: This can also happen if the atoms in your molecules are not
> > Error: sequential (e.g. molecule 1 is atoms 1-4,10-14 and molecule 2
> > Error: is atoms 5-9,15-20). This can be fixed using the 'setMolecules'
> > Error: command in parmed.py.
> > Error: Could not open parm HOH_wat.prmtop
>
> As stated by the error message, this can be caused by either incorrect
> bond information or if atoms are not sequential within your molecules
> (I suspect the latter). Cpptraj requires that atoms in molecules be
> sequential for speed (among other things). Did you try using the
> 'setMolecules' command in parmed.py?
>

As a followup to this comment, this indicates an issue in which the
topology file is broken. While the setMolecules command in ParmEd will fix
the topology file, it does so by re-ordering the atoms inside so that each
molecule is represented continuously in the prmtop.

This invariably 'breaks' compatibility with the generated trajectory file
since the atoms are now in a different order. This issue is currently
caused in complex systems if you use the 'bond' command in tleap to link
together different 'molecules'. The idea is to use setMolecules alongside
a restart file right after the problematic topology is printed out with
tleap so that it reorders both the restart file and the prmtop.

Since this error is very frequently discovered in hindsight, I wrote a
little Python script awhile ago that uses some of the new ParmEd
functionality to reorder the prmtop and the trajectory at the same time.
 The link to the script is on my wiki:
http://jswails.wikidot.com/helpful-scripts#toc21

A couple caveats here:

1) This re-orders the prmtop and the trajectory simultaneously (use the
--help argument to the script to find out how to use it)

2) This script requires an up-to-date AmberTools 13 installation, and
AMBERHOME must be set. Therefore, the input trajectories and input prmtop
must have the same atom ordering.

3) This script will _only_ read NetCDF files (ASCII trajectories would
process WAY too slowly), and will only write out a NetCDF file. You can
provide as many input trajectories as you want, but it will only print out
one output trajectory.

4) This script requires Numpy as well as some version of the NetCDF
bindings, so either install ScientificPython, py-netcdf, or the netCDF4
python bindings.

If you need to fix a trajectory, you can't use ParmEd to fix your prmtop.
 In the future, it may not be a bad idea to use the 'checkValidity' command
in ParmEd for prmtops you create just to check that everything is OK (it
will also catch common mistakes, like not setting disulfide bonds, that may
save you lots of time and wasted simulation).

HTH,
Jason

-- 
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jun 19 2013 - 11:00:04 PDT
Custom Search