Question about pressure calculation with PBCs if using ASE

[asedova]

I am not very familiar with ASE as a simulation engine, but this seems to be the easiest way to get started running MD inference with a NEQUIP model. However, I am not sure about the pressure/stress tensor calculation in ASE. There's a notorious issue with pressure calculation and PBCs. For a pairwise potential like in classical MD, the problem can be solved using the momentum flux approach and summing over pairwise forces (see for instance Haile Appendix B or the GROMACS manual), and for DFT, dE/dV is computed directly, by working it out analytically for each density functional family (like this: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.50.4327). For arbitrary many-body potentials without an explicit analytical derivative, there would be a problem (https://www.sciencedirect.com/science/article/abs/pii/S0009261406001552?via%3Dihub), which was solved by the LAMMPS folks (https://aip.scitation.org/doi/full/10.1063/1.3245303) and is implemented in LAMMPS (https://docs.lammps.org/compute_pressure.html). What about for ASE? Is the pressure/pressure tensor calculated with the LAMMPS method?

[Linux-cpp-lisp]

Hi @asedova ,

but this seems to be the easiest way to get started running MD inference with a NEQUIP model

While there's a little bit of one-time effort to get it compiled, we generally recommend that people use pair_nequip in LAMMPS to run MD with NequIP or Allegro models. ASE is a very useful code, but its not really a production MD engine.

There's a notorious issue with pressure calculation and PBCs.

Regarding the calculation of the stress tensor, we take the directly analytical dE/dV approach using autodifferentiation, see:
https://github.com/mir-group/nequip/blob/main/nequip/nn/_grad_output.py#L218

We provide this stress tensor directly to both ASE, via our calculator, and to LAMMPS, in pair_nequip--- neither computes the stress tensor on its own.

While we are fairly confident that this yields correct results, it is definitely a very subtle point, and if you find anything in reading the code or in your simulations that you think might be wrong please feel free to reach out here or by email (see profile).

Thanks!

[asedova]

Thanks so much for the detailed answer. That's a really great feature.

[Linux-cpp-lisp]

I should just note that the stress support in pair_nequip is still in beta / hasn't been merged yet, and is on the stress branch: https://github.com/mir-group/pair_nequip/tree/stress. Please let use know if you encounter any issues (or if it works well, thats good to know too 😁 )