ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8391 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2012-06-25 16:42:20 +00:00
parent bee842bed7
commit e8ff810ee7
78 changed files with 396 additions and 474 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/CLASS2/angle_class2.cpp
View File

@ -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);

1
src/CLASS2/bond_class2.cpp
View File

@ -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);

4
src/CLASS2/dihedral_class2.cpp
View File

@ -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

6
src/CLASS2/improper_class2.cpp
View File

@ -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

3
src/MOLECULE/angle_charmm.cpp
View File

@ -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);

2
src/MOLECULE/angle_cosine.cpp
View File

@ -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);

2
src/MOLECULE/angle_cosine_delta.cpp
View File

@ -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);

2
src/MOLECULE/angle_cosine_periodic.cpp
View File

@ -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);

2
src/MOLECULE/angle_cosine_squared.cpp
View File

@ -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);

2
src/MOLECULE/angle_harmonic.cpp
View File

@ -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);

2
src/MOLECULE/angle_table.cpp
View File

@ -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);

1
src/MOLECULE/bond_fene.cpp
View File

@ -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

1
src/MOLECULE/bond_fene_expand.cpp
View File

@ -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

1
src/MOLECULE/bond_harmonic.cpp
View File

@ -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);

1
src/MOLECULE/bond_morse.cpp
View File

@ -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);

1
src/MOLECULE/bond_nonlinear.cpp
View File

@ -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);

1
src/MOLECULE/bond_quartic.cpp
View File

@ -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;

1
src/MOLECULE/bond_table.cpp
View File

@ -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);

5
src/MOLECULE/dihedral_charmm.cpp
View File

@ -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;

4
src/MOLECULE/dihedral_harmonic.cpp
View File

@ -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

4
src/MOLECULE/dihedral_helix.cpp
View File

@ -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

4
src/MOLECULE/dihedral_multi_harmonic.cpp
View File

@ -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

4
src/MOLECULE/dihedral_opls.cpp
View File

@ -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

4
src/MOLECULE/improper_cvff.cpp
View File

@ -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

3
src/MOLECULE/improper_harmonic.cpp
View File

@ -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);

3
src/MOLECULE/improper_umbrella.cpp
View File

@ -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

3
src/USER-CG-CMM/angle_cg_cmm.cpp
View File

@ -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);

3
src/USER-CG-CMM/angle_sdk.cpp
View File

@ -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];

2
src/USER-MISC/angle_cosine_shift.cpp
View File

@ -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);

2
src/USER-MISC/angle_cosine_shift_exp.cpp
View File

@ -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);

1
src/USER-MISC/angle_dipole.cpp
View File

@ -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);

4
src/USER-MISC/bond_harmonic_shift.cpp
View File

@ -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

4
src/USER-MISC/bond_harmonic_shift_cut.cpp
View File

@ -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

4
src/USER-MISC/dihedral_cosine_shift_exp.cpp
View File

@ -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

3
src/USER-MISC/improper_cossq.cpp
View File

@ -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);

5
src/USER-MISC/improper_ring.cpp
View File

@ -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]

3
src/USER-OMP/angle_charmm_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_class2_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_delta_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_periodic_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_shift_exp_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_shift_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_cosine_squared_omp.cpp
View File

@ -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);

1
src/USER-OMP/angle_dipole_omp.cpp
View File

@ -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);

2
src/USER-OMP/angle_harmonic_omp.cpp
View File

@ -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);

3
src/USER-OMP/angle_sdk_omp.cpp
View File

@ -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];

2
src/USER-OMP/angle_table_omp.cpp
View File

@ -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);

1
src/USER-OMP/bond_class2_omp.cpp
View File

@ -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);

1
src/USER-OMP/bond_fene_expand_omp.cpp
View File

@ -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);

1
src/USER-OMP/bond_fene_omp.cpp
View File

@ -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];

1
src/USER-OMP/bond_harmonic_omp.cpp
View File

@ -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);

4
src/USER-OMP/bond_harmonic_shift_cut_omp.cpp
View File

@ -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

4
src/USER-OMP/bond_harmonic_shift_omp.cpp
View File

@ -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

1
src/USER-OMP/bond_morse_omp.cpp
View File

@ -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);

1
src/USER-OMP/bond_nonlinear_omp.cpp
View File

@ -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);

1
src/USER-OMP/bond_quartic_omp.cpp
View File

@ -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;

1
src/USER-OMP/bond_table_omp.cpp
View File

@ -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);

5
src/USER-OMP/dihedral_charmm_omp.cpp
View File

@ -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;

4
src/USER-OMP/dihedral_class2_omp.cpp
View File

@ -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

4
src/USER-OMP/dihedral_cosine_shift_exp_omp.cpp
View File

@ -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

4
src/USER-OMP/dihedral_harmonic_omp.cpp
View File

@ -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

4
src/USER-OMP/dihedral_helix_omp.cpp
View File

@ -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

4
src/USER-OMP/dihedral_multi_harmonic_omp.cpp
View File

@ -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

4
src/USER-OMP/dihedral_opls_omp.cpp
View File

@ -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

6
src/USER-OMP/improper_class2_omp.cpp
View File

@ -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

3
src/USER-OMP/improper_cossq_omp.cpp
View File

@ -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);

4
src/USER-OMP/improper_cvff_omp.cpp
View File

@ -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

3
src/USER-OMP/improper_harmonic_omp.cpp
View File

@ -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);

5
src/USER-OMP/improper_ring_omp.cpp
View File

@ -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]

3
src/USER-OMP/improper_umbrella_omp.cpp
View File

@ -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

253
src/atom.cpp
View File

@ -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<int> (natoms/comm->nprocs);
map_nhash = MAX(nper,nmax);
if (map_nhash > natoms) map_nhash = static_cast<int> (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) {

12
src/atom.h
View File

@ -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);
};
}

299
src/atom_map.cpp Normal file
View File

@ -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<int> (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<int> (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;
}

47
src/domain.cpp
View File

@ -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

14
src/domain.h
View File

@ -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
};

3
src/fix_shake.cpp
View File

@ -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");

49
src/neigh_bond.cpp
View File

@ -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++;
}