From 6f4bba7edf244656cd1a552d8765c499646acc73 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Oct 2012 16:24:37 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9006 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_com_molecule.cpp | 16 +-- src/compute_gyration.cpp | 17 ++- src/compute_gyration_molecule.cpp | 57 ++++------ src/compute_msd_molecule.cpp | 22 ++-- src/fix_addforce.cpp | 15 +-- src/fix_momentum.cpp | 18 ++- src/fix_spring_rg.cpp | 15 +-- src/fix_spring_self.cpp | 34 ++---- src/fix_tmd.cpp | 51 +++------ src/group.cpp | 176 +++++++++++------------------- src/velocity.cpp | 15 +-- 11 files changed, 149 insertions(+), 287 deletions(-) diff --git a/src/compute_com_molecule.cpp b/src/compute_com_molecule.cpp index 80e745f056..82a212ca3a 100644 --- a/src/compute_com_molecule.cpp +++ b/src/compute_com_molecule.cpp @@ -96,8 +96,8 @@ void ComputeCOMMolecule::init() void ComputeCOMMolecule::compute_array() { int i,imol; - double xbox,ybox,zbox; double massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -113,23 +113,17 @@ void ComputeCOMMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/compute_gyration.cpp b/src/compute_gyration.cpp index 2e65cf1793..007d56f0c3 100644 --- a/src/compute_gyration.cpp +++ b/src/compute_gyration.cpp @@ -83,25 +83,22 @@ void ComputeGyration::compute_vector() 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 unwrap[3]; double rg[6]; rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; + + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; + rg[0] += dx*dx * massone; rg[1] += dy*dy * massone; rg[2] += dz*dz * massone; diff --git a/src/compute_gyration_molecule.cpp b/src/compute_gyration_molecule.cpp index 8fd9b6751e..1ebc47a9aa 100644 --- a/src/compute_gyration_molecule.cpp +++ b/src/compute_gyration_molecule.cpp @@ -123,8 +123,8 @@ void ComputeGyrationMolecule::init() void ComputeGyrationMolecule::compute_vector() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -141,21 +141,16 @@ void ComputeGyrationMolecule::compute_vector() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - dx = (x[i][0] + xbox*xprd) - comall[imol][0]; - dy = (x[i][1] + ybox*yprd) - comall[imol][1]; - dz = (x[i][2] + zbox*zprd) - comall[imol][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[imol][0]; + dy = unwrap[1] - comall[imol][1]; + dz = unwrap[2] - comall[imol][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg[imol] += (dx*dx + dy*dy + dz*dz) * massone; @@ -171,8 +166,8 @@ void ComputeGyrationMolecule::compute_vector() void ComputeGyrationMolecule::compute_array() { int i,j,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -191,21 +186,15 @@ void ComputeGyrationMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - dx = (x[i][0] + xbox*xprd) - comall[imol][0]; - dy = (x[i][1] + ybox*yprd) - comall[imol][1]; - dz = (x[i][2] + zbox*zprd) - comall[imol][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - comall[imol][0]; + dy = unwrap[1] - comall[imol][1]; + dz = unwrap[2] - comall[imol][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rgt[imol][0] += dx*dx * massone; @@ -233,8 +222,8 @@ void ComputeGyrationMolecule::compute_array() void ComputeGyrationMolecule::molcom() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; for (i = 0; i < nmolecules; i++) com[i][0] = com[i][1] = com[i][2] = 0.0; @@ -248,23 +237,17 @@ void ComputeGyrationMolecule::molcom() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/compute_msd_molecule.cpp b/src/compute_msd_molecule.cpp index c98cca8e1e..270cfad6f8 100644 --- a/src/compute_msd_molecule.cpp +++ b/src/compute_msd_molecule.cpp @@ -111,8 +111,8 @@ void ComputeMSDMolecule::init() void ComputeMSDMolecule::compute_array() { int i,imol; - double xbox,ybox,zbox,dx,dy,dz; - double massone; + double dx,dy,dz,massone; + double unwrap[3]; invoked_array = update->ntimestep; @@ -130,23 +130,17 @@ void ComputeMSDMolecule::compute_array() double *rmass = atom->rmass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; imol = molecule[i]; if (molmap) imol = molmap[imol-idlo]; else imol--; - com[imol][0] += (x[i][0] + xbox*xprd) * massone; - com[imol][1] += (x[i][1] + ybox*yprd) * massone; - com[imol][2] += (x[i][2] + zbox*zprd) * massone; + domain->unmap(x[i],image[i],unwrap); + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + com[imol][0] += unwrap[0] * massone; + com[imol][1] += unwrap[1] * massone; + com[imol][2] += unwrap[2] * massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp index 2523c7cc05..8ec546bb01 100644 --- a/src/fix_addforce.cpp +++ b/src/fix_addforce.cpp @@ -217,12 +217,6 @@ void FixAddForce::min_setup(int vflag) void FixAddForce::post_force(int vflag) { - int xbox,ybox,zbox; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double **x = atom->x; double **f = atom->f; int *mask = atom->mask; @@ -247,18 +241,15 @@ void FixAddForce::post_force(int vflag) // potential energy = - x dot f in unwrapped coords if (varflag == CONSTANT) { + double unwrap[3]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (iregion >= 0 && !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) continue; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - foriginal[0] -= xvalue * (x[i][0]+xbox*xprd) + - yvalue * (x[i][1]+ybox*yprd) + zvalue * (x[i][2]+zbox*zprd); - + domain->unmap(x[i],image[i],unwrap); + foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2]; foriginal[1] += f[i][0]; foriginal[2] += f[i][1]; foriginal[3] += f[i][2]; diff --git a/src/fix_momentum.cpp b/src/fix_momentum.cpp index 69e06aeaed..6f9da4427b 100644 --- a/src/fix_momentum.cpp +++ b/src/fix_momentum.cpp @@ -57,7 +57,8 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) : if (linear) if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 || - zflag < 0 || zflag > 1) error->all(FLERR,"Illegal fix momentum command"); + zflag < 0 || zflag > 1) + error->all(FLERR,"Illegal fix momentum command"); // cannot have 0 atoms in group @@ -121,20 +122,15 @@ void FixMomentum::end_of_step() tagint *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 unwrap[3]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][2] -= omega[0]*dy - omega[1]*dx; diff --git a/src/fix_spring_rg.cpp b/src/fix_spring_rg.cpp index e17aa8966c..0d93c4b628 100644 --- a/src/fix_spring_rg.cpp +++ b/src/fix_spring_rg.cpp @@ -110,20 +110,15 @@ void FixSpringRG::post_force(int vflag) int nlocal = atom->nlocal; double massfrac; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + double unwrap[3]; - int xbox,ybox,zbox; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; term1 = 2.0 * k * (1.0 - rg0/rg); - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; massfrac = mass[type[i]]/masstotal; f[i][0] -= term1*dx*massfrac; f[i][1] -= term1*dy*massfrac; diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index fb2120d4cb..043bd21e2e 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -45,9 +45,10 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command"); xflag = yflag = zflag = 1; + if (narg == 5) { if (strcmp(arg[4],"xyz") == 0) { - ; /* default */ + xflag = yflag = zflag = 1; } else if (strcmp(arg[4],"xy") == 0) { zflag = 0; } else if (strcmp(arg[4],"xz") == 0) { @@ -78,20 +79,9 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : tagint *image = atom->image; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xoriginal[i][0] = x[i][0] + xbox*xprd; - xoriginal[i][1] = x[i][1] + ybox*yprd; - xoriginal[i][2] = x[i][2] + zbox*zprd; - } else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; + if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]); + else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; } espring = 0.0; @@ -161,21 +151,17 @@ void FixSpringSelf::post_force(int vflag) tagint *image = atom->image; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - int xbox,ybox,zbox; double dx,dy,dz; + double unwrap[3]; + espring = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xoriginal[i][0]; - dy = x[i][1] + ybox*yprd - xoriginal[i][1]; - dz = x[i][2] + zbox*zprd - xoriginal[i][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xoriginal[i][0]; + dy = unwrap[1] - xoriginal[i][1]; + dz = unwrap[2] - xoriginal[i][2]; if (!xflag) dx = 0.0; if (!yflag) dy = 0.0; if (!zflag) dz = 0.0; diff --git a/src/fix_tmd.cpp b/src/fix_tmd.cpp index 71efd01bdf..fde131e632 100644 --- a/src/fix_tmd.cpp +++ b/src/fix_tmd.cpp @@ -98,26 +98,17 @@ FixTMD::FixTMD(LAMMPS *lmp, int narg, char **arg) : double *mass = atom->mass; int nlocal = atom->nlocal; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double dx,dy,dz; - int xbox,ybox,zbox; + rho_start = 0.0; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xf[i][0]; - dy = x[i][1] + ybox*yprd - xf[i][1]; - dz = x[i][2] + zbox*zprd - xf[i][2]; + domain->unmap(x[i],image[i],xold[i]); + dx = xold[i][0] - xf[i][0]; + dy = xold[i][1] - xf[i][1]; + dz = xold[i][2] - xf[i][2]; rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz); - xold[i][0] = x[i][0] + xbox*xprd; - xold[i][1] = x[i][1] + ybox*yprd; - xold[i][2] = x[i][2] + zbox*zprd; } else xold[i][0] = xold[i][1] = xold[i][2] = 0.0; } @@ -190,10 +181,8 @@ void FixTMD::initial_integrate(int vflag) double dxold,dyold,dzold,xback,yback,zback; double gamma_forward,gamma_back,gamma_max,lambda; double kt,fr,kttotal,frtotal,dtfm; + double unwrap[3]; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -202,7 +191,6 @@ void FixTMD::initial_integrate(int vflag) int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - int xbox,ybox,zbox; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -216,12 +204,10 @@ void FixTMD::initial_integrate(int vflag) dxold = xold[i][0] - xf[i][0]; dyold = xold[i][1] - xf[i][1]; dzold = xold[i][2] - xf[i][2]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - xf[i][0]; - dy = x[i][1] + ybox*yprd - xf[i][1]; - dz = x[i][2] + zbox*zprd - xf[i][2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xf[i][0]; + dy = unwrap[1] - xf[i][1]; + dz = unwrap[2] - xf[i][2]; a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold); b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold); e += mass[type[i]]*(dx *dx + dy *dy + dz *dz); @@ -258,12 +244,10 @@ void FixTMD::initial_integrate(int vflag) dxold = xold[i][0] - xf[i][0]; dyold = xold[i][1] - xf[i][1]; dzold = xold[i][2] - xf[i][2]; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xback = x[i][0] + xbox*xprd + gamma_back*dxold; - yback = x[i][1] + ybox*yprd + gamma_back*dyold; - zback = x[i][2] + zbox*zprd + gamma_back*dzold; + domain->unmap(x[i],image[i],unwrap); + xback = unwrap[0] + gamma_back*dxold; + yback = unwrap[1] + gamma_back*dyold; + zback = unwrap[2] + gamma_back*dzold; dxkt = xback - xold[i][0]; dykt = yback - xold[i][1]; dzkt = zback - xold[i][2]; @@ -314,12 +298,7 @@ void FixTMD::initial_integrate(int vflag) x[i][2] += gamma_forward*dzold; v[i][2] += gamma_forward*dzold/dtv; f[i][2] += gamma_forward*dzold/dtv/dtfm; - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - xold[i][0] = x[i][0] + xbox*xprd; - xold[i][1] = x[i][1] + ybox*yprd; - xold[i][2] = x[i][2] + zbox*zprd; + domain->unmap(x[i],image[i],xold[i]); } } } diff --git a/src/group.cpp b/src/group.cpp index 21bea8f750..c82d110058 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -766,34 +766,27 @@ void Group::xcm(int igroup, double masstotal, double *cm) 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; + double unwrap[3]; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; 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; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; 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; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } @@ -827,34 +820,27 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion) 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; + double unwrap[3]; 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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; 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; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[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])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; 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; + domain->unmap(x[i],image[i],unwrap); + cmone[0] += unwrap[0] * massone; + cmone[1] += unwrap[1] * massone; + cmone[2] += unwrap[2] * massone; } } @@ -1102,21 +1088,16 @@ double Group::gyration(int igroup, double masstotal, double *cm) 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 unwrap[3]; double rg = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg += (dx*dx + dy*dy + dz*dz) * massone; @@ -1147,21 +1128,16 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion) 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 unwrap[3]; 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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; rg += (dx*dx + dy*dy + dz*dz) * massone; @@ -1192,22 +1168,18 @@ void Group::angmom(int igroup, double *cm, double *lmom) 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 unwrap[3]; + 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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; p[0] += massone * (dy*v[i][2] - dz*v[i][1]); @@ -1238,22 +1210,18 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion) 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 unwrap[3]; + double p[3]; p[0] = p[1] = p[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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; p[0] += massone * (dy*v[i][2] - dz*v[i][1]); @@ -1280,22 +1248,18 @@ void Group::torque(int igroup, double *cm, double *tq) tagint *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 unwrap[3]; + 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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - 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]; @@ -1321,22 +1285,18 @@ void Group::torque(int igroup, double *cm, double *tq, int iregion) tagint *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 unwrap[3]; + 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] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - 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]; @@ -1364,11 +1324,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) 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 unwrap[3]; + double ione[3][3]; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) @@ -1376,12 +1334,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; ione[0][0] += massone * (dy*dy + dz*dz); @@ -1418,11 +1374,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) 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 unwrap[3]; + double ione[3][3]; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) @@ -1430,12 +1384,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - cm[0]; - dy = (x[i][1] + ybox*yprd) - cm[1]; - dz = (x[i][2] + zbox*zprd) - cm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - cm[0]; + dy = unwrap[1] - cm[1]; + dz = unwrap[2] - cm[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; ione[0][0] += massone * (dy*dy + dz*dz); diff --git a/src/velocity.cpp b/src/velocity.cpp index 6063eb46fc..e69b6d67a5 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -718,20 +718,15 @@ void Velocity::zero_rotation() tagint *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 unwrap[3]; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = (x[i][0] + xbox*xprd) - xcm[0]; - dy = (x[i][1] + ybox*yprd) - xcm[1]; - dz = (x[i][2] + zbox*zprd) - xcm[2]; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0]- xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][2] -= omega[0]*dy - omega[1]*dx;