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

src/ASPHERE/pair_line_lj.cpp

@@ -52,7 +52,10 @@ PairLineLJ::~PairLineLJ()
     memory->destroy(setflag);
     memory->destroy(cutsq);
 
+    memory->destroy(subsize);
     memory->destroy(cut);
+    memory->destroy(cutsub);
+    memory->destroy(cutsubsq);
     memory->destroy(epsilon);
     memory->destroy(sigma);
     memory->destroy(lj1);
@@ -66,7 +69,7 @@ PairLineLJ::~PairLineLJ()
 
 void PairLineLJ::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum,itype,jtype;
+  int i,j,ii,jj,inum,jnum,itype,jtype,tmp;
   int ni,nj,npi,npj,ifirst,jfirst;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,r2inv,r6inv,term1,term2,sig,sig3,forcelj;
@@ -130,10 +133,10 @@ void PairLineLJ::compute(int eflag, int vflag)
 
       evdwl = 0.0;
       if (line[i] >= 0 && line[j] >= 0) {
-        if (dnum[i] == 0) discretize(i,sigma[itype][itype]);
+        if (dnum[i] == 0) discretize(i,subsize[itype]);
         npi = dnum[i];
         ifirst = dfirst[i];
-        if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]);
+        if (dnum[j] == 0) discretize(j,subsize[jtype]);
         npj = dnum[j];
         jfirst = dfirst[j];
 
@@ -154,7 +157,11 @@ void PairLineLJ::compute(int eflag, int vflag)
             dely = xi[1] - xj[1];
             rsq = delx*delx + dely*dely;
 
-            sig = 0.5 * (discrete[ifirst+ni].sigma+discrete[jfirst+nj].sigma);
+            // skip this pair of sub-particles if outside sub cutoff
+
+            if (rsq >= cutsubsq[itype][jtype]) continue;
+
+            sig = sigma[itype][jtype];
             sig3 = sig*sig*sig;
             term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
             term1 = 2.0 * term2 * sig3*sig3;
@@ -167,16 +174,15 @@ void PairLineLJ::compute(int eflag, int vflag)
 
             fi[0] = delx*fpair;
             fi[1] = dely*fpair;
+
             f[i][0] += fi[0];
             f[i][1] += fi[1];
             torque[i][2] += dxi*fi[1] - dyi*fi[0];
 
             if (newton_pair || j < nlocal) {
-              fj[0] = -delx*fpair;
+              f[j][0] -= fi[0];
-              fj[1] = -dely*fpair;
+              f[j][1] -= fi[1];
-              f[j][0] += fj[0];
+              torque[j][2] -= dxj*fi[1] - dyj*fi[0];
-              f[j][1] += fj[1];
-              torque[j][2] += dxj*fj[1] - dyj*fj[0];
             }
           }
         }
@@ -185,7 +191,7 @@ void PairLineLJ::compute(int eflag, int vflag)
       // convert line into Np particles based on sigma and line length
 
       } else if (line[i] >= 0) {
-        if (dnum[i] == 0) discretize(i,sigma[itype][itype]);
+        if (dnum[i] == 0) discretize(i,subsize[itype]);
         npi = dnum[i];
         ifirst = dfirst[i];
 
@@ -202,7 +208,11 @@ void PairLineLJ::compute(int eflag, int vflag)
           dely = xi[1] - xj[1];
           rsq = delx*delx + dely*dely;
 
-          sig = 0.5 * (discrete[ifirst+ni].sigma+sigma[jtype][jtype]);
+          // skip this pair of sub-particles if outside sub cutoff
+
+          if (rsq >= cutsubsq[itype][jtype]) continue;
+
+          sig = sigma[itype][jtype];
           sig3 = sig*sig*sig;
           term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
           term1 = 2.0 * term2 * sig3*sig3;
@@ -220,10 +230,8 @@ void PairLineLJ::compute(int eflag, int vflag)
           torque[i][2] += dxi*fi[1] - dyi*fi[0];
 
           if (newton_pair || j < nlocal) {
-            fj[0] = -delx*fpair;
+            f[j][0] -= fi[0];
-            fj[1] = -dely*fpair;
+            f[j][1] -= fi[1];
-            f[j][0] += fj[0];
-            f[j][1] += fj[1];
           }
         }
 
@@ -231,7 +239,7 @@ void PairLineLJ::compute(int eflag, int vflag)
       // convert line into Np particles based on sigma and line length
 
       } else if (line[j] >= 0) {
-        if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]);
+        if (dnum[j] == 0) discretize(j,subsize[jtype]);
         npj = dnum[j];
         jfirst = dfirst[j];
 
@@ -248,7 +256,11 @@ void PairLineLJ::compute(int eflag, int vflag)
           dely = xi[1] - xj[1];
           rsq = delx*delx + dely*dely;
 
-          sig = 0.5 * (sigma[itype][itype]+discrete[jfirst+nj].sigma);
+          // skip this pair of sub-particles if outside sub cutoff
+
+          if (rsq >= cutsubsq[itype][jtype]) continue;
+
+          sig = sigma[itype][jtype];
           sig3 = sig*sig*sig;
           term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
           term1 = 2.0 * term2 * sig3*sig3;
@@ -265,11 +277,9 @@ void PairLineLJ::compute(int eflag, int vflag)
           f[i][1] += fi[1];
 
           if (newton_pair || j < nlocal) {
-            f[j][0] += fj[0];
+            f[j][0] -= fi[0];
-            f[j][1] += fj[1];
+            f[j][1] -= fi[1];
-            fj[0] = -delx*fpair;
+            torque[j][2] -= dxj*fi[1] - dyj*fi[0];
-            fj[1] = -dely*fpair;
-            torque[j][2] += dxj*fj[1] - dyj*fj[0];
           }
         }
 
@@ -318,7 +328,10 @@ void PairLineLJ::allocate()
 
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
 
+  memory->create(subsize,n+1,"pair:subsize");
   memory->create(cut,n+1,n+1,"pair:cut");
+  memory->create(cutsub,n+1,n+1,"pair:cutsub");
+  memory->create(cutsubsq,n+1,n+1,"pair:cutsubsq");
   memory->create(epsilon,n+1,n+1,"pair:epsilon");
   memory->create(sigma,n+1,n+1,"pair:sigma");
   memory->create(lj1,n+1,n+1,"pair:lj1");
@@ -353,7 +366,7 @@ void PairLineLJ::settings(int narg, char **arg)
 
 void PairLineLJ::coeff(int narg, char **arg)
 {
-  if (narg < 4 || narg > 5)
+  if (narg < 7 || narg > 8)
     error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
@@ -361,17 +374,23 @@ void PairLineLJ::coeff(int narg, char **arg)
   force->bounds(arg[0],atom->ntypes,ilo,ihi);
   force->bounds(arg[1],atom->ntypes,jlo,jhi);
 
-  double epsilon_one = force->numeric(FLERR,arg[2]);
+  double size_itype = force->numeric(FLERR,arg[2]);
-  double sigma_one = force->numeric(FLERR,arg[3]);
+  double size_jtype = force->numeric(FLERR,arg[3]);
+  double epsilon_one = force->numeric(FLERR,arg[4]);
+  double sigma_one = force->numeric(FLERR,arg[5]);
+  double cutsub_one = force->numeric(FLERR,arg[6]);
 
   double cut_one = cut_global;
-  if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
+  if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
+      subsize[i] = size_itype;
+      subsize[j] = size_jtype;
       epsilon[i][j] = epsilon_one;
       sigma[i][j] = sigma_one;
+      cutsub[i][j] = cutsub_one;
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -399,12 +418,9 @@ void PairLineLJ::init_style()
 
 double PairLineLJ::init_one(int i, int j)
 {
-  if (setflag[i][j] == 0) {
+  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
-    epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+
-                               sigma[i][i],sigma[j][j]);
+  cutsubsq[i][j] = cutsub[i][j] * cutsub[i][j];
-    sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
-    cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
-  }
 
   lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
   lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
@@ -413,6 +429,8 @@ double PairLineLJ::init_one(int i, int j)
 
   epsilon[j][i] = epsilon[i][j];
   sigma[j][i] = sigma[i][j];
+  cutsubsq[j][i] = cutsubsq[i][j];
+
   lj1[j][i] = lj1[i][j];
   lj2[j][i] = lj2[i][j];
   lj3[j][i] = lj3[i][j];
@@ -422,16 +440,17 @@ double PairLineLJ::init_one(int i, int j)
 }
 
 /* ----------------------------------------------------------------------
-   discretize line segment I into N sub-segments no more than sigma in length
+   discretize line segment I into N sub-particles with <= size separation
-   store new discrete particles in Discrete list
+   store displacement dx,dy of discrete particles in Discrete list
 ------------------------------------------------------------------------- */
 
-void PairLineLJ::discretize(int i, double sigma)
+void PairLineLJ::discretize(int i, double size)
 {
   AtomVecLine::Bonus *bonus = avec->bonus;
   double length = bonus[atom->line[i]].length;
   double theta = bonus[atom->line[i]].theta;
-  int n = static_cast<int> (length/sigma) + 1;
+  int n = static_cast<int> (length/size) + 1;
+  
   dnum[i] = n;
   dfirst[i] = ndiscrete;
 
@@ -441,14 +460,12 @@ void PairLineLJ::discretize(int i, double sigma)
       memory->srealloc(discrete,dmax*sizeof(Discrete),"pair:discrete");
   }
 
-  sigma = length/n;
   double delta;
 
   for (int m = 0; m < n; m++) {
     delta = -0.5 + (2*m+1)/(2.0*n);
     discrete[ndiscrete].dx = delta*length*cos(theta);
     discrete[ndiscrete].dy = delta*length*sin(theta);
-    discrete[ndiscrete].sigma = sigma;
     ndiscrete++;
   }
 }

src/ASPHERE/pair_line_lj.h

@@ -36,14 +36,16 @@ class PairLineLJ : public Pair {
 
  protected:
   double cut_global;
+  double *subsize;
+  double **epsilon,**sigma,**cutsub,**cutsubsq;
   double **cut;
-  double **epsilon,**sigma;
+  double **lj1,**lj2,**lj3,**lj4;     // for sphere/sphere interactions
-  double **lj1,**lj2,**lj3,**lj4;
   class AtomVecLine *avec;
 
+  double *size;     // per-type size of sub-particles to tile line segment
+
   struct Discrete {
     double dx,dy;
-    double sigma;
   };
   Discrete *discrete;           // list of all discretes for all lines
   int ndiscrete;                // number of discretes in list

src/displace_atoms.cpp

@@ -24,11 +24,16 @@
 #include "group.h"
 #include "math_const.h"
 #include "random_park.h"
-#include "error.h"
 #include "force.h"
 #include "input.h"
 #include "variable.h"
+#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
+#include "atom_vec_body.h"
+#include "math_extra.h"
 #include "memory.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -203,8 +208,10 @@ void DisplaceAtoms::command(int narg, char **arg)
   // X = P + C + A cos(theta) + B sin(theta)
 
   if (style == ROTATE) {
-    double axis[3],point[3];
+    double theta_new;
+    double axis[3],point[3],qrotate[4],qnew[4];
     double a[3],b[3],c[3],d[3],disp[3],runit[3];
+    double *quat;
 
     int dim = domain->dimension;
     point[0] = xscale*force->numeric(FLERR,arg[2]);
@@ -224,11 +231,36 @@ void DisplaceAtoms::command(int narg, char **arg)
     runit[1] = axis[1]/len;
     runit[2] = axis[2]/len;
 
-    double sine = sin(MY_PI*theta/180.0);
+    double angle = MY_PI*theta/180.0;
-    double cosine = cos(MY_PI*theta/180.0);
+    double sine = sin(angle);
+    double cosine = cos(angle);
     double ddotr;
 
+    // flags for additional orientation info stored by some atom styles
+
+    int ellipsoid_flag = atom->ellipsoid_flag;
+    int line_flag = atom->line_flag;
+    int tri_flag = atom->tri_flag;
+    int body_flag = atom->body_flag;
+
+    int theta_flag = 0;
+    int quat_flag = 0;
+    if (line_flag) theta_flag = 1;
+    if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
+
+    // AtomVec pointers to retrieve per-atom storage of extra quantities
+    
+    AtomVecEllipsoid *avec_ellipsoid = 
+      (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+    AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
+    AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
+    AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
+
     double **x = atom->x;
+    int *ellipsoid = atom->ellipsoid;
+    int *line = atom->line;
+    int *tri = atom->tri;
+    int *body = atom->body;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
@@ -253,6 +285,36 @@ void DisplaceAtoms::command(int narg, char **arg)
         x[i][0] = point[0] + c[0] + disp[0];
         x[i][1] = point[1] + c[1] + disp[1];
         if (dim == 3) x[i][2] = point[2] + c[2] + disp[2];
+
+        // theta for lines
+          
+        if (theta_flag && line[i] >= 0.0) {
+          theta_new = fmod(avec_line->bonus[line[i]].theta+angle,MY_2PI);
+          avec_line->bonus[atom->line[i]].theta = theta_new;
+        }
+          
+        // quats for ellipsoids, tris, and bodies
+        
+        if (quat_flag) {
+          quat = NULL;
+          if (ellipsoid_flag && ellipsoid[i] >= 0)
+            quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
+          else if (tri_flag && tri[i] >= 0)
+            quat = avec_tri->bonus[tri[i]].quat;
+          else if (body_flag && body[i] >= 0)
+            quat = avec_body->bonus[body[i]].quat;
+          if (quat) {
+            qrotate[0] = cosine;
+            qrotate[1] = runit[0]*sine;
+            qrotate[2] = runit[1]*sine;
+            qrotate[3] = runit[2]*sine;
+            MathExtra::quatquat(qrotate,quat,qnew);
+            quat[0] = qnew[0];
+            quat[1] = qnew[1];
+            quat[2] = qnew[2];
+            quat[3] = qnew[3];
+          }
+        }
       }
     }
   }

src/dump_image.cpp

@@ -18,7 +18,8 @@
 #include "dump_image.h"
 #include "image.h"
 #include "atom.h"
-#include "atom_vec.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
 #include "atom_vec_body.h"
 #include "body.h"
 #include "molecule.h"
@@ -29,6 +30,7 @@
 #include "input.h"
 #include "variable.h"
 #include "math_const.h"
+#include "math_extra.h"
 #include "error.h"
 #include "memory.h"
 
@@ -106,7 +108,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
   // set defaults for optional args
 
   atomflag = YES;
-  bodyflag = NO;
+  lineflag = triflag = bodyflag = NO;
   if (atom->nbondtypes == 0) bondflag = NO;
   else {
     bondflag = YES;
@@ -132,28 +134,19 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 
   int iarg = ioptional;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"adiam") == 0) {
+    if (strcmp(arg[iarg],"atom") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
-      adiam = NUMERIC;
-      adiamvalue = force->numeric(FLERR,arg[iarg+1]);
-      if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
-      iarg += 2;
-
-    } else if (strcmp(arg[iarg],"atom") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES;
       else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO;
       else error->all(FLERR,"Illegal dump image command");
       iarg += 2;
 
-    } else if (strcmp(arg[iarg],"body") == 0) {
+    } else if (strcmp(arg[iarg],"adiam") == 0) {
-      if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
+      if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
-      if (strcmp(arg[iarg+1],"yes") == 0) bodyflag = YES;
+      adiam = NUMERIC;
-      else if (strcmp(arg[iarg+1],"no") == 0) bodyflag = NO;
+      adiamvalue = force->numeric(FLERR,arg[iarg+1]);
-      else error->all(FLERR,"Illegal dump image command");
+      if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
-      bodyflag1 = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 2;
-      bodyflag2 = force->numeric(FLERR,arg[iarg+3]);
-      iarg += 4;
 
     } else if (strcmp(arg[iarg],"bond") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
@@ -174,6 +167,32 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
       else error->all(FLERR,"Illegal dump image command");
       iarg += 3;
 
+
+    } else if (strcmp(arg[iarg],"line") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
+      lineflag = YES;
+      if (strcmp(arg[iarg+1],"type") == 0) lcolor = TYPE;
+      else error->all(FLERR,"Illegal dump image command");
+      ldiam = NUMERIC;
+      ldiamvalue = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+
+    } else if (strcmp(arg[iarg],"tri") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
+      triflag = YES;
+      if (strcmp(arg[iarg+1],"type") == 0) tcolor = TYPE;
+      else error->all(FLERR,"Illegal dump image command");
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"body") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
+      bodyflag = YES;
+      if (strcmp(arg[iarg+1],"type") == 0) bodycolor = TYPE;
+      else error->all(FLERR,"Illegal dump image command");
+      bodyflag1 = force->numeric(FLERR,arg[iarg+2]);
+      bodyflag2 = force->numeric(FLERR,arg[iarg+3]);
+      iarg += 4;
+
     } else if (strcmp(arg[iarg],"size") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
       int width = force->inumeric(FLERR,arg[iarg+1]);
@@ -333,13 +352,26 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
     } else error->all(FLERR,"Illegal dump image command");
   }
 
-  // error check for bodyflag
+  // error checks and setup for lineflag, triflag, bodyflag
 
-  if (bodyflag) {
+  if (lineflag) {
-    AtomVecBody *avec = (AtomVecBody *) atom->style_match("body");
+    avec_line = (AtomVecLine *) atom->style_match("line");
-    if (!avec) error->all(FLERR,"Dump image body yes requires atom style body");
+    if (!avec_line) 
-    bptr = avec->bptr;
+      error->all(FLERR,"Dump image line requires atom style line");
   }
+  if (triflag) {
+    avec_tri = (AtomVecTri *) atom->style_match("tri");
+    if (!avec_tri) 
+      error->all(FLERR,"Dump image tri requires atom style tri");
+  }
+  if (bodyflag) {
+    avec_body = (AtomVecBody *) atom->style_match("body");
+    if (!avec_body) 
+      error->all(FLERR,"Dump image body yes requires atom style body");
+  }
+
+  extraflag = 0;
+  if (lineflag || triflag || bodyflag) extraflag = 1;
 
   // allocate image buffer now that image size is known
 
@@ -679,18 +711,21 @@ void DumpImage::view_params()
 
 void DumpImage::create_image()
 {
-  int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus;
+  int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus,drawflag;
   tagint tagprev;
   double diameter,delx,dely,delz;
   int *bodyvec;
   double **bodyarray;
   double *color,*color1,*color2;
-  double xmid[3];
+  double xmid[3],pt1[3],pt2[3],pt3[3];
+  double mat[3][3];
 
   // render my atoms
 
   if (atomflag) {
     double **x = atom->x;
+    int *line = atom->line;
+    int *tri = atom->tri;
     int *body = atom->body;
 
     m = 0;
@@ -719,9 +754,101 @@ void DumpImage::create_image()
         diameter = buf[m+1];
       }
 
-      if (!body || !bodyflag || body[j] < 0)
+      // do not draw if line,tri,body keywords enabled and atom is one of those
-        image->draw_sphere(x[j],color,diameter);
+
-      else {
+      drawflag = 1;
+      if (extraflag) {
+        if (lineflag && line[j] >= 0) drawflag = 0;
+        if (triflag && tri[j] >= 0) drawflag = 0;
+        if (bodyflag && body[j] >= 0) drawflag = 0;
+      }
+
+      if (drawflag) image->draw_sphere(x[j],color,diameter);
+
+      m += size_one;
+    }
+  }
+
+  // render atoms that are lines
+
+  if (lineflag) {
+    double length,theta,dx,dy;
+    double **x = atom->x;
+    int *line = atom->line;
+    int *type = atom->type;
+
+    for (i = 0; i < nchoose; i++) {
+      j = clist[i];
+      if (line[j] < 0) continue;
+
+      if (lcolor == TYPE) {
+        color = colortype[type[j]];
+      }
+
+      if (ldiam == NUMERIC) {
+        diameter = ldiamvalue;
+      }
+
+      length = avec_line->bonus[line[j]].length;
+      theta = avec_line->bonus[line[j]].theta;
+      dx = 0.5*length*cos(theta);
+      dy = 0.5*length*sin(theta);
+
+      pt1[0] = x[j][0] + dx;
+      pt1[1] = x[j][1] + dy;
+      pt1[2] = 0.0;
+      pt2[0] = x[j][0] - dx;
+      pt2[1] = x[j][1] - dy;
+      pt2[2] = 0.0;
+
+      image->draw_cylinder(pt1,pt2,color,ldiamvalue,3);
+    }
+  }
+
+  // render atoms that are triangles
+
+  if (triflag) {
+    double **x = atom->x;
+    int *tri = atom->tri;
+    int *type = atom->type;
+
+    for (i = 0; i < nchoose; i++) {
+      j = clist[i];
+      if (tri[j] < 0) continue;
+
+      if (tcolor == TYPE) {
+        color = colortype[type[j]];
+      }
+
+      MathExtra::quat_to_mat(avec_tri->bonus[tri[i]].quat,mat);
+      MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c1,pt1);
+      MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c2,pt2);
+      MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c3,pt3);
+      MathExtra::add3(pt1,x[i],pt1);
+      MathExtra::add3(pt2,x[i],pt2);
+      MathExtra::add3(pt3,x[i],pt3);
+
+      image->draw_triangle(pt1,pt2,pt3,color);
+    }
+  }
+
+  // render atoms that are bodies
+
+  if (bodyflag) {
+    Body *bptr = avec_body->bptr;
+    double **x = atom->x;
+    int *body = atom->body;
+
+    m = 0;
+    for (i = 0; i < nchoose; i++) {
+      j = clist[i];
+      if (body[j] < 0) continue;
+
+      if (bodycolor == TYPE) {
+        itype = static_cast<int> (buf[m]);
+        color = colortype[itype];
+      }
+
       ibonus = body[i];
       n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray);
       for (k = 0; k < n; k++) {
@@ -731,7 +858,6 @@ void DumpImage::create_image()
           image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3],
                                color,bodyarray[k][6],3);
       }
-      }
 
       m += size_one;
     }

src/dump_image.h

@@ -37,13 +37,24 @@ class DumpImage : public DumpCustom {
   int filetype;
   enum{PPM,JPG,PNG};
 
+  int atomflag;                    // 0/1 for draw atoms
   int acolor,adiam;                // what determines color/diam of atoms
   double adiamvalue;               // atom diameter value
-  int atomflag,bondflag;           // 0/1 for draw atoms,bonds
+
+  int lineflag;                    // 0/1 for draw atoms as lines
+  int lcolor,ldiam;                // what determines color/diam of lines
+  double ldiamvalue;               // line diameter value
+  int triflag;                     // 0/1 for draw atoms as triangles
+  int tcolor;                      // what determines color of tris
   int bodyflag;                    // 0/1 for draw atoms as bodies
-  double bodyflag1,bodyflag2;      // user params for drawing bodies
+  int bodycolor;                   // what determines color of bodies
+  double bodyflag1,bodyflag2;      // user-specified params for drawing bodies
+  
+  int bondflag;                    // 0/1 for draw bonds
   int bcolor,bdiam;                // what determines color/diam of bonds
   double bdiamvalue;               // bond diameter value
+
+  int extraflag;                   // 0/1 for any of line/tri/body flag set
   char *thetastr,*phistr;          // variables for view theta,phi
   int thetavar,phivar;             // index to theta,phi vars
   int cflag;                       // static/dynamic box center
@@ -64,8 +75,11 @@ class DumpImage : public DumpCustom {
   double *diamtype,*diamelement,*bdiamtype;         // per-type diameters
   double **colortype,**colorelement,**bcolortype;   // per-type colors
 
+  class AtomVecLine *avec_line;    // ptrs to atom style (sub)classes
+  class AtomVecTri *avec_tri;
+  class AtomVecBody *avec_body;
+
   class Image *image;              // class that renders each image
-  class Body *bptr;                // class for Body particles
 
   int *chooseghost;                // extended choose array for comm
   double **bufcopy;                // buffer for communicating bond/atom info

src/fix_move.cpp

@@ -26,7 +26,12 @@
 #include "respa.h"
 #include "input.h"
 #include "variable.h"
+#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
+#include "atom_vec_body.h"
 #include "math_const.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "error.h"
 
@@ -37,6 +42,8 @@ using namespace MathConst;
 enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
 enum{EQUAL,ATOM};
 
+#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
+
 /* ---------------------------------------------------------------------- */
 
 FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
@@ -216,7 +223,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
 
   // set omega_rotate from period
 
-  if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = 2.0*MY_PI / period;
+  if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period;
 
   // runit = unit vector along rotation axis
 
@@ -229,24 +236,54 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     runit[2] = axis[2]/len;
   }
 
-  // set omega_flag if particles store omega
+  // set flags for extra attributes particles may store
+  // relevant extra attributes = omega, angmom, theta, quat
   
   omega_flag = atom->omega_flag;
+  angmom_flag = atom->angmom_flag;
+
+  radius_flag = atom->radius_flag;
+  ellipsoid_flag = atom->ellipsoid_flag;
+  line_flag = atom->line_flag;
+  tri_flag = atom->tri_flag;
+  body_flag = atom->body_flag;
+
+  theta_flag = quat_flag = 0;
+  if (line_flag) theta_flag = 1;
+  if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
+
+  extra_flag = 0;
+  if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1;
     
   // perform initial allocation of atom-based array
   // register with Atom class
 
   xoriginal = NULL;
+  toriginal = NULL;
+  qoriginal = NULL;
   grow_arrays(atom->nmax);
   atom->add_callback(0);
   atom->add_callback(1);
 
   displace = velocity = NULL;
 
+  // AtomVec pointers to retrieve per-atom storage of extra quantities
+
+  avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+  avec_line = (AtomVecLine *) atom->style_match("line");
+  avec_tri = (AtomVecTri *) atom->style_match("tri");
+  avec_body = (AtomVecBody *) atom->style_match("body");
+
   // xoriginal = initial unwrapped positions of atoms
+  // toriginal = initial theta of lines
+  // qoriginal = initial quat of extended particles
 
   double **x = atom->x;
   imageint *image = atom->image;
+  int *ellipsoid = atom->ellipsoid;
+  int *line = atom->line;
+  int *tri = atom->tri;
+  int *body = atom->body;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
@@ -255,6 +292,45 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
     else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
   }
 
+  if (theta_flag) {
+    for (int i = 0; i < nlocal; i++) {
+      if ((mask[i] & groupbit) && line[i] >= 0)
+        toriginal[i] = avec_line->bonus[line[i]].theta;
+      else toriginal[i] = 0.0;
+    }
+  }
+
+  if (quat_flag) {
+    double *quat;
+    for (int i = 0; i < nlocal; i++) {
+      quat = NULL;
+      if (mask[i] & groupbit) {
+        if (ellipsoid_flag && ellipsoid[i] >= 0)
+          quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
+        else if (tri_flag && tri[i] >= 0)
+          quat = avec_tri->bonus[tri[i]].quat;
+        else if (body_flag && body[i] >= 0)
+          quat = avec_body->bonus[body[i]].quat;
+      }
+      if (quat) {
+        qoriginal[i][0] = quat[0];
+        qoriginal[i][1] = quat[1];
+        qoriginal[i][2] = quat[2];
+        qoriginal[i][3] = quat[3];
+      } else qoriginal[i][0] = qoriginal[i][1] = 
+               qoriginal[i][2] = qoriginal[i][3] = 0.0;
+    }
+  }
+
+  // nrestart = size of per-atom restart data
+  // nrestart = 1 + xorig + torig + qorig
+
+  nrestart = 4;
+  if (theta_flag) nrestart++;
+  if (quat_flag) nrestart += 4;
+
+  // time origin for movement = current timestep
+
   time_origin = update->ntimestep;
 }
 
@@ -270,6 +346,8 @@ FixMove::~FixMove()
   // delete locally stored arrays
 
   memory->destroy(xoriginal);
+  memory->destroy(toriginal);
+  memory->destroy(qoriginal);
   memory->destroy(displace);
   memory->destroy(velocity);
 
@@ -381,9 +459,12 @@ void FixMove::init()
 
 void FixMove::initial_integrate(int vflag)
 {
-  double dtfm;
+  int flag;
-  double xold[3],a[3],b[3],c[3],d[3],disp[3];
   double ddotr,dx,dy,dz;
+  double dtfm,theta_new;
+  double xold[3],a[3],b[3],c[3],d[3],disp[3],w[3],ex[3],ey[3],ez[3];
+  double inertia_ellipsoid[3],qrotate[4];
+  double *quat,*inertia,*shape;
 
   double delta = (update->ntimestep - time_origin) * dt;
 
@@ -391,10 +472,17 @@ void FixMove::initial_integrate(int vflag)
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
+  double **angmom = atom->angmom;
+  double *radius = atom->radius;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
+  int *ellipsoid = atom->ellipsoid;
+  int *line = atom->line;
+  int *tri = atom->tri;
+  int *body = atom->body;
   int *mask = atom->mask;
+
   int nlocal = atom->nlocal;
 
   // for linear: X = X0 + V*dt
@@ -553,17 +641,89 @@ void FixMove::initial_integrate(int vflag)
         v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]);
         v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]);
         v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]);
+
+        // set any extra attributes affected by rotation
+
+        if (extra_flag) {
+
+          // omega for spheres and lines
+
           if (omega_flag) {
+            flag = 0;
+            if (radius_flag && radius[i] > 0.0) flag = 1;
+            if (line_flag && line[i] >= 0.0) flag = 1;
+            if (flag) {
               omega[i][0] = omega_rotate*runit[0];
               omega[i][1] = omega_rotate*runit[1];
               omega[i][2] = omega_rotate*runit[2];
             }
+          }
+
+          // angmom for ellipsoids, tris, and bodies
+          
+          if (angmom_flag) {
+            quat = inertia = NULL;
+            if (ellipsoid_flag && ellipsoid[i] >= 0) {
+              quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
+              shape = avec_ellipsoid->bonus[ellipsoid[i]].shape;
+              inertia_ellipsoid[0] = 
+                INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+              inertia_ellipsoid[1] = 
+                INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+              inertia_ellipsoid[2] = 
+                INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+              inertia = inertia_ellipsoid;
+            } else if (tri_flag && tri[i] >= 0) {
+              quat = avec_tri->bonus[tri[i]].quat;
+              inertia = avec_tri->bonus[tri[i]].inertia;
+            } else if (body_flag && body[i] >= 0) {
+              quat = avec_body->bonus[body[i]].quat;
+              inertia = avec_body->bonus[body[i]].inertia;
+            }
+            if (quat) {
+              w[0] = omega_rotate*runit[0];
+              w[1] = omega_rotate*runit[1];
+              w[2] = omega_rotate*runit[2];
+              MathExtra::q_to_exyz(quat,ex,ey,ez);
+              MathExtra::omega_to_angmom(w,ex,ey,ez,inertia,angmom[i]);
+            }
+          }
+
+          // theta for lines
+          
+          if (theta_flag && line[i] >= 0.0) {
+            theta_new = fmod(toriginal[i]+arg,MY_2PI);
+            avec_line->bonus[atom->line[i]].theta = theta_new;
+          }
+          
+          // quats for ellipsoids, tris, and bodies
+          
+          if (quat_flag) {
+            quat = NULL;
+            if (ellipsoid_flag && ellipsoid[i] >= 0)
+              quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
+            else if (tri_flag && tri[i] >= 0)
+              quat = avec_tri->bonus[tri[i]].quat;
+            else if (body_flag && body[i] >= 0)
+              quat = avec_body->bonus[body[i]].quat;
+            if (quat) {
+              qrotate[0] = cosine;
+              qrotate[1] = runit[0]*sine;
+              qrotate[2] = runit[1]*sine;
+              qrotate[3] = runit[2]*sine;
+              MathExtra::quatquat(qrotate,qoriginal[i],quat);
+            }
+          }
+          }
 
         domain->remap_near(x[i],xold);
       }
     }
 
   // for variable: compute x,v from variables
+  // NOTE: also allow for changes to extra attributes?
+  //       omega, angmom, theta, quat
+  //       only necessary if prescribed motion involves rotation
 
   } else if (mstyle == VARIABLE) {
 
@@ -809,6 +969,8 @@ void FixMove::final_integrate_respa(int ilevel, int iloop)
 double FixMove::memory_usage()
 {
   double bytes = atom->nmax*3 * sizeof(double);
+  if (theta_flag) bytes += atom->nmax * sizeof(double);
+  if (quat_flag) bytes += atom->nmax*4 * sizeof(double);
   if (displaceflag) bytes += atom->nmax*3 * sizeof(double);
   if (velocityflag) bytes += atom->nmax*3 * sizeof(double);
   return bytes;
@@ -850,6 +1012,8 @@ void FixMove::restart(char *buf)
 void FixMove::grow_arrays(int nmax)
 {
   memory->grow(xoriginal,nmax,3,"move:xoriginal");
+  if (theta_flag) memory->grow(toriginal,nmax,"move:toriginal");
+  if (quat_flag) memory->grow(qoriginal,nmax,4,"move:qoriginal");
   array_atom = xoriginal;
 }
 
@@ -862,6 +1026,13 @@ void FixMove::copy_arrays(int i, int j, int delflag)
   xoriginal[j][0] = xoriginal[i][0];
   xoriginal[j][1] = xoriginal[i][1];
   xoriginal[j][2] = xoriginal[i][2];
+  if (theta_flag) toriginal[j] = toriginal[i];
+  if (quat_flag) {
+    qoriginal[j][0] = qoriginal[i][0];
+    qoriginal[j][1] = qoriginal[i][1];
+    qoriginal[j][2] = qoriginal[i][2];
+    qoriginal[j][3] = qoriginal[i][3];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -870,8 +1041,15 @@ void FixMove::copy_arrays(int i, int j, int delflag)
 
 void FixMove::set_arrays(int i)
 {
+  double theta;
+  double *quat;
+
   double **x = atom->x;
   imageint *image = atom->image;
+  int *ellipsoid = atom->ellipsoid;
+  int *line = atom->line;
+  int *tri = atom->tri;
+  int *body = atom->body;
   int *mask = atom->mask;
 
   // particle not in group
@@ -932,6 +1110,33 @@ void FixMove::set_arrays(int i)
     xoriginal[i][0] = point[0] + c[0] + disp[0];
     xoriginal[i][1] = point[1] + c[1] + disp[1];
     xoriginal[i][2] = point[2] + c[2] + disp[2];
+    
+    // set theta and quat extra attributes affected by rotation
+
+    if (extra_flag) {
+
+      // theta for lines
+      
+      if (theta_flag && line[i] >= 0.0) {
+        theta = avec_line->bonus[atom->line[i]].theta;
+        toriginal[i] = theta - 0.0;  // NOTE: edit this line
+      }
+      
+      // quats for ellipsoids, tris, and bodies
+      
+      if (quat_flag) {
+        quat = NULL;
+        if (ellipsoid_flag && ellipsoid[i] >= 0)
+          quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
+        else if (tri_flag && tri[i] >= 0)
+          quat = avec_tri->bonus[tri[i]].quat;
+        else if (body_flag && body[i] >= 0)
+          quat = avec_body->bonus[body[i]].quat;
+        if (quat) {
+          // qoriginal = f(quat,-delta);   // NOTE: edit this line
+        }
+      }
+    }
   }
 }
 
@@ -941,10 +1146,18 @@ void FixMove::set_arrays(int i)
 
 int FixMove::pack_exchange(int i, double *buf)
 {
-  buf[0] = xoriginal[i][0];
+  int n = 0;
-  buf[1] = xoriginal[i][1];
+  buf[n++] = xoriginal[i][0];
-  buf[2] = xoriginal[i][2];
+  buf[n++] = xoriginal[i][1];
-  return 3;
+  buf[n++] = xoriginal[i][2];
+  if (theta_flag) buf[n++] = toriginal[i];
+  if (quat_flag) {
+    buf[n++] = qoriginal[i][0];
+    buf[n++] = qoriginal[i][1];
+    buf[n++] = qoriginal[i][2];
+    buf[n++] = qoriginal[i][3];
+  }
+  return n;
 }
 
 /* ----------------------------------------------------------------------
@@ -953,10 +1166,18 @@ int FixMove::pack_exchange(int i, double *buf)
 
 int FixMove::unpack_exchange(int nlocal, double *buf)
 {
-  xoriginal[nlocal][0] = buf[0];
+  int n = 0;
-  xoriginal[nlocal][1] = buf[1];
+  xoriginal[nlocal][0] = buf[n++];
-  xoriginal[nlocal][2] = buf[2];
+  xoriginal[nlocal][1] = buf[n++];
-  return 3;
+  xoriginal[nlocal][2] = buf[n++];
+  if (theta_flag) toriginal[nlocal] = buf[n++];
+  if (quat_flag) {
+    qoriginal[nlocal][0] = buf[n++];
+    qoriginal[nlocal][1] = buf[n++];
+    qoriginal[nlocal][2] = buf[n++];
+    qoriginal[nlocal][3] = buf[n++];
+  }
+  return n;
 }
 
 /* ----------------------------------------------------------------------
@@ -965,11 +1186,19 @@ int FixMove::unpack_exchange(int nlocal, double *buf)
 
 int FixMove::pack_restart(int i, double *buf)
 {
-  buf[0] = 4;
+  int n = 1;
-  buf[1] = xoriginal[i][0];
+  buf[n++] = xoriginal[i][0];
-  buf[2] = xoriginal[i][1];
+  buf[n++] = xoriginal[i][1];
-  buf[3] = xoriginal[i][2];
+  buf[n++] = xoriginal[i][2];
-  return 4;
+  if (theta_flag) buf[n++] = toriginal[i];
+  if (quat_flag) {
+    buf[n++] = qoriginal[i][0];
+    buf[n++] = qoriginal[i][1];
+    buf[n++] = qoriginal[i][2];
+    buf[n++] = qoriginal[i][3];
+  }
+  buf[0] = n;
+  return n;
 }
 
 /* ----------------------------------------------------------------------
@@ -989,6 +1218,13 @@ void FixMove::unpack_restart(int nlocal, int nth)
   xoriginal[nlocal][0] = extra[nlocal][m++];
   xoriginal[nlocal][1] = extra[nlocal][m++];
   xoriginal[nlocal][2] = extra[nlocal][m++];
+  if (theta_flag) toriginal[nlocal] = extra[nlocal][m++];
+  if (quat_flag) {
+    qoriginal[nlocal][0] = extra[nlocal][m++];
+    qoriginal[nlocal][1] = extra[nlocal][m++];
+    qoriginal[nlocal][2] = extra[nlocal][m++];
+    qoriginal[nlocal][3] = extra[nlocal][m++];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -997,7 +1233,7 @@ void FixMove::unpack_restart(int nlocal, int nth)
 
 int FixMove::maxsize_restart()
 {
-  return 4;
+  return nrestart;
 }
 
 /* ----------------------------------------------------------------------
@@ -1006,7 +1242,7 @@ int FixMove::maxsize_restart()
 
 int FixMove::size_restart(int nlocal)
 {
-  return 4;
+  return nrestart;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_move.h

@@ -61,13 +61,23 @@ class FixMove : public Fix {
   double dt,dtv,dtf;
   int xvar,yvar,zvar,vxvar,vyvar,vzvar;
   int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle;
-  int omega_flag,nlevels_respa;
+  int extra_flag,omega_flag,angmom_flag;
+  int radius_flag,ellipsoid_flag,line_flag,tri_flag,body_flag;
+  int theta_flag,quat_flag;
+  int nlevels_respa,nrestart;
   int time_origin;
 
   double **xoriginal;         // original coords of atoms
+  double *toriginal;          // original theta of atoms
+  double **qoriginal;         // original quat of atoms
   int displaceflag,velocityflag;
   int maxatom;
   double **displace,**velocity;
+
+  class AtomVecEllipsoid *avec_ellipsoid;
+  class AtomVecLine *avec_line;
+  class AtomVecTri *avec_tri;
+  class AtomVecBody *avec_body;
 };
 
 }

src/math_extra.cpp

@@ -230,6 +230,7 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq)
    apply evolution operators to quat, quat momentum
    Miller et al., J Chem Phys. 116, 8649-8659 (2002)
 ------------------------------------------------------------------------- */
+
 void no_squish_rotate(int k, double *p, double *q, double *inertia,
                       double dt)
 {

src/math_extra.h

@@ -76,7 +76,8 @@ namespace MathExtra {
   void rotate(double matrix[3][3], int i, int j, int k, int l,
               double s, double tau);
   void richardson(double *q, double *m, double *w, double *moments, double dtq);
-  void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt);
+  void no_squish_rotate(int k, double *p, double *q, double *inertia, 
+                        double dt);
 
   // shape matrix operations
   // upper-triangular 3x3 matrix stored in Voigt notation as 6-vector
@@ -151,7 +152,8 @@ inline void MathExtra::normalize3(const double *v, double *ans)
    scale a vector to length
 ------------------------------------------------------------------------- */
 
-inline void MathExtra::snormalize3(const double length, const double *v, double *ans)
+inline void MathExtra::snormalize3(const double length, const double *v, 
+                                   double *ans)
 {
   double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
   ans[0] = v[0]*scale;
@@ -601,7 +603,8 @@ inline void MathExtra::axisangle_to_quat(const double *v, const double angle,
    Apply principal rotation generator about x to rotation matrix m
 ------------------------------------------------------------------------- */
 
-inline void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3])
+inline void MathExtra::rotation_generator_x(const double m[3][3], 
+                                            double ans[3][3])
 {
   ans[0][0] = 0;
   ans[0][1] = -m[0][2];