Hi all,

I am going to calculate the heat flux vector based on contribution from atoms in the computational domain (something like compute heat/flux) in Tersoff. For example, for two body term, I multiplied 0.5*fpair by velocity of each atom(Please see the highlighted lines). However, the results that I get from LAMMPS is not similar to the code I wrote in Fortran. I finally realized that atom i velocity components (V[i][0], V[i][1],V[i][2]) and some of atoms j velocity components (V[j][0], V[j][1],V[j][2]) is exactly similar to the results of my code but the velocity for some of j atoms (maybe 20%) is zero in lammps but non-zero in my code. I have no idea why the velocity of some j atoms (V[j][0], V[j][1], V[j][2]) is zero in lammps.

Please note that for calculation of convective term of heat flux vector ( sum of eiVi over all atoms). The results of my code is exactly similar to lammps results meaning that velocity calculation for both code is same. However, when it comes to calculation of second term of heat flux (sum of fij*rij*vi). I have this problem.

Thank you very much

Hossein

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];

itag = tag[i];

itype = map[type[i]];

xtmp = x[i][0];

ytmp = x[i][1];

ztmp = x[i][2];

// two-body interactions, skip half of them

jlist = firstneigh[i];

jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {

j = jlist[jj];

j &= NEIGHMASK;

jtag = tag[j];

if (itag > jtag) {

if ((itag+jtag) 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) 2 == 1) continue;

} else {

if (x[j][2] < x[i][2]) continue;

if (x[j][2] == ztmp && x[j][1] < ytmp) continue;

if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;

}

jtype = map[type[j]];

delx = xtmp - x[j][0];

dely = ytmp - x[j][1];

delz = ztmp - x[j][2];

rsq = delx*delx + dely*dely + delz*delz;

iparam_ij = elem2param[itype][jtype][jtype];

if (rsq > params[iparam_ij].cutsq) continue;

repulsive(¶ms[iparam_ij],rsq,fpair,eflag,evdwl);

f[i][0] += delx*fpair;
f[i][1] += dely*fpair;

f[i][2] += delz

*fpair;*

f[j][0] -= delxfpair;

f[j][0] -= delx

f[j][1] -= dely

*fpair;*

f[j][2] -= delzfpair;

f[j][2] -= delz

if (evflag) ev_tally(i,j,nlocal,newton_pair,

evdwl,0.0,fpair,delx,dely,delz);

double **v = atom->v;

energyBit = evdwl;

atomE[i] += 0.5*energyBit;
qdpi = 0.5*fpair*(delx

*v[i][0]+dely*v[i][1]+delz

*v[i][2]);*

q[i][0] += qdpidelx;

q[i][0] += qdpi

q[i][1] += qdpi

*dely;*

q[i][2] += qdpidelz;

q[i][2] += qdpi

atomE[j] += 0.5

*energyBit;*

qdpj = -0.5fpair*(delx

qdpj = -0.5

*v[j][0]+dely*v[j][1]+delz

*v[j][2]);*

q[j][0] -= qdpjdelx;

q[j][0] -= qdpj

q[j][1] -= qdpj

*dely;*

q[j][2] -= qdpjdelz;

q[j][2] -= qdpj