search

dimchris / mdanalysis

Automatically exported from code.google.com/p/mdanalysis
0 stars 0 forks source link

Generating trajectories of Gromacs PCA's #79

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
This is a suggestion for a mdanalysis improvement: to be able to generate 
trajectories of Gromacs PCA's. 

PCA are trr files without forces and velocities (only positions, which 
represent PCA vectors). The mode index is numbered in "step".

These are the first few lines of gmxdump -f pca.trr...

   natoms=      1252  step=         1  time=3.2866406e+00  lambda=         0
   box (3x3):
      box[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   x (1252x3):
      x[    0]={ 2.32158e-02,  1.46866e-03, -7.98529e-03}
      x[    1]={ 2.24531e-02,  3.62402e-03, -8.81686e-03}
      x[    2]={ 1.76530e-02,  5.75658e-03, -6.89898e-03}
      x[    3]={ 1.40843e-02,  1.05042e-02, -5.93055e-03}
      x[    4]={ 7.69775e-03,  1.05042e-02, -4.27620e-03}
      x[    5]={ 5.82995e-03,  1.53628e-02, -2.67087e-03}

Presently, one can generate a trajectory from some given coordinates as:

# iterate over the new frames to write trajectory
for frame,coords in enumerate(newtraj):
        ts._pos[:] = coords     # change coords of current in-memory time step
        ts.frame = frame         # manually change the frame number
        W.write(ts)
W.close() 

But that puts zero velocities and does not include the mode index. That renders 
the trajectory unacceptable to Gromacs as a PCA trr file.

Original issue reported on code.google.com by rcreh...@gmail.com on 21 Sep 2011 at 10:06

GoogleCodeExporter commented 9 years ago 
Thanks, will have a look. Two questions:

- You say that the above recipe (ts.frame = frame) will not set the mode index. 
What does it set instead? Can you gmxdump the resulting trajectory (which still 
has v and f). If the problem is that step starts at 0 but modes start at 1: How 
about using the optional starting offset argument to enumerate():

  for mode, coords in enumerate(newtraj, 1):
      ts.frame = mode
      ...

  i.e. starting the enumeration at 1 instead of 0?

- Is the reported time the eigenvalue or is it some arbitrary number?

Thanks.

Original comment by orbeckst on 21 Sep 2011 at 2:07

GoogleCodeExporter commented 9 years ago 
GROMACS expect the mode number in step variable. 

This way, it works:
for frame,mode in enumerate(modes[4:16]):
    coords = mode.scaledToNorm(1.).array*10 #nm a angstroms
    ts.lmbda = -1
    if frame<=1:
       ts._pos[:] = xav
    else:
       ts._pos[:] = coords
    ts.frame = frame         # manually change the frame number
    ts.step = frame -1
    if frame <= 1:
       ts.time = frame-1
    else:
       ts.time = mode.frequency
    W.write(ts) # convert angstrom to nm for gmx

Original comment by rcreh...@gmail.com on 30 Sep 2011 at 9:42

GoogleCodeExporter commented 9 years ago 
Good that you figured out how to "reformat" the modes.

Do you think that there is anything in MDAnalysis that should be changed/added 
to facilitate this conversion? If so, add it here. Otherwise I'll close the 
report in a few days.

Original comment by orbeckst on 30 Sep 2011 at 3:11

GoogleCodeExporter commented 9 years ago 
I guess that the only thing that needs to be changed is the documentation, 
somehow stating that for a PCA trajectory the PC index goes into the ts.step. 
Just to avoid other people opening new issues :-)

Original comment by rcreh...@gmail.com on 30 Sep 2011 at 3:34

GoogleCodeExporter commented 9 years ago 
Ok, will do that and add your recipe to the docs. 

Good that you found ts.step; I just looked at the source again and saw that 
this is very specific to XTC/TRR and is not automatically set from frame --- 
sorry for misleading you in my original example.

Many thanks again for the feedback and the solution!

Original comment by orbeckst on 30 Sep 2011 at 4:38

GoogleCodeExporter commented 9 years ago 
added the recipe in r906, see 
http://mdanalysis.googlecode.com/svn/trunk/doc/html/documentation_pages/coordina
tes/TRR.html#filling-a-trr-with-pca-modes

Original comment by orbeckst on 5 Oct 2011 at 11:15