DCD Unit cell angles different between MDAnalysis and VMD

[GoogleCodeExporter]

Hello,

As discussed in the mailling list 
(https://groups.google.com/forum/#!msg/mdnalysis-discussion/C5TbkZLB2Eo/vyzcafrF
1AkJ), alpha and gamma angles are swapped between VMD and MDAnalysis when 
opening a DCD file (attached).

MDAnalysis
import MDAnalysis
u = MDAnalysis.Universe("sample.psf", "sample.dcd")
u.trajectory.ts.dimensions
  array([ 38.42659378,  38.39310074,  44.75979996,  60.02891541,  90. ,  90. ], dtype=float32)

VMD
vmd sample.psf sample.dcd
molinfo 0 get {a b c alpha beta gamma}  (Using TkConsole)
  38.426594 38.393101 44.759800 90.000000 90.000000 60.028915

The differences are caused by the indexing of the unitcell array when reading 
DCD files. In lines 732, 734, 738 and 740 of dcd.c 
(https://code.google.com/p/mdanalysis/source/browse/package/src/dcd/dcd.c#732), 
MDAnalysis associates alpha to index 1, while gamma is set to index 4. When 
looking at lines 01013, 01015, 01019 and 01021 of VMD's dcdplugin.c 
(http://www.ks.uiuc.edu/Research/vmd/plugins/doxygen/dcdplugin_8c-source.html#l0
0947), the same indexes are swapped (4 for alpha and 1 for gamma).
This problem became apparent when dealing with a non-orthorhombic box, since 
these exchanges make no difference for orthorhombic boxes.

Thanks

Original issue reported on code.google.com by [email protected] on 7 Jul 2014 at 4:44

Attachments:

• sample.tgz

[kain88-de]

@orbeckst

About the unitcell conversion for CHARMM > 22 do we really need to do the angle conversion in dcd.c. That is currently done for every dcd file.

[orbeckst]

Follow up at #1155 (review)

[kain88-de]

I'm not sure if the current libdcd port fixes the problems described here. I also don't have any experience with CHARM or NAMD to be able to test this. It would be nice if someone with experience in the programs could make a table like in #187 (comment) and add test files for us. Implementing the checks later on should be straight forward.

Also to check what is stored in the DCD the currently unmerged libdcd can be used. It returns the raw unit-cell information. That should make things a bit easier