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@5302 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2010-11-24 14:39:50 +00:00
parent 0d9d109d49
commit 2de8e4d51d
3 changed files with 105 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
src/group.cpp
View File

@ -1240,6 +1240,87 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute the torque T (tq) on group
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void Group::torque(int igroup, double *cm, double *tq)
{
int groupbit = bitmask[igroup];
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[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];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute the torque T (tq) on group of atoms in region
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void Group::torque(int igroup, double *cm, double *tq, int iregion)
{
int groupbit = bitmask[igroup];
Region *region = domain->regions[iregion];
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[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])) {
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];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute moment of inertia tensor around center-of-mass cm of group
must unwrap atoms to compute itensor correctly

2
src/group.h
View File

@ -54,6 +54,8 @@ class Group : protected Pointers {
double gyration(int, double, double *, int);
void angmom(int, double *, double *); // angular momentum of group
void angmom(int, double *, double *, int);
void torque(int, double *, double *); // torque on group
void torque(int, double *, double *, int);
void inertia(int, double *, double [3][3]); // inertia tensor
void inertia(int, double *, double [3][3], int);
void omega(double *, double [3][3], double *); // angular velocity

25
src/variable.cpp
View File

@ -2423,8 +2423,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
customize by adding a group function with optional region arg:
count(group),mass(group),charge(group),
xcm(group,dim),vcm(group,dim),fcm(group,dim),
bound(group,xmin),gyration(group),ke(group),angmom(group),
inertia(group,dim),omega(group,dim)
bound(group,xmin),gyration(group),ke(group),angmom(group,dim),
torque(group,dim),inertia(group,dim),omega(group,dim)
------------------------------------------------------------------------- */
int Variable::group_function(char *word, char *contents, Tree **tree,
@ -2438,7 +2438,8 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
strcmp(word,"vcm") && strcmp(word,"fcm") &&
strcmp(word,"bound") && strcmp(word,"gyration") &&
strcmp(word,"ke") && strcmp(word,"angmom") &&
strcmp(word,"inertia") && strcmp(word,"omega"))
strcmp(word,"torque") && strcmp(word,"inertia") &&
strcmp(word,"omega"))
return 0;
// parse contents for arg1,arg2,arg3 separated by commas
@ -2588,6 +2589,24 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
else if (strcmp(arg2,"z") == 0) value = lmom[2];
else error->all("Invalid group function in variable formula");
} else if (strcmp(word,"torque") == 0) {
atom->check_mass();
double xcm[3],tq[3];
if (narg == 2) {
double masstotal = group->mass(igroup);
group->xcm(igroup,masstotal,xcm);
group->torque(igroup,xcm,tq);
} else if (narg == 3) {
int iregion = region_function(arg3);
double masstotal = group->mass(igroup,iregion);
group->xcm(igroup,masstotal,xcm,iregion);
group->torque(igroup,xcm,tq,iregion);
} else error->all("Invalid group function in variable formula");
if (strcmp(arg2,"x") == 0) value = tq[0];
else if (strcmp(arg2,"y") == 0) value = tq[1];
else if (strcmp(arg2,"z") == 0) value = tq[2];
else error->all("Invalid group function in variable formula");
} else if (strcmp(word,"inertia") == 0) {
atom->check_mass();
double xcm[3],inertia[3][3];