Fitting GAP virials from lammps simulation

[MarcBrueck]

I am trying to fit a GAP potential to simulation data, that I got from a classical MD lammps simulation.

I use the compute_pressure command (https://lammps.sandia.gov/doc/compute_pressure.html), that calculates the pressure tensor:

                           P_ij = sum(m * v_i * v_j) / V + sum(r_i * F_j) / V

Right now I am using only the second part of this sum the so called virial term. Is this assumption correct?

Then I am multiplying this pressure with the volume of my system, convert it to eV and write it into my extended xyz file:
virial_ij = - P_ij * Volume

Sincerely
Marco Brueckner

[gabor1]

This looks correct. Do you have any concerns?

…

-- Gábor

On 28 Jan 2021, at 09:39, Marco Brueckner ***@***.***> wrote:  I am trying to fit a GAP potential to simulation data, that I got from a classical MD lammps simulation. I use the compute_pressure command (https://lammps.sandia.gov/doc/compute_pressure.html), that calculates the pressure tensor: P_ij = sum(m * v_i * v_j) / V + sum(r_i * F_j) / V Right now I am using only the second part of this sum the so called virial term. Is this assumption correct? Then I am multiplying this pressure with the volume of my system, convert it to eV and write it into my extended xyz file: virial_ij = - P_ij * Volume Sincerely Marco Brueckner — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[MarcBrueck]

The energies and forces seem to be predicted correctly by my trained potential and also the off diagonal components of the virial tensor, but the main diagonal component are shifted by a constant. This is important to me since I am trying to get npt simulations to run.
How does GAP calculate the stress of the system?

I get the predicted virial by:

configuration.calc = potential
predicted_virial.append(configuration.get_stress())

where configuration is an ase atoms object and potential is the gap.xml potential file.

Energy_validation.pdf
Force_component_validation.pdf
Virial_main-axis-components_validation.pdf
Virial_off-diagonal-components_validation.pdf

Sincerely Marco Brueckner

[bernstei]

On Feb 9, 2021, at 11:17 AM, Marco Brueckner <[email protected]<mailto:[email protected]>> wrote: The energies and forces seem to be predicted correctly by my trained potential and also the off diagonal components of the virial tensor, but the main diagonal component are shifted by a constant. This is important to me since I am trying to get npt simulations to run. How does GAP calculate the stress of the system? Correctly :) Same as any other gradient - analytical derivatives of the covariance kernels which lead to analytical derivatives of the energy w.r.t. cell deformation. Diagonal stress components off by a constant is the same as getting the right bulk modulus but being off by a constant for the equilibrium volume. Have you compared E(V) curves for DFT and GAP? I get the predicted virial by: configuration.calc = potential predicted_virial.append(configuration.get_stress()) Not that this is not quite right, if you take the variable name seriously, because of factors of -1 and cell volume. But if you just call it "predicted_stress" it'll be consistent with what you're putting in there, and the right thing to compare to an ASE-computed DFT stress, e.g.

[MarcBrueck]

I forgot to mention that I multiply the predicted stress from GAP with the volume. If I wouldn't then the predicted virial and simulated virial would differ by 3 orders of magnitude. The off diagonal components are predicted correctly though suggesting to me that I did the process correctly. The constant offset for the main diagonal components suggests to me that there is a term missing somehow between the lammps compute pressure command and the calculation by GAP.

[bernstei]

On Feb 10, 2021, at 2:48 AM, Marco Brueckner <[email protected]<mailto:[email protected]>> wrote: I forgot to mention that I multiply the predicted stress from GAP with the volume. If you pass that information in the form of the "virial", that's the correct magnitude, but off by a negative sign from the _ASE_ definition of stress, I believe. I'm not sure whether LAMMPS and ASE agree on the sign convention. Note that recent versions of libAtoms/QUIP and GAP (on github) can also read from fields called "stress", and that has the same magnitude (virial/volume), units (eV/A^3), and sign convention as ASE's stress calculation. If I wouldn't then the predicted virial and simulated virial would differ by 3 orders of magnitude. The off diagonal components are predicted correctly though suggesting to me that I did the process correctly. The constant offset for the main diagonal components suggests to me that there is a term missing somehow between the lammps compute pressure command and the calculation by GAP The only difference I can think of is that lammps may be including the kinetic (ideal gas) contribution, and GAP would only report the potential energy gradient contribution. Have you looked at E(V) curves? Noam

[gabor1]

In the original post it is stated that the kinetic term is not used when passing the virial into the XYZ for training. how about in your tests above, are you also only using the second term in your original post?