Stresses, forces, and energy not read from xyz

[harveydevereux]

These quantities don't seem to be present when reading the trajectory with ase.io.read which is apparently a choice made there. Obviously these came be computed but as I have found on #189 these values differ from those computed during the simulation for normal i/o reasons I suppose.

If we specifically write them as new names, @alinelena tried e.g. mace-stress, then they are read into Atom.info I think and so can be accessed.

This came about because If we want to test that the correlator gets the same result as correlating directly from a trajectory file, we find the stresses computed are not the same, hence the correlations are unqeual at ~2-3 decimal place.

[ElliottKasoar]

Could you share an example of an xyz you can't read?

I have one from janus md which includes something like:

8
Lattice="5.64 0.0 0.0 0.0 5.64 0.0 0.0 0.0 5.64" Properties=species:S:1:pos:R:3:momenta:R:3:masses:R:1 spacegroup="P 1" unit_cell=conventional occupancy="_JSON {\"0\": {\"Na\": 1.0}, \"1\": {\"Cl\": 1.0}, \"2\": {\"Na\": 1.0}, \"3\": {\"Cl\": 1.0}, \"4\": {\"Na\": 1.0}, \"5\": {\"Cl\": 1.0}, \"6\": {\"Na\": 1.0}, \"7\": {\"Cl\": 1.0}}" time_fs=0.0 step=0 density=2.1636163899268936 energy=-27.035127799332745 free_energy=-27.035127799332745 stress="-0.004783275999053426 -2.0450471036772005e-19 -2.740042417777791e-19 -2.0450471036772005e-19 -0.004783275999053416 1.0800096722131524e-19 -2.740042417777791e-19 1.0800096722131524e-19 -0.004783275999053422" pbc="T T T"
Na       0.00000000       0.00000000       0.00000000       0.80328659       1.36571868      -0.65792736      22.98976928
Cl       2.82000000       0.00000000       0.00000000       0.48617771      -0.43936593      -0.02838953      35.45000000
Na       0.00000000       2.82000000       2.82000000      -0.29125007      -0.31928102       0.06499317      22.98976928
Cl       2.82000000       2.82000000       2.82000000       0.50929048      -0.10619493       0.60497278      35.45000000
Na       2.82000000       0.00000000       2.82000000       1.28422574       0.65451093      -0.07545326      22.98976928
Cl       0.00000000       0.00000000       2.82000000      -1.46429434       0.18885469       0.13842687      35.45000000
Na       2.82000000       2.82000000       0.00000000      -0.57732809      -0.66852154      -0.12921605      22.98976928
Cl       0.00000000       2.82000000       0.00000000      -0.75010802      -0.67572088       0.08259338      35.45000000

which I can read and get the results from:

atoms = read("test.xyz", index=":")
print(atoms[0].calc.results)

{'energy': -27.035127799332745, 'free_energy': -27.035127799332745, 'stress': array([-4.78327600e-03, -4.78327600e-03, -4.78327600e-03,  1.08000967e-19, -2.74004242e-19, -2.04504710e-19])}

[alinelena]

Could you share an example of an xyz you can't read?

I have one from janus md which includes something like:

8
Lattice="5.64 0.0 0.0 0.0 5.64 0.0 0.0 0.0 5.64" Properties=species:S:1:pos:R:3:momenta:R:3:masses:R:1 spacegroup="P 1" unit_cell=conventional occupancy="_JSON {\"0\": {\"Na\": 1.0}, \"1\": {\"Cl\": 1.0}, \"2\": {\"Na\": 1.0}, \"3\": {\"Cl\": 1.0}, \"4\": {\"Na\": 1.0}, \"5\": {\"Cl\": 1.0}, \"6\": {\"Na\": 1.0}, \"7\": {\"Cl\": 1.0}}" time_fs=0.0 step=0 density=2.1636163899268936 energy=-27.035127799332745 free_energy=-27.035127799332745 stress="-0.004783275999053426 -2.0450471036772005e-19 -2.740042417777791e-19 -2.0450471036772005e-19 -0.004783275999053416 1.0800096722131524e-19 -2.740042417777791e-19 1.0800096722131524e-19 -0.004783275999053422" pbc="T T T"
Na       0.00000000       0.00000000       0.00000000       0.80328659       1.36571868      -0.65792736      22.98976928
Cl       2.82000000       0.00000000       0.00000000       0.48617771      -0.43936593      -0.02838953      35.45000000
Na       0.00000000       2.82000000       2.82000000      -0.29125007      -0.31928102       0.06499317      22.98976928
Cl       2.82000000       2.82000000       2.82000000       0.50929048      -0.10619493       0.60497278      35.45000000
Na       2.82000000       0.00000000       2.82000000       1.28422574       0.65451093      -0.07545326      22.98976928
Cl       0.00000000       0.00000000       2.82000000      -1.46429434       0.18885469       0.13842687      35.45000000
Na       2.82000000       2.82000000       0.00000000      -0.57732809      -0.66852154      -0.12921605      22.98976928
Cl       0.00000000       2.82000000       0.00000000      -0.75010802      -0.67572088       0.08259338      35.45000000

which I can read and get the results from:

atoms = read("test.xyz", index=":")
print(atoms[0].calc.results)

{'energy': -27.035127799332745, 'free_energy': -27.035127799332745, 'stress': array([-4.78327600e-03, -4.78327600e-03, -4.78327600e-03,  1.08000967e-19, -2.74004242e-19, -2.04504710e-19])}

this is practically the ase issue, that I have fixed for single point... new ase silently will not read forces/stress/energy we will need to follow the same as we did in singlepoint and abuse info and arrays

[harveydevereux]

@alinelena I think @ElliottKasoar is correct thought, we just simply did not know about Atoms.calc.results it is there

tests/test_correlator.py::test_md_correlations {'energy': -27.035127799332745, 'free_energy': -27.035127799332745, 'stress': array([-4.78327600e-03, -4.78327600e-03, -4.78327600e-03, -1.90037373e-19,
       -1.08369307e-20, -1.67012180e-19])}

so unless you think otherwise this issue can be closed

[ElliottKasoar]

this is practically the ase issue, that I have fixed for single point... new ase silently will not read forces/stress/energy we will need to follow the same as we did in singlepoint and abuse info and arrays

ASE is reading them in this case though, it just puts them into atoms.calc.results rather than info.

For single point it wasn't too bad to fix, since we call get_potential_energy and write ourselves, but given that the ASE MD functions do both in this case, I'm not sure there's a way to move/copy the results to info for every image without quite a large overhead, so is it definitely worth it?

Do you know if there are plans to tidy up how ASE handles all of this?

Sorry @harveydevereux, your comment hadn't loaded when I wrote this. I think we did discuss this for single point too, and for better generality (e.g. having multiple energies) using info was preferable, and safer from further ASE changes... I'm just not sure whether it's worth it in this case. I may well be forgetting stronger reasons, though.

[ElliottKasoar]

For single point it wasn't too bad to fix, since we call get_potential_energy and write ourselves, but given that the ASE MD functions do both in this case, I'm not sure there's a way to move/copy the results to info for every image without quite a large overhead, so is it definitely worth it?

This is wrong actually, we can modify the _write_traj function e.g.:

for prop in self.dyn.atoms.calc.results:
    self.dyn.atoms.info[f"mace_{prop}"] = self.dyn.atoms.calc.results[prop]
self.dyn.atoms.calc.results = {}

So it wouldn't be too bad, although we'd have to decide whether all of the results (e.g. "node_energy") or just a subset should be saved.

I maintain that the information discussed is there, and readable through ase.io.read, so I think it's a separate issue whether we want to do this to make things neater from a labelling perspective.

[ElliottKasoar]

Closing here in favour of #206, as this isn't a bug as described, but does remain an inconsistency between calculations that we should clean up, as discussed with @alinelena.