Problem in writing NOPBC trajectory

[Vikas Dubey]

Hi everyone, 

I am trying to the extract trajectories of certain atoms and this trajectory should not have any periodic boundary condition effect.

People who use GROMACS, its similar to "gmx trjconv -f test.xtc -pbc nojump -n index.ndx -o nopbc.xtc". I have written the following code to achieve this but still my trajectory seem to have pbc enabled. What could be the possible problem with the code ?

mda.core.flags['use_periodic_selections'] = False

mda.core.flags['use_KDTree_routines'] = True

u=mda.Universe("test.pdb","semifinal.dcd")

phosphate_0=u.atoms.select_atoms("name P")

protein_0=u.atoms.select_atoms("protein and name CA")

phosphate_protein_0=u.atoms.select_atoms("protein and name CA or name P")

phosphate_protein_0.write("output.pdb")

with mda.Writer("output.dcd", phosphate_protein_0.n_atoms) as W:

    for ts in u.trajectory:

        W.write(phosphate_protein_0)

[Max Linke]

Hi

We don't have automatic removal of PBC boundaries yet. Also your code is
just actually doing right now is to write out the coordinates of your
selection unchanged. You have to write the PBC correction yourself. If
you don't want to do that (can be tricky for some geometries) you can
also use MDAnalysis to convert your trajectory into a XTC/TRR and then
use the gromacs tools to remove the PBC jumps (assuming you have set the
correct box for your simulation). Later analysis can then be done in
MDAnalysis.


best Max

> --
> You received this message because you are subscribed to the Google
> Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to mdnalysis-discus...@googlegroups.com
> <mailto:mdnalysis-discus...@googlegroups.com>.
> To post to this group, send email to
> mdnalysis-...@googlegroups.com
> <mailto:mdnalysis-...@googlegroups.com>.
> Visit this group at https://groups.google.com/group/mdnalysis-discussion.
> For more options, visit https://groups.google.com/d/optout.

[Richard Gowers]

The closest we have currently is the function MDAnalysis.lib.mdamath.make_whole https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/lib/mdamath.py#L270

This takes a single AtomGroup (which must represent a single molecule) and removes any jumps within it.  It needs bond information, so you'll need to have used a topology which gives that information.

I think for now Max is right, write it out into a trr/xtc and use gromacs.

Richard

> an email to mdnalysis-discussion+unsub...@googlegroups.com
> <mailto:mdnalysis-discussion+unsub...@googlegroups.com>.

[Vikas Dubey]

Sounds good ! Thanks a lot Max and Richard. I will use gromacs tools in between.