CHARMM Velocity reader

[hesther]

Currently, there is no possibility to read in CHARMM velocity files. These are binaries, very similar to CHARMM dcds, just with a slightly different header, and the data block contains velocities saved in Angstrom/AKMA instead of positions in Angstrom. Would it be possible to add a CHARMM velocity reader?

I suggest to add a new reader class (eg. VelReader) similar to DCDReader (and a new timestep subclass with changed init function to take in velocities instead of positions) in MDAnalysis.coordinates.DCD.py, as well as enabling this new reader in MDAnalysis.coordinates.init.py. I am not quite sure if there are any other changes required.

Currently, I reformat CHARMM velocity files to CHARMM dcd files (which is tedious, and time consuming), create a Universe, and get the velocities as the "positions" times MDAnalysis.units.timeUnit_factor['AKMA'].

I would be very thankful if a CHARMM velocity reader would be included in a future release! Thanks!

[orbeckst]

@hesther , a velocity reader would be a good addition.

However, all developers are really quite busy and I don't see anyone going to implement it any time soon because we have a whole bunch urgent fixes/changes before a 1.0 release that we are focusing our effort on.

Can you try to do it yourself? We are happy to review pull requests and give advice. It might not be too hard if you start with the code in package/MDAnalysis/coordinates/DCD.py. Something along the lines of

class DCDVelocityReader(DCDReader):
    format = 'VDCD'   # we need another format name to distinguish from positions
    flavor = 'CHARMM'
    units = {'time': 'AKMA', 'velocity': 'Angstrom/AKMA'}

    def __init__(self, filename, convert_units=True, dt=None, **kwargs):
         #  need to use different Timestep with velocities --- might be a bit more
         # difficult than just calling the __init__ of the parent class
         super(DCDVelocityReader, self).__init__(filename, convert_units=True, dt=None, **kwargs)

    def _frame_to_ts(...):
            ...
            ts.velocities = frame.xyz
            if self.convert_units:
                self.convert_pos_from_native(ts.dimensions[:3])
                self.convert_pos_from_native(ts.positions)
            return ts

The new reader should be registered as "VDCD".

What I wrote above is not complete but contains some of the ideas. Have a look at the Gromacs XDR (TRR) reader for how we deal with velocities.

[hesther]

Thanks for the reply! If I find the time, I will try, but am very busy right now with finishing my PhD-thesis and preparing for my next PostDoc position ;)

[ayylemao]

It seems that this feature is still not in the package and I still think it would be a good addition to the MDAnalysis package. Are there plans to get this into a future release?

[orbeckst]

There are currently no plans that I’m aware of but MDAnalysis is open source so contributions are welcome.

1. How would you expect MDAnalysis to behave if a VDCD reader existed? Can you write the code that you would like to use, ie, how would you want to load your Universe?
2. Can you provide a link to the format description? Which program produces those files?

MDAnalysis aims to implement standard-compliant file I/O so knowing the exact definition of file formats and how commonly they are used is important.