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

src/GRANULAR/fix_pour.cpp

@@ -117,21 +117,23 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Cannot use fix_pour unless atoms have IDs");
 
   if (mode == MOLECULE) {
-    if (onemol->xflag == 0)
+    for (int i = 0; i < nmol; i++) {
-      error->all(FLERR,"Fix pour molecule must have coordinates");
+      if (onemol[i].xflag == 0)
-    if (onemol->typeflag == 0)
+        error->all(FLERR,"Fix pour molecule must have coordinates");
-      error->all(FLERR,"Fix pour molecule must have atom types");
+      if (onemol[i].typeflag == 0)
-    if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes)
+        error->all(FLERR,"Fix pour molecule must have atom types");
-      error->all(FLERR,"Invalid atom type in fix pour mol command");
+      if (ntype+onemol[i].ntypes <= 0 || ntype+onemol[i].ntypes > atom->ntypes)
+        error->all(FLERR,"Invalid atom type in fix pour mol command");
       
-    if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
+      if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
-      error->all(FLERR,"Fix pour molecule template ID must be same "
+        error->all(FLERR,"Fix pour molecule template ID must be same "
-                 "as atom style template ID");
+                   "as atom style template ID");
-    onemol->check_attributes(0);
+      onemol[i].check_attributes(0);
 
-    // fix pour uses geoemetric center of molecule for insertion
+      // fix pour uses geoemetric center of molecule for insertion
 
-    onemol->compute_center();
+      onemol[i].compute_center();
+    }
   }
 
   if (rigidflag && mode == ATOM)
@@ -144,7 +146,11 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
   // setup of coords and imageflags array
 
   if (mode == ATOM) natom = 1;
-  else natom = onemol->natoms;
+  else {
+    natom = 0;
+    for (int i = 0; i < nmol; i++)
+      natom = MAX(natom,onemol[i].natoms);
+  }
   memory->create(coords,natom,4,"pour:coords");
   memory->create(imageflags,natom,"pour:imageflags");
 
@@ -211,6 +217,12 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
 
   double volume,volume_one;
 
+  molradius_max = 0.0;
+  if (mode == MOLECULE) {
+    for (int i = 0; i < nmol; i++)
+      molradius_max = MAX(molradius_max,onemol[i].molradius);
+  }
+
   if (domain->dimension == 3) {
     if (region_style == 1) {
       double dy = yhi - ylo;
@@ -218,8 +230,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
       volume = (xhi-xlo) * dy * (zhi-zlo);
     } else volume = MY_PI*rc*rc * (zhi-zlo);
     if (mode == MOLECULE) {
-      double molradius = onemol->molradius;
+      volume_one = 4.0/3.0 * MY_PI * molradius_max*molradius_max*molradius_max;
-      volume_one = 4.0/3.0 * MY_PI * molradius*molradius*molradius;
     } else if (dstyle == ONE || dstyle == RANGE) {
       volume_one = 4.0/3.0 * MY_PI * radius_max*radius_max*radius_max;
     } else if (dstyle == POLY) {
@@ -231,8 +242,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
   } else {
     volume = (xhi-xlo) * (yhi-ylo);
     if (mode == MOLECULE) {
-      double molradius = onemol->molradius;
+      volume_one = MY_PI * molradius_max*molradius_max;
-      volume_one = MY_PI * molradius*molradius;
     } else if (dstyle == ONE || dstyle == RANGE) {
       volume_one = MY_PI * radius_max*radius_max;
     } else if (dstyle == POLY) {
@@ -264,6 +274,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
 FixPour::~FixPour()
 {
   delete random;
+  delete [] molfrac;
   delete [] idrigid;
   delete [] idshake;
   delete [] radius_poly;
@@ -356,7 +367,7 @@ void FixPour::init()
 
 void FixPour::pre_exchange()
 {
-  int i,j,m,flag,nlocalprev;
+  int i,j,m,flag,nlocalprev,imol;
   double r[3],rotmat[3][3],quat[4],vnew[3];
   double *newcoord;
 
@@ -478,6 +489,9 @@ void FixPour::pre_exchange()
         imageflags[0] = ((imageint) IMGMAX << IMG2BITS) |
             ((imageint) IMGMAX << IMGBITS) | IMGMAX;
       } else {
+        double rng = random->uniform();
+        imol = 0;
+        while (rng > molfrac[imol]) imol++;
         if (dimension == 3) {
           r[0] = random->uniform() - 0.5;
           r[1] = random->uniform() - 0.5;
@@ -491,7 +505,7 @@ void FixPour::pre_exchange()
         MathExtra::axisangle_to_quat(r,theta,quat);
         MathExtra::quat_to_mat(quat,rotmat);
         for (i = 0; i < natom; i++) {
-          MathExtra::matvec(rotmat,onemol->dx[i],coords[i]);
+          MathExtra::matvec(rotmat,onemol[imol].dx[i],coords[i]);
           coords[i][0] += coord[0];
           coords[i][1] += coord[1];
           coords[i][2] += coord[2];
@@ -500,7 +514,7 @@ void FixPour::pre_exchange()
           // default to 0.5, if radii not defined in Molecule
           //   same as atom->avec->create_atom(), invoked below
 
-          if (onemol->radiusflag) coords[i][3] = onemol->radius[i];
+          if (onemol[imol].radiusflag) coords[i][3] = onemol[imol].radius[i];
           else coords[i][3] = 0.5;
 
           imageflags[i] = ((imageint) IMGMAX << IMG2BITS) |
@@ -589,7 +603,7 @@ void FixPour::pre_exchange()
 
       if (flag) {
         if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]);
-        else atom->avec->create_atom(ntype+onemol->type[m],coords[m]);
+        else atom->avec->create_atom(ntype+onemol[imol].type[m],coords[m]);
         int n = atom->nlocal - 1;
         atom->tag[n] = maxtag_all + m+1;
         if (mode == MOLECULE) {
@@ -608,7 +622,7 @@ void FixPour::pre_exchange()
           radtmp = newcoord[3];
           atom->radius[n] = radtmp;
           atom->rmass[n] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
-        } else atom->add_molecule_atom(onemol,m,n,maxtag_all);
+        } else atom->add_molecule_atom(&onemol[imol],m,n,maxtag_all);
         for (j = 0; j < nfix; j++)
           if (fix[j]->create_attribute) fix[j]->set_arrays(n);
       }
@@ -619,9 +633,9 @@ void FixPour::pre_exchange()
     // FixShake::set_molecule stores shake info for molecule
 
     if (rigidflag)
-      fixrigid->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat);
+      fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
     else if (shakeflag)
-      fixshake->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat);
+      fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
 
     maxtag_all += natom;
     if (mode == MOLECULE && atom->molecule_flag) maxmol_all++;
@@ -645,10 +659,10 @@ void FixPour::pre_exchange()
     if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
       error->all(FLERR,"Too many total atoms");
     if (mode == MOLECULE) {
-      atom->nbonds += onemol->nbonds * ninserted_mols;
+      atom->nbonds += onemol[imol].nbonds * ninserted_mols;
-      atom->nangles += onemol->nangles * ninserted_mols;
+      atom->nangles += onemol[imol].nangles * ninserted_mols;
-      atom->ndihedrals += onemol->ndihedrals * ninserted_mols;
+      atom->ndihedrals += onemol[imol].ndihedrals * ninserted_mols;
-      atom->nimpropers += onemol->nimpropers * ninserted_mols;
+      atom->nimpropers += onemol[imol].nimpropers * ninserted_mols;
     }
     if (maxtag_all >= MAXTAGINT)
       error->all(FLERR,"New atom IDs exceed maximum allowed ID");
@@ -696,7 +710,7 @@ void FixPour::find_maxid()
    check if particle i could overlap with a particle inserted into region
    return 1 if yes, 0 if no
    for ATOM mode, use delta with maximum size for inserted atoms
-   for MOLECULE mode, use delta with radius of inserted molecules
+   for MOLECULE mode, use delta with max radius of inserted molecules
    account for PBC in overlap decision via outside() and minimum_image()
 ------------------------------------------------------------------------- */
 
@@ -704,7 +718,7 @@ int FixPour::overlap(int i)
 {
   double delta;
   if (mode == ATOM) delta = atom->radius[i] + radius_max;
-  else delta = atom->radius[i] + onemol->molradius;
+  else delta = atom->radius[i] + molradius_max;
 
   double *x = atom->x[i];
 
@@ -835,6 +849,7 @@ void FixPour::options(int narg, char **arg)
       iregion = domain->find_region(arg[iarg+1]);
       if (iregion == -1) error->all(FLERR,"Fix pour region ID does not exist");
       iarg += 2;
+
     } else if (strcmp(arg[iarg],"mol") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
       int imol = atom->find_molecule(arg[iarg+1]);
@@ -845,7 +860,26 @@ void FixPour::options(int narg, char **arg)
                        "fix pour has multiple molecules");
       mode = MOLECULE;
       onemol = atom->molecules[imol];
+      nmol = onemol->nset;
+      delete [] molfrac;
+      molfrac = new double[nmol];
+      molfrac[0] = 1.0/nmol;
+      for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol;
+      molfrac[nmol-1] = 1.0;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"molfrac") == 0) {
+      if (mode != MOLECULE) error->all(FLERR,"Illegal fix deposit command");
+      if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit command");
+      molfrac[0] = force->numeric(FLERR,arg[iarg+1]);
+      for (int i = 1; i < nmol; i++) 
+        molfrac[i] = molfrac[i-1] + force->numeric(FLERR,arg[iarg+i+1]);
+      if (molfrac[nmol-1] < 1.0-EPSILON || molfrac[nmol-1] > 1.0+EPSILON) 
+        error->all(FLERR,"Illegal fix deposit command");
+      if (nmol > 1) molfrac[nmol-1] = 1.0 - molfrac[nmol-2];
+      else molfrac[nmol-1] = 1.0;
+      iarg += nmol+1;
+
+
     } else if (strcmp(arg[iarg],"rigid") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
       int n = strlen(arg[iarg+1]) + 1;
@@ -960,21 +994,28 @@ void *FixPour::extract(const char *str, int &itype)
     if (mode == ATOM) {
       if (itype == ntype) oneradius = radius_max;
       else oneradius = 0.0;
+
     } else {
-      double *radius = onemol->radius;
-      int *type = onemol->type;
-      int natoms = onemol->natoms;
 
-      // check radii of matching types in Molecule
+      // find a molecule in template with matching type
-      // default to 0.5, if radii not defined in Molecule
-      //   same as atom->avec->create_atom(), invoked in pre_exchange()
 
-      oneradius = 0.0;
+      for (int i = 0; i < nmol; i++) {
-      for (int i = 0; i < natoms; i++)
+        if (itype-ntype > onemol[i].ntypes) continue;
-        if (type[i] == itype-ntype) {
+        double *radius = onemol[i].radius;
-          if (radius) oneradius = MAX(oneradius,radius[i]);
+        int *type = onemol[i].type;
-          else oneradius = 0.5;
+        int natoms = onemol[i].natoms;
-        }
+
+        // check radii of matching types in Molecule
+        // default to 0.5, if radii not defined in Molecule
+        //   same as atom->avec->create_atom(), invoked in pre_exchange()
+
+        oneradius = 0.0;
+        for (int i = 0; i < natoms; i++)
+          if (type[i] == itype-ntype) {
+            if (radius) oneradius = MAX(oneradius,radius[i]);
+            else oneradius = 0.5;
+          }
+      }
     }
     itype = 0;
     return &oneradius;

src/GRANULAR/fix_pour.h

@@ -52,7 +52,9 @@ class FixPour : public Fix {
   char *idrigid,*idshake;
 
   class Molecule *onemol;
-  int natom;                    // # of atoms per inserted particle
+  int natom,nmol;
+  double molradius_max;
+  double *molfrac;
   double **coords;
   imageint *imageflags;
   class Fix *fixrigid,*fixshake;

src/MISC/fix_deposit.cpp

@@ -38,6 +38,8 @@ using namespace MathConst;
 
 enum{ATOM,MOLECULE};
 
+#define EPSILON 1.0e6
+
 /* ---------------------------------------------------------------------- */
 
 FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
@@ -99,21 +101,23 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Cannot use fix_deposit unless atoms have IDs");
 
   if (mode == MOLECULE) {
-    if (onemol->xflag == 0)
+    for (int i = 0; i < nmol; i++) {
-      error->all(FLERR,"Fix deposit molecule must have coordinates");
+      if (onemol[i].xflag == 0)
-    if (onemol->typeflag == 0)
+        error->all(FLERR,"Fix deposit molecule must have coordinates");
-      error->all(FLERR,"Fix deposit molecule must have atom types");
+      if (onemol[i].typeflag == 0)
-    if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes)
+        error->all(FLERR,"Fix deposit molecule must have atom types");
-      error->all(FLERR,"Invalid atom type in fix deposit mol command");
+      if (ntype+onemol[i].ntypes <= 0 || ntype+onemol[i].ntypes > atom->ntypes)
+        error->all(FLERR,"Invalid atom type in fix deposit mol command");
       
-    if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
+      if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
-      error->all(FLERR,"Fix deposit molecule template ID must be same "
+        error->all(FLERR,"Fix deposit molecule template ID must be same "
-                 "as atom_style template ID");
+                   "as atom_style template ID");
-    onemol->check_attributes(0);
+      onemol[i].check_attributes(0);
 
-    // fix deposit uses geoemetric center of molecule for insertion
+      // fix deposit uses geoemetric center of molecule for insertion
       
-    onemol->compute_center();
+      onemol[i].compute_center();
+    }
   }
 
   if (rigidflag && mode == ATOM)
@@ -126,7 +130,11 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
   // setup of coords and imageflags array
 
   if (mode == ATOM) natom = 1;
-  else natom = onemol->natoms;
+  else {
+    natom = 0;
+    for (int i = 0; i < nmol; i++)
+      natom = MAX(natom,onemol[i].natoms);
+  }
   memory->create(coords,natom,3,"deposit:coords");
   memory->create(imageflags,natom,"deposit:imageflags");
 
@@ -184,6 +192,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
 FixDeposit::~FixDeposit()
 {
   delete random;
+  delete [] molfrac;
   delete [] idrigid;
   delete [] idshake;
   delete [] idregion;
@@ -246,7 +255,7 @@ void FixDeposit::init()
 
 void FixDeposit::pre_exchange()
 {
-  int i,j,m,n,nlocalprev,flag,flagall;
+  int i,j,m,n,nlocalprev,imol,flag,flagall;
   double coord[3],lamda[3],delx,dely,delz,rsq;
   double r[3],vnew[3],rotmat[3][3],quat[4];
   double *newcoord;
@@ -351,6 +360,9 @@ void FixDeposit::pre_exchange()
       imageflags[0] = ((imageint) IMGMAX << IMG2BITS) |
         ((imageint) IMGMAX << IMGBITS) | IMGMAX;
     } else {
+      double rng = random->uniform();
+      imol = 0;
+      while (rng > molfrac[imol]) imol++;
       if (dimension == 3) {
         r[0] = random->uniform() - 0.5;
         r[1] = random->uniform() - 0.5;
@@ -364,7 +376,7 @@ void FixDeposit::pre_exchange()
       MathExtra::axisangle_to_quat(r,theta,quat);
       MathExtra::quat_to_mat(quat,rotmat);
       for (i = 0; i < natom; i++) {
-        MathExtra::matvec(rotmat,onemol->dx[i],coords[i]);
+        MathExtra::matvec(rotmat,onemol[imol].dx[i],coords[i]);
         coords[i][0] += coord[0];
         coords[i][1] += coord[1];
         coords[i][2] += coord[2];
@@ -447,7 +459,7 @@ void FixDeposit::pre_exchange()
 
       if (flag) {
         if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]);
-        else atom->avec->create_atom(ntype+onemol->type[m],coords[m]);
+        else atom->avec->create_atom(ntype+onemol[imol].type[m],coords[m]);
         n = atom->nlocal - 1;
         atom->tag[n] = maxtag_all + m+1;
         if (mode == MOLECULE) {
@@ -462,7 +474,8 @@ void FixDeposit::pre_exchange()
         atom->v[n][0] = vnew[0];
         atom->v[n][1] = vnew[1];
         atom->v[n][2] = vnew[2];
-        if (mode == MOLECULE) atom->add_molecule_atom(onemol,m,n,maxtag_all);
+        if (mode == MOLECULE) 
+          atom->add_molecule_atom(&onemol[imol],m,n,maxtag_all);
         for (j = 0; j < nfix; j++)
           if (fix[j]->create_attribute) fix[j]->set_arrays(n);
       }
@@ -473,9 +486,9 @@ void FixDeposit::pre_exchange()
     // FixShake::set_molecule stores shake info for molecule
 
     if (rigidflag)
-      fixrigid->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat);
+      fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
     else if (shakeflag)
-      fixshake->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat);
+      fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
 
     // old code: unsuccessful if no proc performed insertion of an atom
     // don't think that check is necessary
@@ -504,10 +517,10 @@ void FixDeposit::pre_exchange()
     if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
       error->all(FLERR,"Too many total atoms");
     if (mode == MOLECULE) {
-      atom->nbonds += onemol->nbonds;
+      atom->nbonds += onemol[imol].nbonds;
-      atom->nangles += onemol->nangles;
+      atom->nangles += onemol[imol].nangles;
-      atom->ndihedrals += onemol->ndihedrals;
+      atom->ndihedrals += onemol[imol].ndihedrals;
-      atom->nimpropers += onemol->nimpropers;
+      atom->nimpropers += onemol[imol].nimpropers;
     }
     maxtag_all += natom;
     if (maxtag_all >= MAXTAGINT)
@@ -561,6 +574,7 @@ void FixDeposit::options(int narg, char **arg)
   iregion = -1;
   idregion = NULL;
   mode = ATOM;
+  molfrac = NULL;
   rigidflag = 0;
   idrigid = NULL;
   shakeflag = 0;
@@ -586,17 +600,33 @@ void FixDeposit::options(int narg, char **arg)
       idregion = new char[n];
       strcpy(idregion,arg[iarg+1]);
       iarg += 2;
+
     } else if (strcmp(arg[iarg],"mol") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
       int imol = atom->find_molecule(arg[iarg+1]);
       if (imol == -1)
         error->all(FLERR,"Molecule template ID for fix deposit does not exist");
-      if (atom->molecules[imol]->nset > 1 && comm->me == 0)
-        error->warning(FLERR,"Molecule template for "
-                       "fix deposit has multiple molecules");
       mode = MOLECULE;
       onemol = atom->molecules[imol];
+      nmol = onemol->nset;
+      delete [] molfrac;
+      molfrac = new double[nmol];
+      molfrac[0] = 1.0/nmol;
+      for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol;
+      molfrac[nmol-1] = 1.0;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"molfrac") == 0) {
+      if (mode != MOLECULE) error->all(FLERR,"Illegal fix deposit command");
+      if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit command");
+      molfrac[0] = force->numeric(FLERR,arg[iarg+1]);
+      for (int i = 1; i < nmol; i++) 
+        molfrac[i] = molfrac[i-1] + force->numeric(FLERR,arg[iarg+i+1]);
+      if (molfrac[nmol-1] < 1.0-EPSILON || molfrac[nmol-1] > 1.0+EPSILON) 
+        error->all(FLERR,"Illegal fix deposit command");
+      if (nmol > 1) molfrac[nmol-1] = 1.0 - molfrac[nmol-2];
+      else molfrac[nmol-1] = 1.0;
+      iarg += nmol+1;
+
     } else if (strcmp(arg[iarg],"rigid") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
       int n = strlen(arg[iarg+1]) + 1;
@@ -730,24 +760,32 @@ void *FixDeposit::extract(const char *str, int &itype)
     if (mode == ATOM) {
       if (itype == ntype) oneradius = 0.5;
       else oneradius = 0.0;
+
     } else {
-      double *radius = onemol->radius;
-      int *type = onemol->type;
-      int natoms = onemol->natoms;
 
-      // check radii of matching types in Molecule
+      // find a molecule in template with matching type
-      // default to 0.5, if radii not defined in Molecule
-      //   same as atom->avec->create_atom(), invoked in pre_exchange()
 
-      oneradius = 0.0;
+      for (int i = 0; i < nmol; i++) {
-      for (int i = 0; i < natoms; i++)
+        if (itype-ntype > onemol[i].ntypes) continue;
-        if (type[i] == itype-ntype) {
+        double *radius = onemol[i].radius;
-          if (radius) oneradius = MAX(oneradius,radius[i]);
+        int *type = onemol[i].type;
-          else oneradius = 0.5;
+        int natoms = onemol[i].natoms;
-        }
+
+        // check radii of matching types in Molecule
+        // default to 0.5, if radii not defined in Molecule
+        //   same as atom->avec->create_atom(), invoked in pre_exchange()
+
+        oneradius = 0.0;
+        for (int i = 0; i < natoms; i++)
+          if (type[i] == itype-ntype) {
+            if (radius) oneradius = MAX(oneradius,radius[i]);
+            else oneradius = 0.5;
+          }
+      }
     }
     itype = 0;
     return &oneradius;
   }
+
   return NULL;
 }

src/MISC/fix_deposit.h

@@ -48,7 +48,8 @@ class FixDeposit : public Fix {
   char *idrigid,*idshake;
 
   class Molecule *onemol;
-  int natom;
+  int natom,nmol;
+  double *molfrac;
   double **coords;
   imageint *imageflags;
   class Fix *fixrigid,*fixshake;

src/RIGID/fix_rigid_small.cpp

@@ -178,6 +178,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
         error->warning(FLERR,"Molecule template for "
                        "fix rigid/small has multiple molecules");
       onemol = atom->molecules[imol];
+      nmol = onemol->nset;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"temp") == 0) {
@@ -300,17 +301,19 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   // error check and further setup for Molecule template
 
   if (onemol) {
-    if (onemol->xflag == 0)
+    for (int i = 0; i < nmol; i++) {
-      error->all(FLERR,"Fix rigid/small molecule must have coordinates");
+      if (onemol[i].xflag == 0)
-    if (onemol->typeflag == 0)
+        error->all(FLERR,"Fix rigid/small molecule must have coordinates");
-      error->all(FLERR,"Fix rigid/small molecule must have atom types");
+      if (onemol[i].typeflag == 0)
+        error->all(FLERR,"Fix rigid/small molecule must have atom types");
 
-    // fix rigid/small uses center, masstotal, COM, inertia of molecule
+      // fix rigid/small uses center, masstotal, COM, inertia of molecule
       
-    onemol->compute_center();
+      onemol[i].compute_center();
-    onemol->compute_mass();
+      onemol[i].compute_mass();
-    onemol->compute_com();
+      onemol[i].compute_com();
-    onemol->compute_inertia();
+      onemol[i].compute_inertia();
+    }
   }
 
   // set pstat_flag
@@ -1549,7 +1552,10 @@ void FixRigidSmall::create_bodies()
 
   MPI_Allreduce(&rsqfar,&maxextent,1,MPI_DOUBLE,MPI_MAX,world);
   maxextent = sqrt(maxextent);
-  if (onemol) maxextent = MAX(maxextent,onemol->maxextent);
+  if (onemol) {
+    for (int i = 0; i < nmol; i++)
+      maxextent = MAX(maxextent,onemol[i].maxextent);
+  }
 
   // clean up
 
@@ -1716,8 +1722,17 @@ void FixRigidSmall::setup_bodies_static()
   }
 
   // extended = 1 if using molecule template with finite-size particles
+  // require all molecules in template to have consistent radiusflag
 
-  if (onemol && onemol->radiusflag) extended = 1;
+  if (onemol) {
+    int radiusflag = onemol->radiusflag;
+    for (i = 1; i < nmol; i++) {
+      if (onemol->radiusflag != radiusflag)
+        error->all(FLERR,"Inconsistent use of finite-size particles "
+                   "by molecule template molecules");
+    }
+    if (radiusflag) extended = 1;
+  }
 
   // grow extended arrays and set extended flags for each particle
   // orientflag = 4 if any particle stores ellipsoid or tri orientation
@@ -2579,7 +2594,7 @@ void FixRigidSmall::set_arrays(int i)
           relative to template in Molecule class
 ------------------------------------------------------------------------- */
 
-void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, 
+void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int imol,
                                  double *xgeom, double *vcm, double *quat)
 {
   int m;
@@ -2592,16 +2607,16 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev,
   tagint *tag = atom->tag;
 
   for (int i = nlocalprev; i < nlocal; i++) {
-    bodytag[i] = tagprev + onemol->comatom;
+    bodytag[i] = tagprev + onemol[imol].comatom;
-    if (tag[i]-tagprev == onemol->comatom) bodyown[i] = nlocal_body;
+    if (tag[i]-tagprev == onemol[imol].comatom) bodyown[i] = nlocal_body;
 
     m = tag[i] - tagprev-1;
-    displace[i][0] = onemol->dxbody[m][0];
+    displace[i][0] = onemol[imol].dxbody[m][0];
-    displace[i][1] = onemol->dxbody[m][1];
+    displace[i][1] = onemol[imol].dxbody[m][1];
-    displace[i][2] = onemol->dxbody[m][2];
+    displace[i][2] = onemol[imol].dxbody[m][2];
 
     eflags[i] = 0;
-    if (onemol->radiusflag) {
+    if (onemol[imol].radiusflag) {
       eflags[i] |= SPHERE;
       eflags[i] |= OMEGA;
       eflags[i] |= TORQUE;
@@ -2610,27 +2625,27 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev,
     if (bodyown[i] >= 0) {
       if (nlocal_body == nmax_body) grow_body();
       Body *b = &body[nlocal_body];
-      b->mass = onemol->masstotal;
+      b->mass = onemol[imol].masstotal;
 
-      // new COM = Q (onemol->xcm - onemol->center) + xgeom
+      // new COM = Q (onemol[imol].xcm - onemol[imol].center) + xgeom
       // Q = rotation matrix associated with quat
 
       MathExtra::quat_to_mat(quat,rotmat);
-      MathExtra::sub3(onemol->com,onemol->center,ctr2com);
+      MathExtra::sub3(onemol[imol].com,onemol[imol].center,ctr2com);
       MathExtra::matvec(rotmat,ctr2com,ctr2com_rotate);
       MathExtra::add3(ctr2com_rotate,xgeom,b->xcm);
 
       b->vcm[0] = vcm[0];
       b->vcm[1] = vcm[1];
       b->vcm[2] = vcm[2];
-      b->inertia[0] = onemol->inertia[0];
+      b->inertia[0] = onemol[imol].inertia[0];
-      b->inertia[1] = onemol->inertia[1];
+      b->inertia[1] = onemol[imol].inertia[1];
-      b->inertia[2] = onemol->inertia[2];
+      b->inertia[2] = onemol[imol].inertia[2];
 
       // final quat is product of insertion quat and original quat
       // true even if insertion rotation was not around COM
 
-      MathExtra::quatquat(quat,onemol->quat,b->quat);
+      MathExtra::quatquat(quat,onemol[imol].quat,b->quat);
       MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
 
       b->angmom[0] = b->angmom[1] = b->angmom[2] = 0.0;

src/RIGID/fix_rigid_small.h

@@ -48,7 +48,7 @@ class FixRigidSmall : public Fix {
   void grow_arrays(int);
   void copy_arrays(int, int, int);
   void set_arrays(int);
-  void set_molecule(int, tagint, double *, double *, double *);
+  void set_molecule(int, tagint, int, double *, double *, double *);
 
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
@@ -171,6 +171,7 @@ class FixRigidSmall : public Fix {
   // molecules added on-the-fly as rigid bodies
 
   class Molecule *onemol;
+  int nmol;
 
   // class data used by ring communication callbacks
 

src/RIGID/fix_shake.cpp

@@ -161,14 +161,18 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
         error->warning(FLERR,"Molecule template for "
                        "fix shake has multiple molecules");
       onemol = atom->molecules[imol];
+      nmol = onemol->nset;
       iarg += 2;
     } else error->all(FLERR,"Illegal fix shake command");
   }
 
   // error check for Molecule template
 
-  if (onemol && onemol->shakeflag == 0)
+  if (onemol) {
-    error->all(FLERR,"Fix shake molecule template must have shake info");
+    for (int i = 0; i < nmol; i++)
+      if (onemol[i].shakeflag == 0)
+        error->all(FLERR,"Fix shake molecule template must have shake info");
+  }
 
   // allocate bond and angle distance arrays, indexed from 1 to n
 
@@ -2456,7 +2460,7 @@ void FixShake::update_arrays(int i, int atom_offset)
    xgeom,vcm,quat ignored
 ------------------------------------------------------------------------- */
 
-void FixShake::set_molecule(int nlocalprev, tagint tagprev, 
+void FixShake::set_molecule(int nlocalprev, tagint tagprev, int imol,
                             double *xgeom, double *vcm, double *quat)
 {
   int m,flag;
@@ -2465,13 +2469,13 @@ void FixShake::set_molecule(int nlocalprev, tagint tagprev,
   if (nlocalprev == nlocal) return;
 
   tagint *tag = atom->tag;
-  tagint **mol_shake_atom = onemol->shake_atom;
+  tagint **mol_shake_atom = onemol[imol].shake_atom;
-  int **mol_shake_type = onemol->shake_type;
+  int **mol_shake_type = onemol[imol].shake_type;
 
   for (int i = nlocalprev; i < nlocal; i++) {
     m = tag[i] - tagprev-1;
 
-    flag = shake_flag[i] = onemol->shake_flag[m];
+    flag = shake_flag[i] = onemol[imol].shake_flag[m];
 
     if (flag == 1) {
       shake_atom[i][0] = mol_shake_atom[m][0] + tagprev;

src/RIGID/fix_shake.h

@@ -40,7 +40,7 @@ class FixShake : public Fix {
   void copy_arrays(int, int, int);
   void set_arrays(int);
   void update_arrays(int, int);
-  void set_molecule(int, tagint, double *, double *, double *);
+  void set_molecule(int, tagint, int, double *, double *, double *);
 
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
@@ -105,6 +105,7 @@ class FixShake : public Fix {
 
   class Molecule **onemols;             // atom style template pointer
   class Molecule *onemol;               // molecule added on-the-fly
+  int nmol;
 
   void find_clusters();
   int masscheck(double);

src/fix.h

@@ -125,7 +125,7 @@ class Fix : protected Pointers {
   virtual void copy_arrays(int, int, int) {}
   virtual void set_arrays(int) {}
   virtual void update_arrays(int, int) {}
-  virtual void set_molecule(int, tagint, double *, double *, double *) {}
+  virtual void set_molecule(int, tagint, int, double *, double *, double *) {}
 
   virtual int pack_border(int, int *, double *) {return 0;}
   virtual int unpack_border(int, int, double *) {return 0;}