diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 7c2aad2e85..46c5459e57 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -193,12 +193,18 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : if (nbody == 0) error->all(FLERR,"No rigid bodies defined"); if (me == 0) { - if (screen) fprintf(screen, - "%d rigid bodies with " BIGINT_FORMAT " atoms\n", - nbody,atomall); - if (logfile) fprintf(logfile, - "%d rigid bodies with " BIGINT_FORMAT " atoms\n", - nbody,atomall); + if (screen) { + fprintf(screen,"%d rigid bodies with " BIGINT_FORMAT " atoms\n", + nbody,atomall); + fprintf(screen," %g = max distance from body owner to body atom\n", + maxextent); + } + if (logfile) { + fprintf(logfile,"%d rigid bodies with " BIGINT_FORMAT " atoms\n", + nbody,atomall); + fprintf(logfile," %g = max distance from body owner to body atom\n", + maxextent); + } } // initialize Marsaglia RNG with processor-unique seed @@ -1318,7 +1324,7 @@ void FixRigidSmall::create_bodies() } // idclose = ID of atom in body closest to center pt (smaller ID if tied) - // rsq = distance squared from idclose to center pt + // rsqclose = distance squared from idclose to center pt memory->create(idclose,n,"rigid/small:idclose"); memory->create(rsqclose,n,"rigid/small:rsqclose"); @@ -1341,21 +1347,49 @@ void FixRigidSmall::create_bodies() } // pass buffer around ring of procs - // func = idclose,rsqclose with atom IDs from every proc + // func = update idclose,rsqclose with atom IDs from every proc // when done, have idclose for every rigid body my atoms are part of frsptr = this; comm->ring(m,sizeof(double),buf,2,ring_nearest,NULL); // set bodytag of all owned atoms, based on idclose + // find max value of rsqclose across all procs + double rsqmax = 0.0; for (i = 0; i < nlocal; i++) { bodytag[i] = 0; if (!(mask[i] & groupbit)) continue; m = hash->find(molecule[i])->second; bodytag[i] = idclose[m]; + rsqmax = MAX(rsqmax,rsqclose[m]); } + // pack my atoms into buffer as bodytag of owning atom, unwrapped coords + + m = 0; + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + domain->unmap(x[i],image[i],unwrap); + buf[m++] = bodytag[i]; + buf[m++] = unwrap[0]; + buf[m++] = unwrap[1]; + buf[m++] = unwrap[2]; + } + + // pass buffer around ring of procs + // func = update rsqfar for atoms belonging to bodies I own + // when done, have rsqfar for all atoms in bodies I own + + rsqfar = 0.0; + frsptr = this; + comm->ring(m,sizeof(double),buf,3,ring_farthest,NULL); + + // find maxextent of rsqfar across all procs + + MPI_Allreduce(&rsqfar,&maxextent,1,MPI_DOUBLE,MPI_MAX,world); + maxextent = sqrt(maxextent); + // clean up delete hash; @@ -1437,6 +1471,40 @@ void FixRigidSmall::ring_nearest(int n, char *cbuf) } } +/* ---------------------------------------------------------------------- + process rigid body atoms from another proc + update rsqfar = distance from owning atom to other atom +------------------------------------------------------------------------- */ + +void FixRigidSmall::ring_farthest(int n, char *cbuf) +{ + double **x = frsptr->atom->x; + int *image = frsptr->atom->image; + int nlocal = frsptr->atom->nlocal; + + double *buf = (double *) cbuf; + int ndatums = n/4; + + int itag,iowner; + double delx,dely,delz,rsq; + double *xx; + double unwrap[3]; + + int m = 0; + for (int i = 0; i < ndatums; i++, m += 4) { + itag = static_cast (buf[m]); + iowner = frsptr->atom->map(itag); + if (iowner < 0 || iowner >= nlocal) continue; + frsptr->domain->unmap(x[iowner],image[iowner],unwrap); + xx = &buf[m+1]; + delx = xx[0] - unwrap[0]; + dely = xx[1] - unwrap[1]; + delz = xx[2] - unwrap[2]; + rsq = delx*delx + dely*dely + delz*delz; + frsptr->rsqfar = MAX(frsptr->rsqfar,rsq); + } +} + /* ---------------------------------------------------------------------- one-time initialization of rigid body attributes extended flags, mass, center-of-mass @@ -1908,6 +1976,10 @@ void FixRigidSmall::setup_bodies() fabs(itensor[ibody][5]/norm) > TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments"); } + + // clean up + + memory->destroy(itensor); } /* ---------------------------------------------------------------------- @@ -2457,26 +2529,26 @@ void FixRigidSmall::reset_dt() debug method for sanity checking of atom/body data pointers ------------------------------------------------------------------------- */ +/* void FixRigidSmall::check(int flag) { - /* for (int i = 0; i < atom->nlocal; i++) { if (bodyown[i] >= 0) { if (bodytag[i] != atom->tag[i]) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD AAA"); + errorx->one(FLERR,"BAD AAA"); } if (bodyown[i] < 0 || bodyown[i] >= nlocal_body) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD BBB"); + errorx->one(FLERR,"BAD BBB"); } if (atom2body[i] != bodyown[i]) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD CCC"); + errorx->one(FLERR,"BAD CCC"); } if (body[bodyown[i]].ilocal != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD DDD"); + errorx->one(FLERR,"BAD DDD"); } } } @@ -2485,11 +2557,11 @@ void FixRigidSmall::check(int flag) if (bodyown[i] < 0 && bodytag[i] > 0) { if (atom2body[i] < 0 || atom2body[i] >= nlocal_body+nghost_body) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD EEE"); + errorx->one(FLERR,"BAD EEE"); } if (bodytag[i] != atom->tag[body[atom2body[i]].ilocal]) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD FFF"); + errorx->one(FLERR,"BAD FFF"); } } } @@ -2501,11 +2573,11 @@ void FixRigidSmall::check(int flag) printf("Values %d %d: %d %d %d\n", i,atom->tag[i],bodyown[i],nlocal_body,nghost_body); printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD GGG"); + errorx->one(FLERR,"BAD GGG"); } if (body[bodyown[i]].ilocal != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD HHH"); + errorx->one(FLERR,"BAD HHH"); } } } @@ -2513,12 +2585,12 @@ void FixRigidSmall::check(int flag) for (int i = 0; i < nlocal_body; i++) { if (body[i].ilocal < 0 || body[i].ilocal >= atom->nlocal) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD III"); + errorx->one(FLERR,"BAD III"); } if (bodytag[body[i].ilocal] != atom->tag[body[i].ilocal] || bodyown[body[i].ilocal] != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD JJJ"); + errorx->one(FLERR,"BAD JJJ"); } } @@ -2526,12 +2598,12 @@ void FixRigidSmall::check(int flag) if (body[i].ilocal < atom->nlocal || body[i].ilocal >= atom->nlocal + atom->nghost) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD KKK"); + errorx->one(FLERR,"BAD KKK"); } if (bodyown[body[i].ilocal] != i) { printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag); - error->one(FLERR,"BAD LLL"); + errorx->one(FLERR,"BAD LLL"); } } - */ } +*/ diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index abf3bcd0f9..2e5e9f0bc4 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -71,6 +71,7 @@ class FixRigidSmall : public Fix { int firstflag; // 1 for first-time setup of rigid bodies int commflag; // various modes of forward/reverse comm int nbody; // total # of rigid bodies + double maxextent; // furthest distance from body owner to body atom struct Body { double mass; // total mass of body @@ -141,6 +142,7 @@ class FixRigidSmall : public Fix { double **ctr; int *idclose; double *rsqclose; + double rsqfar; void set_xv(); void set_v(); @@ -153,10 +155,11 @@ class FixRigidSmall : public Fix { static void ring_bbox(int, char *); static void ring_nearest(int, char *); + static void ring_farthest(int, char *); // debug - void check(int); + //void check(int); }; } diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 6699ab7478..078e935ef7 100644 --- a/src/delete_bonds.cpp +++ b/src/delete_bonds.cpp @@ -86,12 +86,14 @@ void DeleteBonds::command(int narg, char **arg) int undo_flag = 0; int remove_flag = 0; int special_flag = 0; + int induce_flag = 0; while (iarg < narg) { if (strcmp(arg[iarg],"any") == 0) any_flag = 1; else if (strcmp(arg[iarg],"undo") == 0) undo_flag = 1; else if (strcmp(arg[iarg],"remove") == 0) remove_flag = 1; else if (strcmp(arg[iarg],"special") == 0) special_flag = 1; + else if (strcmp(arg[iarg],"induce") == 0) induce_flag = 1; else error->all(FLERR,"Illegal delete_bonds command"); iarg++; } @@ -109,9 +111,9 @@ void DeleteBonds::command(int narg, char **arg) // set topology interactions either off or on // criteria for an interaction to potentially be changed (set flag = 1) - // all atoms or one atom in interaction must be in group, based on any_flag - // for style = MULTI, no other criteria - // for style = ATOM, at least one atom is specified type + // all atoms or any atom in interaction must be in group, based on any_flag + // for style = MULTI, all bond/angle/dihedral/improper, no other criteria + // for style = ATOM, same as MULTI, at least one atom is specified type // for style = BOND/ANGLE/DIHEDRAL/IMPROPER, interaction is specified type // for style = STATS only compute stats, flag is always 0 // if flag = 1 @@ -125,7 +127,8 @@ void DeleteBonds::command(int narg, char **arg) int i,m,n,consider,flag; int atom1,atom2,atom3,atom4; - if (atom->avec->bonds_allow) { + if (atom->avec->bonds_allow && + (style == BOND || style == MULTI || style == ATOM)) { int *num_bond = atom->num_bond; int **bond_type = atom->bond_type; @@ -155,7 +158,8 @@ void DeleteBonds::command(int narg, char **arg) } } - if (atom->avec->angles_allow) { + if (atom->avec->angles_allow && + (style == ANGLE || style == MULTI || style == ATOM)) { int *num_angle = atom->num_angle; int **angle_type = atom->angle_type; @@ -189,7 +193,8 @@ void DeleteBonds::command(int narg, char **arg) } } - if (atom->avec->dihedrals_allow) { + if (atom->avec->dihedrals_allow && + (style == DIHEDRAL || style == MULTI || style == ATOM)) { int *num_dihedral = atom->num_dihedral; int **dihedral_type = atom->dihedral_type; @@ -225,7 +230,8 @@ void DeleteBonds::command(int narg, char **arg) } } - if (atom->avec->impropers_allow) { + if (atom->avec->impropers_allow && + (style == IMPROPER || style == MULTI || style == ATOM)) { int *num_improper = atom->num_improper; int **improper_type = atom->improper_type; @@ -261,8 +267,20 @@ void DeleteBonds::command(int narg, char **arg) } } + // induce turn off of angles, dihedral, impropers due to turned off bonds + // induce turn off of dihedrals due to turned off angles + // all atoms or any atom in interaction must be in group, based on any_flag + + if (induce_flag) { + + // circulate list of turned off bonds around ring of procs + + // circulate list of turned off angles around ring of procs + + } + // remove interactions if requested - // only if all atoms in bond, angle, etc are in the delete_bonds group + // all atoms or any atom in interaction must be in group, based on any_flag if (remove_flag) { @@ -272,7 +290,12 @@ void DeleteBonds::command(int narg, char **arg) while (m < atom->num_bond[i]) { if (atom->bond_type[i][m] <= 0) { atom1 = atom->map(atom->bond_atom[i][m]); - if (mask[i] & groupbit && mask[atom1] & groupbit) { + flag = 0; + if (!any_flag && mask[i] & groupbit && mask[atom1] & groupbit) + flag = 1; + if (any_flag && (mask[i] & groupbit || mask[atom1] & groupbit)) + flag = 1; + if (flag) { n = atom->num_bond[i]; atom->bond_type[i][m] = atom->bond_type[i][n-1]; atom->bond_atom[i][m] = atom->bond_atom[i][n-1]; @@ -291,8 +314,12 @@ void DeleteBonds::command(int narg, char **arg) atom1 = atom->map(atom->angle_atom1[i][m]); atom2 = atom->map(atom->angle_atom2[i][m]); atom3 = atom->map(atom->angle_atom3[i][m]); - if (mask[atom1] & groupbit && mask[atom2] & groupbit && - mask[atom3] & groupbit) { + flag = 0; + if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit && + mask[atom3] & groupbit) flag = 1; + if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit || + mask[atom3] & groupbit)) flag = 1; + if (flag) { n = atom->num_angle[i]; atom->angle_type[i][m] = atom->angle_type[i][n-1]; atom->angle_atom1[i][m] = atom->angle_atom1[i][n-1]; @@ -314,8 +341,13 @@ void DeleteBonds::command(int narg, char **arg) atom2 = atom->map(atom->dihedral_atom2[i][m]); atom3 = atom->map(atom->dihedral_atom3[i][m]); atom4 = atom->map(atom->dihedral_atom4[i][m]); - if (mask[atom1] & groupbit && mask[atom2] & groupbit && - mask[atom3] & groupbit && mask[atom4] & groupbit) { + flag = 0; + if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit && + mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1; + if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit || + mask[atom3] & groupbit || mask[atom4] & groupbit)) + flag = 1; + if (flag) { n = atom->num_dihedral[i]; atom->dihedral_type[i][m] = atom->dihedral_type[i][n-1]; atom->dihedral_atom1[i][m] = atom->dihedral_atom1[i][n-1]; @@ -338,8 +370,13 @@ void DeleteBonds::command(int narg, char **arg) atom2 = atom->map(atom->improper_atom2[i][m]); atom3 = atom->map(atom->improper_atom3[i][m]); atom4 = atom->map(atom->improper_atom4[i][m]); - if (mask[atom1] & groupbit && mask[atom2] & groupbit && - mask[atom3] & groupbit && mask[atom4] & groupbit) { + flag = 0; + if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit && + mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1; + if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit || + mask[atom3] & groupbit || mask[atom4] & groupbit)) + flag = 1; + if (flag) { n = atom->num_improper[i]; atom->improper_type[i][m] = atom->improper_type[i][n-1]; atom->improper_atom1[i][m] = atom->improper_atom1[i][n-1];