diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 53b10ddba5..482795edbe 100755 --- a/src/MC/fix_bond_break.cpp +++ b/src/MC/fix_bond_break.cpp @@ -128,13 +128,14 @@ int FixBondBreak::setmask() void FixBondBreak::init() { - // require special bonds = 0,1,1 + // require special bonds = *,1,1 + // [0] can be anything b/c I,J are removed from each other's special list + // [1],[2] must be 1.0 b/c only special lists of I,J are updated when + // bond I-J is broken, not special lists of neighbors of I,J,etc int flag = 0; - if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || - force->special_lj[3] != 1.0) flag = 1; - if (force->special_coul[1] != 0.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0) flag = 1; + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) flag = 1; + if (force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) flag = 1; if (flag) error->all(FLERR,"Fix bond/break requires special_bonds = 0,1,1"); // warn if angles, dihedrals, impropers are being used diff --git a/src/MC/fix_bond_create.cpp b/src/MC/fix_bond_create.cpp index ee7cfa2767..3e9ab83518 100755 --- a/src/MC/fix_bond_create.cpp +++ b/src/MC/fix_bond_create.cpp @@ -180,16 +180,17 @@ void FixBondCreate::init() if (force->pair == NULL || cutsq > force->pair->cutsq[iatomtype][jatomtype]) error->all(FLERR,"Fix bond/create cutoff is longer than pairwise cutoff"); - // require special bonds = 0,1,1 + // require special bonds = *,1,1 + // [0] can be anything b/c duplicate bond is checked for + // [1],[2] must be 1 b/c only special lists of I,J are updated when + // bond I-J is created, not special lists of neighbors of I,J,etc - if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || - force->special_lj[3] != 1.0) - error->all(FLERR,"Fix bond/create requires special_bonds lj = 0,1,1"); + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) + error->all(FLERR,"Fix bond/create requires special_bonds lj = *,1,1"); if (atom->q_flag) - if (force->special_coul[1] != 0.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0) - error->all(FLERR,"Fix bond/create requires special_bonds coul = 0,1,1"); + if (force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) + error->all(FLERR,"Fix bond/create requires special_bonds coul = *,1,1"); // warn if angles, dihedrals, impropers are being used @@ -267,7 +268,7 @@ void FixBondCreate::setup(int vflag) void FixBondCreate::post_integrate() { - int i,j,m,ii,jj,inum,jnum,itype,jtype,n1,n3,possible; + int i,j,k,m,ii,jj,inum,jnum,itype,jtype,n1,n3,possible; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *ilist,*jlist,*numneigh,**firstneigh; tagint *slist; @@ -309,6 +310,8 @@ void FixBondCreate::post_integrate() double **x = atom->x; tagint *tag = atom->tag; + tagint **bond_atom = atom->bond_atom; + int *num_bond = atom->num_bond; int *mask = atom->mask; int *type = atom->type; @@ -345,6 +348,16 @@ void FixBondCreate::post_integrate() } if (!possible) continue; + // do not allow a duplicate bond to be created + // check existing bonds of both I and J + + for (k = 0; k < num_bond[i]; k++) + if (bond_atom[i][k] == tag[j]) possible = 0; + if (j < nlocal) + for (k = 0; k < num_bond[j]; k++) + if (bond_atom[j][k] == tag[i]) possible = 0; + if (!possible) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -386,8 +399,6 @@ void FixBondCreate::post_integrate() // and probability constraint is satisfied int **bond_type = atom->bond_type; - tagint **bond_atom = atom->bond_atom; - int *num_bond = atom->num_bond; int **nspecial = atom->nspecial; tagint **special = atom->special; int newton_bond = force->newton_bond; @@ -437,11 +448,10 @@ void FixBondCreate::post_integrate() // increment bondcount, convert atom to new type if limit reached bondcount[i]++; - if (type[i] == iatomtype) { + if (type[i] == iatomtype) if (bondcount[i] == imaxbond) type[i] = inewtype; - } else { + else if (bondcount[i] == jmaxbond) type[i] = jnewtype; - } // count the created bond once diff --git a/src/fix.h b/src/fix.h index 6f11272be6..4be3991f16 100644 --- a/src/fix.h +++ b/src/fix.h @@ -196,6 +196,10 @@ class Fix : protected Pointers { void v_setup(int); void v_tally(int, int *, double, double *); + inline int sbmask(int j) { + return j >> SBBITS & 3; + } + // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation diff --git a/src/read_data.cpp b/src/read_data.cpp index 628f9f1e59..4c4a6ad995 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -139,10 +139,9 @@ void ReadData::command(int narg, char **arg) } else error->all(FLERR,"Illegal read_data command"); } - // perform 1-pass read if no molecular topoogy in file - // perform 2-pass read if molecular topology - // 1st pass calculates max topology/atom + // perform 2-pass read if molecular topology, + // first pass calculates max topology/atom int atomflag,topoflag; int bondflag,angleflag,dihedralflag,improperflag; @@ -180,7 +179,6 @@ void ReadData::command(int narg, char **arg) atom->allocate_type_arrays(); atom->avec->grow(n); - n = atom->nmax; domain->print_box(" "); domain->set_initial_box(); @@ -428,16 +426,19 @@ void ReadData::command(int narg, char **arg) if (!topoflag) break; firstpass = 0; - // reallocate bond,angle,diehdral,improper arrays via grow(), - // using new bond,angle,dihedral,improper per-atom values from 1st pass + // reallocate bond,angle,diehdral,improper arrays via grow() + // use new bond,angle,dihedral,improper per-atom values from 1st pass // should leave other atom arrays unchanged, since already nmax in length + // if bonds/etc not in data file, initialize per-atom size + // with extra settings before grow() of these topology arrays if (bondflag) { memory->destroy(atom->bond_type); memory->destroy(atom->bond_atom); atom->bond_type = NULL; atom->bond_atom = NULL; - } + } else atom->bond_per_atom = atom->extra_bond_per_atom; + if (angleflag) { memory->destroy(atom->angle_type); memory->destroy(atom->angle_atom1); @@ -445,7 +446,8 @@ void ReadData::command(int narg, char **arg) memory->destroy(atom->angle_atom3); atom->angle_type = NULL; atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; - } + } else atom->angle_per_atom = atom->extra_angle_per_atom; + if (dihedralflag) { memory->destroy(atom->dihedral_type); memory->destroy(atom->dihedral_atom1); @@ -455,7 +457,8 @@ void ReadData::command(int narg, char **arg) atom->dihedral_type = NULL; atom->dihedral_atom1 = atom->dihedral_atom2 = atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; - } + } else atom->dihedral_per_atom = atom->extra_dihedral_per_atom; + if (improperflag) { memory->destroy(atom->improper_type); memory->destroy(atom->improper_atom1); @@ -465,7 +468,7 @@ void ReadData::command(int narg, char **arg) atom->improper_type = NULL; atom->improper_atom1 = atom->improper_atom2 = atom->improper_atom3 = atom->improper_atom4 = NULL; - } + } else atom->improper_per_atom = atom->extra_improper_per_atom; atom->avec->grow(atom->nmax); }