Inconsistent pressure in NVT runs between LAMMPS and ASE - which one to trust?

[casteln]

Hello,

When running NVT trajectories of the same system and NequIP model (learned with stress) with ASE and LAMMPS, I get different average results for the pressure.
I find a systematic shift of roughly 200 MPa, with the pressure outputted by LAMMPS always higher.

I've created a minimal reproducible example available here, which details the thermostat and versions used, as well as how the pressure is read from the output files.
It also contains a figure that shows that this shift was observed across models (i.e. that it's not an issue with the units).

Would you know where this is coming from? And which value should I trust?

[Linux-cpp-lisp]

Hi @casteln ,

Thanks very much for the detailed report!!

Which LAMMPS branch? One of them lacks stress support and will only give the kinetic term. Sorry about this, I've been meaning to clean up that branch situation for a while now and have not yet have the time, I expect to very soon.

[casteln]

Any idea on why the pressure output may still be different in the stress branch of pair_nequip compared to ASE?

[ssrokyz]

Does the stress function of the ASE NequIP calculator also include the kinetic term?
LAMMPS calculates the stress as the sum of kinetic and viral terms [check this out].
The LAMMPS kinetic term is naturally included by LAMMPS itself, while I don't know in the ASE case.

[casteln]

Thank you @ssrokyz for your insight!
From this answer by @Linux-cpp-lisp, it would appear that the stress computed by ASE is only the gradient of the potential energy, and thus do not include the kinetic term.

My understanding is therefore that the difference may come from this kinetic term.
Do you know if it's easy to compute this kinetic term with ASE?

[Linux-cpp-lisp]

Hi all, sorry for my late response here and thanks @ssrokyz for chiming in!!

Yes, I would suspect that this is the case. You can see how in the pair_nequip unit tests we ask LAMMPS for the stress tensor (not the virial, but still without the temperature / kinetic term) which is what we compare to the output from ASE:
https://github.com/mir-group/pair_nequip/blob/stress/tests/test_python_repro.py#L148-L152
and
https://github.com/mir-group/pair_nequip/blob/stress/tests/test_python_repro.py#L317-L344

[casteln]

Thank you for answer!
I've considered LAMMPS stress output after disabling the kinetic term, and found that it gave the same average pressure value than ASE.

[Linux-cpp-lisp]

Great!