1

I would like to extract one chain from my molecular dynamics trajectory (xtc file) using MDAnalysis. I expected it to be very simple, but an error occurred and I am not sure why I am getting it. Here is the code:

import MDAnalysis as mda

u = mda.Universe('trajectory1.xtc')
protein = u.select_atoms('segid A')
protein.write('trajectory1-A.xtc')

And it returns error: AttributeError: AtomGroup has no attribute segids

On MDAnalysis page, "segid" is used with ".select_atoms". Is the problem that I have trajectory file and not classic PDB file?

eyllanesc
  • 235,170
  • 19
  • 170
  • 241

2 Answers2

1

It is supposed that in the Universe creation (I feel like God) you have to provide the topology and the trajectory. As far as I remember, the xtc file only contains the trajectory (triplets of numbers, the coordinates for the atoms), but not the topology. With no topology, there are neither residues nor segments.

Check if you have some segment with:

u.atoms.segments
list(u.atoms.segments)

Probably it will return an empty set.

In that case, you should provide to the Universe creation a suitable topology with the segments you are interested in.

Something like:

import MDAnalysis as mda

u = mda.Universe('mytopology.pdb', 'trajectory1.xtc')
protein = u.select_atoms('segid A')

Now, you can traverse/process your trajectory in any way you see fit.

Poshi
  • 5,332
  • 3
  • 15
  • 32
  • Oh yes! Thanks! This was the problem, topology file. Another problem now is that it extracted one chain from trajectory, but only as 1 frame and not full trajectory of that chain. – HungryMolecule Feb 24 '20 at 09:43
  • Read the documentation thoroughly. A selection returns an `AtomGroup`, which is a group of atoms, in a single snapshot. If you want the next frame in the trajectory, you should advance to the next one: `universe.trajectory.next()`. – Poshi Feb 24 '20 at 10:13
  • Can you amend your answer with what @HungryMolecule has in their solution, i.e., show how to load a `Universe` with topology and trajectory? In particular, chains will only be available in a PDB file. IIRC, Gromacs TPR files or GRO files do not contain this information. – orbeckst Feb 25 '20 at 17:30
1

So, I revised the code because I found out that using trajectory file it is needed to have slightly modified code, because otherwise it will extract only one frame of the trajectory, but with this it will extract the whole trajectory of the one chain:

import MDAnalysis as mda

u = mda.Universe('mytopology.pdb', 'trajectory1.xtc')
protein = u.select_atoms('segid A')
with mda.Writer("traj1-Achain.xtc", protein.n_atoms) as W:
    for ts in u.trajectory:
        W.write(protein)

Thanks for help!

orbeckst
  • 3,189
  • 1
  • 17
  • 14