From 1904459e45f7f54e52debdcbab4dcbcf6d418d0f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 28 Apr 2009 19:46:52 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2776 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/group.cpp | 356 ++++++++++++++++++++++++++++++++++++++++++++++- src/group.h | 9 ++ src/variable.cpp | 125 ++++++++++++----- src/variable.h | 1 + 4 files changed, 447 insertions(+), 44 deletions(-) diff --git a/src/group.cpp b/src/group.cpp index 1cdf092048..d75de55e9d 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -501,6 +501,30 @@ double Group::count(int igroup) return nall; } +/* ---------------------------------------------------------------------- + count atoms in group and region + compute in double precision in case system is huge +------------------------------------------------------------------------- */ + +double Group::count(int igroup, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int n = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) 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 @@ -510,10 +534,10 @@ 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 *mask = atom->mask; + int *type = atom->type; int nlocal = atom->nlocal; double one = 0.0; @@ -531,6 +555,40 @@ double Group::mass(int igroup) return all; } +/* ---------------------------------------------------------------------- + compute the total mass of group of atoms in region + use either per-type mass or per-atom rmass +------------------------------------------------------------------------- */ + +double Group::mass(int igroup, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + double one = 0.0; + + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) + one += rmass[i]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) + 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 ------------------------------------------------------------------------- */ @@ -539,8 +597,8 @@ double Group::charge(int igroup) { int groupbit = bitmask[igroup]; - int *mask = atom->mask; double *q = atom->q; + int *mask = atom->mask; int nlocal = atom->nlocal; double qone = 0.0; @@ -552,6 +610,30 @@ double Group::charge(int igroup) return qall; } +/* ---------------------------------------------------------------------- + compute the total charge of group of atoms in region +------------------------------------------------------------------------- */ + +double Group::charge(int igroup, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double qone = 0.0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) + 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 @@ -596,7 +678,52 @@ void Group::bounds(int igroup, double *minmax) } /* ---------------------------------------------------------------------- - compute the center-of-mass coords of group with total mass = masstotal + compute the coordinate bounds of the group of atoms in region + periodic images are not considered, so atoms are NOT unwrapped +------------------------------------------------------------------------- */ + +void Group::bounds(int igroup, double *minmax, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 of atoms + masstotal = total mass return center-of-mass coords in cm[] must unwrap atoms to compute center-of-mass correctly ------------------------------------------------------------------------- */ @@ -653,7 +780,67 @@ void Group::xcm(int igroup, double masstotal, double *cm) } /* ---------------------------------------------------------------------- - compute the center-of-mass velocity of group with total mass = masstotal + compute the center-of-mass coords of group of atoms in region + mastotal = total mass + 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 iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 of atoms + masstotal = total mass return center-of-mass velocity in cm[] ------------------------------------------------------------------------- */ @@ -695,6 +882,52 @@ void Group::vcm(int igroup, double masstotal, double *cm) cm[2] /= masstotal; } +/* ---------------------------------------------------------------------- + compute the center-of-mass velocity of group of atoms in region + masstotal = total mass + return center-of-mass velocity in cm[] +------------------------------------------------------------------------- */ + +void Group::vcm(int igroup, double masstotal, double *cm, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 ------------------------------------------------------------------------- */ @@ -721,7 +954,34 @@ void Group::fcm(int igroup, double *cm) } /* ---------------------------------------------------------------------- - compute the total kinetic energy of group and return it + compute the total force on group of atoms in region +------------------------------------------------------------------------- */ + +void Group::fcm(int igroup, double *cm, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 of atoms and return it ------------------------------------------------------------------------- */ double Group::ke(int igroup) @@ -756,7 +1016,45 @@ double Group::ke(int igroup) } /* ---------------------------------------------------------------------- - compute the radius-of-gyration of group around center-of-mass cm + compute the total kinetic energy of group of atoms in region and return it +------------------------------------------------------------------------- */ + +double Group::ke(int igroup, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + double **x = atom->x; + 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 && region->match(x[i][0],x[i][1],x[i][2])) + 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 && region->match(x[i][0],x[i][1],x[i][2])) + 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 of atoms + around center-of-mass cm must unwrap atoms to compute Rg correctly ------------------------------------------------------------------------- */ @@ -797,6 +1095,50 @@ double Group::gyration(int igroup, double masstotal, double *cm) return sqrt(rg_all/masstotal); } +/* ---------------------------------------------------------------------- + compute the radius-of-gyration of group of atoms in region + around center-of-mass cm + must unwrap atoms to compute Rg correctly +------------------------------------------------------------------------- */ + +double Group::gyration(int igroup, double masstotal, double *cm, int iregion) +{ + int groupbit = bitmask[igroup]; + Region *region = domain->regions[iregion]; + + 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 && region->match(x[i][0],x[i][1],x[i][2])) { + 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 diff --git a/src/group.h b/src/group.h index 74565a3d68..0e4e5a6d47 100644 --- a/src/group.h +++ b/src/group.h @@ -33,14 +33,23 @@ class Group : protected Pointers { void read_restart(FILE *); double count(int); // count atoms in group + double count(int,int); // count atoms in group & region double mass(int); // total mass of atoms in group + double mass(int,int); double charge(int); // total charge of atoms in group + double charge(int,int); void bounds(int, double *); // bounds of atoms in group + void bounds(int, double *, int); void xcm(int, double, double *); // center-of-mass coords of group + void xcm(int, double, double *, int); void vcm(int, double, double *); // center-of-mass velocity of group + void vcm(int, double, double *, int); void fcm(int, double *); // total force on group + void fcm(int, double *, int); double ke(int); // kinetic energy of group + double ke(int, int); double gyration(int, double, double *); // radius-of-gyration of group + double gyration(int, double, double *, int); void angmom(int, double *, double *); // angular momentum of group void inertia(int, double *, double [3][3]); // inertia tensor void omega(double *, double [3][3], double *); // angular velocity diff --git a/src/variable.cpp b/src/variable.cpp index 39f2e32fa7..231d3de53a 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1408,35 +1408,46 @@ int Variable::math_function(char *word, char *contents, Tree **tree, } /* ---------------------------------------------------------------------- - process a group function in formula + process a group function in formula with optional region arg push result onto tree or arg stack word = group function - contents = str bewteen parentheses with one or two args + contents = str bewteen parentheses with one,two,three args return 0 if not a match, 1 if successfully processed customize by adding a group function: count(group),mass(group),charge(group), xcm(group,dim),vcm(group,dim),fcm(group,dim), - bound(group,xmin),gyration(group) + bound(group,xmin),gyration(group),ke(group) ------------------------------------------------------------------------- */ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **treestack, int &ntreestack, double *argstack, int &nargstack) { - // parse contents for arg1,arg2 separated by comma + // parse contents for arg1,arg2,arg3 separated by commas - char *arg1,*arg2; - char *ptr = strchr(contents,','); - if (ptr) *ptr = '\0'; + char *ptr1 = strchr(contents,','); + if (ptr1) *ptr1 = '\0'; + char *ptr2 = strchr(ptr1+1,','); + if (ptr2) *ptr2 = '\0'; + + char *arg1,*arg2,*arg3; int n = strlen(contents) + 1; arg1 = new char[n]; strcpy(arg1,contents); - if (ptr) { - n = strlen(ptr+1) + 1; + int narg = 1; + if (ptr1) { + n = strlen(ptr1+1) + 1; arg2 = new char[n]; - strcpy(arg2,ptr+1); + strcpy(arg2,ptr1+1); + narg = 2; } else arg2 = NULL; - + if (ptr2) { + n = strlen(ptr2+1) + 1; + arg3 = new char[n]; + strcpy(arg3,ptr1+1); + narg = 3; + } else arg3 = NULL; + int igroup = group->find(arg1); if (igroup == -1) error->all("Group ID in variable formula does not exist"); @@ -1454,59 +1465,67 @@ int Variable::group_function(char *word, char *contents, Tree **tree, // match word to group function if (strcmp(word,"count") == 0) { - if (arg2) - error->all("Invalid group function in variable formula"); - value = group->count(igroup); + if (narg == 1) value = group->count(igroup); + else if (narg == 2) value = group->count(igroup,region_function(arg2)); + else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"mass") == 0) { - if (arg2) - error->all("Invalid group function in variable formula"); - value = group->mass(igroup); + if (narg == 1) value = group->mass(igroup); + else if (narg == 2) value = group->mass(igroup,region_function(arg2)); + else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"charge") == 0) { - if (arg2) - error->all("Invalid group function in variable formula"); - value = group->charge(igroup); + if (narg == 1) value = group->charge(igroup); + else if (narg == 2) value = group->charge(igroup,region_function(arg2)); + else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"xcm") == 0) { - if (!arg2) - error->all("Invalid group function in variable formula"); atom->check_mass(); - double masstotal = group->mass(igroup); double xcm[3]; - group->xcm(igroup,masstotal,xcm); + if (narg == 2) { + double masstotal = group->mass(igroup); + group->xcm(igroup,masstotal,xcm); + } else if (narg == 3) { + int iregion = region_function(arg3); + double masstotal = group->mass(igroup,iregion); + group->xcm(igroup,masstotal,xcm,iregion); + } else error->all("Invalid group function in variable formula"); if (strcmp(arg2,"x") == 0) value = xcm[0]; else if (strcmp(arg2,"y") == 0) value = xcm[1]; else if (strcmp(arg2,"z") == 0) value = xcm[2]; else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"vcm") == 0) { - if (!arg2) - error->all("Invalid group function in variable formula"); atom->check_mass(); - double masstotal = group->mass(igroup); double vcm[3]; - group->vcm(igroup,masstotal,vcm); + if (narg == 2) { + double masstotal = group->mass(igroup); + group->vcm(igroup,masstotal,vcm); + } else if (narg == 3) { + int iregion = region_function(arg3); + double masstotal = group->mass(igroup,iregion); + group->vcm(igroup,masstotal,vcm,iregion); + } else error->all("Invalid group function in variable formula"); if (strcmp(arg2,"x") == 0) value = vcm[0]; else if (strcmp(arg2,"y") == 0) value = vcm[1]; else if (strcmp(arg2,"z") == 0) value = vcm[2]; else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"fcm") == 0) { - if (!arg2) - error->all("Invalid group function in variable formula"); double fcm[3]; - group->fcm(igroup,fcm); + if (narg == 2) group->fcm(igroup,fcm); + else if (narg == 3) group->fcm(igroup,fcm,region_function(arg3)); + else error->all("Invalid group function in variable formula"); if (strcmp(arg2,"x") == 0) value = fcm[0]; else if (strcmp(arg2,"y") == 0) value = fcm[1]; else if (strcmp(arg2,"z") == 0) value = fcm[2]; else error->all("Invalid group function in variable formula"); } else if (strcmp(word,"bound") == 0) { - if (!arg2) - error->all("Invalid group function in variable formula"); double minmax[6]; - group->bounds(igroup,minmax); + if (narg == 2) group->bounds(igroup,minmax); + else if (narg == 3) group->bounds(igroup,minmax,region_function(arg3)); + else error->all("Invalid group function in variable formula"); if (strcmp(arg2,"xmin") == 0) value = minmax[0]; else if (strcmp(arg2,"xmax") == 0) value = minmax[1]; else if (strcmp(arg2,"ymin") == 0) value = minmax[2]; @@ -1517,10 +1536,22 @@ int Variable::group_function(char *word, char *contents, Tree **tree, } else if (strcmp(word,"gyration") == 0) { atom->check_mass(); - double masstotal = group->mass(igroup); double xcm[3]; - group->xcm(igroup,masstotal,xcm); - value = group->gyration(igroup,masstotal,xcm); + if (narg == 1) { + double masstotal = group->mass(igroup); + group->xcm(igroup,masstotal,xcm); + value = group->gyration(igroup,masstotal,xcm); + } else if (narg == 2) { + int iregion = region_function(arg2); + double masstotal = group->mass(igroup,iregion); + group->xcm(igroup,masstotal,xcm,iregion); + value = group->gyration(igroup,masstotal,xcm,iregion); + } else error->all("Invalid group function in variable formula"); + + } else if (strcmp(word,"ke") == 0) { + if (narg == 1) value = group->ke(igroup); + else if (narg == 2) value = group->ke(igroup,region_function(arg2)); + else error->all("Invalid group function in variable formula"); } else return 0; @@ -1533,6 +1564,26 @@ int Variable::group_function(char *word, char *contents, Tree **tree, return 1; } +/* ---------------------------------------------------------------------- + process a group function in formula with optional region arg + push result onto tree or arg stack + word = group function + contents = str bewteen parentheses with one,two,three args + return 0 if not a match, 1 if successfully processed + customize by adding a group function: + count(group),mass(group),charge(group), + xcm(group,dim),vcm(group,dim),fcm(group,dim), + bound(group,xmin),gyration(group) +------------------------------------------------------------------------- */ + +int Variable::region_function(char *id) +{ + int iregion = domain->find_region(id); + if (iregion == -1) + error->all("Invalid region in group function in variable formula"); + return iregion; +} + /* ---------------------------------------------------------------------- extract a global value from a per-atom quantity in a formula flag = 0 -> word is an atom vector diff --git a/src/variable.h b/src/variable.h index 2c60f80430..60da3aeb40 100644 --- a/src/variable.h +++ b/src/variable.h @@ -61,6 +61,7 @@ class Variable : protected Pointers { int int_between_brackets(char *, int, int &, int); int math_function(char *, char *, Tree **, Tree **, int &, double *, int &); int group_function(char *, char *, Tree **, Tree **, int &, double *, int &); + int region_function(char *); void peratom2global(int, char *, double *, int, int, Tree **, Tree **, int &, double *, int &); void atom_vector(char *, Tree **, Tree **, int &);