Energy moment needs ghost atoms' position, forces and velocities

Xiaoliang_Zhang1 · January 28, 2016, 1:22pm

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~

akohlmey · January 28, 2016, 1:36pm

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?

there is a more fundamental problem here. you don't seem to understand
the nature of domain decomposition and the fact that with your scheme
there is a lot of double counting going on, and that is changing with
communication cutoff and number of processors.

if all you need to compute is a \sum_i from some expression over only
i-indexed properties for all atoms in the system, then you don't need
to worry about ghost atoms at all. just loop over the nlocal atoms
owned by each subdomain and then do the reduction over all MPI ranks
to combine all subdomains into the total.

you have not mentioned, however, what determines the "unique N-body
interaction groups".

furthermore, since you are not using neighbor lists, i don't
understand the purpose of using the "includegroup" feature. on the
other hand, computes typically operate on a group, but there is no
checking of the group bitmask in your code.

axel.

Xiaoliang_Zhang1 · January 28, 2016, 2:21pm

Dear Axel,

Thanks for your kind help.

I think the summation sum_i ( r_i * (f_i . v_i) ) is over all local and ghost atoms.

If i is a ghost atom then the position of the ghost atom is used, not the corresponding local atom. And the force f_i is the force on the ghost atom i due to all interactions.

Is it possible to realize this? Thanks.

Best regards,
Xiaoliang

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?

there is a more fundamental problem here. you don't seem to understand
the nature of domain decomposition and the fact that with your scheme
there is a lot of double counting going on, and that is changing with
communication cutoff and number of processors.

if all you need to compute is a \sum_i from some expression over only
i-indexed properties for all atoms in the system, then you don't need
to worry about ghost atoms at all. just loop over the nlocal atoms
owned by each subdomain and then do the reduction over all MPI ranks
to combine all subdomains into the total.

you have not mentioned, however, what determines the "unique N-body
interaction groups".

furthermore, since you are not using neighbor lists, i don't
understand the purpose of using the "includegroup" feature. on the
other hand, computes typically operate on a group, but there is no
checking of the group bitmask in your code.

I just copied the codes from pair.cpp for the void Pair::virial_fdotr_compute() subroutine, and modified the virial codes to be my codes. Yes, maybe I don't need the "includegroup" feature, because what I need is the energy moment for all atoms.

akohlmey · January 28, 2016, 2:38pm

Dear Axel,

Thanks for your kind help.

I think the summation sum_i ( r_i * (f_i . v_i) ) is over all local and
ghost atoms.

why do you ask for help and then don't pay any attention to what you get told?

If i is a ghost atom then the position of the ghost atom is used, not the
corresponding local atom. And the force f_i is the force on the ghost atom i
due to all interactions.

as i already mentioned, that makes no sense. *every* ghost atom, is a
local atom *somewhere*.

Is it possible to realize this? Thanks.

you may compute anything you want. but don't blame me, if your results are crap.

sjplimp · January 28, 2016, 3:34pm

A compute can acquire info about ghost atoms

by invoking comm->forward_comm_compute()

and providing the needed pack() and unpack() methods.

However, it looks like you doing a variant of sum (local+ghost) r dot F

which I don’t think will work at the end of the timestep

when your compute is likely to be invoked. The typical

reason to do that is include ghost atoms to avoid accounting

for periodic boundaries. At the end of the timestep the

forces on local atoms have already been summed from

ghost atoms, so the summation will not be correct.

There is a relatively new mechanism in the code to

add hooks to the force calculation and compute new

quantities like virial variants which happens at the

correct point during the timestep to avoid this

problem. Axel can elaborate.

Steve

akohlmey · January 29, 2016, 10:45am

A compute can acquire info about ghost atoms
by invoking comm->forward_comm_compute()
and providing the needed pack() and unpack() methods.

However, it looks like you doing a variant of sum (local+ghost) r dot F
which I don't think will work at the end of the timestep
when your compute is likely to be invoked. The typical
reason to do that is include ghost atoms to avoid accounting
for periodic boundaries. At the end of the timestep the
forces on local atoms have already been summed from
ghost atoms, so the summation will not be correct.

There is a relatively new mechanism in the code to
add hooks to the force calculation and compute new
quantities like virial variants which happens at the
correct point during the timestep to avoid this
problem. Axel can elaborate.

steve refers here to the USER-TALLY package, which comes with several
example computes.

http://lammps.sandia.gov/doc/compute_tally.html

there is a special complication for computes that need to access
velocities. during the force compute, that is when the tally callback
is executed, the velocities are not fully updated, so their
contributions must be computed (or corrected) later. compute
heat/flux/tally is an example for that.

axel.

athomps · February 2, 2016, 10:46pm

I looked at the paper[1], but found it to be a little short on detail. I can see the advantages of using an Einstein formulation, and this is already well established for viscosity calculations, although it is not available in LAMMPS [2] The proposed method for the time-integrated heat flux (energy moment) seems to have some similarity to the way that LAMMPS calculates the stress tensor for arbitrary manybody potentials using only partial forces accumulated on local and ghost atoms before ghost atom forces are applied to the local atoms by reverse communication. [3] Getting LAMMPS to do the same thing for the time-integrated heat flux would be non-trivial. Before embarking on an effort to implement this in LAMMPS, you should first try the methods already available in LAMMPS, under examples/KAPPA.

Aidan

[1] Kinaci, Haskins, Cagin, J.Chem. Phys. 137, 014106 (2012)
[2] Mondello, M Grest, GS, JChemPhys, 106, 9327 (1997)
[3] Thompson, Plimpton, Mattson, J Chem Phys, 131, 154107 (2009)

*They make an unfortunate choice of nomenclature, using the letter ‘k’ and the word “kinetic” to refer to what is normally referred to as the virial, force, or mechanical work contribution to the heat flux, while using the letter ‘p’ and “potential” to refer to what is normally called the convective or kinetic contribution.