DCDReader calculates wrong time (istart is in integrator steps not frames)

[orbeckst]

Expected behaviour

The DCDReader should correctly set the time step of a DCD file.

Actual behaviour

The starting frame is not zero but a large number (actually, it is nsavc to larger, see below)

Code to reproduce the behaviour

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD
u = mda.Universe(PSF,DCD)

times = [ts.time for ts in u.trajectory]
print(times[:10], "...", times[-10:])

gives

[999.9999119200186, 1000.9999118319386, 1001.9999117438587, 1002.9999116557786, 1003.9999115676986, 1004.9999114796186, 1005.9999113915387, 1006.9999113034587, 1007.9999112153787, 1008.9999111272988] ... [1087.9999041689803, 1088.9999040809003, 1089.9999039928202, 1090.9999039047402, 1091.9999038166602, 1092.9999037285804, 1093.9999036405004, 1094.9999035524204, 1095.9999034643404, 1096.9999033762604]

but it should be

[0.0 1.0 ... 97.0]

Note the header

print(u.trajectory._file.header)

contains istart = 1000

{'natoms': 3341, 'istart': 1000, 'nsavc': 1000, 'delta': 0.020454827696084976, 'is_periodic': False, 'remarks': '* DIMS ADK SEQUENCE FOR PORE PROGRAM                                            * WRITTEN BY LIZ DENNING (6.2008)                                               *  DATE:     6/ 6/ 8     17:23:56      CREATED BY USER: denniej0                '}

which is wrongly interpreted as frames instead of integrator steps in various places:

• wrong time calculation

  mdanalysis/package/MDAnalysis/coordinates/DCD.py

  Line 219 in e94b1c4

  ts.time = (ts.frame + self._file.header['istart']) * self.ts.dt

• docs

  mdanalysis/package/MDAnalysis/coordinates/DCD.py

  Line 384 in e94b1c4

  istart : int (optional)

• docs:

  mdanalysis/package/MDAnalysis/lib/formats/libdcd.pyx

  Line 428 in e94b1c4

  istart: starting frame number

Background

From readdcd.h:

/*
 * Read the header information from a dcd file.
 * Input: fd - a file struct opened for binary reading.
 * Output: 0 on success, negative error code on failure.
 * Side effects: *natoms set to number of atoms per frame
 *               *nsets set to number of frames in dcd file
 *               *istart set to starting timestep of dcd file
 *               *nsavc set to timesteps between dcd saves
 *               *delta set to value of trajectory timestep
 *               *nfixed set to number of fixed atoms 
 *               *freeind may be set to heap-allocated space
 *               *reverse set to one if reverse-endian, zero if not.
 *               *charmm set to internal code for handling charmm data.
 */

https://www.charmm.org/charmm/documentation/by-version/c42b1/params/doc/dynamc/

4. Query a trajectory file
   TRAJectory QUERy UNIT integer

The query command rewinds an open trajectory file and then reads the header information from this trajectory file. It prints a summary and sets the following command line substitution parameters:

• 'NFILE' - Number of frames in the trajectory file
• 'START' - Step number for the first frame
• 'SKIP' - Frequency at which frames were saved (NSTEP=NFILE*SKIP when not using restart files)
• 'NSTEP' - Total number of steps from the simulation
• 'NDEGF' - Number of degrees of freedom from the simulation (Can be use to get the temperature with velocity files).
• 'DELTA' - The dynamics step length (in picoseconds).

Currently version of MDAnalysis:

(run python -c "import MDAnalysis as mda; print(mda.__version__)")

develop

[orbeckst]

mdanalysis/package/MDAnalysis/coordinates/DCD.py

Line 219 in e94b1c4

ts.time = (ts.frame + self._file.header['istart']) * self.ts.dt

should be changed to

ts.time = (ts.frame + self._file.header['istart']/self._file.header['nsavc']) * self.ts.dt

[tylerjereddy]

Can we safely ignore the heterogeneity between DCD formats in applying the fix & writing appropriate tests?

[kain88-de]

That was my oversight. @orbeckst do you know if this is stored similarly in CHARMM and NAMD?

[orbeckst]

As far as I can tell, NAMD uses the same convention (NAMD 2.9: DCD). The comment in readdcd.h quoted above #1819 (comment) (which is from VMD/catdcd) also says

*istart set to starting timestep of dcd file
*nsavc set to timesteps between dcd saves

so given that nsavc is definitely integrator timesteps it follows that istart should also be measured in (integrator) timesteps.

My verdict is that it should also apply to NAMD DCDs.

Note that our "NAMD" trajectory ("PSF_NAMD_TRICLINIC", "DCD_NAMD_TRICLINIC") seems to have wrong header information (maybe because it is just a single frame or because it was written by VMD's DCD plugin?):

{'delta': 1.0,
 'is_periodic': True,
 'istart': 0,
 'natoms': 5545,
 'nsavc': 1,
 'remarks': 'Created by DCD pluginREMARKS Created 06 July, 2014 at 17:29Y5~CORD,'}

whereas a recent NAMD trajectory that I managed to look at looks as expected:

{'delta': 0.020454829558730125,
 'is_periodic': False,
 'istart': 1000,
 'natoms': 8822,
 'nsavc': 1000,
 'remarks': 'REMARKS FILENAME=protein.dcd CREATED BY NAMD                                  REMARKS DATE: 02/28/18 CREATED BY USER: unknown                                '}

Do we have a proper test NAMD trajectory in our collection? E.g., the same AdK system that we have for Gromacs but ran in NAMD? (@sseyler do you have one??)

[sseyler]

I have three TMD trajectories for the AdK closed-to-open transition that I produced in NAMD. @orbeckst you can find them on NFS here:

/nfs/homes4/sseyler/Simulations/adk/namd/TMD/4-run

with the names tmd0001.dcd, tmd0002.dcd, and tmd0003.dcd. Maybe one of these can be incorporated as a test trajectory in MDAnalysis?

[orbeckst]

@sseyler thanks, I found a ~4 MB trajectory that is native NAMD (the tmd0001.dcd etc were already processed with MDAnalysis):

>>> import MDAnalysis as mda
>>> PSFNAMD = "adk_closed.psf"
>>> DCDNAMD = "adk_gbis_tmd-fast1.dcd'"
>>> u = mda.Universe(PSFNAMD, DCDNAMD)
>>> u.trajectory._file.header
{'delta': 0.020454829558730125,
 'is_periodic': False,
 'istart': 10,
 'natoms': 3341,
 'nsavc': 10,
 'remarks': 'REMARKS FILENAME=adk_gbis_tmd-fast1.dcd CREATED BY NAMD                         REMARKS DATE: 07/21/15 CREATED BY USER: sseyler                                 '}

where the delta is 1 fs

>>> mda.units.convert(u.trajectory._file.header['delta'], 'AKMA', 'ps')
0.0010000000029814057

All of this agrees with the NAMD input parameters (for NAMD 2.10 for Linux-x86_64-multicore-CUDA)

# Integrator Parameters
timestep            1.0  ;# 1fs/step

# Output
dcdfreq                 10

and the reported output from NAMD

Info: TIMESTEP               1
Info: DCD FILENAME           adk_gbis_tmd-fast1.dcd
Info: DCD FREQUENCY          10
Info: DCD FIRST STEP         10

We could add trajectory and PSF as native NAMD trajectories.

[orbeckst]

I added tests and the fix; some of the other tests are failing now: need to look at them in more detail and check if the reference needs to be changed or code logic somewhere.