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

src/molecule.cpp

@@ -16,6 +16,7 @@
 #include "molecule.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "atom_vec_body.h"
 #include "force.h"
 #include "comm.h"
 #include "domain.h"
@@ -163,7 +164,8 @@ Molecule::~Molecule()
    compute center = geometric center of molecule
    also compute:
      dx = displacement of each atom from center
-     molradius = radius of molecule from center including finite-size particles
+     molradius = radius of molecule from center 
+       including finite-size particles or body particles
 ------------------------------------------------------------------------- */
 
 void Molecule::compute_center()
@@ -446,6 +448,13 @@ void Molecule::read(int flag)
       itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale;
       itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale;
     }
+    else if (strstr(line,"body")) {
+      bodyflag = 1;
+      avec_body = (AtomVecBody *) atom->style_match("body");
+      if (!avec_body) 
+        error->all(FLERR,"Molecule file requires atom style body");
+      sscanf(line,"%d %d",&nibody,&ndbody);
+    }
 
     else break;
   }
@@ -542,6 +551,19 @@ void Molecule::read(int flag)
       if (flag) shaketype_read(line);
       else skip_lines(natoms,line);
 
+    } else if (strcmp(keyword,"Body Integers") == 0) {
+      if (bodyflag == 0 || nibody == 0)
+	error->all(FLERR,"Molecule file has body params "
+                   "but no setting for them");
+      ibodyflag = 1;
+      body(flag,0,line);
+    } else if (strcmp(keyword,"Body Doubles") == 0) {
+      if (bodyflag == 0 || ndbody == 0)
+	error->all(FLERR,"Molecule file has body params "
+                   "but no setting for them");
+      dbodyflag = 1;
+      body(flag,1,line);
+
     } else error->one(FLERR,"Unknown section in molecule file");
 
     parse_keyword(1,line,keyword);
@@ -560,6 +582,10 @@ void Molecule::read(int flag)
       error->all(FLERR,"Molecule file has special flags but no bonds");
     if ((shakeflagflag || shakeatomflag || shaketypeflag) && !shakeflag)
       error->all(FLERR,"Molecule file shake info is incomplete");
+    if (bodyflag && nibody && ibodyflag == 0)
+      error->all(FLERR,"Molecule file has no Body Integers section");
+    if (bodyflag && ndbody && dbodyflag == 0)
+      error->all(FLERR,"Molecule file has no Body Doubles section");
   }
 
   // auto-generate special bonds
@@ -570,6 +596,21 @@ void Molecule::read(int flag)
     maxspecial = atom->maxspecial;
     if (flag) special_generate();
   }
+
+  // body particle must have natom = 1
+  // set radius by having body class compute its own radius
+  
+  if (bodyflag) {
+    radiusflag = 1;
+    if (natoms != 1) 
+      error->all(FLERR,"Molecule natoms must be 1 for body particle");
+    if (sizescale != 1.0)
+      error->all(FLERR,"Molecule sizescale must be 1.0 for body particle");
+    if (flag) {
+      radius[0] = avec_body->radius_body(nibody,ndbody,ibodyparams,dbodyparams);
+      maxradius = radius[0];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1244,6 +1285,44 @@ void Molecule::shaketype_read(char *line)
   }
 }
 
+/* ----------------------------------------------------------------------
+   read body params from file
+   pflag = 0/1 for integer/double params
+------------------------------------------------------------------------- */
+
+void Molecule::body(int flag, int pflag, char *line)
+{
+  int i,ncount;
+
+  int nparam = nibody;
+  if (pflag) nparam = ndbody;
+
+  int nword = 0;
+  while (nword < nparam) {
+    readline(line);
+
+    ncount = atom->count_words(line);
+    if (ncount == 0)
+      error->one(FLERR,"Too few values in body section of molecule file");
+    if (nword+ncount > nparam) 
+      error->all(FLERR,"Too many values in body section of molecule file");
+    
+    if (flag) {
+      if (pflag == 0) {
+        ibodyparams[nword++] = force->inumeric(FLERR,strtok(line," \t\n\r\f"));
+        for (i = 1; i < ncount; i++)
+          ibodyparams[nword++] = 
+            force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
+      } else {
+        dbodyparams[nword++] = force->numeric(FLERR,strtok(line," \t\n\r\f"));
+        for (i = 1; i < ncount; i++)
+          dbodyparams[nword++] = 
+            force->numeric(FLERR,strtok(NULL," \t\n\r\f"));
+      }
+    } else nword += ncount;
+  }
+}
+
 /* ----------------------------------------------------------------------
    error check molecule attributes and topology against system settings
    flag = 0, just check this molecule
@@ -1318,6 +1397,7 @@ void Molecule::initialize()
   nbonds = nangles = ndihedrals = nimpropers = 0;
   ntypes = 0;
   nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
+  nibody = ndbody = 0;
 
   bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
   maxspecial = 0;
@@ -1326,6 +1406,7 @@ void Molecule::initialize()
   bondflag = angleflag = dihedralflag = improperflag = 0;
   nspecialflag = specialflag = 0;
   shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0;
+  bodyflag = ibodyflag = dbodyflag = 0;
 
   centerflag = massflag = comflag = inertiaflag = 0;
   tag_require = 0;
@@ -1359,6 +1440,9 @@ void Molecule::initialize()
   shake_atom = NULL;
   shake_type = NULL;
 
+  ibodyparams = NULL;
+  dbodyparams = NULL;
+
   dx = NULL;
   dxcom = NULL;
   dxbody = NULL;
@@ -1445,6 +1529,11 @@ void Molecule::allocate()
     memory->create(shake_atom,natoms,4,"molecule:shake_flag");
     memory->create(shake_type,natoms,3,"molecule:shake_flag");
   }
+
+  if (bodyflag) {
+    if (nibody) memory->create(ibodyparams,nibody,"molecule:ibodyparams");
+    if (ndbody) memory->create(dbodyparams,ndbody,"molecule:dbodyparams");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1493,6 +1582,9 @@ void Molecule::deallocate()
   memory->destroy(dx);
   memory->destroy(dxcom);
   memory->destroy(dxbody);
+
+  memory->destroy(ibodyparams);
+  memory->destroy(dbodyparams);
 }
 
 /* ----------------------------------------------------------------------