Velocity in the code

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 fijrijvi). 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 = delxdelx + delydely + delz*delz;

iparam_ij = elem2param[itype][jtype][jtype];
if (rsq > params[iparam_ij].cutsq) continue;

repulsive(&params[iparam_ij],rsq,fpair,eflag,evdwl);

f[i][0] += delxfpair;
f[i][1] += dely
fpair;
f[i][2] += delzfpair;
f[j][0] -= delx
fpair;
f[j][1] -= delyfpair;
f[j][2] -= delz
fpair;

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.5energyBit;
qdpi = 0.5
fpair*(delxv[i][0]+delyv[i][1]+delzv[i][2]);
q[i][0] += qdpi
delx;
q[i][1] += qdpidely;
q[i][2] += qdpi
delz;
atomE[j] += 0.5energyBit;
qdpj = -0.5
fpair*(delxv[j][0]+delyv[j][1]+delzv[j][2]);
q[j][0] -= qdpj
delx;
q[j][1] -= qdpjdely;
q[j][2] -= qdpj
delz;

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.

you do this calculation in the wrong place. it would be better to do
this as a compute (or fix).
what you are discarding with your code is the fact that some of the
"j" atoms in your calculation
are "ghost atoms". please also note that the i,j loop is not over all
pairs of atoms but over all i,j pairs within the neighborlist cutoff.

to have the velocity of ghost atoms available, you also have to use the command
comm_modify vel yes

if you need all pairs, and not just the neighbor list, you need to do
some programming with a ring buffer to have all information of atoms,
that are on other processors pass by each processor.
otherwise, you can probably save some computational effort by using a
half neighbor list.

axel.

You do know that compute heat/flux works fine
already with many-body potentials such as Tersoff?
It includes the energy contribution from the 2-body
and 3-body terms. If you are trying to get only
the 2-body contribution, then I suggest you use
the “rerun” command, to read a dump file
(of positions and velocities), comment out the
3-body energy tallying in pair_tersoff.cpp, then
use compute heat/flux in the rerun script.

Steve