I have implemented a many-body potential based on neural networks and the 4D bispectrum (i.e., Spectral Neighborhood Analysis, SNA in LAMMPS). The potential currently produces correct energies and forces, but the pressure calculation is significantly off of my reference calculations. The pressure calculated within LAMMPS also does not correspond to the pressure calculated by numerical differentiation of the energy with respect to the cell volume. If I understand the article in JCP 131 (154107 (2009), https://aip.scitation.org/doi/10.1063/1.3245303) correctly, then the virial calculation in LAMMPS should be valid for arbitrary many-body potentials and “all I need” should be correct forces and positions plus the correct (local) summation cells. As with other many-body potentials in LAMMPS, I thought a call to Pair::virial_fdotr_compute at the end of Pair::compute should thus be sufficient.
Is there anything that I am missing from a technical point of view and that might interfere with the validity of the virial calculation? What are the pre-requisites for making Pair::virial_fdotr_compute work correctly (technically)? I could imagine that there might be something from a neighborlist-setup or something in relation to the fact that newton_pair/force->newton needs to be on in my case because of the SNA(D) calculation.
I would appreciate your help and will of course provide more detailed explanations and the code, if required. Some details are listed at the end of the mail.
Best regards,
Frank
More Details:
Only the inference is done within LAMMPS using a small C++ neural network class. The training is done based on data extracted from LAMMPS (SNA and SNAD components) using a python-based framework called pytorch. The potential (i.e., Pair::compute) is then just implemented as a loop over nlocal atoms, where in which in each iteration the neural network energy/forces are evaluated according to the SNA(D) components. The resulting energies and forces correspond very well to the reference calculations, both from the evaluations in the training code and the original reference calculations.
Everything you say seems to be correct. However you don’t say anything about ghost atoms, so that may be where the problem is. The fdotr method relies on each force on a neighbor j due to the energy of atom i being applied to one and only one periodic image of j. Eventually, all the contributions on different ghost atoms are accumulated to the corresponding local atom via reverse communication, but the fdotr calculation is performed before that step. I think if you look at how force and virial are handled in pair_style snap, this will be clearer to you.
thank you for your time and response. Now it took me a while to get
back to it and the response grew accordingly.
I didn't do anything with the ghosts at that time and I think I now
better see how that is problematic. Re-reading your 2009-paper on this
topic, I gather that, in a simplified manner, you calculate the virial
in separate groups that at that moment don't see the rest of the
system (i.e., before reverse communication). Those groups are
effectively isolated from the rest of the system by only accounting
for energy contributions from (local and ghost) atoms within the local
(minimum-image) cell.
Is that more or less correct?
I am trying to figure out if that is applicable to my artificial
neural network (ANN) potentials as well. The potential energy in an
ANN with non-linear activation functions cannot be separated into
group contributions that sum up to the correct total. Hence,
Pair:virial_fdotr_compute would not be applicable at all. Correct?
As a simple test, I used the snav/atom compute and multiplied its
per-atom values with the ANN energy-derivative with respect to the
descriptor. With that, I get close to finite difference calculations
of the pressure. However, I think the former is not completely right
because of the energy separation issue.
The energy of the atomic ANNs was deduced from all of the neighbors of
an atom. Which should be okay for local atoms, but not for ghost
atoms.
Any insight would again be appreciated. A couple more separate
questions arose while thinking about this:
Looking at your SNAP implementation, how is double counting of force
contributions avoided? With full neighbor lists and newton_pair on,
forces on ghost atoms from local atoms also appear as forces on local
atoms from ghost atoms (for a different local cell), don't they?
Btw, in the neighbor lists are all neighbors of the ghost atoms also included?
The potential energy in an
ANN with non-linear activation functions cannot be separated into
group contributions that sum up to the correct total.
This depends on how your define your model. Most ANN models that I have seen do the same thing SNAP does, which is write the potential energy as a sum over single-atom contributions, each of which can now be handled as an isolated group. The processor that owns the central atom calculates the group energy and associated forces. That is how we avoid double counting forces in SNAP.
Indeed, most ANN models represent the energy as a sum of atomic contributions.
However, the individual atomic energies depend on all atoms (typically
within a certain cutoff) in a non-additive fashion.
After calculating the atomic energies from the whole system, then they
are separable.
I.e., if the energy of an atom is calculated in parts with an
incomplete set of neighbors, then the resulting addition will be
wrong.
This, however, appears not to be a problem for the virial calculation
when using the virial_fdotr_compute method.
In the force calculation, in the loop over j neighbors of i, one only
needs the derivative of the i'th atomic neural network for the
respective force contributions to i and j.
And for this network, all neighbors of i within the respective cutoff
are accounted for.
I implemented the neural network potential following your SNAP
implementation and forces and pressure both correspond well to finite
difference calculations, respectively.