diff --git a/src/neigh_bond.cpp b/src/neigh_bond.cpp index 4840c19639..2cbda027c3 100644 --- a/src/neigh_bond.cpp +++ b/src/neigh_bond.cpp @@ -66,6 +66,8 @@ void Neighbor::bond_all() nbondlist++; } } + + if (cluster_check) bond_check(); } /* ---------------------------------------------------------------------- */ @@ -106,6 +108,33 @@ void Neighbor::bond_partial() nbondlist++; } } + + if (cluster_check) bond_check(); +} + +/* ---------------------------------------------------------------------- */ + +void Neighbor::bond_check() +{ + int i,j; + double dx,dy,dz,dxstart,dystart,dzstart; + + double **x = atom->x; + int flag = 0; + + for (int m = 0; m < nbondlist; m++) { + i = bondlist[m][0]; + j = bondlist[m][1]; + dxstart = dx = x[i][0] - x[j][0]; + dystart = dy = x[i][1] - x[j][1]; + dzstart = dz = x[i][2] - x[j][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + } + + int flag_all; + MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); + if (flag_all) error->all(FLERR,"Bond extent > half of periodic box length"); } /* ---------------------------------------------------------------------- */ @@ -153,6 +182,8 @@ void Neighbor::angle_all() nanglelist++; } } + + if (cluster_check) angle_check(); } /* ---------------------------------------------------------------------- */ @@ -201,6 +232,39 @@ void Neighbor::angle_partial() nanglelist++; } } + + if (cluster_check) angle_check(); +} + +/* ---------------------------------------------------------------------- */ + +void Neighbor::angle_check() +{ + int i,j,k; + double dx,dy,dz,dxstart,dystart,dzstart; + + double **x = atom->x; + int flag = 0; + + for (int m = 0; m < nbondlist; m++) { + i = anglelist[m][0]; + j = anglelist[m][1]; + k = anglelist[m][1]; + dxstart = dx = x[i][0] - x[j][0]; + dystart = dy = x[i][1] - x[j][1]; + dzstart = dz = x[i][2] - x[j][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + dxstart = dx = x[i][0] - x[k][0]; + dystart = dy = x[i][1] - x[k][1]; + dzstart = dz = x[i][2] - x[k][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + } + + int flag_all; + MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); + if (flag_all) error->all(FLERR,"Angle extent > half of periodic box length"); } /* ---------------------------------------------------------------------- */ @@ -254,6 +318,8 @@ void Neighbor::dihedral_all() ndihedrallist++; } } + + if (cluster_check) dihedral_check(dihedrallist); } /* ---------------------------------------------------------------------- */ @@ -308,6 +374,46 @@ void Neighbor::dihedral_partial() ndihedrallist++; } } + + if (cluster_check) dihedral_check(dihedrallist); +} + +/* ---------------------------------------------------------------------- */ + +void Neighbor::dihedral_check(int **list) +{ + int i,j,k,l; + double dx,dy,dz,dxstart,dystart,dzstart; + + double **x = atom->x; + int flag = 0; + + for (int m = 0; m < nbondlist; m++) { + i = list[m][0]; + j = list[m][1]; + k = list[m][1]; + l = list[m][1]; + dxstart = dx = x[i][0] - x[j][0]; + dystart = dy = x[i][1] - x[j][1]; + dzstart = dz = x[i][2] - x[j][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + dxstart = dx = x[i][0] - x[k][0]; + dystart = dy = x[i][1] - x[k][1]; + dzstart = dz = x[i][2] - x[k][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + dxstart = dx = x[i][0] - x[l][0]; + dystart = dy = x[i][1] - x[l][1]; + dzstart = dz = x[i][2] - x[l][2]; + domain->minimum_image(dx,dy,dz); + if (dx != dxstart || dy != dystart || dz != dzstart) flag = 1; + } + + int flag_all; + MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); + if (flag_all) + error->all(FLERR,"Dihedral/improper extent > half of periodic box length"); } /* ---------------------------------------------------------------------- */ @@ -346,7 +452,7 @@ void Neighbor::improper_all() atom1 = domain->closest_image(i,atom1); atom2 = domain->closest_image(i,atom2); atom3 = domain->closest_image(i,atom3); - atom4 = domain->closest_image(i,atom4); + atom4 = domain-> closest_image(i,atom4); if (newton_bond || (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)) { if (nimproperlist == maximproper) { @@ -361,6 +467,8 @@ void Neighbor::improper_all() nimproperlist++; } } + + if (cluster_check) dihedral_check(improperlist); } /* ---------------------------------------------------------------------- */ @@ -415,4 +523,6 @@ void Neighbor::improper_partial() nimproperlist++; } } + + if (cluster_check) dihedral_check(improperlist); } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 26dc09ccc0..676688fef9 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -68,6 +68,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) oneatom = 2000; binsizeflag = 0; build_once = 0; + cluster_check = 0; cutneighsq = NULL; cutneighghostsq = NULL; @@ -1672,6 +1673,13 @@ void Neighbor::modify_params(int narg, char **arg) if (binsize_user <= 0.0) binsizeflag = 0; else binsizeflag = 1; iarg += 2; + } else if (strcmp(arg[iarg],"cluster") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command"); + if (strcmp(arg[iarg+1],"yes") == 0) cluster_check = 1; + else if (strcmp(arg[iarg+1],"no") == 0) cluster_check = 0; + else error->all(FLERR,"Illegal neigh_modify command"); + iarg += 2; + } else if (strcmp(arg[iarg],"include") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command"); includegroup = group->find(arg[iarg+1]); @@ -1682,6 +1690,7 @@ void Neighbor::modify_params(int narg, char **arg) error->all(FLERR, "Neigh_modify include group != atom_modify first group"); iarg += 2; + } else if (strcmp(arg[iarg],"exclude") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command"); @@ -1729,6 +1738,7 @@ void Neighbor::modify_params(int narg, char **arg) } else if (strcmp(arg[iarg+1],"none") == 0) { nex_type = nex_group = nex_mol = 0; iarg += 2; + } else error->all(FLERR,"Illegal neigh_modify command"); } else error->all(FLERR,"Illegal neigh_modify command"); diff --git a/src/neighbor.h b/src/neighbor.h index 6bcb38ae55..3c8d20205c 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -97,6 +97,7 @@ class Neighbor : protected Pointers { double *cuttypesq; // cuttype squared double triggersq; // trigger = build when atom moves this dist + int cluster_check; // 1 if check bond/angle/etc satisfies minimg double **xhold; // atom coords at last neighbor build int maxhold; // size of xhold array @@ -255,14 +256,17 @@ class Neighbor : protected Pointers { BondPtr bond_build; // ptr to bond list functions void bond_all(); // bond list with all bonds void bond_partial(); // exclude certain bonds + void bond_check(); BondPtr angle_build; // ptr to angle list functions void angle_all(); // angle list with all angles void angle_partial(); // exclude certain angles + void angle_check(); BondPtr dihedral_build; // ptr to dihedral list functions void dihedral_all(); // dihedral list with all dihedrals void dihedral_partial(); // exclude certain dihedrals + void dihedral_check(int **); BondPtr improper_build; // ptr to improper list functions void improper_all(); // improper list with all impropers