Can not reproduce virial calculation

[FromColleWithLove]

Hi,

I'm trying to use ACEpotentials on a dataset I created, and I cannot find the definition of virial used by the code.
My dataset uses the LAMMPS' definition, as in https://docs.lammps.org/compute_pressure.html , and it sums f_i * r_j for over the atoms.
When I try to do the same with the Si dataset provided in the examples, I obtain different values.
What is the stress calculation performed for training? Maybe I missed e.g. periodic images contributions and volume tries?

[wcwitt]

Hi, can you send a code snippet showing your approach for the Si dataset? I would expect they are the same up to a sign (although always be careful with units in LAMMPS).

[cortner]

and it sums f_i * r_j for over the atomsJ

I'm confused by this expression. If this is indeed the sum of forces * positions, then this is not the virial. E.g. in a one-atom unit cell this would always be zero. So you probably mean something else?

[FromColleWithLove]

Yes. I agree on that.
This expression is taken from
https://docs.lammps.org/compute_pressure.html
I was wondering how you computed stresses because I have an ab initio trajectory for other uses (FLARE) and it adds a stress calculator. I wished to make the same trajectory compatible, therefore was not sure whether the LAMMPS default calculator was correct

[jameskermode]

I think the issue is that that expression is valid only for a pair potential. In general, the virial components are computed as the derivative of the potential energy with respect to the cell.

[cortner]

Ok. If that’s your question then I suggest you just look into the source code. For site potentials it’s a simple formula.

ACEpotentials.jl/src/models/calculators.jl

Line 105 in 6b95848

_site_virial(dV::AbstractVector{SVector{3, T1}},

[FromColleWithLove]

@wcwitt

import numpy as np
from ase.io import read

atoms=read("Si_dataset.xyz",index=200,format="extxyz")

def virial(forces,positions,i,j):
    virial=0
    for f,x in zip(forces[:,i],positions[:,j]):
        virial+=f*x
    return virial

pos=atoms.get_positions()
vel=atoms.get_velocities()
forces=atoms.arrays['dft_force']
data_vir= atoms.info['dft_virial']
print(data_vir)
data_vir=data_vir[0]
print("As in dict  , computed")
my_vir=virial(forces,pos,0,0)/3
print(f"{data_vir} eV(?) , {my_vir} eV")
for i in range(3):
    for j in range(3):
        print(i,j,virial(forces,pos,i,j))

Output:

[ 7.10295495e+01 -6.36740730e-01 -1.38160460e-01 -6.36740730e-01
  6.87016022e+01 -3.68170500e-02 -1.38160460e-01 -3.68170500e-02
  6.35478722e+01]
As in dict  , computed
71.02954952 eV(?) , 0.6831223959085334 eV
0 0 2.0493671877256
0 1 5.9200526851071995
0 2 -1.3249875238165005
1 0 3.0732016549816015
1 1 2.208839692488698
1 2 2.9895363949582996
2 0 -0.07916525015140063
2 1 -0.41241244087150153
2 2 -4.312182497595102

This assumes that units are eV and Angstroms, I see virials are in eV so I assumed they are not normalized by volume.
The fact that my results are not symmetrical confusea me. How are the virials in the dataset computed? Am I missing something crucial?
And importantly, the first configurations of the dataset (sh) with null forces should have a null virial, shouldn't them? Unless volume tries are accounted for.

[cortner]

See my comment below. This is not the virial. Not even up to bugs or signs. I’m not sure what a good source is to recommend to learn about this. Maybe Chuck has a recommendation?

[tjjarvinen]

Here is alternative resource for for virial definition https://en.wikipedia.org/wiki/Virial_stress that gives a different view.

The main difference is that you are using total force to calculate virial - You cannot do this!

The way to calculate virial is from the individual force components before you sum up the forces. This is also how it is done in our source code, which you can see by looking the code location Christoph linked.

This means that you cannot recreate the virials from the training set data alone.