From 7c30111fa06e88bad0381c748d169ce4f97c497c Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 5 May 2025 11:50:59 -0600 Subject: [PATCH 1/7] Fixing passed arg from type to index --- src/RHEO/fix_rheo_pressure.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 121473f7e3..6088c44348 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -191,14 +191,13 @@ void FixRHEOPressure::setup_pre_force(int /*vflag*/) void FixRHEOPressure::pre_force(int /*vflag*/) { int *mask = atom->mask; - int *type = atom->type; double *rho = atom->rho; double *pressure = atom->pressure; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) pressure[i] = calc_pressure(rho[i], type[i]); + if (mask[i] & groupbit) pressure[i] = calc_pressure(rho[i], i); if (comm_forward) comm->forward_comm(this); } From dc07a1471ec938aaa62f452942fd46c3adf76a3b Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 5 May 2025 12:45:35 -0600 Subject: [PATCH 2/7] Properly processing >1 fix set force --- src/RHEO/compute_rheo_interface.cpp | 41 +++++++++++------------------ 1 file changed, 16 insertions(+), 25 deletions(-) diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 0cd6e46417..c2440100ae 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -344,31 +344,22 @@ void ComputeRHEOInterface::store_forces() // If forces are overwritten by a fix, there are no pressure forces // so just normalize auto fixlist = modify->get_fix_by_style("setforce"); - if (fixlist.size() != 0) { - for (const auto &fix : fixlist) { - for (int i = 0; i < atom->nlocal; i++) { - if (rmass) - minv = 1.0 / rmass[i]; - else - minv = 1.0 / mass[type[i]]; - if (mask[i] & fix->groupbit) - for (int a = 0; a < 3; a++) fp_store[i][a] = f[i][a] * minv; - else - for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; - } - } - } else { - if (rmass) { - for (int i = 0; i < atom->nlocal; i++) { - minv = 1.0 / rmass[i]; - for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; - } - } else { - for (int i = 0; i < atom->nlocal; i++) { - minv = 1.0 / mass[type[i]]; - for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; - } - } + int skip_from_setforce; + for (int i = 0; i < atom->nlocal; i++) { + if (rmass) + minv = 1.0 / rmass[i]; + else + minv = 1.0 / mass[type[i]]; + + skip_from_setforce = 0; + for (const auto &fix : fixlist) + if (mask[i] & fix->groupbit) + skip_from_setforce = 1; + + if (skip_from_setforce) + for (int a = 0; a < 3; a++) fp_store[i][a] = f[i][a] * minv; + else + for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } // Forward comm forces From b94c41c05a95459d2a04a83f43768db4b1066052 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 23 Jun 2025 20:35:36 -0600 Subject: [PATCH 3/7] Fixing some static code analysis issues in RHEO --- src/RHEO/bond_rheo_shell.cpp | 16 +++++++--------- src/RHEO/compute_rheo_kernel.cpp | 13 ++++++++----- src/RHEO/compute_rheo_vshift.cpp | 2 ++ src/RHEO/fix_rheo.cpp | 5 ++--- src/RHEO/fix_rheo_thermal.cpp | 2 +- src/RHEO/pair_rheo.cpp | 12 +++++------- 6 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index d5b2dab1b0..4aecb4a1fe 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -217,12 +217,6 @@ void BondRHEOShell::compute(int eflag, int vflag) i2 = itmp; } - delx = x[i1][0] - x[i2][0]; - dely = x[i1][1] - x[i2][1]; - delz = x[i1][2] - x[i2][2]; - rsq = delx * delx + dely * dely + delz * delz; - r = sqrt(rsq); - // If bond hasn't been set - zero data if (t < EPSILON || std::isnan(t)) t = store_bond(n, i1, i2); @@ -472,7 +466,7 @@ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; - for (i = first; i < last; i++) { buf[m++] = dbond[i]; } + for (i = first; i < last; i++) { buf[m++] = ubuf(dbond[i]).d; } return m; } @@ -485,7 +479,7 @@ void BondRHEOShell::unpack_reverse_comm(int n, int *list, double *buf) m = 0; for (i = 0; i < n; i++) { j = list[i]; - dbond[j] += buf[m++]; + dbond[j] += (int) ubuf(buf[m++]).i; } } @@ -495,7 +489,8 @@ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce) { if (type <= 0) return 0.0; - double r0, t; + double r0 = -1; + double t = -1; for (int n = 0; n < atom->num_bond[i]; n++) { if (atom->bond_atom[i][n] == atom->tag[j]) { r0 = fix_bond_history->get_atom_value(i, n, 0); @@ -503,6 +498,9 @@ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce) } } + if (r0 == -1) + error->one(FLERR, "Could not find bond"); + svector[1] = t; if (t < tform) return 0.0; diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index e44e356133..16ecca00e9 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -251,13 +251,14 @@ double ComputeRHEOKernel::calc_w_quintic(double r) double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s; s = r * 3.0 * cutinv; - if (s > 3.0) { w = 0.0; } - if (s <= 3.0) { tmp3 = 3.0 - s; tmp3sq = tmp3 * tmp3; w = tmp3sq * tmp3sq * tmp3; + } else { + w = 0.0; } + if (s <= 2.0) { tmp2 = 2.0 - s; tmp2sq = tmp2 * tmp2; @@ -285,12 +286,14 @@ double ComputeRHEOKernel::calc_dw_scalar_quintic(double r) s = r * 3.0 * cutinv; - if (s > 3.0) { wp = 0.0; } if (s <= 3.0) { tmp3 = 3.0 - s; tmp3sq = tmp3 * tmp3; wp = -5.0 * tmp3sq * tmp3sq; + } else { + wp = 0.0; } + if (s <= 2.0) { tmp2 = 2.0 - s; tmp2sq = tmp2 * tmp2; @@ -855,7 +858,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pb for (int i = 0; i < n; i++) { int j = list[i]; if (comm_stage == 0) { - buf[m++] = coordination[j]; + buf[m++] = ubuf(coordination[j]).d; } else { if (kernel_style == RK0) { buf[m++] = C0[j]; @@ -878,7 +881,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) for (int i = first; i < last; i++) { if (comm_stage == 0) { - coordination[i] = buf[m++]; + coordination[i] = (int) ubuf(buf[m++]).i; } else { if (kernel_style == RK0) { C0[i] = buf[m++]; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 5fbc780e91..03dc001ad2 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -271,6 +271,8 @@ void ComputeRHEOVShift::correct_surfaces() int dim = domain->dimension; double nx, ny, nz, vx, vy, vz, dot; + nz = 0.0; + vz = 0.0; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index c2ff2d8b09..ffaedc2637 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -110,7 +110,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : kernel_style = RK2; } else error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]); - zmin_kernel = utils::numeric(FLERR, arg[5], false, lmp); + zmin_kernel = utils::inumeric(FLERR, arg[5], false, lmp); int iarg = 6; while (iarg < narg) { @@ -342,8 +342,6 @@ void FixRHEO::initial_integrate(int /*vflag*/) double *rmass = atom->rmass; double **gradr = compute_grad->gradr; double **gradv = compute_grad->gradv; - double **vshift; - if (shift_flag) vshift = compute_vshift->vshift; int nlocal = atom->nlocal; int rmass_flag = atom->rmass_flag; @@ -400,6 +398,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Shifting atoms if (shift_flag) { + double vshift = compute_vshift->vshift; for (i = 0; i < nlocal; i++) { if (status[i] & STATUS_NO_SHIFT) continue; if (status[i] & PHASECHECK) continue; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index e9842e19e2..1d3eaf6f19 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -189,7 +189,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg], "react") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal react", error); cut_bond = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - btype = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + btype = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); comm_forward = 4; if (cut_bond <= 0.0) error->all(FLERR, "Illegal max bond length must be greater than zero"); if ((btype < 1) || (btype > atom->nbondtypes)) diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index db83b01539..af9802e07c 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -80,7 +80,7 @@ void PairRHEO::compute(int eflag, int vflag) int i, j, a, b, ii, jj, inum, jnum, itype, jtype; int pair_force_flag, pair_rho_flag, pair_avisc_flag; int fluidi, fluidj; - double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave; + double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, cs_ave; double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, csqi, csqj; double eta_ave, kappa_ave, dT_prefactor; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; @@ -189,11 +189,10 @@ void PairRHEO::compute(int eflag, int vflag) if (!variable_csq) { cs_ave = 0.5 * (cs[itype] + cs[jtype]); } else { - csqi = fix_pressure->calc_csq(rhoi, i); - csqj = fix_pressure->calc_csq(rhoj, j); + csqi = fix_pressure->calc_csq(rho[i], i); + csqj = fix_pressure->calc_csq(rho[j], j); cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj)); } - csq_ave = cs_ave * cs_ave; pair_rho_flag = 0; pair_force_flag = 0; @@ -231,7 +230,7 @@ void PairRHEO::compute(int eflag, int vflag) Pj = fix_pressure->calc_pressure(rhoj, j); if ((chi[j] > 0.9) && (r < (cutk * 0.5))) - fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * csq_ave * cutk * rinv; + fmag = (chi[j] - 0.9) * (cutk * 0.5 - r) * rho0j * cs_ave * cs_ave * cutk * rinv; } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vi, vj, i, j); @@ -239,7 +238,7 @@ void PairRHEO::compute(int eflag, int vflag) Pi = fix_pressure->calc_pressure(rhoi, i); if (chi[i] > 0.9 && r < (cutk * 0.5)) - fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * csq_ave * cutk * rinv; + fmag = (chi[i] - 0.9) * (cutk * 0.5 - r) * rho0i * cs_ave * cs_ave * cutk * rinv; } else if ((!fluidi) && (!fluidj)) { rhoi = rho0i; @@ -251,7 +250,6 @@ void PairRHEO::compute(int eflag, int vflag) csqi = fix_pressure->calc_csq(rhoi, i); csqj = fix_pressure->calc_csq(rhoj, j); cs_ave = 0.5 * (sqrt(csqi) + sqrt(csqj)); - csq_ave = cs_ave * cs_ave; } } From 727c0b251a27a91e3a5730ee76d30b1f64d143c8 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 23 Jun 2025 20:35:53 -0600 Subject: [PATCH 4/7] Clarifying exception in fix deform doc --- doc/src/fix_deform.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_deform.rst b/doc/src/fix_deform.rst index 829e4947eb..93f4136744 100644 --- a/doc/src/fix_deform.rst +++ b/doc/src/fix_deform.rst @@ -298,7 +298,8 @@ Note that it should return the "change" in box length, not the absolute box length. This means it should evaluate to 0.0 when invoked on the initial timestep of the run following the definition of fix deform. It should evaluate to a value > 0.0 to dilate the box at -future times, or a value < 0.0 to compress the box. +future times, or a value < 0.0 to compress the box. The exception +would be if the run command uses the *pre no* option. The variable *name2* must also be an :doc:`equal-style variable ` and should calculate the rate of box length change, in From 7e3d6026d45762501f4ebe27c79b5866dc6f8085 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 23 Jun 2025 20:53:51 -0600 Subject: [PATCH 5/7] Restoring rheo options to set --- src/set.cpp | 48 +++++++++++++++++++++++++++++++++++++++++++++++- src/set.h | 2 ++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/set.cpp b/src/set.cpp index f42dd175ab..f1b1acdd63 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -48,7 +48,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{ANGLE,ANGMOM,BOND,CC,CHARGE,DENSITY,DIAMETER,DIHEDRAL,DIPOLE, DIPOLE_RANDOM,DPD_THETA,EDPD_CV,EDPD_TEMP,EPSILON,IMAGE,IMPROPER,LENGTH, - MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,SHAPE, + MASS,MOLECULE,OMEGA,QUAT,QUAT_RANDOM,RADIUS_ELECTRON,RHEO_STATUS,SHAPE, SMD_CONTACT_RADIUS,SMD_MASS_DENSITY,SPH_CV,SPH_E,SPH_RHO, SPIN_ATOM,SPIN_ATOM_RANDOM,SPIN_ELECTRON,TEMPERATURE,THETA,THETA_RANDOM, TRI,TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,VOLUME,VX,VY,VZ,X,Y,Z, @@ -292,6 +292,14 @@ void Set::process_args(int caller_flag, int narg, char **arg) action->keyword = RADIUS_ELECTRON; process_radius_election(iarg,narg,arg,action); invoke_choice[naction++] = &Set::invoke_radius_election; + } else if (strcmp(arg[iarg],"rheo/rho") == 0) { + action->keyword = SPH_RHO; + process_sph_rho(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_sph_rho; + } else if (strcmp(arg[iarg],"rheo/status") == 0) { + action->keyword = RHEO_STATUS; + process_rheo_status(iarg,narg,arg,action); + invoke_choice[naction++] = &Set::invoke_rheo_status; } else if (strcmp(arg[iarg],"shape") == 0) { action->keyword = SHAPE; process_shape(iarg,narg,arg,action); @@ -1988,6 +1996,44 @@ void Set::invoke_radius_election(Action *action) /* ---------------------------------------------------------------------- */ +void Set::process_rheo_status(int &iarg, int narg, char **arg, Action *action) +{ + if (!atom->rheo_status_flag) + error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style()); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error); + + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1,action); + else { + action->ivalue1 = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (action->ivalue1 != 0 && action->ivalue1 != 1) + error->one(FLERR,"Invalid rheo/status {} in set command", action->ivalue1); + } + + iarg += 2; +} + +void Set::invoke_rheo_status(Action *action) +{ + int nlocal = atom->nlocal; + int *status = atom->rheo_status; + + int varflag = action->varflag; + int rheo_status; + if (!action->varflag1) rheo_status = action->ivalue1; + + for (int i = 0; i < nlocal; i++) { + if (!select[i]) continue; + if (varflag) { + rheo_status = static_cast (vec1[i]); + if (rheo_status != 0 && rheo_status != 1) + error->one(FLERR,"Invalid rheo/status in set command"); + } + status[i] = rheo_status; + } +} + +/* ---------------------------------------------------------------------- */ + void Set::process_shape(int &iarg, int narg, char **arg, Action *action) { if (!atom->ellipsoid_flag) diff --git a/src/set.h b/src/set.h index 65489d9932..3e3e1f3a48 100644 --- a/src/set.h +++ b/src/set.h @@ -110,6 +110,7 @@ class Set : public Command { void process_quat(int &, int, char **, Action *); void process_quat_random(int &, int, char **, Action *); void process_radius_election(int &, int, char **, Action *); + void process_rheo_status(int &, int, char **, Action *); void process_shape(int &, int, char **, Action *); void process_smd_contact_radius(int &, int, char **, Action *); void process_smd_mass_density(int &, int, char **, Action *); @@ -162,6 +163,7 @@ class Set : public Command { void invoke_quat(Action *); void invoke_quat_random(Action *); void invoke_radius_election(Action *); + void invoke_rheo_status(Action *); void invoke_shape(Action *); void invoke_smd_contact_radius(Action *); void invoke_smd_mass_density(Action *); From 08d81a5ac27c91cb8f401af4b26ec0bfc43f8327 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 24 Jun 2025 12:08:54 -0600 Subject: [PATCH 6/7] Correcting simple variable type mistake --- src/RHEO/fix_rheo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index ffaedc2637..423847c715 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -398,7 +398,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Shifting atoms if (shift_flag) { - double vshift = compute_vshift->vshift; + double **vshift = compute_vshift->vshift; for (i = 0; i < nlocal; i++) { if (status[i] & STATUS_NO_SHIFT) continue; if (status[i] & PHASECHECK) continue; From 03611c95aea6f41bfc03c0bd02092d6eeebaa2dc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 24 Jun 2025 20:38:29 -0400 Subject: [PATCH 7/7] improve error messages --- src/RHEO/bond_rheo_shell.cpp | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index 4aecb4a1fe..b74a969cf2 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -147,7 +147,8 @@ void BondRHEOShell::store_data() // map to find index n j = atom->map(atom->bond_atom[i][m]); - if (j == -1) error->one(FLERR, "Atom missing in BPM bond"); + if (j == -1) error->one(FLERR, Error::NOLASTLINE, + "Atom {} missing in BPM bond", atom->bond_atom[i][m]); fix_bond_history->update_atom_value(i, m, 0, 0.0); fix_bond_history->update_atom_value(i, m, 1, 0.0); @@ -352,19 +353,20 @@ void BondRHEOShell::init_style() BondBPM::init_style(); if (comm->ghost_velocity == 0) - error->all(FLERR, "Bond rheo/shell requires ghost atoms store velocity"); + error->all(FLERR, Error::NOLASTLINE, "Bond rheo/shell requires ghost atoms store velocity"); auto fixes = modify->get_fix_by_style("^rheo$"); - if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell"); + if (fixes.size() == 0) + error->all(FLERR, Error::NOLASTLINE, "Need to define fix rheo to use bond rheo/shell"); class FixRHEO *fix_rheo = dynamic_cast(fixes[0]); if (!fix_rheo->surface_flag) - error->all(FLERR, "Bond rheo/shell requires surface calculation in fix rheo"); + error->all(FLERR, Error::NOLASTLINE, "Bond rheo/shell requires surface calculation in fix rheo"); compute_surface = fix_rheo->compute_surface; fixes = modify->get_fix_by_style("^rheo/oxidation$"); if (fixes.size() == 0) - error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell"); + error->all(FLERR, Error::NOLASTLINE, "Need to define fix rheo/oxidation to use bond rheo/shell"); class FixRHEOOxidation *fix_rheo_oxidation = dynamic_cast(fixes[0]); rsurf = fix_rheo_oxidation->rsurf; @@ -383,14 +385,13 @@ void BondRHEOShell::settings(int narg, char **arg) if (strcmp(arg[iarg], "t/form") == 0) { if (iarg + 1 > narg) utils::missing_cmd_args(FLERR, "bond rheo/shell t/form", error); tform = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (tform < 0.0) + error->all(FLERR, iarg + 1, "Illegal bond rheo/shell command, must specify positive formation time"); i += 1; } else { - error->all(FLERR, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]); + error->all(FLERR, iarg, "Illegal bond rheo/shell command, invalid argument {}", arg[iarg]); } } - - if (tform < 0.0) - error->all(FLERR, "Illegal bond rheo/shell command, must specify positive formation time"); } /* ---------------------------------------------------------------------- @@ -491,15 +492,17 @@ double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce) double r0 = -1; double t = -1; + const auto *tag = atom->tag; for (int n = 0; n < atom->num_bond[i]; n++) { - if (atom->bond_atom[i][n] == atom->tag[j]) { + if (atom->bond_atom[i][n] == tag[j]) { r0 = fix_bond_history->get_atom_value(i, n, 0); t = fix_bond_history->get_atom_value(i, n, 1); } } if (r0 == -1) - error->one(FLERR, "Could not find bond"); + error->one(FLERR, Error::NOLASTLINE, + "Could not find bond between atoms {} and {}", tag[i], tag[j]); svector[1] = t; if (t < tform) return 0.0;