ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2013-02-19 16:13:19 +00:00
parent bb9a46af64
commit 28cd1e1af6
3 changed files with 150 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
src/RIGID/fix_rigid_small.cpp
View File

@ -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<int> (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");
}
}
*/
}
*/

5
src/RIGID/fix_rigid_small.h
View File

@ -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);
};
}

67
src/delete_bonds.cpp
View File

@ -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];