Reading TRR files

[ericsmoll]

Hello MDAnalysis users,

I am new to MDAnalysis and am trying to evaluate how easily I can read and write gromacs trr data.

from MDAnalysis.coordinates.TRR import *

a = TRRReader("A.trr")

last_frame = a[-1]

last_frame.has_x # => True

last_frame.has_v # => False

last_frame.has_f # => True

b[0] gives me a three element array of type float32 (I ran a single precision calculation) with the coordinates of the first atom in my simulation.

How do I access the forces? I examined the method available to "b":

copy

has_f

numatoms

time

copy_slice

has_v

prec

volume

dimensions

has_x

status

frame

lmbda

step

None of these seem to be able to access the forces.

Also, I assume that the coordinates and forces will be in double precision if I run a double precision simulation, correct? There is no TRRReader that I am missing, correct?

Best,

Eric

[Richard Gowers]

Hi Eric

To read a trajectory you want to use a Universe object, which itself handles the reader for you. You'll also need to give it some sort of topology file, like a gro or pdb file

import MDAnalysis as mda

u = mda.Universe(topologyfile, 'yourfile.trr')

You can then skip to the last frame using

u.trajectory[-1]

And then access the forces which are attached to the AtomGroup objects. 

u.atoms.forces

Whenever you query information like positions, velocities etc you usually do it via the atoms, ie select some atoms -> make an AtomGroup, for this AtomGroup, what are the forces.  The AtomGroup of everything is u.atoms, but you'll likely want to make subselections of this.  There's a good tutorial that's just been wrote:

http://orbeckst.github.io/MDAnalysisTutorial/basics.html#universe-and-atomgroup


Hope that helps
Richard

[Oliver Beckstein]

Hi Eric,

On 9 Nov, 2014, at 15:03, ericsmoll wrote:

I am new to MDAnalysis and am trying to evaluate how easily I can read and write gromacs trr data.

Richard already covered the important point: normally you access everything through a Universe, which makes trajectory data available in the context of atoms (see eg http://orbeckst.github.io/MDAnalysisTutorial/trajectories.html ) 

Just to give an example: calculate the average force on a protein atom and a solvent atom (using some of the test trajectories that contain forces and velocities): 

import MDAnalysis

from MDAnalysis.tests.datafiles import PDB_xvf, TRR_xvf

u = MDAnalysis.Universe(PDB_xvf, TRR_xvf)

p = u.selectAtoms("protein")

s = u.selectAtoms("not protein")

for ts in u.trajectory:

f_protein = p.forces.mean(axis=0)  

f_solvent = s.forces.mean(axis=0)

print("{0} {1} ps:  <F_prot> = {2} kJ/mol/A   <F_solv> = {3} kJ/mol/A".format(ts.frame, u.trajectory.time, f_protein, f_solvent))

For me this prints

1 50.0 ps:  <F_prot> = [ 0.40944228 -0.25576955 -0.04409491] kJ/mol/A   <F_solv> = [-0.02028263  0.01273486  0.00221454] kJ/mol/A

2 100.0 ps:  <F_prot> = [ 0.48914304 -0.07478244  0.01495386] kJ/mol/A   <F_solv> = [-0.02445253  0.00366654 -0.00072841] kJ/mol/A

3 150.0 ps:  <F_prot> = [ 0.13971081  0.14064492 -0.28849533] kJ/mol/A   <F_solv> = [-0.00698189 -0.00691999  0.01438792] kJ/mol/A

Note that MDAnalysis converts to its native units (i.e. kJ/(mol Å) for forces and ps for time).

If you REALLY only want the trajectory you CAN use a TRRReader on its own:

import MDAnalysis

from MDAnalysis.tests.datafiles import PDB_xvf, TRR_xvf

trr = MDAnalysis.coordinates.TRR.TRRReader(TRR_xvf)

forces = trr.ts._forces

You get the forces from the private _forces attribute of a Timestep object. You then have to slice your forces array using your own atom indices, e.g. for the forces of the atoms 0-3: forces[0:4]. The units are the original Gromacs units.

Also, I assume that the coordinates and forces will be in double precision if I run a double precision simulation, correct?

No, the reader will always use float32. If you need float64 you will have to add additional C-code together with Python code (and the appropriate SWIG interface) that does everything in double. It's not really hard but a bit tedious and will take a while if you're new to the code.

There is no TRRReader that I am missing, correct?

Yes, you're not missing one.

Ask if you have further questions!

Oliver

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst  * orbe...@gmail.com

[Eric Smoll]

Hello Oliver,

I see. I would like double precision coordinates, forces, etc. I started doing some digging and it seems that gromacs has a C library with the functions necessary to read and write trr files. When I look at the headers, the exported functions ("read_trr" and "write_trr") both return floats for the coordinates, forces, etc. Is this what you mean when you say I will have to add additional C-code? I need to implement versions of these functions that returns double precision values and then write a python wrapper using (for example) the ctypes module?

Best,

Eric

--
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.
To post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at http://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

[Oliver Beckstein]

Hi Eric,

On 11 Nov, 2014, at 20:43, Eric Smoll wrote:

> I see. I would like double precision coordinates, forces, etc. I started doing some digging and it seems that gromacs has a C library with the functions necessary to read and write trr files.

We based our XTC and TRR reader on libxdr but Manel added some significant improvements to what we now call libxdr2 (namely, fast seeking!). The code is under package/src/xdrfile in the distribution and you can browse it at https://code.google.com/p/mdanalysis/source/browse/package#package%2Fsrc%2Fxdrfile

> When I look at the headers, the exported functions ("read_trr" and "write_trr") both return floats for the coordinates, forces, etc. Is this what you mean when you say I will have to add additional C-code? I need to implement versions of these functions that returns double precision values and then write a python wrapper using (for example) the ctypes module?

If you look at xdrfile_trr.h you see that coordinates are of type *rvec. In xdrfile.h,

typedef float rvec[DIM]

You would probably want to turn this into

typedef double rvec[DIM]

and also check that the rest of the code is compatible. You could then duplicate the exported routines and export them via the SWIG interface https://code.google.com/p/mdanalysis/source/browse/package/src/xdrfile/libxdrfile2.i (which in my opinion makes it easy to create very pythonic wrappers).

Within MDAnalysis we actually assume that coordinate arrays are float32 or at least other cython/C-based code will only accept those but if you do not plan on using e.g. distance_array() then you should be fine.

If you come up with a solution that can be seamlessly integrated into MDAnalysis we're happy to consider it for a future release!

[Richard Gowers]

Changing the distance subroutines to allow double precision would be pretty trivial, I could get that done if it was required.

Maybe for 1.1 we could implement double precision everywhere...

[ericsmoll]

Hello Richard and Oliver,

Thank you for the guidance. When I load a gro / trr pair I can access a matrix of coordinates / velocities / forces with
u.atoms.get_positions() / u.atoms.get_velocities() / u.atoms.get_forces().

Is the cell size loaded as well? I cannot find it in the methods of u.atoms.

Best,
Eric

[ericsmoll]

Hello Richard and Oliver,

I believe I found it. The dimensions method has the cell information.

I have another question.
if I run len(list(u.atoms[:])) I get a series of entries like this:
< Atom 19184: name 'C9' of type 'C' of resname 'MOL', resid 1280 and segid 'SYSTEM'>

u.atoms.names() will yield an array of the name entries. There are similar methods for most of the other printed info except "type." How do I access the type parameter? Is the atomic-number stored anywhere?

Best,
Eric


On Thursday, November 13, 2014 1:05:32 PM UTC-7, Richard Gowers wrote:

[Oliver Beckstein]

Hi Eric,

On 16 Nov, 2014, at 13:44, ericsmoll wrote:

> I believe I found it. The dimensions method has the cell information.

yes, Timestep.dimensions()

>
> I have another question.
> if I run len(list(u.atoms[:])) I get a series of entries like this:
> < Atom 19184: name 'C9' of type 'C' of resname 'MOL', resid 1280 and segid 'SYSTEM'>
>
> u.atoms.names() will yield an array of the name entries. There are similar methods for most of the other printed info except "type."

u.atoms.types()

... is going to be added in the next release - it's in already in devel ;-)

In the meantime:

> How do I access the type parameter? Is the atomic-number stored anywhere?

Yes, Atom.type.

E.g.

atypes = numpy.array([a.type for a in u.atoms])

[Eric Smoll]

Hello Oliver,

Thanks again for the guidance. I followed your suggestions and made the xdr library changes necessary to read and write trr files in double precision. I tested it and it appears to work. The numbers loaded into python match a full precision dump of the trr file using the double precision gromacs utility "g_traj_d -fp"

Can you provide some input on appending to a trr file? Is this possible in MDAnalysis? If not, does the underlying xdrfile library support appending to a trr file? If not, what would need to be done?

Best,

Eric

[Oliver Beckstein]

Hi Eric,

On 17 Nov, 2014, at 08:56, Eric Smoll wrote:

> Thanks again for the guidance. I followed your suggestions and made the xdr library changes necessary to read and write trr files in double precision. I tested it and it appears to work. The numbers loaded into python match a full precision dump of the trr file using the double precision gromacs utility "g_traj_d -fp"

Great.

> Can you provide some input on appending to a trr file? Is this possible in MDAnalysis? If not, does the underlying xdrfile library support appending to a trr file? If not, what would need to be done?

I don't think the xdrlibrary supports this but I am not 100% sure. Maybe Manel can chime in. Of the top of my head I don't know what needs to be done but I assume it's possible because frames are independent in a TRR. You'd have to figure out how appending works in XDR files (which is the TRR file format).

I assume you don't just want to concatenate files?

At least I should mention that in MDAnalysis you can also read a sequence of files and concatenate them on the fly:

u = MDAnalysis.Universe("system.tpr", ["01.trr", "02.trr", "03.trr"])

(or even u = MDAnalysis.Universe("system.tpr", "01.trr", "02.trr", "03.trr") )