Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Drop TRZ reading forces, and change DL_Poly units #2393

Open
lilyminium opened this issue Nov 6, 2019 · 11 comments · May be fixed by #4983
Open

Drop TRZ reading forces, and change DL_Poly units #2393

lilyminium opened this issue Nov 6, 2019 · 11 comments · May be fixed by #4983

Comments

@lilyminium
Copy link
Member

lilyminium commented Nov 6, 2019

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.

if has_forces:
forces = np.array(forces, dtype=np.float32, order='F')

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:

forces=self.has_force,

if self.has_force:
ts._forces[:, 0] = data['fx']
ts._forces[:, 1] = data['fy']
ts._forces[:, 2] = data['fz']

if not self.has_force:
frame_contents += [('pad7', '<i4')]
else:
frame_contents += [
('pad7', '<2i4'),
('fx', readarg),
('pad8', '<2i4'),
('fy', readarg),
('pad9', '<2i4'),
('fz', readarg),
('pad10', '<i4')]

@orbeckst
Copy link
Member

orbeckst commented Nov 7, 2019

@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!

@orbeckst orbeckst changed the title Are these the correct units? Are these the correct units for DL_POLY and TRZ (IBisCO and YASP)? Nov 7, 2019
@richardjgowers
Copy link
Member

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
Copy link
Member

orbeckst commented Nov 18, 2019

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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

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

  • 7.1.2.1 The CONFIG File Format (input coordinates): uses
    • Å for cartesian position coordinates (and cell box vectors)
    • Å/ps for velocities
    • Å * Dalton/ps2 for forces

    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.

  • REVCON (restart file) has the same format as the CONFIG file.
  • 7.2.1 The HISTORY File (trajectory output)
    • time: ps
    • cell vectors: Å
    • atomic mass: amu (Dalton)
    • atomic charge: e
    • position: ? (assume Å – would be nice to state it even it it is "obvious")
    • velocity: Å/ps
    • force: Å * Dalton / ps2

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
Copy link
Member Author

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
Copy link
Member

Yes, I think they should.

@lilyminium lilyminium changed the title Are these the correct units for DL_POLY and TRZ (IBisCO and YASP)? Drop TRZ writing forces and change DL_Poly units Feb 21, 2020
@lilyminium lilyminium changed the title Drop TRZ writing forces and change DL_Poly units Drop TRZ writing forces, fix TRZ reading force units, and change DL_Poly units Feb 21, 2020
@lilyminium lilyminium changed the title Drop TRZ writing forces, fix TRZ reading force units, and change DL_Poly units Fix TRZ reading force units, and change DL_Poly units Feb 21, 2020
@lilyminium lilyminium changed the title Fix TRZ reading force units, and change DL_Poly units Drop TRZ reading forces, and change DL_Poly units Feb 22, 2020
@lexi-x
Copy link
Contributor

lexi-x commented Mar 20, 2025

GSOC Contributor Inquiry

Hello!

I'm a prospective GSOC contributor hoping to work on this issue. I made the aforementioned edits to add unit conversions for force in the DL_Poly file with ConfigReader and HistoryReader and dropped the force-reading elements in TRZ in the highlighted locations. However, I have some questions before submitting a pull request.

  1. I removed the error checks for forces in the _read_trz_header method, but I was wondering if it is necessary to remove force from the TRZ trajectory header as well, as I am unsure if this will disrupt later processing.

  2. If I understood correctly, I need to add test cases to the respective files in the test suite. I was thinking of adding checks to see whether DL_Poly produces the right units and to see if the removal of force considerations in TRZ causes errors. I would love any recommendations on implementation or more tests that need to be written.

I am relatively new to open source, so any feedback is appreciated!

@orbeckst
Copy link
Member

@lexi-x can you put up your PR (and mention this issue number in it so that it gets linked)? Looking at code — even incomplete one — is a lot easier than talking about it in the abstract. Once the tests run in GitHub actions, you'll also get a coverage report that shows which parts of your code changes need to be tested.

I also should say that this is a fairly old issue and things have been changing. For instance, we are going to drop TRZ support #4311 so we are probably not going to be very concerned with TRZ. For DL_Poly we'll have to look. If you can work with DL_Poly yourself then that would be really useful, I think, because it's not a program that I have any personal experience with.

@lexi-x lexi-x linked a pull request Mar 20, 2025 that will close this issue
5 tasks
@lexi-x
Copy link
Contributor

lexi-x commented Mar 20, 2025

@orbeckst Thank you for your response! I just created a draft pull request with what I have so far, I kept the TRZ edits in the pull request just in case but noted for the future that TRZ is not a priority. Meanwhile, I will keep exploring what I can do with DL_Poly. Please let me know any suggestions you have based on the pull request, thank you!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants