diff --git a/src/molecule.cpp b/src/molecule.cpp index 8e772218c1..8ac2c037ee 100644 --- a/src/molecule.cpp +++ b/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); } /* ----------------------------------------------------------------------