diff --git a/src/domain.cpp b/src/domain.cpp index cece653486..c2b700ea2e 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -678,7 +678,7 @@ void Domain::remap(double *x, int &image) /* ---------------------------------------------------------------------- unmap the point via image flags - don't reset image flag + x overwritten with result, don't reset image flag for triclinic, use h[] to add in tilt factors in other dims as needed ------------------------------------------------------------------------- */ @@ -699,6 +699,29 @@ void Domain::unmap(double *x, int image) } } +/* ---------------------------------------------------------------------- + unmap the point via image flags + result returned in y, don't reset image flag + for triclinic, use h[] to add in tilt factors in other dims as needed +------------------------------------------------------------------------- */ + +void Domain::unmap(double *x, int image, double *y) +{ + int xbox = (image & 1023) - 512; + int ybox = (image >> 10 & 1023) - 512; + int zbox = (image >> 20) - 512; + + if (triclinic == 0) { + y[0] = x[0] + xbox*xprd; + y[1] = x[1] + ybox*yprd; + y[2] = x[2] + zbox*zprd; + } else { + y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; + y[1] = x[1] + h[1]*ybox + h[3]*zbox; + y[2] = x[2] + h[2]*zbox; + } +} + /* ---------------------------------------------------------------------- create a lattice delete it if style = none diff --git a/src/domain.h b/src/domain.h index a0bd61eec5..318a01fe06 100644 --- a/src/domain.h +++ b/src/domain.h @@ -91,6 +91,7 @@ class Domain : protected Pointers { void pbc(); void remap(double *, int &); void unmap(double *, int); + void unmap(double *, int, double *); void minimum_image(double &, double &, double &); void minimum_image(double *); void set_lattice(int, char **); diff --git a/src/fix_msd.cpp b/src/fix_msd.cpp index 90b74d8d06..dd2b92a901 100644 --- a/src/fix_msd.cpp +++ b/src/fix_msd.cpp @@ -72,14 +72,8 @@ FixMSD::FixMSD(LAMMPS *lmp, int narg, char **arg) : int xbox,ybox,zbox; 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; - 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; } // nmsd = # of atoms in group @@ -149,6 +143,7 @@ void FixMSD::end_of_step() int *image = atom->image; int nlocal = atom->nlocal; + double *h = domain->h; double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; @@ -157,19 +152,39 @@ void FixMSD::end_of_step() double msd[4]; msd[0] = msd[1] = msd[2] = msd[3] = 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 - xoriginal[i][0]; - dy = x[i][1] + ybox*yprd - xoriginal[i][1]; - dz = x[i][2] + zbox*zprd - xoriginal[i][2]; - msd[0] += dx*dx; - msd[1] += dy*dy; - msd[2] += dz*dz; - msd[3] += dx*dx + dy*dy + dz*dz; - } + // dx,dy,dz = displacement of atom from original position + // for triclinic, need to unwrap current atom coord via h matrix + + if (domain->triclinic == 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 - xoriginal[i][0]; + dy = x[i][1] + ybox*yprd - xoriginal[i][1]; + dz = x[i][2] + zbox*zprd - xoriginal[i][2]; + msd[0] += dx*dx; + msd[1] += dy*dy; + msd[2] += dz*dz; + msd[3] += dx*dx + dy*dy + dz*dz; + } + + } else { + 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] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0]; + dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1]; + dz = x[i][2] + h[2]*zbox - xoriginal[i][2]; + msd[0] += dx*dx; + msd[1] += dy*dy; + msd[2] += dz*dz; + msd[3] += dx*dx + dy*dy + dz*dz; + } + } double msd_all[4]; MPI_Allreduce(msd,msd_all,4,MPI_DOUBLE,MPI_SUM,world);