git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3590 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute.cpp

@@ -15,18 +15,21 @@
 #include "stdlib.h"
 #include "string.h"
 #include "ctype.h"
-#include "comm.h"
 #include "compute.h"
+#include "atom.h"
 #include "domain.h"
+#include "comm.h"
 #include "group.h"
-#include "modify.h"
-#include "lattice.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 #define DELTA 4
+#define BIG 2000000000
+
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
 
 /* ---------------------------------------------------------------------- */
 
@@ -77,6 +80,10 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   
   ntime = maxtime = 0;
   tlist = NULL;
+
+  // setup map for molecule IDs
+
+  molmap = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -87,6 +94,7 @@ Compute::~Compute()
   delete [] style;
 
   memory->sfree(tlist);
+  memory->sfree(molmap);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -189,3 +197,99 @@ void Compute::clearstep()
 {
   ntime = 0;
 }
+
+/* ----------------------------------------------------------------------
+   identify molecule IDs with atoms in group
+   warn if any atom in group has molecule ID = 0
+   warn if any molecule has only some atoms in group
+   return Ncount = # of molecules with atoms in group
+   set molmap to NULL if molecule IDs include all in range from 1 to Ncount
+   else: molecule IDs range from idlo to idhi
+         set molmap to vector of length idhi-idlo+1
+	 molmap[id-idlo] = index from 0 to Ncount-1
+         return idlo and idhi
+------------------------------------------------------------------------- */
+
+int Compute::molecules_in_group(int &idlo, int &idhi)
+{
+  int i;
+
+  memory->sfree(molmap);
+  molmap = NULL;
+
+  // find lo/hi molecule ID for any atom in group
+  // warn if atom in group has ID = 0
+
+  int *molecule = atom->molecule;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int lo = BIG;
+  int hi = -BIG;
+  int flag = 0;
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (molecule[i] == 0) {
+	flag = 1;
+	continue;
+      }
+      lo = MIN(lo,molecule[i]);
+      hi = MAX(hi,molecule[i]);
+    }
+
+  int flagall;
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+  if (flagall && comm->me == 0)
+    error->warning("Atom with molecule ID = 0 included in "
+		   "compute molecule group");
+
+  MPI_Allreduce(&lo,&idlo,1,MPI_INT,MPI_MIN,world);
+  MPI_Allreduce(&hi,&idhi,1,MPI_INT,MPI_MAX,world);
+  if (idlo == BIG) return 0;
+
+  // molmap = vector of length nlen
+  // set to 1 for IDs that appear in group across all procs, else 0
+
+  int nlen = idhi-idlo+1;
+  molmap = (int *) memory->smalloc(nlen*sizeof(int),"compute:molmap");
+  for (i = 0; i < nlen; i++) molmap[i] = 0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      molmap[molecule[i]-idlo] = 1;
+
+  int *molmapall = 
+    (int *) memory->smalloc(nlen*sizeof(int),"compute:molmapall");
+  MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world);
+
+  // nmolecules = # of non-zero IDs in molmap
+  // molmap[i] = index of molecule, skipping molecules not in group with -1
+
+  int nmolecules = 0;
+  for (i = 0; i < nlen; i++)
+    if (molmapall[i]) molmap[i] = nmolecules++;
+    else molmap[i] = -1;
+  memory->sfree(molmapall);
+
+  // warn if any molecule has some atoms in group and some not in group
+
+  flag = 0;
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) continue;
+    if (molecule[i] == 0) continue;
+    if (molecule[i] < idlo || molecule[i] > idhi) continue;
+    if (molmap[molecule[i]-idlo] >= 0) flag = 1;
+  }
+
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+  if (flagall && comm->me == 0)
+    error->warning("One or more compute molecules has atoms not in group");
+
+  // if molmap simply stores 1 to Nmolecules, then free it
+
+  if (nmolecules < nlen) return nmolecules;
+  if (idlo > 1) return nmolecules;
+  memory->sfree(molmap);
+  molmap = NULL;
+  return nmolecules;
+}

src/compute.h

@@ -118,6 +118,10 @@ class Compute : protected Pointers {
   double vbias[3];             // stored velocity bias for one atom
   double **vbiasall;           // stored velocity bias for all atoms
   int maxbias;                 // size of vbiasall array
+
+  int *molmap;                 // convert molecule ID to local index
+
+  int molecules_in_group(int &, int &);
 };
 
 }

src/compute_com_molecule.cpp

@@ -0,0 +1,156 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_com_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeCOMMolecule::ComputeCOMMolecule(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all("Illegal compute com/molecule command");
+
+  if (atom->molecular == 0)
+    error->all("Compute com/molecule requires molecular atom style");
+
+  array_flag = 1;
+  size_array_cols = 3;
+  extarray = 0;
+
+  // setup molecule-based data
+
+  nmolecules = molecules_in_group(idlo,idhi);
+  size_array_rows = nmolecules;
+
+  massproc = (double *) memory->smalloc(nmolecules*sizeof(double),
+					"com/molecule:massproc");
+  masstotal = (double *) memory->smalloc(nmolecules*sizeof(double),
+					 "com/molecule:masstotal");
+  com = memory->create_2d_double_array(nmolecules,3,"com/molecule:com");
+  comall = memory->create_2d_double_array(nmolecules,3,"com/molecule:comall");
+  array = comall;
+
+  // compute masstotal for each molecule
+
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  int i,imol;
+  double massone;
+
+  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      massproc[imol] += massone;
+    }
+
+  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeCOMMolecule::~ComputeCOMMolecule()
+{
+  memory->sfree(massproc);
+  memory->sfree(masstotal);
+  memory->destroy_2d_double_array(com);
+  memory->destroy_2d_double_array(comall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeCOMMolecule::init()
+{
+  int ntmp = molecules_in_group(idlo,idhi);
+  if (ntmp != nmolecules)
+    error->all("Molecule count changed in compute com/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeCOMMolecule::compute_array()
+{
+  int i,imol;
+  double xbox,ybox,zbox;
+  double massone;
+
+  invoked_array = update->ntimestep;
+
+  for (i = 0; i < nmolecules; i++)
+    com[i][0] = com[i][1] = com[i][2] = 0.0;
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      com[imol][0] += (x[i][0] + xbox*xprd) * massone;
+      com[imol][1] += (x[i][1] + ybox*yprd) * massone;
+      com[imol][2] += (x[i][2] + zbox*zprd) * massone;
+    }
+
+  MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
+		MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < nmolecules; i++) {
+    comall[i][0] /= masstotal[i];
+    comall[i][1] /= masstotal[i];
+    comall[i][2] /= masstotal[i];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeCOMMolecule::memory_usage()
+{
+  double bytes = 2*nmolecules * sizeof(double);
+  if (molmap) bytes += nmolecules * sizeof(int);
+  bytes += 2*nmolecules*3 * sizeof(double);
+  return bytes;
+}

src/compute_com_molecule.h

@@ -0,0 +1,39 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_COM_MOLECULE_H
+#define COMPUTE_COM_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeCOMMolecule : public Compute {
+ public:
+  ComputeCOMMolecule(class LAMMPS *, int, char **);
+  ~ComputeCOMMolecule();
+  void init();
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int nmolecules;
+  int idlo,idhi;
+
+  double *massproc,*masstotal;
+  double **com,**comall;
+};
+
+}
+
+#endif

src/compute_gyration_molecule.cpp

@@ -0,0 +1,186 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "compute_gyration_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGyrationMolecule::ComputeGyrationMolecule(LAMMPS *lmp, 
+						 int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all("Illegal compute gyration/molecule command");
+
+  if (atom->molecular == 0)
+    error->all("Compute gyration/molecule requires molecular atom style");
+
+  vector_flag = 1;
+  extvector = 0;
+
+  // setup molecule-based data
+
+  nmolecules = molecules_in_group(idlo,idhi);
+  size_vector = nmolecules;
+
+  massproc = (double *) memory->smalloc(nmolecules*sizeof(double),
+					"gyration/molecule:massproc");
+  masstotal = (double *) memory->smalloc(nmolecules*sizeof(double),
+					 "gyration/molecule:masstotal");
+  com = memory->create_2d_double_array(nmolecules,3,"gyration/molecule:com");
+  comall = memory->create_2d_double_array(nmolecules,3,
+					  "gyration/molecule:comall");
+  rg = (double *) memory->smalloc(nmolecules*sizeof(double),
+				  "gyration/molecule:rg");
+  rgall = (double *) memory->smalloc(nmolecules*sizeof(double),
+				     "gyration/molecule:rgall");
+  vector = rgall;
+
+  // compute masstotal for each molecule
+
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  int i,imol;
+  double massone;
+
+  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      massproc[imol] += massone;
+    }
+
+  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGyrationMolecule::~ComputeGyrationMolecule()
+{
+  memory->sfree(massproc);
+  memory->sfree(masstotal);
+  memory->destroy_2d_double_array(com);
+  memory->destroy_2d_double_array(comall);
+  memory->sfree(rg);
+  memory->sfree(rgall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGyrationMolecule::init()
+{
+  int ntmp = molecules_in_group(idlo,idhi);
+  if (ntmp != nmolecules)
+    error->all("Molecule count changed in compute gyration/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGyrationMolecule::compute_vector()
+{
+  int i,imol;
+  double xbox,ybox,zbox,dx,dy,dz;
+  double massone;
+
+  invoked_array = update->ntimestep;
+
+  for (i = 0; i < nmolecules; i++)
+    com[i][0] = com[i][1] = com[i][2] = 0.0;
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      com[imol][0] += (x[i][0] + xbox*xprd) * massone;
+      com[imol][1] += (x[i][1] + ybox*yprd) * massone;
+      com[imol][2] += (x[i][2] + zbox*zprd) * massone;
+    }
+
+  MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
+		MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < nmolecules; i++) {
+    comall[i][0] /= masstotal[i];
+    comall[i][1] /= masstotal[i];
+    comall[i][2] /= masstotal[i];
+  }
+
+  for (i = 0; i < nmolecules; i++) rg[i] = 0.0;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      dx = (x[i][0] + xbox*xprd) - comall[imol][0];
+      dy = (x[i][1] + ybox*yprd) - comall[imol][1];
+      dz = (x[i][2] + zbox*zprd) - comall[imol][2];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      rg[imol] += (dx*dx + dy*dy + dz*dz) * massone;
+    }
+
+  MPI_Allreduce(rg,rgall,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+
+  for (i = 0; i < nmolecules; i++) rgall[i] = sqrt(rgall[i]/masstotal[i]);
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeGyrationMolecule::memory_usage()
+{
+  double bytes = 4*nmolecules * sizeof(double);
+  if (molmap) bytes += nmolecules * sizeof(int);
+  bytes += 2*nmolecules*3 * sizeof(double);
+  return bytes;
+}

src/compute_gyration_molecule.h

@@ -0,0 +1,40 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_GYRATION_MOLECULE_H
+#define COMPUTE_GYRATION_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeGyrationMolecule : public Compute {
+ public:
+  ComputeGyrationMolecule(class LAMMPS *, int, char **);
+  ~ComputeGyrationMolecule();
+  void init();
+  void compute_vector();
+  double memory_usage();
+
+ private:
+  int nmolecules;
+  int idlo,idhi;
+
+  double *massproc,*masstotal;
+  double **com,**comall;
+  double *rg,*rgall;
+};
+
+}
+
+#endif

src/compute_msd_molecule.cpp

@@ -0,0 +1,190 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_msd_molecule.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeMSDMolecule::ComputeMSDMolecule(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all("Illegal compute msd/molecule command");
+
+  if (atom->molecular == 0)
+    error->all("Compute msd/molecule requires molecular atom style");
+
+  array_flag = 1;
+  size_array_cols = 4;
+  extarray = 0;
+
+  // setup molecule-based data and initial COM positions
+
+  nmolecules = molecules_in_group(idlo,idhi);
+  size_array_rows = nmolecules;
+
+  massproc = (double *) memory->smalloc(nmolecules*sizeof(double),
+					"msd/molecule:massproc");
+  masstotal = (double *) memory->smalloc(nmolecules*sizeof(double),
+					 "msd/molecule:masstotal");
+  com = memory->create_2d_double_array(nmolecules,3,"msd/molecule:com");
+  comall = memory->create_2d_double_array(nmolecules,3,"msd/molecule:comall");
+  cominit = memory->create_2d_double_array(nmolecules,3,
+					   "msd/molecule:cominit");
+  msd = memory->create_2d_double_array(nmolecules,4,"msd/molecule:msd");
+  array = msd;
+
+  // compute masstotal for each molecule
+
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  int i,imol;
+  double massone;
+
+  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      massproc[imol] += massone;
+    }
+
+  MPI_Allreduce(massproc,masstotal,nmolecules,MPI_DOUBLE,MPI_SUM,world);
+
+  // compute initial COM for each molecule
+
+  firstflag = 1;
+  compute_array();
+  for (i = 0; i < nmolecules; i++) {
+    cominit[i][0] = comall[i][0];
+    cominit[i][1] = comall[i][1];
+    cominit[i][2] = comall[i][2];
+  }
+  firstflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeMSDMolecule::~ComputeMSDMolecule()
+{
+  memory->sfree(massproc);
+  memory->sfree(masstotal);
+  memory->destroy_2d_double_array(com);
+  memory->destroy_2d_double_array(comall);
+  memory->destroy_2d_double_array(cominit);
+  memory->destroy_2d_double_array(msd);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeMSDMolecule::init()
+{
+  int ntmp = molecules_in_group(idlo,idhi);
+  if (ntmp != nmolecules)
+    error->all("Molecule count changed in compute msd/molecule");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeMSDMolecule::compute_array()
+{
+  int i,imol;
+  double xbox,ybox,zbox,dx,dy,dz;
+  double massone;
+
+  invoked_array = update->ntimestep;
+
+  // compute current COM positions
+
+  for (i = 0; i < nmolecules; i++)
+    com[i][0] = com[i][1] = com[i][2] = 0.0;
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *molecule = atom->molecule;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = (image[i] & 1023) - 512;
+      ybox = (image[i] >> 10 & 1023) - 512;
+      zbox = (image[i] >> 20) - 512;
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      imol = molecule[i];
+      if (molmap) imol = molmap[imol-idlo];
+      else imol--;
+      com[imol][0] += (x[i][0] + xbox*xprd) * massone;
+      com[imol][1] += (x[i][1] + ybox*yprd) * massone;
+      com[imol][2] += (x[i][2] + zbox*zprd) * massone;
+    }
+
+  MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
+		MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < nmolecules; i++) {
+    comall[i][0] /= masstotal[i];
+    comall[i][1] /= masstotal[i];
+    comall[i][2] /= masstotal[i];
+  }
+
+  // MSD is difference between current and initial COM
+  // cominit does not yet exist when called from constructor
+
+  if (firstflag) return;
+
+  for (i = 0; i < nmolecules; i++) {
+    dx = comall[i][0] - cominit[i][0];
+    dy = comall[i][1] - cominit[i][1];
+    dz = comall[i][2] - cominit[i][2];
+    msd[i][0] = dx*dx;
+    msd[i][1] = dy*dy;
+    msd[i][2] = dz*dz;
+    msd[i][3] = dx*dx + dy*dy + dz*dz;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeMSDMolecule::memory_usage()
+{
+  double bytes = 2*nmolecules * sizeof(double);
+  if (molmap) bytes += nmolecules * sizeof(int);
+  bytes += 2*nmolecules*3 * sizeof(double);
+  bytes += nmolecules*4 * sizeof(double);
+  return bytes;
+}

src/compute_msd_molecule.h

@@ -0,0 +1,41 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_MSD_MOLECULE_H
+#define COMPUTE_MSD_MOLECULE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeMSDMolecule : public Compute {
+ public:
+  ComputeMSDMolecule(class LAMMPS *, int, char **);
+  ~ComputeMSDMolecule();
+  void init();
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int nmolecules;
+  int idlo,idhi;
+  int firstflag;
+
+  double *massproc,*masstotal;
+  double **com,**comall,**cominit;
+  double **msd;
+};
+
+}
+
+#endif

src/style.h

@@ -81,16 +81,19 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_centro_atom.h"
 #include "compute_cna_atom.h"
 #include "compute_com.h"
+#include "compute_com_molecule.h"
 #include "compute_coord_atom.h"
 #include "compute_dihedral_local.h"
 #include "compute_displace_atom.h"
 #include "compute_group_group.h"
 #include "compute_gyration.h"
+#include "compute_gyration_molecule.h"
 #include "compute_heat_flux.h"
 #include "compute_improper_local.h"
 #include "compute_ke.h"
 #include "compute_ke_atom.h"
 #include "compute_msd.h"
+#include "compute_msd_molecule.h"
 #include "compute_pe.h"
 #include "compute_pe_atom.h"
 #include "compute_pressure.h"
@@ -117,16 +120,19 @@ ComputeStyle(bond/local,ComputeBondLocal)
 ComputeStyle(centro/atom,ComputeCentroAtom)
 ComputeStyle(cna/atom,ComputeCNAAtom)
 ComputeStyle(com,ComputeCOM)
+ComputeStyle(com/molecule,ComputeCOMMolecule)
 ComputeStyle(coord/atom,ComputeCoordAtom)
 ComputeStyle(dihedral/local,ComputeDihedralLocal)
 ComputeStyle(displace/atom,ComputeDisplaceAtom)
 ComputeStyle(group/group,ComputeGroupGroup)
 ComputeStyle(gyration,ComputeGyration)
+ComputeStyle(gyration/molecule,ComputeGyrationMolecule)
 ComputeStyle(heat/flux,ComputeHeatFlux)
 ComputeStyle(improper/local,ComputeImproperLocal)
 ComputeStyle(ke,ComputeKE)
 ComputeStyle(ke/atom,ComputeKEAtom)
 ComputeStyle(msd,ComputeMSD)
+ComputeStyle(msd/molecule,ComputeMSDMolecule)
 ComputeStyle(pe,ComputePE)
 ComputeStyle(pe/atom,ComputePEAtom)
 ComputeStyle(pressure,ComputePressure)