Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86892091
group.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Oct 9, 06:05
Size
24 KB
Mime Type
text/x-c
Expires
Fri, Oct 11, 06:05 (2 d)
Engine
blob
Format
Raw Data
Handle
21501187
Attached To
rLAMMPS lammps
group.cpp
View Options
/* ----------------------------------------------------------------------
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 "force.h"
#include "region.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAX_GROUP 32
enum{TYPE,MOLECULE,ID};
enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN};
#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 = (char *) "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
// add to group if atom is in region
if (strcmp(arg[1],"region") == 0) {
if (narg != 3) error->all("Illegal group command");
int iregion = domain->find_region(arg[2]);
if (iregion == -1) error->all("Group region ID does not exist");
for (i = 0; i < nlocal; i++)
if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
mask[i] |= bit;
// style = logical condition
} else if (narg >= 3 &&
(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 = LT;
else if (strcmp(arg[2],"<=") == 0) condition = LE;
else if (strcmp(arg[2],">") == 0) condition = GT;
else if (strcmp(arg[2],">=") == 0) condition = GE;
else if (strcmp(arg[2],"==") == 0) condition = EQ;
else if (strcmp(arg[2],"!=") == 0) condition = NEQ;
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 == LT) {
for (i = 0; i < nlocal; i++) if (attribute[i] < bound1) mask[i] |= bit;
} else if (condition == LE) {
for (i = 0; i < nlocal; i++) if (attribute[i] <= bound1) mask[i] |= bit;
} else if (condition == GT) {
for (i = 0; i < nlocal; i++) if (attribute[i] > bound1) mask[i] |= bit;
} else if (condition == GE) {
for (i = 0; i < nlocal; i++) if (attribute[i] >= bound1) mask[i] |= bit;
} else if (condition == EQ) {
for (i = 0; i < nlocal; i++) if (attribute[i] == bound1) mask[i] |= bit;
} else if (condition == NEQ) {
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;
n = 0;
for (i = 0; i < nlocal; i++) if (mask[i] & bit) n++;
double rlocal = n;
double all;
MPI_Allreduce(&rlocal,&all,1,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
if (screen) fprintf(screen,"%.15g atoms in group %s\n",all,names[igroup]);
if (logfile)
fprintf(logfile,"%.15g 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(const 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 one = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) one += rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) one += mass[type[i]];
}
double all;
MPI_Allreduce(&one,&all,1,MPI_DOUBLE,MPI_SUM,world);
return all;
}
/* ----------------------------------------------------------------------
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);
minmax[0] = -minmax[0];
minmax[2] = -minmax[2];
minmax[4] = -minmax[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 (rmass) {
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;
}
} 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 = 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;
}
}
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 (rmass) {
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;
}
} else {
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;
}
}
MPI_Allreduce(p,cm,3,MPI_DOUBLE,MPI_SUM,world);
cm[0] /= masstotal;
cm[1] /= masstotal;
cm[2] /= masstotal;
}
/* ----------------------------------------------------------------------
compute the total force on group of atoms
------------------------------------------------------------------------- */
void Group::fcm(int igroup, double *cm)
{
int groupbit = bitmask[igroup];
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double flocal[3];
flocal[0] = flocal[1] = flocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
flocal[0] += f[i][0];
flocal[1] += f[i][1];
flocal[2] += f[i][2];
}
MPI_Allreduce(flocal,cm,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute the total kinetic energy of group and return it
------------------------------------------------------------------------- */
double Group::ke(int igroup)
{
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 one = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
}
double all;
MPI_Allreduce(&one,&all,1,MPI_DOUBLE,MPI_SUM,world);
all *= 0.5 * force->mvv2e;
return all;
}
/* ----------------------------------------------------------------------
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 (rmass) massone = rmass[i];
else massone = mass[type[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 (rmass) massone = rmass[i];
else massone = mass[type[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 (rmass) massone = rmass[i];
else massone = mass[type[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];
}
Event Timeline
Log In to Comment