VASP forces from OUTCAR

[tfrederiksen]

Nothing fancy, but this PR enables users to read atomic positions and forces from VASP (eg., in a relaxation run) from the OUTCAR file

• Closes #xxxx
• Tests added
• Documentation for functionality
• User visible changes (including notable bug fixes) are documented in CHANGELOG

[zerothi]

Thanks for this @tfrederiksen!

I think we need to conform to the idea that read_forces should only return forces, not the geometry as well.

Perhaps we should have a method read_md|steps|trajectory|path|? to return MD stuff (more than geometry|forces), this would fall in line with the read_scf method which also contains such information.

[tfrederiksen]

OK, for me it would be fine to rename the current function to read_trajectory. And we could also have a read_forces (relying on the current function). Is this an acceptable approach?

Note that the relevant block(s) in the OUTCAR contains both positions and forces, so it is basically for free to return both.

[zerothi]

OK, for me it would be fine to rename the current function to read_trajectory.

Yes, lets do something like this. What should we name it?
trajectory I now think is advocating only coordinates.

Ideas:

• steps
• time_steps

In any case I think we should return a PropertyDict with fields like:

steps = PropertyDict()
steps.geometries = GeometryCollection(geometries)
steps.forces = Collection(forces)
return steps

The advantage is that future additions can be done without disrupting code, i.e. some codes can return
energy information as well as forces etc.

And we could also have a read_forces (relying on the current function). Is this an acceptable approach?
Agreed, that would be ok, I think :)
It would be great if they could use the @sile_read_multiple decorator for automatic adding arguments.

Note that the relevant block(s) in the OUTCAR contains both positions and forces, so it is basically for free to return both.
Ok.

[tfrederiksen]

I prefer the wording trajectory here and think can imply more than just coordinates: in a general sense to system motion in some phase space. In an MD simulation the trajectory could be positions+velocities+cell and in a geometry relaxation to positions+forces+cell.

I like the suggestion of a propertyDict. However, the chemical species info is not available in OUTCAR, so maybe just something like this?

def read_trajectory(all=False):
    ...
    trajectory = PropertyDict()
    if all:
        trajectory.xyz = Collection(xyz)
        trajectory.force = Collection(force)
        trajectory.cell = Collection(cell)
    else:
        trajectory.xyz = xyz
        trajectory.force = force
        trajectory.cell = cell
    return trajectory

I note that sisl generally uses the singular form (read_geometry etc) so I would adopt this here too.

[zerothi]

Ok, lets do that then. Later when/if we can parse the data for the chemical species, we can add a geometry field.

[tfrederiksen]

It would of course be easy with the OUTCAR basepath to also read a POSCAR for the chemical species. But opening another file internally in this function seems a bit messy, or?

[zerothi]

It would of course be easy with the OUTCAR basepath to also read a POSCAR for the chemical species. But opening another file internally in this function seems a bit messy, or?

I have done this for e.g. the fdf file. If they are present, then it would be ok. Better to have consistent data. But then we shouldn't read the geometry, only the atomic information.

[tfrederiksen]

OK, I'll give it a try when I find time.

[tfrederiksen]

Maybe this has been discussed elsewhere, but couldn't the geometry class not optionally hold info about forces, strain, velocities, constraints, etc?

[zerothi]

Maybe this has been discussed elsewhere, but couldn't the geometry class not optionally hold info about forces, strain, velocities, constraints, etc?

This issue is tracked in #11 I have yet to find a stable solution that also works when manipulating the geometry. For instance, what would happen if one translates a subset of the coordinates, does it make sense to remove the forces (etc) or what should it do. But we should continue the discussion in #11

[codecov]

Codecov Report

Merging #573 (517dfe1) into main (adb10c5) will increase coverage by 0.00%.
The diff coverage is 95.73%.

❗ Current head 517dfe1 differs from pull request most recent head a659fcf. Consider uploading reports for the commit a659fcf to get more accurate results

@@           Coverage Diff           @@
##             main     #573   +/-   ##
=======================================
  Coverage   87.28%   87.29%           
=======================================
  Files         375      375           
  Lines       47564    47615   +51     
=======================================
+ Hits        41518    41565   +47     
- Misses       6046     6050    +4     

Impacted Files	Coverage Δ
src/sisl/io/__init__.py	100.00% <ø> (ø)
src/sisl/geometry.py	85.77% <25.00%> (-0.13%)	⬇️
src/sisl/io/_multiple.py	88.37% <88.37%> (+1.55%)	⬆️
src/sisl/io/sile.py	84.82% <96.15%> (ø)
src/sisl/io/vasp/stdout.py	95.95% <98.24%> (-0.30%)	⬇️
src/sisl/io/orca/stdout.py	97.61% <99.16%> (ø)
src/sisl/io/openmx/md.py	100.00% <100.00%> (ø)
src/sisl/io/orca/tests/test_stdout.py	100.00% <100.00%> (ø)
src/sisl/io/siesta/ani.py	100.00% <100.00%> (ø)
src/sisl/io/siesta/tests/test_ani.py	100.00% <100.00%> (ø)
... and 5 more

... and 16 files with indirect coverage changes

[zerothi]

Thanks @tfrederiksen! Some comments and improvements are needed.

[tfrederiksen]

@zerothi , the errors exposed above in the test runs seem unrelated to this branch.

[zerothi]

yeah, I need to check exactly why, I couldn't reproduce on my box, so I'll have another look when I have the time.

[zerothi]

minor details