Different potential energy calculation results

Dear LAMMPS users,

I am referring the function single () defined in pair.cpp to calculate the pairwise force coefficient fpair and the potential energy.
eng=pair->single (i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
The current potential I am using is EAM potential.

The weird observation I have is that my summation of the potential energy (eng value of each pair) from all the pairs does not match with the output from log file. The value I get is -21936114 and the value from log file is -63005273, the difference is close to 3 times.

The calculation of pressure based on fpair seems to be consistent with the value from log file, but the difference in potential energy is hard to understand.

I will appreciate it so much if anyone can give me any idea about this issue.
Thanks!

Dear LAMMPS users,

I am referring the function *single ()* defined in pair.cpp to calculate
the pairwise force coefficient *fpair* and the *potential energy*.
eng=pair->single (i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
The current potential I am using is EAM potential.

The weird observation I have is that my summation of the potential energy
(eng value of each pair) from all the pairs does not match with the output
from log file. The value I get is -21936114 and the value from log file is
-63005273, the difference is close to 3 times.

The calculation of pressure based on fpair seems to be consistent with the
value from log file, but the difference in potential energy is hard to
understand.

​your e-mail is missing two important pieces of information.

1) which LAMMPS version you are using exactly. what you are reporting
sounds a lot like a bug that was fixed not so long ago.
2) an example input deck that allows to quickly reproduce your issue (in
case this is a real bug and still exists with the latest patch release).

axel.​

Dear Dr. Kohlmeyer,

Sorry for the missing information.
I am using the latest stable version of lammps, 11Aug2017.
I tried another case with Morse potential, which generates consistent potential value to the log file.
While for EAM potential, it does not.

I am not sure how to quickly reproduce the simulation, because I do the modification in the cpp code.
This is how I did in the code.
I put this code segments within the end_of_step part in fix_ttm_mod.cpp and output the summation of all the potential energy as potential_all[0],
which is used for comparison with the log file.

int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, eng, r, dr, dexp, factor_coul, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double vir[3];

double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;

inum = force->pair->list->inum;
ilist = force->pair->list->ilist;
numneigh = force->pair->list->numneigh;
firstneigh = force->pair->list->firstneigh;

potential_single[0] = 0.0;
potential_all[0] = 0.0;

for (ii = 0; ii < inum; ii++)
{
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype])
{

eng = pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
potential_single[0] += eng;

}
}
}
MPI_Allreduce(&potential_single[0], &potential_all[0], 1, MPI_DOUBLE, MPI_SUM, world);

I am not sure why different results can be generated for different potential, Morse and EAM,
and is there any difference to refer the same function single() by using different potentials.

Thanks!

Dear Dr. Kohlmeyer,

Sorry for the missing information.
I am using the latest stable version of lammps, 11Aug2017.
I tried another case with *Morse potential*, which generates consistent
potential value to the log file.
While for EAM potential, it does not.

​[...]​

I am not sure why different results can be generated for different
potential, Morse and EAM,
and is there any difference to refer the same function single() by using
different potentials.

​there is a very subtle difference between morse and eam. the former is a
true pair-wise additive potential, the latter is not.
the PairEAM::single() function is dependent on position dependent
properties in data structures, that are filled with suitable data during
the original force computation.​

​do you get the same "wrong" number, when you are using compute group/group
(which should be doing something similar to what you are doing with your
modification?​

​axel.​

Dear Dr. Kohlmeyer,

I am not sure about the results from compute group/group.
If that is the case, does it mean that the single() function cannot be used directly for the calculation the atom potential in EAM potential.
I feel by referring this function, the coeff of fpair is calculated correctly.
Is there other way to refer the value of single pair value?

Thanks!

Dear Dr. Kohlmeyer,

In the code of pair_eam.cpp, the overall potential energy seems to be also summed from pairwise potential energy as

recip = 1.0/r;
phi = z2recip;
phip = z2p
recip - phirecip;
psip = fp[i]rhojp + fp[j]rhoip + phip;
fpair = -psip
recip;
f[i][0] += delx
fpair;
f[i][1] += dely
fpair;
f[i][2] += delzfpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx
fpair;
f[j][1] -= delyfpair;
f[j][2] -= delz
fpair;
}
if (eflag) evdwl = phi;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);

the value of phi should be the same output from single () function.

This part seems to be similar with pairwise potential like Morse

if (eflag) evdwl = phi;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);

So the I still have the curiosity about the difference in calculation of potential.
Thanks!

You are only looking at the pair term from the EAM potential, but EAM in fact includes another 3-body term called the “embedding” term.

Just several lines above the snippet of code you pasted, the code first calculates the embedding density (rho) for each of the atoms (I-J loop), then calculates the embedding energy (phi) using an I-loop. Then the variable phi is reused to calculate the pair potential in the next I-J loop (the snippet you pasted). By looking only at the second I-J loop, you are missing the embedding term.

The single() function only calculates the pair term, not including the 3-body embedding term.

Ray

Dear Dr. Shan,

Thanks for your explanation.
I may want to use the atom potential value in external fix.cpp,
if this is the case, is there any way to include this embedding term into the potential energy by referring single() function or other method?

Thanks!

Dear Dr. Shan,

Thanks for your explanation.
I may want to use the atom potential value in external fix.cpp,
if this is the case, is there any way to include this embedding term into
the potential energy by referring single() function or other method?

​you cannot compute the embedding energy directly from Pair::single(),
since it is not a pairwise term, but a per-atom energy (with contributions
from its neighbors).
that embedding energy term isn't stored, so it cannot be accessed. you
would have to modify pair_eam.cpp to make it accessible.

if you just need the total (pair style) potential energy term​, you can
access it via the thermo keyword epair. if you need the per-atom
decomposition, you can define compute pe/atom and access it this way.

​in short, you need to provide more detail, what exactly it is that you
need to use as input in your fix.

axel.​

Dear Dr. Kohlmeyer,

Thanks a lot for your explanation.
I understand now the difference between these different potential from their nature and definition.
I will try to solve this problem by access to the variable pair->eatom, together with the potential from pair->single().

Thanks!

Dear Dr. Kohlmeyer,

Thanks a lot for your explanation.
I understand now the difference between these different potential from
their nature and definition.
I will try to solve this problem by access to the variable pair->eatom,
together with the potential from pair->single().

​that will count the pairwise contributions twice.

axel.

Dear Dr. Kohlmeyer,

Thanks a lot for your remind.
I also notice this so I just use eatom to count the total potential energy.
Still there is little difference like -62655360.55 vs -63005273, this difference is much smaller,
is this due to calculation issue or anything else?

Thanks!

Dear Dr. Kohlmeyer,

Thanks a lot for your remind.
I also notice this so I just use eatom to count the total potential
energy.
Still there is little difference like -62655360.55 vs -63005273, this
difference is much smaller,
is this due to calculation issue or anything else?

​i don't know. i don't have a crystal ball.

axel.​

Dear Dr. Kohlmeyer,

The value of -62655360.55 is from the summation of eatom[i] for i<nlocal,
if the nlocal is changed to nall=nlocal+nghost, the potential value will be exactly same as 63005273,
so what is the reason to include the potential energy from ghost atoms on each processor but not only the local atoms.

Thanks!