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,
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 *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;
potential_all = 0.0;
for (ii = 0; ii < inum; ii++)
i = ilist[ii];
xtmp = x[i];
ytmp = x[i];
ztmp = x[i];
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];
dely = ytmp - x[j];
delz = ztmp - x[j];
rsq = delx*delx + dely*dely + 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 += eng;
MPI_Allreduce(&potential_single, &potential_all, 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.