From 0f4be10d4a06eb2432c4c4d91c80fc43cab8d19e Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 14 Jun 2012 22:17:17 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8289 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/domain.cpp | 67 ++++++++++++++++++++++++++++++++++++++++ src/domain.h | 3 +- src/neigh_full.cpp | 26 +++++++++++++--- src/neigh_half_bin.cpp | 23 +++++++++++--- src/neigh_half_multi.cpp | 21 ++++++++++--- src/neigh_half_nsq.cpp | 11 +++++-- src/neigh_respa.cpp | 48 ++++++++++++++++++++++++---- 7 files changed, 177 insertions(+), 22 deletions(-) diff --git a/src/domain.cpp b/src/domain.cpp index 96f3772c21..472da91734 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -41,6 +41,7 @@ using namespace MathConst; #define BIG 1.0e20 #define SMALL 1.0e-4 #define DELTA 1 +#define BONDSTRETCH 1.1 enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp @@ -97,6 +98,12 @@ Domain::~Domain() void Domain::init() { + // check for too small a periodic box for molecular system + + if (atom->molecular && box_too_small()) + error->all(FLERR,"Bond/angle/dihedral extent must be < " + "half of periodic box dimension"); + // set box_change if box dimensions/shape ever changes // due to shrink-wrapping, fixes that change volume (npt, vol/rescale, etc) @@ -515,6 +522,66 @@ void Domain::pbc() } } +/* ---------------------------------------------------------------------- + check that no pair of atoms in a bonded interaction + are further apart than half a periodic box length + return 1 if any pair is, else 0 +------------------------------------------------------------------------- */ + +int Domain::box_too_small() +{ + int i,j,k; + + // only need to check if some dimension is periodic + + if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return 0; + + // maxbondall = longest current bond length + // NOTE: if box is tiny (less than 2 * bond-length), + // the check itself may compute bad bond lengths + // not sure how to account for that extreme case + + int *num_bond = atom->num_bond; + int **bond_atom = atom->bond_atom; + double **x = atom->x; + int nlocal = atom->nlocal; + + double delx,dely,delz,rsq,r; + double maxbondme = 0.0; + + for (i = 0; i < nlocal; i++) + for (j = 0; j < num_bond[i]; j++) { + k = atom->map(bond_atom[i][j]); + if (k < 0) error->one(FLERR,"Bond atom missing in box size check"); + delx = x[i][0] - x[k][0]; + dely = x[i][1] - x[k][1]; + delz = x[i][2] - x[k][2]; + domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + maxbondme = MAX(maxbondme,rsq); + } + + double maxbondall; + MPI_Allreduce(&maxbondme,&maxbondall,1,MPI_DOUBLE,MPI_MAX,world); + maxbondall = sqrt(maxbondall); + + // maxdelta = furthest apart 2 atoms in a bonded interaction can be + // include BONDSTRETCH factor to account for dynamics + + double maxdelta = maxbondall * BONDSTRETCH; + if (atom->nangles) maxdelta = 2.0 * maxbondall * BONDSTRETCH; + if (atom->ndihedrals) maxdelta = 3.0 * maxbondall * BONDSTRETCH; + + // maxdelta cannot be more than half a periodic box length, + // else when use minimg() in bond/angle/dihdral compute, + // could calculate incorrect distance between 2 atoms + + if (xperiodic && maxdelta > xprd_half) return 1; + if (yperiodic && maxdelta > yprd_half) return 1; + if (dimension == 3 && zperiodic && maxdelta > zprd_half) return 1; + return 0; +} + /* ---------------------------------------------------------------------- minimum image convention check return 1 if any distance > 1/2 of box size diff --git a/src/domain.h b/src/domain.h index 5922908071..58303527c0 100644 --- a/src/domain.h +++ b/src/domain.h @@ -70,7 +70,7 @@ class Domain : protected Pointers { // triclinic box double xy,xz,yz; // 3 tilt factors - double h[6],h_inv[6]; // shape matrix in Voigt notation + double h[6],h_inv[6]; // shape matrix in Voigt notation double h_rate[6],h_ratelo[3]; // rate of box size/shape change int box_change; // 1 if box bounds ever change, 0 if fixed @@ -93,6 +93,7 @@ class Domain : protected Pointers { virtual void set_local_box(); virtual void reset_box(); virtual void pbc(); + int box_too_small(); int minimum_image_check(double, double, double); void minimum_image(double &, double &, double &); void minimum_image(double *); diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index 15ca926340..98812d077f 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -14,6 +14,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" #include "group.h" #include "error.h" @@ -89,7 +90,10 @@ void Neighbor::full_nsq(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -173,7 +177,10 @@ void Neighbor::full_nsq_ghost(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -281,7 +288,10 @@ void Neighbor::full_bin(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -379,7 +389,10 @@ void Neighbor::full_bin_ghost(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -507,7 +520,10 @@ void Neighbor::full_multi(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index 09a5ccab92..aa748c5f7e 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -14,6 +14,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" #include "error.h" using namespace LAMMPS_NS; @@ -98,7 +99,10 @@ void Neighbor::half_bin_no_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -196,7 +200,11 @@ void Neighbor::half_bin_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -217,7 +225,11 @@ void Neighbor::half_bin_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -320,7 +332,10 @@ void Neighbor::half_bin_newton_tri(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp index 917699c444..bee0bfffcd 100644 --- a/src/neigh_half_multi.cpp +++ b/src/neigh_half_multi.cpp @@ -14,6 +14,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" #include "error.h" using namespace LAMMPS_NS; @@ -106,7 +107,10 @@ void Neighbor::half_multi_no_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -207,7 +211,10 @@ void Neighbor::half_multi_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -235,7 +242,10 @@ void Neighbor::half_multi_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -348,7 +358,10 @@ void Neighbor::half_multi_newton_tri(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp index 22f0952d15..98a8893c64 100644 --- a/src/neigh_half_nsq.cpp +++ b/src/neigh_half_nsq.cpp @@ -14,6 +14,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" #include "group.h" #include "error.h" @@ -89,7 +90,10 @@ void Neighbor::half_nsq_no_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } @@ -193,7 +197,10 @@ void Neighbor::half_nsq_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; } } diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp index 6770b26a51..c152eba2b8 100644 --- a/src/neigh_respa.cpp +++ b/src/neigh_respa.cpp @@ -14,6 +14,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "atom.h" +#include "domain.h" #include "group.h" #include "error.h" @@ -81,6 +82,7 @@ void Neighbor::respa_nsq_no_newton(NeighList *list) int npnt_middle = 0; int which = 0; + int minchange = 0; for (i = 0; i < nlocal; i++) { @@ -132,16 +134,21 @@ void Neighbor::respa_nsq_no_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); } @@ -242,6 +249,7 @@ void Neighbor::respa_nsq_newton(NeighList *list) int npnt_middle = 0; int which = 0; + int minchange = 0; for (i = 0; i < nlocal; i++) { @@ -310,17 +318,22 @@ void Neighbor::respa_nsq_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); } @@ -423,6 +436,7 @@ void Neighbor::respa_bin_no_newton(NeighList *list) int npnt_middle = 0; int which = 0; + int minchange = 0; for (i = 0; i < nlocal; i++) { @@ -480,11 +494,15 @@ void Neighbor::respa_bin_no_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } @@ -492,6 +510,7 @@ void Neighbor::respa_bin_no_newton(NeighList *list) if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); } @@ -594,6 +613,7 @@ void Neighbor::respa_bin_newton(NeighList *list) int npnt_middle = 0; int which = 0; + int minchange = 0; for (i = 0; i < nlocal; i++) { @@ -654,17 +674,22 @@ void Neighbor::respa_bin_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); } @@ -687,11 +712,15 @@ void Neighbor::respa_bin_newton(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } @@ -699,6 +728,7 @@ void Neighbor::respa_bin_newton(NeighList *list) if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); } @@ -801,6 +831,7 @@ void Neighbor::respa_bin_newton_tri(NeighList *list) int npnt_middle = 0; int which = 0; + int minchange = 0; for (i = 0; i < nlocal; i++) { @@ -866,11 +897,15 @@ void Neighbor::respa_bin_newton_tri(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { which = find_special(special[i],nspecial[i],tag[j]); - if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + if (which == 0) neighptr[n++] = j; + else if (minchange = domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; + else if (minchange) neighptr_inner[n_inner++] = j; else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } @@ -878,6 +913,7 @@ void Neighbor::respa_bin_newton_tri(NeighList *list) if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; + else if (minchange) neighptr_middle[n_middle++] = j; else if (which > 0) neighptr_middle[n_middle++] = j ^ (which << SBBITS); }