fix DCDReader istart (#1819) #1832

orbeckst · 2018-03-15T07:25:52Z

Changes made in this Pull Request:

add native NAMD 2.10 test DCD trajectory with PSF
add tests
add fix

PR Checklist

Tests?
Docs?
CHANGELOG updated?
Issue raised/referenced?

orbeckst · 2018-03-15T08:10:39Z

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.

Help welcome.

kain88-de · 2018-03-15T08:29:50Z

package/MDAnalysis/coordinates/DCD.py

        """convert a dcd-frame to a :class:`TimeStep`"""
        ts.frame = self._frame
-        ts.time = (ts.frame + self._file.header['istart']) * self.ts.dt
+        ts.time = (ts.frame + self._file.header['istart']/self._file.header['nsavc']) * self.ts.dt


I thought there where more places where we need changes. I'm pretty sure the writer needs to be updated too.

I am pretty sure that you're right... I just committed what I had before falling asleep.

I think I fixed all occurences of istart/nsavc in the reader and writer. I also added a test for the writer to check the defaults that we are setting. Namely, istart=0 is not what CHARMM would do: it would set it to nsavc which is what MDAnalysis now does for istart=None. However, I kept the MDAnalysis value (which we had for a while) as the default.

orbeckst · 2018-03-15T16:42:11Z

#1830 needs to be fixed before the tests will be useful

kain88-de · 2018-03-15T16:43:25Z

you can rebase. We already applied a band aid. The fix for #1830 needs some more time

orbeckst · 2018-03-15T16:58:34Z

ok

kain88-de · 2018-03-18T08:37:26Z

testsuite/MDAnalysisTests/analysis/test_rms.py

    @pytest.fixture()
    def correct_values(self):
-        return [[0, 1000, 0], [49, 1049, 4.68953]]
+        return [[0, 1, 0], [49, 50, 4.68953]]


@orbeckst this makes the tests pass. I'm unsure why we have to adjust the step by 1 position now.

Do you mean, adjusting the step by 1 position instead of 1000? The reason for this is that istart is measured in integrator time steps. 1 frame contains nsavc integrator time steps. If a DCD wants to state that the trajectory starts "at frame 1" then it has to set istart = nsavc. In 0.17.0 we wrongly interpreted istart as frames instead of integrator time steps.

If your question is "why do our test DCDs now have a starting time "at frame 1", e.g.

>>> import MDAnalysis as mda ... from MDAnalysisTests.datafiles import PSF, DCD, DCD_NAMD_GBIS, PSF_NAMD_GBIS, TPR, XTC ... >>> # DCD time starts at 1*dt >>> u = mda.Universe(PSF, DCD) >>> print(u.trajectory.frame, u.trajectory.time, u.trajectory.n_frames, u.trajectory.totaltime) 0 0.9999999119200186 98 96.9999914562418 >>> u = mda.Universe(PSF_NAMD_GBIS, DCD_NAMD_GBIS) >>> print(u.trajectory.frame, u.trajectory.time, u.trajectory.n_frames, u.trajectory.totaltime) 0 0.010000000029814058 100 0.9900000029515917 >>> # XTC starts at 0*dt >>> u = mda.Universe(TPR, XTC) >>> print(u.trajectory.frame, u.trajectory.time, u.trajectory.n_frames, u.trajectory.totaltime) 0 0.0 10 900.0000686645508

then the answer is because this is how we interpret the meaning of istart at the moment: as an offset of istart/nsavc frames from 0 (based on the docs that I cited in #1819). I am admittedly not 100% sure if this is the best way to do this; it seems to imply that the starting frame is not saved in a trajectory and the first frame that was saved would be the one istart integrator time steps after the begin of the simulation. The CHARMM docs (Example under TRAJECTORY command) imply that the DYNAmics command does indeed start saving after the coordinates with the initial conditions:

Example. Assume that the three files traj1.trj, traj2.trj and traj3.trj were created using the following dyanmics commands:
dyna start nstep 1000 nsavc 100 ! saves 10 frames (100, 200, ...1000)
dyna restart nstep 10000 nsavc 100 ! saves 100 frames (1100, 1200, ...10100)
dyna restart nstep 10000 nsavc 100 ! saves 100 frames (10200, 10300, ...20100)

I might have to play around with the CHARMM TRAJectory command to get a better idea of what CHARMM really thinks of a trajectory.

However, I'd like to fix istart first (which is clearly wrong) and then discuss what is more of a discussion about how to interpret quirks in the DCD format and how to make it work seamlessly in MDAnalysis.

The change of 1000 to 1 is OK. We made this change in the current development cycle. I'm curious why we needed to change from step 49 to 50. My first try to fix the istart behavior was in #1767. The change in the iteration introduced by your changes is new to this pr. Because such a change would break user code I'm hesitant to merge this.

There are some NAMD experts in my group. I will ask them what behavior is best for them.

Can you explain with a simple example what the change "from 49 to 50" means?

I quickly browsed through the discussion in #1767. The choice istart=1 is now istart=None or manually setting istart=nsavc. Setting istart=0 as the default for the writer is retained.

If we wanted to make clearer that we want the user to be able to set frame offsets then we should introduce something like frame_offset. We already have offset for the Timestep.

I haven't understood yet where changes in frames come from. It might well be that something needs to be fixed but I'd like to have a test case for this in the dcd tests (and not primarily in the rms code).

I will have some time during the week to have a closer look at this.

I'm curious why we needed to change from step 49 to 50.

I think the same reason: before we had 1049, which was calculated as istart + frames = 1000 + 49. Now we have istart/nsavc + frames = 1000/1000 + 49 = 1 + 49 = 50.

Sorry I didn't have time yet to look at this. But your guess is likely right. This means from now on all reported time will start from 1? except when istart is 0.

I'm cool with that. This will affect users by showing later times now. I would change the CHANGELOG message to say that we used to report to early timestamps and now read times correctly and users should expect that times are shifted to later values. I like to be a bit more exact then just we did it wrong

Otherwise looks good to me. It can be merged.

orbeckst · 2018-03-18T18:17:54Z

Last commit harmonizes use of decimal and should make the lone failure on travis pass now.

orbeckst · 2018-03-18T18:19:47Z

@kain88-de I am travelling for most of the day and will be busy over the next few days. Please move forward with this PR as you see fit.

orbeckst · 2018-03-18T18:25:17Z

testsuite/MDAnalysisTests/coordinates/test_dcd.py

+@pytest.mark.parametrize("variable, default", (("istart", 0), ("nsavc", 1)))
+def test_DCDWriter_default(tmpdir, variable, default):
+    outfile = str(tmpdir.join("test.dcd"))
+    with mda.Writer(outfile, n_atoms=10) as W:


In principle, writing the empty DCD and returning the header could be made a module-level fixture. I am out of time to do this right now.

orbeckst · 2018-03-18T18:25:46Z

testsuite/MDAnalysisTests/coordinates/test_dcd.py

                            ts_dcd.positions,
                            3)
+
+@pytest.fixture(params=[(PSF, DCD),


Should this be declared with scope="module"?

orbeckst

comments inline – don't have time to address them

@sseyler

- 100 steps of TMD in GBIS implicit solvent from closed AdK to open AdK; dt = 1fs; output frequency 10, first step 10 - NAMD version 2.10 - created by @sseyler - see #1819 (comment) for details

- fixes #1819 - add tests for the istart and nsavc defaults in DCDWriter - update docs for nsavc and istart - adjust test_writer_dt() so that we get the same times for dcd and xtc (which is a bit tricky) - update CHANGELOG

orbeckst · 2018-03-19T00:11:59Z

For some reason, the minimal tests failed with something related to duecredit (https://travis-ci.org/MDAnalysis/mdanalysis/jobs/355076569). I rebased this branch against develop to make sure that everything is in here and force-pushed.

Apologies @kain88-de, you will have to make sure that your branch pulls properly.

kain88-de · 2018-03-19T08:17:45Z

Apologies @kain88-de, you will have to make sure that your branch pulls properly.

git reset --hard upstream/branch-name works like a charm everytime ;-)

orbeckst · 2018-03-22T17:42:48Z

@kain88-de any more comments?

tylerjereddy · 2018-03-22T22:34:23Z

I tried resolving the merge conflict in the github GUI since it was a pretty simple thing in the changelog.

orbeckst · 2018-03-23T00:00:28Z

still need to modify CHANGELOG following #1832 (comment)

- fixes #1819 - see PR #1832 and #1767 for discussions on DCD and istart; see also https://github.com/MDAnalysis/mdanalysis/wiki/FileFormats

orbeckst added the Work in progress label Mar 15, 2018

orbeckst added the help wanted label Mar 15, 2018

kain88-de requested changes Mar 15, 2018

View reviewed changes

orbeckst force-pushed the issue-1819-dcd-istart branch 2 times, most recently from 9fd0736 to 29f89f5 Compare March 18, 2018 07:13

orbeckst removed Work in progress help wanted labels Mar 18, 2018

kain88-de reviewed Mar 18, 2018

View reviewed changes

orbeckst commented Mar 18, 2018

View reviewed changes

orbeckst and others added 6 commits March 18, 2018 19:10

add native NAMD 2.10 PSF and DCD

2ea22fd

- 100 steps of TMD in GBIS implicit solvent from closed AdK to open AdK; dt = 1fs; output frequency 10, first step 10 - NAMD version 2.10 - created by @sseyler - see #1819 (comment) for details

add tests for #1819

b23e688

correctly calculate start of time in DCD

81a44aa

fixing istart/nsavc for writers

486b27b

- fixes #1819 - add tests for the istart and nsavc defaults in DCDWriter - update docs for nsavc and istart - adjust test_writer_dt() so that we get the same times for dcd and xtc (which is a bit tricky) - update CHANGELOG

fix rms tests

83648d2

DCD test: use decimal=5 for time tests

4e7d4aa

orbeckst force-pushed the issue-1819-dcd-istart branch from ea153c9 to 4e7d4aa Compare March 19, 2018 00:11

kain88-de approved these changes Mar 22, 2018

View reviewed changes

Merge branch 'develop' into issue-1819-dcd-istart

994f08c

updated CHANGELOG for #1819

6412516

- fixes #1819 - see PR #1832 and #1767 for discussions on DCD and istart; see also https://github.com/MDAnalysis/mdanalysis/wiki/FileFormats

orbeckst self-assigned this Mar 23, 2018

orbeckst merged commit 5a232a3 into develop Mar 23, 2018

orbeckst deleted the issue-1819-dcd-istart branch March 23, 2018 06:04

fix DCDReader istart (#1819) #1832

fix DCDReader istart (#1819) #1832

Uh oh!

Conversation

orbeckst commented Mar 15, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR Checklist

Uh oh!

orbeckst commented Mar 15, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

orbeckst commented Mar 15, 2018

Uh oh!

kain88-de commented Mar 15, 2018

Uh oh!

orbeckst commented Mar 15, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kain88-de Mar 18, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

orbeckst commented Mar 18, 2018

Uh oh!

orbeckst commented Mar 18, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

orbeckst left a comment

Choose a reason for hiding this comment

Uh oh!

orbeckst commented Mar 19, 2018

Uh oh!

kain88-de commented Mar 19, 2018

Uh oh!

orbeckst commented Mar 22, 2018

Uh oh!

tylerjereddy commented Mar 22, 2018

Uh oh!

orbeckst commented Mar 23, 2018

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

orbeckst commented Mar 15, 2018 •

edited

Loading

kain88-de Mar 18, 2018 •

edited

Loading