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

src/atom.cpp

@@ -31,14 +31,18 @@
 #include "update.h"
 #include "domain.h"
 #include "group.h"
+#include "molecule.h"
 #include "accelerator_cuda.h"
 #include "atom_masks.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define DELTA 1
+#define DELTA_MOLECULE 1
 #define DELTA_MEMSTR 1024
 #define EPSILON 1.0e-6
 #define CUDA_CHUNK 3000
@@ -106,6 +110,11 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   improper_type = improper_atom1 = improper_atom2 = NULL;
   improper_atom3 = improper_atom4 = NULL;
 
+  // user-defined molecules
+
+  nmolecule = maxmol = 0;
+  molecules = NULL;
+
   // custom atom arrays
 
   nivector = ndvector = 0;
@@ -230,6 +239,11 @@ Atom::~Atom()
   memory->destroy(improper_atom3);
   memory->destroy(improper_atom4);
 
+  // delete user-defined molecules
+
+  for (int i = 0; i < nmolecule; i++) delete molecules[i];
+  memory->sfree(molecules);
+
   // delete custom atom arrays
 
   for (int i = 0; i < nivector; i++) {
@@ -1152,6 +1166,99 @@ int Atom::shape_consistency(int itype,
   return 1;
 }
 
+/* ----------------------------------------------------------------------
+   add a new molecule template
+------------------------------------------------------------------------- */
+
+void Atom::add_molecule(int narg, char **arg)
+{
+  if (narg < 1) error->all(FLERR,"Illegal molecule command");
+  if (find_molecule(arg[0]) >= 0) error->all(FLERR,"Reuse of molecule ID");
+
+  // extend molecule list if necessary
+
+  if (nmolecule == maxmol) {
+    maxmol += DELTA_MOLECULE;
+    molecules = (Molecule **)
+      memory->srealloc(molecules,maxmol*sizeof(Molecule *),"atom::molecules");
+  }
+
+  molecules[nmolecule++] = new Molecule(lmp,narg,arg);
+}
+
+/* ----------------------------------------------------------------------
+   find which molecule has molecule ID
+   return -1 if does not exist
+------------------------------------------------------------------------- */
+
+int Atom::find_molecule(char *id)
+{
+  int imol;
+  for (imol = 0; imol < nmolecule; imol++)
+    if (strcmp(id,molecules[imol]->id) == 0) return imol;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   add info for iatom from molecule template onemol to current atom ilocal
+------------------------------------------------------------------------- */
+
+void Atom::add_molecule_atom(Molecule *onemol, int iatom,
+                             int ilocal, int offset)
+{
+  if (onemol->qflag) q[ilocal] = onemol->q[iatom];
+  if (onemol->radiusflag) radius[ilocal] = onemol->radius[iatom];
+  if (onemol->rmassflag) rmass[ilocal] = onemol->rmass[iatom];
+  else if (rmass_flag) 
+    rmass[ilocal] = 4.0*MY_PI/3.0 *
+      radius[ilocal]*radius[ilocal]*radius[ilocal];
+
+  if (onemol->bondflag) {
+    num_bond[ilocal] = onemol->num_bond[iatom];
+    for (int i = 0; i < num_bond[ilocal]; i++) {
+      bond_type[ilocal][i] = onemol->bond_type[iatom][i];
+      bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
+    }
+  }
+  if (onemol->angleflag) {
+    num_angle[ilocal] = onemol->num_angle[iatom];
+    for (int i = 0; i < num_angle[ilocal]; i++) {
+      angle_type[ilocal][i] = onemol->angle_type[iatom][i];
+      angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
+      angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
+      angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
+    }
+  }
+  if (onemol->dihedralflag) {
+    num_dihedral[ilocal] = onemol->num_dihedral[iatom];
+    for (int i = 0; i < num_dihedral[ilocal]; i++) {
+      dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
+      dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset;
+      dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset;
+      dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
+      dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
+    }
+  }
+  if (onemol->improperflag) {
+    num_improper[ilocal] = onemol->num_improper[iatom];
+    for (int i = 0; i < num_improper[ilocal]; i++) {
+      improper_type[ilocal][i] = onemol->improper_type[iatom][i];
+      improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset;
+      improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset;
+      improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
+      improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
+    }
+  }
+
+  if (onemol->specialflag) {
+    nspecial[ilocal][0] = onemol->nspecial[iatom][0];
+    nspecial[ilocal][1] = onemol->nspecial[iatom][1];
+    int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];
+    for (int i = 0; i < n; i++)
+      special[ilocal][i] = onemol->special[iatom][i] + offset;
+  }
+}
+
 /* ----------------------------------------------------------------------
    reorder owned atoms so those in firstgroup appear first
    called by comm->exchange() if atom_modify first group is set

src/atom.h

@@ -109,6 +109,11 @@ class Atom : protected Pointers {
   int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag;
   int rho_flag,e_flag,cv_flag,vest_flag;
 
+  // molecules
+
+  int nmolecule,maxmol;
+  class Molecule **molecules;
+
   // extra peratom info in restart file destined for fix & diag
 
   double **extra;
@@ -178,6 +183,10 @@ class Atom : protected Pointers {
   int radius_consistency(int, double &);
   int shape_consistency(int, double &, double &, double &);
 
+  void add_molecule(int, char **);
+  int find_molecule(char *);
+  void add_molecule_atom(class Molecule *, int, int, int);
+
   void first_reorder();
   void sort();
 

src/fix_deposit.cpp

@@ -17,6 +17,7 @@
 #include "fix_deposit.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "molecule.h"
 #include "force.h"
 #include "update.h"
 #include "modify.h"
@@ -26,11 +27,16 @@
 #include "lattice.h"
 #include "region.h"
 #include "random_park.h"
+#include "math_extra.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
+using namespace MathConst;
+
+enum{ATOM,MOLECULE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -55,6 +61,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
 
   iregion = -1;
   idregion = NULL;
+  mode = ATOM;
   idnext = 0;
   globalflag = localflag = 0;
   lo = hi = deltasq = 0.0;
@@ -96,6 +103,24 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
       error->all(FLERR,"Deposition region extends outside simulation box");
   }
 
+  // error check and further setup for mode = MOLECULE
+
+  if (mode == MOLECULE) {
+    if (onemol->xflag == 0)
+      error->all(FLERR,"Fix deposit molecule must have coordinates");
+    if (onemol->typeflag == 0)
+      error->all(FLERR,"Fix deposit molecule must have atom types");
+
+    onemol->compute_center();
+  }
+
+  // setup of coords and imagesflags array
+
+  if (mode == ATOM) natom = 1;
+  else natom = onemol->natoms;
+  memory->create(coords,natom,3,"deposit:coords");
+  memory->create(imageflags,natom,"deposit:imageflags");
+
   // setup scaling
 
   double xscale,yscale,zscale;
@@ -129,16 +154,9 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
   ty *= yscale;
   tz *= zscale;
 
-  // maxtag_all = current max tag for all atoms
+  // find current max atom and molecule IDs if necessary
 
-  if (idnext) {
+  if (idnext) find_maxid();
-    int *tag = atom->tag;
-    int nlocal = atom->nlocal;
-
-    int maxtag = 0;
-    for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
-    MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world);
-  }
 
   // random number generator, same for all procs
 
@@ -158,6 +176,8 @@ FixDeposit::~FixDeposit()
 {
   delete random;
   delete [] idregion;
+  memory->destroy(coords);
+  memory->destroy(imageflags);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -186,15 +206,19 @@ void FixDeposit::init()
 
 void FixDeposit::pre_exchange()
 {
-  int i,j;
+  int i,j,m;
   int flag,flagall;
   double coord[3],lamda[3],delx,dely,delz,rsq;
-  double *newcoord;
+  double alpha,beta,gamma;
+  double r[3],rotmat[3][3],quat[4];
+  double *newcoord,*center;
 
   // just return if should not be called on this timestep
 
   if (next_reneighbor != update->ntimestep) return;
 
+  int dimension = domain->dimension;
+
   // compute current offset = bottom of insertion volume
 
   double offset = 0.0;
@@ -209,6 +233,10 @@ void FixDeposit::pre_exchange()
     subhi = domain->subhi_lamda;
   }
 
+  // find current max atom and molecule IDs if necessary
+
+  if (!idnext) find_maxid();
+
   // attempt an insertion until successful
 
   int nfix = modify->nfix;
@@ -219,7 +247,7 @@ void FixDeposit::pre_exchange()
   while (attempt < maxattempt) {
     attempt++;
 
-    // choose random position for new atom within region
+    // choose random position for new particle within region
 
     coord[0] = xlo + random->uniform() * (xhi-xlo);
     coord[1] = ylo + random->uniform() * (yhi-ylo);
@@ -232,18 +260,19 @@ void FixDeposit::pre_exchange()
 
     // adjust vertical coord by offset
 
-    if (domain->dimension == 2) coord[1] += offset;
+    if (dimension == 2) coord[1] += offset;
     else coord[2] += offset;
 
     // if global, reset vertical coord to be lo-hi above highest atom
     // if local, reset vertical coord to be lo-hi above highest "nearby" atom
     // local computation computes lateral distance between 2 particles w/ PBC
+    // when done, have final coord of atom or center pt of molecule
 
     if (globalflag || localflag) {
       int dim;
       double max,maxall,delx,dely,delz,rsq;
 
-      if (domain->dimension == 2) {
+      if (dimension == 2) {
         dim = 1;
         max = domain->boxlo[1];
       } else {
@@ -259,7 +288,7 @@ void FixDeposit::pre_exchange()
           dely = coord[1] - x[i][1];
           delz = 0.0;
           domain->minimum_image(delx,dely,delz);
-          if (domain->dimension == 2) rsq = delx*delx;
+          if (dimension == 2) rsq = delx*delx;
           else rsq = delx*delx + dely*dely;
           if (rsq > deltasq) continue;
         }
@@ -267,38 +296,76 @@ void FixDeposit::pre_exchange()
       }
 
       MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
-      if (domain->dimension == 2)
+      if (dimension == 2)
         coord[1] = maxall + lo + random->uniform()*(hi-lo);
       else
         coord[2] = maxall + lo + random->uniform()*(hi-lo);
     }
 
-    // now have final coord
+    // coords = coords of all atoms
-    // if distance to any atom is less than near, try again
+    // for molecule, perform random rotation around center pt
+    // apply PBC so final coords are inside box
+    // also store image flag modified due to PBC
+
+    if (mode == ATOM) {
+      coords[0][0] = coord[0];
+      coords[0][1] = coord[1];
+      coords[0][2] = coord[2];
+    } else {
+      if (dimension == 3) {
+        r[0] = random->uniform() - 0.5;
+        r[1] = random->uniform() - 0.5;
+        r[2] = random->uniform() - 0.5;
+      } else {
+        r[0] = r[1] = 0.0;
+        r[2] = 1.0;
+      }
+      double theta = random->uniform() * MY_2PI;
+      MathExtra::norm3(r);
+      MathExtra::axisangle_to_quat(r,theta,quat);
+      MathExtra::quat_to_mat(quat,rotmat);
+      center = onemol->center;
+      for (i = 0; i < natom; i++) {
+        MathExtra::matvec(rotmat,onemol->dx[i],coords[i]);
+        coords[i][0] += center[0] + coord[0];
+        coords[i][1] += center[1] + coord[1];
+        coords[i][2] += center[2] + coord[2];
+
+        imageflags[i] = ((tagint) IMGMAX << IMG2BITS) |
+          ((tagint) IMGMAX << IMGBITS) | IMGMAX;
+        domain->remap(coords[i],imageflags[i]);
+      }
+    }
+
+    // if distance to any inserted atom is less than near, try again
 
     double **x = atom->x;
     int nlocal = atom->nlocal;
 
     flag = 0;
-    for (i = 0; i < nlocal; i++) {
+    for (m = 0; m < natom; m++) {
-      delx = coord[0] - x[i][0];
+      for (i = 0; i < nlocal; i++) {
-      dely = coord[1] - x[i][1];
+        delx = coords[m][0] - x[i][0];
-      delz = coord[2] - x[i][2];
+        dely = coords[m][1] - x[i][1];
-      domain->minimum_image(delx,dely,delz);
+        delz = coords[m][2] - x[i][2];
-      rsq = delx*delx + dely*dely + delz*delz;
+        domain->minimum_image(delx,dely,delz);
-      if (rsq < nearsq) flag = 1;
+        rsq = delx*delx + dely*dely + delz*delz;
+        if (rsq < nearsq) flag = 1;
+      }
     }
     MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
     if (flagall) continue;
 
-    // insertion will proceed
+    // proceed with insertion
-    // choose random velocity for new atom
+    // choose random velocity for new particle
+    // used for every atom in molecule
 
     double vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
     double vytmp = vylo + random->uniform() * (vyhi-vylo);
     double vztmp = vzlo + random->uniform() * (vzhi-vzlo);
 
-    // if we have a sputter target change velocity vector accordingly
+    // if target specified, change velocity vector accordingly
+
     if (targetflag) {
       double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp);
       delx = tx - coord[0];
@@ -313,67 +380,84 @@ void FixDeposit::pre_exchange()
       }
     }
 
-    // check if new atom is in my sub-box or above it if I'm highest proc
+    // check if new atoms are in my sub-box or above it if I am highest proc
-    // if so, add to my list via create_atom()
+    // if so, add atom to my list via create_atom()
-    // initialize info about the atoms
+    // initialize additional info about the atoms
     // set group mask to "all" plus fix group
 
-    if (domain->triclinic) {
+    for (m = 0; m < natom; m++) {
-      domain->x2lamda(coord,lamda);
+      if (domain->triclinic) {
-      newcoord = lamda;
+        domain->x2lamda(coords[m],lamda);
-    } else newcoord = coord;
+        newcoord = lamda;
+      } else newcoord = coords[m];
+      
+      flag = 0;
+      if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+          newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
+          newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
+      else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] &&
+               comm->myloc[2] == comm->procgrid[2]-1 &&
+               newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+               newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
+      else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
+               comm->myloc[1] == comm->procgrid[1]-1 &&
+               newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
 
-    flag = 0;
+      if (flag) {
-    if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+        if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]);
-        newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
+        else atom->avec->create_atom(onemol->type[m],coords[m]);
-        newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
+        int n = atom->nlocal - 1;
-    else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] &&
+        atom->tag[n] = maxtag_all + m+1;
-             comm->myloc[2] == comm->procgrid[2]-1 &&
+        if (atom->molecule_flag) atom->molecule[n] = maxmol_all;
-             newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
+        atom->mask[n] = 1 | groupbit;
-             newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
+        atom->image[n] = imageflags[m];
-    else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
+        atom->v[n][0] = vxtmp;
-             comm->myloc[1] == comm->procgrid[1]-1 &&
+        atom->v[n][1] = vytmp;
-             newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
+        atom->v[n][2] = vztmp;
-
+        if (mode == MOLECULE) atom->add_molecule_atom(onemol,m,n,maxtag_all);
-    if (flag) {
+        for (j = 0; j < nfix; j++)
-      atom->avec->create_atom(ntype,coord);
+          if (fix[j]->create_attribute) fix[j]->set_arrays(n);
-      int m = atom->nlocal - 1;
+      }
-      atom->type[m] = ntype;
-      atom->mask[m] = 1 | groupbit;
-      atom->v[m][0] = vxtmp;
-      atom->v[m][1] = vytmp;
-      atom->v[m][2] = vztmp;
-      for (j = 0; j < nfix; j++)
-        if (fix[j]->create_attribute) fix[j]->set_arrays(m);
     }
-    MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
+
+    // old code: unsuccessful if no proc performed insertion of an atom
+    // don't think that check is necessary
+    // if get this far, should always be succesful
+    // would be hard to undo partial insertion for a molecule
+    // better to check how many atoms could be inserted (w/out inserting)
+    //   then sum to insure all are inserted, before doing actual insertion
+    // MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
+
+    success = 1;
     break;
   }
 
-  // warn if not successful b/c too many attempts or no proc owned particle
+  // warn if not successful b/c too many attempts
 
   if (!success && comm->me == 0)
     error->warning(FLERR,"Particle deposition was unsuccessful",0);
 
-  // reset global natoms
+  // reset global natoms,nbonds,etc
-  // if idnext, set new atom ID to incremented maxtag_all
+  // increment maxtag_all and maxmol_all if necessary
-  // else set new atom ID to value beyond all current atoms
   // if global map exists, reset it now instead of waiting for comm
   // since adding an atom messes up ghosts
 
   if (success) {
-    atom->natoms += 1;
+    atom->natoms += natom;
-    if (atom->tag_enable) {
+    if (mode == MOLECULE) {
-      if (idnext) {
+      atom->nbonds += onemol->nbonds;
-        maxtag_all++;
+      atom->nangles += onemol->nangles;
-        if (atom->nlocal && atom->tag[atom->nlocal-1] == 0) 
+      atom->ndihedrals += onemol->ndihedrals;
-          atom->tag[atom->nlocal-1] = maxtag_all;
+      atom->nimpropers += onemol->nimpropers;
-      } else atom->tag_extend();
+    }
-      if (atom->map_style) {
+    if (idnext) {
-        atom->nghost = 0;
+      maxtag_all += natom;
-        atom->map_init();
+      if (mode == MOLECULE) maxmol_all++;
-        atom->map_set();
+    }
-      }
+    if (atom->map_style) {
+      atom->nghost = 0;
+      atom->map_init();
+      atom->map_set();
     }
   }
 
@@ -385,6 +469,28 @@ void FixDeposit::pre_exchange()
   else next_reneighbor = 0;
 }
 
+/* ----------------------------------------------------------------------
+   maxtag_all = current max atom ID for all atoms
+   maxmol_all = current max molecule ID for all atoms
+------------------------------------------------------------------------- */
+
+void FixDeposit::find_maxid()
+{
+  int *tag = atom->tag;
+  int *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+
+  int max = 0;
+  for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
+  MPI_Allreduce(&max,&maxtag_all,1,MPI_INT,MPI_MAX,world);
+
+  if (mode == MOLECULE) {
+    max = 0;
+    for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
+    MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world);
+  }
+}
+
 /* ----------------------------------------------------------------------
    parse optional parameters at end of input line
 ------------------------------------------------------------------------- */
@@ -404,6 +510,14 @@ 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 ID for fix deposit does not exist");
+      mode = MOLECULE;
+      onemol = atom->molecules[imol];
+     iarg += 2;
     } else if (strcmp(arg[iarg],"id") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
       if (strcmp(arg[iarg+1],"max") == 0) idnext = 0;
@@ -423,11 +537,13 @@ void FixDeposit::options(int narg, char **arg)
       globalflag = 0;
       lo = force->numeric(FLERR,arg[iarg+1]);
       hi = force->numeric(FLERR,arg[iarg+2]);
-      deltasq = force->numeric(FLERR,arg[iarg+3])*force->numeric(FLERR,arg[iarg+3]);
+      deltasq = force->numeric(FLERR,arg[iarg+3]) * 
+        force->numeric(FLERR,arg[iarg+3]);
       iarg += 4;
     } else if (strcmp(arg[iarg],"near") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
-      nearsq = force->numeric(FLERR,arg[iarg+1])*force->numeric(FLERR,arg[iarg+1]);
+      nearsq = force->numeric(FLERR,arg[iarg+1]) * 
+        force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg],"attempt") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");

src/fix_deposit.h

@@ -40,15 +40,24 @@ class FixDeposit : public Fix {
  private:
   int ninsert,nfreq,seed;
   int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag;
-  char *idregion;
+  int mode;
   double lo,hi,deltasq,nearsq,rate;
   double vxlo,vxhi,vylo,vyhi,vzlo,vzhi;
   double xlo,xhi,ylo,yhi,zlo,zhi;
   double tx,ty,tz;
+  char *idregion;
+  class Molecule *onemol;
+
+  int natom;
+  double **coords;
+  int *imageflags;
+
   int nfirst,ninserted;
-  int idnext,maxtag_all;
+  int idnext;
+  int maxtag_all,maxmol_all;
   class RanPark *random;
 
+  void find_maxid();
   void options(int, char **);
 };
 

src/input.cpp

@@ -564,6 +564,7 @@ int Input::execute_command()
   else if (!strcmp(command,"mass")) mass();
   else if (!strcmp(command,"min_modify")) min_modify();
   else if (!strcmp(command,"min_style")) min_style();
+  else if (!strcmp(command,"molecule")) molecule();
   else if (!strcmp(command,"neigh_modify")) neigh_modify();
   else if (!strcmp(command,"neighbor")) neighbor_command();
   else if (!strcmp(command,"newton")) newton();
@@ -1311,6 +1312,13 @@ void Input::min_style()
 
 /* ---------------------------------------------------------------------- */
 
+void Input::molecule()
+{
+  atom->add_molecule(narg,arg);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Input::neigh_modify()
 {
   neighbor->modify_params(narg,arg);

src/input.h

@@ -103,6 +103,7 @@ class Input : protected Pointers {
   void mass();
   void min_modify();
   void min_style();
+  void molecule();
   void neigh_modify();
   void neighbor_command();
   void newton();

src/molecule.h

@@ -0,0 +1,118 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_ONE_MOLECULE_H
+#define LMP_ONE_MOLEUCULE_H
+
+#include "pointers.h"
+
+namespace LAMMPS_NS {
+
+class Molecule : protected Pointers {
+ public:
+  char *id;
+
+  // number of atoms,bonds,etc in molecule
+
+  int natoms;
+  int nbonds,nangles,ndihedrals,nimpropers;
+
+  // max bond,angle,etc per atom
+
+  int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
+  int maxspecial;
+
+  // 1 if attribute defined in file, 0 if not
+
+  int xflag,typeflag,qflag,radiusflag,rmassflag;
+  int bondflag,angleflag,dihedralflag,improperflag;
+  int nspecialflag,specialflag;
+
+  // 1 if attribute defined or computed, 0 if not
+
+  int centerflag,comflag,inertiaflag;
+
+  // attributes
+
+  double **x;          // displacement of each atom from origin
+  int *type;           // type of each atom
+  double *q;           // charge on each atom
+  double *radius;      // radius of each atom
+  double *rmass;       // mass of each atom
+
+  int *num_bond;       // bonds, angles, dihedrals, impropers for each atom
+  int **bond_type;
+  int **bond_atom;
+  
+  int *num_angle;
+  int **angle_type;
+  int **angle_atom1,**angle_atom2,**angle_atom3;
+  
+  int *num_dihedral;
+  int **dihedral_type;
+  int **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
+  
+  int *num_improper;
+  int **improper_type;
+  int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
+
+  int **nspecial;
+  int **special;
+
+  double center[3];    // geometric center of molecule
+  double xcm[3];       // center of mass of molecule
+  double inertia[6];   // moments of inertia of molecule
+
+  double **dx;         // displacement of each atom relative to center
+  double **dxcom;      // displacement of each atom relative to COM
+
+  Molecule(class LAMMPS *, int, char **);
+  ~Molecule();
+  void compute_center();
+  void compute_xcm();
+  void compute_inertia();
+
+ private:
+  int me;
+  FILE *fp;
+  int *count;
+
+  void read(int);
+  void coords(char *);
+  void types(char *);
+  void charges(char *);
+  void diameters(char *);
+  void masses(char *);
+  void bonds(int, char *);
+  void angles(int, char *);
+  void dihedrals(int, char *);
+  void impropers(int, char *);
+  void nspecial_read(int, char *);
+  void special_read(char *);
+
+  void initialize();
+  void allocate();
+  void deallocate();
+
+  void open(char *);
+  void readline(char *);
+  void parse_keyword(int, char *, char *);
+  void skip_lines(int, char *);
+  int parse(char *, char **, int);
+
+  // void print();
+};
+
+}
+
+#endif