search

MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.3k stars 648 forks source link

Drop TRZ reading forces, and change DL_Poly units #2393

Open lilyminium opened 4 years ago

lilyminium commented 4 years ago

MDAnalysis is supposed to have force units of kJ/mol.A.

DL_Poly

The DLPoly format seems to work in 10 J/mol.A :

Last but not least, it is worth pointing out that composite entities, such as velocities and forces, have their units expressed as composites of the default DL POLY units as shown in Section 1.3.7.

(manual: ftp://ftp.dl.ac.uk/ccp5/DL_POLY/DL_POLY_4.0/DOCUMENTS/USRMAN4.pdf)

Section 1.3.7:

Internally all DL POLY 4 subroutines and functions assume the use of the following defined molecular units: • The unit of time (to) is 1 × 10−12 seconds (i.e. picoseconds) • The unit of length (lo) is 1 × 10−10 metres (i.e. ̊Angstroms) • The unit of mass (mo) is 1.6605402 × 10−27 kilograms (i.e. Daltons - atomic mass units) • The unit of charge (qo) is 1.60217733 × 10−19 Coulombs (i.e. electrons - units of proton charge) • The unit of energy (Eo = mo(lo/to)2) is 1.6605402 × 10−23 Joules (10 J mol−1) • The unit of pressure (P = E l−3) is 1.6605402 × 107 Pascals (163.882576 atmospheres) ooo • Planck’s constant ( ̄h) which is 6.350780668 × Eoto .

But the ConfigReader and HistoryReader don't appear to convert units.

https://github.com/MDAnalysis/mdanalysis/blob/3e8ea8e0744f3163ddbf7b4f5f66874dd26796a7/package/MDAnalysis/coordinates/DLPoly.py#L114-L115

We probably need to change this to 

     forces = np.array(forces, dtype=np.float32, order='F') /100

and write according tests.

TRZ (IBIsCO, YASP)

Similarly, IBIsCO trajectories have force units of kJ/mol.nm , but TRZParser doesn't seem to convert units. Edit: From discussion below, force-reading should be dropped. So we should remove considerations of force in the TRZReader, including below:

https://github.com/MDAnalysis/mdanalysis/blob/3e8ea8e0744f3163ddbf7b4f5f66874dd26796a7/package/MDAnalysis/coordinates/TRZ.py#L179

https://github.com/MDAnalysis/mdanalysis/blob/3e8ea8e0744f3163ddbf7b4f5f66874dd26796a7/package/MDAnalysis/coordinates/TRZ.py#L267-L270

https://github.com/MDAnalysis/mdanalysis/blob/3e8ea8e0744f3163ddbf7b4f5f66874dd26796a7/package/MDAnalysis/coordinates/TRZ.py#L213-L223

orbeckst commented 4 years ago

@richardjgowers wrote these readers IIRC. I have no experience with them, but if they specify their formats then we better adhere to them!

We might have to issue a warning with these readers for at least one release cycle, saying that they now convert!

richardjgowers commented 4 years ago

So I hacked in writing ibisco force to the output, so I'm not sure that there's anyone else creating these files... And I can't remember what units these are in.

Not sure about the dl poly units, if they are specified we should probably convert? Tbh a better solution would be to keep track of units (maybe even in something like an xarray array) and do conversion properly (peek into astropy and see how they handle it?)

orbeckst commented 4 years ago

So I hacked in writing ibisco force to the output, so I'm not sure that there's anyone else creating these files... And I can't remember what units these are in.

How can we find out? Are there any docs?

Or if this is an MDA extension, then we need to document it here and whatever we say is the rule.

Or we drop it.

Not sure about the dl poly units, if they are specified we should probably convert?

Yes! — Are there docs for the file formats available somewhere? EDIT: sorry, DL_Poly docs were linked in the post...

Most file formats that we encounter come with defined units, so the units are part of the file format definition (the semantics). We should always be able to convert to our unit system and convert back without having to track units. (I imagine that astropy spans much larger length and time scales so it makes sense to be flexible with units.)

richardjgowers commented 4 years ago

Yeah we can drop TRZ and nobody will likely notice.

The DLP stuff looks confusing, section 1.3.7 says those are the internal units, which could well be different to the input/output units...

orbeckst commented 4 years ago

Yeah we can drop TRZ and nobody will likely notice.

Do you mean dropping support for forces in TRZ (because only your own custom code was writing forces)? The IBisCO docs for the II. coordinate file format only mention positions and velocities

This file contains the coordinates, the velocities and the connectivity table.

so by dropping force support we would simply be implementing the published format. I would be ok with this approach and dropping non-standard features.

orbeckst commented 4 years ago

DL_Poly looks complicated. Cited from the manual

1.3.7 Units

[...]

Note: In the DL POLY 4 OUTPUT file, the print out of pressure is in units of katms (kilo-atmospheres) at all times. The unit of energy is either DL POLY units specified above, or in other units specified by the user at run time (see Section 7.1.3). The default is the DL POLY unit. Externally, DL POLY 4 accepts information in its own specific formatting as described in Section 7.1. Irrespective of formatting rules, all values provided to define input entities are read in DL POLY units (except otherwise specified as in the case of energy units) or their composite mixture representing the corresponding entity physically, i.e. velocities’ components are in Ångstroms/picosecond.

It looks as if one can set units pretty freely in the input file (in particular the CONTROL and FIELD (force field) files).

However, it does not look as if the output files automatically carry the input file units. From the docs it looks as if the output files have well-defined, fixed units:

Here are some things I found

I think we should start by assuming that files are always in the documented DL_POLY units and add a note to the reader/writer stating this explicitly.

lilyminium commented 4 years ago

So should DLPoly.ConfigReader and DLPoly.HistoryReader change to convert units when reading in forces? They don't seem to do that at present.

orbeckst commented 4 years ago

Yes, I think they should.