Dear all,

I'm interested in calculating the energy moment which is related to Rk = sum_i ( r_i * (f_i . v_i) ), where i is a summation over all particles (real and ghost atoms) that make up unique N-body interaction groups (J. Chem. Phys. 137, 014106 (2012)).

I write some codes to realize this:

void ComputeEnergyMoment::compute_vector()

{

invoked_vector = update->ntimestep;

double **x = atom->x;

double **f = atom->f;

double **v = atom->v;

// compute energy moment

// heat flux vector = Rk[3]

// Rk[alpha] = sum_i( x_i[alpha]*(f_i[0]*v_i[0]+f_i[1]*v_i[1]+f_i[2]*v_i[2]) )

double Rk[3] = {0.0,0.0,0.0};

// sum over force on all particles including ghosts

if (neighbor->includegroup == 0) {

int nall = atom->nlocal + atom->nghost;

for (int i = 0; i < nall; i++) {

Rk[0] += x[i][0] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[1] += x[i][1] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[2] += x[i][2] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

}

// neighbor includegroup flag is set

// sum over force on initial nfirst particles and ghosts

} else {

int nall = atom->nfirst;

for (int i = 0; i < nall; i++) {

Rk[0] += x[i][0] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[1] += x[i][1] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[2] += x[i][2] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

}

nall = atom->nlocal + atom->nghost;

for (int i = atom->nlocal; i < nall; i++) {

Rk[0] += x[i][0] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[1] += x[i][1] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

Rk[2] += x[i][2] * (f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]);

}

}

MPI_Allreduce(Rk,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);

}

But it seems that it's wrong, and I don't know how to achieve the correct results because the formula Rk = sum_i ( r_i * (f_i . v_i) ) needs not only the positions, forces and velocities of the local atoms but also the ghost atoms. Can anyone give me some help about how to write the codes to realize the formula?

Thanks a lot and have a nice day~