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

src/compute_chunk_atom.cpp

@@ -191,6 +191,7 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
   maxflag[1] = UPPER;
   maxflag[2] = UPPER;
   scaleflag = LATTICE;
+  pbcflag = 0;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"region") == 0) {
@@ -269,6 +270,12 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED;
       else error->all(FLERR,"Illegal compute chunk/atom command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"pbc") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
+      if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1;
+      else error->all(FLERR,"Illegal compute chunk/atom command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal compute chunk/atom command");
   }
 
@@ -1265,7 +1272,7 @@ int ComputeChunkAtom::setup_xyz_bins()
 int ComputeChunkAtom::setup_sphere_bins()
 {
   // convert sorigin_user to sorigin
-  // sorigin is always in box units, for orthogonal or triclinic domains
+  // sorigin,srad are always in box units, for orthogonal or triclinic domains
   // lamda2x works for either orthogonal or triclinic
 
   if (scaleflag == REDUCED) {
@@ -1280,6 +1287,23 @@ int ComputeChunkAtom::setup_sphere_bins()
     sradmax = sradmax_user;
   }
 
+  // if pbcflag set, sradmax must be < 1/2 box in any periodic dim
+  // treat orthongonal and triclinic the same
+  // check every time bins are created
+
+  if (pbcflag) {
+    double *prd_half = domain->prd_half;
+    int *periodicity = domain->periodicity;
+    int flag = 0;
+    if (periodicity[0] && sradmax > prd_half[0]) flag = 1;
+    if (periodicity[1] && sradmax > prd_half[1]) flag = 1;
+    if (domain->dimension == 3 && 
+        periodicity[2] && sradmax > prd_half[2]) flag = 1;
+    if (flag) 
+      error->all(FLERR,"Compute chunk/atom bin/sphere radius "
+                 "is too large for periodic box");
+  }
+
   sinvrad = nsbin / (sradmax-sradmin);
 
   // allocate and set bin coordinates
@@ -1328,6 +1352,21 @@ int ComputeChunkAtom::setup_cylinder_bins()
     cradmax = cradmax_user;
   }
 
+  // if pbcflag set, sradmax must be < 1/2 box in any periodic non-axis dim
+  // treat orthongonal and triclinic the same
+  // check every time bins are created
+
+  if (pbcflag) {
+    double *prd_half = domain->prd_half;
+    int *periodicity = domain->periodicity;
+    int flag = 0;
+    if (periodicity[cdim1] && sradmax > prd_half[cdim1]) flag = 1;
+    if (periodicity[cdim2] && sradmax > prd_half[cdim2]) flag = 1;
+    if (flag) 
+      error->all(FLERR,"Compute chunk/atom bin/cylinder radius "
+                 "is too large for periodic box");
+  }
+
   cinvrad = ncbin / (cradmax-cradmin);
 
   // allocate and set radial bin coordinates
@@ -1743,6 +1782,7 @@ void ComputeChunkAtom::atom2binsphere()
   double *boxlo = domain->boxlo;
   double *boxhi = domain->boxhi;
   double *prd = domain->prd;
+  double *prd_half = domain->prd_half;
   int *periodicity = domain->periodicity;
 
   // remap each atom's relevant coords back into box via PBC if necessary
@@ -1773,6 +1813,33 @@ void ComputeChunkAtom::atom2binsphere()
     dx = xremap - sorigin[0];
     dy = yremap - sorigin[1];
     dz = zremap - sorigin[2];
+
+    // if requested, apply PBC to distance from sphere center
+    // treat orthogonal and triclinic the same
+    //   with dx,dy,dz = lengths independent of each other
+    // so do not use domain->minimum_image() which couples for triclinic
+
+    if (pbcflag) {
+      if (periodicity[0]) {
+        if (fabs(dx) > prd_half[0]) {
+          if (dx < 0.0) dx += prd[0];
+          else dx -= prd[0];
+        }
+      }
+      if (periodicity[1]) {
+        if (fabs(dy) > prd_half[1]) {
+          if (dy < 0.0) dy += prd[1];
+          else dy -= prd[1];
+        }
+      }
+      if (periodicity[2]) {
+        if (fabs(dz) > prd_half[2]) {
+          if (dz < 0.0) dz += prd[2];
+          else dz -= prd[2];
+        }
+      }
+    }
+
     r = sqrt(dx*dx + dy*dy + dz*dz);
 
     ibin = static_cast<int> ((r - sradmin) * sinvrad);
@@ -1792,7 +1859,6 @@ void ComputeChunkAtom::atom2binsphere()
 
 /* ----------------------------------------------------------------------
    assign each atom to a cylindrical bin
-   
 ------------------------------------------------------------------------- */
 
 void ComputeChunkAtom::atom2bincylinder()
@@ -1812,6 +1878,7 @@ void ComputeChunkAtom::atom2bincylinder()
   double *boxlo = domain->boxlo;
   double *boxhi = domain->boxhi;
   double *prd = domain->prd;
+  double *prd_half = domain->prd_half;
   int *periodicity = domain->periodicity;
 
   // remap each atom's relevant coords back into box via PBC if necessary
@@ -1837,6 +1904,26 @@ void ComputeChunkAtom::atom2bincylinder()
 
     d1 = remap1 - corigin[cdim1];
     d2 = remap2 - corigin[cdim2];
+
+    // if requested, apply PBC to distance from cylinder axis
+    // treat orthogonal and triclinic the same
+    //   with d1,d2 = lengths independent of each other
+
+    if (pbcflag) {
+      if (periodicity[cdim1]) {
+        if (fabs(d1) > prd_half[cdim1]) {
+          if (d1 < 0.0) d1 += prd[cdim1];
+          else d1 -= prd[cdim1];
+        }
+      }
+      if (periodicity[cdim2]) {
+        if (fabs(d2) > prd_half[cdim2]) {
+          if (d2 < 0.0) d2 += prd[cdim2];
+          else d2 -= prd[cdim2];
+        }
+      }
+    }
+
     r = sqrt(d1*d1 + d2*d2);
 
     rbin = static_cast<int> ((r - cradmin) * cinvrad);

src/compute_chunk_atom.h

@@ -51,7 +51,7 @@ class ComputeChunkAtom : public Compute {
   int which,binflag;
   int regionflag,nchunksetflag,nchunkflag,discard;
   int limit,limitstyle,limitfirst;
-  int scaleflag;
+  int scaleflag,pbcflag;
   double xscale,yscale,zscale;
   int argindex;
   char *cfvid;

src/fix_store.cpp

@@ -15,9 +15,9 @@
 #include <string.h>
 #include "fix_store.h"
 #include "atom.h"
+#include "force.h"
 #include "memory.h"
 #include "error.h"
-#include "force.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -80,16 +80,6 @@ int FixStore::setmask()
   return mask;
 }
 
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based array
-------------------------------------------------------------------------- */
-
-double FixStore::memory_usage()
-{
-  double bytes = atom->nmax*nvalues * sizeof(double);
-  return bytes;
-}
-
 /* ----------------------------------------------------------------------
    allocate atom-based array
 ------------------------------------------------------------------------- */
@@ -189,3 +179,14 @@ int FixStore::size_restart(int nlocal)
 {
   return nvalues+1;
 }
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixStore::memory_usage()
+{
+  double bytes = atom->nmax*nvalues * sizeof(double);
+  return bytes;
+}
+

src/fix_store.h

@@ -33,7 +33,6 @@ class FixStore : public Fix {
   ~FixStore();
   int setmask();
 
-  double memory_usage();
   void grow_arrays(int);
   void copy_arrays(int, int, int);
   int pack_exchange(int, double *);
@@ -43,6 +42,8 @@ class FixStore : public Fix {
   int size_restart(int);
   int maxsize_restart();
 
+  double memory_usage();
+
  private:
   int nvalues;                  // total # of values per atom
   int vecflag;                  // 1 if nvalues = 1

src/fix_store_state.cpp

@@ -105,6 +105,19 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
       if (domain->triclinic)
         pack_choice[nvalues++] = &FixStoreState::pack_zu_triclinic;
       else pack_choice[nvalues++] = &FixStoreState::pack_zu;
+    } else if (strcmp(arg[iarg],"xsu") == 0) {
+      if (domain->triclinic)
+        pack_choice[nvalues++] = &FixStoreState::pack_xsu_triclinic;
+      else pack_choice[nvalues++] = &FixStoreState::pack_xsu;
+    } else if (strcmp(arg[iarg],"ysu") == 0) {
+      if (domain->triclinic)
+        pack_choice[nvalues++] = &FixStoreState::pack_ysu_triclinic;
+      else pack_choice[nvalues++] = &FixStoreState::pack_ysu;
+    } else if (strcmp(arg[iarg],"zsu") == 0) {
+      if (domain->triclinic)
+        pack_choice[nvalues++] = &FixStoreState::pack_zsu_triclinic;
+      else pack_choice[nvalues++] = &FixStoreState::pack_zsu;
+
     } else if (strcmp(arg[iarg],"ix") == 0) {
       pack_choice[nvalues++] = &FixStoreState::pack_ix;
     } else if (strcmp(arg[iarg],"iy") == 0) {
@@ -145,12 +158,22 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,
                    "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_muz;
+    } else if (strcmp(arg[iarg],"mu") == 0) {
+      if (!atom->mu_flag)
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
+      pack_choice[nvalues++] = &FixStoreState::pack_mu;
 
     } else if (strcmp(arg[iarg],"radius") == 0) {
       if (!atom->radius_flag)
         error->all(FLERR,
                    "Fix store/state for atom property that isn't allocated");
       pack_choice[nvalues++] = &FixStoreState::pack_radius;
+    } else if (strcmp(arg[iarg],"diameter") == 0) {
+      if (!atom->radius_flag)
+        error->all(FLERR,
+                   "Fix store/state for atom property that isn't allocated");
+      pack_choice[nvalues++] = &FixStoreState::pack_diameter;
     } else if (strcmp(arg[iarg],"omegax") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,
@@ -989,6 +1012,128 @@ void FixStoreState::pack_zu_triclinic(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void FixStoreState::pack_xsu(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double boxxlo = domain->boxlo[0];
+  double invxprd = 1.0/domain->xprd;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      vbuf[n] = (x[i][0]-boxxlo)*invxprd + ((image[i] & IMGMASK) - IMGMAX);
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixStoreState::pack_ysu(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double boxylo = domain->boxlo[1];
+  double invyprd = 1.0/domain->yprd;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      vbuf[n] = (x[i][1]-boxylo)*invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixStoreState::pack_zsu(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double boxzlo = domain->boxlo[2];
+  double invzprd = 1.0/domain->zprd;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      vbuf[n] = (x[i][2]-boxzlo)*invzprd + (image[i] >> IMG2BITS) - IMGMAX;
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_xsu_triclinic(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double *boxlo = domain->boxlo;
+  double *h_inv = domain->h_inv;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      vbuf[n] = h_inv[0]*(x[i][0]-boxlo[0]) + h_inv[5]*(x[i][1]-boxlo[1]) +
+	h_inv[4]*(x[i][2]-boxlo[2]) + (image[i] & IMGMASK) - IMGMAX;
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_ysu_triclinic(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double *boxlo = domain->boxlo;
+  double *h_inv = domain->h_inv;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      vbuf[n] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]) +
+	(image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_zsu_triclinic(int n)
+{
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double *boxlo = domain->boxlo;
+  double *h_inv = domain->h_inv;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) 
+      buf[n] = h_inv[2]*(x[i][2]-boxlo[2]) + (image[i] >> IMG2BITS) - IMGMAX;
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixStoreState::pack_ix(int n)
 {
   imageint *image = atom->image;
@@ -1184,6 +1329,21 @@ void FixStoreState::pack_muz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void FixStoreState::pack_mu(int n)
+{
+  double **mu = atom->mu;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) vbuf[n] = mu[i][3];
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixStoreState::pack_radius(int n)
 {
   double *radius = atom->radius;
@@ -1199,6 +1359,21 @@ void FixStoreState::pack_radius(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void FixStoreState::pack_diameter(int n)
+{
+  double *radius = atom->radius;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) vbuf[n] = 2.0*radius[i];
+    else vbuf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixStoreState::pack_omegax(int n)
 {
   double **omega = atom->omega;

src/fix_store_state.h

@@ -79,6 +79,13 @@ class FixStoreState : public Fix {
   void pack_xu_triclinic(int);
   void pack_yu_triclinic(int);
   void pack_zu_triclinic(int);
+  void pack_xsu(int);
+  void pack_ysu(int);
+  void pack_zsu(int);
+  void pack_xsu_triclinic(int);
+  void pack_ysu_triclinic(int);
+  void pack_zsu_triclinic(int);
+
   void pack_ix(int);
   void pack_iy(int);
   void pack_iz(int);
@@ -93,7 +100,9 @@ class FixStoreState : public Fix {
   void pack_mux(int);
   void pack_muy(int);
   void pack_muz(int);
+  void pack_mu(int);
   void pack_radius(int);
+  void pack_diameter(int);
   void pack_omegax(int);
   void pack_omegay(int);
   void pack_omegaz(int);