diff --git a/src/CLASS2/angle_class2.cpp b/src/CLASS2/angle_class2.cpp index 4129216ed3..8b0ba3c232 100644 --- a/src/CLASS2/angle_class2.cpp +++ b/src/CLASS2/angle_class2.cpp @@ -97,7 +97,6 @@ void AngleClass2::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -107,7 +106,6 @@ void AngleClass2::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/CLASS2/bond_class2.cpp b/src/CLASS2/bond_class2.cpp index 2234227abe..08c9ee3e90 100644 --- a/src/CLASS2/bond_class2.cpp +++ b/src/CLASS2/bond_class2.cpp @@ -72,7 +72,6 @@ void BondClass2::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/CLASS2/dihedral_class2.cpp b/src/CLASS2/dihedral_class2.cpp index dcf204a37d..ea13d09bcc 100644 --- a/src/CLASS2/dihedral_class2.cpp +++ b/src/CLASS2/dihedral_class2.cpp @@ -135,26 +135,22 @@ void DihedralClass2::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // distances diff --git a/src/CLASS2/improper_class2.cpp b/src/CLASS2/improper_class2.cpp index 5c05493840..4d06ba9333 100644 --- a/src/CLASS2/improper_class2.cpp +++ b/src/CLASS2/improper_class2.cpp @@ -118,17 +118,14 @@ void ImproperClass2::compute(int eflag, int vflag) delr[0][0] = x[i1][0] - x[i2][0]; delr[0][1] = x[i1][1] - x[i2][1]; delr[0][2] = x[i1][2] - x[i2][2]; - domain->minimum_image(delr[0]); delr[1][0] = x[i3][0] - x[i2][0]; delr[1][1] = x[i3][1] - x[i2][1]; delr[1][2] = x[i3][2] - x[i2][2]; - domain->minimum_image(delr[1]); delr[2][0] = x[i4][0] - x[i2][0]; delr[2][1] = x[i4][1] - x[i2][1]; delr[2][2] = x[i4][2] - x[i2][2]; - domain->minimum_image(delr[2]); // bond lengths and associated values @@ -661,17 +658,14 @@ void ImproperClass2::angleangle(int eflag, int vflag) delxAB = x[i1][0] - x[i2][0]; delyAB = x[i1][1] - x[i2][1]; delzAB = x[i1][2] - x[i2][2]; - domain->minimum_image(delxAB,delyAB,delzAB); delxBC = x[i3][0] - x[i2][0]; delyBC = x[i3][1] - x[i2][1]; delzBC = x[i3][2] - x[i2][2]; - domain->minimum_image(delxBC,delyBC,delzBC); delxBD = x[i4][0] - x[i2][0]; delyBD = x[i4][1] - x[i2][1]; delzBD = x[i4][2] - x[i2][2]; - domain->minimum_image(delxBD,delyBD,delzBD); // bond lengths diff --git a/src/MOLECULE/angle_charmm.cpp b/src/MOLECULE/angle_charmm.cpp index 8d086838fd..9c9a86179d 100644 --- a/src/MOLECULE/angle_charmm.cpp +++ b/src/MOLECULE/angle_charmm.cpp @@ -82,7 +82,6 @@ void AngleCharmm::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -92,7 +91,6 @@ void AngleCharmm::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -102,7 +100,6 @@ void AngleCharmm::compute(int eflag, int vflag) delxUB = x[i3][0] - x[i1][0]; delyUB = x[i3][1] - x[i1][1]; delzUB = x[i3][2] - x[i1][2]; - domain->minimum_image(delxUB,delyUB,delzUB); rsqUB = delxUB*delxUB + delyUB*delyUB + delzUB*delzUB; rUB = sqrt(rsqUB); diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp index d6bd3965d7..6ff6985277 100644 --- a/src/MOLECULE/angle_cosine.cpp +++ b/src/MOLECULE/angle_cosine.cpp @@ -73,7 +73,6 @@ void AngleCosine::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -83,7 +82,6 @@ void AngleCosine::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/angle_cosine_delta.cpp b/src/MOLECULE/angle_cosine_delta.cpp index b7cceb6186..08124142e5 100644 --- a/src/MOLECULE/angle_cosine_delta.cpp +++ b/src/MOLECULE/angle_cosine_delta.cpp @@ -65,7 +65,6 @@ void AngleCosineDelta::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -75,7 +74,6 @@ void AngleCosineDelta::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp index 7c676d309b..39ab78114f 100644 --- a/src/MOLECULE/angle_cosine_periodic.cpp +++ b/src/MOLECULE/angle_cosine_periodic.cpp @@ -80,7 +80,6 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -90,7 +89,6 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/angle_cosine_squared.cpp b/src/MOLECULE/angle_cosine_squared.cpp index f9275589e8..9eaa4c31af 100644 --- a/src/MOLECULE/angle_cosine_squared.cpp +++ b/src/MOLECULE/angle_cosine_squared.cpp @@ -79,7 +79,6 @@ void AngleCosineSquared::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -89,7 +88,6 @@ void AngleCosineSquared::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp index a1d36aca6e..87a733ba14 100644 --- a/src/MOLECULE/angle_harmonic.cpp +++ b/src/MOLECULE/angle_harmonic.cpp @@ -75,7 +75,6 @@ void AngleHarmonic::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -85,7 +84,6 @@ void AngleHarmonic::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index a5d16910d1..0be817acc7 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -91,7 +91,6 @@ void AngleTable::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -101,7 +100,6 @@ void AngleTable::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp index 27880ee1e3..41d81f5e17 100644 --- a/src/MOLECULE/bond_fene.cpp +++ b/src/MOLECULE/bond_fene.cpp @@ -72,7 +72,6 @@ void BondFENE::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); // force from log term diff --git a/src/MOLECULE/bond_fene_expand.cpp b/src/MOLECULE/bond_fene_expand.cpp index 6b80478e51..a8051468ac 100644 --- a/src/MOLECULE/bond_fene_expand.cpp +++ b/src/MOLECULE/bond_fene_expand.cpp @@ -74,7 +74,6 @@ void BondFENEExpand::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); // force from log term diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp index ab67e8bec6..4fb18d13ca 100644 --- a/src/MOLECULE/bond_harmonic.cpp +++ b/src/MOLECULE/bond_harmonic.cpp @@ -66,7 +66,6 @@ void BondHarmonic::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/MOLECULE/bond_morse.cpp b/src/MOLECULE/bond_morse.cpp index 070c301136..22d2909244 100644 --- a/src/MOLECULE/bond_morse.cpp +++ b/src/MOLECULE/bond_morse.cpp @@ -71,7 +71,6 @@ void BondMorse::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/MOLECULE/bond_nonlinear.cpp b/src/MOLECULE/bond_nonlinear.cpp index 5bc12e650d..6577d4380c 100644 --- a/src/MOLECULE/bond_nonlinear.cpp +++ b/src/MOLECULE/bond_nonlinear.cpp @@ -67,7 +67,6 @@ void BondNonlinear::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index 8b5264bacd..9c49a5e510 100755 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -89,7 +89,6 @@ void BondQuartic::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index 19b3f66402..1f83c7980f 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -83,7 +83,6 @@ void BondTable::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/MOLECULE/dihedral_charmm.cpp b/src/MOLECULE/dihedral_charmm.cpp index 644655f589..6c07195939 100644 --- a/src/MOLECULE/dihedral_charmm.cpp +++ b/src/MOLECULE/dihedral_charmm.cpp @@ -101,26 +101,22 @@ void DihedralCharmm::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); ax = vb1y*vb2zm - vb1z*vb2ym; ay = vb1z*vb2xm - vb1x*vb2zm; @@ -269,7 +265,6 @@ void DihedralCharmm::compute(int eflag, int vflag) delx = x[i1][0] - x[i4][0]; dely = x[i1][1] - x[i4][1]; delz = x[i1][2] - x[i4][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; diff --git a/src/MOLECULE/dihedral_harmonic.cpp b/src/MOLECULE/dihedral_harmonic.cpp index a951f3816b..6f2e887f74 100644 --- a/src/MOLECULE/dihedral_harmonic.cpp +++ b/src/MOLECULE/dihedral_harmonic.cpp @@ -87,26 +87,22 @@ void DihedralHarmonic::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c,s calculation diff --git a/src/MOLECULE/dihedral_helix.cpp b/src/MOLECULE/dihedral_helix.cpp index e9c17392d4..04bc7b029c 100644 --- a/src/MOLECULE/dihedral_helix.cpp +++ b/src/MOLECULE/dihedral_helix.cpp @@ -90,26 +90,22 @@ void DihedralHelix::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/MOLECULE/dihedral_multi_harmonic.cpp b/src/MOLECULE/dihedral_multi_harmonic.cpp index cff19447f7..d004a4cd94 100644 --- a/src/MOLECULE/dihedral_multi_harmonic.cpp +++ b/src/MOLECULE/dihedral_multi_harmonic.cpp @@ -86,26 +86,22 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/MOLECULE/dihedral_opls.cpp b/src/MOLECULE/dihedral_opls.cpp index 676434a061..0e40843fa9 100644 --- a/src/MOLECULE/dihedral_opls.cpp +++ b/src/MOLECULE/dihedral_opls.cpp @@ -86,26 +86,22 @@ void DihedralOPLS::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/MOLECULE/improper_cvff.cpp b/src/MOLECULE/improper_cvff.cpp index a75d6c077d..a0aa12832e 100644 --- a/src/MOLECULE/improper_cvff.cpp +++ b/src/MOLECULE/improper_cvff.cpp @@ -81,26 +81,22 @@ void ImproperCvff::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/MOLECULE/improper_harmonic.cpp b/src/MOLECULE/improper_harmonic.cpp index dae70277c9..da6f50bd84 100644 --- a/src/MOLECULE/improper_harmonic.cpp +++ b/src/MOLECULE/improper_harmonic.cpp @@ -81,17 +81,14 @@ void ImproperHarmonic::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp index 61055a9598..939566a24b 100644 --- a/src/MOLECULE/improper_umbrella.cpp +++ b/src/MOLECULE/improper_umbrella.cpp @@ -84,21 +84,18 @@ void ImproperUmbrella::compute(int eflag, int vflag) vb1x = x[i2][0] - x[i1][0]; vb1y = x[i2][1] - x[i1][1]; vb1z = x[i2][2] - x[i1][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i1][0]; vb2y = x[i3][1] - x[i1][1]; vb2z = x[i3][2] - x[i1][2]; - domain->minimum_image(vb2x,vb2y,vb2z); // 3rd bond vb3x = x[i4][0] - x[i1][0]; vb3y = x[i4][1] - x[i1][1]; vb3z = x[i4][2] - x[i1][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation // A = vb1 X vb2 is perpendicular to IJK plane diff --git a/src/USER-CG-CMM/angle_cg_cmm.cpp b/src/USER-CG-CMM/angle_cg_cmm.cpp index f147a87237..dbeb04f50c 100644 --- a/src/USER-CG-CMM/angle_cg_cmm.cpp +++ b/src/USER-CG-CMM/angle_cg_cmm.cpp @@ -166,7 +166,6 @@ void AngleCGCMM::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -176,7 +175,6 @@ void AngleCGCMM::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -200,7 +198,6 @@ void AngleCGCMM::compute(int eflag, int vflag) delx3 = x[i1][0] - x[i3][0]; dely3 = x[i1][1] - x[i3][1]; delz3 = x[i1][2] - x[i3][2]; - domain->minimum_image(delx3,dely3,delz3); rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3; r3 = sqrt(rsq3); diff --git a/src/USER-CG-CMM/angle_sdk.cpp b/src/USER-CG-CMM/angle_sdk.cpp index 35c4416af1..fef80ea832 100644 --- a/src/USER-CG-CMM/angle_sdk.cpp +++ b/src/USER-CG-CMM/angle_sdk.cpp @@ -88,7 +88,6 @@ void AngleSDK::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -98,7 +97,6 @@ void AngleSDK::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -128,7 +126,6 @@ void AngleSDK::compute(int eflag, int vflag) delx3 = x[i1][0] - x[i3][0]; dely3 = x[i1][1] - x[i3][1]; delz3 = x[i1][2] - x[i3][2]; - domain->minimum_image(delx3,dely3,delz3); rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3; const int type1 = atom->type[i1]; diff --git a/src/USER-MISC/angle_cosine_shift.cpp b/src/USER-MISC/angle_cosine_shift.cpp index 11b4cd070c..edeececffb 100644 --- a/src/USER-MISC/angle_cosine_shift.cpp +++ b/src/USER-MISC/angle_cosine_shift.cpp @@ -78,7 +78,6 @@ void AngleCosineShift::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -88,7 +87,6 @@ void AngleCosineShift::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-MISC/angle_cosine_shift_exp.cpp b/src/USER-MISC/angle_cosine_shift_exp.cpp index 0c607b4c43..db75fb0b23 100644 --- a/src/USER-MISC/angle_cosine_shift_exp.cpp +++ b/src/USER-MISC/angle_cosine_shift_exp.cpp @@ -82,7 +82,6 @@ void AngleCosineShiftExp::compute(int eflag, int vflag) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -92,7 +91,6 @@ void AngleCosineShiftExp::compute(int eflag, int vflag) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-MISC/angle_dipole.cpp b/src/USER-MISC/angle_dipole.cpp index 4020842413..b80007e50d 100644 --- a/src/USER-MISC/angle_dipole.cpp +++ b/src/USER-MISC/angle_dipole.cpp @@ -78,7 +78,6 @@ void AngleDipole::compute(int eflag, int vflag) delx = x[iRef][0] - x[iDip][0]; dely = x[iRef][1] - x[iDip][1]; delz = x[iRef][2] - x[iDip][2]; - domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); diff --git a/src/USER-MISC/bond_harmonic_shift.cpp b/src/USER-MISC/bond_harmonic_shift.cpp index e140892141..21b0ff40ab 100644 --- a/src/USER-MISC/bond_harmonic_shift.cpp +++ b/src/USER-MISC/bond_harmonic_shift.cpp @@ -71,7 +71,6 @@ void BondHarmonicShift::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); @@ -84,7 +83,8 @@ void BondHarmonicShift::compute(int eflag, int vflag) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (eflag) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); + if (eflag) + ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-MISC/bond_harmonic_shift_cut.cpp b/src/USER-MISC/bond_harmonic_shift_cut.cpp index 8ebb820e0b..d867b5e0cf 100644 --- a/src/USER-MISC/bond_harmonic_shift_cut.cpp +++ b/src/USER-MISC/bond_harmonic_shift_cut.cpp @@ -71,7 +71,6 @@ void BondHarmonicShiftCut::compute(int eflag, int vflag) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); @@ -86,7 +85,8 @@ void BondHarmonicShiftCut::compute(int eflag, int vflag) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (eflag) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); + if (eflag) + ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type])); // apply force to each of 2 atoms diff --git a/src/USER-MISC/dihedral_cosine_shift_exp.cpp b/src/USER-MISC/dihedral_cosine_shift_exp.cpp index d607051d9f..0491f87c06 100644 --- a/src/USER-MISC/dihedral_cosine_shift_exp.cpp +++ b/src/USER-MISC/dihedral_cosine_shift_exp.cpp @@ -90,26 +90,22 @@ void DihedralCosineShiftExp::compute(int eflag, int vflag) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c,s calculation diff --git a/src/USER-MISC/improper_cossq.cpp b/src/USER-MISC/improper_cossq.cpp index 9bd5be609f..1614deaed1 100644 --- a/src/USER-MISC/improper_cossq.cpp +++ b/src/USER-MISC/improper_cossq.cpp @@ -86,7 +86,6 @@ void ImproperCossq::compute(int eflag, int vflag) vb1x = x[i2][0] - x[i1][0]; vb1y = x[i2][1] - x[i1][1]; vb1z = x[i2][2] - x[i1][2]; - domain->minimum_image(vb1x,vb1y,vb1z); rjisq = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z ; rji = sqrt(rjisq); @@ -94,13 +93,11 @@ void ImproperCossq::compute(int eflag, int vflag) vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); /* separation vector between i3 and i4, (i4-i3) */ vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); rlksq = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z ; rlk = sqrt(rlksq); diff --git a/src/USER-MISC/improper_ring.cpp b/src/USER-MISC/improper_ring.cpp index 0afbb197ca..fbbe9dd75b 100644 --- a/src/USER-MISC/improper_ring.cpp +++ b/src/USER-MISC/improper_ring.cpp @@ -123,13 +123,10 @@ void ImproperRing::compute(int eflag, int vflag) Although, they are irrelevant to the calculation of the potential, we keep them for maximal compatibility. */ vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); /* Pass the atom tags to form the necessary combinations. */ @@ -148,7 +145,6 @@ void ImproperRing::compute(int eflag, int vflag) bvec1x[icomb] = x[at2[icomb]][0] - x[at1[icomb]][0]; bvec1y[icomb] = x[at2[icomb]][1] - x[at1[icomb]][1]; bvec1z[icomb] = x[at2[icomb]][2] - x[at1[icomb]][2]; - domain -> minimum_image(bvec1x[icomb], bvec1y[icomb], bvec1z[icomb]); /* also calculate the norm of the vector: */ bvec1n[icomb] = sqrt( bvec1x[icomb]*bvec1x[icomb] + bvec1y[icomb]*bvec1y[icomb] @@ -157,7 +153,6 @@ void ImproperRing::compute(int eflag, int vflag) bvec2x[icomb] = x[at3[icomb]][0] - x[at2[icomb]][0]; bvec2y[icomb] = x[at3[icomb]][1] - x[at2[icomb]][1]; bvec2z[icomb] = x[at3[icomb]][2] - x[at2[icomb]][2]; - domain -> minimum_image(bvec2x[icomb], bvec2y[icomb], bvec2z[icomb]); /* also calculate the norm of the vector: */ bvec2n[icomb] = sqrt( bvec2x[icomb]*bvec2x[icomb] + bvec2y[icomb]*bvec2y[icomb] diff --git a/src/USER-OMP/angle_charmm_omp.cpp b/src/USER-OMP/angle_charmm_omp.cpp index 1efe62d4d3..4cc86dd11c 100644 --- a/src/USER-OMP/angle_charmm_omp.cpp +++ b/src/USER-OMP/angle_charmm_omp.cpp @@ -106,7 +106,6 @@ void AngleCharmmOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -116,7 +115,6 @@ void AngleCharmmOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -126,7 +124,6 @@ void AngleCharmmOMP::eval(int nfrom, int nto, ThrData * const thr) delxUB = x[i3][0] - x[i1][0]; delyUB = x[i3][1] - x[i1][1]; delzUB = x[i3][2] - x[i1][2]; - domain->minimum_image(delxUB,delyUB,delzUB); rsqUB = delxUB*delxUB + delyUB*delyUB + delzUB*delzUB; rUB = sqrt(rsqUB); diff --git a/src/USER-OMP/angle_class2_omp.cpp b/src/USER-OMP/angle_class2_omp.cpp index 0dbd0cbfda..ef887a93e2 100644 --- a/src/USER-OMP/angle_class2_omp.cpp +++ b/src/USER-OMP/angle_class2_omp.cpp @@ -107,7 +107,6 @@ void AngleClass2OMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -117,7 +116,6 @@ void AngleClass2OMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_delta_omp.cpp b/src/USER-OMP/angle_cosine_delta_omp.cpp index 5ba57d20b5..caf259b15c 100644 --- a/src/USER-OMP/angle_cosine_delta_omp.cpp +++ b/src/USER-OMP/angle_cosine_delta_omp.cpp @@ -104,7 +104,6 @@ void AngleCosineDeltaOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -114,7 +113,6 @@ void AngleCosineDeltaOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_omp.cpp b/src/USER-OMP/angle_cosine_omp.cpp index a773f94c50..51d3494a4f 100644 --- a/src/USER-OMP/angle_cosine_omp.cpp +++ b/src/USER-OMP/angle_cosine_omp.cpp @@ -104,7 +104,6 @@ void AngleCosineOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -114,7 +113,6 @@ void AngleCosineOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_periodic_omp.cpp b/src/USER-OMP/angle_cosine_periodic_omp.cpp index a86b68611f..6c5932e70c 100644 --- a/src/USER-OMP/angle_cosine_periodic_omp.cpp +++ b/src/USER-OMP/angle_cosine_periodic_omp.cpp @@ -105,7 +105,6 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_shift_exp_omp.cpp b/src/USER-OMP/angle_cosine_shift_exp_omp.cpp index 81cd8ecad3..18af91cf7d 100644 --- a/src/USER-OMP/angle_cosine_shift_exp_omp.cpp +++ b/src/USER-OMP/angle_cosine_shift_exp_omp.cpp @@ -105,7 +105,6 @@ void AngleCosineShiftExpOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleCosineShiftExpOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_shift_omp.cpp b/src/USER-OMP/angle_cosine_shift_omp.cpp index 2c830cf393..a7f40fee8a 100644 --- a/src/USER-OMP/angle_cosine_shift_omp.cpp +++ b/src/USER-OMP/angle_cosine_shift_omp.cpp @@ -105,7 +105,6 @@ void AngleCosineShiftOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleCosineShiftOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_cosine_squared_omp.cpp b/src/USER-OMP/angle_cosine_squared_omp.cpp index ccdaf2bdbc..e8eaff6b2a 100644 --- a/src/USER-OMP/angle_cosine_squared_omp.cpp +++ b/src/USER-OMP/angle_cosine_squared_omp.cpp @@ -105,7 +105,6 @@ void AngleCosineSquaredOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleCosineSquaredOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_dipole_omp.cpp b/src/USER-OMP/angle_dipole_omp.cpp index f9cbc6f51d..025ffcd16c 100644 --- a/src/USER-OMP/angle_dipole_omp.cpp +++ b/src/USER-OMP/angle_dipole_omp.cpp @@ -102,7 +102,6 @@ void AngleDipoleOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[iRef][0] - x[iDip][0]; dely = x[iRef][1] - x[iDip][1]; delz = x[iRef][2] - x[iDip][2]; - domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); diff --git a/src/USER-OMP/angle_harmonic_omp.cpp b/src/USER-OMP/angle_harmonic_omp.cpp index b6ff1a1ddd..c30f830b77 100644 --- a/src/USER-OMP/angle_harmonic_omp.cpp +++ b/src/USER-OMP/angle_harmonic_omp.cpp @@ -105,7 +105,6 @@ void AngleHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/angle_sdk_omp.cpp b/src/USER-OMP/angle_sdk_omp.cpp index 43e0049e16..bd780e0fe0 100644 --- a/src/USER-OMP/angle_sdk_omp.cpp +++ b/src/USER-OMP/angle_sdk_omp.cpp @@ -107,7 +107,6 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -117,7 +116,6 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); @@ -147,7 +145,6 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr) delx3 = x[i1][0] - x[i3][0]; dely3 = x[i1][1] - x[i3][1]; delz3 = x[i1][2] - x[i3][2]; - domain->minimum_image(delx3,dely3,delz3); rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3; const int type1 = atom->type[i1]; diff --git a/src/USER-OMP/angle_table_omp.cpp b/src/USER-OMP/angle_table_omp.cpp index 8eb73b4c2b..15c32c8943 100644 --- a/src/USER-OMP/angle_table_omp.cpp +++ b/src/USER-OMP/angle_table_omp.cpp @@ -105,7 +105,6 @@ void AngleTableOMP::eval(int nfrom, int nto, ThrData * const thr) delx1 = x[i1][0] - x[i2][0]; dely1 = x[i1][1] - x[i2][1]; delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; r1 = sqrt(rsq1); @@ -115,7 +114,6 @@ void AngleTableOMP::eval(int nfrom, int nto, ThrData * const thr) delx2 = x[i3][0] - x[i2][0]; dely2 = x[i3][1] - x[i2][1]; delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; r2 = sqrt(rsq2); diff --git a/src/USER-OMP/bond_class2_omp.cpp b/src/USER-OMP/bond_class2_omp.cpp index a31ac2e3aa..580d9b9804 100644 --- a/src/USER-OMP/bond_class2_omp.cpp +++ b/src/USER-OMP/bond_class2_omp.cpp @@ -96,7 +96,6 @@ void BondClass2OMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/bond_fene_expand_omp.cpp b/src/USER-OMP/bond_fene_expand_omp.cpp index 39c43d978e..4706ef8232 100644 --- a/src/USER-OMP/bond_fene_expand_omp.cpp +++ b/src/USER-OMP/bond_fene_expand_omp.cpp @@ -99,7 +99,6 @@ void BondFENEExpandOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/bond_fene_omp.cpp b/src/USER-OMP/bond_fene_omp.cpp index 72596d1864..34f90ea3a0 100644 --- a/src/USER-OMP/bond_fene_omp.cpp +++ b/src/USER-OMP/bond_fene_omp.cpp @@ -98,7 +98,6 @@ void BondFENEOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r0sq = r0[type] * r0[type]; diff --git a/src/USER-OMP/bond_harmonic_omp.cpp b/src/USER-OMP/bond_harmonic_omp.cpp index 9c8c6d2c71..dc8bab1c54 100644 --- a/src/USER-OMP/bond_harmonic_omp.cpp +++ b/src/USER-OMP/bond_harmonic_omp.cpp @@ -95,7 +95,6 @@ void BondHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp index de2c72e890..6ea5f6c2c1 100644 --- a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp +++ b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp @@ -95,7 +95,6 @@ void BondHarmonicShiftCutOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); @@ -110,7 +109,8 @@ void BondHarmonicShiftCutOMP::eval(int nfrom, int nto, ThrData * const thr) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (EFLAG) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); + if (EFLAG) + ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-OMP/bond_harmonic_shift_omp.cpp b/src/USER-OMP/bond_harmonic_shift_omp.cpp index 2d566f771b..2f6f1d37b9 100644 --- a/src/USER-OMP/bond_harmonic_shift_omp.cpp +++ b/src/USER-OMP/bond_harmonic_shift_omp.cpp @@ -95,7 +95,6 @@ void BondHarmonicShiftOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); @@ -107,7 +106,8 @@ void BondHarmonicShiftOMP::eval(int nfrom, int nto, ThrData * const thr) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - if (EFLAG) ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); + if (EFLAG) + ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) ); // apply force to each of 2 atoms diff --git a/src/USER-OMP/bond_morse_omp.cpp b/src/USER-OMP/bond_morse_omp.cpp index 0a91857c8e..21063c636a 100644 --- a/src/USER-OMP/bond_morse_omp.cpp +++ b/src/USER-OMP/bond_morse_omp.cpp @@ -95,7 +95,6 @@ void BondMorseOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/bond_nonlinear_omp.cpp b/src/USER-OMP/bond_nonlinear_omp.cpp index 46eed47604..9e61ca3215 100644 --- a/src/USER-OMP/bond_nonlinear_omp.cpp +++ b/src/USER-OMP/bond_nonlinear_omp.cpp @@ -95,7 +95,6 @@ void BondNonlinearOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/bond_quartic_omp.cpp b/src/USER-OMP/bond_quartic_omp.cpp index d07ee9aca6..79d7621a61 100644 --- a/src/USER-OMP/bond_quartic_omp.cpp +++ b/src/USER-OMP/bond_quartic_omp.cpp @@ -109,7 +109,6 @@ void BondQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; diff --git a/src/USER-OMP/bond_table_omp.cpp b/src/USER-OMP/bond_table_omp.cpp index fefadfc44a..3184991315 100644 --- a/src/USER-OMP/bond_table_omp.cpp +++ b/src/USER-OMP/bond_table_omp.cpp @@ -96,7 +96,6 @@ void BondTableOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; delz = x[i1][2] - x[i2][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r = sqrt(rsq); diff --git a/src/USER-OMP/dihedral_charmm_omp.cpp b/src/USER-OMP/dihedral_charmm_omp.cpp index fbebfbaaac..b069af43c5 100644 --- a/src/USER-OMP/dihedral_charmm_omp.cpp +++ b/src/USER-OMP/dihedral_charmm_omp.cpp @@ -124,26 +124,22 @@ void DihedralCharmmOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c,s calculation @@ -293,7 +289,6 @@ void DihedralCharmmOMP::eval(int nfrom, int nto, ThrData * const thr) delx = x[i1][0] - x[i4][0]; dely = x[i1][1] - x[i4][1]; delz = x[i1][2] - x[i4][2]; - domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; diff --git a/src/USER-OMP/dihedral_class2_omp.cpp b/src/USER-OMP/dihedral_class2_omp.cpp index 0d150851c8..d21c2869e9 100644 --- a/src/USER-OMP/dihedral_class2_omp.cpp +++ b/src/USER-OMP/dihedral_class2_omp.cpp @@ -116,26 +116,22 @@ void DihedralClass2OMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // distances diff --git a/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp b/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp index 15f4dd97c1..ecb1a85fa4 100644 --- a/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp +++ b/src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp @@ -111,26 +111,22 @@ void DihedralCosineShiftExpOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c,s calculation diff --git a/src/USER-OMP/dihedral_harmonic_omp.cpp b/src/USER-OMP/dihedral_harmonic_omp.cpp index 11023ee36c..e3b94c8f9a 100644 --- a/src/USER-OMP/dihedral_harmonic_omp.cpp +++ b/src/USER-OMP/dihedral_harmonic_omp.cpp @@ -110,26 +110,22 @@ void DihedralHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c,s calculation diff --git a/src/USER-OMP/dihedral_helix_omp.cpp b/src/USER-OMP/dihedral_helix_omp.cpp index a869367bdb..fd7025cdb4 100644 --- a/src/USER-OMP/dihedral_helix_omp.cpp +++ b/src/USER-OMP/dihedral_helix_omp.cpp @@ -114,26 +114,22 @@ void DihedralHelixOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/USER-OMP/dihedral_multi_harmonic_omp.cpp b/src/USER-OMP/dihedral_multi_harmonic_omp.cpp index a40bdba955..a68bebfdf6 100644 --- a/src/USER-OMP/dihedral_multi_harmonic_omp.cpp +++ b/src/USER-OMP/dihedral_multi_harmonic_omp.cpp @@ -111,26 +111,22 @@ void DihedralMultiHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/USER-OMP/dihedral_opls_omp.cpp b/src/USER-OMP/dihedral_opls_omp.cpp index 30b805611b..3810dc1985 100644 --- a/src/USER-OMP/dihedral_opls_omp.cpp +++ b/src/USER-OMP/dihedral_opls_omp.cpp @@ -112,26 +112,22 @@ void DihedralOPLSOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/USER-OMP/improper_class2_omp.cpp b/src/USER-OMP/improper_class2_omp.cpp index 0a4758f55e..08f4a21e00 100644 --- a/src/USER-OMP/improper_class2_omp.cpp +++ b/src/USER-OMP/improper_class2_omp.cpp @@ -133,17 +133,14 @@ void ImproperClass2OMP::eval(int nfrom, int nto, ThrData * const thr) delr[0][0] = x[i1][0] - x[i2][0]; delr[0][1] = x[i1][1] - x[i2][1]; delr[0][2] = x[i1][2] - x[i2][2]; - domain->minimum_image(delr[0]); delr[1][0] = x[i3][0] - x[i2][0]; delr[1][1] = x[i3][1] - x[i2][1]; delr[1][2] = x[i3][2] - x[i2][2]; - domain->minimum_image(delr[1]); delr[2][0] = x[i4][0] - x[i2][0]; delr[2][1] = x[i4][1] - x[i2][1]; delr[2][2] = x[i4][2] - x[i2][2]; - domain->minimum_image(delr[2]); // bond lengths and associated values @@ -542,17 +539,14 @@ void ImproperClass2OMP::angleangle_thr(int nfrom, int nto, ThrData * const thr) delxAB = x[i1][0] - x[i2][0]; delyAB = x[i1][1] - x[i2][1]; delzAB = x[i1][2] - x[i2][2]; - domain->minimum_image(delxAB,delyAB,delzAB); delxBC = x[i3][0] - x[i2][0]; delyBC = x[i3][1] - x[i2][1]; delzBC = x[i3][2] - x[i2][2]; - domain->minimum_image(delxBC,delyBC,delzBC); delxBD = x[i4][0] - x[i2][0]; delyBD = x[i4][1] - x[i2][1]; delzBD = x[i4][2] - x[i2][2]; - domain->minimum_image(delxBD,delyBD,delzBD); // bond lengths diff --git a/src/USER-OMP/improper_cossq_omp.cpp b/src/USER-OMP/improper_cossq_omp.cpp index 1e1b53eaf7..c49cd0b3ce 100644 --- a/src/USER-OMP/improper_cossq_omp.cpp +++ b/src/USER-OMP/improper_cossq_omp.cpp @@ -106,7 +106,6 @@ void ImproperCossqOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i2][0] - x[i1][0]; vb1y = x[i2][1] - x[i1][1]; vb1z = x[i2][2] - x[i1][2]; - domain->minimum_image(vb1x,vb1y,vb1z); rjisq = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z ; rji = sqrt(rjisq); @@ -114,13 +113,11 @@ void ImproperCossqOMP::eval(int nfrom, int nto, ThrData * const thr) vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); /* separation vector between i3 and i4, (i4-i3) */ vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); rlksq = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z ; rlk = sqrt(rlksq); diff --git a/src/USER-OMP/improper_cvff_omp.cpp b/src/USER-OMP/improper_cvff_omp.cpp index fd8b4aca1d..819a22b66e 100644 --- a/src/USER-OMP/improper_cvff_omp.cpp +++ b/src/USER-OMP/improper_cvff_omp.cpp @@ -109,26 +109,22 @@ void ImproperCvffOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); // 3rd bond vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation diff --git a/src/USER-OMP/improper_harmonic_omp.cpp b/src/USER-OMP/improper_harmonic_omp.cpp index 77f1a1c596..4ac8351a80 100644 --- a/src/USER-OMP/improper_harmonic_omp.cpp +++ b/src/USER-OMP/improper_harmonic_omp.cpp @@ -108,17 +108,14 @@ void ImproperHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); diff --git a/src/USER-OMP/improper_ring_omp.cpp b/src/USER-OMP/improper_ring_omp.cpp index 6f8a73c41c..78b15a0453 100644 --- a/src/USER-OMP/improper_ring_omp.cpp +++ b/src/USER-OMP/improper_ring_omp.cpp @@ -122,13 +122,10 @@ void ImproperRingOMP::eval(int nfrom, int nto, ThrData * const thr) Although, they are irrelevant to the calculation of the potential, we keep them for maximal compatibility. */ vb1x = x[i1][0] - x[i2][0]; vb1y = x[i1][1] - x[i2][1]; vb1z = x[i1][2] - x[i2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); vb2x = x[i3][0] - x[i2][0]; vb2y = x[i3][1] - x[i2][1]; vb2z = x[i3][2] - x[i2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); vb3x = x[i4][0] - x[i3][0]; vb3y = x[i4][1] - x[i3][1]; vb3z = x[i4][2] - x[i3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); /* Pass the atom tags to form the necessary combinations. */ @@ -146,7 +143,6 @@ void ImproperRingOMP::eval(int nfrom, int nto, ThrData * const thr) bvec1x[icomb] = x[at2[icomb]][0] - x[at1[icomb]][0]; bvec1y[icomb] = x[at2[icomb]][1] - x[at1[icomb]][1]; bvec1z[icomb] = x[at2[icomb]][2] - x[at1[icomb]][2]; - domain -> minimum_image(bvec1x[icomb], bvec1y[icomb], bvec1z[icomb]); /* also calculate the norm of the vector: */ bvec1n[icomb] = sqrt( bvec1x[icomb]*bvec1x[icomb] + bvec1y[icomb]*bvec1y[icomb] @@ -155,7 +151,6 @@ void ImproperRingOMP::eval(int nfrom, int nto, ThrData * const thr) bvec2x[icomb] = x[at3[icomb]][0] - x[at2[icomb]][0]; bvec2y[icomb] = x[at3[icomb]][1] - x[at2[icomb]][1]; bvec2z[icomb] = x[at3[icomb]][2] - x[at2[icomb]][2]; - domain -> minimum_image(bvec2x[icomb], bvec2y[icomb], bvec2z[icomb]); /* also calculate the norm of the vector: */ bvec2n[icomb] = sqrt( bvec2x[icomb]*bvec2x[icomb] + bvec2y[icomb]*bvec2y[icomb] diff --git a/src/USER-OMP/improper_umbrella_omp.cpp b/src/USER-OMP/improper_umbrella_omp.cpp index ebc1fc299f..9a64da1fb1 100644 --- a/src/USER-OMP/improper_umbrella_omp.cpp +++ b/src/USER-OMP/improper_umbrella_omp.cpp @@ -107,21 +107,18 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr) vb1x = x[i2][0] - x[i1][0]; vb1y = x[i2][1] - x[i1][1]; vb1z = x[i2][2] - x[i1][2]; - domain->minimum_image(vb1x,vb1y,vb1z); // 2nd bond vb2x = x[i3][0] - x[i1][0]; vb2y = x[i3][1] - x[i1][1]; vb2z = x[i3][2] - x[i1][2]; - domain->minimum_image(vb2x,vb2y,vb2z); // 3rd bond vb3x = x[i4][0] - x[i1][0]; vb3y = x[i4][1] - x[i1][1]; vb3z = x[i4][2] - x[i1][2]; - domain->minimum_image(vb3x,vb3y,vb3z); // c0 calculation // A = vb1 X vb2 is perpendicular to IJK plane diff --git a/src/atom.cpp b/src/atom.cpp index 7d4564359c..760a22db7f 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -129,21 +129,17 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) nextra_store = 0; extra = NULL; - // default mapping values and hash table primes + // default mapping values tag_enable = 1; map_style = 0; map_tag_max = 0; map_nhash = 0; - - nprimes = 38; - primes = new int[nprimes]; - int plist[] = {5041,10007,20011,30011,40009,50021,60013,70001,80021, - 90001,100003,110017,120011,130003,140009,150001,160001, - 170003,180001,190027,200003,210011,220009,230003,240007, - 250007,260003,270001,280001,290011,300007,310019,320009, - 330017,340007,350003,362881,3628801}; - for (int i = 0; i < nprimes; i++) primes[i] = plist[i]; + + sametag = NULL; + map_array = NULL; + map_bucket = NULL; + map_hash = NULL; // default atom style = atomic @@ -236,7 +232,6 @@ Atom::~Atom() // delete mapping data structures map_delete(); - delete [] primes; } /* ---------------------------------------------------------------------- @@ -431,243 +426,6 @@ void Atom::modify_params(int narg, char **arg) } } -/* ---------------------------------------------------------------------- - allocate and initialize array or hash table for global -> local map - set map_tag_max = largest atom ID (may be larger than natoms) - for array option: - array length = 1 to largest tag of any atom - set entire array to -1 as initial values - for hash option: - map_nhash = length of hash table - map_nbucket = # of hash buckets, prime larger than map_nhash - so buckets will only be filled with 0 or 1 atoms on average -------------------------------------------------------------------------- */ - -void Atom::map_init() -{ - map_delete(); - - if (tag_enable == 0) - error->all(FLERR,"Cannot create an atom map unless atoms have IDs"); - - int max = 0; - for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); - MPI_Allreduce(&max,&map_tag_max,1,MPI_INT,MPI_MAX,world); - - if (map_style == 1) { - memory->create(map_array,map_tag_max+1,"atom:map_array"); - for (int i = 0; i <= map_tag_max; i++) map_array[i] = -1; - - } else { - - // map_nhash = max of atoms/proc or total atoms, times 2, at least 1000 - - int nper = static_cast (natoms/comm->nprocs); - map_nhash = MAX(nper,nmax); - if (map_nhash > natoms) map_nhash = static_cast (natoms); - if (comm->nprocs > 1) map_nhash *= 2; - map_nhash = MAX(map_nhash,1000); - - // map_nbucket = prime just larger than map_nhash - - int n = map_nhash/10000; - n = MIN(n,nprimes-1); - map_nbucket = primes[n]; - if (map_nbucket < map_nhash && n < nprimes-1) map_nbucket = primes[n+1]; - - // set all buckets to empty - // set hash to map_nhash in length - // put all hash entries in free list and point them to each other - - map_bucket = new int[map_nbucket]; - for (int i = 0; i < map_nbucket; i++) map_bucket[i] = -1; - - map_hash = new HashElem[map_nhash]; - map_nused = 0; - map_free = 0; - for (int i = 0; i < map_nhash; i++) map_hash[i].next = i+1; - map_hash[map_nhash-1].next = -1; - } -} - -/* ---------------------------------------------------------------------- - clear global -> local map for all of my own and ghost atoms - for hash table option: - global ID may not be in table if image atom was already cleared -------------------------------------------------------------------------- */ - -void Atom::map_clear() -{ - if (map_style == 1) { - int nall = nlocal + nghost; - for (int i = 0; i < nall; i++) map_array[tag[i]] = -1; - - } else { - int previous,global,ibucket,index; - int nall = nlocal + nghost; - for (int i = 0; i < nall; i++) { - - // search for key - // if don't find it, done - - previous = -1; - global = tag[i]; - ibucket = global % map_nbucket; - index = map_bucket[ibucket]; - while (index > -1) { - if (map_hash[index].global == global) break; - previous = index; - index = map_hash[index].next; - } - if (index == -1) continue; - - // delete the hash entry and add it to free list - // special logic if entry is 1st in the bucket - - if (previous == -1) map_bucket[ibucket] = map_hash[index].next; - else map_hash[previous].next = map_hash[index].next; - - map_hash[index].next = map_free; - map_free = index; - map_nused--; - } - } -} - -/* ---------------------------------------------------------------------- - set global -> local map for all of my own and ghost atoms - loop in reverse order so that nearby images take precedence over far ones - and owned atoms take precedence over images - this enables valid lookups of bond topology atoms - for hash table option: - if hash table too small, re-init - global ID may already be in table if image atom was set -------------------------------------------------------------------------- */ - -void Atom::map_set() -{ - if (map_style == 1) { - int nall = nlocal + nghost; - for (int i = nall-1; i >= 0 ; i--) map_array[tag[i]] = i; - - } else { - int previous,global,ibucket,index; - int nall = nlocal + nghost; - if (nall > map_nhash) map_init(); - - for (int i = nall-1; i >= 0 ; i--) { - - // search for key - // if found it, just overwrite local value with index - - previous = -1; - global = tag[i]; - ibucket = global % map_nbucket; - index = map_bucket[ibucket]; - while (index > -1) { - if (map_hash[index].global == global) break; - previous = index; - index = map_hash[index].next; - } - if (index > -1) { - map_hash[index].local = i; - continue; - } - - // take one entry from free list - // add the new global/local pair as entry at end of bucket list - // special logic if this entry is 1st in bucket - - index = map_free; - map_free = map_hash[map_free].next; - if (previous == -1) map_bucket[ibucket] = index; - else map_hash[previous].next = index; - map_hash[index].global = global; - map_hash[index].local = i; - map_hash[index].next = -1; - map_nused++; - } - } -} - -/* ---------------------------------------------------------------------- - set global to local map for one atom - for hash table option: - global ID may already be in table if atom was already set -------------------------------------------------------------------------- */ - -void Atom::map_one(int global, int local) -{ - if (map_style == 1) map_array[global] = local; - - else { - // search for key - // if found it, just overwrite local value with index - - int previous = -1; - int ibucket = global % map_nbucket; - int index = map_bucket[ibucket]; - while (index > -1) { - if (map_hash[index].global == global) break; - previous = index; - index = map_hash[index].next; - } - if (index > -1) { - map_hash[index].local = local; - return; - } - - // take one entry from free list - // add the new global/local pair as entry at end of bucket list - // special logic if this entry is 1st in bucket - - index = map_free; - map_free = map_hash[map_free].next; - if (previous == -1) map_bucket[ibucket] = index; - else map_hash[previous].next = index; - map_hash[index].global = global; - map_hash[index].local = local; - map_hash[index].next = -1; - map_nused++; - } -} - -/* ---------------------------------------------------------------------- - free the array or hash table for global to local mapping -------------------------------------------------------------------------- */ - -void Atom::map_delete() -{ - if (map_style == 1) { - if (map_tag_max) memory->destroy(map_array); - } else { - if (map_nhash) { - delete [] map_bucket; - delete [] map_hash; - } - map_nhash = 0; - } - map_tag_max = 0; -} - -/* ---------------------------------------------------------------------- - lookup global ID in hash table, return local index -------------------------------------------------------------------------- */ - -int Atom::map_find_hash(int global) -{ - int local = -1; - int index = map_bucket[global % map_nbucket]; - while (index > -1) { - if (map_hash[index].global == global) { - local = map_hash[index].local; - break; - } - index = map_hash[index].next; - } - return local; -} - /* ---------------------------------------------------------------------- add unique tags to any atoms with tag = 0 new tags are grouped by proc and start after max current tag @@ -1644,6 +1402,7 @@ bigint Atom::memory_usage() bigint bytes = avec->memory_usage(); memory->destroy(memstr); + bytes += smax*sizeof(int); if (map_style == 1) bytes += memory->usage(map_array,map_tag_max+1); else if (map_style == 2) { diff --git a/src/atom.h b/src/atom.h index 9a93e620d1..83779da061 100644 --- a/src/atom.h +++ b/src/atom.h @@ -118,6 +118,10 @@ class Atom : protected Pointers { int sortfreq; // sort atoms every this many steps, 0 = off bigint nextsort; // next timestep to sort on + // indices of atoms with same ID + + int *sametag; // sametag[I] = next atom with same ID, -1 if no more + // functions Atom(class LAMMPS *); @@ -190,8 +194,9 @@ class Atom : protected Pointers { // global to local ID mapping - int map_tag_max; - int *map_array; + int map_tag_max; // size of map_array + int *map_array; // direct map of length max atom ID + 1 + int smax; // max size of sametag struct HashElem { int global; // key to search on = global ID @@ -204,8 +209,6 @@ class Atom : protected Pointers { int map_nbucket; // # of hash buckets int *map_bucket; // ptr to 1st entry in each bucket HashElem *map_hash; // hash table - int *primes; // table of prime #s for hashing - int nprimes; // # of primes // spatial sorting of atoms @@ -224,6 +227,7 @@ class Atom : protected Pointers { char *memstr; // string of array names already counted void setup_sort_bins(); + int next_prime(int); }; } diff --git a/src/atom_map.cpp b/src/atom_map.cpp new file mode 100644 index 0000000000..aede10c535 --- /dev/null +++ b/src/atom_map.cpp @@ -0,0 +1,299 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "atom.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define EXTRA 1000 + +/* ---------------------------------------------------------------------- + allocate and initialize array or hash table for global -> local map + set map_tag_max = largest atom ID (may be larger than natoms) + for array option: + array length = 1 to largest tag of any atom + set entire array to -1 as initial values + for hash option: + map_nhash = length of hash table + map_nbucket = # of hash buckets, prime larger than map_nhash + so buckets will only be filled with 0 or 1 atoms on average +------------------------------------------------------------------------- */ + +void Atom::map_init() +{ + map_delete(); + + if (tag_enable == 0) + error->all(FLERR,"Cannot create an atom map unless atoms have IDs"); + + int max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); + MPI_Allreduce(&max,&map_tag_max,1,MPI_INT,MPI_MAX,world); + + smax = nlocal + nghost + EXTRA; + memory->create(sametag,smax,"atom:sametag"); + + if (map_style == 1) { + memory->create(map_array,map_tag_max+1,"atom:map_array"); + for (int i = 0; i <= map_tag_max; i++) map_array[i] = -1; + + } else { + + // map_nhash = max # of atoms that can be hashed on this proc + // set to max of ave atoms/proc or atoms I can store + // multiply by 2, require at least 1000 + // doubling means hash table will be re-init only rarely + + int nper = static_cast (natoms/comm->nprocs); + map_nhash = MAX(nper,nmax); + map_nhash *= 2; + map_nhash = MAX(map_nhash,1000); + + // map_nbucket = prime just larger than map_nhash + + map_nbucket = next_prime(map_nhash); + + // set all buckets to empty + // set hash to map_nhash in length + // put all hash entries in free list and point them to each other + + map_bucket = new int[map_nbucket]; + for (int i = 0; i < map_nbucket; i++) map_bucket[i] = -1; + + map_hash = new HashElem[map_nhash]; + map_nused = 0; + map_free = 0; + for (int i = 0; i < map_nhash; i++) map_hash[i].next = i+1; + map_hash[map_nhash-1].next = -1; + } +} + +/* ---------------------------------------------------------------------- + clear global -> local map for all of my own and ghost atoms + for hash table option: + global ID may not be in table if image atom was already cleared +------------------------------------------------------------------------- */ + +void Atom::map_clear() +{ + if (map_style == 1) { + int nall = nlocal + nghost; + for (int i = 0; i < nall; i++) { + sametag[i] = -1; + map_array[tag[i]] = -1; + } + + } else { + int previous,global,ibucket,index; + int nall = nlocal + nghost; + for (int i = 0; i < nall; i++) { + sametag[i] = -1; + + // search for key + // if don't find it, done + + previous = -1; + global = tag[i]; + ibucket = global % map_nbucket; + index = map_bucket[ibucket]; + while (index > -1) { + if (map_hash[index].global == global) break; + previous = index; + index = map_hash[index].next; + } + if (index == -1) continue; + + // delete the hash entry and add it to free list + // special logic if entry is 1st in the bucket + + if (previous == -1) map_bucket[ibucket] = map_hash[index].next; + else map_hash[previous].next = map_hash[index].next; + + map_hash[index].next = map_free; + map_free = index; + map_nused--; + } + } +} + +/* ---------------------------------------------------------------------- + set global -> local map for all of my own and ghost atoms + loop in reverse order so that nearby images take precedence over far ones + and owned atoms take precedence over images + this enables valid lookups of bond topology atoms + for hash table option: + if hash table too small, re-init + global ID may already be in table if image atom was set +------------------------------------------------------------------------- */ + +void Atom::map_set() +{ + int nall = nlocal + nghost; + if (nall > smax) { + smax = nall + EXTRA; + memory->destroy(sametag); + memory->create(sametag,smax,"atom:sametag"); + } + + if (map_style == 1) { + for (int i = nall-1; i >= 0 ; i--) { + sametag[i] = map_array[tag[i]]; + map_array[tag[i]] = i; + } + + } else { + int previous,global,ibucket,index; + if (nall > map_nhash) map_init(); + + for (int i = nall-1; i >= 0 ; i--) { + sametag[i] = map_find_hash(tag[i]); + + // search for key + // if found it, just overwrite local value with index + + previous = -1; + global = tag[i]; + ibucket = global % map_nbucket; + index = map_bucket[ibucket]; + while (index > -1) { + if (map_hash[index].global == global) break; + previous = index; + index = map_hash[index].next; + } + if (index > -1) { + map_hash[index].local = i; + continue; + } + + // take one entry from free list + // add the new global/local pair as entry at end of bucket list + // special logic if this entry is 1st in bucket + + index = map_free; + map_free = map_hash[map_free].next; + if (previous == -1) map_bucket[ibucket] = index; + else map_hash[previous].next = index; + map_hash[index].global = global; + map_hash[index].local = i; + map_hash[index].next = -1; + map_nused++; + } + } +} + +/* ---------------------------------------------------------------------- + set global to local map for one atom + for hash table option: + global ID may already be in table if atom was already set + called by Special class +------------------------------------------------------------------------- */ + +void Atom::map_one(int global, int local) +{ + if (map_style == 1) map_array[global] = local; + else { + // search for key + // if found it, just overwrite local value with index + + int previous = -1; + int ibucket = global % map_nbucket; + int index = map_bucket[ibucket]; + while (index > -1) { + if (map_hash[index].global == global) break; + previous = index; + index = map_hash[index].next; + } + if (index > -1) { + map_hash[index].local = local; + return; + } + + // take one entry from free list + // add the new global/local pair as entry at end of bucket list + // special logic if this entry is 1st in bucket + + index = map_free; + map_free = map_hash[map_free].next; + if (previous == -1) map_bucket[ibucket] = index; + else map_hash[previous].next = index; + map_hash[index].global = global; + map_hash[index].local = local; + map_hash[index].next = -1; + map_nused++; + } +} + +/* ---------------------------------------------------------------------- + free the array or hash table for global to local mapping +------------------------------------------------------------------------- */ + +void Atom::map_delete() +{ + memory->destroy(sametag); + + if (map_style == 1) { + if (map_tag_max) memory->destroy(map_array); + } else { + if (map_nhash) { + delete [] map_bucket; + delete [] map_hash; + } + map_nhash = 0; + } + map_tag_max = 0; +} + +/* ---------------------------------------------------------------------- + lookup global ID in hash table, return local index + called by map() in atom.h +------------------------------------------------------------------------- */ + +int Atom::map_find_hash(int global) +{ + int local = -1; + int index = map_bucket[global % map_nbucket]; + while (index > -1) { + if (map_hash[index].global == global) { + local = map_hash[index].local; + break; + } + index = map_hash[index].next; + } + return local; +} + +/* ---------------------------------------------------------------------- + return next prime larger than n +------------------------------------------------------------------------- */ + +int Atom::next_prime(int n) +{ + int factor; + + int nprime = n+1; + if (nprime % 2 == 0) nprime++; + int root = static_cast (sqrt(1.0*n)) + 2; + + while (nprime <= MAXSMALLINT) { + for (factor = 3; factor < root; factor++) + if (nprime % factor == 0) break; + if (factor == root) return nprime; + nprime += 2; + } + + return MAXSMALLINT; +} diff --git a/src/domain.cpp b/src/domain.cpp index c5a4b72d2b..140f08053e 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -580,19 +580,6 @@ void Domain::box_too_small_check() "Bond/angle/dihedral extent > half of periodic box length"); } -/* ---------------------------------------------------------------------- - minimum image convention check - return 1 if any distance > 1/2 of box size -------------------------------------------------------------------------- */ - -int Domain::minimum_image_check(double dx, double dy, double dz) -{ - if (xperiodic && fabs(dx) > xprd_half) return 1; - if (yperiodic && fabs(dy) > yprd_half) return 1; - if (zperiodic && fabs(dz) > zprd_half) return 1; - return 0; -} - /* ---------------------------------------------------------------------- minimum image convention use 1/2 of box size as test @@ -717,6 +704,40 @@ void Domain::minimum_image(double *delta) } } +/* ---------------------------------------------------------------------- + return local index of atom J or any of its images that is closest to atom I + if J is not a valid index like -1, just return it +------------------------------------------------------------------------- */ + +int Domain::closest_image(int i, int j) +{ + if (j < 0) return j; + + int *sametag = atom->sametag; + double **x = atom->x; + double *xi = x[i]; + + int closest = j; + double delx = xi[0] - x[j][0]; + double dely = xi[1] - x[j][1]; + double delz = xi[2] - x[j][2]; + double rsqmin = delx*delx + dely*dely + delz*delz; + double rsq; + + while (sametag[j] >= 0) { + j = sametag[j]; + delx = xi[0] - x[j][0]; + dely = xi[1] - x[j][1]; + delz = xi[2] - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < rsqmin) { + rsqmin = rsq; + closest = j; + } + } + return closest; +} + /* ---------------------------------------------------------------------- find and return Xj image = periodic image of Xj that is closest to Xi for triclinic, add/subtract tilt factors in other dims as needed diff --git a/src/domain.h b/src/domain.h index 517fc565b6..13cd6a44c2 100644 --- a/src/domain.h +++ b/src/domain.h @@ -14,6 +14,7 @@ #ifndef LMP_DOMAIN_H #define LMP_DOMAIN_H +#include "math.h" #include "pointers.h" namespace LAMMPS_NS { @@ -94,9 +95,9 @@ class Domain : protected Pointers { virtual void reset_box(); virtual void pbc(); void box_too_small_check(); - int minimum_image_check(double, double, double); void minimum_image(double &, double &, double &); void minimum_image(double *); + int closest_image(int, int); void closest_image(const double * const, const double * const, double * const); void remap(double *, int &); @@ -121,6 +122,17 @@ class Domain : protected Pointers { void bbox(double *, double *, double *, double *); void box_corners(); + // minimum image convention check + // return 1 if any distance > 1/2 of box size + // inline since called from neighbor build inner loop + + inline int minimum_image_check(double dx, double dy, double dz) { + if (xperiodic && fabs(dx) > xprd_half) return 1; + if (yperiodic && fabs(dy) > yprd_half) return 1; + if (zperiodic && fabs(dz) > zprd_half) return 1; + return 0; + } + private: double small[3]; // fractions of box lengths }; diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index a4ec615b17..554b96c989 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -126,7 +126,8 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : } else if (mode == 'm') { double massone = atof(arg[next]); if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake"); - if (nmass == atom->ntypes) error->all(FLERR,"Too many masses for fix shake"); + if (nmass == atom->ntypes) + error->all(FLERR,"Too many masses for fix shake"); mass_list[nmass++] = massone; } else error->all(FLERR,"Illegal fix shake command"); diff --git a/src/neigh_bond.cpp b/src/neigh_bond.cpp index df95e272cb..327d301797 100644 --- a/src/neigh_bond.cpp +++ b/src/neigh_bond.cpp @@ -15,6 +15,7 @@ #include "atom.h" #include "force.h" #include "update.h" +#include "domain.h" #include "memory.h" #include "error.h" @@ -53,7 +54,7 @@ void Neighbor::bond_all() memory->grow(bondlist,maxbond,3,"neighbor:bondlist"); } bondlist[nbondlist][0] = i; - bondlist[nbondlist][1] = atom1; + bondlist[nbondlist][1] = domain->closest_image(i,atom1); bondlist[nbondlist][2] = bond_type[i][m]; nbondlist++; } @@ -92,7 +93,7 @@ void Neighbor::bond_partial() memory->grow(bondlist,maxbond,3,"neighbor:bondlist"); } bondlist[nbondlist][0] = i; - bondlist[nbondlist][1] = atom1; + bondlist[nbondlist][1] = domain->closest_image(i,atom1); bondlist[nbondlist][2] = bond_type[i][m]; nbondlist++; } @@ -134,9 +135,9 @@ void Neighbor::angle_all() maxangle += BONDDELTA; memory->grow(anglelist,maxangle,4,"neighbor:anglelist"); } - anglelist[nanglelist][0] = atom1; - anglelist[nanglelist][1] = atom2; - anglelist[nanglelist][2] = atom3; + anglelist[nanglelist][0] = domain->closest_image(i,atom1); + anglelist[nanglelist][1] = domain->closest_image(i,atom2); + anglelist[nanglelist][2] = domain->closest_image(i,atom3); anglelist[nanglelist][3] = angle_type[i][m]; nanglelist++; } @@ -179,9 +180,9 @@ void Neighbor::angle_partial() maxangle += BONDDELTA; memory->grow(anglelist,maxangle,4,"neighbor:anglelist"); } - anglelist[nanglelist][0] = atom1; - anglelist[nanglelist][1] = atom2; - anglelist[nanglelist][2] = atom3; + anglelist[nanglelist][0] = domain->closest_image(i,atom1); + anglelist[nanglelist][1] = domain->closest_image(i,atom2); + anglelist[nanglelist][2] = domain->closest_image(i,atom3); anglelist[nanglelist][3] = angle_type[i][m]; nanglelist++; } @@ -227,10 +228,10 @@ void Neighbor::dihedral_all() maxdihedral += BONDDELTA; memory->grow(dihedrallist,maxdihedral,5,"neighbor:dihedrallist"); } - dihedrallist[ndihedrallist][0] = atom1; - dihedrallist[ndihedrallist][1] = atom2; - dihedrallist[ndihedrallist][2] = atom3; - dihedrallist[ndihedrallist][3] = atom4; + dihedrallist[ndihedrallist][0] = domain->closest_image(i,atom1); + dihedrallist[ndihedrallist][1] = domain->closest_image(i,atom2); + dihedrallist[ndihedrallist][2] = domain->closest_image(i,atom3); + dihedrallist[ndihedrallist][3] = domain->closest_image(i,atom4); dihedrallist[ndihedrallist][4] = dihedral_type[i][m]; ndihedrallist++; } @@ -277,10 +278,10 @@ void Neighbor::dihedral_partial() maxdihedral += BONDDELTA; memory->grow(dihedrallist,maxdihedral,5,"neighbor:dihedrallist"); } - dihedrallist[ndihedrallist][0] = atom1; - dihedrallist[ndihedrallist][1] = atom2; - dihedrallist[ndihedrallist][2] = atom3; - dihedrallist[ndihedrallist][3] = atom4; + dihedrallist[ndihedrallist][0] = domain->closest_image(i,atom1); + dihedrallist[ndihedrallist][1] = domain->closest_image(i,atom2); + dihedrallist[ndihedrallist][2] = domain->closest_image(i,atom3); + dihedrallist[ndihedrallist][3] = domain->closest_image(i,atom4); dihedrallist[ndihedrallist][4] = dihedral_type[i][m]; ndihedrallist++; } @@ -326,10 +327,10 @@ void Neighbor::improper_all() maximproper += BONDDELTA; memory->grow(improperlist,maximproper,5,"neighbor:improperlist"); } - improperlist[nimproperlist][0] = atom1; - improperlist[nimproperlist][1] = atom2; - improperlist[nimproperlist][2] = atom3; - improperlist[nimproperlist][3] = atom4; + improperlist[nimproperlist][0] = domain->closest_image(i,atom1); + improperlist[nimproperlist][1] = domain->closest_image(i,atom2); + improperlist[nimproperlist][2] = domain->closest_image(i,atom3); + improperlist[nimproperlist][3] = domain->closest_image(i,atom4); improperlist[nimproperlist][4] = improper_type[i][m]; nimproperlist++; } @@ -376,10 +377,10 @@ void Neighbor::improper_partial() maximproper += BONDDELTA; memory->grow(improperlist,maximproper,5,"neighbor:improperlist"); } - improperlist[nimproperlist][0] = atom1; - improperlist[nimproperlist][1] = atom2; - improperlist[nimproperlist][2] = atom3; - improperlist[nimproperlist][3] = atom4; + improperlist[nimproperlist][0] = domain->closest_image(i,atom1); + improperlist[nimproperlist][1] = domain->closest_image(i,atom2); + improperlist[nimproperlist][2] = domain->closest_image(i,atom3); + improperlist[nimproperlist][3] = domain->closest_image(i,atom4); improperlist[nimproperlist][4] = improper_type[i][m]; nimproperlist++; }