diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 0703171e480d097b419a4230d096d58531875491..33a4b441fda4f9f3eefb3039db46c4ea2724a454 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -137,9 +137,10 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : int custom_flag = strcmp(arg[3],"custom") == 0; if (custom_flag) { if (narg < 5) error->all(FLERR,"Illegal fix rigid command"); + // determine whether atom-style variable or atom property is used. if (strstr(arg[4],"i_") == arg[4]) { - int is_double; + int is_double=0; int custom_index = atom->find_custom(arg[4]+2,is_double); if (custom_index == -1) error->all(FLERR,"Fix rigid custom requires previously defined property/atom"); @@ -153,13 +154,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world); molecule = new tagint[nlocal]; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) molecule[i] = (tagint)(value[i] - minval + 1); + if (mask[i] & groupbit) + molecule[i] = (tagint)(value[i] - minval + 1); + else + molecule[i] = 0; + } else if (strstr(arg[4],"v_") == arg[4]) { int ivariable = input->variable->find(arg[4]+2); if (ivariable < 0) error->all(FLERR,"Variable name for fix rigid custom does not exist"); if (input->variable->atomstyle(ivariable) == 0) - error->all(FLERR,"Fix rigid custom variable is not atom-style variable"); + error->all(FLERR,"Fix rigid custom variable is no atom-style variable"); double *value = new double[nlocal]; input->variable->compute_atom(ivariable,0,value,1,0); int minval = INT_MAX; @@ -169,7 +174,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world); molecule = new tagint[nlocal]; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) molecule[i] = (tagint)((int)value[i] - minval + 1); + if (mask[i] & groupbit) molecule[i] = (tagint)((tagint)value[i] - minval + 1); delete[] value; } else error->all(FLERR,"Unsupported fix rigid custom property"); } else { diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 1404c3bf58da22f20739dbf521c8aa2adcfea723..13208d4c4c434f4be2602a1447650dced9dca3f2 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -29,7 +29,9 @@ #include "group.h" #include "comm.h" #include "force.h" +#include "input.h" #include "output.h" +#include "variable.h" #include "random_mars.h" #include "math_const.h" #include "memory.h" @@ -64,11 +66,12 @@ enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), step_respa(NULL), - infile(NULL), body(NULL), bodyown(NULL), bodytag(NULL), atom2body(NULL), - xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL), - avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL), - itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL), id_dilate(NULL), - onemols(NULL), hash(NULL), bbox(NULL), ctr(NULL), idclose(NULL), rsqclose(NULL) + infile(NULL), body(NULL), bodyown(NULL), bodytag(NULL), atom2body(NULL), + xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL), + avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL), + itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL), + id_dilate(NULL), onemols(NULL), hash(NULL), bbox(NULL), ctr(NULL), + idclose(NULL), rsqclose(NULL) { int i; @@ -89,7 +92,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : // perform initial allocation of atom-based arrays // register with Atom class - extended = orientflag = dorientflag = 0; + extended = orientflag = dorientflag = customflag = 0; bodyown = NULL; bodytag = NULL; atom2body = NULL; @@ -103,24 +106,70 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : // parse args for rigid body specification + int *mask = atom->mask; + tagint *bodyid = NULL; + int nlocal = atom->nlocal; + if (narg < 4) error->all(FLERR,"Illegal fix rigid/small command"); - if (strcmp(arg[3],"molecule") != 0) - error->all(FLERR,"Illegal fix rigid/small command"); + if (strcmp(arg[3],"molecule") == 0) { + if (atom->molecule_flag == 0) + error->all(FLERR,"Fix rigid/small requires atom attribute molecule"); + + } else if (strcmp(arg[3],"custom") == 0) { + if (narg < 5) error->all(FLERR,"Illegal fix rigid/small command"); + bodyid = new tagint[nlocal]; + customflag = 1; + + // determine whether atom-style variable or atom property is used. + if (strstr(arg[4],"i_") == arg[4]) { + int is_double=0; + int custom_index = atom->find_custom(arg[4]+2,is_double); + if (custom_index == -1) + error->all(FLERR,"Fix rigid/small custom requires previously defined property/atom"); + else if (is_double) + error->all(FLERR,"Fix rigid/small custom requires integer-valued property/atom"); + + int minval = INT_MAX; + int *value = atom->ivector[custom_index]; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) minval = MIN(minval,value[i]); + int vmin = minval; + MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world); + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + bodyid[i] = (tagint)(value[i] - minval + 1); + else bodyid[i] = 0; + + } else if (strstr(arg[4],"v_") == arg[4]) { + int ivariable = input->variable->find(arg[4]+2); + if (ivariable < 0) + error->all(FLERR,"Variable name for fix rigid/small custom does not exist"); + if (input->variable->atomstyle(ivariable) == 0) + error->all(FLERR,"Fix rigid/small custom variable is no atom-style variable"); + double *value = new double[nlocal]; + input->variable->compute_atom(ivariable,0,value,1,0); + int minval = INT_MAX; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) minval = MIN(minval,(int)value[i]); + int vmin = minval; + MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world); + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + bodyid[i] = (tagint)((tagint)value[i] - minval + 1); + else bodyid[0] = 0; + delete[] value; + } else error->all(FLERR,"Unsupported fix rigid custom property"); + } else error->all(FLERR,"Illegal fix rigid/small command"); - if (atom->molecule_flag == 0) - error->all(FLERR,"Fix rigid/small requires atom attribute molecule"); if (atom->map_style == 0) error->all(FLERR,"Fix rigid/small requires an atom map, see atom_modify"); - // maxmol = largest molecule # - - int *mask = atom->mask; - tagint *molecule = atom->molecule; - int nlocal = atom->nlocal; - + // maxmol = largest bodyid # maxmol = -1; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]); + if (mask[i] & groupbit) maxmol = MAX(maxmol,bodyid[i]); tagint itmp; MPI_Allreduce(&maxmol,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world); @@ -155,6 +204,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : } int iarg = 4; + if (customflag) ++iarg; + while (iarg < narg) { if (strcmp(arg[iarg],"langevin") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal fix rigid/small command"); @@ -344,11 +395,12 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : if (pcouple == XYZ || (domain->dimension == 2 && pcouple == XY)) pstyle = ISO; else pstyle = ANISO; - // create rigid bodies based on molecule ID + // create rigid bodies based on molecule or custom ID // sets bodytag for owned atoms // body attributes are computed later by setup_bodies() - create_bodies(); + create_bodies(bodyid); + if (customflag) delete [] bodyid; // set nlocal_body and allocate bodies I own @@ -1424,7 +1476,7 @@ void FixRigidSmall::set_v() set bodytag for all owned atoms ------------------------------------------------------------------------- */ -void FixRigidSmall::create_bodies() +void FixRigidSmall::create_bodies(tagint *bodyid) { int i,m,n; double unwrap[3]; @@ -1464,8 +1516,8 @@ void FixRigidSmall::create_bodies() double *buf; memory->create(buf,ncount*percount,"rigid/small:buf"); - // create map hash for storing unique molecule IDs of my atoms - // key = molecule ID + // create map hash for storing unique body IDs of my atoms + // key = body ID // value = index into per-body data structure // n = # of entries in hash @@ -1477,12 +1529,10 @@ void FixRigidSmall::create_bodies() // value = index into N-length data structure // n = count of unique bodies my atoms are part of - tagint *molecule = atom->molecule; - n = 0; for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - if (hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = n++; + if (hash->find(bodyid[i]) == hash->end()) (*hash)[bodyid[i]] = n++; } // bbox = bounding box of each rigid body my atoms are part of @@ -1494,7 +1544,7 @@ void FixRigidSmall::create_bodies() bbox[i][1] = bbox[i][3] = bbox[i][5] = -BIG; } - // pack my atoms into buffer as molecule ID, unwrapped coords + // pack my atoms into buffer as body ID, unwrapped coords double **x = atom->x; @@ -1502,7 +1552,7 @@ void FixRigidSmall::create_bodies() for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; domain->unmap(x[i],image[i],unwrap); - buf[m++] = molecule[i]; + buf[m++] = bodyid[i]; buf[m++] = unwrap[0]; buf[m++] = unwrap[1]; buf[m++] = unwrap[2]; @@ -1542,7 +1592,7 @@ void FixRigidSmall::create_bodies() for (i = 0; i < n; i++) rsqclose[i] = BIG; - // pack my atoms into buffer as molecule ID, atom ID, unwrapped coords + // pack my atoms into buffer as body ID, atom ID, unwrapped coords tagint *tag = atom->tag; @@ -1550,7 +1600,7 @@ void FixRigidSmall::create_bodies() for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; domain->unmap(x[i],image[i],unwrap); - buf[m++] = molecule[i]; + buf[m++] = bodyid[i]; buf[m++] = ubuf(tag[i]).d; buf[m++] = unwrap[0]; buf[m++] = unwrap[1]; @@ -1570,7 +1620,7 @@ void FixRigidSmall::create_bodies() for (i = 0; i < nlocal; i++) { bodytag[i] = 0; if (!(mask[i] & groupbit)) continue; - m = hash->find(molecule[i])->second; + m = hash->find(bodyid[i])->second; bodytag[i] = idclose[m]; rsqmax = MAX(rsqmax,rsqclose[m]); } diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index b07dea4f333ea3f27b7919fb74a4b2c52fc9585e..22f9b0c16c31e7ba8910a558dc1802fff6d4b128 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -79,6 +79,7 @@ class FixRigidSmall : public Fix { char *infile; // file to read rigid body attributes from int setupflag; // 1 if body properties are setup, else 0 int commflag; // various modes of forward/reverse comm + int customflag; // 1 if custom property/variable define bodies int nbody; // total # of rigid bodies int nlinear; // total # of linear rigid bodies tagint maxmol; // max mol-ID @@ -187,7 +188,7 @@ class FixRigidSmall : public Fix { void image_shift(); void set_xv(); void set_v(); - void create_bodies(); + void create_bodies(tagint *); void setup_bodies_static(); void setup_bodies_dynamic(); void readfile(int, double **, int *);