diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 7da19241e9..ddafd33d84 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -718,6 +718,8 @@ void FixGCMC::attempt_atomic_translation() coord[2] = x[i][2] + displace*rz; } } + if (!domain->inside(coord)) + error->one(FLERR,"Fix gcmc put atom outside box"); double energy_after = energy(i,ngcmc_type,-1,coord); if (random_unequal->uniform() < @@ -805,6 +807,9 @@ void FixGCMC::attempt_atomic_insertion() coord[2] = zlo + random_equal->uniform() * (zhi-zlo); } + if (!domain->inside(coord)) + error->one(FLERR,"Fix gcmc put atom outside box"); + int proc_flag = 0; if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && @@ -919,6 +924,8 @@ void FixGCMC::attempt_molecule_translation() coord[0] = x[i][0] + com_displace[0]; coord[1] = x[i][1] + com_displace[1]; coord[2] = x[i][2] + com_displace[2]; + if (!domain->inside(coord)) + error->one(FLERR,"Fix gcmc put atom outside box"); energy_after += energy(i,atom->type[i],translation_molecule,coord); } } @@ -1002,6 +1009,8 @@ void FixGCMC::attempt_molecule_rotation() xtmp[1] = atom_coord[n][1]; xtmp[2] = atom_coord[n][2]; domain->remap(xtmp); + if (!domain->inside(xtmp)) + error->one(FLERR,"Fix gcmc put atom outside box"); energy_after += energy(i,atom->type[i],rotation_molecule,xtmp); n++; } @@ -1121,6 +1130,8 @@ void FixGCMC::attempt_molecule_insertion() xtmp[1] = atom_coord[i][1]; xtmp[2] = atom_coord[i][2]; domain->remap(xtmp); + if (!domain->inside(xtmp)) + error->one(FLERR,"Fix gcmc put atom outside box"); procflag[i] = false; if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] && @@ -1313,6 +1324,8 @@ void FixGCMC::attempt_atomic_translation_full() coord[2] = x[i][2] + displace*rz; } } + if (!domain->inside(coord)) + error->one(FLERR,"Fix gcmc put atom outside box"); xtmp[0] = x[i][0]; xtmp[1] = x[i][1]; xtmp[2] = x[i][2]; @@ -1426,6 +1439,9 @@ void FixGCMC::attempt_atomic_insertion_full() coord[2] = zlo + random_equal->uniform() * (zhi-zlo); } + if (!domain->inside(coord)) + error->one(FLERR,"Fix gcmc put atom outside box"); + int proc_flag = 0; if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && @@ -1540,6 +1556,8 @@ void FixGCMC::attempt_molecule_translation_full() x[i][0] += com_displace[0]; x[i][1] += com_displace[1]; x[i][2] += com_displace[2]; + if (!domain->inside(x[i])) + error->one(FLERR,"Fix gcmc put atom outside box"); } } @@ -1621,6 +1639,8 @@ void FixGCMC::attempt_molecule_rotation_full() x[i][2] += com[2]; image[i] = imagetmp; domain->remap(x[i],image[i]); + if (!domain->inside(x[i])) + error->one(FLERR,"Fix gcmc put atom outside box"); n++; } } @@ -1787,6 +1807,8 @@ void FixGCMC::attempt_molecule_insertion_full() xtmp[2] += com_coord[2]; domain->remap(xtmp); + if (!domain->inside(xtmp)) + error->one(FLERR,"Fix gcmc put atom outside box"); if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] && xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] && diff --git a/src/domain.cpp b/src/domain.cpp index da93f50720..134065a506 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -583,6 +583,40 @@ void Domain::pbc() } } +/* ---------------------------------------------------------------------- + check that point is inside box boundaries, in [lo,hi) sense + return 1 if true, 0 if false +------------------------------------------------------------------------- */ + +int Domain::inside(double* x) +{ + double *lo,*hi,*period; + double delta[3]; + + if (triclinic == 0) { + lo = boxlo; + hi = boxhi; + period = prd; + } else { + lo = boxlo_lamda; + hi = boxhi_lamda; + period = prd_lamda; + + delta[0] = x[0] - boxlo[0]; + delta[1] = x[1] - boxlo[1]; + delta[2] = x[2] - boxlo[2]; + + x[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; + x[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; + x[2] = h_inv[2]*delta[2]; + } + + if (x[0] < lo[0] || x[0] >= hi[0] || + x[1] < lo[1] || x[1] >= hi[1] || + x[2] < lo[2] || x[2] >= hi[2]) return 0; + else return 1; +} + /* ---------------------------------------------------------------------- warn if image flags of any bonded atoms are inconsistent could be a problem when using replicate or fix rigid diff --git a/src/domain.h b/src/domain.h index 12eb9aacef..d6373ecfb1 100644 --- a/src/domain.h +++ b/src/domain.h @@ -126,6 +126,7 @@ class Domain : protected Pointers { virtual void x2lamda(int); virtual void lamda2x(double *, double *); virtual void x2lamda(double *, double *); + int inside(double *); void x2lamda(double *, double *, double *, double *); void bbox(double *, double *, double *, double *); void box_corners();