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

This commit is contained in:
sjplimp 2012-10-25 16:28:17 +00:00
parent 43b161eb33
commit e9fa893b1c
2 changed files with 30 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/USER-MISC/compute_temp_rotate.cpp
View File

@ -89,11 +89,8 @@ double ComputeTempRotate::compute_scalar()
{ {
double vthermal[3]; double vthermal[3];
double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3];
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
invoked_scalar = update->ntimestep; invoked_scalar = update->ntimestep;
@ -123,17 +120,13 @@ double ComputeTempRotate::compute_scalar()
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],unwrap)
xbox = (image[i] & IMGMASK) - IMGMAX; dx = unwrap[0] - xcm[0];
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dy = unwrap[1] - xcm[1];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dz = unwrap[2] - xcm[2];
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2];
vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0];
vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1];
vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[0] = v[i][0] - vbiasall[i][0];
vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[1] = v[i][1] - vbiasall[i][1];
vthermal[2] = v[i][2] - vbiasall[i][2]; vthermal[2] = v[i][2] - vbiasall[i][2];
@ -155,14 +148,10 @@ double ComputeTempRotate::compute_scalar()
void ComputeTempRotate::compute_vector() void ComputeTempRotate::compute_vector()
{ {
int i;
double vthermal[3]; double vthermal[3];
double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3];
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
invoked_vector = update->ntimestep; invoked_vector = update->ntimestep;
@ -189,25 +178,20 @@ void ComputeTempRotate::compute_vector()
} }
double massone,t[6]; double massone,t[6];
for (i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < 6; i++) t[i] = 0.0;
for (i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],unwrap)
xbox = (image[i] & IMGMASK) - IMGMAX; dx = unwrap[0] - xcm[0];
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dy = unwrap[1] - xcm[1];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dz = unwrap[2] - xcm[2];
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2];
vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0];
vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1];
vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[0] = v[i][0] - vbiasall[i][0];
vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[1] = v[i][1] - vbiasall[i][1];
vthermal[2] = v[i][2] - vbiasall[i][2]; vthermal[2] = v[i][2] - vbiasall[i][2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0]; t[0] += massone * vthermal[0]*vthermal[0];
@ -219,7 +203,7 @@ void ComputeTempRotate::compute_vector()
} }
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

37
src/USER-MISC/fix_addtorque.cpp
View File

@ -109,19 +109,22 @@ void FixAddTorque::init()
if (xstr) { if (xstr) {
xvar = input->variable->find(xstr); xvar = input->variable->find(xstr);
if (xvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); if (xvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL; if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style"); else error->all(FLERR,"Variable for fix addtorque is invalid style");
} }
if (ystr) { if (ystr) {
yvar = input->variable->find(ystr); yvar = input->variable->find(ystr);
if (yvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); if (yvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL; if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style"); else error->all(FLERR,"Variable for fix addtorque is invalid style");
} }
if (zstr) { if (zstr) {
zvar = input->variable->find(zstr); zvar = input->variable->find(zstr);
if (zvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist"); if (zvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL; if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style"); else error->all(FLERR,"Variable for fix addtorque is invalid style");
} }
@ -168,16 +171,14 @@ void FixAddTorque::post_force(int vflag)
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double mvv2e = force->mvv2e; double mvv2e = force->mvv2e;
int xbox,ybox,zbox;
double dx,dy,dz,vx,vy,vz,fx,fy,fz,massone,omegadotr; double dx,dy,dz,vx,vy,vz,fx,fy,fz,massone,omegadotr;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tcm[3],xcm[3],angmom[3],omega[3],itorque[3],domegadt[3],tlocal[3]; double tcm[3],xcm[3],angmom[3],omega[3],itorque[3],domegadt[3],tlocal[3];
double inertia[3][3]; double inertia[3][3];
double unwrap[3];
// foriginal[0] = "potential energy" for added force // foriginal[0] = "potential energy" for added force
// foriginal[123] = torque on atoms before extra force added // foriginal[123] = torque on atoms before extra force added
foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0; foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
force_flag = 0; force_flag = 0;
@ -200,12 +201,10 @@ void FixAddTorque::post_force(int vflag)
tlocal[0] = tlocal[1] = tlocal[2] = 0.0; tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap)
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - xcm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xcm[1];
dx = (x[i][0] + xbox*xprd) - xcm[0]; dz = unwrap[2] - xcm[2];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
omegadotr = omega[0]*dx+omega[1]*dy+omega[2]*dz; omegadotr = omega[0]*dx+omega[1]*dy+omega[2]*dz;
@ -222,12 +221,10 @@ void FixAddTorque::post_force(int vflag)
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap)
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - xcm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xcm[1];
dx = (x[i][0] + xbox*xprd) - xcm[0]; dz = unwrap[2] - xcm[2];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
vx = mvv2e*(dz*omega[1]-dy*omega[2]); vx = mvv2e*(dz*omega[1]-dy*omega[2]);
vy = mvv2e*(dx*omega[2]-dz*omega[0]); vy = mvv2e*(dx*omega[2]-dz*omega[0]);
vz = mvv2e*(dy*omega[0]-dx*omega[1]); vz = mvv2e*(dy*omega[0]-dx*omega[1]);