Incorrect box dimensions when writing .rst7 from .dcd

[tuckerburgin]

I'm getting some weird behavior from pytraj when attempting to write out a frame from a .dcd trajectory produced by CP2K to an .rst7 file. This is best demonstrated by the following comparison to the behavior from mdtraj, which is correct:

tburgin$ python
>>> import pytraj
>>> pytraj.__version__
'2.0.6'
>>> import mdtraj
>>> ptraj = pytraj.load('2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd', 'Cl-CCCl.prmtop')
>>> mtraj = mdtraj.load('2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd', top='Cl-CCCl.prmtop')
>>> pytraj.write_traj('pytraj.rst7', ptraj, frame_indices=[0])
>>> mtraj[0].save_amberrst7('mdtraj.rst7')
>>> quit()
tburgin$ tail -n 1 pytraj.rst7.1 
 129.7251139 137.7776043 129.5068450  32.7895050  41.8777299  32.3208509
tburgin$ tail -n 1 mdtraj.rst7 
  25.0720005  52.7510033  23.9169998  90.0000000  90.0000000  90.0000000

I'm not sure what's going on here, but the box info produced by mdtraj is the same as the box dimensions from the simulation that produced the trajectory. Here are the relevant files:

pytraj_report.tar.gz

Thanks!

[tuckerburgin]

Maybe this is related to #1604? For practical reasons I am still using a standalone pytraj release, I'm having trouble linking the ambertools version into my workflow (for reasons that are more to do with me and with that workflow than with pytraj, I think).

[drroe]

Hi. The problem here is the DCD trajectory (which based on the info in the DCD itself appears to have been written by CP2K) you have does not properly conform to the Charmm DCD standard. DCD trajectories from Charmm version 22 and up are supposed to store symmetric shape matrix information, not unit cell parameters. You can see that cpptraj (the engine underlying pytraj) expects this:

> trajin 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz 1 1
  [trajin 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz 1 1]
        Reading '2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz' as Charmm DCD
Warning: 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz: Reported number of frames in DCD file is 0,
Warning:        actual number of frames is 191. Only reading 191 frames.
        Version >= 22; assuming shape matrix is stored.

(Also note that CP2K is also not properly setting the number of frames in the DCD header).

CP2K is setting the DCD version bit in the header to 24, so it should be converting the unit cell parameters to a symmetric shape matrix, but it actually seems to be storing unit cell parameters instead. Since cpptraj expects symmetric shape matrix info in the DCD based on the version, the box ends up skewed and you get this warning:

Warning: Box is too skewed to perform accurate imaging.
Warning:  Images and imaged distances may not be the absolute minimum.

If you force reading unit cell information with the ucell keyword, things now work:

> trajin 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz 1 1 ucell
  [trajin 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz 1 1 ucell]
        Reading '2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz' as Charmm DCD
Warning: CHARMM version is >= 22 but 'ucell' specified.
Warning: Assuming box info is stored as unit cell and not shape matrix.
Warning: 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz: Reported number of frames in DCD file is 0,
Warning:        actual number of frames is 191. Only reading 191 frames.
        Version < 22; assuming X-aligned cell.

The following input:

parm Cl-CCCl.prmtop.gz
trajin 2.1.heat_interface_frame_2032.pdb_0_ts_guess.dcd.gz 1 1 ucell
trajout temp.rst7

Results in "good" box info:

$ tail -n 1 temp.rst7
  25.0720000  52.7510023  23.9170003  90.0000000  90.0000000  90.0000000

So you need to force cpptraj/pytraj to read unit cell info instead of symmetric shape info; @hainm is there a way to pass trajin-type arguments to pytraj.io?

Also, you should contact the CP2K devs and have them fix their usage of the DCD format. Probably the easiest thing would be for them to change the version bit in the DCD header to something less than 22 (in cpptraj I just use 21).

[tuckerburgin]

I see! Thanks so much for the comprehensive answer. I have alerted a CP2K dev to this issue, hopefully it can be fixed on their end.

[hainm]

@hainm is there a way to pass trajin-type arguments to pytraj.io?

hi @drroe: pytraj has not supported yet. I've filed #1621 to keep track.