diff --git a/src/atom.cpp b/src/atom.cpp index 24362be4d4..9f822a6071 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1318,45 +1318,51 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset) void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body, tagint id_offset) { - int j,m,tagdata,ninteger,ndouble; + int j,m,nvalues,tagdata,ninteger,ndouble; int maxint = 0; int maxdouble = 0; - char **ivalues = NULL; - char **dvalues = NULL; + int *ivalues = NULL; + double *dvalues = NULL; // loop over lines of body data - // tokenize the lines into ivalues and dvalues - // if I own atom tag, unpack its values + // if I own atom tag, tokenize lines into ivalues/dvalues, call data_body() + // else skip values for (int i = 0; i < n; i++) { if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset; else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset; - ninteger = atoi(strtok(NULL," \t\n\r\f")); - ndouble = atoi(strtok(NULL," \t\n\r\f")); - - if (ninteger > maxint) { - delete [] ivalues; - maxint = ninteger; - ivalues = new char*[maxint]; - } - if (ndouble > maxdouble) { - delete [] dvalues; - maxdouble = ndouble; - dvalues = new char*[maxdouble]; - } - - for (j = 0; j < ninteger; j++) - ivalues[j] = strtok(NULL," \t\n\r\f"); - for (j = 0; j < ndouble; j++) - dvalues[j] = strtok(NULL," \t\n\r\f"); - if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bodies section of data file"); + + if ((m = map(tagdata)) >= 0) { + ninteger = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); + ndouble = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); - if ((m = map(tagdata)) >= 0) + if (ninteger > maxint) { + delete [] ivalues; + maxint = ninteger; + ivalues = new int[maxint]; + } + if (ndouble > maxdouble) { + delete [] dvalues; + maxdouble = ndouble; + dvalues = new double[maxdouble]; + } + + for (j = 0; j < ninteger; j++) + ivalues[j] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); + for (j = 0; j < ndouble; j++) + dvalues[j] = force->numeric(FLERR,strtok(NULL," \t\n\r\f")); + avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues); + + } else { + nvalues = 2 + ninteger + ndouble; // number of values to skip + for (j = 0; j < nvalues; j++) + strtok(NULL," \t\n\r\f"); + } } delete [] ivalues; @@ -1588,7 +1594,13 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom, else if (rmass_flag) rmass[ilocal] = 4.0*MY_PI/3.0 * radius[ilocal]*radius[ilocal]*radius[ilocal]; - + if (onemol->bodyflag) { + body[ilocal] = 0; // as if a body read from data file + onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody, + onemol->ibodyparams,onemol->dbodyparams); + onemol->avec_body->set_quat(ilocal,onemol->quat_external); + } + if (molecular != 1) return; // add bond topology info diff --git a/src/set.cpp b/src/set.cpp index 33f6b77dc7..8ee5d7654b 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -21,6 +21,7 @@ #include "atom_vec_ellipsoid.h" #include "atom_vec_line.h" #include "atom_vec_tri.h" +#include "atom_vec_body.h" #include "domain.h" #include "region.h" #include "group.h" @@ -41,7 +42,7 @@ using namespace MathConst; enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI, - DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,ANGMOM,OMEGA, + DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA, DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,INAME,DNAME}; @@ -221,7 +222,7 @@ void Set::command(int narg, char **arg) else zvalue = force->numeric(FLERR,arg[iarg+3]); if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) varparse(arg[iarg+4],4); else wvalue = force->numeric(FLERR,arg[iarg+4]); - if (!atom->ellipsoid_flag && !atom->tri_flag) + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); set(QUAT); iarg += 5; @@ -229,7 +230,7 @@ void Set::command(int narg, char **arg) } else if (strcmp(arg[iarg],"quat/random") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); ivalue = force->inumeric(FLERR,arg[iarg+1]); - if (!atom->ellipsoid_flag && !atom->tri_flag) + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); if (ivalue <= 0) error->all(FLERR,"Invalid random number seed in set command"); @@ -248,6 +249,16 @@ void Set::command(int narg, char **arg) set(THETA); iarg += 2; + } else if (strcmp(arg[iarg],"theta/random") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); + ivalue = force->inumeric(FLERR,arg[iarg+1]); + if (!atom->line_flag) + error->all(FLERR,"Cannot set this attribute for this atom style"); + if (ivalue <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + set(THETA_RANDOM); + iarg += 2; + } else if (strcmp(arg[iarg],"angmom") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal set command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); @@ -256,7 +267,7 @@ void Set::command(int narg, char **arg) else yvalue = force->numeric(FLERR,arg[iarg+2]); if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3); else zvalue = force->numeric(FLERR,arg[iarg+3]); - if (!atom->ellipsoid_flag && !atom->tri_flag) + if (!atom->angmom_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); set(ANGMOM); iarg += 4; @@ -269,7 +280,7 @@ void Set::command(int narg, char **arg) else yvalue = force->numeric(FLERR,arg[iarg+2]); if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3); else zvalue = force->numeric(FLERR,arg[iarg+3]); - if (!atom->sphere_flag) + if (!atom->omega_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); set(OMEGA); iarg += 4; @@ -558,6 +569,7 @@ void Set::set(int keyword) (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"); int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { @@ -604,8 +616,9 @@ void Set::set(int keyword) else if (keyword == MESO_E) atom->e[i] = dvalue; else if (keyword == MESO_CV) atom->cv[i] = dvalue; else if (keyword == MESO_RHO) atom->rho[i] = dvalue; - else if (keyword == SMD_MASS_DENSITY) { // set mass from volume and supplied mass density - atom->rmass[i] = atom->vfrac[i] * dvalue; + else if (keyword == SMD_MASS_DENSITY) { + // set mass from volume and supplied mass density + atom->rmass[i] = atom->vfrac[i] * dvalue; } else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue; @@ -678,17 +691,23 @@ void Set::set(int keyword) mu[i][2]*mu[i][2]); } - // set quaternion orientation of ellipsoid or tri particle - + // set quaternion orientation of ellipsoid or tri or body particle + // enforce quat rotation vector in z dir for 2d systems + else if (keyword == QUAT) { double *quat; if (avec_ellipsoid && atom->ellipsoid[i] >= 0) quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; else if (avec_tri && atom->tri[i] >= 0) quat = avec_tri->bonus[atom->tri[i]].quat; + else if (avec_body && atom->body[i] >= 0) + quat = avec_body->bonus[atom->body[i]].quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); - + if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0)) + error->one(FLERR,"Cannot set quaternion with xy components " + "for 2d system"); + double theta2 = MY_PI2 * wvalue/180.0; double sintheta2 = sin(theta2); quat[0] = cos(theta2); @@ -706,7 +725,7 @@ void Set::set(int keyword) avec_line->bonus[atom->line[i]].theta = dvalue; } - // set angmom of ellipsoidal or tri particle + // set angmom or omega of particle else if (keyword == ANGMOM) { atom->angmom[i][0] = xvalue; @@ -720,7 +739,6 @@ void Set::set(int keyword) atom->omega[i][2] = zvalue; } - // reset any or all of 3 image flags else if (keyword == IMAGE) { @@ -768,7 +786,9 @@ void Set::setrandom(int keyword) 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"); RanPark *random = new RanPark(lmp,1); double **x = atom->x; @@ -828,7 +848,7 @@ void Set::setrandom(int keyword) } } - // set quaternions to random orientations in 3d or 2d + // set quaternions to random orientations in 3d and 2d } else if (keyword == QUAT_RANDOM) { int nlocal = atom->nlocal; @@ -842,6 +862,8 @@ void Set::setrandom(int keyword) quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; else if (avec_tri && atom->tri[i] >= 0) quat = avec_tri->bonus[atom->tri[i]].quat; + else if (avec_body && atom->body[i] >= 0) + quat = avec_body->bonus[atom->body[i]].quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); @@ -864,6 +886,8 @@ void Set::setrandom(int keyword) if (select[i]) { if (avec_ellipsoid && atom->ellipsoid[i] >= 0) quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; + else if (avec_body && atom->body[i] >= 0) + quat = avec_body->bonus[atom->body[i]].quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); @@ -876,6 +900,21 @@ void Set::setrandom(int keyword) count++; } } + + // set theta to random orientation in 2d + + } else if (keyword == QUAT_RANDOM) { + int nlocal = atom->nlocal; + double theta; + for (i = 0; i < nlocal; i++) { + if (select[i]) { + if (atom->line[i] < 0) + error->one(FLERR,"Cannot set theta for atom that is not a line"); + random->reset(seed,x[i]); + avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform(); + count++; + } + } } delete random;