From 4e6dffc7cd1c1851a8223a6520709df94eed382e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Apr 2023 11:45:40 -0400 Subject: [PATCH 01/20] silence compiler warning, reformat beginning of file. --- src/RIGID/fix_shake.cpp | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index b9156a0747..6cd624c755 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -31,15 +31,15 @@ #include "respa.h" #include "update.h" -#include #include +#include #include using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -#define RVOUS 1 // 0 for irregular, 1 for all2all +#define RVOUS 1 // 0 for irregular, 1 for all2all static constexpr double BIG = 1.0e20; static constexpr double MASSDELTA = 0.1; @@ -47,15 +47,15 @@ static constexpr double MASSDELTA = 0.1; /* ---------------------------------------------------------------------- */ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), bond_flag(nullptr), angle_flag(nullptr), - type_flag(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), - loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), - vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), - shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), - b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr), - b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), - a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr), - a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr), closest_list(nullptr) + Fix(lmp, narg, arg), bond_flag(nullptr), angle_flag(nullptr), type_flag(nullptr), + mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), loop_respa(nullptr), + step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr), + mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), shake_atom(nullptr), + shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), closest_list(nullptr), + b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr), + b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), + a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr), + a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr) { energy_global_flag = energy_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1; @@ -76,7 +76,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : molecular = atom->molecular; if (molecular == Atom::ATOMIC) - error->all(FLERR,"Cannot use fix {} with non-molecular system", style); + error->all(FLERR, "Cannot use fix {} with non-molecular system", style); // perform initial allocation of atom-based arrays // register with Atom class @@ -97,13 +97,13 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : comm_forward = 3; // parse SHAKE args - auto mystyle = fmt::format("fix {}",style); + auto mystyle = fmt::format("fix {}", style); - if (narg < 8) utils::missing_cmd_args(FLERR,mystyle, error); + if (narg < 8) utils::missing_cmd_args(FLERR, mystyle, error); - tolerance = utils::numeric(FLERR,arg[3],false,lmp); - max_iter = utils::inumeric(FLERR,arg[4],false,lmp); - output_every = utils::inumeric(FLERR,arg[5],false,lmp); + tolerance = utils::numeric(FLERR, arg[3], false, lmp); + max_iter = utils::inumeric(FLERR, arg[4], false, lmp); + output_every = utils::inumeric(FLERR, arg[5], false, lmp); // parse SHAKE args for bond and angle types // will be used by find_clusters @@ -111,11 +111,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : // store args for "m" in list of length nmass for looping over // for "m" verify that atom masses have been set - bond_flag = new int[atom->nbondtypes+1]; + bond_flag = new int[atom->nbondtypes + 1]; for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0; - angle_flag = new int[atom->nangletypes+1]; + angle_flag = new int[atom->nangletypes + 1]; for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0; - type_flag = new int[atom->ntypes+1]; + type_flag = new int[atom->ntypes + 1]; for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0; mass_list = new double[atom->ntypes]; nmass = 0; From 27127a46ccd4e524a1e03a052f39badbd491e8e0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Apr 2023 21:41:09 -0400 Subject: [PATCH 02/20] enable and apply clang-format --- src/region_block.cpp | 357 +++++++++++++++++---------------- src/region_cone.cpp | 335 ++++++++++++++++--------------- src/region_cylinder.cpp | 428 ++++++++++++++++++++++------------------ src/region_prism.cpp | 339 +++++++++++++++++-------------- src/region_sphere.cpp | 119 ++++++----- 5 files changed, 844 insertions(+), 734 deletions(-) diff --git a/src/region_block.cpp b/src/region_block.cpp index c20d3abb17..8b50521d7c 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -24,100 +23,125 @@ using namespace LAMMPS_NS; -enum{CONSTANT,VARIABLE}; +enum { CONSTANT, VARIABLE }; -#define BIG 1.0e20 +static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), zlostr(nullptr), zhistr(nullptr) + Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), + zlostr(nullptr), zhistr(nullptr) { - options(narg-8,&arg[8]); + options(narg - 8, &arg[8]); xlostyle = CONSTANT; - if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { + if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[2],"INF") == 0) xlo = -BIG; - else if (domain->triclinic == 0) xlo = domain->boxlo[0]; - else xlo = domain->boxlo_bound[0]; - } else if (utils::strmatch(arg[2],"^v_")) { - xlostr = utils::strdup(arg[2]+2); - xlo = 0.0; - xlostyle = VARIABLE; - varshape = 1; - } else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[2], "INF") == 0) + xlo = -BIG; + else if (domain->triclinic == 0) + xlo = domain->boxlo[0]; + else + xlo = domain->boxlo_bound[0]; + } else if (utils::strmatch(arg[2], "^v_")) { + xlostr = utils::strdup(arg[2] + 2); + xlo = 0.0; + xlostyle = VARIABLE; + varshape = 1; + } else + xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp); xhistyle = CONSTANT; - if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { + if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[3],"INF") == 0) xhi = BIG; - else if (domain->triclinic == 0) xhi = domain->boxhi[0]; - else xhi = domain->boxhi_bound[0]; - } else if (utils::strmatch(arg[3],"^v_")) { - xhistr = utils::strdup(arg[3]+2); - xhi = 0.0; - xhistyle = VARIABLE; - varshape = 1; - } else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[3], "INF") == 0) + xhi = BIG; + else if (domain->triclinic == 0) + xhi = domain->boxhi[0]; + else + xhi = domain->boxhi_bound[0]; + } else if (utils::strmatch(arg[3], "^v_")) { + xhistr = utils::strdup(arg[3] + 2); + xhi = 0.0; + xhistyle = VARIABLE; + varshape = 1; + } else + xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp); ylostyle = CONSTANT; - if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { + if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[4],"INF") == 0) ylo = -BIG; - else if (domain->triclinic == 0) ylo = domain->boxlo[1]; - else ylo = domain->boxlo_bound[1]; - } else if (utils::strmatch(arg[4],"^v_")) { - ylostr = utils::strdup(arg[4]+2); - ylo = 0.0; - ylostyle = VARIABLE; - varshape = 1; - } else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[4], "INF") == 0) + ylo = -BIG; + else if (domain->triclinic == 0) + ylo = domain->boxlo[1]; + else + ylo = domain->boxlo_bound[1]; + } else if (utils::strmatch(arg[4], "^v_")) { + ylostr = utils::strdup(arg[4] + 2); + ylo = 0.0; + ylostyle = VARIABLE; + varshape = 1; + } else + ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp); yhistyle = CONSTANT; - if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { + if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[5],"INF") == 0) yhi = BIG; - else if (domain->triclinic == 0) yhi = domain->boxhi[1]; - else yhi = domain->boxhi_bound[1]; - } else if (utils::strmatch(arg[5],"^v_")) { - yhistr = utils::strdup(arg[5]+2); - yhi = 0.0; - yhistyle = VARIABLE; - varshape = 1; - } else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[5], "INF") == 0) + yhi = BIG; + else if (domain->triclinic == 0) + yhi = domain->boxhi[1]; + else + yhi = domain->boxhi_bound[1]; + } else if (utils::strmatch(arg[5], "^v_")) { + yhistr = utils::strdup(arg[5] + 2); + yhi = 0.0; + yhistyle = VARIABLE; + varshape = 1; + } else + yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp); zlostyle = CONSTANT; - if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { + if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[6],"INF") == 0) zlo = -BIG; - else if (domain->triclinic == 0) zlo = domain->boxlo[2]; - else zlo = domain->boxlo_bound[2]; - } else if (utils::strmatch(arg[6],"^v_")) { - zlostr = utils::strdup(arg[6]+2); - zlo = 0.0; - zlostyle = VARIABLE; - varshape = 1; - } else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[6], "INF") == 0) + zlo = -BIG; + else if (domain->triclinic == 0) + zlo = domain->boxlo[2]; + else + zlo = domain->boxlo_bound[2]; + } else if (utils::strmatch(arg[6], "^v_")) { + zlostr = utils::strdup(arg[6] + 2); + zlo = 0.0; + zlostyle = VARIABLE; + varshape = 1; + } else + zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp); zhistyle = CONSTANT; - if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { + if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[7],"INF") == 0) zhi = BIG; - else if (domain->triclinic == 0) zhi = domain->boxhi[2]; - else zhi = domain->boxhi_bound[2]; - } else if (utils::strmatch(arg[7],"^v_")) { - zhistr = utils::strdup(arg[7]+2); - zhi = 0.0; - zhistyle = VARIABLE; - varshape = 1; - } else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[7], "INF") == 0) + zhi = BIG; + else if (domain->triclinic == 0) + zhi = domain->boxhi[2]; + else + zhi = domain->boxhi_bound[2]; + } else if (utils::strmatch(arg[7], "^v_")) { + zhistr = utils::strdup(arg[7] + 2); + zhi = 0.0; + zhistyle = VARIABLE; + varshape = 1; + } else + zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp); if (varshape) { variable_check(); @@ -126,9 +150,9 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : // error check - if (xlo > xhi) error->all(FLERR,"Illegal region block xlo: {} >= xhi: {}", xlo, xhi); - if (ylo > yhi) error->all(FLERR,"Illegal region block ylo: {} >= yhi: {}", ylo, yhi); - if (zlo > zhi) error->all(FLERR,"Illegal region block zlo: {} >= zhi: {}", zlo, zhi); + if (xlo > xhi) error->all(FLERR, "Illegal region block xlo: {} >= xhi: {}", xlo, xhi); + if (ylo > yhi) error->all(FLERR, "Illegal region block ylo: {} >= yhi: {}", ylo, yhi); + if (zlo > zhi) error->all(FLERR, "Illegal region block zlo: {} >= zhi: {}", zlo, zhi); // extent of block @@ -140,15 +164,18 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : extent_yhi = yhi; extent_zlo = zlo; extent_zhi = zhi; - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to all 6 planes // particle can only touch 3 planes cmax = 6; contact = new Contact[cmax]; - if (interior) tmax = 3; - else tmax = 1; + if (interior) + tmax = 3; + else + tmax = 1; // open face data structs @@ -235,13 +262,13 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : RegBlock::~RegBlock() { if (copymode) return; - delete [] xlostr; - delete [] xhistr; - delete [] ylostr; - delete [] yhistr; - delete [] zlostr; - delete [] zhistr; - delete [] contact; + delete[] xlostr; + delete[] xhistr; + delete[] ylostr; + delete[] yhistr; + delete[] zlostr; + delete[] zhistr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -259,8 +286,7 @@ void RegBlock::init() int RegBlock::inside(double x, double y, double z) { - if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) - return 1; + if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) return 1; return 0; } @@ -277,8 +303,7 @@ int RegBlock::surface_interior(double *x, double cutoff) // x is exterior to block - if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || - x[2] < zlo || x[2] > zhi) return 0; + if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0; // x is interior to block or on its surface @@ -352,17 +377,16 @@ int RegBlock::surface_interior(double *x, double cutoff) int RegBlock::surface_exterior(double *x, double cutoff) { - double xp,yp,zp; - double xc,yc,zc,dist,mindist; + double xp, yp, zp; + double xc, yc, zc, dist, mindist; // x is far enough from block that there is no contact // x is interior to block - if (x[0] <= xlo-cutoff || x[0] >= xhi+cutoff || - x[1] <= ylo-cutoff || x[1] >= yhi+cutoff || - x[2] <= zlo-cutoff || x[2] >= zhi+cutoff) return 0; - if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && - x[2] > zlo && x[2] < zhi) return 0; + if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff || + x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff) + return 0; + if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0; // x is exterior to block or on its surface // xp,yp,zp = point on surface of block that x is closest to @@ -370,20 +394,29 @@ int RegBlock::surface_exterior(double *x, double cutoff) // do not add contact point if r >= cutoff if (!openflag) { - if (x[0] < xlo) xp = xlo; - else if (x[0] > xhi) xp = xhi; - else xp = x[0]; - if (x[1] < ylo) yp = ylo; - else if (x[1] > yhi) yp = yhi; - else yp = x[1]; - if (x[2] < zlo) zp = zlo; - else if (x[2] > zhi) zp = zhi; - else zp = x[2]; + if (x[0] < xlo) + xp = xlo; + else if (x[0] > xhi) + xp = xhi; + else + xp = x[0]; + if (x[1] < ylo) + yp = ylo; + else if (x[1] > yhi) + yp = yhi; + else + yp = x[1]; + if (x[2] < zlo) + zp = zlo; + else if (x[2] > zhi) + zp = zhi; + else + zp = x[2]; } else { mindist = BIG; for (int i = 0; i < 6; i++) { if (open_faces[i]) continue; - dist = find_closest_point(i,x,xc,yc,zc); + dist = find_closest_point(i, x, xc, yc, zc); if (dist < mindist) { xp = xc; yp = yc; @@ -393,7 +426,7 @@ int RegBlock::surface_exterior(double *x, double cutoff) } } - add_contact(0,x,xp,yp,zp); + add_contact(0, x, xp, yp, zp); contact[0].iwall = 0; if (contact[0].r < cutoff) return 1; return 0; @@ -403,28 +436,17 @@ int RegBlock::surface_exterior(double *x, double cutoff) change region shape via variable evaluation ------------------------------------------------------------------------- */ -void RegBlock::shape_update() // addition +void RegBlock::shape_update() // addition { - if (xlostyle == VARIABLE) - xlo = xscale * input->variable->compute_equal(xlovar); - - if (xhistyle == VARIABLE) - xhi = xscale * input->variable->compute_equal(xhivar); - - if (ylostyle == VARIABLE) - ylo = yscale * input->variable->compute_equal(ylovar); - - if (yhistyle == VARIABLE) - yhi = yscale * input->variable->compute_equal(yhivar); - - if (zlostyle == VARIABLE) - zlo = zscale * input->variable->compute_equal(zlovar); - - if (zhistyle == VARIABLE) - zhi = zscale * input->variable->compute_equal(zhivar); + if (xlostyle == VARIABLE) xlo = xscale * input->variable->compute_equal(xlovar); + if (xhistyle == VARIABLE) xhi = xscale * input->variable->compute_equal(xhivar); + if (ylostyle == VARIABLE) ylo = yscale * input->variable->compute_equal(ylovar); + if (yhistyle == VARIABLE) yhi = yscale * input->variable->compute_equal(yhivar); + if (zlostyle == VARIABLE) zlo = zscale * input->variable->compute_equal(zlovar); + if (zhistyle == VARIABLE) zhi = zscale * input->variable->compute_equal(zhivar); if (xlo > xhi || ylo > yhi || zlo > zhi) - error->one(FLERR,"Variable evaluation in region gave bad value"); + error->one(FLERR, "Variable evaluation in region gave bad value"); // face[0] @@ -489,54 +511,48 @@ void RegBlock::shape_update() // addition error check on existence of variable ------------------------------------------------------------------------- */ -void RegBlock::variable_check() // addition +void RegBlock::variable_check() // addition { if (xlostyle == VARIABLE) { xlovar = input->variable->find(xlostr); - if (xlovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (xlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(xlovar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } if (xhistyle == VARIABLE) { xhivar = input->variable->find(xhistr); - if (xhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (xhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(xhivar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } if (ylostyle == VARIABLE) { ylovar = input->variable->find(ylostr); - if (ylovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (ylovar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(ylovar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } if (yhistyle == VARIABLE) { yhivar = input->variable->find(yhistr); - if (yhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (yhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(yhivar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } if (zlostyle == VARIABLE) { zlovar = input->variable->find(zlostr); - if (zlovar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (zlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(zlovar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } if (zhistyle == VARIABLE) { zhivar = input->variable->find(zhistr); - if (zhivar < 0) - error->all(FLERR,"Variable name for region block does not exist"); + if (zhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); if (!input->variable->equalstyle(zhivar)) - error->all(FLERR,"Variable for region block is invalid style"); + error->all(FLERR, "Variable for region block is invalid style"); } } @@ -545,36 +561,35 @@ void RegBlock::variable_check() // addition store closest point in xc,yc,zc --------------------------------------------------------------------------*/ -double RegBlock::find_closest_point(int i, double *x, - double &xc, double &yc, double &zc) +double RegBlock::find_closest_point(int i, double *x, double &xc, double &yc, double &zc) { - double dot,d2,d2min; - double xr[3],xproj[3],p[3]; + double dot, d2, d2min; + double xr[3], xproj[3], p[3]; xr[0] = x[0] - corners[i][0][0]; xr[1] = x[1] - corners[i][0][1]; xr[2] = x[2] - corners[i][0][2]; - dot = face[i][0]*xr[0] + face[i][1]*xr[1] + face[i][2]*xr[2]; - xproj[0] = xr[0] - dot*face[i][0]; - xproj[1] = xr[1] - dot*face[i][1]; - xproj[2] = xr[2] - dot*face[i][2]; + dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2]; + xproj[0] = xr[0] - dot * face[i][0]; + xproj[1] = xr[1] - dot * face[i][1]; + xproj[2] = xr[2] - dot * face[i][2]; d2min = BIG; // check if point projects inside of face if (inside_face(xproj, i)) { - d2 = d2min = dot*dot; + d2 = d2min = dot * dot; xc = xproj[0] + corners[i][0][0]; yc = xproj[1] + corners[i][0][1]; zc = xproj[2] + corners[i][0][2]; - // check each edge + // check each edge } else { - point_on_line_segment(corners[i][0],corners[i][1],x,p); - d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + - (p[2]-x[2])*(p[2]-x[2]); + point_on_line_segment(corners[i][0], corners[i][1], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -582,9 +597,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - point_on_line_segment(corners[i][1],corners[i][2],x,p); - d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + - (p[2]-x[2])*(p[2]-x[2]); + point_on_line_segment(corners[i][1], corners[i][2], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -592,9 +607,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - point_on_line_segment(corners[i][2],corners[i][3],x,p); - d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + - (p[2]-x[2])*(p[2]-x[2]); + point_on_line_segment(corners[i][2], corners[i][3], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -602,9 +617,9 @@ double RegBlock::find_closest_point(int i, double *x, zc = p[2]; } - point_on_line_segment(corners[i][3],corners[i][0],x,p); - d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + - (p[2]-x[2])*(p[2]-x[2]); + point_on_line_segment(corners[i][3], corners[i][0], x, p); + d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) + + (p[2] - x[2]) * (p[2] - x[2]); if (d2 < d2min) { d2min = d2; xc = p[0]; @@ -623,14 +638,12 @@ double RegBlock::find_closest_point(int i, double *x, int RegBlock::inside_face(double *xproj, int iface) { if (iface < 2) { - if (xproj[1] > 0 && (xproj[1] < yhi-ylo) && - xproj[2] > 0 && (xproj[2] < zhi-zlo)) return 1; + if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; } else if (iface < 4) { - if (xproj[0] > 0 && (xproj[0] < (xhi-xlo)) && - xproj[2] > 0 && (xproj[2] < (zhi-zlo))) return 1; + if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo))) + return 1; } else { - if (xproj[0] > 0 && xproj[0] < (xhi-xlo) && - xproj[1] > 0 && xproj[1] < (yhi-ylo)) return 1; + if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1; } return 0; diff --git a/src/region_cone.cpp b/src/region_cone.cpp index e71a6a72a3..76bedf4c92 100644 --- a/src/region_cone.cpp +++ b/src/region_cone.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -26,102 +25,120 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 +static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ -RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), lo(0.0), hi(0.0) +RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), lo(0.0), hi(0.0) { - options(narg-9,&arg[9]); + options(narg - 9, &arg[9]); // check open face settings if (openflag) - for (int i=3; i<6; i++) - if (open_faces[i]) error->all(FLERR,"Illegal region cone open face: {}", i+1); + for (int i = 3; i < 6; i++) + if (open_faces[i]) error->all(FLERR, "Illegal region cone open face: {}", i + 1); - if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) - error->all(FLERR,"Illegal region cone axis: {}", arg[2]); + if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0) + error->all(FLERR, "Illegal region cone axis: {}", arg[2]); axis = arg[2][0]; if (axis == 'x') { - c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = yscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = yscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = yscale * utils::numeric(FLERR, arg[5], false, lmp); + radiushi = yscale * utils::numeric(FLERR, arg[6], false, lmp); } else if (axis == 'y') { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp); + radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp); } else if (axis == 'z') { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); - c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); - radiuslo = xscale*utils::numeric(FLERR,arg[5],false,lmp); - radiushi = xscale*utils::numeric(FLERR,arg[6],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); + c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp); + radiuslo = xscale * utils::numeric(FLERR, arg[5], false, lmp); + radiushi = xscale * utils::numeric(FLERR, arg[6], false, lmp); } - if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { + if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); if (axis == 'x') { - if (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[0]; - else lo = domain->boxlo_bound[0]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[0]; + else + lo = domain->boxlo_bound[0]; } if (axis == 'y') { - if (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[1]; - else lo = domain->boxlo_bound[1]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[1]; + else + lo = domain->boxlo_bound[1]; } if (axis == 'z') { - if (strcmp(arg[7],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[2]; - else lo = domain->boxlo_bound[2]; + if (strcmp(arg[7], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[2]; + else + lo = domain->boxlo_bound[2]; } } else { - if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[7],false,lmp); - if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[7],false,lmp); - if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[7],false,lmp); + if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[7], false, lmp); + if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[7], false, lmp); + if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[7], false, lmp); } - if (strcmp(arg[8],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { + if (strcmp(arg[8], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); if (axis == 'x') { - if (strcmp(arg[8],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[0]; - else hi = domain->boxhi_bound[0]; + if (strcmp(arg[8], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[0]; + else + hi = domain->boxhi_bound[0]; } if (axis == 'y') { - if (strcmp(arg[8],"INF") == 0) hi = BIG; - if (domain->triclinic == 0) hi = domain->boxhi[1]; - else hi = domain->boxhi_bound[1]; + if (strcmp(arg[8], "INF") == 0) hi = BIG; + if (domain->triclinic == 0) + hi = domain->boxhi[1]; + else + hi = domain->boxhi_bound[1]; } if (axis == 'z') { - if (strcmp(arg[8],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[2]; - else hi = domain->boxhi_bound[2]; + if (strcmp(arg[8], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[2]; + else + hi = domain->boxhi_bound[2]; } } else { - if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[8],false,lmp); - if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[8],false,lmp); - if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[8],false,lmp); + if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[8], false, lmp); + if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[8], false, lmp); + if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[8], false, lmp); } // error check - if (radiuslo < 0.0) error->all(FLERR,"Illegal radius in region cone command"); - if (radiushi < 0.0) error->all(FLERR,"Illegal radius in region cone command"); + if (radiuslo < 0.0) error->all(FLERR, "Illegal radius in region cone command"); + if (radiushi < 0.0) error->all(FLERR, "Illegal radius in region cone command"); if (radiuslo == 0.0 && radiushi == 0.0) - error->all(FLERR,"Illegal radius in region cone command"); - if (hi == lo) error->all(FLERR,"Illegal cone length in region cone command"); + error->all(FLERR, "Illegal radius in region cone command"); + if (hi == lo) error->all(FLERR, "Illegal cone length in region cone command"); // extent of cone - if (radiuslo > radiushi) maxradius = radiuslo; - else maxradius = radiushi; + if (radiuslo > radiushi) + maxradius = radiuslo; + else + maxradius = radiushi; if (interior) { bboxflag = 1; @@ -149,22 +166,25 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) : extent_zlo = lo; extent_zhi = hi; } - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to cone surface and 2 ends // particle can only touch surface and 1 end cmax = 3; contact = new Contact[cmax]; - if (interior) tmax = 2; - else tmax = 1; + if (interior) + tmax = 2; + else + tmax = 1; } /* ---------------------------------------------------------------------- */ RegCone::~RegCone() { - delete [] contact; + delete[] contact; } /* ---------------------------------------------------------------------- @@ -174,30 +194,36 @@ RegCone::~RegCone() int RegCone::inside(double x, double y, double z) { - double del1,del2,dist; + double del1, del2, dist; double currentradius; if (axis == 'x') { del1 = y - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && x >= lo && x <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && x >= lo && x <= hi) + return 1; + else + return 0; } else if (axis == 'y') { del1 = x - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (y-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && y >= lo && y <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (y - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && y >= lo && y <= hi) + return 1; + else + return 0; } else if (axis == 'z') { del1 = x - c1; del2 = y - c2; - dist = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (z-lo)*(radiushi-radiuslo)/(hi-lo); - if (dist <= currentradius && z >= lo && z <= hi) return 1; - else return 0; + dist = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (z - lo) * (radiushi - radiuslo) / (hi - lo); + if (dist <= currentradius && z >= lo && z <= hi) + return 1; + else + return 0; } return 0; @@ -213,16 +239,16 @@ int RegCone::inside(double x, double y, double z) int RegCone::surface_interior(double *x, double cutoff) { - double del1,del2,r,currentradius,delx,dely,delz,dist,delta; - double surflo[3],surfhi[3],xs[3]; + double del1, del2, r, currentradius, delx, dely, delz, dist, delta; + double surflo[3], surfhi[3], xs[3]; int n = 0; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo); // x is exterior to cone @@ -234,23 +260,22 @@ int RegCone::surface_interior(double *x, double cutoff) if (r > 0.0 && !open_faces[2]) { surflo[0] = lo; - surflo[1] = c1 + del1*radiuslo/r; - surflo[2] = c2 + del2*radiuslo/r; + surflo[1] = c1 + del1 * radiuslo / r; + surflo[2] = c2 + del2 * radiuslo / r; surfhi[0] = hi; - surfhi[1] = c1 + del1*radiushi/r; - surfhi[2] = c2 + del2*radiushi/r; - point_on_line_segment(surflo,surfhi,x,xs); + surfhi[1] = c1 + del1 * radiushi / r; + surfhi[2] = c2 + del2 * radiushi / r; + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; - contact[n].radius = -2.0*(radiuslo + (xs[0]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[0] - lo) * (radiushi - radiuslo) / (hi - lo)); contact[n].iwall = 2; n++; } @@ -278,9 +303,8 @@ int RegCone::surface_interior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[1]-lo)* - (radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo); // y is exterior to cone @@ -291,25 +315,24 @@ int RegCone::surface_interior(double *x, double cutoff) // surfhi = pt on outer circle of top end plane, same dir as y vs axis if (r > 0.0 && !open_faces[2]) { - surflo[0] = c1 + del1*radiuslo/r; + surflo[0] = c1 + del1 * radiuslo / r; surflo[1] = lo; - surflo[2] = c2 + del2*radiuslo/r; - surfhi[0] = c1 + del1*radiushi/r; + surflo[2] = c2 + del2 * radiuslo / r; + surfhi[0] = c1 + del1 * radiushi / r; surfhi[1] = hi; - surfhi[2] = c2 + del2*radiushi/r; - point_on_line_segment(surflo,surfhi,x,xs); + surfhi[2] = c2 + del2 * radiushi / r; + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; contact[n].iwall = 2; - contact[n].radius = -2.0*(radiuslo + (xs[1]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[1] - lo) * (radiushi - radiuslo) / (hi - lo)); n++; } } @@ -336,8 +359,8 @@ int RegCone::surface_interior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo); // z is exterior to cone @@ -348,31 +371,30 @@ int RegCone::surface_interior(double *x, double cutoff) // surfhi = pt on outer circle of top end plane, same dir as z vs axis if (r > 0.0 && !open_faces[2]) { - surflo[0] = c1 + del1*radiuslo/r; - surflo[1] = c2 + del2*radiuslo/r; + surflo[0] = c1 + del1 * radiuslo / r; + surflo[1] = c2 + del2 * radiuslo / r; surflo[2] = lo; - surfhi[0] = c1 + del1*radiushi/r; - surfhi[1] = c2 + del2*radiushi/r; + surfhi[0] = c1 + del1 * radiushi / r; + surfhi[1] = c2 + del2 * radiushi / r; surfhi[2] = hi; - point_on_line_segment(surflo,surfhi,x,xs); + point_on_line_segment(surflo, surfhi, x, xs); delx = x[0] - xs[0]; dely = x[1] - xs[1]; delz = x[2] - xs[2]; - dist = sqrt(delx*delx + dely*dely + delz*delz); + dist = sqrt(delx * delx + dely * dely + delz * delz); if (dist < cutoff) { contact[n].r = dist; contact[n].delx = delx; contact[n].dely = dely; contact[n].delz = delz; contact[n].iwall = 2; - contact[n].radius = -2.0*(radiuslo + (xs[2]-lo)* - (radiushi-radiuslo)/(hi-lo)); + contact[n].radius = -2.0 * (radiuslo + (xs[2] - lo) * (radiushi - radiuslo) / (hi - lo)); n++; } } delta = x[2] - lo; - if (delta < cutoff && !open_faces[0]) { + if (delta < cutoff && !open_faces[0]) { contact[n].r = delta; contact[n].delz = delta; contact[n].delx = contact[n].dely = 0.0; @@ -381,7 +403,7 @@ int RegCone::surface_interior(double *x, double cutoff) n++; } delta = hi - x[2]; - if (delta < cutoff && !open_faces[1]) { + if (delta < cutoff && !open_faces[1]) { contact[n].r = delta; contact[n].delz = -delta; contact[n].delx = contact[n].dely = 0.0; @@ -402,14 +424,14 @@ int RegCone::surface_interior(double *x, double cutoff) int RegCone::surface_exterior(double *x, double cutoff) { - double del1,del2,r,currentradius,distsq,distsqprev,crad; - double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3]; + double del1, del2, r, currentradius, distsq, distsqprev, crad; + double corner1[3], corner2[3], corner3[3], corner4[3], xp[3], nearest[3]; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[0] - lo) * (radiushi - radiuslo) / (hi - lo); // radius of curvature, only used for granular walls @@ -418,8 +440,7 @@ int RegCone::surface_exterior(double *x, double cutoff) // x is far enough from cone that there is no contact // x is interior to cone - if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) - return 0; + if (r >= maxradius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0; if (r < currentradius && x[0] > lo && x[0] < hi) return 0; // x is exterior to cone or on its surface @@ -430,11 +451,11 @@ int RegCone::surface_exterior(double *x, double cutoff) // do not add contact point if r >= cutoff corner1[0] = lo; - corner1[1] = c1 + del1*radiuslo/r; - corner1[2] = c2 + del2*radiuslo/r; + corner1[1] = c1 + del1 * radiuslo / r; + corner1[2] = c2 + del2 * radiuslo / r; corner2[0] = hi; - corner2[1] = c1 + del1*radiushi/r; - corner2[2] = c2 + del2*radiushi/r; + corner2[1] = c1 + del1 * radiushi / r; + corner2[2] = c2 + del2 * radiushi / r; corner3[0] = lo; corner3[1] = c1; corner3[2] = c2; @@ -445,27 +466,27 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[0]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[0] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0.0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0.0; } if (distsq == BIG) return 0; - add_contact(0,x,nearest[0],nearest[1],nearest[2]); + add_contact(0, x, nearest[0], nearest[1], nearest[2]); contact[0].radius = crad; contact[0].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -474,8 +495,8 @@ int RegCone::surface_exterior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[1] - lo) * (radiushi - radiuslo) / (hi - lo); // radius of curvature, only used for granular walls @@ -484,8 +505,7 @@ int RegCone::surface_exterior(double *x, double cutoff) // y is far enough from cone that there is no contact // y is interior to cone - if (r >= maxradius+cutoff || - x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0; + if (r >= maxradius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0; if (r < currentradius && x[1] > lo && x[1] < hi) return 0; // y is exterior to cone or on its surface @@ -495,12 +515,12 @@ int RegCone::surface_exterior(double *x, double cutoff) // could be edge of cone // do not add contact point if r >= cutoff - corner1[0] = c1 + del1*radiuslo/r; + corner1[0] = c1 + del1 * radiuslo / r; corner1[1] = lo; - corner1[2] = c2 + del2*radiuslo/r; - corner2[0] = c1 + del1*radiushi/r; + corner1[2] = c2 + del2 * radiuslo / r; + corner2[0] = c1 + del1 * radiushi / r; corner2[1] = hi; - corner2[2] = c2 + del2*radiushi/r; + corner2[2] = c2 + del2 * radiushi / r; corner3[0] = c1; corner3[1] = lo; corner3[2] = c2; @@ -511,26 +531,26 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[1]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[1] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } - add_contact(0,x,nearest[0],nearest[1],nearest[2]); + add_contact(0, x, nearest[0], nearest[1], nearest[2]); contact[0].radius = crad; contact[0].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -539,8 +559,8 @@ int RegCone::surface_exterior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); - currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo); + r = sqrt(del1 * del1 + del2 * del2); + currentradius = radiuslo + (x[2] - lo) * (radiushi - radiuslo) / (hi - lo); // radius of curvature, only used for granular walls @@ -549,8 +569,7 @@ int RegCone::surface_exterior(double *x, double cutoff) // z is far enough from cone that there is no contact // z is interior to cone - if (r >= maxradius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) - return 0; + if (r >= maxradius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0; if (r < currentradius && x[2] > lo && x[2] < hi) return 0; // z is exterior to cone or on its surface @@ -560,11 +579,11 @@ int RegCone::surface_exterior(double *x, double cutoff) // could be edge of cone // do not add contact point if r >= cutoff - corner1[0] = c1 + del1*radiuslo/r; - corner1[1] = c2 + del2*radiuslo/r; + corner1[0] = c1 + del1 * radiuslo / r; + corner1[1] = c2 + del2 * radiuslo / r; corner1[2] = lo; - corner2[0] = c1 + del1*radiushi/r; - corner2[1] = c2 + del2*radiushi/r; + corner2[0] = c1 + del1 * radiushi / r; + corner2[1] = c2 + del2 * radiushi / r; corner2[2] = hi; corner3[0] = c1; corner3[1] = c2; @@ -576,26 +595,26 @@ int RegCone::surface_exterior(double *x, double cutoff) distsq = BIG; if (!open_faces[2]) { - point_on_line_segment(corner1,corner2,x,xp); - distsq = closest(x,xp,nearest,distsq); - crad = -2.0*(radiuslo + (nearest[2]-lo)*(radiushi-radiuslo)/(hi-lo)); + point_on_line_segment(corner1, corner2, x, xp); + distsq = closest(x, xp, nearest, distsq); + crad = -2.0 * (radiuslo + (nearest[2] - lo) * (radiushi - radiuslo) / (hi - lo)); } if (!open_faces[0]) { - point_on_line_segment(corner1,corner3,x,xp); + point_on_line_segment(corner1, corner3, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } if (!open_faces[1]) { - point_on_line_segment(corner2,corner4,x,xp); + point_on_line_segment(corner2, corner4, x, xp); distsqprev = distsq; - distsq = closest(x,xp,nearest,distsq); + distsq = closest(x, xp, nearest, distsq); if (distsq < distsqprev) crad = 0; } - add_contact(0,x,nearest[0],nearest[1],nearest[2]); + add_contact(0, x, nearest[0], nearest[1], nearest[2]); contact[0].radius = crad; contact[0].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -610,7 +629,7 @@ double RegCone::closest(double *x, double *near, double *nearest, double dsq) double delx = x[0] - near[0]; double dely = x[1] - near[1]; double delz = x[2] - near[2]; - double rsq = delx*delx + dely*dely + delz*delz; + double rsq = delx * delx + dely * dely + delz * delz; if (rsq >= dsq) return dsq; nearest[0] = near[0]; diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index 259b855a66..c6119c0c3b 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -25,94 +24,97 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 -enum{CONSTANT,VARIABLE}; +static constexpr double BIG = 1.0e20; + +enum { CONSTANT, VARIABLE }; /* ---------------------------------------------------------------------- */ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr) + Region(lmp, narg, arg), c1str(nullptr), c2str(nullptr), rstr(nullptr) { - options(narg-8,&arg[8]); + options(narg - 8, &arg[8]); // check open face settings if (openflag) - for (int i=3; i<6; i++) - if (open_faces[i]) error->all(FLERR,"Illegal region cylinder open face: {}", i+1); + for (int i = 3; i < 6; i++) + if (open_faces[i]) error->all(FLERR, "Illegal region cylinder open face: {}", i + 1); - if (strcmp(arg[2],"x") != 0 && strcmp(arg[2],"y") != 0 && strcmp(arg[2],"z") != 0) - error->all(FLERR,"Illegal region cylinder axis: {}", arg[2]); + if (strcmp(arg[2], "x") != 0 && strcmp(arg[2], "y") != 0 && strcmp(arg[2], "z") != 0) + error->all(FLERR, "Illegal region cylinder axis: {}", arg[2]); axis = arg[2][0]; if (axis == 'x') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = yscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = yscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } else if (axis == 'y') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = zscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = zscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } else if (axis == 'z') { - if (utils::strmatch(arg[3],"^v_")) { - c1str = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + c1str = utils::strdup(arg[3] + 2); c1 = 0.0; c1style = VARIABLE; varshape = 1; } else { - c1 = xscale*utils::numeric(FLERR,arg[3],false,lmp); + c1 = xscale * utils::numeric(FLERR, arg[3], false, lmp); c1style = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - c2str = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + c2str = utils::strdup(arg[4] + 2); c2 = 0.0; c2style = VARIABLE; varshape = 1; } else { - c2 = yscale*utils::numeric(FLERR,arg[4],false,lmp); + c2 = yscale * utils::numeric(FLERR, arg[4], false, lmp); c2style = CONSTANT; } } - if (utils::strmatch(arg[5],"^v_")) { - rstr = utils::strdup(arg[5]+2); + if (utils::strmatch(arg[5], "^v_")) { + rstr = utils::strdup(arg[5] + 2); radius = 0.0; rstyle = VARIABLE; varshape = 1; } else { - radius = utils::numeric(FLERR,arg[5],false,lmp); - if (axis == 'x') radius *= yscale; - else radius *= xscale; + radius = utils::numeric(FLERR, arg[5], false, lmp); + if (axis == 'x') + radius *= yscale; + else + radius *= xscale; rstyle = CONSTANT; } @@ -121,57 +123,75 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : RegCylinder::shape_update(); } - if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { + if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); if (axis == 'x') { - if (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[0]; - else lo = domain->boxlo_bound[0]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[0]; + else + lo = domain->boxlo_bound[0]; } if (axis == 'y') { - if (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[1]; - else lo = domain->boxlo_bound[1]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[1]; + else + lo = domain->boxlo_bound[1]; } if (axis == 'z') { - if (strcmp(arg[6],"INF") == 0) lo = -BIG; - else if (domain->triclinic == 0) lo = domain->boxlo[2]; - else lo = domain->boxlo_bound[2]; + if (strcmp(arg[6], "INF") == 0) + lo = -BIG; + else if (domain->triclinic == 0) + lo = domain->boxlo[2]; + else + lo = domain->boxlo_bound[2]; } } else { - if (axis == 'x') lo = xscale*utils::numeric(FLERR,arg[6],false,lmp); - if (axis == 'y') lo = yscale*utils::numeric(FLERR,arg[6],false,lmp); - if (axis == 'z') lo = zscale*utils::numeric(FLERR,arg[6],false,lmp); + if (axis == 'x') lo = xscale * utils::numeric(FLERR, arg[6], false, lmp); + if (axis == 'y') lo = yscale * utils::numeric(FLERR, arg[6], false, lmp); + if (axis == 'z') lo = zscale * utils::numeric(FLERR, arg[6], false, lmp); } - if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { + if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); if (axis == 'x') { - if (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[0]; - else hi = domain->boxhi_bound[0]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[0]; + else + hi = domain->boxhi_bound[0]; } if (axis == 'y') { - if (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[1]; - else hi = domain->boxhi_bound[1]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[1]; + else + hi = domain->boxhi_bound[1]; } if (axis == 'z') { - if (strcmp(arg[7],"INF") == 0) hi = BIG; - else if (domain->triclinic == 0) hi = domain->boxhi[2]; - else hi = domain->boxhi_bound[2]; + if (strcmp(arg[7], "INF") == 0) + hi = BIG; + else if (domain->triclinic == 0) + hi = domain->boxhi[2]; + else + hi = domain->boxhi_bound[2]; } } else { - if (axis == 'x') hi = xscale*utils::numeric(FLERR,arg[7],false,lmp); - if (axis == 'y') hi = yscale*utils::numeric(FLERR,arg[7],false,lmp); - if (axis == 'z') hi = zscale*utils::numeric(FLERR,arg[7],false,lmp); + if (axis == 'x') hi = xscale * utils::numeric(FLERR, arg[7], false, lmp); + if (axis == 'y') hi = yscale * utils::numeric(FLERR, arg[7], false, lmp); + if (axis == 'z') hi = zscale * utils::numeric(FLERR, arg[7], false, lmp); } // error check - if (radius <= 0.0) error->all(FLERR,"Illegal radius {} in region cylinder command", radius); + if (radius <= 0.0) error->all(FLERR, "Illegal radius {} in region cylinder command", radius); // extent of cylinder // for variable radius, uses initial radius @@ -202,25 +222,28 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : extent_zlo = lo; extent_zhi = hi; } - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to cylinder surface and 2 ends // particle can only touch surface and 1 end cmax = 3; contact = new Contact[cmax]; - if (interior) tmax = 2; - else tmax = 1; + if (interior) + tmax = 2; + else + tmax = 1; } /* ---------------------------------------------------------------------- */ RegCylinder::~RegCylinder() { - delete [] c1str; - delete [] c2str; - delete [] rstr; - delete [] contact; + delete[] c1str; + delete[] c2str; + delete[] rstr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -238,27 +261,33 @@ void RegCylinder::init() int RegCylinder::inside(double x, double y, double z) { - double del1,del2,dist; + double del1, del2, dist; int inside; if (axis == 'x') { del1 = y - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && x >= lo && x <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && x >= lo && x <= hi) + inside = 1; + else + inside = 0; } else if (axis == 'y') { del1 = x - c1; del2 = z - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && y >= lo && y <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && y >= lo && y <= hi) + inside = 1; + else + inside = 0; } else { del1 = x - c1; del2 = y - c2; - dist = sqrt(del1*del1 + del2*del2); - if (dist <= radius && z >= lo && z <= hi) inside = 1; - else inside = 0; + dist = sqrt(del1 * del1 + del2 * del2); + if (dist <= radius && z >= lo && z <= hi) + inside = 1; + else + inside = 0; } return inside; @@ -274,14 +303,14 @@ int RegCylinder::inside(double x, double y, double z) int RegCylinder::surface_interior(double *x, double cutoff) { - double del1,del2,r,delta; + double del1, del2, r, delta; int n = 0; if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // x is exterior to cylinder @@ -293,9 +322,9 @@ int RegCylinder::surface_interior(double *x, double cutoff) if (delta < cutoff && r > 0.0 && !open_faces[2]) { contact[n].r = delta; contact[n].delx = 0.0; - contact[n].dely = del1*(1.0-radius/r); - contact[n].delz = del2*(1.0-radius/r); - contact[n].radius = -2.0*radius; + contact[n].dely = del1 * (1.0 - radius / r); + contact[n].delz = del2 * (1.0 - radius / r); + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -324,7 +353,7 @@ int RegCylinder::surface_interior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // y is exterior to cylinder @@ -335,10 +364,10 @@ int RegCylinder::surface_interior(double *x, double cutoff) delta = radius - r; if (delta < cutoff && r > 0.0 && !open_faces[2]) { contact[n].r = delta; - contact[n].delx = del1*(1.0-radius/r); + contact[n].delx = del1 * (1.0 - radius / r); contact[n].dely = 0.0; - contact[n].delz = del2*(1.0-radius/r); - contact[n].radius = -2.0*radius; + contact[n].delz = del2 * (1.0 - radius / r); + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -367,7 +396,7 @@ int RegCylinder::surface_interior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // z is exterior to cylinder @@ -378,10 +407,10 @@ int RegCylinder::surface_interior(double *x, double cutoff) delta = radius - r; if (delta < cutoff && r > 0.0 && !open_faces[2]) { contact[n].r = delta; - contact[n].delx = del1*(1.0-radius/r); - contact[n].dely = del2*(1.0-radius/r); + contact[n].delx = del1 * (1.0 - radius / r); + contact[n].dely = del2 * (1.0 - radius / r); contact[n].delz = 0.0; - contact[n].radius = -2.0*radius; + contact[n].radius = -2.0 * radius; contact[n].iwall = 2; contact[n].varflag = 1; n++; @@ -419,12 +448,12 @@ int RegCylinder::surface_interior(double *x, double cutoff) int RegCylinder::surface_exterior(double *x, double cutoff) { - double del1,del2,r; - double xp,yp,zp; + double del1, del2, r; + double xp, yp, zp; double dx, dr, dr2, d2, d2prev; - // radius of curvature for granular - // 0 for flat surfaces (infinite case), 2*radius for curved portion + // radius of curvature for granular + // 0 for flat surfaces (infinite case), 2*radius for curved portion double crad = 0.0; int varflag = 0; @@ -432,12 +461,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (axis == 'x') { del1 = x[1] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // x is far enough from cylinder that there is no contact // x is interior to cylinder - if (r >= radius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0; + if (r >= radius + cutoff || x[0] <= lo - cutoff || x[0] >= hi + cutoff) return 0; if (r < radius && x[0] > lo && x[0] < hi) return 0; // x is exterior to cylinder or on its surface @@ -448,48 +477,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - yp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; - crad = 2.0*radius; + yp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { yp = x[1]; zp = x[2]; } - if (x[0] < lo) xp = lo; - else if (x[0] > hi) xp = hi; - else xp = x[0]; + if (x[0] < lo) + xp = lo; + else if (x[0] > hi) + xp = hi; + else + xp = x[0]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - yp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; + yp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; if (x[0] < lo) { - dx = lo-x[0]; + dx = lo - x[0]; xp = lo; - } - else if (x[0] > hi) { - dx = x[0]-hi; + } else if (x[0] > hi) { + dx = x[0] - hi; xp = hi; - } - else { + } else { dx = 0; xp = x[0]; } - d2 = d2prev = dr2 + dx*dx; + d2 = d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[0]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { xp = lo; if (r < radius) { @@ -504,8 +536,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[0]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { xp = hi; if (r < radius) { @@ -516,7 +550,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } } - add_contact(0,x,xp,yp,zp); + add_contact(0, x, xp, yp, zp); contact[0].radius = crad; contact[0].varflag = varflag; contact[0].iwall = 0; @@ -526,12 +560,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } else if (axis == 'y') { del1 = x[0] - c1; del2 = x[2] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // y is far enough from cylinder that there is no contact // y is interior to cylinder - if (r >= radius+cutoff || x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0; + if (r >= radius + cutoff || x[1] <= lo - cutoff || x[1] >= hi + cutoff) return 0; if (r < radius && x[1] > lo && x[1] < hi) return 0; // y is exterior to cylinder or on its surface @@ -542,48 +576,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - xp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; - crad = 2.0*radius; + xp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { xp = x[0]; zp = x[2]; } - if (x[1] < lo) yp = lo; - else if (x[1] > hi) yp = hi; - else yp = x[1]; + if (x[1] < lo) + yp = lo; + else if (x[1] > hi) + yp = hi; + else + yp = x[1]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - xp = c1 + del1*radius/r; - zp = c2 + del2*radius/r; + xp = c1 + del1 * radius / r; + zp = c2 + del2 * radius / r; if (x[1] < lo) { - dx = lo-x[1]; + dx = lo - x[1]; yp = lo; - } - else if (x[1] > hi) { - dx = x[1]-hi; + } else if (x[1] > hi) { + dx = x[1] - hi; yp = hi; - } - else { + } else { dx = 0; yp = x[1]; } - d2 = d2prev = dr2 + dx*dx; + d2 = d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[1]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { yp = lo; if (r < radius) { @@ -598,8 +635,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[1]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { yp = hi; if (r < radius) { @@ -610,7 +649,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } } - add_contact(0,x,xp,yp,zp); + add_contact(0, x, xp, yp, zp); contact[0].radius = crad; contact[0].varflag = varflag; contact[0].iwall = 0; @@ -620,12 +659,12 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } else { del1 = x[0] - c1; del2 = x[1] - c2; - r = sqrt(del1*del1 + del2*del2); + r = sqrt(del1 * del1 + del2 * del2); // z is far enough from cylinder that there is no contact // z is interior to cylinder - if (r >= radius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0; + if (r >= radius + cutoff || x[2] <= lo - cutoff || x[2] >= hi + cutoff) return 0; if (r < radius && x[2] > lo && x[2] < hi) return 0; // z is exterior to cylinder or on its surface @@ -636,46 +675,51 @@ int RegCylinder::surface_exterior(double *x, double cutoff) d2prev = BIG; if (!openflag) { if (r > radius) { - xp = c1 + del1*radius/r; - yp = c2 + del2*radius/r; - crad = 2.0*radius; + xp = c1 + del1 * radius / r; + yp = c2 + del2 * radius / r; + crad = 2.0 * radius; varflag = 1; } else { xp = x[0]; yp = x[1]; } - if (x[2] < lo) zp = lo; - else if (x[2] > hi) zp = hi; - else zp = x[2]; + if (x[2] < lo) + zp = lo; + else if (x[2] > hi) + zp = hi; + else + zp = x[2]; } else { - // closest point on curved surface + // closest point on curved surface dr = r - radius; - dr2 = dr*dr; + dr2 = dr * dr; if (!open_faces[2]) { - xp = c1 + del1*radius/r; - yp = c2 + del2*radius/r; + xp = c1 + del1 * radius / r; + yp = c2 + del2 * radius / r; if (x[2] < lo) { - dx = lo-x[2]; + dx = lo - x[2]; zp = lo; } else if (x[2] > hi) { - dx = x[2]-hi; + dx = x[2] - hi; zp = hi; } else { dx = 0; zp = x[2]; } - d2prev = dr2 + dx*dx; + d2prev = dr2 + dx * dx; } // closest point on bottom cap if (!open_faces[0]) { dx = lo - x[2]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { zp = lo; if (r < radius) { @@ -690,8 +734,10 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (!open_faces[1]) { dx = hi - x[2]; - if (r < radius) d2 = dx*dx; - else d2 = dr2 + dx*dx; + if (r < radius) + d2 = dx * dx; + else + d2 = dr2 + dx * dx; if (d2 < d2prev) { zp = hi; if (r < radius) { @@ -702,7 +748,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff) } } - add_contact(0,x,xp,yp,zp); + add_contact(0, x, xp, yp, zp); contact[0].radius = crad; contact[0].varflag = varflag; contact[0].iwall = 0; @@ -721,22 +767,21 @@ void RegCylinder::shape_update() if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var); if (rstyle == VARIABLE) { radius = input->variable->compute_equal(rvar); - if (radius < 0.0) - error->one(FLERR,"Variable evaluation in region gave bad value"); + if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value"); } if (axis == 'x') { if (c1style == VARIABLE) c1 *= yscale; if (c2style == VARIABLE) c2 *= zscale; - if (rstyle == VARIABLE) radius *= yscale; + if (rstyle == VARIABLE) radius *= yscale; } else if (axis == 'y') { if (c1style == VARIABLE) c1 *= xscale; if (c2style == VARIABLE) c2 *= zscale; - if (rstyle == VARIABLE) radius *= xscale; - } else { // axis == 'z' + if (rstyle == VARIABLE) radius *= xscale; + } else { // axis == 'z' if (c1style == VARIABLE) c1 *= xscale; if (c2style == VARIABLE) c2 *= yscale; - if (rstyle == VARIABLE) radius *= xscale; + if (rstyle == VARIABLE) radius *= xscale; } } @@ -748,30 +793,26 @@ void RegCylinder::variable_check() { if (c1style == VARIABLE) { c1var = input->variable->find(c1str); - if (c1var < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (c1var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); if (!input->variable->equalstyle(c1var)) - error->all(FLERR,"Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable for region cylinder is invalid style"); } if (c2style == VARIABLE) { c2var = input->variable->find(c2str); - if (c2var < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (c2var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); if (!input->variable->equalstyle(c2var)) - error->all(FLERR,"Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable for region cylinder is invalid style"); } if (rstyle == VARIABLE) { rvar = input->variable->find(rstr); - if (rvar < 0) - error->all(FLERR,"Variable name for region cylinder does not exist"); + if (rvar < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); if (!input->variable->equalstyle(rvar)) - error->all(FLERR,"Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable for region cylinder is invalid style"); } } - /* ---------------------------------------------------------------------- Set values needed to calculate velocity due to shape changes. These values do not depend on the contact, so this function is @@ -795,35 +836,34 @@ void RegCylinder::set_velocity_shape() xcenter[2] = 0; } forward_transform(xcenter[0], xcenter[1], xcenter[2]); - if (update->ntimestep > 0) rprev = prev[4]; - else rprev = radius; + if (update->ntimestep > 0) + rprev = prev[4]; + else + rprev = radius; prev[4] = radius; } - - /* ---------------------------------------------------------------------- add velocity due to shape change to wall velocity ------------------------------------------------------------------------- */ void RegCylinder::velocity_contact_shape(double *vwall, double *xc) { - double delx, dely, delz; // Displacement of contact point in x,y,z + double delx, dely, delz; // Displacement of contact point in x,y,z if (axis == 'x') { delx = 0; - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); } else if (axis == 'y') { - delx = (xc[0] - xcenter[0])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); dely = 0; - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); } else { - delx = (xc[0] - xcenter[0])*(1 - rprev/radius); - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); delz = 0; } - vwall[0] += delx/update->dt; - vwall[1] += dely/update->dt; - vwall[2] += delz/update->dt; + vwall[0] += delx / update->dt; + vwall[1] += dely / update->dt; + vwall[2] += delz / update->dt; } - diff --git a/src/region_prism.cpp b/src/region_prism.cpp index 7c06c0aa7d..da2048cce0 100644 --- a/src/region_prism.cpp +++ b/src/region_prism.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -26,105 +25,126 @@ using namespace LAMMPS_NS; -#define BIG 1.0e20 +static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) { - options(narg-11,&arg[11]); + options(narg - 11, &arg[11]); - if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { + if (strcmp(arg[2], "INF") == 0 || strcmp(arg[2], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[2],"INF") == 0) xlo = -BIG; - else xlo = domain->boxlo[0]; - } else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[2], "INF") == 0) + xlo = -BIG; + else + xlo = domain->boxlo[0]; + } else + xlo = xscale * utils::numeric(FLERR, arg[2], false, lmp); - if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) { + if (strcmp(arg[3], "INF") == 0 || strcmp(arg[3], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[3],"INF") == 0) xhi = BIG; - else xhi = domain->boxhi[0]; - } else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[3], "INF") == 0) + xhi = BIG; + else + xhi = domain->boxhi[0]; + } else + xhi = xscale * utils::numeric(FLERR, arg[3], false, lmp); - if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) { + if (strcmp(arg[4], "INF") == 0 || strcmp(arg[4], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[4],"INF") == 0) ylo = -BIG; - else ylo = domain->boxlo[1]; - } else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[4], "INF") == 0) + ylo = -BIG; + else + ylo = domain->boxlo[1]; + } else + ylo = yscale * utils::numeric(FLERR, arg[4], false, lmp); - if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) { + if (strcmp(arg[5], "INF") == 0 || strcmp(arg[5], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[5],"INF") == 0) yhi = BIG; - else yhi = domain->boxhi[1]; - } else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[5], "INF") == 0) + yhi = BIG; + else + yhi = domain->boxhi[1]; + } else + yhi = yscale * utils::numeric(FLERR, arg[5], false, lmp); - if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { + if (strcmp(arg[6], "INF") == 0 || strcmp(arg[6], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[6],"INF") == 0) zlo = -BIG; - else zlo = domain->boxlo[2]; - } else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[6], "INF") == 0) + zlo = -BIG; + else + zlo = domain->boxlo[2]; + } else + zlo = zscale * utils::numeric(FLERR, arg[6], false, lmp); - if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) { + if (strcmp(arg[7], "INF") == 0 || strcmp(arg[7], "EDGE") == 0) { if (domain->box_exist == 0) - error->all(FLERR,"Cannot use region INF or EDGE when box does not exist"); - if (strcmp(arg[7],"INF") == 0) zhi = BIG; - else zhi = domain->boxhi[2]; - } else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp); + error->all(FLERR, "Cannot use region INF or EDGE when box does not exist"); + if (strcmp(arg[7], "INF") == 0) + zhi = BIG; + else + zhi = domain->boxhi[2]; + } else + zhi = zscale * utils::numeric(FLERR, arg[7], false, lmp); - xy = xscale*utils::numeric(FLERR,arg[8],false,lmp); - xz = xscale*utils::numeric(FLERR,arg[9],false,lmp); - yz = yscale*utils::numeric(FLERR,arg[10],false,lmp); + xy = xscale * utils::numeric(FLERR, arg[8], false, lmp); + xz = xscale * utils::numeric(FLERR, arg[9], false, lmp); + yz = yscale * utils::numeric(FLERR, arg[10], false, lmp); // error check // prism cannot be 0 thickness in any dim, else inverse blows up // non-zero tilt values cannot be used if either dim is INF on both ends - if (xlo >= xhi) error->all(FLERR,"Illegal region prism xlo: {} >= xhi: {}", xlo, xhi); - if (ylo >= yhi) error->all(FLERR,"Illegal region prism ylo: {} >= yhi: {}", ylo, yhi); - if (zlo >= zhi) error->all(FLERR,"Illegal region prism zlo: {} >= zhi: {}", zlo ,zhi); + if (xlo >= xhi) error->all(FLERR, "Illegal region prism xlo: {} >= xhi: {}", xlo, xhi); + if (ylo >= yhi) error->all(FLERR, "Illegal region prism ylo: {} >= yhi: {}", ylo, yhi); + if (zlo >= zhi) error->all(FLERR, "Illegal region prism zlo: {} >= zhi: {}", zlo, zhi); if (xy != 0.0 && xlo == -BIG && xhi == BIG) - error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite x size"); + error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite x size"); if (xy != 0.0 && ylo == -BIG && yhi == BIG) - error->all(FLERR,"Illegal region prism non-zero xy tilt with infinite y size"); + error->all(FLERR, "Illegal region prism non-zero xy tilt with infinite y size"); if (xz != 0.0 && xlo == -BIG && xhi == BIG) - error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite x size"); + error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite x size"); if (xz != 0.0 && zlo == -BIG && zhi == BIG) - error->all(FLERR,"Illegal region prism non-zero xz tilt with infinite z size"); + error->all(FLERR, "Illegal region prism non-zero xz tilt with infinite z size"); if (yz != 0.0 && ylo == -BIG && yhi == BIG) - error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite y size"); + error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite y size"); if (yz != 0.0 && zlo == -BIG && zhi == BIG) - error->all(FLERR,"Illegal region prism non-zero yz tilt with infinite z size"); + error->all(FLERR, "Illegal region prism non-zero yz tilt with infinite z size"); // extent of prism if (interior) { bboxflag = 1; - extent_xlo = MIN(xlo,xlo+xy); - extent_xlo = MIN(extent_xlo,extent_xlo+xz); - extent_ylo = MIN(ylo,ylo+yz); + extent_xlo = MIN(xlo, xlo + xy); + extent_xlo = MIN(extent_xlo, extent_xlo + xz); + extent_ylo = MIN(ylo, ylo + yz); extent_zlo = zlo; - extent_xhi = MAX(xhi,xhi+xy); - extent_xhi = MAX(extent_xhi,extent_xhi+xz); - extent_yhi = MAX(yhi,yhi+yz); + extent_xhi = MAX(xhi, xhi + xy); + extent_xhi = MAX(extent_xhi, extent_xhi + xz); + extent_yhi = MAX(yhi, yhi + yz); extent_zhi = zhi; - } else bboxflag = 0; + } else + bboxflag = 0; // particle could be close to all 6 planes // particle can only touch 3 planes cmax = 6; contact = new Contact[cmax]; - if (interior) tmax = 3; - else tmax = 1; + if (interior) + tmax = 3; + else + tmax = 1; // h = transformation matrix from tilt coords (0-1) to box coords (xyz) // columns of h are edge vectors of tilted box @@ -140,26 +160,26 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) h[1][2] = yz; h[2][2] = zhi - zlo; - hinv[0][0] = 1.0/h[0][0]; - hinv[0][1] = -h[0][1] / (h[0][0]*h[1][1]); - hinv[0][2] = (h[0][1]*h[1][2] - h[0][2]*h[1][1]) / (h[0][0]*h[1][1]*h[2][2]); - hinv[1][1] = 1.0/h[1][1]; - hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]); - hinv[2][2] = 1.0/h[2][2]; + hinv[0][0] = 1.0 / h[0][0]; + hinv[0][1] = -h[0][1] / (h[0][0] * h[1][1]); + hinv[0][2] = (h[0][1] * h[1][2] - h[0][2] * h[1][1]) / (h[0][0] * h[1][1] * h[2][2]); + hinv[1][1] = 1.0 / h[1][1]; + hinv[1][2] = -h[1][2] / (h[1][1] * h[2][2]); + hinv[2][2] = 1.0 / h[2][2]; // corners = 8 corner points of prism // order = x varies fastest, then y, finally z // clo/chi = lo and hi corner pts of prism - a[0] = xhi-xlo; + a[0] = xhi - xlo; a[1] = 0.0; a[2] = 0.0; b[0] = xy; - b[1] = yhi-ylo; + b[1] = yhi - ylo; b[2] = 0.0; c[0] = xz; c[1] = yz; - c[2] = zhi-zlo; + c[2] = zhi - zlo; clo[0] = corners[0][0] = xlo; clo[1] = corners[0][1] = ylo; @@ -196,19 +216,18 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) // face = 6 inward-facing unit normals to prism faces // order = xy plane, xz plane, yz plane - MathExtra::cross3(a,b,face[0]); - MathExtra::cross3(b,a,face[1]); - MathExtra::cross3(c,a,face[2]); - MathExtra::cross3(a,c,face[3]); - MathExtra::cross3(b,c,face[4]); - MathExtra::cross3(c,b,face[5]); + MathExtra::cross3(a, b, face[0]); + MathExtra::cross3(b, a, face[1]); + MathExtra::cross3(c, a, face[2]); + MathExtra::cross3(a, c, face[3]); + MathExtra::cross3(b, c, face[4]); + MathExtra::cross3(c, b, face[5]); // remap open face indices to be consistent if (openflag) { int temp[6]; - for (int i = 0; i < 6; i++) - temp[i] = open_faces[i]; + for (int i = 0; i < 6; i++) temp[i] = open_faces[i]; open_faces[0] = temp[4]; open_faces[1] = temp[5]; open_faces[2] = temp[2]; @@ -223,27 +242,51 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) // verts in each tri are ordered so that right-hand rule gives inward norm // order = xy plane, xz plane, yz plane - tri[0][0] = 0; tri[0][1] = 1; tri[0][2] = 3; - tri[1][0] = 0; tri[1][1] = 3; tri[1][2] = 2; - tri[2][0] = 4; tri[2][1] = 7; tri[2][2] = 5; - tri[3][0] = 4; tri[3][1] = 6; tri[3][2] = 7; + tri[0][0] = 0; + tri[0][1] = 1; + tri[0][2] = 3; + tri[1][0] = 0; + tri[1][1] = 3; + tri[1][2] = 2; + tri[2][0] = 4; + tri[2][1] = 7; + tri[2][2] = 5; + tri[3][0] = 4; + tri[3][1] = 6; + tri[3][2] = 7; - tri[4][0] = 0; tri[4][1] = 4; tri[4][2] = 5; - tri[5][0] = 0; tri[5][1] = 5; tri[5][2] = 1; - tri[6][0] = 2; tri[6][1] = 7; tri[6][2] = 6; - tri[7][0] = 2; tri[7][1] = 3; tri[7][2] = 7; + tri[4][0] = 0; + tri[4][1] = 4; + tri[4][2] = 5; + tri[5][0] = 0; + tri[5][1] = 5; + tri[5][2] = 1; + tri[6][0] = 2; + tri[6][1] = 7; + tri[6][2] = 6; + tri[7][0] = 2; + tri[7][1] = 3; + tri[7][2] = 7; - tri[8][0] = 2; tri[8][1] = 6; tri[8][2] = 4; - tri[9][0] = 2; tri[9][1] = 4; tri[9][2] = 0; - tri[10][0] = 1; tri[10][1] = 5; tri[10][2] = 7; - tri[11][0] = 1; tri[11][1] = 7; tri[11][2] = 3; + tri[8][0] = 2; + tri[8][1] = 6; + tri[8][2] = 4; + tri[9][0] = 2; + tri[9][1] = 4; + tri[9][2] = 0; + tri[10][0] = 1; + tri[10][1] = 5; + tri[10][2] = 7; + tri[11][0] = 1; + tri[11][1] = 7; + tri[11][2] = 3; } /* ---------------------------------------------------------------------- */ RegPrism::~RegPrism() { - delete [] contact; + delete[] contact; } /* ---------------------------------------------------------------------- @@ -258,12 +301,11 @@ RegPrism::~RegPrism() int RegPrism::inside(double x, double y, double z) { - double a = hinv[0][0]*(x-xlo) + hinv[0][1]*(y-ylo) + hinv[0][2]*(z-zlo); - double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo); - double c = hinv[2][2]*(z-zlo); + double a = hinv[0][0] * (x - xlo) + hinv[0][1] * (y - ylo) + hinv[0][2] * (z - zlo); + double b = hinv[1][1] * (y - ylo) + hinv[1][2] * (z - zlo); + double c = hinv[2][2] * (z - zlo); - if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) - return 1; + if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) return 1; return 0; } @@ -283,10 +325,12 @@ int RegPrism::surface_interior(double *x, double cutoff) // x is exterior to prism for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else corner = clo; - dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + - (x[2]-corner[2])*face[i][2]; + if (i % 2) + corner = chi; + else + corner = clo; + dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] + + (x[2] - corner[2]) * face[i][2]; if (dot < 0.0) return 0; } @@ -296,15 +340,17 @@ int RegPrism::surface_interior(double *x, double cutoff) for (i = 0; i < 6; i++) { if (open_faces[i]) continue; - if (i % 2) corner = chi; - else corner = clo; - dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + - (x[2]-corner[2])*face[i][2]; + if (i % 2) + corner = chi; + else + corner = clo; + dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] + + (x[2] - corner[2]) * face[i][2]; if (dot < cutoff) { contact[n].r = dot; - contact[n].delx = dot*face[i][0]; - contact[n].dely = dot*face[i][1]; - contact[n].delz = dot*face[i][2]; + contact[n].delx = dot * face[i][0]; + contact[n].dely = dot * face[i][1]; + contact[n].delz = dot * face[i][2]; contact[n].radius = 0; contact[n].iwall = i; n++; @@ -325,25 +371,29 @@ int RegPrism::surface_exterior(double *x, double cutoff) int i; double dot; double *corner; - double xp,yp,zp; + double xp, yp, zp; // x is far enough from prism that there is no contact for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else corner = clo; - dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + - (x[2]-corner[2])*face[i][2]; + if (i % 2) + corner = chi; + else + corner = clo; + dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] + + (x[2] - corner[2]) * face[i][2]; if (dot <= -cutoff) return 0; } // x is interior to prism for (i = 0; i < 6; i++) { - if (i % 2) corner = chi; - else corner = clo; - dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + - (x[2]-corner[2])*face[i][2]; + if (i % 2) + corner = chi; + else + corner = clo; + dot = (x[0] - corner[0]) * face[i][0] + (x[1] - corner[1]) * face[i][1] + + (x[2] - corner[2]) * face[i][2]; if (dot <= 0.0) break; } @@ -354,8 +404,8 @@ int RegPrism::surface_exterior(double *x, double cutoff) // could be edge or corner pt of prism // do not add contact point if r >= cutoff - find_nearest(x,xp,yp,zp); - add_contact(0,x,xp,yp,zp); + find_nearest(x, xp, yp, zp); + add_contact(0, x, xp, yp, zp); contact[0].radius = 0; contact[0].iwall = 0; if (contact[0].r < cutoff) return 1; @@ -369,8 +419,8 @@ int RegPrism::surface_exterior(double *x, double cutoff) void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp) { - int i,j,k,iface; - double xproj[3],xline[3],nearest[3]; + int i, j, k, iface; + double xproj[3], xline[3], nearest[3]; double dot; // generate successive xnear points, one nearest to x is (xp,yp,zp) @@ -384,27 +434,25 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp) double distsq = BIG; for (int itri = 0; itri < 12; itri++) { - iface = itri/2; + iface = itri / 2; if (open_faces[iface]) continue; i = tri[itri][0]; j = tri[itri][1]; k = tri[itri][2]; - dot = (x[0]-corners[i][0])*face[iface][0] + - (x[1]-corners[i][1])*face[iface][1] + - (x[2]-corners[i][2])*face[iface][2]; - xproj[0] = x[0] - dot*face[iface][0]; - xproj[1] = x[1] - dot*face[iface][1]; - xproj[2] = x[2] - dot*face[iface][2]; - if (inside_tri(xproj,corners[i],corners[j],corners[k],face[iface])) { - distsq = closest(x,xproj,nearest,distsq); - } - else { - point_on_line_segment(corners[i],corners[j],xproj,xline); - distsq = closest(x,xline,nearest,distsq); - point_on_line_segment(corners[j],corners[k],xproj,xline); - distsq = closest(x,xline,nearest,distsq); - point_on_line_segment(corners[i],corners[k],xproj,xline); - distsq = closest(x,xline,nearest,distsq); + dot = (x[0] - corners[i][0]) * face[iface][0] + (x[1] - corners[i][1]) * face[iface][1] + + (x[2] - corners[i][2]) * face[iface][2]; + xproj[0] = x[0] - dot * face[iface][0]; + xproj[1] = x[1] - dot * face[iface][1]; + xproj[2] = x[2] - dot * face[iface][2]; + if (inside_tri(xproj, corners[i], corners[j], corners[k], face[iface])) { + distsq = closest(x, xproj, nearest, distsq); + } else { + point_on_line_segment(corners[i], corners[j], xproj, xline); + distsq = closest(x, xline, nearest, distsq); + point_on_line_segment(corners[j], corners[k], xproj, xline); + distsq = closest(x, xline, nearest, distsq); + point_on_line_segment(corners[i], corners[k], xproj, xline); + distsq = closest(x, xline, nearest, distsq); } } @@ -422,25 +470,24 @@ void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp) if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle ------------------------------------------------------------------------- */ -int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, - double *norm) +int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3, double *norm) { - double edge[3],pvec[3],xproduct[3]; + double edge[3], pvec[3], xproduct[3]; - MathExtra::sub3(v2,v1,edge); - MathExtra::sub3(x,v1,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v2, v1, edge); + MathExtra::sub3(x, v1, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; - MathExtra::sub3(v3,v2,edge); - MathExtra::sub3(x,v2,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v3, v2, edge); + MathExtra::sub3(x, v2, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; - MathExtra::sub3(v1,v3,edge); - MathExtra::sub3(x,v3,pvec); - MathExtra::cross3(edge,pvec,xproduct); - if (MathExtra::dot3(xproduct,norm) < 0.0) return 0; + MathExtra::sub3(v1, v3, edge); + MathExtra::sub3(x, v3, pvec); + MathExtra::cross3(edge, pvec, xproduct); + if (MathExtra::dot3(xproduct, norm) < 0.0) return 0; return 1; } @@ -452,7 +499,7 @@ double RegPrism::closest(double *x, double *near, double *nearest, double dsq) double delx = x[0] - near[0]; double dely = x[1] - near[1]; double delz = x[2] - near[2]; - double rsq = delx*delx + dely*dely + delz*delz; + double rsq = delx * delx + dely * dely + delz * delz; if (rsq >= dsq) return dsq; nearest[0] = near[0]; diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index 03eb762266..b67e805ac4 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,52 +22,52 @@ using namespace LAMMPS_NS; -enum{CONSTANT,VARIABLE}; +enum { CONSTANT, VARIABLE }; /* ---------------------------------------------------------------------- */ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr) + Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr) { - options(narg-6,&arg[6]); + options(narg - 6, &arg[6]); - if (utils::strmatch(arg[2],"^v_")) { - xstr = utils::strdup(arg[2]+2); + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); xc = 0.0; xstyle = VARIABLE; varshape = 1; } else { - xc = xscale*utils::numeric(FLERR,arg[2],false,lmp); + xc = xscale * utils::numeric(FLERR, arg[2], false, lmp); xstyle = CONSTANT; } - if (utils::strmatch(arg[3],"^v_")) { - ystr = utils::strdup(arg[3]+2); + if (utils::strmatch(arg[3], "^v_")) { + ystr = utils::strdup(arg[3] + 2); yc = 0.0; ystyle = VARIABLE; varshape = 1; } else { - yc = yscale*utils::numeric(FLERR,arg[3],false,lmp); + yc = yscale * utils::numeric(FLERR, arg[3], false, lmp); ystyle = CONSTANT; } - if (utils::strmatch(arg[4],"^v_")) { - zstr = utils::strdup(arg[4]+2); + if (utils::strmatch(arg[4], "^v_")) { + zstr = utils::strdup(arg[4] + 2); zc = 0.0; zstyle = VARIABLE; varshape = 1; } else { - zc = zscale*utils::numeric(FLERR,arg[4],false,lmp); + zc = zscale * utils::numeric(FLERR, arg[4], false, lmp); zstyle = CONSTANT; } - if (utils::strmatch(arg[5],"^v_")) { - rstr = utils::strdup(arg[5]+2); + if (utils::strmatch(arg[5], "^v_")) { + rstr = utils::strdup(arg[5] + 2); radius = 0.0; rstyle = VARIABLE; varshape = 1; } else { - radius = xscale*utils::numeric(FLERR,arg[5],false,lmp); + radius = xscale * utils::numeric(FLERR, arg[5], false, lmp); rstyle = CONSTANT; } @@ -79,7 +78,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : // error check - if (radius < 0.0) error->all(FLERR,"Illegal region sphere radius: {}", radius); + if (radius < 0.0) error->all(FLERR, "Illegal region sphere radius: {}", radius); // extent of sphere // for variable radius, uses initial radius and origin for variable center @@ -92,7 +91,8 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : extent_yhi = yc + radius; extent_zlo = zc - radius; extent_zhi = zc + radius; - } else bboxflag = 0; + } else + bboxflag = 0; cmax = 1; contact = new Contact[cmax]; @@ -103,11 +103,11 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::~RegSphere() { - delete [] xstr; - delete [] ystr; - delete [] zstr; - delete [] rstr; - delete [] contact; + delete[] xstr; + delete[] ystr; + delete[] zstr; + delete[] rstr; + delete[] contact; } /* ---------------------------------------------------------------------- */ @@ -128,7 +128,7 @@ int RegSphere::inside(double x, double y, double z) double delx = x - xc; double dely = y - yc; double delz = z - zc; - double r = sqrt(delx*delx + dely*dely + delz*delz); + double r = sqrt(delx * delx + dely * dely + delz * delz); if (r <= radius) return 1; return 0; @@ -146,15 +146,15 @@ int RegSphere::surface_interior(double *x, double cutoff) double delx = x[0] - xc; double dely = x[1] - yc; double delz = x[2] - zc; - double r = sqrt(delx*delx + dely*dely + delz*delz); + double r = sqrt(delx * delx + dely * dely + delz * delz); if (r > radius || r == 0.0) return 0; double delta = radius - r; if (delta < cutoff) { contact[0].r = delta; - contact[0].delx = delx*(1.0-radius/r); - contact[0].dely = dely*(1.0-radius/r); - contact[0].delz = delz*(1.0-radius/r); + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); contact[0].radius = -radius; contact[0].iwall = 0; contact[0].varflag = 1; @@ -174,15 +174,15 @@ int RegSphere::surface_exterior(double *x, double cutoff) double delx = x[0] - xc; double dely = x[1] - yc; double delz = x[2] - zc; - double r = sqrt(delx*delx + dely*dely + delz*delz); + double r = sqrt(delx * delx + dely * dely + delz * delz); if (r < radius) return 0; double delta = r - radius; if (delta < cutoff) { contact[0].r = delta; - contact[0].delx = delx*(1.0-radius/r); - contact[0].dely = dely*(1.0-radius/r); - contact[0].delz = delz*(1.0-radius/r); + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); contact[0].radius = radius; contact[0].iwall = 0; contact[0].varflag = 1; @@ -197,19 +197,15 @@ int RegSphere::surface_exterior(double *x, double cutoff) void RegSphere::shape_update() { - if (xstyle == VARIABLE) - xc = xscale * input->variable->compute_equal(xvar); + if (xstyle == VARIABLE) xc = xscale * input->variable->compute_equal(xvar); - if (ystyle == VARIABLE) - yc = yscale * input->variable->compute_equal(yvar); + if (ystyle == VARIABLE) yc = yscale * input->variable->compute_equal(yvar); - if (zstyle == VARIABLE) - zc = zscale * input->variable->compute_equal(zvar); + if (zstyle == VARIABLE) zc = zscale * input->variable->compute_equal(zvar); if (rstyle == VARIABLE) { radius = xscale * input->variable->compute_equal(rvar); - if (radius < 0.0) - error->one(FLERR,"Variable evaluation in region gave bad value"); + if (radius < 0.0) error->one(FLERR, "Variable evaluation in region gave bad value"); } } @@ -221,34 +217,30 @@ void RegSphere::variable_check() { if (xstyle == VARIABLE) { xvar = input->variable->find(xstr); - if (xvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (xvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); if (!input->variable->equalstyle(xvar)) - error->all(FLERR,"Variable for region sphere is invalid style"); + error->all(FLERR, "Variable for region sphere is invalid style"); } if (ystyle == VARIABLE) { yvar = input->variable->find(ystr); - if (yvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (yvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); if (!input->variable->equalstyle(yvar)) - error->all(FLERR,"Variable for region sphere is invalid style"); + error->all(FLERR, "Variable for region sphere is invalid style"); } if (zstyle == VARIABLE) { zvar = input->variable->find(zstr); - if (zvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (zvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); if (!input->variable->equalstyle(zvar)) - error->all(FLERR,"Variable for region sphere is invalid style"); + error->all(FLERR, "Variable for region sphere is invalid style"); } if (rstyle == VARIABLE) { rvar = input->variable->find(rstr); - if (rvar < 0) - error->all(FLERR,"Variable name for region sphere does not exist"); + if (rvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); if (!input->variable->equalstyle(rvar)) - error->all(FLERR,"Variable for region sphere is invalid style"); + error->all(FLERR, "Variable for region sphere is invalid style"); } } @@ -265,27 +257,26 @@ void RegSphere::set_velocity_shape() xcenter[1] = yc; xcenter[2] = zc; forward_transform(xcenter[0], xcenter[1], xcenter[2]); - if (update->ntimestep > 0) rprev = prev[4]; - else rprev = radius; + if (update->ntimestep > 0) + rprev = prev[4]; + else + rprev = radius; prev[4] = radius; } - - /* ---------------------------------------------------------------------- add velocity due to shape change to wall velocity ------------------------------------------------------------------------- */ void RegSphere::velocity_contact_shape(double *vwall, double *xc) { - double delx, dely, delz; // Displacement of contact point in x,y,z + double delx, dely, delz; // Displacement of contact point in x,y,z - delx = (xc[0] - xcenter[0])*(1 - rprev/radius); - dely = (xc[1] - xcenter[1])*(1 - rprev/radius); - delz = (xc[2] - xcenter[2])*(1 - rprev/radius); + delx = (xc[0] - xcenter[0]) * (1 - rprev / radius); + dely = (xc[1] - xcenter[1]) * (1 - rprev / radius); + delz = (xc[2] - xcenter[2]) * (1 - rprev / radius); - vwall[0] += delx/update->dt; - vwall[1] += dely/update->dt; - vwall[2] += delz/update->dt; + vwall[0] += delx / update->dt; + vwall[1] += dely / update->dt; + vwall[2] += delz / update->dt; } - From dfed9bf854c907054608a1a067d8ed1cf236058a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Apr 2023 21:49:43 -0400 Subject: [PATCH 03/20] improve error messages --- src/region_block.cpp | 27 ++++++++++++--------------- src/region_cylinder.cpp | 12 ++++++------ src/region_sphere.cpp | 16 ++++++++-------- 3 files changed, 26 insertions(+), 29 deletions(-) diff --git a/src/region_block.cpp b/src/region_block.cpp index 8b50521d7c..e674cfc61c 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -515,44 +515,41 @@ void RegBlock::variable_check() // addition { if (xlostyle == VARIABLE) { xlovar = input->variable->find(xlostr); - if (xlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (xlovar < 0) error->all(FLERR, "Variable {} for region block does not exist", xlostr); if (!input->variable->equalstyle(xlovar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", xlostr); } - if (xhistyle == VARIABLE) { xhivar = input->variable->find(xhistr); - if (xhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (xhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", xhistr); if (!input->variable->equalstyle(xhivar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", xhistr); } if (ylostyle == VARIABLE) { ylovar = input->variable->find(ylostr); - if (ylovar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (ylovar < 0) error->all(FLERR, "Variable {} for region block does not exist", ylostr); if (!input->variable->equalstyle(ylovar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", ylostr); } - if (yhistyle == VARIABLE) { yhivar = input->variable->find(yhistr); - if (yhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (yhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", yhistr); if (!input->variable->equalstyle(yhivar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", yhistr); } if (zlostyle == VARIABLE) { zlovar = input->variable->find(zlostr); - if (zlovar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (zlovar < 0) error->all(FLERR, "Variable {} for region block does not exist", zlostr); if (!input->variable->equalstyle(zlovar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", zlostr); } - if (zhistyle == VARIABLE) { zhivar = input->variable->find(zhistr); - if (zhivar < 0) error->all(FLERR, "Variable name for region block does not exist"); + if (zhivar < 0) error->all(FLERR, "Variable {} for region block does not exist", zhistr); if (!input->variable->equalstyle(zhivar)) - error->all(FLERR, "Variable for region block is invalid style"); + error->all(FLERR, "Variable {} for region block is invalid style", zhistr); } } diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index c6119c0c3b..601f70e66b 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -793,23 +793,23 @@ void RegCylinder::variable_check() { if (c1style == VARIABLE) { c1var = input->variable->find(c1str); - if (c1var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); + if (c1var < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", c1str); if (!input->variable->equalstyle(c1var)) - error->all(FLERR, "Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable {} for region cylinder is invalid style", c1str); } if (c2style == VARIABLE) { c2var = input->variable->find(c2str); - if (c2var < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); + if (c2var < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", c2str); if (!input->variable->equalstyle(c2var)) - error->all(FLERR, "Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable {} for region cylinder is invalid style", c2str); } if (rstyle == VARIABLE) { rvar = input->variable->find(rstr); - if (rvar < 0) error->all(FLERR, "Variable name for region cylinder does not exist"); + if (rvar < 0) error->all(FLERR, "Variable {} for region cylinder does not exist", rstr); if (!input->variable->equalstyle(rvar)) - error->all(FLERR, "Variable for region cylinder is invalid style"); + error->all(FLERR, "Variable {} for region cylinder is invalid style", rstr); } } diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index b67e805ac4..cd20a697d4 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -217,30 +217,30 @@ void RegSphere::variable_check() { if (xstyle == VARIABLE) { xvar = input->variable->find(xstr); - if (xvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); + if (xvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", xstr); if (!input->variable->equalstyle(xvar)) - error->all(FLERR, "Variable for region sphere is invalid style"); + error->all(FLERR, "Variable {} for region sphere is invalid style", xstr); } if (ystyle == VARIABLE) { yvar = input->variable->find(ystr); - if (yvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); + if (yvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", ystr); if (!input->variable->equalstyle(yvar)) - error->all(FLERR, "Variable for region sphere is invalid style"); + error->all(FLERR, "Variable {} for region sphere is invalid style", ystr); } if (zstyle == VARIABLE) { zvar = input->variable->find(zstr); - if (zvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); + if (zvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", zstr); if (!input->variable->equalstyle(zvar)) - error->all(FLERR, "Variable for region sphere is invalid style"); + error->all(FLERR, "Variable {} for region sphere is invalid style", zstr); } if (rstyle == VARIABLE) { rvar = input->variable->find(rstr); - if (rvar < 0) error->all(FLERR, "Variable name for region sphere does not exist"); + if (rvar < 0) error->all(FLERR, "Variable {} for region sphere does not exist", rstr); if (!input->variable->equalstyle(rvar)) - error->all(FLERR, "Variable for region sphere is invalid style"); + error->all(FLERR, "Variable {} for region sphere is invalid style", rstr); } } From e59c9d0f678f6bbca127ab282eefdbc9135c9b84 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 3 Apr 2023 21:49:53 -0400 Subject: [PATCH 04/20] silence compiler warning --- src/region_block.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/region_block.cpp b/src/region_block.cpp index e674cfc61c..6f97bab9c5 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -30,8 +30,8 @@ static constexpr double BIG = 1.0e20; /* ---------------------------------------------------------------------- */ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : - Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), - zlostr(nullptr), zhistr(nullptr) + Region(lmp, narg, arg), xlostr(nullptr), ylostr(nullptr), zlostr(nullptr), xhistr(nullptr), + yhistr(nullptr), zhistr(nullptr) { options(narg - 8, &arg[8]); From 106f029941b951f6dc29f9c12ae686729e66f511 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Apr 2023 06:05:03 -0400 Subject: [PATCH 05/20] improve Linux distribution detection for recent Fedora versions --- cmake/Modules/LAMMPSUtils.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/LAMMPSUtils.cmake b/cmake/Modules/LAMMPSUtils.cmake index 8a118e7a3b..1ceec7e06e 100644 --- a/cmake/Modules/LAMMPSUtils.cmake +++ b/cmake/Modules/LAMMPSUtils.cmake @@ -152,7 +152,7 @@ endfunction(FetchPotentials) # set CMAKE_LINUX_DISTRO and CMAKE_DISTRO_VERSION on Linux if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (EXISTS /etc/os-release)) file(STRINGS /etc/os-release distro REGEX "^NAME=") - string(REGEX REPLACE "NAME=\"?([^\"]*)\"?" "\\1" distro "${distro}") + string(REGEX REPLACE "NAME=\"?([^ ]+).*\"?" "\\1" distro "${distro}") file(STRINGS /etc/os-release disversion REGEX "^VERSION_ID=") string(REGEX REPLACE "VERSION_ID=\"?([^\"]*)\"?" "\\1" disversion "${disversion}") set(CMAKE_LINUX_DISTRO ${distro}) From c9605e1cba5519c9941da30918437c27b4b4fc24 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Apr 2023 06:05:46 -0400 Subject: [PATCH 06/20] update custom linker support also for "mold" (even faster than lld) --- cmake/Modules/Testing.cmake | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/cmake/Modules/Testing.cmake b/cmake/Modules/Testing.cmake index c82114a1dd..5345211178 100644 --- a/cmake/Modules/Testing.cmake +++ b/cmake/Modules/Testing.cmake @@ -18,29 +18,33 @@ if(ENABLE_TESTING) # we need to build and link a LOT of tester executables, so it is worth checking if # a faster linker is available. requires GNU or Clang compiler, newer CMake. - # also only verified with Fedora Linux > 30 and Ubuntu <= 18.04 (Ubuntu 20.04 fails) + # also only verified with Fedora Linux > 30 and Ubuntu 18.04 or 22.04+(Ubuntu 20.04 fails) if((CMAKE_SYSTEM_NAME STREQUAL "Linux") AND (CMAKE_VERSION VERSION_GREATER_EQUAL 3.13) - AND ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") - OR (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))) - if(((CMAKE_LINUX_DISTRO STREQUAL "Ubuntu") AND (CMAKE_DISTRO_VERSION VERSION_LESS_EQUAL 18.04)) + AND ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") OR (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))) + if(((CMAKE_LINUX_DISTRO STREQUAL "Ubuntu") AND + ((CMAKE_DISTRO_VERSION VERSION_LESS_EQUAL 18.04) OR (CMAKE_DISTRO_VERSION VERSION_GREATER_EQUAL 22.04))) OR ((CMAKE_LINUX_DISTRO STREQUAL "Fedora") AND (CMAKE_DISTRO_VERSION VERSION_GREATER 30))) include(CheckCXXCompilerFlag) set(CMAKE_CUSTOM_LINKER_DEFAULT default) + check_cxx_compiler_flag(-fuse-ld=mold HAVE_MOLD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=lld HAVE_LLD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=gold HAVE_GOLD_LINKER_FLAG) check_cxx_compiler_flag(-fuse-ld=bfd HAVE_BFD_LINKER_FLAG) + find_program(HAVE_MOLD_LINKER_BIN ld.mold) find_program(HAVE_LLD_LINKER_BIN lld ld.lld) find_program(HAVE_GOLD_LINKER_BIN ld.gold) find_program(HAVE_BFD_LINKER_BIN ld.bfd) - if(HAVE_LLD_LINKER_FLAG AND HAVE_LLD_LINKER_BIN) + if(HAVE_MOLD_LINKER_FLAG AND HAVE_MOLD_LINKER_BIN) + set(CMAKE_CUSTOM_LINKER_DEFAULT mold) + elseif(HAVE_LLD_LINKER_FLAG AND HAVE_LLD_LINKER_BIN) set(CMAKE_CUSTOM_LINKER_DEFAULT lld) elseif(HAVE_GOLD_LINKER_FLAG AND HAVE_GOLD_LINKER_BIN) set(CMAKE_CUSTOM_LINKER_DEFAULT gold) elseif(HAVE_BFD_LINKER_FLAG AND HAVE_BFD_LINKER_BIN) set(CMAKE_CUSTOM_LINKER_DEFAULT bfd) endif() - set(CMAKE_CUSTOM_LINKER_VALUES lld gold bfd default) - set(CMAKE_CUSTOM_LINKER ${CMAKE_CUSTOM_LINKER_DEFAULT} CACHE STRING "Choose a custom linker for faster linking (lld, gold, bfd, default)") + set(CMAKE_CUSTOM_LINKER_VALUES mold lld gold bfd default) + set(CMAKE_CUSTOM_LINKER ${CMAKE_CUSTOM_LINKER_DEFAULT} CACHE STRING "Choose a custom linker for faster linking (mold, lld, gold, bfd, default)") validate_option(CMAKE_CUSTOM_LINKER CMAKE_CUSTOM_LINKER_VALUES) mark_as_advanced(CMAKE_CUSTOM_LINKER) if(NOT "${CMAKE_CUSTOM_LINKER}" STREQUAL "default") From e670a94b8a9dbddffec0fae5c2380ee6f0299297 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Apr 2023 08:09:47 -0400 Subject: [PATCH 07/20] cosmetic --- src/MANYBODY/pair_bop.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp index f8c119c436..518a6dad79 100644 --- a/src/MANYBODY/pair_bop.cpp +++ b/src/MANYBODY/pair_bop.cpp @@ -1873,8 +1873,7 @@ void PairBOP::read_table(char *filename) reader = new PotentialFileReader(lmp, filename, "BOP"); bop_types = reader->next_int(); if (bop_types <= 0) - error->one(FLERR,fmt::format("BOP potential file with {} " - "elements",bop_types)); + error->one(FLERR,fmt::format("BOP potential file with {} elements",bop_types)); bop_elements = new char*[bop_types]; bop_masses = new double[bop_types]; From 4e48ddb9757e9c72765b8d510866e371cfdca9bb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 6 Apr 2023 07:21:37 -0400 Subject: [PATCH 08/20] No need to use nvcc_wrapper globally when configuring with CMake --- doc/src/Build_extras.rst | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index a2125244ba..0ecf54f744 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -683,20 +683,11 @@ This list was last updated for version 3.7.1 of the Kokkos library. -D Kokkos_ARCH_GPUARCH=yes # GPUARCH = GPU from list above -D Kokkos_ENABLE_CUDA=yes -D Kokkos_ENABLE_OPENMP=yes - -D CMAKE_CXX_COMPILER=wrapper # wrapper = full path to Cuda nvcc wrapper This will also enable executing FFTs on the GPU, either via the internal KISSFFT library, or - by preference - with the cuFFT library bundled with the CUDA toolkit, depending on whether CMake - can identify its location. The *wrapper* value for - ``CMAKE_CXX_COMPILER`` variable is the path to the CUDA nvcc - compiler wrapper provided in the Kokkos library: - ``lib/kokkos/bin/nvcc_wrapper``\ . The setting should include the - full path name to the wrapper, e.g. - - .. code-block:: bash - - -D CMAKE_CXX_COMPILER=${HOME}/lammps/lib/kokkos/bin/nvcc_wrapper + can identify its location. For AMD or NVIDIA GPUs using HIP, set these variables: From e44aa7740379e02bdefe78c27360cbb81caaadde Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 6 Apr 2023 07:23:49 -0400 Subject: [PATCH 09/20] fix copy-n-paste error --- doc/src/fix_wall_gran.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 0006e1750b..5435b6d5b8 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -71,9 +71,9 @@ Examples fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 - fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox - fix 4 all wall/gran/region granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone - fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0 + fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox + fix 4 all wall/gran granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone + fix 5 all wall/gran granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0 fix 6 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts Description From ebcb4432376b2458c8e5c958358d0e6d76dd6e35 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 6 Apr 2023 11:42:19 -0400 Subject: [PATCH 10/20] correct link --- doc/src/pair_colloid.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_colloid.rst b/doc/src/pair_colloid.rst index 1e75ae1c78..26e7fa7da4 100644 --- a/doc/src/pair_colloid.rst +++ b/doc/src/pair_colloid.rst @@ -93,7 +93,7 @@ with :math:`A_{ss}` set appropriately, which results from letting both particle sizes go to zero. When used in combination with :doc:`pair_style yukawa/colloid -`, the two terms become the so-called DLVO potential, +`, the two terms become the so-called DLVO potential, which combines electrostatic repulsion and van der Waals attraction. The following coefficients must be defined for each pair of atoms From 4c403e5b71d23366ed49d2e97477b229f0f3097e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 6 Apr 2023 12:58:30 -0400 Subject: [PATCH 11/20] close and finalize LAMMPS instance in MLIAP example python scripts --- examples/mliap/mliap_numpy_Ta06A.py | 6 ++-- examples/mliap/mliap_pytorch_Ta06A.py | 4 ++- examples/mliap/mliap_pytorch_Ta06A_kokkos.py | 4 ++- examples/mliap/mliap_unified_lj_Ar.py | 32 +++++++++++--------- 4 files changed, 26 insertions(+), 20 deletions(-) diff --git a/examples/mliap/mliap_numpy_Ta06A.py b/examples/mliap/mliap_numpy_Ta06A.py index 25c928cd89..180d221d35 100644 --- a/examples/mliap/mliap_numpy_Ta06A.py +++ b/examples/mliap/mliap_numpy_Ta06A.py @@ -85,8 +85,7 @@ class LinearModel(): def __call__(self,elems,bispectrum,beta,energy): energy[:] = bispectrum @ self.weights + self.bias beta[:] = self.weights - - + mymodel = LinearModel("Ta06A.mliap.model") import lammps @@ -98,4 +97,5 @@ lammps.mliap.activate_mliappy(lmp) lmp.commands_string(before_loading) lammps.mliap.load_model(mymodel) lmp.commands_string(after_loading) - +lmp.close() +lmp.finalize() diff --git a/examples/mliap/mliap_pytorch_Ta06A.py b/examples/mliap/mliap_pytorch_Ta06A.py index 589c9661b8..7c8949d8d6 100644 --- a/examples/mliap/mliap_pytorch_Ta06A.py +++ b/examples/mliap/mliap_pytorch_Ta06A.py @@ -103,6 +103,8 @@ model = torch.load(torch_model) # Connect the PyTorch model to the mliap pair style. lammps.mliap.load_model(model) - + # run the simulation with the mliap pair style lmp.commands_string(after_loading) +lmp.close() +lmp.finalize() diff --git a/examples/mliap/mliap_pytorch_Ta06A_kokkos.py b/examples/mliap/mliap_pytorch_Ta06A_kokkos.py index f0d2b9bf3e..8afa3f0df3 100644 --- a/examples/mliap/mliap_pytorch_Ta06A_kokkos.py +++ b/examples/mliap/mliap_pytorch_Ta06A_kokkos.py @@ -103,6 +103,8 @@ model = torch.load(torch_model) # Connect the PyTorch model to the mliap pair style. lammps.mliap.load_model_kokkos(model) - + # run the simulation with the mliap pair style lmp.commands_string(after_loading) +lmp.close() +lmp.finalize() diff --git a/examples/mliap/mliap_unified_lj_Ar.py b/examples/mliap/mliap_unified_lj_Ar.py index 2b404bd98e..d48bdd385a 100644 --- a/examples/mliap/mliap_unified_lj_Ar.py +++ b/examples/mliap/mliap_unified_lj_Ar.py @@ -6,30 +6,30 @@ before_loading =\ """# 3d Lennard-Jones melt -units lj -atom_style atomic +units lj +atom_style atomic -lattice fcc 0.8442 -region box block 0 10 0 10 0 10 -create_box 1 box -create_atoms 1 box -mass 1 1.0 +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 -velocity all create 3.0 87287 loop geom +velocity all create 3.0 87287 loop geom """ after_loading =\ """ -pair_style mliap unified EXISTS -pair_coeff * * Ar +pair_style mliap unified EXISTS +pair_coeff * * Ar -neighbor 0.3 bin -neigh_modify every 20 delay 0 check no +neighbor 0.3 bin +neigh_modify every 20 delay 0 check no -fix 1 all nve +fix 1 all nve -thermo 50 -run 250 +thermo 50 +run 250 """ import lammps @@ -63,3 +63,5 @@ lammps.mliap.load_unified(unified) # Run the simulation with the mliap unified pair style # Use pre-loaded model by specifying model filename as "EXISTS" lmp.commands_string(after_loading) +lmp.close() +lmp.finalize() From a98a77041ef3cd1c677e5bbea28d176df7e96c45 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 15:05:46 -0400 Subject: [PATCH 12/20] add option to include PLUMED into cross-compiled Windows binaries --- cmake/Modules/Packages/PLUMED.cmake | 196 +++++++++++++++------------- 1 file changed, 107 insertions(+), 89 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index 9a4a9556ee..22d992bcaf 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -1,106 +1,124 @@ -set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library") -set(PLUMED_MODE_VALUES static shared runtime) -set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES}) -validate_option(PLUMED_MODE PLUMED_MODE_VALUES) -string(TOUPPER ${PLUMED_MODE} PLUMED_MODE) +# Plumed2 support for PLUMED package -set(PLUMED_LINK_LIBS) -if(PLUMED_MODE STREQUAL "STATIC") - find_package(LAPACK REQUIRED) - find_package(BLAS REQUIRED) - find_package(GSL REQUIRED) - list(APPEND PLUMED_LINK_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} GSL::gsl) - find_package(ZLIB QUIET) - if(ZLIB_FOUND) - list(APPEND PLUMED_LINK_LIBS ZLIB::ZLIB) - endif() - find_package(FFTW3 QUIET) - if(FFTW3_FOUND) - list(APPEND PLUMED_LINK_LIBS FFTW3::FFTW3) - endif() -endif() +set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") +set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") -find_package(PkgConfig QUIET) -set(DOWNLOAD_PLUMED_DEFAULT ON) -if(PKG_CONFIG_FOUND) - pkg_check_modules(PLUMED QUIET plumed) - if(PLUMED_FOUND) - set(DOWNLOAD_PLUMED_DEFAULT OFF) - endif() -endif() +mark_as_advanced(PLUMED_URL) +mark_as_advanced(PLUMED_MD5) +GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) -option(DOWNLOAD_PLUMED "Download Plumed package instead of using an already installed one" ${DOWNLOAD_PLUMED_DEFAULT}) -if(DOWNLOAD_PLUMED) - if(BUILD_MPI) - set(PLUMED_CONFIG_MPI "--enable-mpi") - set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER}) - set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER}) - else() - set(PLUMED_CONFIG_MPI "--disable-mpi") - set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER}) - set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER}) +if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) + if(NOT PLUMED_BUILD_DIR) + message(FATAL_ERROR "Must set PLUMED_BUILD_DIR when cross-compiling for Windows") endif() - if(BUILD_OMP) - set(PLUMED_CONFIG_OMP "--enable-openmp") - else() - set(PLUMED_CONFIG_OMP "--disable-openmp") - endif() - message(STATUS "PLUMED download requested - we will build our own") + add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) + set_target_properties(LAMMPS::PLUMED PROPERTIES + IMPORTED_LOCATION ${PLUMED_BUILD_DIR}/lib/install/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} + INTERFACE_LINK_LIBRARIES "-Wl,--image-base -Wl,0x10000000 -lfftw3 -lz -fstack-protector -lssp -fopenmp" + INTERFACE_INCLUDE_DIRECTORIES ${PLUMED_BUILD_DIR}/include) + target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) + target_compile_options(lammps PRIVATE -mcmodel=medium) +else() + set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library") + set(PLUMED_MODE_VALUES static shared runtime) + set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES}) + validate_option(PLUMED_MODE PLUMED_MODE_VALUES) + string(TOUPPER ${PLUMED_MODE} PLUMED_MODE) + + set(PLUMED_LINK_LIBS) if(PLUMED_MODE STREQUAL "STATIC") - set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}") - elseif(PLUMED_MODE STREQUAL "SHARED") - set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX};/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") - elseif(PLUMED_MODE STREQUAL "RUNTIME") - set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_PREFIX}") + find_package(LAPACK REQUIRED) + find_package(BLAS REQUIRED) + find_package(GSL REQUIRED) + list(APPEND PLUMED_LINK_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} GSL::gsl) + find_package(ZLIB QUIET) + if(ZLIB_FOUND) + list(APPEND PLUMED_LINK_LIBS ZLIB::ZLIB) + endif() + find_package(FFTW3 QUIET) + if(FFTW3_FOUND) + list(APPEND PLUMED_LINK_LIBS FFTW3::FFTW3) + endif() endif() - set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") - set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") + find_package(PkgConfig QUIET) + set(DOWNLOAD_PLUMED_DEFAULT ON) + if(PKG_CONFIG_FOUND) + pkg_check_modules(PLUMED QUIET plumed) + if(PLUMED_FOUND) + set(DOWNLOAD_PLUMED_DEFAULT OFF) + endif() + endif() - mark_as_advanced(PLUMED_URL) - mark_as_advanced(PLUMED_MD5) - GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) + option(DOWNLOAD_PLUMED "Download Plumed package instead of using an already installed one" ${DOWNLOAD_PLUMED_DEFAULT}) + if(DOWNLOAD_PLUMED) + if(BUILD_MPI) + set(PLUMED_CONFIG_MPI "--enable-mpi") + set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER}) + set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER}) + else() + set(PLUMED_CONFIG_MPI "--disable-mpi") + set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER}) + set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER}) + endif() + if(BUILD_OMP) + set(PLUMED_CONFIG_OMP "--enable-openmp") + else() + set(PLUMED_CONFIG_OMP "--disable-openmp") + endif() + message(STATUS "PLUMED download requested - we will build our own") + if(PLUMED_MODE STREQUAL "STATIC") + set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}") + elseif(PLUMED_MODE STREQUAL "SHARED") + set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX};/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") + elseif(PLUMED_MODE STREQUAL "RUNTIME") + set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_PREFIX}") + endif() - include(ExternalProject) - ExternalProject_Add(plumed_build - URL ${PLUMED_URL} ${PLUMED_FALLBACK} - URL_MD5 ${PLUMED_MD5} - BUILD_IN_SOURCE 1 - CONFIGURE_COMMAND /configure --prefix= + include(ExternalProject) + ExternalProject_Add(plumed_build + URL ${PLUMED_URL} ${PLUMED_FALLBACK} + URL_MD5 ${PLUMED_MD5} + BUILD_IN_SOURCE 1 + CONFIGURE_COMMAND /configure --prefix= ${CONFIGURE_REQUEST_PIC} --enable-modules=all + --enable-cxx=11 + --disable-python + --disable-doc ${PLUMED_CONFIG_MPI} ${PLUMED_CONFIG_OMP} CXX=${PLUMED_CONFIG_CXX} CC=${PLUMED_CONFIG_CC} - BUILD_BYPRODUCTS ${PLUMED_BUILD_BYPRODUCTS} - ) - ExternalProject_get_property(plumed_build INSTALL_DIR) - add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) - add_dependencies(LAMMPS::PLUMED plumed_build) - if(PLUMED_MODE STREQUAL "STATIC") - set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${PLUMED_LINK_LIBS};${CMAKE_DL_LIBS}") - elseif(PLUMED_MODE STREQUAL "SHARED") - set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX};${CMAKE_DL_LIBS}") - elseif(PLUMED_MODE STREQUAL "RUNTIME") - set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") - set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${CMAKE_DL_LIBS}") + BUILD_BYPRODUCTS ${PLUMED_BUILD_BYPRODUCTS} + ) + ExternalProject_get_property(plumed_build INSTALL_DIR) + add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) + add_dependencies(LAMMPS::PLUMED plumed_build) + if(PLUMED_MODE STREQUAL "STATIC") + set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${PLUMED_LINK_LIBS};${CMAKE_DL_LIBS}") + elseif(PLUMED_MODE STREQUAL "SHARED") + set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${CMAKE_SHARED_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX};${CMAKE_DL_LIBS}") + elseif(PLUMED_MODE STREQUAL "RUNTIME") + set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${INSTALL_DIR}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") + set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION ${INSTALL_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_SUFFIX} INTERFACE_LINK_LIBRARIES "${CMAKE_DL_LIBS}") + endif() + set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include) + file(MAKE_DIRECTORY ${INSTALL_DIR}/include) + else() + find_package(PkgConfig REQUIRED) + pkg_check_modules(PLUMED REQUIRED plumed) + add_library(LAMMPS::PLUMED INTERFACE IMPORTED) + if(PLUMED_MODE STREQUAL "STATIC") + include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static) + elseif(PLUMED_MODE STREQUAL "SHARED") + include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared) + elseif(PLUMED_MODE STREQUAL "RUNTIME") + set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") + include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime) + endif() + set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}") + set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}") endif() - set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES ${INSTALL_DIR}/include) - file(MAKE_DIRECTORY ${INSTALL_DIR}/include) -else() - find_package(PkgConfig REQUIRED) - pkg_check_modules(PLUMED REQUIRED plumed) - add_library(LAMMPS::PLUMED INTERFACE IMPORTED) - if(PLUMED_MODE STREQUAL "STATIC") - include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static) - elseif(PLUMED_MODE STREQUAL "SHARED") - include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared) - elseif(PLUMED_MODE STREQUAL "RUNTIME") - set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}") - include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime) - endif() - set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}") - set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}") + target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) endif() -target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) From a6a8f2c451902230c5db2ca968efb387e457fe5f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 15:28:18 -0400 Subject: [PATCH 13/20] more tweaks for cross-compiling plumed for windows --- cmake/Modules/Packages/PLUMED.cmake | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index 22d992bcaf..337c475d20 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -1,24 +1,31 @@ # Plumed2 support for PLUMED package -set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") -set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") - -mark_as_advanced(PLUMED_URL) -mark_as_advanced(PLUMED_MD5) -GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) - if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) + # special case for cross-compiling to windows with externally cross-compiled plumed tree if(NOT PLUMED_BUILD_DIR) message(FATAL_ERROR "Must set PLUMED_BUILD_DIR when cross-compiling for Windows") + else() + set(PLUMED_INSTALL_DIR "${PLUMED_BUILD_DIR/src/lib/install") endif() add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) set_target_properties(LAMMPS::PLUMED PROPERTIES - IMPORTED_LOCATION ${PLUMED_BUILD_DIR}/lib/install/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX} + IMPORTED_LOCATION "${PLUMED_INSTALL_DIR}/libplumed.a" INTERFACE_LINK_LIBRARIES "-Wl,--image-base -Wl,0x10000000 -lfftw3 -lz -fstack-protector -lssp -fopenmp" - INTERFACE_INCLUDE_DIRECTORIES ${PLUMED_BUILD_DIR}/include) + INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_BUILD_DIR}/src/include") target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) - target_compile_options(lammps PRIVATE -mcmodel=medium) + add_custom_command(OUTPUT ${CMAKE_BINARY_DIR}/plumed.exe + COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PLUMED_INSTALL_DIR}/plumed.exe) + add_custom_command(OUTPUT ${CMAKE_BINARY_DIR}/plumed_patches + COMMAND ${CMAKE_COMMAND} -E copy_directory_if_different ${PLUMED_BUILD_DIR}/patches) else() + + set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") + set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") + + mark_as_advanced(PLUMED_URL) + mark_as_advanced(PLUMED_MD5) + GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) + set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library") set(PLUMED_MODE_VALUES static shared runtime) set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES}) From d5680f0c6d63b8c1c283f7cfb5e74cbb3354e691 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 15:44:56 -0400 Subject: [PATCH 14/20] fix typo --- cmake/Modules/Packages/PLUMED.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index 337c475d20..68a0a390a6 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -5,7 +5,7 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) if(NOT PLUMED_BUILD_DIR) message(FATAL_ERROR "Must set PLUMED_BUILD_DIR when cross-compiling for Windows") else() - set(PLUMED_INSTALL_DIR "${PLUMED_BUILD_DIR/src/lib/install") + set(PLUMED_INSTALL_DIR "${PLUMED_BUILD_DIR}/src/lib/install") endif() add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) set_target_properties(LAMMPS::PLUMED PROPERTIES From d71de7cc387395415656114a61d423346986ea9c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 16:36:52 -0400 Subject: [PATCH 15/20] handle exceptions in destructors --- src/MC/fix_charge_regulation.cpp | 7 ++++++- src/MC/fix_gcmc.cpp | 14 ++++++++++++-- src/MC/fix_widom.cpp | 14 ++++++++++++-- 3 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/MC/fix_charge_regulation.cpp b/src/MC/fix_charge_regulation.cpp index 0626762faf..a828f276ea 100644 --- a/src/MC/fix_charge_regulation.cpp +++ b/src/MC/fix_charge_regulation.cpp @@ -160,7 +160,12 @@ FixChargeRegulation::~FixChargeRegulation() { if (exclusion_group_bit && group) { auto group_id = std::string("FixChargeRegulation:gcmc_exclusion_group:") + id; - group->assign(group_id + " delete"); + try { + group->assign(group_id + " delete"); + } catch (std::exception &e) { + if (comm->me == 0) + fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what()); + } } if (group) { diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index f8193f6235..634b512936 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -427,12 +427,22 @@ FixGCMC::~FixGCMC() if (exclusion_group_bit && group) { auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id; - group->assign(group_id + " delete"); + try { + group->assign(group_id + " delete"); + } catch (std::exception &e) { + if (comm->me == 0) + fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what()); + } } if (molecule_group_bit && group) { auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id; - group->assign(group_id + " delete"); + try { + group->assign(group_id + " delete"); + } catch (std::exception &e) { + if (comm->me == 0) + fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what()); + } } if (full_flag && group && neighbor) { diff --git a/src/MC/fix_widom.cpp b/src/MC/fix_widom.cpp index d1f2b95740..a98f29da5e 100644 --- a/src/MC/fix_widom.cpp +++ b/src/MC/fix_widom.cpp @@ -255,12 +255,22 @@ FixWidom::~FixWidom() if (exclusion_group_bit && group) { auto group_id = std::string("FixWidom:widom_exclusion_group:") + id; - group->assign(group_id + " delete"); + try { + group->assign(group_id + " delete"); + } catch (std::exception &e) { + if (comm->me == 0) + fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what()); + } } if (molecule_group_bit && group) { auto group_id = std::string("FixWidom:rotation_gas_atoms:") + id; - group->assign(group_id + " delete"); + try { + group->assign(group_id + " delete"); + } catch (std::exception &e) { + if (comm->me == 0) + fprintf(stderr, "Error deleting group %s: %s\n", group_id.c_str(), e.what()); + } } if (full_flag && group && neighbor) { From 891c284754ed45c6ab94c67a5c86cf0848a434d0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 16:38:39 -0400 Subject: [PATCH 16/20] avoid static code analysis warnings --- src/EXTRA-PAIR/pair_born_gauss.cpp | 2 ++ src/region_block.cpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/EXTRA-PAIR/pair_born_gauss.cpp b/src/EXTRA-PAIR/pair_born_gauss.cpp index cf64bf2ef1..f60cc4dc6f 100644 --- a/src/EXTRA-PAIR/pair_born_gauss.cpp +++ b/src/EXTRA-PAIR/pair_born_gauss.cpp @@ -38,6 +38,8 @@ PairBornGauss::PairBornGauss(LAMMPS *lmp) : single_enable = 1; respa_enable = 0; writedata = 1; + + cut_global = 0.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/region_block.cpp b/src/region_block.cpp index 6f97bab9c5..efa3d8ca6a 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -33,6 +33,8 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg), xlostr(nullptr), ylostr(nullptr), zlostr(nullptr), xhistr(nullptr), yhistr(nullptr), zhistr(nullptr) { + xlovar = xhivar = ylovar = yhivar = zlovar = zhivar = -1; + options(narg - 8, &arg[8]); xlostyle = CONSTANT; From 71700b8765cc7ae0e9bd2c990b089713745fef35 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 8 Apr 2023 18:07:03 -0400 Subject: [PATCH 17/20] use explicit target to update/copy precompiled plumed files --- cmake/Modules/Packages/PLUMED.cmake | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index 68a0a390a6..a5592cff70 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -13,10 +13,14 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) INTERFACE_LINK_LIBRARIES "-Wl,--image-base -Wl,0x10000000 -lfftw3 -lz -fstack-protector -lssp -fopenmp" INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_BUILD_DIR}/src/include") target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) - add_custom_command(OUTPUT ${CMAKE_BINARY_DIR}/plumed.exe - COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PLUMED_INSTALL_DIR}/plumed.exe) - add_custom_command(OUTPUT ${CMAKE_BINARY_DIR}/plumed_patches - COMMAND ${CMAKE_COMMAND} -E copy_directory_if_different ${PLUMED_BUILD_DIR}/patches) + + add_custom_target(copy_plumed ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches + COMMAND ${CMAKE_COMMAND} -E copy ${PLUMED_INSTALL_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed.exe + COMMAND ${CMAKE_COMMAND} -E copy_directory ${PLUMED_BUILD_DIR}/patches ${CMAKE_BINARY_DIR}/plumed_patches + BYPRODUCTS ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches + COMMENT "Copying Plumed files" + ) + else() set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") From f84a31dfdf12ca4852ee397dea9ff2a7f768a1c4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 9 Apr 2023 01:46:58 -0400 Subject: [PATCH 18/20] change folder for patches, so we can set PLUMED_ROOT accordigly --- cmake/Modules/Packages/PLUMED.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index a5592cff70..1d5b7b9d02 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -16,8 +16,8 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) add_custom_target(copy_plumed ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches COMMAND ${CMAKE_COMMAND} -E copy ${PLUMED_INSTALL_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed.exe - COMMAND ${CMAKE_COMMAND} -E copy_directory ${PLUMED_BUILD_DIR}/patches ${CMAKE_BINARY_DIR}/plumed_patches - BYPRODUCTS ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches + COMMAND ${CMAKE_COMMAND} -E copy_directory ${PLUMED_BUILD_DIR}/patches ${CMAKE_BINARY_DIR}/patches + BYPRODUCTS ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/patches COMMENT "Copying Plumed files" ) From bc4d664f2b6b009d4f668a5b827f90a2c031beb0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 10 Apr 2023 07:14:07 -0400 Subject: [PATCH 19/20] Fully integrate cross-compiling Plumed2 lib into CMake build system --- cmake/Modules/Packages/PLUMED.cmake | 87 ++++++++++++++++++++--------- 1 file changed, 60 insertions(+), 27 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index 1d5b7b9d02..bb31da6bdf 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -1,12 +1,66 @@ # Plumed2 support for PLUMED package +if(BUILD_MPI) + set(PLUMED_CONFIG_MPI "--enable-mpi") + set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER}) + set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER}) + set(PLUMED_CONFIG_CPP "-I ${MPI_CXX_INCLUDE_PATH}") + set(PLUMED_CONFIG_LIB "${MPI_CXX_LIBRARIES}") + set(PLUMED_CONFIG_DEP "mpi4win_build") +else() + set(PLUMED_CONFIG_MPI "--disable-mpi") + set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER}) + set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER}) + set(PLUMED_CONFIG_CPP "") + set(PLUMED_CONFIG_LIB "") + set(PLUMED_CONFIG_DEP "") +endif() +if(BUILD_OMP) + set(PLUMED_CONFIG_OMP "--enable-openmp") +else() + set(PLUMED_CONFIG_OMP "--disable-openmp") +endif() + +set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" + CACHE STRING "URL for PLUMED tarball") +set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") + +mark_as_advanced(PLUMED_URL) +mark_as_advanced(PLUMED_MD5) +GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) + if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) - # special case for cross-compiling to windows with externally cross-compiled plumed tree - if(NOT PLUMED_BUILD_DIR) - message(FATAL_ERROR "Must set PLUMED_BUILD_DIR when cross-compiling for Windows") + if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64") + set(CROSS_CONFIGURE mingw64-configure) + elseif(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86") + set(CROSS_CONFIGURE mingw32-configure) else() - set(PLUMED_INSTALL_DIR "${PLUMED_BUILD_DIR}/src/lib/install") + message(FATAL_ERROR "Unsupported target system: ${CMAKE_SYSTEM_NAME}/${CMAKE_SYSTEM_PROCESSOR}") endif() + message(STATUS "Downloading and cross-compiling Plumed2 for ${CMAKE_SYSTEM_NAME}/${CMAKE_SYSTEM_PROCESSOR} with ${CROSS_CONFIGURE}") + include(ExternalProject) + ExternalProject_Add(plumed_build + URL ${PLUMED_URL} ${PLUMED_FALLBACK} + URL_MD5 ${PLUMED_MD5} + BUILD_IN_SOURCE 1 + CONFIGURE_COMMAND ${CROSS_CONFIGURE} --disable-shared --disable-bsymbolic + --disable-python --enable-cxx=11 + --enable-modules=-adjmat:+crystallization:-dimred:+drr:+eds:-fisst:+funnel:+logmfd:+manyrestraints:+maze:+opes:+multicolvar:-pamm:-piv:+s2cm:-sasa:-ves + ${PLUMED_CONFIG_OMP} + ${PLUMED_CONFIG_MPI} + CXX=${PLUMED_CONFIG_CXX} + CC=${PLUMED_CONFIG_CC} + CPPFLAGS=${PLUMED_CONFIG_CPP} + LIBS=${PLUMED_CONFIG_LIB} + INSTALL_COMMAND "" + BUILD_BYPRODUCTS "/src/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}" "/src/lib/install/plumed{EXE_SUFFIX}" + DEPENDS "${PLUMED_MPI_CONFIG_DEP}" + ) + ExternalProject_Get_Property(plumed_build SOURCE_DIR) + set(PLUMED_BUILD_DIR ${SOURCE_DIR}) + set(PLUMED_INSTALL_DIR ${PLUMED_BUILD_DIR}/src/lib/install) + file(MAKE_DIRECTORY ${PLUMED_BUILD_DIR}/src/include) + add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION "${PLUMED_INSTALL_DIR}/libplumed.a" @@ -14,22 +68,16 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_BUILD_DIR}/src/include") target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) - add_custom_target(copy_plumed ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches + add_custom_target(plumed_copy ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches COMMAND ${CMAKE_COMMAND} -E copy ${PLUMED_INSTALL_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed.exe COMMAND ${CMAKE_COMMAND} -E copy_directory ${PLUMED_BUILD_DIR}/patches ${CMAKE_BINARY_DIR}/patches BYPRODUCTS ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/patches + DEPENDS plumed_build COMMENT "Copying Plumed files" ) else() - set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball") - set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball") - - mark_as_advanced(PLUMED_URL) - mark_as_advanced(PLUMED_MD5) - GetFallbackURL(PLUMED_URL PLUMED_FALLBACK) - set(PLUMED_MODE "static" CACHE STRING "Linkage mode for Plumed2 library") set(PLUMED_MODE_VALUES static shared runtime) set_property(CACHE PLUMED_MODE PROPERTY STRINGS ${PLUMED_MODE_VALUES}) @@ -63,20 +111,6 @@ else() option(DOWNLOAD_PLUMED "Download Plumed package instead of using an already installed one" ${DOWNLOAD_PLUMED_DEFAULT}) if(DOWNLOAD_PLUMED) - if(BUILD_MPI) - set(PLUMED_CONFIG_MPI "--enable-mpi") - set(PLUMED_CONFIG_CC ${CMAKE_MPI_C_COMPILER}) - set(PLUMED_CONFIG_CXX ${CMAKE_MPI_CXX_COMPILER}) - else() - set(PLUMED_CONFIG_MPI "--disable-mpi") - set(PLUMED_CONFIG_CC ${CMAKE_C_COMPILER}) - set(PLUMED_CONFIG_CXX ${CMAKE_CXX_COMPILER}) - endif() - if(BUILD_OMP) - set(PLUMED_CONFIG_OMP "--enable-openmp") - else() - set(PLUMED_CONFIG_OMP "--disable-openmp") - endif() message(STATUS "PLUMED download requested - we will build our own") if(PLUMED_MODE STREQUAL "STATIC") set(PLUMED_BUILD_BYPRODUCTS "/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}") @@ -96,7 +130,6 @@ else() --enable-modules=all --enable-cxx=11 --disable-python - --disable-doc ${PLUMED_CONFIG_MPI} ${PLUMED_CONFIG_OMP} CXX=${PLUMED_CONFIG_CXX} From 415be03f6c2de7d7344ac6ebc59e83c062a9e74d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 10 Apr 2023 09:16:31 -0400 Subject: [PATCH 20/20] improve dependency processing --- cmake/Modules/Packages/PLUMED.cmake | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cmake/Modules/Packages/PLUMED.cmake b/cmake/Modules/Packages/PLUMED.cmake index bb31da6bdf..b90bff4b9c 100644 --- a/cmake/Modules/Packages/PLUMED.cmake +++ b/cmake/Modules/Packages/PLUMED.cmake @@ -53,7 +53,7 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) CPPFLAGS=${PLUMED_CONFIG_CPP} LIBS=${PLUMED_CONFIG_LIB} INSTALL_COMMAND "" - BUILD_BYPRODUCTS "/src/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumed${CMAKE_STATIC_LIBRARY_SUFFIX}" "/src/lib/install/plumed{EXE_SUFFIX}" + BUILD_BYPRODUCTS "/src/lib/install/libplumed.a" "/src/lib/install/plumed.exe" DEPENDS "${PLUMED_MPI_CONFIG_DEP}" ) ExternalProject_Get_Property(plumed_build SOURCE_DIR) @@ -62,11 +62,11 @@ if((CMAKE_SYSTEM_NAME STREQUAL "Windows") AND (CMAKE_CROSSCOMPILING)) file(MAKE_DIRECTORY ${PLUMED_BUILD_DIR}/src/include) add_library(LAMMPS::PLUMED UNKNOWN IMPORTED) + add_dependencies(LAMMPS::PLUMED plumed_build) set_target_properties(LAMMPS::PLUMED PROPERTIES IMPORTED_LOCATION "${PLUMED_INSTALL_DIR}/libplumed.a" INTERFACE_LINK_LIBRARIES "-Wl,--image-base -Wl,0x10000000 -lfftw3 -lz -fstack-protector -lssp -fopenmp" INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_BUILD_DIR}/src/include") - target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) add_custom_target(plumed_copy ALL ${CMAKE_COMMAND} -E rm -rf ${CMAKE_BINARY_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed_patches COMMAND ${CMAKE_COMMAND} -E copy ${PLUMED_INSTALL_DIR}/plumed.exe ${CMAKE_BINARY_DIR}/plumed.exe @@ -164,5 +164,6 @@ else() set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}") set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}") endif() - target_link_libraries(lammps PRIVATE LAMMPS::PLUMED) endif() + +target_link_libraries(lammps PRIVATE LAMMPS::PLUMED)