Create a sub-universe out of a gromacs trajectory

[mverissi]

Dear MDAnalysers,

I have four large GROMACS trajectories, each composed of 30,000 frames, with three chains of a protein solvated in water. There are roughly 13000 protein atoms, and the rest is water. This was necessary because it's a globular protein interacting with a coiled-coil. I'm interested in working only with a small subset of the protein atoms, such as

"chainID A and (resid 63-71 or resid 205-208 or resid 250-256 or resid 370 or resid 417-435)"
and
"(chainID B and resid 476-504) or (chainID C and resid 601-624)"

and analyse the hydrogen bonds between only these atoms. Those are at the interface between chain A, chain B and chain C. In principle, it would be possible to select the atoms directly with Gromacs, but I'm not sure how to create a new working TPR file that MDAnalysis could read.

To save time, I would like to consider the universe as a whole, including its entire trajectory, and create a new “small universe” comprising only these atoms, in which I would then analyse hydrogen bonds. Would this be possible? I would greatly appreciate any help.

Best regards,

Marcos Verissimo Alves
Universidade Federal Fluminense, Brazil.

[orbeckst]

You can write a subset of your system to new trajectory:

u = mda.Universe(TPR, "big.xtc")
subsystem = u.select_atoms("(chainID A and (resid 63-71 or resid 205-208 or resid 250-256 or resid 370 or resid 417-435)) or ((chainID B and resid 476-504) or (chainID C and resid 601-624))")

subsystem.atoms.write("subsystem.xtc", frames="all")

# also write a PDB as a topology file
subsystem.atoms.write("subsystem.pdb")

You can then load your system as Universe("subsystem.pdb", "subsystem.xtc").

MDAnalysis cannot write TPR files.

If you really need the bond information then you can do the following trick: load the whole system with bond information and then create a smaller universe to which you attach the trajectory:

big_u = mda.Universe(TPR, "big.xtc")
subsystem = big_u.select_atoms("(chainID A and (resid 63-71 or resid 205-208 or resid 250-256 or resid 370 or resid 417-435)) or ((chainID B and resid 476-504) or (chainID C and resid 601-624))")

u = mda.Merge(subsystem)
u.load_new("subsystem.xtc")

(I didn't try the above code so please check and use it as a starting point.)