/* ----------------------------------------------------------------------
   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 "mpi.h"
#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "group.h"
#include "domain.h"
#include "atom.h"
#include "atom_vec_hybrid.h"
#include "region.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define MAX_GROUP 32

#define TYPE      1
#define MOLECULE  2
#define ID        3

#define LESS_THAN      1
#define GREATER_THAN   2
#define LESS_EQUAL     3
#define GREATER_EQUAL  4
#define BETWEEN        5

#define BIG 1.0e20

#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))

/* ----------------------------------------------------------------------
   initialize group memory
------------------------------------------------------------------------- */

Group::Group(LAMMPS *lmp) : Pointers(lmp)
{
  MPI_Comm_rank(world,&me);

  names = (char **) memory->smalloc(MAX_GROUP*sizeof(char *),"group:names");
  bitmask = new int[MAX_GROUP];
  inversemask = new int[MAX_GROUP];
  for (int i = 0; i < MAX_GROUP; i++) bitmask[i] = 1 << i;
  for (int i = 0; i < MAX_GROUP; i++) inversemask[i] = bitmask[i] ^ ~0;

  // create "all" group

  char *str = "all";
  int n = strlen(str) + 1;
  names[0] = (char *) memory->smalloc(n*sizeof(char),"group:names[]");
  strcpy(names[0],str);
  ngroup = 1;
}

/* ----------------------------------------------------------------------
   free all memory
------------------------------------------------------------------------- */

Group::~Group()
{
  for (int i = 0; i < ngroup; i++) memory->sfree(names[i]);
  memory->sfree(names);
  delete [] bitmask;
  delete [] inversemask;
}

/* ----------------------------------------------------------------------
   assign atoms to a new or existing group
------------------------------------------------------------------------- */

void Group::assign(int narg, char **arg)
{
  int i;

  if (domain->box_exist == 0) 
    error->all("Group command before simulation box is defined");
  if (narg < 2) error->all("Illegal group command");

  // find group in existing list
  // igroup = -1 is a new group name, add it

  int igroup = find(arg[0]);

  if (igroup == -1) {
    if (ngroup == MAX_GROUP) error->all("Too many groups");
    igroup = ngroup;
    ngroup++;
    int n = strlen(arg[0]) + 1;
    names[igroup] = (char *) memory->smalloc(n*sizeof(char),"group:names[]");
    strcpy(names[igroup],arg[0]);
  }

  double **x = atom->x;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  int bit = bitmask[igroup];
    
  // style = region

  if (strcmp(arg[1],"region") == 0) {

    if (narg != 3) error->all("Illegal group command");
    
    int iregion;
    for (iregion = 0; iregion < domain->nregion; iregion++) 
      if (strcmp(arg[2],domain->regions[iregion]->id) == 0) break;
    if (iregion == domain->nregion)
      error->all("Group region ID does not exist");

    // add to group if atom is in region

    for (i = 0; i < nlocal; i++)
      if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
	mask[i] |= bit;

  // style = hybrid

  } else if (strcmp(arg[1],"hybrid") == 0) {

    if (narg != 3) error->all("Illegal group command");
    if (strcmp(atom->atom_style,"hybrid") != 0)
      error->all("Group hybrid command requires atom style hybrid");
    
    AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) atom->avec;
    int ihybrid;
    for (ihybrid = 0; ihybrid < avec_hybrid->nstyles; ihybrid++)
      if (strcmp(avec_hybrid->keywords[ihybrid],arg[2]) == 0) break;
    if (ihybrid == avec_hybrid->nstyles)
      error->all("Group hybrid sub-style does not exist");

    // add to group if atom matches hybrid sub-style

    int *hybrid = atom->hybrid;

    for (i = 0; i < nlocal; i++)
      if (hybrid[i] == ihybrid)	mask[i] |= bit;

  // style = logical condition

  } else if (strcmp(arg[2],"<") == 0 || strcmp(arg[2],">") == 0 || 
	     strcmp(arg[2],"<=") == 0 || strcmp(arg[2],">=") == 0 || 
	     strcmp(arg[2],"<>") == 0) {

    if (narg < 4 || narg > 5) error->all("Illegal group command");
    int category,condition,bound1,bound2;

    if (strcmp(arg[1],"type") == 0) category = TYPE;
    else if (strcmp(arg[1],"molecule") == 0) category = MOLECULE;
    else if (strcmp(arg[1],"id") == 0) category = ID;
    else error->all("Illegal group command");
    
    if (strcmp(arg[2],"<") == 0) condition = LESS_THAN;
    else if (strcmp(arg[2],">") == 0) condition = GREATER_THAN;
    else if (strcmp(arg[2],"<=") == 0) condition = LESS_EQUAL;
    else if (strcmp(arg[2],">=") == 0) condition = GREATER_EQUAL;
    else if (strcmp(arg[2],"<>") == 0) condition = BETWEEN;
    else error->all("Illegal group command");

    bound1 = atoi(arg[3]);
    bound2 = -1;

    if (condition == BETWEEN) {
      if (narg != 5) error->all("Illegal group command");
      bound2 = atoi(arg[4]);
    }

    int *attribute;
    if (category == TYPE) attribute = atom->type;
    else if (category == MOLECULE) attribute = atom->molecule;
    else if (category == ID) attribute = atom->tag;

    // add to group if meets condition

    if (condition == LESS_THAN) {
      for (i = 0; i < nlocal; i++) if (attribute[i] < bound1) mask[i] |= bit;
    } else if (condition == GREATER_THAN) {
      for (i = 0; i < nlocal; i++) if (attribute[i] > bound1) mask[i] |= bit;
    } else if (condition == LESS_EQUAL) {
      for (i = 0; i < nlocal; i++) if (attribute[i] <= bound1) mask[i] |= bit;
    } else if (condition == GREATER_EQUAL) {
      for (i = 0; i < nlocal; i++) if (attribute[i] >= bound1) mask[i] |= bit;
    } else if (condition == BETWEEN) {
      for (i = 0; i < nlocal; i++)
	if (attribute[i] >= bound1 && attribute[i] <= bound2) mask[i] |= bit;
    }

  // style = list of values

  } else if (strcmp(arg[1],"type") == 0 || strcmp(arg[1],"molecule") == 0 ||
	     strcmp(arg[1],"id") == 0) {

    if (narg < 3) error->all("Illegal group command");

    int length = narg-2;
    int *list = new int[length];

    int category;
    if (strcmp(arg[1],"type") == 0) category = TYPE;
    else if (strcmp(arg[1],"molecule") == 0) category = MOLECULE;
    else if (strcmp(arg[1],"id") == 0) category = ID;
    else error->all("Illegal group command");
    
    length = narg - 2;
    for (int iarg = 2; iarg < narg; iarg++) list[iarg-2] = atoi(arg[iarg]);

    int *attribute;
    if (category == TYPE) attribute = atom->type;
    else if (category == MOLECULE) attribute = atom->molecule;
    else if (category == ID) attribute = atom->tag;

    // add to group if attribute is any in list

    for (int ilist = 0; ilist < length; ilist++)
      for (i = 0; i < nlocal; i++)
	if (attribute[i] == list[ilist]) mask[i] |= bit;

    delete [] list;
    
  // style = subtract

  } else if (strcmp(arg[1],"subtract") == 0) {

    if (narg < 4) error->all("Illegal group command");
    
    int length = narg-2;
    int *list = new int[length];

    int jgroup;
    for (int iarg = 2; iarg < narg; iarg++) {
      jgroup = find(arg[iarg]);
      if (jgroup == -1) error->all("Group ID does not exist");
      list[iarg-2] = jgroup;
    }

    // add to group if in 1st group in list

    int otherbit = bitmask[list[0]];

    for (i = 0; i < nlocal; i++) 
      if (mask[i] & otherbit) mask[i] |= bit;

    // remove atoms if they are in any of the other groups
    // AND with inverse mask removes the atom from group

    int inverse = inversemask[igroup];

    for (int ilist = 1; ilist < length; ilist++) {
      otherbit = bitmask[list[ilist]];
      for (i = 0; i < nlocal; i++) 
	if (mask[i] & otherbit) mask[i] &= inverse;
    }

    delete [] list;

  // style = union

  } else if (strcmp(arg[1],"union") == 0) {

    if (narg < 3) error->all("Illegal group command");
    
    int length = narg-2;
    int *list = new int[length];

    int jgroup;
    for (int iarg = 2; iarg < narg; iarg++) {
      jgroup = find(arg[iarg]);
      if (jgroup == -1) error->all("Group ID does not exist");
      list[iarg-2] = jgroup;
    }

    // add to group if in any other group in list

    int otherbit;

    for (int ilist = 0; ilist < length; ilist++) {
      otherbit = bitmask[list[ilist]];
      for (i = 0; i < nlocal; i++) 
	if (mask[i] & otherbit) mask[i] |= bit;
    }

    delete [] list;

  // style = intersect

  } else if (strcmp(arg[1],"intersect") == 0) {

    if (narg < 4) error->all("Illegal group command");
    
    int length = narg-2;
    int *list = new int[length];

    int jgroup;
    for (int iarg = 2; iarg < narg; iarg++) {
      jgroup = find(arg[iarg]);
      if (jgroup == -1) error->all("Group ID does not exist");
      list[iarg-2] = jgroup;
    }

    // add to group if in all groups in list

    int otherbit,ok,ilist;

    for (i = 0; i < nlocal; i++) {
      ok = 1;
      for (ilist = 0; ilist < length; ilist++) {
	otherbit = bitmask[list[ilist]];
	if ((mask[i] & otherbit) == 0) ok = 0;
      }
      if (ok) mask[i] |= bit;
    }

    delete [] list;

  // not a valid group style

  } else error->all("Illegal group command");

  // print stats for changed group

  int n,all;
  n = 0;
  for (i = 0; i < nlocal; i++) if (mask[i] & bit) n++;
  MPI_Allreduce(&n,&all,1,MPI_INT,MPI_SUM,world);

  if (me == 0) {
    if (screen) fprintf(screen,"%d atoms in group %s\n",all,names[igroup]);
    if (logfile) fprintf(logfile,"%d atoms in group %s\n",all,names[igroup]);
  }
}
 
/* ----------------------------------------------------------------------
   add flagged atoms to a new or existing group
------------------------------------------------------------------------- */

void Group::create(char *name, int *flag)
{
  int i;

  // find group in existing list
  // igroup = -1 is a new group name, add it

  int igroup = find(name);

  if (igroup == -1) {
    if (ngroup == MAX_GROUP) error->all("Too many groups");
    igroup = ngroup;
    ngroup++;
    int n = strlen(name) + 1;
    names[igroup] = (char *) memory->smalloc(n*sizeof(char),"group:names[]");
    strcpy(names[igroup],name);
  }

  // add atom to group if flag is set

  int *mask = atom->mask;
  int nlocal = atom->nlocal;
  int bit = bitmask[igroup];

  for (i = 0; i < nlocal; i++) 
    if (flag[i]) mask[i] |= bit;
}

/* ----------------------------------------------------------------------
   return group index if name matches existing group, -1 if no such group
------------------------------------------------------------------------- */

int Group::find(char *name)
{
  for (int igroup = 0; igroup < ngroup; igroup++)
    if (strcmp(name,names[igroup]) == 0) return igroup;
  return -1;
}

/* ----------------------------------------------------------------------
   write group info to a restart file
   only called by proc 0
------------------------------------------------------------------------- */

void Group::write_restart(FILE *fp)
{
  fwrite(&ngroup,sizeof(int),1,fp);

  int n;
  for (int i = 0; i < ngroup; i++) {
    n = strlen(names[i]) + 1;
    fwrite(&n,sizeof(int),1,fp);
    fwrite(names[i],sizeof(char),n,fp);
  }
}

/* ----------------------------------------------------------------------
   read group info from a restart file
   proc 0 reads, bcast to all procs
------------------------------------------------------------------------- */

void Group::read_restart(FILE *fp)
{
  int i,n;

  for (i = 0; i < ngroup; i++) memory->sfree(names[i]);

  if (me == 0) fread(&ngroup,sizeof(int),1,fp);
  MPI_Bcast(&ngroup,1,MPI_INT,0,world);
  
  for (i = 0; i < ngroup; i++) {
    if (me == 0) fread(&n,sizeof(int),1,fp);
    MPI_Bcast(&n,1,MPI_INT,0,world);
    names[i] = (char *) memory->smalloc(n*sizeof(char),"group:names[]");
    if (me == 0) fread(names[i],sizeof(char),n,fp);
    MPI_Bcast(names[i],n,MPI_CHAR,0,world);
  }
}

// ----------------------------------------------------------------------
// computations on a group of atoms
// ----------------------------------------------------------------------

/* ----------------------------------------------------------------------
   count atoms in group
   compute in double precision in case system is huge
------------------------------------------------------------------------- */

double Group::count(int igroup)
{
  int groupbit = bitmask[igroup];

  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  int n = 0;
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) n++;

  double nsingle = n;
  double nall;
  MPI_Allreduce(&nsingle,&nall,1,MPI_DOUBLE,MPI_SUM,world);
  return nall;
}

/* ----------------------------------------------------------------------
   compute the total mass of group of atoms
   use either per-type mass or per-atom rmass
------------------------------------------------------------------------- */

double Group::mass(int igroup)
{
  int groupbit = bitmask[igroup];

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

  double m = 0.0;

  if (mass) {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) m += mass[type[i]];
  } else {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) m += rmass[i];
  }

  double mall;
  MPI_Allreduce(&m,&mall,1,MPI_DOUBLE,MPI_SUM,world);
  return mall;
}

/* ----------------------------------------------------------------------
   compute the total charge of group of atoms
------------------------------------------------------------------------- */

double Group::charge(int igroup)
{
  int groupbit = bitmask[igroup];

  int *mask = atom->mask;
  double *q = atom->q;
  int nlocal = atom->nlocal;

  double qone = 0.0;
  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) qone += q[i];

  double qall;
  MPI_Allreduce(&qone,&qall,1,MPI_DOUBLE,MPI_SUM,world);
  return qall;
}

/* ----------------------------------------------------------------------
   compute the coordinate bounds of the group of atoms
   periodic images are not considered, so atoms are NOT unwrapped
------------------------------------------------------------------------- */

void Group::bounds(int igroup, double *minmax)
{
  int groupbit = bitmask[igroup];

  double extent[6];
  extent[0] = extent[2] = extent[4] = BIG;
  extent[1] = extent[3] = extent[5] = -BIG;
  
  double **x = atom->x;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++) {
    if (mask[i] & groupbit) {
      extent[0] = MIN(extent[0],x[i][0]);
      extent[1] = MAX(extent[1],x[i][0]);
      extent[2] = MIN(extent[2],x[i][1]);
      extent[3] = MAX(extent[3],x[i][1]);
      extent[4] = MIN(extent[4],x[i][2]);
      extent[5] = MAX(extent[5],x[i][2]);
    }
  }
  
  // compute extent across all procs
  // flip sign of MIN to do it in one Allreduce MAX
  // set box by extent in shrink-wrapped dims
  
  extent[0] = -extent[0];
  extent[2] = -extent[2];
  extent[4] = -extent[4];
  
  MPI_Allreduce(extent,minmax,6,MPI_DOUBLE,MPI_MAX,world);

  extent[0] = -extent[0];
  extent[2] = -extent[2];
  extent[4] = -extent[4];
}

/* ----------------------------------------------------------------------
   compute the center-of-mass coords of group with total mass = masstotal
   return center-of-mass coords in cm[]
   must unwrap atoms to compute center-of-mass correctly
------------------------------------------------------------------------- */

void Group::xcm(int igroup, double masstotal, double *cm)
{
  int groupbit = bitmask[igroup];

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

  double cmone[3];
  cmone[0] = cmone[1] = cmone[2] = 0.0;
  double xprd = domain->xprd;
  double yprd = domain->yprd;
  double zprd = domain->zprd;

  int xbox,ybox,zbox;
  double massone;

  if (mass) {
    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;
	massone = mass[type[i]];
	cmone[0] += (x[i][0] + xbox*xprd) * massone;
	cmone[1] += (x[i][1] + ybox*yprd) * massone;
	cmone[2] += (x[i][2] + zbox*zprd) * massone;
      }
  } else {
    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;
	massone = rmass[i];
	cmone[0] += (x[i][0] + xbox*xprd) * massone;
	cmone[1] += (x[i][1] + ybox*yprd) * massone;
	cmone[2] += (x[i][2] + zbox*zprd) * massone;
      }
  }

  MPI_Allreduce(cmone,cm,3,MPI_DOUBLE,MPI_SUM,world);
  cm[0] /= masstotal;
  cm[1] /= masstotal;
  cm[2] /= masstotal;
}

/* ----------------------------------------------------------------------
   compute the center-of-mass velocity of group with total mass = masstotal
   return center-of-mass velocity in cm[]
------------------------------------------------------------------------- */

void Group::vcm(int igroup, double masstotal, double *cm)
{
  int groupbit = bitmask[igroup];

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

  double p[3],massone;
  p[0] = p[1] = p[2] = 0.0;

  if (mass) {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) {
	massone = mass[type[i]];
	p[0] += v[i][0]*massone;
	p[1] += v[i][1]*massone;
	p[2] += v[i][2]*massone;
      }
  } else {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) {
	massone = rmass[i];
	p[0] += v[i][0]*massone;
	p[1] += v[i][1]*massone;
	p[2] += v[i][2]*massone;
      }
  }

  MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world);
  cm[0] /= masstotal;
  cm[1] /= masstotal;
  cm[2] /= masstotal;
}

/* ----------------------------------------------------------------------
   compute the radius-of-gyration of group around center-of-mass cm
   must unwrap atoms to compute Rg correctly
------------------------------------------------------------------------- */

double Group::gyration(int igroup, double masstotal, double *cm)
{
  int groupbit = bitmask[igroup];

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

  int xbox,ybox,zbox;
  double dx,dy,dz,massone;
  double xprd = domain->xprd;
  double yprd = domain->yprd;
  double zprd = domain->zprd;
  double rg = 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;
      dx = (x[i][0] + xbox*xprd) - cm[0];
      dy = (x[i][1] + ybox*yprd) - cm[1];
      dz = (x[i][2] + zbox*zprd) - cm[2];
      if (mass) massone = mass[type[i]];
      else massone = rmass[i];
      rg += (dx*dx + dy*dy + dz*dz) * massone;
    }
  double rg_all;
  MPI_Allreduce(&rg,&rg_all,1,MPI_DOUBLE,MPI_SUM,world);
  
  return sqrt(rg_all/masstotal);
}

/* ----------------------------------------------------------------------
   compute the angular momentum L (lmom) of group around center-of-mass cm
   must unwrap atoms to compute L correctly
------------------------------------------------------------------------- */

void Group::angmom(int igroup, double *cm, double *lmom)
{
  int groupbit = bitmask[igroup];

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

  int xbox,ybox,zbox;
  double dx,dy,dz,massone;
  double xprd = domain->xprd;
  double yprd = domain->yprd;
  double zprd = domain->zprd;
  double p[3];
  p[0] = p[1] = p[2] = 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;
      dx = (x[i][0] + xbox*xprd) - cm[0];
      dy = (x[i][1] + ybox*yprd) - cm[1];
      dz = (x[i][2] + zbox*zprd) - cm[2];
      if (mass) massone = mass[type[i]];
      else massone = rmass[i];
      p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
      p[1] += massone * (dz*v[i][0] - dx*v[i][2]);
      p[2] += massone * (dx*v[i][1] - dy*v[i][0]);
    }

  MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
}

/* ----------------------------------------------------------------------
   compute moment of inertia tensor around center-of-mass cm
   must unwrap atoms to compute itensor correctly
------------------------------------------------------------------------- */

void Group::inertia(int igroup, double *cm, double itensor[3][3])
{
  int i,j;

  int groupbit = bitmask[igroup];

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

  int xbox,ybox,zbox;
  double dx,dy,dz,massone;
  double xprd = domain->xprd;
  double yprd = domain->yprd;
  double zprd = domain->zprd;
  double ione[3][3];
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      ione[i][j] = 0.0;

  for (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;
      dx = (x[i][0] + xbox*xprd) - cm[0];
      dy = (x[i][1] + ybox*yprd) - cm[1];
      dz = (x[i][2] + zbox*zprd) - cm[2];
      if (mass) massone = mass[type[i]];
      else massone = rmass[i];
      ione[0][0] += massone * (dy*dy + dz*dz);
      ione[1][1] += massone * (dx*dx + dz*dz);
      ione[2][2] += massone * (dx*dx + dy*dy);
      ione[0][1] -= massone * dx*dy;
      ione[1][2] -= massone * dy*dz;
      ione[0][2] -= massone * dx*dz;
    }
  ione[1][0] = ione[0][1];
  ione[2][1] = ione[1][2];
  ione[2][0] = ione[0][2];

  MPI_Allreduce(&ione[0][0],&itensor[0][0],9,MPI_DOUBLE,MPI_SUM,world);
}

/* ----------------------------------------------------------------------
   compute angular velocity omega from L = Iw, inverting I to solve for w
   really not a group operation, but L and I were computed for a group
------------------------------------------------------------------------- */

void Group::omega(double *angmom, double inertia[3][3], double *w)
{
  double inverse[3][3];

  inverse[0][0] = inertia[1][1]*inertia[2][2] - inertia[1][2]*inertia[2][1];
  inverse[0][1] = -(inertia[0][1]*inertia[2][2] - inertia[0][2]*inertia[2][1]);
  inverse[0][2] = inertia[0][1]*inertia[1][2] - inertia[0][2]*inertia[1][1];

  inverse[1][0] = -(inertia[1][0]*inertia[2][2] - inertia[1][2]*inertia[2][0]);
  inverse[1][1] = inertia[0][0]*inertia[2][2] - inertia[0][2]*inertia[2][0];
  inverse[1][2] = -(inertia[0][0]*inertia[1][2] - inertia[0][2]*inertia[1][0]);

  inverse[2][0] = inertia[1][0]*inertia[2][1] - inertia[1][1]*inertia[2][0];
  inverse[2][1] = -(inertia[0][0]*inertia[2][1] - inertia[0][1]*inertia[2][0]);
  inverse[2][2] = inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0];

  double determinant = inertia[0][0]*inertia[1][1]*inertia[2][2] + 
    inertia[0][1]*inertia[1][2]*inertia[2][0] + 
    inertia[0][2]*inertia[1][0]*inertia[2][1] - 
    inertia[0][0]*inertia[1][2]*inertia[2][1] -
    inertia[0][1]*inertia[1][0]*inertia[2][2] - 
    inertia[2][0]*inertia[1][1]*inertia[0][2]; 

  int i,j;
  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      inverse[i][j] /= determinant;

  w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] + 
    inverse[0][2]*angmom[2];
  w[1] = inverse[1][0]*angmom[0] + inverse[1][1]*angmom[1] + 
    inverse[1][2]*angmom[2];
  w[2] = inverse[2][0]*angmom[0] + inverse[2][1]*angmom[1] + 
    inverse[2][2]*angmom[2];
}