[lammps-users] new compute command

LAMMPS LAMMPS Mailing List Mirror
1

Dear Steve,

I am trying to modify the simple compute_ke.cpp code for writing a new
compute command later.
So I replace the v[i][0],v[i][1] and v[i][2] with x[i][0],x[i][1] and
x[i][2] and then I compare the result with 0.5*m*(x^2 + y^2 +z^2) from
dump file. When I use a single processing, they are the same and when I
use multiprocessing, they are different. Would you please tell me how I
can get the same values? Here’s my modified code.

Regards,
Ali

#include "mpi.h"
#include "compute_ke.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "group.h"
#include "error.h"

using namespace LAMMPS_NS;

#define INVOKED_SCALAR 1

/* ---------------------------------------------------------------------- */

ComputeKE_new::ComputeKE_new(LAMMPS *lmp, int narg, char **arg) :
  Compute(lmp, narg, arg)
{
  if (narg != 3) error->all("Illegal compute ke command");

  scalar_flag = 1;
  extscalar = 1;
}

/* ---------------------------------------------------------------------- */

void ComputeKE_new::init()
{
  pfactor = 0.5 * force->mvv2e;
}

/* ---------------------------------------------------------------------- */

double ComputeKE_new::compute_scalar()
{
  invoked |= INVOKED_SCALAR;

  double **v = atom->v;
  double *rmass = atom->rmass;
  double *mass = atom->mass;
  int *mask = atom->mask;
  int *type = atom->type;
  int nlocal = atom->nlocal;

  double ke_new = 0.0;

  if (mass) {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit)
  ke_new += mass[type[i]] *
    (x[i][0]*x[i][0] + x[i][1]*x[i][1] + x[i][2]*x[i][2]);

  } else {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit)
  ke_new += rmass[i] * (x[i][0]*x[i][0] + x[i][1]*x[i][1] +
x[i][2]*x[i][2]);
  }

  MPI_Allreduce(&ke_new,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
  scalar *= pfactor;
  return scalar;
}

2

When you say compare to the dump file, you mean you are summing
up all the per-atom values in the dump file and comparing
to what this compute outputs? If so, it seems like they
would give the same answer.

Steve

3

Thanks for your reply. I made some mistakes and you are right.
Now I am trying to write a new compute command for computing rij^2 over a
few particles.
I have modified the compute_ke.cpp as following but the answer is wrong
with parallel processing.
Would you please tell me how I can write a parallel compute command for
this computation?

#include "mpi.h"
#include "compute_rij.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "group.h"
#include "error.h"

using namespace LAMMPS_NS;

#define INVOKED_SCALAR 1

/* ---------------------------------------------------------------------- */

Compute_rij::Compute_rij(LAMMPS *lmp, int narg, char **arg) :
  Compute(lmp, narg, arg)
{
  if (narg != 3) error->all("Illegal compute rij command");

  scalar_flag = 1;
  extscalar = 1;
}

/* ---------------------------------------------------------------------- */

void Compute_rij::init()
{
  pfactor = 0.5* force->mvv2e;
}

/* ---------------------------------------------------------------------- */

double Compute_rij::compute_scalar()
{
  invoked |= INVOKED_SCALAR;

  double **v = atom->v;
  double **x = atom->x;
  double **f = atom->f;
  double *rmass = atom->rmass;
  double *mass = atom->mass;
  int *mask = atom->mask;
  int *type = atom->type;
  int nlocal = atom->nlocal;

  double rij = 0.0;

for (int j=0; j < nlocal; j++){
   if (mask[j] & groupbit)
     for ( int i=0; i < nlocal; i++){
         if (mask[i] & groupbit)
            rij += (x[i][0]-x[j][0])*(x[i][0]-x[j][0])
                     +(x[i][1]-x[j][1])*(x[i][1]-x[j][1])
                     +(x[i][2]-x[j][2])*(x[i][2]-x[j][2])
    }
}

  MPI_Allreduce(rij,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
  scalar *= pfactor;
  return scalar;
}

When you say compare to the dump file, you mean you are summing up all
the per-atom values in the dump file and comparing
to what this compute outputs? If so, it seems like they
would give the same answer.

Steve

Dear Steve,
I am trying to modify the simple compute_ke.cpp code for writing a new

compute command later.

So I replace the v[i][0],v[i][1] and v[i][2] with x[i][0],x[i][1] and

x[i][2] and then I compare the >>>result with 0.5*m*(x^2 + y^2 +z^2)
from dump file. When I use a single processing, they are the

same and when I
use multiprocessing, they are different. Would you please tell me how

I can get the same values? >>>Here's my modified code.

4

ali,

Thanks for your reply. I made some mistakes and you are right.

Now I am trying to write a new compute command for computing rij^2
over a few particles. I have modified the compute_ke.cpp as
following but the answer is wrong with parallel processing.

yes. it _has_ to be. you don't factor in that lammps
does domain decomposition and thus not all nodes have
access to all coordinates and that particularly the
number of "local" coordinates can be different after
each rebuild of the neighbor lists.

the compute_ke.cpp code is not a very good template for
your calculation, since it sums over a 'per atom' property,
whereas you are after a 'per atom pair' property.

you only compute the pairs that are local to each node.
that works in serial but fails in parallel, since you
are missing the 'between nodes' pairs.

Would you please tell me how I can write a parallel compute command for
this computation?

what you want to do is really best done in postprocessing.
to understand why, i suggest you first have a look at the
papers describing how lammps is parallelized particularly
the domain decomposition part.

cheers,
    axel.

5

I can't help you debug your code. I'd run it on a 1 atom system,
then a 2 atom system, etc and figure out what's wrong.

Steve