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

Merge remote-tracking branch 'devel/whitespace-cleanup' into collected-small-fixes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2018-12-03 15:06:23 -05:00
parent 4af7e706d7 1c795a6be3
commit e1aff27cd4
69 changed files with 390 additions and 390 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/BODY/body_rounded_polygon.cpp
View File

@ -170,7 +170,7 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
// nentries = number of double entries to be read from Body section:
// 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
int nentries = 6 + 3*nsub + 1;
int nentries = 6 + 3*nsub + 1;
if (ndouble != nentries)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
@ -279,7 +279,7 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 1;
k += 2;
}
}
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
@ -341,7 +341,7 @@ double BodyRoundedPolygon::radius_body(int /*ninteger*/, int ndouble,
double maxrad = 0.0;
double delta[3];
int offset = 6;
int offset = 6;
for (int i = 0; i < nsub; i++) {
delta[0] = dfile[offset];
delta[1] = dfile[offset+1];
@ -350,9 +350,9 @@ double BodyRoundedPolygon::radius_body(int /*ninteger*/, int ndouble,
onerad = MathExtra::len3(delta);
maxrad = MAX(maxrad,onerad);
}
// add in radius of rounded corners
return maxrad + 0.5*dfile[offset];
}
@ -401,7 +401,7 @@ int BodyRoundedPolygon::image(int ibonus, double flag1, double /*flag2*/,
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int n = bonus->ivalue[0];
if (n == 1) {
for (int i = 0; i < n; i++) {
imflag[i] = SPHERE;
@ -418,7 +418,7 @@ int BodyRoundedPolygon::image(int ibonus, double flag1, double /*flag2*/,
}
} else {
// first end pt of each line
for (int i = 0; i < n; i++) {
@ -436,7 +436,7 @@ int BodyRoundedPolygon::image(int ibonus, double flag1, double /*flag2*/,
}
// second end pt of each line
for (int i = 0; i < n; i++) {
j = i+1;
if (j == n) j = 0;

12
src/BODY/body_rounded_polyhedron.cpp
View File

@ -132,7 +132,7 @@ double BodyRoundedPolyhedron::enclosing_radius(struct AtomVecBody::Bonus *bonus)
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2);
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
MAX_FACE_SIZE*nfaces(bonus));
}
@ -143,7 +143,7 @@ double BodyRoundedPolyhedron::rounded_radius(struct AtomVecBody::Bonus *bonus)
int nvertices = bonus->ivalue[0];
if (nvertices == 1 || nvertices == 2)
return *(bonus->dvalue+3*nsub(bonus)+2+1);
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
return *(bonus->dvalue+3*nsub(bonus) + 2*nedges(bonus) +
MAX_FACE_SIZE*nfaces(bonus)+1);
}
@ -207,7 +207,7 @@ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
// nsub == 1 || nsub == 2 || nsub == 3:
// 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
// nsub > 3:
// 6 for inertia + 3*nsub for vertex coords + 2*nsub for edges +
// 6 for inertia + 3*nsub for vertex coords + 2*nsub for edges +
// 3*nfaces + 1 for rounded radius
int nedges,nentries;
@ -327,7 +327,7 @@ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
bonus->dvalue[k] = 0;
*(&bonus->dvalue[k]+1) = 1;
k += 2;
}
}
erad = sqrt(erad2);
bonus->dvalue[k] = erad;
@ -407,7 +407,7 @@ double BodyRoundedPolyhedron::radius_body(int /*ninteger*/, int ndouble,
double maxrad = 0.0;
double delta[3];
int offset = 6;
int offset = 6;
for (int i = 0; i < nsub; i++) {
delta[0] = dfile[offset];
delta[1] = dfile[offset+1];
@ -420,7 +420,7 @@ double BodyRoundedPolyhedron::radius_body(int /*ninteger*/, int ndouble,
if (nsub > 2) offset += (2*nedges+MAX_FACE_SIZE*nfac);
// add in radius of rounded corners
return maxrad + 0.5*dfile[offset];
}

4
src/BODY/fix_wall_body_polygon.cpp
View File

@ -263,7 +263,7 @@ void FixWallBodyPolygon::post_force(int /*vflag*/)
}
ndiscrete = nedge = 0;
for (i = 0; i < nlocal; i++)
for (i = 0; i < nlocal; i++)
dnum[i] = ednum[i] = 0;
for (i = 0; i < nlocal; i++) {
@ -675,7 +675,7 @@ void FixWallBodyPolygon::contact_forces(Contact& contact, double j_a,
AtomVecBody::Bonus *bonus;
ibody = contact.ibody;
// compute the velocity of the vertex in the space-fixed frame
ibonus = atom->body[ibody];

18
src/BODY/fix_wall_body_polyhedron.cpp
View File

@ -102,7 +102,7 @@ FixWallBodyPolyhedron::FixWallBodyPolyhedron(LAMMPS *lmp, int narg, char **arg)
if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
else hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
}
}
// check for trailing keyword/values
@ -274,7 +274,7 @@ void FixWallBodyPolyhedron::post_force(int /*vflag*/)
}
ndiscrete = nedge = nface = 0;
for (i = 0; i < nlocal; i++)
for (i = 0; i < nlocal; i++)
dnum[i] = ednum[i] = facnum[i] = 0;
for (i = 0; i < nlocal; i++) {
@ -320,7 +320,7 @@ void FixWallBodyPolyhedron::post_force(int /*vflag*/)
wall_pos = whi;
side = ZHI;
}
}
}
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) continue;
@ -502,7 +502,7 @@ int FixWallBodyPolyhedron::sphere_against_wall(int i, double wall_pos,
hi[0] = x[i][0];
hi[1] = x[i][1];
hi[2] = wall_pos;
}
}
distance(hi, x[i], d);
@ -536,7 +536,7 @@ int FixWallBodyPolyhedron::sphere_against_wall(int i, double wall_pos,
Output:
contact_list = list of contacts between i and the wall
num_contacts = number of contacts between i's vertices and the wall
Return:
Return:
number of contacts of the edge to the wall (0, 1 or 2)
---------------------------------------------------------------------- */
@ -578,7 +578,7 @@ int FixWallBodyPolyhedron::edge_against_wall(int i, double wall_pos,
------------------------------------------------------------------------- */
int FixWallBodyPolyhedron::compute_distance_to_wall(int ibody, int edge_index,
double *xmi, double rounded_radius_i, double wall_pos,
double *xmi, double rounded_radius_i, double wall_pos,
int /*side*/, double* vwall, int &contact)
{
int mode,ifirst,iefirst,npi1,npi2;
@ -623,7 +623,7 @@ int FixWallBodyPolyhedron::compute_distance_to_wall(int ibody, int edge_index,
hi[0] = xpi1[0];
hi[1] = xpi1[1];
hi[2] = wall_pos;
}
}
distance(hi, xpi1, d1);
@ -658,7 +658,7 @@ int FixWallBodyPolyhedron::compute_distance_to_wall(int ibody, int edge_index,
hi[0] = xpi2[0];
hi[1] = xpi2[1];
hi[2] = wall_pos;
}
}
distance(hi, xpi2, d2);
@ -790,7 +790,7 @@ void FixWallBodyPolyhedron::contact_forces(Contact& contact, double j_a,
AtomVecBody::Bonus *bonus;
ibody = contact.ibody;
// compute the velocity of the vertex in the space-fixed frame
ibonus = atom->body[ibody];

2
src/BODY/fix_wall_body_polyhedron.h
View File

@ -90,7 +90,7 @@ class FixWallBodyPolyhedron : public Fix {
double rounded_radius_i, double wall_pos, int side,
double* vwall, int &contact);
double contact_separation(const Contact& c1, const Contact& c2);
void contact_forces(int ibody, double j_a, double *xi, double *xj,
void contact_forces(int ibody, double j_a, double *xi, double *xj,
double delx, double dely, double delz,
double fx, double fy, double fz, double** x, double** v,
double** angmom, double** f, double** torque, double* vwall);

22
src/BODY/pair_body_rounded_polygon.cpp
View File

@ -275,9 +275,9 @@ void PairBodyRoundedPolygon::compute(int eflag, int vflag)
// scale the force at both contacts
contact_forces(contact_list[m], j_a, x, v, angmom, f, torque,
contact_forces(contact_list[m], j_a, x, v, angmom, f, torque,
evdwl, facc);
contact_forces(contact_list[n], j_a, x, v, angmom, f, torque,
contact_forces(contact_list[n], j_a, x, v, angmom, f, torque,
evdwl, facc);
done = 1;
@ -288,15 +288,15 @@ void PairBodyRoundedPolygon::compute(int eflag, int vflag)
"xv = %f %f %f; xe = %f %f %f\n",
m, contact_list[m].vertex, contact_list[m].ibody,
contact_list[m].edge, contact_list[m].jbody,
contact_list[m].xv[0], contact_list[m].xv[1],
contact_list[m].xv[2], contact_list[m].xe[0],
contact_list[m].xv[0], contact_list[m].xv[1],
contact_list[m].xv[2], contact_list[m].xe[0],
contact_list[m].xe[1], contact_list[m].xe[2]);
printf(" %d: vertex %d of body %d and edge %d of body %d; "
"xv = %f %f %f; xe = %f %f %f\n",
n, contact_list[n].vertex, contact_list[n].ibody,
contact_list[n].edge, contact_list[n].jbody,
contact_list[n].xv[0], contact_list[n].xv[1],
contact_list[n].xv[2], contact_list[n].xe[0],
contact_list[n].xv[0], contact_list[n].xv[1],
contact_list[n].xv[2], contact_list[n].xe[0],
contact_list[n].xe[1], contact_list[n].xe[2]);
#endif
@ -420,7 +420,7 @@ void PairBodyRoundedPolygon::coeff(int narg, char **arg)
void PairBodyRoundedPolygon::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
if (!avec)
if (!avec)
error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
error->all(FLERR,"Pair body/rounded/polygon requires "
@ -501,7 +501,7 @@ void PairBodyRoundedPolygon::init_style()
if (dnum[i] == 0) body2space(i);
eradi = enclosing_radius[i];
if (eradi > merad[itype]) merad[itype] = eradi;
} else
} else
merad[itype] = 0;
}
@ -669,7 +669,7 @@ void PairBodyRoundedPolygon::sphere_against_sphere(int i, int j,
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
@ -860,7 +860,7 @@ int PairBodyRoundedPolygon::vertex_against_edge(int i, int j,
#ifdef _CONVEX_POLYGON
// done with the edges from body j,
// given that vertex ni interacts with only one vertex
// given that vertex ni interacts with only one vertex
// from one edge of body j
break;
#endif
@ -1124,7 +1124,7 @@ int PairBodyRoundedPolygon::compute_distance_to_vertex(int ibody,
mode = EDGE;
if (d < contact_dist + EPSILON)
contact = 1;
} else { // t < 0 || t > 1: closer to either vertices of the edge
if (t < 0) {

72
src/BODY/pair_body_rounded_polyhedron.cpp
View File

@ -283,7 +283,7 @@ void PairBodyRoundedPolyhedron::compute(int eflag, int vflag)
// check interaction between i's edges and j' faces
#ifdef _POLYHEDRON_DEBUG
printf("INTERACTION between edges of %d vs. faces of %d:\n", i, j);
#endif
#endif
interact = edge_against_face(i, j, itype, jtype, x, contact_list,
num_contacts, evdwl, facc);
@ -297,7 +297,7 @@ void PairBodyRoundedPolyhedron::compute(int eflag, int vflag)
// check interaction between i's edges and j' edges
#ifdef _POLYHEDRON_DEBUG
printf("INTERACTION between edges of %d vs. edges of %d:\n", i, j);
#endif
#endif
interact = edge_against_edge(i, j, itype, jtype, x, contact_list,
num_contacts, evdwl, facc);
@ -478,7 +478,7 @@ void PairBodyRoundedPolyhedron::init_style()
if (dnum[i] == 0) body2space(i);
eradi = enclosing_radius[i];
if (eradi > merad[itype]) merad[itype] = eradi;
} else
} else
merad[itype] = 0;
}
@ -630,7 +630,7 @@ void PairBodyRoundedPolyhedron::sphere_against_sphere(int ibody, int jbody,
vr3 = v[ibody][2] - v[jbody][2];
// normal component
rsqinv = 1.0/rsq;
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr * rsqinv;
@ -664,7 +664,7 @@ void PairBodyRoundedPolyhedron::sphere_against_sphere(int ibody, int jbody,
f[ibody][0] += fx;
f[ibody][1] += fy;
f[ibody][2] += fz;
if (newton_pair || jbody < nlocal) {
f[jbody][0] -= fx;
f[jbody][1] -= fy;
@ -797,7 +797,7 @@ void PairBodyRoundedPolyhedron::sphere_against_edge(int ibody, int jbody,
fn[1] = -c_n * vn2;
fn[2] = -c_n * vn3;
// tangential friction term at contact,
// tangential friction term at contact,
// excluding the tangential deformation term
ft[0] = -c_t * vt1;
@ -873,7 +873,7 @@ void PairBodyRoundedPolyhedron::sphere_against_face(int ibody, int jbody,
xi3[2] = x[ibody][2] + discrete[ifirst+npi3][2];
// find the normal unit vector of the face
MathExtra::sub3(xi2, xi1, ui);
MathExtra::sub3(xi3, xi1, vi);
MathExtra::cross3(ui, vi, n);
@ -1076,7 +1076,7 @@ int PairBodyRoundedPolyhedron::edge_against_face(int ibody, int jbody,
itype, jtype, cut_inner,
contact_list, num_contacts,
energy, facc);
}
}
} // end for looping through the edges of body i
@ -1171,7 +1171,7 @@ int PairBodyRoundedPolyhedron::interaction_edge_to_edge(int ibody,
printf(" edge npi1 = %d (%f %f %f); npi2 = %d (%f %f %f) vs."
" edge npj1 = %d (%f %f %f); npj2 = %d (%f %f %f): "
"t1 = %f; t2 = %f; r = %f; dot = %f\n",
npi1, xi1[0], xi1[1], xi1[2], npi2, xi2[0], xi2[1], xi2[2],
npi1, xi1[0], xi1[1], xi1[2], npi2, xi2[0], xi2[1], xi2[2],
npj1, xpj1[0], xpj1[1], xpj1[2], npj2, xpj2[0], xpj2[1], xpj2[2],
t1, t2, r, dot);
#endif
@ -1276,7 +1276,7 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
xi3[2] = xmi[2] + discrete[ifirst+npi3][2];
// find the normal unit vector of the face, ensure it point outward of the body
MathExtra::sub3(xi2, xi1, ui);
MathExtra::sub3(xi3, xi1, vi);
MathExtra::cross3(ui, vi, n);
@ -1305,7 +1305,7 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
xpj2[1] = xmj[1] + discrete[jfirst+npj2][1];
xpj2[2] = xmj[2] + discrete[jfirst+npj2][2];
// no interaction if two ends of the edge
// no interaction if two ends of the edge
// are on the same side with the COM wrt the face
if (opposite_sides(n, xi1, xmi, xpj1) == 0 &&
@ -1340,9 +1340,9 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
int jflag = 1;
#ifdef _POLYHEDRON_DEBUG
if (interact == EF_SAME_SIDE_OF_FACE)
if (interact == EF_SAME_SIDE_OF_FACE)
printf(" - same side of face\n");
else if (interact == EF_PARALLEL)
else if (interact == EF_PARALLEL)
printf(" - parallel\n");
printf(" face: xi1 (%f %f %f) xi2 (%f %f %f) xi3 (%f %f %f)\n",
xi1[0], xi1[1], xi1[2], xi2[0], xi2[1], xi2[2], xi3[0], xi3[1], xi3[2]);
@ -1387,9 +1387,9 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
} else {
num_outside++;
}
}
}
// xpj2 is in the interaction zone
// xpj2 is in the interaction zone
// and its projection on the face is inside the triangle
// compute vertex-face interaction and accumulate force/torque to both bodies
@ -1401,7 +1401,7 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
jflag, energy, facc);
#ifdef _POLYHEDRON_DEBUG
printf(" - compute pair force between vertex %d from edge %d of body %d "
"with face %d of body %d: d2 = %f\n",
"with face %d of body %d: d2 = %f\n",
npj2, edge_index, jbody, face_index, ibody, d2);
#endif
@ -1448,7 +1448,7 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
// compute interaction between the edge with the three edges of the face
#ifdef _POLYHEDRON_DEBUG
printf(" - intersect outside triangle\n");
printf(" - intersect outside triangle\n");
printf(" - compute pair force between edge %d of body %d "
"with face %d of body %d\n", edge_index, jbody, face_index, ibody);
printf(" face: xi1 (%f %f %f) xi2 (%f %f %f) xi3 (%f %f %f)\n",
@ -1693,7 +1693,7 @@ void PairBodyRoundedPolyhedron::rescale_cohesive_forces(double** x,
xc[0] /= (double)num_unique_contacts;
xc[1] /= (double)num_unique_contacts;
xc[2] /= (double)num_unique_contacts;
contact_area = 0.0;
for (int m = 0; m < num_contacts; m++) {
if (contact_list[m].unique == 0) continue;
@ -1847,12 +1847,12 @@ int PairBodyRoundedPolyhedron::edge_face_intersect(double* x1, double* x2,
if (t < 0 || t > 1) {
interact = EF_SAME_SIDE_OF_FACE;
} else {
if (inside == 1)
if (inside == 1)
interact = EF_INTERSECT_INSIDE;
else
interact = EF_INTERSECT_OUTSIDE;
}
return interact;
}
@ -1998,7 +1998,7 @@ void PairBodyRoundedPolyhedron::project_pt_plane(const double* q,
q_proj[2] = q[2] + n[2] * t;
// check if the projection point is inside the triangle
// exclude the edges and vertices
// exclude the edges and vertices
// edge-sphere and sphere-sphere interactions are handled separately
inside = 0;
@ -2022,7 +2022,7 @@ void PairBodyRoundedPolyhedron::project_pt_line(const double* q,
MathExtra::sub3(xi2, xi1, u);
MathExtra::norm3(u);
MathExtra::sub3(q, xi1, v);
s = MathExtra::dot3(u, v);
h[0] = xi1[0] + s * u[0];
h[1] = xi1[1] + s * u[1];
@ -2039,13 +2039,13 @@ void PairBodyRoundedPolyhedron::project_pt_line(const double* q,
t = (h[2] - xi1[2])/(xi2[2] - xi1[2]);
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
compute the shortest distance between two edges (line segments)
x1, x2: two endpoints of the first edge
x3, x4: two endpoints of the second edge
h1: the end point of the shortest segment perpendicular to both edges
h1: the end point of the shortest segment perpendicular to both edges
on the line (x1;x2)
h2: the end point of the shortest segment perpendicular to both edges
h2: the end point of the shortest segment perpendicular to both edges
on the line (x3;x4)
t1: fraction of h1 in the segment (x1,x2)
t2: fraction of h2 in the segment (x3,x4)
@ -2080,7 +2080,7 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
double s1,s2,x13[3],x23[3],x13h[3];
double t13,t23,t31,t41,x31[3],x41[3];
t13=t23=t31=t41=0.0;
MathExtra::sub3(x1,x3,x13); // x13 = x1 - x3
MathExtra::sub3(x2,x3,x23); // x23 = x2 - x3
@ -2089,7 +2089,7 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
x13h[1] = x13[1] - s1*v[1];
x13h[2] = x13[2] - s1*v[2];
r = MathExtra::len3(x13h);
// x13 is the projection of x1 on x3-x4
x13[0] = x3[0] + s1*v[0];
@ -2102,7 +2102,7 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
x23[0] = x3[0] + s2*v[0];
x23[1] = x3[1] + s2*v[1];
x23[2] = x3[2] + s2*v[2];
// find the fraction of the projection points on the edges
if (fabs(x4[0] - x3[0]) > 0)
@ -2177,8 +2177,8 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
h1[0] = (x1[0]+x2[0])/2;
h1[1] = (x1[1]+x2[1])/2;
h1[2] = (x1[2]+x2[2])/2;
h2[0] = (x13[0]+x23[0])/2;
h2[1] = (x13[1]+x23[1])/2;
h2[0] = (x13[0]+x23[0])/2;
h2[1] = (x13[1]+x23[1])/2;
h2[2] = (x13[2]+x23[2])/2;
t1 = 0.5;
t2 = (t13+t23)/2;
@ -2216,8 +2216,8 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
h1[0] = (x31[0]+x41[0])/2;
h1[1] = (x31[1]+x41[1])/2;
h1[2] = (x31[2]+x41[2])/2;
h2[0] = (x3[0]+x4[0])/2;
h2[1] = (x3[1]+x4[1])/2;
h2[0] = (x3[0]+x4[0])/2;
h2[1] = (x3[1]+x4[1])/2;
h2[2] = (x3[2]+x4[2])/2;
t1 = (t31+t41)/2;
t2 = 0.5;
@ -2225,7 +2225,7 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
n2++;
}
}
}
}
// if n1 == 0 and n2 == 0 at this point,
// which means no overlapping segments bt two parallel edges,
@ -2233,10 +2233,10 @@ void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
return;
}
}
// find the vector n perpendicular to both edges
MathExtra::cross3(u, v, n);
MathExtra::norm3(n);
@ -2326,7 +2326,7 @@ double PairBodyRoundedPolyhedron::contact_separation(const Contact& c1,
find the number of unique contacts
------------------------------------------------------------------------- */
void PairBodyRoundedPolyhedron::find_unique_contacts(Contact* contact_list,
void PairBodyRoundedPolyhedron::find_unique_contacts(Contact* contact_list,
int& num_contacts)
{
int n = num_contacts;

6
src/BODY/pair_body_rounded_polyhedron.h
View File

@ -155,13 +155,13 @@ class PairBodyRoundedPolyhedron : public Pair {
int& inside_a, int& inside_b);
// helper functions
int opposite_sides(double* n, double* x0, double* a, double* b);
void project_pt_plane(const double* q, const double* p,
void project_pt_plane(const double* q, const double* p,
const double* n, double* q_proj, double &d);
void project_pt_plane(const double* q, const double* x1, const double* x2,
void project_pt_plane(const double* q, const double* x1, const double* x2,
const double* x3, double* q_proj, double &d, int& inside);
void project_pt_line(const double* q, const double* xi1, const double* xi2,
double* h, double& d, double& t);
void inside_polygon(int ibody, int face_index, double* xmi,
void inside_polygon(int ibody, int face_index, double* xmi,
const double* q1, const double* q2, int& inside1, int& inside2);
void distance_bt_edges(const double* x1, const double* x2,

4
src/KOKKOS/atom_vec_sphere_kokkos.cpp
View File

@ -792,7 +792,7 @@ int AtomVecSphereKokkos::pack_comm_self(
const int &n, const DAT::tdual_int_2d &list, const int &iswap,
const int nfirst, const int &pbc_flag, const int* const pbc) {
// Fallback to AtomVecKokkos if radvary == 0
if (radvary == 0)
if (radvary == 0)
return AtomVecKokkos::pack_comm_self(n,list,iswap,nfirst,pbc_flag,pbc);
if(commKK->forward_comm_on_host) {
sync(Host,X_MASK|RADIUS_MASK|RMASS_MASK);
@ -1462,7 +1462,7 @@ struct AtomVecSphereKokkos_PackBorder {
const int maxsend = (buf.extent(0)*buf.extent(1))/elements;
_buf = typename ArrayTypes<DeviceType>::t_xfloat_2d_um(buf.data(),maxsend,elements);
}
KOKKOS_INLINE_FUNCTION
void operator() (const int& i) const {
const int j = _list(_iswap,i);

2
src/KOKKOS/comm_kokkos.cpp
View File

@ -210,7 +210,7 @@ void CommKokkos::forward_comm_device(int dummy)
space,X_MASK);
}
} else if (ghost_velocity) {
if (size_forward_recv[iswap]) {
if (size_forward_recv[iswap]) {
MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
size_forward_recv[iswap],MPI_DOUBLE,
recvproc[iswap],0,world,&request);

2
src/KOKKOS/fix_freeze_kokkos.cpp
View File

@ -75,7 +75,7 @@ void FixFreezeKokkos<DeviceType>::post_force(int vflag)
OriginalForce original;
Kokkos::parallel_reduce(nlocal, *this, original);
copymode = 0;
foriginal[0] = original.values[0];
foriginal[1] = original.values[1];
foriginal[2] = original.values[2];

4
src/KOKKOS/fix_freeze_kokkos.h
View File

@ -26,7 +26,7 @@ FixStyle(freeze/kk/host,FixFreezeKokkos<LMPHostType>)
#include "kokkos_type.h"
namespace LAMMPS_NS {
template<class DeviceType>
class FixFreezeKokkos : public FixFreeze {
public:
@ -66,7 +66,7 @@ class FixFreezeKokkos : public FixFreeze {
KOKKOS_INLINE_FUNCTION
void operator()(const int i, OriginalForce &original) const;
private:
typename ArrayTypes<DeviceType>::t_f_array f;
typename ArrayTypes<DeviceType>::t_f_array torque;

2
src/KOKKOS/fix_gravity_kokkos.h
View File

@ -50,7 +50,7 @@ class FixGravityKokkos : public FixGravity {
typename ArrayTypes<DeviceType>::t_int_1d type;
typename ArrayTypes<DeviceType>::t_int_1d mask;
};
} // namespace LAMMPS_NS
#endif // LMP_FIX_GRAVITY_KOKKOS_H

8
src/KOKKOS/fix_neigh_history_kokkos.cpp
View File

@ -77,7 +77,7 @@ template <class DeviceType>
void FixNeighHistoryKokkos<DeviceType>::pre_exchange()
{
copymode = 1;
h_resize() = 1;
while (h_resize() > 0) {
FixNeighHistoryKokkosZeroPartnerCountFunctor<DeviceType> zero(this);
@ -156,7 +156,7 @@ template <class DeviceType>
void FixNeighHistoryKokkos<DeviceType>::post_neighbor()
{
tag = atomKK->k_tag.view<DeviceType>();
int inum = pair->list->inum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(pair->list);
d_numneigh = k_list->d_numneigh;
@ -179,7 +179,7 @@ void FixNeighHistoryKokkos<DeviceType>::post_neighbor()
}
copymode = 1;
FixNeighHistoryKokkosPostNeighborFunctor<DeviceType> f(this);
Kokkos::parallel_for(inum,f);
@ -251,7 +251,7 @@ void FixNeighHistoryKokkos<DeviceType>::grow_arrays(int nmax)
k_npartner.template sync<LMPHostType>(); // force reallocation on host
k_partner.template sync<LMPHostType>();
k_valuepartner.template sync<LMPHostType>();
memoryKK->grow_kokkos(k_npartner,npartner,nmax,"neighbor_history:npartner");
memoryKK->grow_kokkos(k_partner,partner,nmax,maxpartner,"neighbor_history:partner");
memoryKK->grow_kokkos(k_valuepartner,valuepartner,nmax,dnum*maxpartner,"neighbor_history:valuepartner");

4
src/KOKKOS/fix_neigh_history_kokkos.h
View File

@ -61,11 +61,11 @@ class FixNeighHistoryKokkos : public FixNeighHistory {
typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
typename ArrayTypes<DeviceType>::t_tagint_1d tag;
typename ArrayTypes<DeviceType>::t_int_1d d_npartner;
typename ArrayTypes<DeviceType>::t_tagint_2d d_partner;
typename ArrayTypes<DeviceType>::t_float_2d d_valuepartner;
typename ArrayTypes<DeviceType>::t_float_2d d_valuepartner;
typename ArrayTypes<DeviceType>::t_int_scalar d_resize;
typename ArrayTypes<LMPHostType>::t_int_scalar h_resize;

6
src/KOKKOS/fix_nve_sphere_kokkos.cpp
View File

@ -75,7 +75,7 @@ void FixNVESphereKokkos<DeviceType>::initial_integrate(int vflag)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
FixNVESphereKokkosInitialIntegrateFunctor<DeviceType> f(this);
FixNVESphereKokkosInitialIntegrateFunctor<DeviceType> f(this);
Kokkos::parallel_for(nlocal,f);
}
@ -95,7 +95,7 @@ void FixNVESphereKokkos<DeviceType>::initial_integrate_item(const int i) const
x(i,0) += dtv * v(i,0);
x(i,1) += dtv * v(i,1);
x(i,2) += dtv * v(i,2);
const double dtirotate = dtfrotate / (radius(i)*radius(i)*rmass(i));
omega(i,0) += dtirotate * torque(i,0);
omega(i,1) += dtirotate * torque(i,1);
@ -139,7 +139,7 @@ void FixNVESphereKokkos<DeviceType>::final_integrate_item(const int i) const
v(i,0) += dtfm * f(i,0);
v(i,1) += dtfm * f(i,1);
v(i,2) += dtfm * f(i,2);
const double dtirotate = dtfrotate / (radius(i)*radius(i)*rmass(i));
omega(i,0) += dtirotate * torque(i,0);
omega(i,1) += dtirotate * torque(i,1);

4
src/KOKKOS/fix_nve_sphere_kokkos.h
View File

@ -26,7 +26,7 @@ FixStyle(nve/sphere/kk/host,FixNVESphereKokkos<LMPHostType>)
#include "kokkos_type.h"
namespace LAMMPS_NS {
template<class DeviceType>
class FixNVESphereKokkos : public FixNVESphere {
public:
@ -36,7 +36,7 @@ class FixNVESphereKokkos : public FixNVESphere {
void init();
void initial_integrate(int);
void final_integrate();
KOKKOS_INLINE_FUNCTION
void initial_integrate_item(const int i) const;
KOKKOS_INLINE_FUNCTION

4
src/KOKKOS/kokkos.h
View File

@ -46,10 +46,10 @@ class KokkosLMP : protected Pointers {
int need_dup()
{
int value = 0;
if (neighflag == HALFTHREAD)
value = NeedDup<HALFTHREAD,DeviceType>::value;
return value;
}

14
src/KOKKOS/pair_gran_hooke_history_kokkos.cpp
View File

@ -331,7 +331,7 @@ void PairGranHookeHistoryKokkos<DeviceType>::operator()(TagPairGranHookeHistoryC
const LMP_FLOAT imass = rmass[i];
const LMP_FLOAT irad = radius[i];
const int jnum = d_numneigh_touch[i];
F_FLOAT fx_i = 0.0;
F_FLOAT fy_i = 0.0;
F_FLOAT fz_i = 0.0;
@ -343,7 +343,7 @@ void PairGranHookeHistoryKokkos<DeviceType>::operator()(TagPairGranHookeHistoryC
for (int jj = 0; jj < jnum; jj++) {
const int m = d_neighbors_touch(i, jj);
const int j = d_neighbors(i, m) & NEIGHMASK;
const X_FLOAT delx = xtmp - x(j,0);
const X_FLOAT dely = ytmp - x(j,1);
const X_FLOAT delz = ztmp - x(j,2);
@ -363,16 +363,16 @@ void PairGranHookeHistoryKokkos<DeviceType>::operator()(TagPairGranHookeHistoryC
V_FLOAT vr1 = v(i,0) - v(j,0);
V_FLOAT vr2 = v(i,1) - v(j,1);
V_FLOAT vr3 = v(i,2) - v(j,2);
// normal component
V_FLOAT vnnr = vr1*delx + vr2*dely + vr3*delz;
V_FLOAT vn1 = delx*vnnr * rsqinv;
V_FLOAT vn2 = dely*vnnr * rsqinv;
V_FLOAT vn3 = delz*vnnr * rsqinv;
// tangential component
V_FLOAT vt1 = vr1 - vn1;
V_FLOAT vt2 = vr2 - vn2;
V_FLOAT vt3 = vr3 - vn3;
@ -467,7 +467,7 @@ void PairGranHookeHistoryKokkos<DeviceType>::operator()(TagPairGranHookeHistoryC
torquex_i -= irad*tor1;
torquey_i -= irad*tor2;
torquez_i -= irad*tor3;
if (NEWTON_PAIR || j < nlocal) {
a_f(j,0) -= fx;
a_f(j,1) -= fy;
@ -550,7 +550,7 @@ void PairGranHookeHistoryKokkos<DeviceType>::ev_tally_xyz_atom(EV_FLOAT &ev, int
F_FLOAT fx, F_FLOAT fy, F_FLOAT fz,
X_FLOAT delx, X_FLOAT dely, X_FLOAT delz) const
{
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
F_FLOAT v[6];

6
src/KOKKOS/pair_gran_hooke_history_kokkos.h
View File

@ -30,7 +30,7 @@ namespace LAMMPS_NS {
template <class DeviceType>
class FixNeighHistoryKokkos;
template<int NEIGHFLAG, int NEWTON_PAIR, int EVFLAG, int SHEARUPDATE>
struct TagPairGranHookeHistoryCompute {};
@ -68,7 +68,7 @@ class PairGranHookeHistoryKokkos : public PairGranHookeHistory {
void ev_tally_xyz_atom(EV_FLOAT &ev, int i, int j,
F_FLOAT fx, F_FLOAT fy, F_FLOAT fz,
X_FLOAT delx, X_FLOAT dely, X_FLOAT delz) const;
protected:
typename AT::t_x_array_randomread x;
typename AT::t_x_array c_x;
@ -96,7 +96,7 @@ class PairGranHookeHistoryKokkos : public PairGranHookeHistory {
typename AT::t_neighbors_2d d_neighbors_touch;
typename AT::t_int_1d d_numneigh_touch;
int newton_pair;
double special_lj[4];

2
src/KOKKOS/pair_reaxc_kokkos.cpp
View File

@ -854,7 +854,7 @@ void PairReaxCKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
//dup_Cdbo = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_Cdbo);
//dup_Cdbopi = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_Cdbopi);
//dup_Cdbopi2 = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_Cdbopi2);
} else {
} else {
ndup_CdDelta = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_CdDelta);
//ndup_Cdbo = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_Cdbo);
//ndup_Cdbopi = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated>(d_Cdbopi);

2
src/KSPACE/ewald.cpp
View File

@ -72,7 +72,7 @@ void Ewald::settings(int narg, char **arg)
if (narg != 1) error->all(FLERR,"Illegal kspace_style ewald command");
accuracy_relative = fabs(force->numeric(FLERR,arg[0]));
}
}
/* ----------------------------------------------------------------------
free all memory

2
src/LATTE/fix_latte.cpp
View File

@ -299,7 +299,7 @@ void FixLatte::post_force(int vflag)
if (coulomb) forces = &flatte[0][0];
else forces = &atom->f[0][0];
int maxiter = -1;
latte(flags,&natoms,coords,type,&ntypes,mass,boxlo,boxhi,&domain->xy,
&domain->xz,&domain->yz,forces,&maxiter,&latte_energy,
&atom->v[0][0],&update->dt,virial,&newsystem,&latteerror);

14
src/MANYBODY/pair_atm.cpp
View File

@ -272,7 +272,7 @@ double PairATM::init_one(int i, int j)
int ntypes = atom->ntypes;
for (int k = j; k <= ntypes; k++)
nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] =
nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] =
nu[i][j][k];
return cut_global;
@ -290,8 +290,8 @@ void PairATM::write_restart(FILE *fp)
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j])
for (k = j; k <= atom->ntypes; k++)
if (setflag[i][j])
for (k = j; k <= atom->ntypes; k++)
fwrite(&nu[i][j][k],sizeof(double),1,fp);
}
}
@ -361,12 +361,12 @@ void PairATM::interaction_ddd(double nu, double r6,
rrk = rjk[0]*rik[0] + rjk[1]*rik[1] + rjk[2]*rik[2];
rrr = 5.0*rri*rrj*rrk;
for (int i = 0; i < 3; i++) {
fj[i] = rrj*(rrk - rri)*rik[i] -
(rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] +
fj[i] = rrj*(rrk - rri)*rik[i] -
(rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] +
(rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
fj[i] *= 3.0*r5inv;
fk[i] = rrk*(rri + rrj)*rij[i] +
(rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] +
fk[i] = rrk*(rri + rrj)*rij[i] +
(rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] +
(rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
fk[i] *= 3.0*r5inv;
}

2
src/MANYBODY/pair_atm.h
View File

@ -43,7 +43,7 @@ class PairATM : public Pair {
double ***nu;
void allocate();
void interaction_ddd(double, double, double, double, double, double *,
void interaction_ddd(double, double, double, double, double, double *,
double *, double *, double *, double *, int, double &);
};

2
src/MANYBODY/pair_gw.cpp
View File

@ -521,7 +521,7 @@ void PairGW::setup_params()
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement &&
k == params[m].kelement) {
if (n >= 0)
if (n >= 0)
error->all(FLERR,"Potential file has duplicate entry");
n = m;
}

6
src/MESSAGE/fix_client_md.cpp
View File

@ -40,11 +40,11 @@ enum{FORCES=1,ENERGY,PRESSURE,ERROR};
FixClientMD::FixClientMD(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (lmp->clientserver != 1)
if (lmp->clientserver != 1)
error->all(FLERR,"Fix client/md requires LAMMPS be running as a client");
if (!atom->map_style) error->all(FLERR,"Fix client/md requires atom map");
if (sizeof(tagint) != 4)
if (sizeof(tagint) != 4)
error->all(FLERR,"Fix client/md requires 4-byte atom IDs");
if (strcmp(update->unit_style,"real") == 0) units = REAL;
@ -306,7 +306,7 @@ void FixClientMD::receive_fev(int vflag)
}
eng = econvert * cs->unpack_double(ENERGY);
if (vflag) {
double *v = (double *) cs->unpack(PRESSURE);

12
src/MESSAGE/message.cpp
View File

@ -51,13 +51,13 @@ void Message::command(int narg, char **arg)
strcmp(arg[2],"mpi/two") == 0) {
if (narg != 4) error->all(FLERR,"Illegal message command");
lmp->cslib = new CSlib(clientserver-1,arg[2],arg[3],&world);
} else if (strcmp(arg[2],"mpi/one") == 0) {
if (narg != 3) error->all(FLERR,"Illegal message command");
if (!lmp->cscomm)
if (!lmp->cscomm)
error->all(FLERR,"Message mpi/one mode, but -mpi cmdline arg not used");
lmp->cslib = new CSlib(clientserver-1,arg[2],&lmp->cscomm,&world);
} else error->all(FLERR,"Illegal message command");
// perform initial handshake between client and server
@ -75,16 +75,16 @@ void Message::command(int narg, char **arg)
int *fieldID,*fieldtype,*fieldlen;
int msgID = cs->recv(nfield,fieldID,fieldtype,fieldlen);
if (msgID != 0) error->one(FLERR,"Bad initial client/server handshake");
} else {
int nfield;
int *fieldID,*fieldtype,*fieldlen;
int msgID = cs->recv(nfield,fieldID,fieldtype,fieldlen);
if (msgID != 0) error->one(FLERR,"Bad initial client/server handshake");
char *pstr = cs->unpack_string(1);
if (strcmp(pstr,arg[1]) != 0)
if (strcmp(pstr,arg[1]) != 0)
error->one(FLERR,"Mismatch in client/server protocol");
cs->send(0,0);
}
}

2
src/MESSAGE/server_mc.cpp
View File

@ -128,7 +128,7 @@ void ServerMC::loop()
update->nsteps = nsteps;
update->firststep = update->ntimestep;
update->laststep = update->ntimestep + nsteps;
update->integrate->setup(1);
update->integrate->run(nsteps);

16
src/MESSAGE/server_md.cpp
View File

@ -100,7 +100,7 @@ void ServerMD::loop()
if (msgID < 0) break;
// SETUP receive at beginning of each run
// required fields: DIM, PERIODICTY, ORIGIN, BOX,
// required fields: DIM, PERIODICTY, ORIGIN, BOX,
// NATOMS, NTYPES, TYPES, COORDS
// optional fields: others in enum above
@ -152,7 +152,7 @@ void ServerMD::loop()
} else error->all(FLERR,"Server md setup field unknown");
}
if (dim == 0 || !periodicity || !origin || !box ||
if (dim == 0 || !periodicity || !origin || !box ||
natoms < 0 || ntypes < 0 || !types || !coords)
error->all(FLERR,"Required server md setup field not received");
@ -164,8 +164,8 @@ void ServerMD::loop()
// reset box, global and local
// reset proc decomposition
if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) &&
if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) &&
domain->triclinic == 0)
error->all(FLERR,"Server md is not initialized for a triclinic box");
@ -194,7 +194,7 @@ void ServerMD::loop()
int ntotal;
MPI_Allreduce(&atom->nlocal,&ntotal,1,MPI_INT,MPI_SUM,world);
if (ntotal != natoms)
if (ntotal != natoms)
error->all(FLERR,"Server md atom count does not match client");
atom->map_init();
@ -247,7 +247,7 @@ void ServerMD::loop()
// reset global/local box like FixDeform at end_of_step()
if (origin && box) {
if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) &&
if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) &&
domain->triclinic == 0)
error->all(FLERR,"Server md is not initialized for a triclinic box");
box_change(origin,box);
@ -348,7 +348,7 @@ void ServerMD::send_fev(int msgID)
CSlib *cs = (CSlib *) lmp->cslib;
cs->send(msgID,3);
double *forces = NULL;
if (atom->nlocal) {
if (units != REAL) forces = &atom->f[0][0];
@ -370,7 +370,7 @@ void ServerMD::send_fev(int msgID)
MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
engall *= econvert;
cs->pack_double(ENERGY,engall);
double v[6],vall[6];
for (int i = 0; i < 6; i++)
v[i] = force->pair->virial[i];

2
src/POEMS/fix_poems.cpp
View File

@ -366,7 +366,7 @@ void FixPOEMS::init()
int pflag = 0;
for (i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"poems") == 0) pflag = 1;
if (pflag && (modify->fmask[i] & POST_FORCE) &&
if (pflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag) {
char str[128];
snprintf(str,128,"Fix %s alters forces after fix poems",modify->fix[i]->id);

8
src/REPLICA/fix_hyper_global.cpp
View File

@ -44,7 +44,7 @@ using namespace FixConst;
FixHyperGlobal::FixHyperGlobal(LAMMPS *lmp, int narg, char **arg) :
FixHyper(lmp, narg, arg), blist(NULL), xold(NULL), tagold(NULL)
{
if (atom->map_style == 0)
if (atom->map_style == 0)
error->all(FLERR,"Fix hyper/global command requires atom map");
if (narg != 7) error->all(FLERR,"Illegal fix hyper/global command");
@ -392,7 +392,7 @@ void FixHyperGlobal::grow_bond()
maxbond += DELTA;
if (maxbond < 0 || maxbond > MAXSMALLINT)
error->one(FLERR,"Fix hyper/local per-processor bond count is too big");
blist = (OneBond *)
blist = (OneBond *)
memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/local:blist");
}
@ -429,7 +429,7 @@ double FixHyperGlobal::compute_vector(int i)
// i = 7 = max bond length during this run
// i = 8 = cummulative hyper time since fix created
// i = 9 = cummulative # of event timesteps since fix created
// i = 9 = cummulative # of event timesteps since fix created
// i = 10 = cummulative # of atoms in events since fix created
if (i == 0) return outvec[1];
@ -478,7 +478,7 @@ double FixHyperGlobal::query(int i)
{
if (i == 1) return compute_vector(8); // cummulative hyper time
if (i == 2) return compute_vector(9); // nevent
if (i == 3) return compute_vector(10); // nevent_atom
if (i == 3) return compute_vector(10); // nevent_atom
if (i == 4) return compute_vector(4); // ave bonds/atom
if (i == 5) return compute_vector(6); // maxdrift
if (i == 6) return compute_vector(7); // maxbondlen

46
src/REPLICA/fix_hyper_local.cpp
View File

@ -54,7 +54,7 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
// local index vs global ID in same variable
// maybe need to declare them all tagint, not int
if (atom->map_style == 0)
if (atom->map_style == 0)
error->all(FLERR,"Fix hyper/local command requires atom map");
if (sizeof(tagint) != sizeof(int))
@ -81,7 +81,7 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
alpha_user = force->numeric(FLERR,arg[8]);
boosttarget = force->numeric(FLERR,arg[9]);
if (cutbond < 0.0 || qfactor < 0.0 || vmax < 0.0 ||
if (cutbond < 0.0 || qfactor < 0.0 || vmax < 0.0 ||
tequil <= 0.0 || dcut <= 0.0 || alpha_user <= 0.0 || boosttarget < 1.0)
error->all(FLERR,"Illegal fix hyper/local command");
@ -107,7 +107,7 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
histo_count = force->inumeric(FLERR,arg[iarg+2]);
histo_delta = force->numeric(FLERR,arg[iarg+3]);
histo_print = force->inumeric(FLERR,arg[iarg+4]);
if (histo_every <= 0 || histo_count % 2 ||
if (histo_every <= 0 || histo_count % 2 ||
histo_delta <= 0.0 || histo_print <= 0)
error->all(FLERR,"Illegal fix hyper/local command");
iarg += 5;
@ -276,16 +276,16 @@ void FixHyperLocal::init()
// warn if no drift distance added to cutghost
if (firstflag) {
double cutghost;
if (force->pair)
double cutghost;
if (force->pair)
cutghost = MAX(force->pair->cutforce+neighbor->skin,comm->cutghostuser);
else
else
cutghost = comm->cutghostuser;
if (cutghost < dcut)
if (cutghost < dcut)
error->all(FLERR,"Fix hyper/local bond cutoff exceeds ghost atom range - "
"use comm_modify cutoff command");
if (cutghost < dcut+cutbond/2.0 && me == 0)
if (cutghost < dcut+cutbond/2.0 && me == 0)
error->warning(FLERR,"Fix hyper/local ghost atom range "
"may not allow for atom drift between events");
}
@ -371,8 +371,8 @@ void FixHyperLocal::pre_neighbor()
missing_coeff += bonds[i][m].boostcoeff;
if (lostbond != IGNORE) {
char str[128];
sprintf(str,"Fix hyper/local bond info missing for bond "
TAGINT_FORMAT "," TAGINT_FORMAT
sprintf(str,"Fix hyper/local bond info missing for bond "
TAGINT_FORMAT "," TAGINT_FORMAT
" with coeff %g at step " BIGINT_FORMAT,
atom->tag[i],bonds[i][m].jtag,bonds[i][m].boostcoeff,
update->ntimestep);
@ -520,7 +520,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
i = old2now[iold];
emax = maxstrain[i];
for (jj = 0; jj < jnum; jj++) {
jold = jlist[jj];
j = old2now[jold];
@ -722,13 +722,13 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
if (me == 0) {
if (screen) {
fprintf(screen,"Histogram of bias coeffs:\n");
for (i = 0; i < histo_count+2; i++)
for (i = 0; i < histo_count+2; i++)
fprintf(screen," %g",1.0*allhisto[i]/total);
fprintf(screen,"\n");
}
if (logfile) {
fprintf(logfile,"Histogram of bias coeffs:\n");
for (i = 0; i < histo_count+2; i++)
for (i = 0; i < histo_count+2; i++)
fprintf(logfile," %g",1.0*allhisto[i]/total);
fprintf(logfile,"\n");
}
@ -790,7 +790,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
if (checkcoeff && update->ntimestep % checkcoeff_every == 0) {
int jb,jbonds;
for (i = 0; i < nlocal; i++) {
nbond = numbond[i];
for (m = 0; m < nbond; m++) {
@ -802,9 +802,9 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
jbonds = numbond[j];
for (jb = 0; jb < jbonds; jb++)
if (bonds[j][jb].jtag == itag) break;
if (jb == jbonds)
if (jb == jbonds)
error->one(FLERR,"Fix hyper/local could not find duplicate bond");
if (bonds[i][m].boostcoeff != bonds[j][jb].boostcoeff)
if (bonds[i][m].boostcoeff != bonds[j][jb].boostcoeff)
checkcoeff_count++;
}
}
@ -885,7 +885,7 @@ void FixHyperLocal::build_bond_list(int natom)
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
while (1) {
bonds = (OneBond **) memory->create(bonds,nmax,maxbondperatom,
"hyper/local:bonds");
@ -941,7 +941,7 @@ void FixHyperLocal::build_bond_list(int natom)
nbond++;
continue;
}
bonds[i][nbond].r0 = sqrt(rsq);
bonds[i][nbond].jtag = tag[j];
bonds[i][nbond].j = j;
@ -998,7 +998,7 @@ void FixHyperLocal::build_bond_list(int natom)
/* ---------------------------------------------------------------------- */
int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf,
int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf,
int /* pbc_flag */, int * /* pbc */)
{
int i,j,m;
@ -1233,7 +1233,7 @@ double FixHyperLocal::compute_vector(int i)
// i = 9 = average bias potential for all bonds during this run
// i = 10 = max bias potential for any bond during this run
// i = 11 = min bias potential for any bond during this run
// i = 12 = max dist from my box of any ghost atom with
// i = 12 = max dist from my box of any ghost atom with
// maxstain < qfactor during this run
// i = 13 = max dist from my box of any ghost atom with
// any maxstrain during this run
@ -1245,7 +1245,7 @@ double FixHyperLocal::compute_vector(int i)
// i = 18 = count of non-matching bias coefficients found during this run
// i = 19 = cummulative hyper time
// i = 20 = cummulative # of event timesteps since fix created
// i = 20 = cummulative # of event timesteps since fix created
// i = 21 = cummulative # of atoms in events since fix created
// i = 22 = cummulative # of new bonds formed since fix created
@ -1394,7 +1394,7 @@ double FixHyperLocal::query(int i)
{
if (i == 1) return compute_vector(19); // cummulative hyper time
if (i == 2) return compute_vector(20); // nevent
if (i == 3) return compute_vector(21); // nevent_atom
if (i == 3) return compute_vector(21); // nevent_atom
if (i == 4) return compute_vector(3); // ave bonds/atom
if (i == 5) return compute_vector(6); // maxdrift
if (i == 6) return compute_vector(7); // maxbondlen

4
src/REPLICA/fix_hyper_local.h
View File

@ -82,7 +82,7 @@ class FixHyperLocal : public FixHyper {
double allboost; // sum of boostcoeff on all bonds on this step
int nnewbond; // running tally of number of new bonds created
int maxbondperatom; // max # of bonds any atom ever has
int maxbondperatom; // max # of bonds any atom ever has
int commflag; // flag for communication mode
int nevent; // # of events that trigger bond rebuild
int nevent_atom; // # of atoms that experienced an event
@ -115,7 +115,7 @@ class FixHyperLocal : public FixHyper {
double **xold; // coords of owned+ghost atoms when bonds created
tagint *tagold; // global IDs of owned+ghost atoms when b created
int maxold; // allocated size of old2now
int maxbond; // allocated size of bonds
int old_nall; // nlocal+nghost when old2now was last setup

10
src/REPLICA/hyper.cpp
View File

@ -40,7 +40,7 @@ enum{NOHYPER,GLOBAL,LOCAL};
/* ---------------------------------------------------------------------- */
Hyper::Hyper(LAMMPS *lmp) :
Hyper::Hyper(LAMMPS *lmp) :
Pointers(lmp), dumplist(NULL)
{}
@ -73,7 +73,7 @@ void Hyper::command(int narg, char **arg)
// total # of timesteps must be multiple of t_event
if (t_event <= 0)
if (t_event <= 0)
error->all(FLERR,"Invalid t_event in hyper command");
if (nsteps % t_event)
error->all(FLERR,"Hyper nsteps must be multiple of t_event");
@ -209,7 +209,7 @@ void Hyper::command(int narg, char **arg)
dynamics(t_event,time_dynamics);
fix_event->store_state_quench();
quench(0);
ecount = compute_event->all_events();
if (ecount) {
@ -226,7 +226,7 @@ void Hyper::command(int narg, char **arg)
fix_event->store_event();
if (hyperenable) fix_hyper->build_bond_list(ecount);
}
fix_event->restore_state_quench();
istep = update->ntimestep - update->beginstep;
}
@ -519,7 +519,7 @@ void Hyper::options(int narg, char **arg)
if (iarg+2 > narg) error->all(FLERR,"Illegal hyper command");
dumpflag = 1;
int idump = output->find_dump(arg[iarg+1]);
if (idump < 0)
if (idump < 0)
error->all(FLERR,"Dump ID in hyper command does not exist");
memory->grow(dumplist,ndump+1,"hyper:dumplist");
dumplist[ndump++] = idump;

6
src/RIGID/fix_ehex.h
View File

@ -75,7 +75,7 @@ class FixEHEX : public Fix {
/* ERROR/WARNING messages:
E: Illegal fix ehex command: wrong number of parameters
E: Illegal fix ehex command: wrong number of parameters
UNDOCUMENTED
@ -87,11 +87,11 @@ E: Region ID for fix ehex does not exist
Self-explanatory.
E: Illegal fix ehex keyword
E: Illegal fix ehex keyword
UNDOCUMENTED
E: You can only use the keyword 'com' together with the keyword 'constrain'
E: You can only use the keyword 'com' together with the keyword 'constrain'
UNDOCUMENTED

8
src/RIGID/fix_rattle.h
View File

@ -79,7 +79,7 @@ class FixRattle : public FixShake {
/* ERROR/WARNING messages:
W: Fix rattle should come after all other integration fixes
W: Fix rattle should come after all other integration fixes
UNDOCUMENTED
@ -88,15 +88,15 @@ E: Rattle determinant = 0.0
The determinant of the matrix being solved for a single cluster
specified by the fix rattle command is numerically invalid.
E: Rattle failed
E: Rattle failed
UNDOCUMENTED
E: Coordinate constraints are not satisfied up to desired tolerance
E: Coordinate constraints are not satisfied up to desired tolerance
UNDOCUMENTED
E: Velocity constraints are not satisfied up to desired tolerance
E: Velocity constraints are not satisfied up to desired tolerance
UNDOCUMENTED

6
src/RIGID/fix_rigid.cpp
View File

@ -174,7 +174,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
MPI_Allreduce(&vmin,&minval,1,MPI_INT,MPI_MIN,world);
molecule = new tagint[nlocal];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mask[i] & groupbit)
molecule[i] = (tagint)((tagint)value[i] - minval + 1);
delete[] value;
} else error->all(FLERR,"Unsupported fix rigid custom property");
@ -727,7 +727,7 @@ void FixRigid::init()
int rflag = 0;
for (i = 0; i < modify->nfix; i++) {
if (modify->fix[i]->rigid_flag) rflag = 1;
if (rflag && (modify->fmask[i] & POST_FORCE) &&
if (rflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag) {
char str[128];
snprintf(str,128,"Fix %s alters forces after fix rigid",modify->fix[i]->id);
@ -2639,7 +2639,7 @@ int FixRigid::modify_param(int narg, char **arg)
else error->all(FLERR,"Illegal fix_modify command");
// reset fix mask
// must do here and not in init,
// must do here and not in init,
// since modify.cpp::init() uses fix masks before calling fix::init()
for (int i = 0; i < modify->nfix; i++)

4
src/RIGID/fix_rigid_small.cpp
View File

@ -566,7 +566,7 @@ void FixRigidSmall::init()
int rflag = 0;
for (i = 0; i < modify->nfix; i++) {
if (modify->fix[i]->rigid_flag) rflag = 1;
if (rflag && (modify->fmask[i] & POST_FORCE) &&
if (rflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag) {
char str[128];
snprintf(str,128,"Fix %s alters forces after fix rigid",modify->fix[i]->id);
@ -3430,7 +3430,7 @@ int FixRigidSmall::modify_param(int narg, char **arg)
else error->all(FLERR,"Illegal fix_modify command");
// reset fix mask
// must do here and not in init,
// must do here and not in init,
// since modify.cpp::init() uses fix masks before calling fix::init()
for (int i = 0; i < modify->nfix; i++)

2
src/SPIN/pair_spin_dmi.cpp
View File

@ -414,7 +414,7 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
compute the mechanical force due to the dmi interaction between atom i and atom j
------------------------------------------------------------------------- */
void PairSpinDmi::compute_dmi_mech(int i, int j, double rsq, double /*eij*/[3],
void PairSpinDmi::compute_dmi_mech(int i, int j, double rsq, double /*eij*/[3],
double fi[3], double spi[3], double spj[3])
{
int *type = atom->type;

4
src/SPIN/pair_spin_exchange.cpp
View File

@ -134,7 +134,7 @@ void PairSpinExchange::coeff(int narg, char **arg)
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
}
@ -395,7 +395,7 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
compute the mechanical force due to the exchange interaction between atom i and atom j
------------------------------------------------------------------------- */
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double eij[3],
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double eij[3],
double fi[3], double spi[3], double spj[3])
{
int *type = atom->type;

2
src/USER-INTEL/fix_intel.cpp
View File

@ -456,7 +456,7 @@ void FixIntel::pair_init_check(const bool cdmessage)
" exclusions with Intel");
}
#endif
int need_tag = 0;
if (atom->molecular) need_tag = 1;

12
src/USER-INTEL/fix_intel.h
View File

@ -83,7 +83,7 @@ class FixIntel : public Fix {
}
inline void set_reduce_flag() { if (_nthreads > 1) _need_reduce = 1; }
inline int lrt() {
if (force->kspace_match("pppm/intel", 0) && update->whichflag == 1)
if (force->kspace_match("pppm/intel", 0) && update->whichflag == 1)
return _lrt;
else return 0;
}
@ -104,7 +104,7 @@ class FixIntel : public Fix {
int _pair_intel_count, _pair_hybrid_flag;
// These should be removed in subsequent update w/ simpler hybrid arch
int _pair_hybrid_zero, _hybrid_nonpair, _zero_master;
public:
inline int* get_overflow_flag() { return _overflow_flag; }
inline int* get_off_overflow_flag() { return _off_overflow_flag; }
@ -215,7 +215,7 @@ class FixIntel : public Fix {
_alignvar(double _stopwatch_offload_pair[1],64);
void _sync_main_arrays(const int prereverse);
template <class ft>
void reduce_results(ft * _noalias const f_in);
@ -512,13 +512,13 @@ issues. Please use 14.0.1.106 or 15.1.133 or later.
E: Currently, cannot offload more than one intel style with hybrid.
Currently, when using offload, hybrid pair styles can only use the intel
Currently, when using offload, hybrid pair styles can only use the intel
suffix for one of the pair styles.
E: Cannot yet use hybrid styles with Intel offload.
The hybrid pair style configuration is not yet supported when using offload
within the Intel package. Support is limited to hybrid/overlay or a hybrid
within the Intel package. Support is limited to hybrid/overlay or a hybrid
style that does not require a skip list.
W: Leaving a core/node free can improve performance for offload
@ -564,7 +564,7 @@ atoms throughout the simulation.
E: Intel package requires fdotr virial with newton on.
This error can occur with a hybrid pair style that mixes styles that are
incompatible with the newton pair setting turned on. Try turning the
incompatible with the newton pair setting turned on. Try turning the
newton pair setting off.
E: Add -DLMP_INTEL_NBOR_COMPAT to build for special_bond exclusions with Intel

4
src/USER-INTEL/intel_buffers.cpp
View File

@ -255,7 +255,7 @@ void IntelBuffers<flt_t, acc_t>::free_list_local()
#endif
lmp->memory->destroy(cnumneigh);
}
#ifdef _LMP_INTEL_OFFLOAD
if (_off_map_ilist != NULL) {
const int * ilist = _off_map_ilist;
@ -295,7 +295,7 @@ void IntelBuffers<flt_t, acc_t>::grow_data3(NeighList *list,
{
const int size = list->get_maxlocal();
int list_num;
for (list_num = 0; list_num < _n_list_ptrs; list_num++)
for (list_num = 0; list_num < _n_list_ptrs; list_num++)
if (_neigh_list_ptrs[list_num].list_ptr == (void*)list) break;
if (list_num == _n_list_ptrs) {
if (_n_list_ptrs == _max_list_ptrs) {

2
src/USER-INTEL/intel_preprocess.h
View File

@ -146,7 +146,7 @@ enum {TIME_PACK, TIME_HOST_NEIGHBOR, TIME_HOST_PAIR, TIME_OFFLOAD_NEIGHBOR,
}
#else
#define IP_PRE_edge_align(n, esize) \
#endif

4
src/USER-INTEL/npair_full_bin_ghost_intel.cpp
View File

@ -482,7 +482,7 @@ void NPairFullBinGhostIntel::fbi(const int offload, NeighList * list,
#endif
}
}
#ifndef _LMP_INTEL_OFFLOAD
if (exclude) {
int alln = n;
@ -515,7 +515,7 @@ void NPairFullBinGhostIntel::fbi(const int offload, NeighList * list,
}
}
#endif
int ns = n - maxnbors;
int alln = n;
atombin[i] = ns;

4
src/USER-INTEL/npair_halffull_newton_intel.cpp
View File

@ -61,7 +61,7 @@ void NPairHalffullNewtonIntel::build_t(NeighList *list,
const int * _noalias const numneigh_full = list->listfull->numneigh;
const int ** _noalias const firstneigh_full =
(const int ** const)list->listfull->firstneigh;
#if defined(_OPENMP)
#pragma omp parallel
#endif
@ -148,7 +148,7 @@ void NPairHalffullNewtonIntel::build_t3(NeighList *list, int *numhalf)
const int * _noalias const numneigh_full = numhalf;
const int ** _noalias const firstneigh_full =
(const int ** const)list->listfull->firstneigh;
int packthreads = 1;
if (comm->nthreads > INTEL_HTHREADS) packthreads = comm->nthreads;

16
src/USER-INTEL/npair_intel.cpp
View File

@ -72,7 +72,7 @@ void NPairIntel::copy_cutsq_info(IntelBuffers<flt_t,acc_t> *buffers) {
if (cutneighghostsq)
use_ghost_cut = 1;
buffers->set_ntypes(tp1, use_ghost_cut);
flt_t **cutneighsqb = buffers->get_cutneighsq();
for (int i = 1; i <= atom->ntypes; i++)
for (int j = 1; j <= atom->ntypes; j++)
@ -116,7 +116,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
#ifdef _LMP_INTEL_OFFLOAD
if (offload_noghost && offload) nall_t = atom->nlocal;
#endif
const int pack_width = _fix->nbor_pack_width();
const ATOM_T * _noalias const x = buffers->get_x();
@ -146,7 +146,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
int ** _noalias const firstneigh = list->firstneigh;
int * _noalias const numneigh = list->numneigh;
int * _noalias const cnumneigh = buffers->cnumneigh();
const int nstencil = this->nstencil;
const int * _noalias const stencil = this->stencil;
const flt_t * _noalias const cutneighsq = buffers->get_cutneighsq()[0];
@ -204,7 +204,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
}
}
const int special_bound = sb;
#ifdef _LMP_INTEL_OFFLOAD
const int * _noalias const binhead = this->binhead;
const int * _noalias const bins = this->bins;
@ -547,7 +547,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
j = -j - 1;
} else
ofind_special(which, special, nspecial, i, tag[j]);
if (which) {
j = j ^ (which << SBBITS);
if (which < special_bound) addme = 0;
@ -562,7 +562,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
if (THREE) {
alln = n2;
n2 = pack_offset + maxnbors;
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#ifdef LMP_INTEL_NBOR_COMPAT
@ -592,7 +592,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
}
}
}
#ifndef _LMP_INTEL_OFFLOAD
if (exclude) {
neighptr2 = neighptr;
@ -643,7 +643,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
n += pack_width;
}
#endif
for (int u = pack_offset + maxnbors; u < n2; u++) {
#ifdef LMP_INTEL_3BODY_FAST
neighptr[n] = neighptr2[u];

4
src/USER-INTEL/npair_intel.h
View File

@ -76,7 +76,7 @@ class NPairIntel : public NPair {
NPairIntel(class LAMMPS *);
~NPairIntel();
virtual void copy_neighbor_info();
#ifdef _LMP_INTEL_OFFLOAD
void grow_stencil();
#endif
@ -86,7 +86,7 @@ class NPairIntel : public NPair {
template <class flt_t, class acc_t>
void copy_cutsq_info(IntelBuffers<flt_t,acc_t> *);
template <class flt_t, class acc_t, int, int, int, int, int>
void bin_newton(const int, NeighList *, IntelBuffers<flt_t,acc_t> *,
const int, const int, const int offload_end = 0);

4
src/USER-INTEL/npair_skip_intel.cpp
View File

@ -109,7 +109,7 @@ void NPairSkipIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
int my_inum = ifrom;
_inum_starts[tid] = ifrom;
// loop over parent full list
for (int ii = ifrom; ii < ito; ii++) {
const int i = ilist_skip[ii];
@ -131,7 +131,7 @@ void NPairSkipIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
const int j = joriginal & NEIGHMASK;
if (!ijskip[itype][type[j]]) neighptr[n++] = joriginal;
}
numhalf[my_inum] = n;
numhalf[my_inum] = n;
for (int jj = jnumhalf; jj < jnum; jj++) {
const int joriginal = jlist[jj];

2
src/USER-INTEL/npair_skip_intel.h
View File

@ -46,7 +46,7 @@ class NPairSkipIntel : public NPair {
~NPairSkipIntel();
virtual void copy_neighbor_info();
void build(class NeighList *);
protected:
FixIntel *_fix;
int *_inum_starts, *_inum_counts, *_full_props;

2
src/USER-INTEL/pair_airebo_intel.h
View File

@ -110,7 +110,7 @@ path and name are correct.
E: Cannot yet use airebo/intel with hybrid.
Pair style airebo/intel cannot currently be used as part of a hybrid
pair style (with the exception of hybrid/overlay).
pair style (with the exception of hybrid/overlay).
*/

2
src/USER-INTEL/pair_buck_coul_cut_intel.cpp
View File

@ -204,7 +204,7 @@ void PairBuckCoulCutIntel::eval(const int offload, const int vflag,
acc_t oevdwl, oecoul, ov0, ov1, ov2, ov3, ov4, ov5;
if (EFLAG || vflag)
oevdwl = oecoul = ov0 = ov1 = ov2 = ov3 = ov4 = ov5 = (acc_t)0;
if (NEWTON_PAIR == 0 && inum != nlocal)
if (NEWTON_PAIR == 0 && inum != nlocal)
memset(f_start, 0, f_stride * sizeof(FORCE_T));
// loop over neighbors of my atoms

2
src/USER-INTEL/pair_dpd_intel.cpp
View File

@ -521,7 +521,7 @@ void PairDPDIntel::pack_force_const(ForceConst<flt_t> &fc,
}
}
if (mytypes > 1 || atom->molecular) _onetype = 0;
for (int i = 0; i < 4; i++) {
fc.special_lj[i] = force->special_lj[i];
fc.special_lj[0] = 1.0;

4
src/USER-INTEL/pppm_intel.cpp
View File

@ -586,7 +586,7 @@ void PPPMIntel::fieldforce_ik(IntelBuffers<flt_t,acc_t> *buffers)
if (force->newton_pair) zl += atom->nghost;
memset(f, 0, zl * sizeof(FORCE_T));
}
#if defined(_OPENMP)
#pragma omp parallel default(none) \
shared(nlocal, nthr) if(!_use_lrt)
@ -737,7 +737,7 @@ void PPPMIntel::fieldforce_ad(IntelBuffers<flt_t,acc_t> *buffers)
if (force->newton_pair) zl += atom->nghost;
memset(f, 0, zl * sizeof(FORCE_T));
}
#if defined(_OPENMP)
#pragma omp parallel default(none) \
shared(nlocal, nthr) if(!_use_lrt)

8
src/USER-MEAMC/meam_impl.cpp
View File

@ -34,12 +34,12 @@ MEAM::MEAM(Memory* mem)
maxneigh = 0;
scrfcn = dscrfcn = fcpair = NULL;
neltypes = 0;
for (int i = 0; i < maxelt; i++) {
Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
beta1_meam[i]= beta2_meam[i] = beta3_meam[i] =
t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] =
Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
beta1_meam[i]= beta2_meam[i] = beta3_meam[i] =
t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] =
rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;
for (int j = 0; j < maxelt; j++) {
lattce_meam[i][j] = FCC;

38
src/USER-MISC/compute_stress_mop.h
View File

@ -56,47 +56,47 @@ namespace LAMMPS_NS {
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute stress/mop incompatible with simulation dimension
Compute stress/mop only works with 3D simulations.
E: Compute stress/mop incompatible with triclinic simulation box
Self-explanatory.
E: Compute stress/mop requires a fixed simulation box
Compute stress/mop is not compatible with any change of volume or shape
or boundary conditions of the simulation box.
E: No pair style is defined for compute stress/mop
Self-explanatory. Compute stress/mop requires the definition of a pair style.
E: Pair style does not support compute stress/mop
The pair style does not have a single() function, so it can
not be invoked by compute stress/mop.
W: compute stress/mop does not account for bond potentials
W: compute stress/mop does not account for angle potentials
W: compute stress/mop does not account for dihedral potentials
W: compute stress/mop does not account for improper potentials
W: compute stress/mop does not account for kspace contributions
Compute stress/mop only accounts for pairwise additive interactions for
the computation of local stress tensor components.
*/

36
src/USER-MISC/compute_stress_mop_profile.h
View File

@ -70,43 +70,43 @@ namespace LAMMPS_NS {
command-line option when running LAMMPS to see the offending line.
E: Compute stress/mop/profile incompatible with simulation dimension
Compute stress/mop/profile only works with 3D simulations.
E: Compute stress/mop/profile incompatible with triclinic simulation box
Self-explanatory.
E: Compute stress/mop/profile requires a fixed simulation box
Compute stress/mop/profile is not compatible with any change of volume or shape
or boundary conditions of the simulation box.
E: No pair style is defined for compute stress/mop/profile
Self-explanatory. Compute stress/mop/profile requires the definition of a pair style.
E: Pair style does not support compute stress/mop/profile
The pair style does not have a single() function, so it can
not be invoked by compute stress/mop/profile.
E: Origin of bins for compute stress/mop/profile is out of bounds
Self-explanatory.
W: compute stress/mop/profile does not account for bond potentials
W: compute stress/mop/profile does not account for angle potentials
W: compute stress/mop/profile does not account for dihedral potentials
W: compute stress/mop/profile does not account for improper potentials
W: compute stress/mop/profile does not account for kspace contributions
Compute stress/mop/profile only accounts for pairwise additive interactions for
the computation of local stress tensor components.
*/

8
src/USER-MISC/dihedral_table_cut.cpp
View File

@ -454,7 +454,7 @@ DihedralTableCut::~DihedralTableCut()
void DihedralTableCut::compute(int eflag, int vflag)
{
int i1,i2,i3,i4,i,j,k,n,type;
double edihedral;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
@ -704,12 +704,12 @@ void DihedralTableCut::compute(int eflag, int vflag)
double gptt = 0;
if ( acos(costh12) > aat_theta0_1[type]) {
gt *= 1-da1*da1/dtheta/dtheta;
gt *= 1-da1*da1/dtheta/dtheta;
gpt = -aat_k[type]*2*da1/dtheta/dtheta;
}
if ( acos(costh23) > aat_theta0_1[type]) {
gtt *= 1-da2*da2/dtheta/dtheta;
gtt *= 1-da2*da2/dtheta/dtheta;
gptt = -aat_k[type]*2*da2/dtheta/dtheta;
}
@ -717,7 +717,7 @@ void DihedralTableCut::compute(int eflag, int vflag)
for (i = 0; i < 4; i++)
for (j = 0; j < 3; j++)
fabcd[i][j] -= - gt*gtt*fpphi*dphidr[i][j]
fabcd[i][j] -= - gt*gtt*fpphi*dphidr[i][j]
- gt*gptt*fphi*dthetadr[1][i][j] + gpt*gtt*fphi*dthetadr[0][i][j];
// apply force to each of 4 atoms

14
src/USER-MISC/pair_kolmogorov_crespi_z.cpp
View File

@ -85,7 +85,7 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
double rsq,r,rhosq,exp1,exp2,r6,r8;
double frho,sumC,sumC2,sumCff,fsum,rdsq;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
@ -121,7 +121,7 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
// rho^2 = r^2 - (n,r) = r^2 - z^2
rhosq = delx*delx + dely*dely;
rsq = rhosq + delz*delz;
if (rsq < cutsq[itype][jtype]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
@ -137,7 +137,7 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
exp2 = exp(-rdsq);
// note that f(rho_ij) equals f(rho_ji) as normals are all along z
sumC = p.C0+p.C2*rdsq+p.C4*rdsq*rdsq;
sumC = p.C0+p.C2*rdsq+p.C4*rdsq*rdsq;
sumC2 = (2*p.C2+4*p.C4*rdsq)*p.delta2inv;
frho = exp2*sumC;
sumCff = p.C + 2*frho;
@ -221,9 +221,9 @@ void PairKolmogorovCrespiZ::settings(int narg, char **arg)
void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
{
int i,j,n;
int i,j,n;
if (narg != 3 + atom->ntypes)
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
@ -262,7 +262,7 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
read_file(arg[2]);
double cut_one = cut_global;
int count = 0;
@ -407,7 +407,7 @@ void PairKolmogorovCrespiZ::read_file(char *filename)
params[nparams].S = atof(words[10]);
// energies in meV further scaled by S
double meV = 1.0e-3*params[nparams].S;
double meV = 1.0e-3*params[nparams].S;
params[nparams].C *= meV;
params[nparams].A *= meV;
params[nparams].C0 *= meV;

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

@ -82,7 +82,7 @@ void PairLJMDF::compute(int eflag, int vflag)
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -109,7 +109,7 @@ void PairLJMDF::compute(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_inner_sq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
rr = sqrt(rsq);
dp = (cut[itype][jtype] - cut_inner[itype][jtype]);

38
src/USER-MOFFF/pair_buck6d_coul_gauss_dsf.cpp
View File

@ -142,7 +142,7 @@ void PairBuck6dCoulGaussDSF::compute(int eflag, int vflag)
term3 = term2*term2;
term4 = 1.0/(1.0 + term2);
term5 = 1.0/(1.0 + 2.0*term2 + term3);
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
@ -150,11 +150,11 @@ void PairBuck6dCoulGaussDSF::compute(int eflag, int vflag)
if (rsq > rsmooth_sq[itype][jtype]) {
rcu = r*rsq;
rqu = rsq*rsq;
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
// forcebuck6d is -dE/dr*r
// forcebuck6d is -dE/dr*r
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
ebuck6d *= sme;
}
@ -167,10 +167,10 @@ void PairBuck6dCoulGaussDSF::compute(int eflag, int vflag)
arg = alpha_ij[itype][jtype]*r;
erfcd = MathSpecial::expmsq(arg);
erfcc = 1 - (MathSpecial::my_erfcx(arg) * erfcd);
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
r*f_shift_ij[itype][jtype]) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else forcecoul = 0.0;
@ -191,7 +191,7 @@ void PairBuck6dCoulGaussDSF::compute(int eflag, int vflag)
} else evdwl = 0.0;
if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
ecoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
rsq*f_shift_ij[itype][jtype]);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
@ -317,7 +317,7 @@ void PairBuck6dCoulGaussDSF::init_style()
neighbor->request(this,instance_me);
cut_coulsq = cut_coul * cut_coul;
cut_coulsq = cut_coul * cut_coul;
}
/* ----------------------------------------------------------------------
@ -330,7 +330,7 @@ double PairBuck6dCoulGaussDSF::init_one(int i, int j)
double cut = MAX(cut_lj[i][j],cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
//calculation of smoothing coefficients c0-c5
c0[i][j] = c1[i][j] = c2[i][j] = c3[i][j] = c4[i][j] = c5[i][j] = 0.0;
rsmooth_sq[i][j] = cut_ljsq[i][j];
@ -347,7 +347,7 @@ double PairBuck6dCoulGaussDSF::init_one(int i, int j)
c5[i][j] = -6.0/denom;
rsmooth_sq[i][j] = rsm_sq;
}
// if offset_flag, shift is only invoked if there is not already smoothing
if (offset_flag && vdwl_smooth >= 1.0) {
double term1 = buck6d3[i][j]/pow(cut_lj[i][j],6.0);
@ -355,12 +355,12 @@ double PairBuck6dCoulGaussDSF::init_one(int i, int j)
double rexp = exp(-cut_lj[i][j]*buck6d2[i][j]);
offset[i][j] = buck6d1[i][j]*rexp - term1*term4;
} else offset[i][j] = 0.0;
double erfcd_c = exp(-alpha_ij[i][j]*alpha_ij[i][j]*cut_coul*cut_coul);
double erfcc_c = erf(alpha_ij[i][j]*cut_coul);
f_shift_ij[i][j] = -erfcc_c/cut_coulsq + 2.0/MY_PIS*alpha_ij[i][j]*erfcd_c/cut_coul;
e_shift_ij[i][j] = erfcc_c/cut_coul - f_shift_ij[i][j]*cut_coul;
cut_ljsq[j][i] = cut_ljsq[i][j];
alpha_ij[j][i] = alpha_ij[i][j];
f_shift_ij[j][i] = f_shift_ij[i][j];
@ -520,18 +520,18 @@ double PairBuck6dCoulGaussDSF::single(int i, int j, int itype, int jtype, double
term3 = term2*term2;
term4 = 1.0/(1.0 + term2);
term5 = 1.0/(1.0 + 2.0*term2 + term3);
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
// smoothing term
if (rsq > rsmooth_sq[itype][jtype]) {
rcu = r*rsq;
rqu = rsq*rsq;
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
// forcebuck6d is -dE/dr*r
// forcebuck6d is -dE/dr*r
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
ebuck6d *= sme;
}
@ -542,7 +542,7 @@ double PairBuck6dCoulGaussDSF::single(int i, int j, int itype, int jtype, double
arg = alpha_ij[itype][jtype]*r;
erfcd = MathSpecial::expmsq(arg);
erfcc = 1 - (MathSpecial::my_erfcx(arg) * erfcd);
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
r*f_shift_ij[itype][jtype]) * r;
} else forcecoul = 0.0;
@ -555,7 +555,7 @@ double PairBuck6dCoulGaussDSF::single(int i, int j, int itype, int jtype, double
}
if (rsq < cut_coulsq) {
phicoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
phicoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
rsq*f_shift_ij[itype][jtype]);
eng += phicoul;
}

36
src/USER-MOFFF/pair_buck6d_coul_gauss_long.cpp
View File

@ -144,17 +144,17 @@ void PairBuck6dCoulGaussLong::compute(int eflag, int vflag)
term3 = term2*term2;
term4 = 1.0/(1.0 + term2);
term5 = 1.0/(1.0 + 2.0*term2 + term3);
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
// smoothing term
if (rsq > rsmooth_sq[itype][jtype]) {
rcu = r*rsq;
rqu = rsq*rsq;
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
// forcebuck6d is -dE/dr*r
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r;
@ -172,15 +172,15 @@ void PairBuck6dCoulGaussLong::compute(int eflag, int vflag)
arg = alpha_ij[itype][jtype]*r;
expa = MathSpecial::expmsq(arg);
erfa = 1 - (MathSpecial::my_erfcx(arg) * expa);
prefactor = qqrd2e*qtmp*q[j]/r;
falpha = erfa - EWALD_F*arg*expa;
forcecoul = prefactor * (falpha - erf + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*falpha;
// (q*q/r) * (gauss(alpha_ij) - gauss(alpha_long)
ealpha = prefactor * (erfa-erf);
// smoothing term - NOTE: ingnored in special_bonds correction
// smoothing term - NOTE: ingnored in special_bonds correction
// since likely rsmooth_sq_c >> d(special)
if (rsq > rsmooth_sq_c) {
rcu = r*rsq;
@ -191,7 +191,7 @@ void PairBuck6dCoulGaussLong::compute(int eflag, int vflag)
ealpha *= sme;
}
} else forcecoul = 0.0;
fpair = (forcecoul + factor_lj*forcebuck6d) * r2inv;
f[i][0] += delx*fpair;
@ -340,8 +340,8 @@ void PairBuck6dCoulGaussLong::init_style()
neighbor->request(this,instance_me);
cut_coulsq = cut_coul * cut_coul;
cut_coulsq = cut_coul * cut_coul;
//calculation of smoothing coefficients c0_c-c5_c for coulomb smoothing
c0_c = c1_c = c2_c = c3_c = c4_c = c5_c = 0.0;
rsmooth_sq_c = cut_coulsq;
@ -370,7 +370,7 @@ double PairBuck6dCoulGaussLong::init_one(int i, int j)
double cut = MAX(cut_lj[i][j],cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
//calculation of smoothing coefficients c0-c5
c0[i][j] = c1[i][j] = c2[i][j] = c3[i][j] = c4[i][j] = c5[i][j] = 0.0;
rsmooth_sq[i][j] = cut_ljsq[i][j];
@ -395,7 +395,7 @@ double PairBuck6dCoulGaussLong::init_one(int i, int j)
double rexp = exp(-cut_lj[i][j]*buck6d2[i][j]);
offset[i][j] = buck6d1[i][j]*rexp - term1*term4;
} else offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[i][j];
alpha_ij[j][i] = alpha_ij[i][j];
buck6d1[j][i] = buck6d1[i][j];
@ -557,18 +557,18 @@ double PairBuck6dCoulGaussLong::single(int i, int j, int itype, int jtype, doubl
term3 = term2*term2;
term4 = 1.0/(1.0 + term2);
term5 = 1.0/(1.0 + 2.0*term2 + term3);
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
// smoothing term
if (rsq > rsmooth_sq[itype][jtype]) {
rcu = r*rsq;
rqu = rsq*rsq;
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
// forcebuck6d is -dE/dr*r
// forcebuck6d is -dE/dr*r
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
ebuck6d *= sme;
}
@ -584,12 +584,12 @@ double PairBuck6dCoulGaussLong::single(int i, int j, int itype, int jtype, doubl
arg = alpha_ij[itype][jtype]*r;
expa = MathSpecial::expmsq(arg);
erfa = 1 - (MathSpecial::my_erfcx(arg) * expa);
prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
falpha = erfa - EWALD_F*arg*expa;
forcecoul = prefactor * (falpha - erf + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*falpha;
ealpha = prefactor * (erfa-erf);
// smoothing term
if (rsq > rsmooth_sq_c) {

14
src/USER-PTM/ptm_polar.cpp
View File

@ -1,5 +1,5 @@
/*******************************************************************************
* -/_|:|_|_\-
* -/_|:|_|_\-
*
* This code is a modification of D.L. Theobald's QCP rotation code.
* It has been adapted to calculate the polar decomposition of a 3x3 matrix
@ -14,7 +14,7 @@
* USA
*
* dtheobald@brandeis.edu
*
*
* Pu Liu
* Johnson & Johnson Pharmaceutical Research and Development, L.L.C.
* 665 Stockton Drive
@ -22,7 +22,7 @@
* USA
*
* pliu24@its.jnj.com
*
*
*
* If you use this QCP rotation calculation method in a publication, please
* reference:
@ -33,7 +33,7 @@
* Acta Crystallographica A 61(4):478-480.
*
* Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald (2009)
* "Fast determination of the optimal rotational matrix for macromolecular
* "Fast determination of the optimal rotational matrix for macromolecular
* superpositions."
* Journal of Computational Chemistry 31(7):1561-1563.
*
@ -63,7 +63,7 @@
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* Source: started anew.
*
@ -82,7 +82,7 @@
*
* 2016/05/29 QCP method adapted for polar decomposition of a 3x3 matrix,
* for use in Polyhedral Template Matching.
*
*
******************************************************************************/
#include <cmath>
@ -327,7 +327,7 @@ void InnerProduct(double *A, int num, const double (*coords1)[3], double (*coord
A[6] += z1 * x2;
A[7] += z1 * y2;
A[8] += z1 * z2;
A[8] += z1 * z2;
}
}

90
src/USER-SCAFACOS/scafacos.cpp
View File

@ -65,10 +65,10 @@ void Scafacos::settings(int narg, char **arg)
if (strcmp(method,"fmm") == 0) {
tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
fmm_tuning_flag = 0;
} else if (strcmp(method,"p3m") == 0 ||
strcmp(method,"p2nfft") == 0 ||
} else if (strcmp(method,"p3m") == 0 ||
strcmp(method,"p2nfft") == 0 ||
strcmp(method,"ewald") == 0) {
tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
} else if (strcmp(method,"direct") == 0) {
; // direct summation has no tolerance type
} else {
@ -100,7 +100,7 @@ void Scafacos::init()
if (logfile && me == 0) fprintf(logfile,
"Setting up ScaFaCoS with solver %s ...\n",method);
if (!atom->q_flag)
if (!atom->q_flag)
error->all(FLERR,"Kspace style requires atom attribute q");
if (domain->dimension == 2)
@ -112,7 +112,7 @@ void Scafacos::init()
if (atom->natoms > INT_MAX && sizeof(int) != 8)
error->all(FLERR,"Scafacos atom count exceeds 2B");
if (atom->molecular > 0)
if (atom->molecular > 0)
error->all(FLERR,
"Cannot use Scafacos with molecular charged systems yet");
@ -140,7 +140,7 @@ void Scafacos::init()
}
double **x = atom->x;
double *q = atom->q;
double *q = atom->q;
int nlocal = atom->nlocal;
if (strcmp(method,"fmm") == 0)
@ -164,7 +164,7 @@ void Scafacos::init()
result = fcs_tune((FCS)fcs,nlocal,&x[0][0],q);
check_result((void*)&result);
// more useful here, since the parameters should be tuned now
// more useful here, since the parameters should be tuned now
if (me == 0) fcs_print_parameters((FCS)fcs);
}
@ -194,12 +194,12 @@ void Scafacos::compute(int eflag, int vflag)
}
if (eflag || vflag) ev_setup(eflag,vflag);
else
else
{
eflag_atom = 0;
vflag_global = 0;
}
// grow xpbc, epot, efield if necessary
if (nlocal > maxatom || maxatom == 0) {
@ -270,7 +270,7 @@ void Scafacos::compute(int eflag, int vflag)
f[i][1] += qone * efield[i][1];
f[i][2] += qone * efield[i][2];
myeng += 0.5 * qone * epot[i];
}
}
if (eflag_atom) {
for (int i = 0; i < nlocal; i++)
@ -296,50 +296,50 @@ int Scafacos::modify_param(int narg, char **arg)
if (narg < 3) error->all(FLERR,
"Illegal kspace_modify command (tolerance)");
if (strcmp(arg[2],"energy") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
else if (strcmp(arg[2],"energy_rel") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_ENERGY_REL;
tolerance_type = FCS_TOLERANCE_TYPE_ENERGY_REL;
else if (strcmp(arg[2],"field") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
else if (strcmp(arg[2],"field_rel") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_FIELD_REL;
tolerance_type = FCS_TOLERANCE_TYPE_FIELD_REL;
else if (strcmp(arg[2],"potential") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL;
tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL;
else if (strcmp(arg[2],"potential_rel") == 0)
tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL_REL;
tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL_REL;
else error->all(FLERR,
"Illegal kspace_modify command (tolerance argument)");
// check if method is compatatible to chosen tolerance type
if(
(
strcmp(method,"fmm") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_ENERGY &&
strcmp(method,"fmm") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_ENERGY &&
tolerance_type != FCS_TOLERANCE_TYPE_ENERGY_REL
)
) ||
(
strcmp(method,"p2nfft") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD &&
tolerance_type != FCS_TOLERANCE_TYPE_POTENTIAL
)
) ||
(
strcmp(method,"p3m") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD
)
) ||
(
strcmp(method,"ewald") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD
)
)
) ||
(
strcmp(method,"p2nfft") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD &&
tolerance_type != FCS_TOLERANCE_TYPE_POTENTIAL
)
) ||
(
strcmp(method,"p3m") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD
)
) ||
(
strcmp(method,"ewald") == 0 &&
(
tolerance_type != FCS_TOLERANCE_TYPE_FIELD
)
)
)
error->all(FLERR,"Illegal kspace_modify command \
(invalid tolerance / method combination)");
(invalid tolerance / method combination)");
return 3;
}
@ -375,7 +375,7 @@ double Scafacos::memory_usage()
}
/* ----------------------------------------------------------------------
setup of ScaFaCoS handle with common parameters
setup of ScaFaCoS handle with common parameters
------------------------------------------------------------------------- */
void Scafacos::setup_handle()
@ -443,7 +443,7 @@ bool Scafacos::box_has_changed()
int *periodicity = domain->periodicity;
double *prd = domain->prd;
bool changed =
bool changed =
(periodicity[0] != old_periodicity[0]) ||
(periodicity[1] != old_periodicity[1]) ||
(periodicity[2] != old_periodicity[2]) ||
@ -462,14 +462,14 @@ bool Scafacos::box_has_changed()
check ScaFaCoS result for error condition
------------------------------------------------------------------------- */
void Scafacos::check_result(void* result_p)
void Scafacos::check_result(void* result_p)
{
FCSResult result = *(FCSResult*)result_p;
if (!result) return;
std::stringstream ss;
ss << "ScaFaCoS: " << fcs_result_get_function(result) << "\n"
ss << "ScaFaCoS: " << fcs_result_get_function(result) << "\n"
<< fcs_result_get_message(result) << "\n";
fcs_result_destroy(result);
std::string err_msg = ss.str();
@ -477,4 +477,4 @@ void Scafacos::check_result(void* result_p)
error->one(FLERR,str);
}

2
src/USER-SDPD/fix_meso_move.cpp
View File

@ -360,7 +360,7 @@ void FixMesoMove::init () {
}
void FixMesoMove::setup_pre_force (int /*vflag*/) {
// set vest equal to v
// set vest equal to v
double **v = atom->v;
double **vest = atom->vest;
int *mask = atom->mask;

18
src/USER-SDPD/fix_rigid_meso.cpp
View File

@ -20,7 +20,7 @@
/* ----------------------------------------------------------------------
Contributing author:
Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com
This is an extension of fix/rigid/nve to SPH/SDPD particles
You can see the original copyright notice of fix/rigid authors above
Note that the Kamberaj paper was related to the nvt variant
@ -53,10 +53,10 @@ FixRigid (lmp, narg, arg) {
if ((atom->e_flag != 1) || (atom->rho_flag != 1))
error->all (FLERR, "fix rigid/meso command requires atom_style with"
" both energy and density");
if (langflag || tstat_flag)
error->all (FLERR,"Can not use thermostat with fix rigid/meso");
if (pstat_flag)
error->all (FLERR,"Can not use barostat with fix rigid/meso");
@ -277,11 +277,11 @@ void FixRigidMeso::set_xv () {
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
// half-step update of particle internal energy and density
e[i] += dtf * de[i];
rho[i] += dtf * drho[i];
ibody = body[i];
xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
@ -301,7 +301,7 @@ void FixRigidMeso::set_xv () {
x2 = x[i][2] + zbox*zprd;
}
}
v0 = v[i][0];
v1 = v[i][1];
v2 = v[i][2];
@ -319,7 +319,7 @@ void FixRigidMeso::set_xv () {
vcm[ibody][1];
v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] +
vcm[ibody][2];
vest[i][0] = 2*v[i][0] - v0;
vest[i][1] = 2*v[i][1] - v1;
vest[i][2] = 2*v[i][2] - v2;
@ -406,11 +406,11 @@ void FixRigidMeso::set_v () {
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
// half-step update of particle internal energy and density
e[i] += dtf * de[i];
rho[i] += dtf * drho[i];
const int ibody = body[i];
MathExtra::matvec (ex_space[ibody],ey_space[ibody],

24
src/USER-SDPD/pair_sdpd_taitwater_isothermal.cpp
View File

@ -14,7 +14,7 @@
/* ----------------------------------------------------------------------
Contributing author:
Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com
references: Espanol and Revenga, Phys Rev E 67, 026705 (2003)
------------------------------------------------------------------------- */
@ -169,7 +169,7 @@ void PairSDPDTaitwaterIsothermal::compute (int eflag, int vflag) {
wiener[0][0] = gaussian (generator);
wiener[1][1] = gaussian (generator);
wiener[2][2] = gaussian (generator);
wiener[0][1] = wiener[1][0] = sqrt_2_inv * gaussian (generator);
wiener[0][2] = wiener[2][0] = sqrt_2_inv * gaussian (generator);
wiener[1][2] = wiener[2][1] = sqrt_2_inv * gaussian (generator);
@ -177,18 +177,18 @@ void PairSDPDTaitwaterIsothermal::compute (int eflag, int vflag) {
wiener[0][0] = random->gaussian ();
wiener[1][1] = random->gaussian ();
wiener[2][2] = random->gaussian ();
wiener[0][1] = wiener[1][0] = sqrt_2_inv * random->gaussian ();
wiener[0][2] = wiener[2][0] = sqrt_2_inv * random->gaussian ();
wiener[1][2] = wiener[2][1] = sqrt_2_inv * random->gaussian ();
#endif
prefactor = sqrt (-4. * kBoltzmann*temperature * fvisc * dtinv) / r;
f_random[0] = prefactor * (wiener[0][0]*delx + wiener[0][1]*dely + wiener[0][2]*delz);
f_random[1] = prefactor * (wiener[1][0]*delx + wiener[1][1]*dely + wiener[1][2]*delz);
f_random[2] = prefactor * (wiener[2][0]*delx + wiener[2][1]*dely + wiener[2][2]*delz);
f[i][0] += delx * fpair + (velx + delx * delVdotDelR / rsq) * fvisc + f_random[0];
f[i][1] += dely * fpair + (vely + dely * delVdotDelR / rsq) * fvisc + f_random[1];
f[i][2] += delz * fpair + (velz + delz * delVdotDelR / rsq) * fvisc + f_random[2];
@ -241,13 +241,13 @@ void PairSDPDTaitwaterIsothermal::settings (int narg, char **arg) {
if (narg != 2 && narg != 3)
error->all (FLERR, "Illegal number of arguments for "
"pair_style sdpd/taitwater/morris/isothermal");
temperature = force->numeric (FLERR, arg[0]);
viscosity = force->numeric (FLERR, arg[1]);
if (temperature <= 0) error->all (FLERR, "Temperature must be positive");
if (viscosity <= 0) error->all (FLERR, "Viscosity must be positive");
// seed is immune to underflow/overflow because it is unsigned
seed = comm->nprocs + comm->me + atom->nlocal;
if (narg == 3) seed += force->inumeric (FLERR, arg[2]);
@ -266,7 +266,7 @@ void PairSDPDTaitwaterIsothermal::coeff (int narg, char **arg) {
if (narg != 5)
error->all (FLERR, "Incorrect args for pair_style "
"sph/taitwater/morris coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
@ -277,7 +277,7 @@ void PairSDPDTaitwaterIsothermal::coeff (int narg, char **arg) {
double soundspeed_one = force->numeric (FLERR,arg[3]);
double cut_one = force->numeric (FLERR,arg[4]);
double B_one = soundspeed_one * soundspeed_one * rho0_one / 7.0;
if (rho0_one <= 0) error->all (FLERR, "Density must be positive");
if (soundspeed_one <= 0) error->all (FLERR, "Sound speed must be positive");
if (cut_one <= 0) error->all (FLERR, "Cutoff must be positive");
@ -304,7 +304,7 @@ void PairSDPDTaitwaterIsothermal::coeff (int narg, char **arg) {
------------------------------------------------------------------------- */
double PairSDPDTaitwaterIsothermal::init_one (int i, int j) {
if (setflag[i][j] == 0)
if (setflag[i][j] == 0)
error->all(FLERR,"Not all pair sph/taitwater/morris coeffs are set");
cut[j][i] = cut[i][j];

14
src/USER-UEF/uef_utils.cpp
View File

@ -150,7 +150,7 @@ void UEFBox::step_deform(const double ex, const double ey)
theta[0] +=winv[0][0]*ex + winv[0][1]*ey;
theta[1] +=winv[1][0]*ex + winv[1][1]*ey;
// deformation of the box. reduce() needs to be called regularly or
// deformation of the box. reduce() needs to be called regularly or
// calculation will become unstable
double eps[3];
@ -169,7 +169,7 @@ void UEFBox::step_deform(const double ex, const double ey)
------------------------------------------------------------------------- */
bool UEFBox::reduce()
{
// determine how many times to apply the automorphisms and find new theta
// determine how many times to apply the automorphisms and find new theta
// values
int f1 = round(theta[0]);
@ -184,7 +184,7 @@ bool UEFBox::reduce()
for (int j=0;j<3;j++)
r0[k][j]=r[k][j];
// this modifies the old change basis matrix to handle the case where the
// this modifies the old change basis matrix to handle the case where the
// automorphism transforms the box but the reduced basis doesn't change
// (r0 should still equal r at the end)
@ -417,7 +417,7 @@ void make_unique(double b[3][3], int r[3][3], int ri[3][3])
if (fabs(b[0][0]) < fabs(b[0][1])) {
col_swap(b,0,1);
col_swap(r,0,1);
col_swap(ri,0,1);
col_swap(ri,0,1);
}
if (fabs(b[0][0]) < fabs(b[0][2])) {
col_swap(b,0,2);
@ -433,7 +433,7 @@ void make_unique(double b[3][3], int r[3][3], int ri[3][3])
if (b[0][0] < 0) {
neg_col(b,0);
neg_col(r,0);
neg_col(ri,0);
neg_col(ri,0);
}
if (b[1][1] < 0) {
neg_col(b,1);
@ -442,8 +442,8 @@ void make_unique(double b[3][3], int r[3][3], int ri[3][3])
}
if (det(b) < 0) {
neg_col(b,2);
neg_col(r,2);
neg_col(ri,2);
neg_col(r,2);
neg_col(ri,2);
}
}
}}