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 = 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[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!