In the documentation I saw the Virial term is computed following Thompson, Plimpton, Mattson JCP 2009. I imagine Lammps uses equation 25 in that paper (the “atom” form):
W=sum_i sum_n r_in . F’_in
Where F’_in is “the partial force on the atom located at ri+Hn due to all the groups associated with the local cell”. My concern is about the “groups associated with the local cell” part.
If a local atom has an image that belong to a group associated with the local cell, then that local atom will be part of the image of that group, which is NOT associated with the local cel (e.g. atom 2 in FIG 1 of the paper). Therefore, virial terms related to that local atom and neighboring atoms of that group image should not be take into account for the virial calculation.
Looking at PairLJCubic::compute routine and the corresponding PairLJCubic::compute I can not figure out how the virial terms arising from the images of the group associated with the cell are filter out.
I would appreciate your help to understand how this work.
Thank you a lot,
Alexis
The last line of the compute() method is this:
if (vflag_fdotr) virial_fdotr_compute();
which is the case for nearly all the pairwise potentials in LAMMPS.
It compute the virial as sum of (F dot R) for
all owned and ghost atoms (which includes atoms
across periodic boundaries). The ghost atoms
have coords which are unwrapped by the periodicity.
That’s all there is to it.
Steve
Thank you. So, if lammps do that, I might have a problem to understand the paper. I will try to explain me better.
Consider a 1 D system with a diatomic molecule centered just in the limit of the box with lenght L. So the system including ghost atoms looks like:
0 L
-------2—|---1------------------2—|---1--------
_____/ group
Where I also indicated the “group associated with the local cell” of this system. Using equation 25 of the paper virial will be:
W=x1*f1 + (x2-L)*f2
I did not include the term x2*f2 and (x1+L)*f1 since atoms at x2 and x1+L do not interact with any group associated with the local cell.
If lammps add all the local and ghost atoms terms, x2*f2 and (x1+L)*f1 will be there. Should not these terms be avoided?
PS: Using eq 22 of Bekker et al. Molecular Simulation 14 (1995) 3 I get W’=(x1+L)f1 + x2f2 for this system, which is equivalent to W (also half of the terms are out). I also get a similar result using equation B.29 of Haile book “Molecular dynamic simulations”.
Thank you. So, if lammps do that, I might have a problem to understand the paper. I will try to explain me better.
Consider a 1 D system with a diatomic molecule centered just in the limit of the box with lenght L. So the system including ghost atoms looks like:
0 L
-------2—|---1------------------2—|---1--------
_____/ group
Where I also indicated the “group associated with the local cell” of this system. Using equation 25 of the paper virial will be:
W=x1*f1 + (x2-L)*f2
I did not include the term x2*f2 and (x1+L)*f1 since atoms at x2 and x1+L do not interact with any group associated with the local cell.
If lammps add all the local and ghost atoms terms, x2*f2 and (x1+L)*f1 will be there. Should not these terms be avoided?
PS: Using eq 22 of Bekker et al. Molecular Simulation 14 (1995) 3 I get W’=(x1+L)f1 + x2f2 for this system, which is equivalent to W (also half of the terms are out). I also get a similar result using equation B.29 of Haile book “Molecular dynamic simulations”.
Thank you. So, if lammps do that, I might have a problem to understand the
paper. I will try to explain me better.
Consider a 1 D system with a diatomic molecule centered just in the limit
of the box with lenght L. So the system including ghost atoms looks like:
0 L
-------2---|---1------------------2---|---1--------
\_____/ group
Where I also indicated the "group associated with the local cell" of this
system. Using equation 25 of the paper virial will be:
W=x1*f1 + (x2-L)*f2
I did not include the term x2*f2 and (x1+L)*f1 since atoms at x2 and x1+L
do not interact with any group associated with the local cell.
If lammps add all the local and ghost atoms terms, x2*f2 and (x1+L)*f1
will be there. Should not these terms be avoided?
this term will not be there with the typical setup (i.e. newton's 3rd law
enabled for interactions crossing subdomain boundaries). in this case, each
subdomain crossing interaction is listed only once in the neighbor lists.
when newton's 3rd law is disabled, no forces are accumulated on ghost atoms
and the reverse communication will not be needed.
axel.