jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2009-04-28 19:46:52 +00:00
parent 2dacebfcd0
commit 1904459e45
4 changed files with 447 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

356
src/group.cpp
View File

@ -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

9
src/group.h
View File

@ -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

125
src/variable.cpp
View File

@ -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

1
src/variable.h
View File

@ -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 &);