search

dimchris / mdanalysis

Automatically exported from code.google.com/p/mdanalysis
0 stars 0 forks source link

DCD info: trajectory.delta wrong, trajectory.dt in wrong units #84

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
What steps will reproduce the problem?
1. u = Universe("brown.psf", "brown.dcd")
2. print u.trajectory.delta
   9.99999974738e-05
3. print u.trajectory.dt
  0.0048888208765
4. print u.trajectory.units["time"]
  AKMA
5.print 
u.trajectory.dt*MDAnalysis.core.units.timeUnit_factor[u.trajectory.units["time"]
]/u.trajectory.skip_timestep
  0.0999999974738

What is the expected output? What do you see instead?

for 2:
Expected: 0.0001
Instead: 9.99999974738e-05
Relative Error: 2.52621248364e-08

for 3:
Expected: 0.1
Instead: 0.0048888208765

for 4:
I do not know what I should expect. No units?

for 5:
Expected: 0.1
Instead: 0.0999999974738
Relative Error: 2.52621247254e-08

What version of the product are you using? On what operating system?
MDAnalysis 0.7.5 and MDAnalysis 0.7.4
Linux ubuntu 11.04 AMD 64 2.6.38-13
python 2.7

Please provide any additional information below.

Concerning 4:
The .lam file is in lj units, so the output should be in dimensionless
time units.
Here, the system seems to believes that the data is in AKMA units and converts 
internally to some other time unit.
But even when only considering the values for .delta, I still have a relative 
error of 2.e-8 even for simple values such as 0.0001.
Note that the error is zero for values of delta = 1.

This error is a pain in the ass, for long timeseries this adds up to more than 
one frame.

Original issue reported on code.google.com by bjrnfrd...@gmail.com on 14 Nov 2011 at 3:36

Attachments:

GoogleCodeExporter commented 9 years ago 
Thanks; I think this is related to Issue 67.

There seem to be two issues:

1) Rounding errors.

2) Assumption that any DCD records the time in AKMA.

Ad 1) I need to check the conversion factors. Assuming there's no mistake 
there, would it be acceptable to round the dt to, say, 6 decimal digits? 

Ad 2) The fact is that CHARMM "DCD" trajectories are defined to be in AKMA 
units and it is really LAMMPS fault if it uses the format but records a 
different time unit. So we'll need a way to make MDAnalysis pick up on this 
behaviour. We can try to sniff the comment field of the LAMMPS DCD (but this 
will still not work for LJ reduced units). We'll probably have to use the 
'format' keyword argument for Universe to un-ambiguously specify that this is a 
"LAMMPS_DCD" or "LJ_DCD" or something similar. Any suggestions?

Original comment by orbeckst on 14 Nov 2011 at 4:15

GoogleCodeExporter commented 9 years ago 
ROUNDING --- comments welcome!

(addressing (2) and (5) in the report)

I can't find a problem with the conversion factors in units, and in any case, 
the problem seems to be already that a delta of 0.0001 is expected to be read 
from the DCD but 9.99999974738e-05 is returned. delta is whatever is found in 
the dcd_header:

  print u.trajectory._dcd_header()

  {'charmm': 5,
   'delta': 9.9999997473787516e-05,
   'file_desc': 28,
   'first': 0,
   'fixedcoords_ptr': 0,
   'freeind_ptr': 0,
   'header_size': 276,
   'istart': 0,
   'natoms': 1,
   'nfixed': 0,
   'nsavc': 1000,
   'nsets': 31,
   'reverse': 0,
   'setsread': 21,
   'with_unitcell': 1}

DCDs store the dt as a 32-bit float so I am not entirely surprised to see 
something slightly different from 0.0001 (although I suppose that a relative 
error of 2.5e-8 is a little bit higher than what I'd naively expect, 2^-31 = 
4.6e-10).

Maybe this is due to some conversions between the C-code and numpy; if anyone 
can shed more light on this I am happy to fix it.

The number of time steps (in trajectory.frame or ts.frame) is an integer and 
will be exact. You could always compute the time as

  round(u.trajectory.delta, 6) * u.trajectory.skip_timestep * u.trajectory.frame

At some point we did round dt to 6 decimals or so but that didn't really fix 
the problem: once you have a timestep that's not near e.g. 0.001 you again make 
things worse.  

Would it be helpful if one could tell a trajectory reader that it should round 
dt to make it *appear* more exact?

For the rounding errors I cannot think of a way to really fix this, mainly 
because that's an inherent problem with floating point arithmetic, see 
http://docs.python.org/tutorial/floatingpoint.html as far as I understand it.

Comments and discussion welcome!

Original comment by orbeckst on 17 Nov 2011 at 3:30

GoogleCodeExporter commented 9 years ago 
DIFFERENT TIME UNIT in LAMMPS DCD

(addressing (3) in the report)

I am attaching a patch for a "LAMMPS DCD" Reader/Writer. By default it assumes 
that the dcd's timestep is in ps. However, this can be changed with the 
timeunit keyword (there's also lengthunit):

  >>> u1 = MDAnalysis.Universe("brown.psf", "brown.dcd", format="LAMMPS")
  >>> print u1.trajectory.delta,  u1.trajectory.dt   
  9.99999974738e-05 0.0999999974738     # ps, note that skip_timestep = 1000

  >>> u2 = MDAnalysis.Universe("brown.psf", "brown.dcd", format="LAMMPS", timeunit="fs")
  >>> print u2.trajectory.delta,  u2.trajectory.dt  
  9.99999974738e-05 9.99999974738e-05

  >>> u3 = MDAnalysis.Universe("brown.psf", "brown.dcd", format="LAMMPS", timeunit="ns")
  >>> print u3.trajectory.delta,  u3.trajectory.dt  
  9.99999974738e-05 99.9999974738

When reading as a normal (i.e. CHARMM) DCD then the previous behaviour is 
retained:

  >>> u = MDAnalysis.Universe("brown.psf", "brown.dcd")
  >>> print u.trajectory.delta,  u.trajectory.dt
  9.99999974738e-05 0.0048888208765

Would this be a useful addition?

Original comment by orbeckst on 17 Nov 2011 at 3:51

Attachments:

GoogleCodeExporter commented 9 years ago 
> I already have one question:
> Say I have an external source that knows the delta value exactly, can specify 
it after initialization of the Universe?
> Like:
> u.trajectory.delta=known_delta

I think that this should work; the header is only read once.

> If I specify it afterwards, is there a way to re-initialize the internal 
timestep > (dt) such that this reflects the changes in delta?
> Something like 
> u.trajectory.reset_time_values_with_correct_delta(known_delta)

dt = skip_timestep * delta

(or rather convert_time_from_native(self.delta)) so I would think that won't be 
necessary. Try it out by setting delta. If there's still a problem we might 
need to look closely at convert_time_from_native().

Original comment by orbeckst on 18 Nov 2011 at 3:57

GoogleCodeExporter commented 9 years ago 
I added the patch from Comment 3 in r925: set the format argument of Universe: 

  u = Universe(PSF, 'lammps.dcd', format='LAMMPS', timeunit='ps') 

to force the time unit. Other supported time units are fs, ns, AKMA.

(Note that there is also a new keyword 'lengthunit' to tell the LAMMPS DCD 
Reader how interpret coordinates, e.g. lengthunit='nm'.)

I am still not sure how to address the other concerns. I will need some more 
feedback on this.

Original comment by orbeckst on 28 Nov 2011 at 6:35

GoogleCodeExporter commented 9 years ago 
Has anyone got any more input on this issue?

Bascially, without more feedback I don't really know what else to do here 
except to close the issue even though not all points have been fully addressed.

Original comment by orbeckst on 13 Feb 2012 at 3:19

GoogleCodeExporter commented 9 years ago 
The workaround with setting 

  u.trajectory.delta = known_delta 

works anytime. As soon as you query the time step with

  print u.trajectory.dt

it will use the updated delta to compute dt. 

Example:

  >>> import MDAnalysis; from MDAnalysis.tests.datafiles import *
  >>> u = MDAnalysis.Universe(PSF,DCD)
  >>> print u.trajectory.delta
  0.0204548276961             # time step in native AKMA units
  >>> u.trajectory[10]
  < Timestep 11 with unit cell dimensions array([  0.,   0.,   0.,  90.,  90.,  90.], dtype=float32) >
  >>> print u.trajectory.dt   # time step converted from AKMA to ps
  0.99999991192

  # Now change the time step
  >>> u.trajectory.delta *= 1000
  >>> print u.trajectory.dt   # new time step (converted from AKMA to ps)
  999.99991192  

I think that this is pretty transparent and works the way one would expect it 
to work hence I don't see the need to implement an additional reset method.

Therefore I am closing this issue (using Fixed instead of WontFix because half 
of it resulted in code changes.) Please feel free to reopen it.

Thanks,
Oliver

Original comment by orbeckst on 29 Feb 2012 at 12:46