From a6cd4a935ebda9328a5c8fa2d134423f5885d707 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Sat, 27 Jun 2020 01:17:46 -0500 Subject: [PATCH 01/31] Replace accumulated displacement by accumulated force for tangential force in styles mindlin and mindlin_rescale. Change documentation accordingly --- doc/src/pair_granular.rst | 65 ++++++++++++++++++++------------ src/GRANULAR/pair_granular.cpp | 69 +++++++++++++++++++++++----------- 2 files changed, 89 insertions(+), 45 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index ee90d1537a..a8509a6019 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -93,7 +93,7 @@ on particle *i* due to contact with particle *j* is given by: .. math:: - \mathbf{F}_{ne, Hooke} = k_N \delta_{ij} \mathbf{n} + \mathbf{F}_{ne, Hooke} = k_n \delta_{ij} \mathbf{n} Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two @@ -106,7 +106,7 @@ For the *hertz* model, the normal component of force is given by: .. math:: - \mathbf{F}_{ne, Hertz} = k_N R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n} + \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n} Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective radius, denoted for simplicity as *R* from here on. For *hertz*\ , the @@ -123,7 +123,7 @@ Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\r modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of types *i* and *j*\ . Note that if the elastic modulus and the shear modulus of the two particles are the same, the *hertz/material* model -is equivalent to the *hertz* model with :math:`k_N = 4/3 E_{eff}` +is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}` The *dmt* model corresponds to the :ref:`(Derjaguin-Muller-Toporov) ` cohesive model, where the force @@ -268,7 +268,7 @@ coefficient, and :math:`k_t` is the tangential stiffness coefficient. For *tangential linear_nohistory*, a simple velocity-dependent Coulomb friction criterion is used, which mimics the behavior of the *pair -gran/hooke* style. The tangential force (\mathbf{F}_t\) is given by: +gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by: .. math:: @@ -294,8 +294,8 @@ keyword also affects the tangential damping. The parameter literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall `, :ref:`Tsuji et al `, :ref:`Silbert et al `). The relative tangential velocity at the point of contact is given by -:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}{n}`, -:math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i`. +:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\mathbf{n}`, +:math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i` . The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` . The normal force value :math:`F_{n0}` used to compute the critical force @@ -319,10 +319,11 @@ form: Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and :math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ . -The remaining tangential options all use accumulated tangential -displacement (i.e. contact history). This is discussed below in the -context of the *linear_history* option, but the same treatment of the -accumulated displacement applies to the other options as well. +The remaining tangential options all use accumulated tangential displacement, +or accumulated tangential force (i.e. contact history). This is discussed +in details below for accumulated tangential displacement in the context +of the *linear_history* option. The same treatment of the accumulated +displacement, or accumulated force, applies to the other options as well. For *tangential linear_history*, the tangential force is given by: @@ -372,7 +373,7 @@ discussion): .. math:: - \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}\right) + \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} - \mathbf{F}_{t,damp}\right) The tangential force is added to the total normal force (elastic plus damping) to produce the total force on the particle. The tangential @@ -387,35 +388,51 @@ overlap region) to induce a torque on each particle according to: \mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t -For *tangential mindlin*\ , the :ref:`Mindlin ` no-slip solution is used, which differs from the *linear_history* -option by an additional factor of *a*\ , the radius of the contact region. The tangential force is given by: +For *tangential mindlin*\ , the :ref:`Mindlin ` no-slip solution +is used which differs from the *linear_history* option by an additional factor +of *a*\ , the radius of the contact region. The tangential stiffness depends +on the radius of the contact region and the force must therefore be computed +incrementally. The accumulated tangential force characterizes the contact history. +The increment of the elastic component of the tangential force :math:`\mathbf{F}_{te}` +is given by: .. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} + \mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau + +The tangential force is given by: + +.. math:: + + \mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}` for all normal contact models, except for *jkr*\ , where it is given implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see -discussion above. To match the Mindlin solution, one should set :math:`k_t = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus -:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This -can also be achieved by specifying *NULL* for :math:`k_t`, in which case a -normal contact model that specifies material parameters :math:`E` and -:math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this -case, mixing of the shear modulus for different particle types *i* and -*j* is done according to: +discussion above. To match the Mindlin solution, one should set +:math:`k_t = 8G_{eff}`, where :math:`G` is the effective shear modulus given by: .. math:: - 1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j + G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1} + +where :math:`G` is the shear modulus, related to Young's modulus :math:`E` +and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be +achieved by specifying *NULL* for :math:`k_t`, in which case a +normal contact model that specifies material parameters :math:`E` and +:math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this +case, mixing of the shear modulus for different particle types *i* and +*j* is done according to the formula above. + + The *mindlin_rescale* option uses the same form as *mindlin*\ , but the -magnitude of the tangential displacement is re-scaled as the contact +magnitude of the tangential force is re-scaled as the contact unloads, i.e. if :math:`a < a_{t_{n-1}}`: .. math:: - \mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}} + \mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}} Here, :math:`t_{n-1}` indicates the value at the previous time step. This rescaling accounts for the fact that a decrease in the diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index fef8ded5f7..ece1e9e211 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -377,6 +377,9 @@ void PairGranular::compute(int eflag, int vflag) // tangential force, including history effects //**************************************** + // For linear forces, history = cumulative tangential displacement + // For nonlinear forces, history = cumulative elastic tangential force + // tangential component vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -424,7 +427,7 @@ void PairGranular::compute(int eflag, int vflag) } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { k_tangential *= a; - // on unloading, rescale the shear displacements + // on unloading, rescale the shear force if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; @@ -432,11 +435,11 @@ void PairGranular::compute(int eflag, int vflag) history[2] *= factor; } } - // rotate and update displacements. + // rotate and update displacements / force // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; + if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON still valid when history is a force??? if (rsht > 0) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -451,18 +454,30 @@ void PairGranular::compute(int eflag, int vflag) history[1] *= scalefac; history[2] *= scalefac; } - // update history - history[0] += vtr1*dt; - history[1] += vtr2*dt; - history[2] += vtr3*dt; - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) - history[3] = a; + // update tangential displacement / force history + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; + } else { // Eq. (18) Thornton et al., 2013 (http://dx.doi.org/10.1016/j.powtec.2012.08.012) + history[0] -= k_tangential*vtr1*dt; + history[1] -= k_tangential*vtr2*dt; + history[2] -= k_tangential*vtr3*dt; + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) + history[3] = a; + } } // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + } else { + fs1 = history[0] - damp_tangential*vtr1; + fs2 = history[1] - damp_tangential*vtr2; + fs3 = history[2] - damp_tangential*vtr3; + } // rescale frictional displacements and forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); @@ -470,12 +485,18 @@ void PairGranular::compute(int eflag, int vflag) shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + - damp_tangential*vtr1); - history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + - damp_tangential*vtr2); - history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + - damp_tangential*vtr3); + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + + damp_tangential*vtr1); + history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + + damp_tangential*vtr2); + history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + + damp_tangential*vtr3); + } else { + history[0] = Fscrit*fs1/fs + damp_tangential*vtr1; + history[1] = Fscrit*fs2/fs + damp_tangential*vtr2; + history[2] = Fscrit*fs3/fs + damp_tangential*vtr3; + } fs1 *= Fscrit/fs; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; @@ -1534,9 +1555,15 @@ double PairGranular::single(int i, int j, int itype, int jtype, history[2]*history[2]); // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + } else { + fs1 = history[0] - damp_tangential*vtr1; + fs2 = history[1] - damp_tangential*vtr2; + fs3 = history[2] - damp_tangential*vtr3; + } // rescale frictional forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); From 57f639c0e535b781738bd394abcae45f2a0d8064 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 12:42:47 -0600 Subject: [PATCH 02/31] bond/react:reset_mol_ids keyword --- src/USER-REACTION/fix_bond_react.cpp | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 2c6c20c844..c908526067 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -144,7 +144,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : int iarg = 3; stabilization_flag = 0; - int num_common_keywords = 1; + reset_mol_ids_flag = 1; + int num_common_keywords = 2; for (int m = 0; m < num_common_keywords; m++) { if (strcmp(arg[iarg],"stabilization") == 0) { if (strcmp(arg[iarg+1],"no") == 0) { @@ -162,6 +163,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : nve_limit_xmax = arg[iarg+3]; iarg += 4; } + } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " + "'reset_mol_ids' keyword has too few arguments"); + if (strcmp(arg[iarg+1],"yes") == 0) iarg += 2; //default + if (strcmp(arg[iarg+1],"no") == 0) { + reset_mol_ids_flag = 0; + iarg += 2; + } } else if (strcmp(arg[iarg],"react") == 0) { break; } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); @@ -2503,11 +2512,14 @@ void FixBondReact::ghost_glovecast() } /* ---------------------------------------------------------------------- -update charges, types, special lists and all topology +update molecule IDs, charges, types, special lists and all topology ------------------------------------------------------------------------- */ void FixBondReact::update_everything() { + if (reset_mol_ids_flag) + input->one("reset_mol_ids " + std::string(group->names[igroup])); + int *type = atom->type; int nlocal = atom->nlocal; From 371a5c5b6156de97ba992b46304b84be2888d47b Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 12:44:34 -0600 Subject: [PATCH 03/31] bond/react: reset_mol_ids docs --- doc/src/fix_bond_react.rst | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index d5f917bfdb..d1671d7c78 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -14,19 +14,22 @@ Syntax react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ... ... -* ID, group-ID are documented in :doc:`fix ` command. Group-ID is ignored. +* ID, group-ID are documented in :doc:`fix ` command. * bond/react = style name of this fix command * the common keyword/values may be appended directly after 'bond/react' * this applies to all reaction specifications (below) -* common_keyword = *stabilization* +* common_keyword = *stabilization* or *reset_mol_ids* .. parsed-literal:: *stabilization* values = *no* or *yes* *group-ID* *xmax* - *no* = no reaction site stabilization + *no* = no reaction site stabilization (default) *yes* = perform reaction site stabilization *group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction *xmax* = xmax value that is used by an internally-created :doc:`nve/limit ` integrator + *reset_mol_ids* values = *yes* or *no* + *yes* = update molecule IDs based on new global topology (default) + *no* = do not update molecule IDs * react = mandatory argument indicating new reaction specification * react-ID = user-assigned name for the reaction @@ -50,9 +53,9 @@ Syntax *stabilize_steps* value = timesteps timesteps = number of timesteps to apply the internally-created :doc:`nve/limit ` fix to reacting atoms *update_edges* value = *none* or *charges* or *custom* - none = do not update topology near the edges of reaction templates - charges = update atomic charges of all atoms in reaction templates - custom = force the update of user-specified atomic charges + *none* = do not update topology near the edges of reaction templates + *charges* = update atomic charges of all atoms in reaction templates + *custom* = force the update of user-specified atomic charges Examples """""""" @@ -154,6 +157,13 @@ due to the internal dynamic grouping performed by fix bond/react. If the group-ID is an existing static group, react-group-IDs should also be specified as this static group, or a subset. +The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids ` +command after a reaction occurs, to ensure that molecule IDs are +consistent with the new bond topology. The group-ID used for +:doc:`reset_mol_ids ` is the group-ID for this fix. +Resetting molecule IDs is necessarily a global operation, and so can +be slow for very large systems. + The following comments pertain to each *react* argument (in other words, can be customized for each reaction, or reaction step): @@ -553,7 +563,7 @@ Default """"""" The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, -update_edges = none +reset_mol_ids = yes, update_edges = none ---------- From da91f81d4062b4acad75c25af8523ea41edb91de Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 13:42:47 -0600 Subject: [PATCH 04/31] bond/react:doc clarification --- doc/src/fix_bond_react.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index d1671d7c78..82ab6beaa3 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -213,9 +213,9 @@ surrounding topology. As described below, the bonding atom pairs of the pre-reacted template are specified by atom ID in the map file. The pre-reacted molecule template should contain as few atoms as possible while still completely describing the topology of all atoms affected -by the reaction. For example, if the force field contains dihedrals, -the pre-reacted template should contain any atom within three bonds of -reacting atoms. +by the reaction (which includes all atoms that change atom type). For +example, if the force field contains dihedrals, the pre-reacted +template should contain any atom within three bonds of reacting atoms. Some atoms in the pre-reacted template that are not reacting may have missing topology with respect to the simulation. For example, the From 6272b7d2bf55e9b8637ae33591d9a0b19d29c5e7 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 13:52:13 -0600 Subject: [PATCH 05/31] add (undocumented) verbosity option to reset_mol_ids --- src/reset_mol_ids.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index e3cc877548..7f71aee0b4 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -56,6 +56,7 @@ void ResetMolIDs::command(int narg, char **arg) int compressflag = 1; int singleflag = 0; tagint offset = -1; + int verbose = 1; int iarg = 1; while (iarg < narg) { @@ -76,6 +77,10 @@ void ResetMolIDs::command(int narg, char **arg) offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command"); iarg += 2; + } else if (strcmp(arg[iarg],"verbose") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); + verbose = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); + iarg += 2; } else error->all(FLERR,"Illegal reset_mol_ids command"); } @@ -203,7 +208,7 @@ void ResetMolIDs::command(int narg, char **arg) MPI_Barrier(world); - if (comm->me == 0) { + if (verbose == 1 && comm->me == 0) { if (nchunk < 0) utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n")); else From a2547701e6bec294d362685442b6252b033ddd5f Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 13:59:30 -0600 Subject: [PATCH 06/31] fix verbose reset_mol_ids --- src/reset_mol_ids.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 7f71aee0b4..f62d1bf402 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -84,7 +84,7 @@ void ResetMolIDs::command(int narg, char **arg) } else error->all(FLERR,"Illegal reset_mol_ids command"); } - if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); + if (verbose == 1 && comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); // record wall time for resetting molecule IDs From e00b0e96f6a3162518ee05f48716dba3a3afe89e Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 14:00:46 -0600 Subject: [PATCH 07/31] bond/react: prevent reset_mol_ids printing --- src/USER-REACTION/fix_bond_react.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index c908526067..ea64fb5447 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -2518,7 +2518,7 @@ update molecule IDs, charges, types, special lists and all topology void FixBondReact::update_everything() { if (reset_mol_ids_flag) - input->one("reset_mol_ids " + std::string(group->names[igroup])); + input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose 0"); int *type = atom->type; From 9ec5708f2fb878a2bfc8cc7b1cd02e2997a96430 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 14:21:10 -0600 Subject: [PATCH 08/31] Update reset_mol_ids.cpp --- src/reset_mol_ids.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index f62d1bf402..ff663ddf24 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -79,7 +79,7 @@ void ResetMolIDs::command(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"verbose") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); - verbose = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); + verbose = utils::inumeric(FLERR,arg[iarg+1],1,lmp); iarg += 2; } else error->all(FLERR,"Illegal reset_mol_ids command"); } From fe6efe8861e287510931835b64dba97bc010946e Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 14:29:39 -0600 Subject: [PATCH 09/31] need header file! --- src/USER-REACTION/fix_bond_react.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index e8a4253ce7..01e8a5ee96 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -61,6 +61,7 @@ class FixBondReact : public Fix { int *max_rxn,*nlocalskips,*nghostlyskips; tagint lastcheck; int stabilization_flag; + int reset_mol_ids_flag; int custom_exclude_flag; int *stabilize_steps_flag; int *update_edges_flag; From edd3fb7108682a853dee3335f57bccf712af4c24 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 20:51:14 -0600 Subject: [PATCH 10/31] reset_mol_ids: documented verbose option --- doc/src/reset_mol_ids.rst | 9 +++++++-- src/USER-REACTION/fix_bond_react.cpp | 2 +- src/reset_mol_ids.cpp | 8 +++++--- 3 files changed, 13 insertions(+), 6 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 0d6063b3ef..9554ac86c2 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -19,6 +19,7 @@ Syntax *compress* value = *yes* or *no* *offset* value = *Noffset* >= -1 *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules + *verbose* value = *yes* or *no* Examples """""""" @@ -95,6 +96,10 @@ is not *all* there may be collisions with the molecule IDs of other atoms. does not perform a full update of the bond topology data structures within LAMMPS. +The *verbose* keyword determines if this command outputs run-time +information to the screen and log file, including the number of new +molecule IDs and CPU time used by the command. + Restrictions """""""""""" none @@ -112,5 +117,5 @@ Related commands Default """"""" -The default keyword settings are compress = yes, single = no, and -offset = -1. +The default keyword settings are compress = yes, single = no, +verbose = yes, and offset = -1. diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index ea64fb5447..54e9720ea7 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -2518,7 +2518,7 @@ update molecule IDs, charges, types, special lists and all topology void FixBondReact::update_everything() { if (reset_mol_ids_flag) - input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose 0"); + input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no"); int *type = atom->type; diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index ff663ddf24..7caa8f1e77 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -79,12 +79,14 @@ void ResetMolIDs::command(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"verbose") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); - verbose = utils::inumeric(FLERR,arg[iarg+1],1,lmp); + if (strcmp(arg[iarg+1],"yes") == 0) verbose = 1; + else if (strcmp(arg[iarg+1],"no") == 0) verbose = 0; + else error->all(FLERR,"Illegal reset_mol_ids command"); iarg += 2; } else error->all(FLERR,"Illegal reset_mol_ids command"); } - if (verbose == 1 && comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); + if (verbose && (comm->me == 0)) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); // record wall time for resetting molecule IDs @@ -208,7 +210,7 @@ void ResetMolIDs::command(int narg, char **arg) MPI_Barrier(world); - if (verbose == 1 && comm->me == 0) { + if (verbose && (comm->me == 0)) { if (nchunk < 0) utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n")); else From 845b9185010b12b8a655c487f9d134d5acc45f52 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 20:59:17 -0600 Subject: [PATCH 11/31] probably better reset_mol_id doc version --- doc/src/reset_mol_ids.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 9554ac86c2..d0e8bc179a 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -97,8 +97,10 @@ is not *all* there may be collisions with the molecule IDs of other atoms. within LAMMPS. The *verbose* keyword determines if this command outputs run-time -information to the screen and log file, including the number of new -molecule IDs and CPU time used by the command. +information to the screen and log file. If the setting is *yes* +(the default), the command prints out information including the number +of new molecule IDs and CPU time used by the command. If the setting +is *no*, no information is printed. Restrictions """""""""""" From 295d75f230d231f5a696fc5d9400a90a4f9f3883 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 18 Jul 2020 21:14:36 -0600 Subject: [PATCH 12/31] guess should reset_mol_IDs *after* updating bonds --- src/USER-REACTION/fix_bond_react.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 54e9720ea7..6ce19bf22a 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -2517,9 +2517,6 @@ update molecule IDs, charges, types, special lists and all topology void FixBondReact::update_everything() { - if (reset_mol_ids_flag) - input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no"); - int *type = atom->type; int nlocal = atom->nlocal; @@ -3054,6 +3051,10 @@ void FixBondReact::update_everything() atom->natoms -= ndel; // done deleting atoms + // reset mol ids + if (reset_mol_ids_flag) + input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no"); + // something to think about: this could done much more concisely if // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG From cf83ce674548f4424166d095816f66f47ed17c30 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sat, 25 Jul 2020 00:52:39 -0600 Subject: [PATCH 13/31] reset_mol_ids->reset() version --- doc/src/reset_mol_ids.rst | 11 +--- src/USER-REACTION/fix_bond_react.cpp | 13 ++++- src/USER-REACTION/fix_bond_react.h | 1 + src/reset_mol_ids.cpp | 84 +++++++++++++++------------- src/reset_mol_ids.h | 7 +++ 5 files changed, 66 insertions(+), 50 deletions(-) diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index d0e8bc179a..0d6063b3ef 100644 --- a/doc/src/reset_mol_ids.rst +++ b/doc/src/reset_mol_ids.rst @@ -19,7 +19,6 @@ Syntax *compress* value = *yes* or *no* *offset* value = *Noffset* >= -1 *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules - *verbose* value = *yes* or *no* Examples """""""" @@ -96,12 +95,6 @@ is not *all* there may be collisions with the molecule IDs of other atoms. does not perform a full update of the bond topology data structures within LAMMPS. -The *verbose* keyword determines if this command outputs run-time -information to the screen and log file. If the setting is *yes* -(the default), the command prints out information including the number -of new molecule IDs and CPU time used by the command. If the setting -is *no*, no information is printed. - Restrictions """""""""""" none @@ -119,5 +112,5 @@ Related commands Default """"""" -The default keyword settings are compress = yes, single = no, -verbose = yes, and offset = -1. +The default keyword settings are compress = yes, single = no, and +offset = -1. diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 6ce19bf22a..16a9ca29f1 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -33,6 +33,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) #include "neigh_list.h" #include "neigh_request.h" #include "random_mars.h" +#include "reset_mol_ids.h" #include "molecule.h" #include "group.h" #include "citeme.h" @@ -92,6 +93,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : fix1 = NULL; fix2 = NULL; fix3 = NULL; + reset_mol_ids = NULL; if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: " "too few arguments"); @@ -166,7 +168,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " "'reset_mol_ids' keyword has too few arguments"); - if (strcmp(arg[iarg+1],"yes") == 0) iarg += 2; //default + if (strcmp(arg[iarg+1],"yes") == 0) { // default + delete reset_mol_ids; + reset_mol_ids = new ResetMolIDs(lmp); + iarg += 2; + } if (strcmp(arg[iarg+1],"no") == 0) { reset_mol_ids_flag = 0; iarg += 2; @@ -508,6 +514,8 @@ FixBondReact::~FixBondReact() } delete [] random; + delete reset_mol_ids; + memory->destroy(partner); memory->destroy(finalpartner); memory->destroy(ncreate); @@ -3052,8 +3060,7 @@ void FixBondReact::update_everything() // done deleting atoms // reset mol ids - if (reset_mol_ids_flag) - input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no"); + if (reset_mol_ids_flag) reset_mol_ids->reset(group->names[igroup]); // something to think about: this could done much more concisely if // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index 01e8a5ee96..1b98614064 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -94,6 +94,7 @@ class FixBondReact : public Fix { class RanMars **random; // random number for 'prob' keyword class RanMars **rrhandom; // random number for Arrhenius constraint class NeighList *list; + class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs int *reacted_mol,*unreacted_mol; int *limit_duration; // indicates how long to relax diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 7caa8f1e77..e6add7df9b 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -33,9 +33,15 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {} +ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) { + compressflag = 1; + singleflag = 0; + offset = -1; +} -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + called as reset_mol_ids command in input script +------------------------------------------------------------------------- */ void ResetMolIDs::command(int narg, char **arg) { @@ -49,14 +55,7 @@ void ResetMolIDs::command(int narg, char **arg) // process args if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command"); - int igroup = group->find(arg[0]); - if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); - int groupbit = group->bitmask[igroup]; - - int compressflag = 1; - int singleflag = 0; - tagint offset = -1; - int verbose = 1; + char *groupid = arg[0]; int iarg = 1; while (iarg < narg) { @@ -77,39 +76,61 @@ void ResetMolIDs::command(int narg, char **arg) offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp); if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command"); iarg += 2; - } else if (strcmp(arg[iarg],"verbose") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command"); - if (strcmp(arg[iarg+1],"yes") == 0) verbose = 1; - else if (strcmp(arg[iarg+1],"no") == 0) verbose = 0; - else error->all(FLERR,"Illegal reset_mol_ids command"); - iarg += 2; } else error->all(FLERR,"Illegal reset_mol_ids command"); } - if (verbose && (comm->me == 0)) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); + if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n"); // record wall time for resetting molecule IDs MPI_Barrier(world); double time1 = MPI_Wtime(); + // initialize system since comm->borders() will be invoked + + lmp->init(); + + // reset molecule IDs + + reset(groupid); + + // total time + + MPI_Barrier(world); + + if (comm->me == 0) { + if (nchunk < 0) + utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n")); + else + utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); + utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", + MPI_Wtime()-time1)); + } +} + +/* ---------------------------------------------------------------------- + called from command() and directly from fixes that update molecules +------------------------------------------------------------------------- */ + +void ResetMolIDs::reset(char *groupid) +{ + int igroup = group->find(groupid); + if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); + int groupbit = group->bitmask[igroup]; + // create instances of compute fragment/atom, compute reduce (if needed), // and compute chunk/atom. all use the group-ID for this command const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; if (singleflag) - modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0])); + modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid)); else - modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0])); + modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid)); const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; if (compressflag) modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", - idchunk,arg[0])); - - // initialize system since comm->borders() will be invoked - - lmp->init(); + idchunk,groupid)); // setup domain, communication // exchange will clear map, borders will reset @@ -145,7 +166,7 @@ void ResetMolIDs::command(int narg, char **arg) // if compressflag = 0, done // set nchunk = -1 since cannot easily determine # of new molecule IDs - int nchunk = -1; + nchunk = -1; // if compressflag = 1, invoke peratom method of compute chunk/atom // will compress new molecule IDs to be contiguous 1 to Nmol @@ -205,17 +226,4 @@ void ResetMolIDs::command(int narg, char **arg) modify->delete_compute(idfrag); if (compressflag) modify->delete_compute(idchunk); - - // total time - - MPI_Barrier(world); - - if (verbose && (comm->me == 0)) { - if (nchunk < 0) - utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n")); - else - utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk)); - utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n", - MPI_Wtime()-time1)); - } } diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h index 83c0b5d80f..cb7c472724 100644 --- a/src/reset_mol_ids.h +++ b/src/reset_mol_ids.h @@ -28,6 +28,13 @@ class ResetMolIDs : protected Pointers { public: ResetMolIDs(class LAMMPS *); void command(int, char **); + void reset(char *); + +private: + int nchunk; + int compressflag; // 1 = contiguous values for new IDs + int singleflag; // 0 = mol IDs of single atoms set to 0 + tagint offset; // offset for contiguous mol ID values }; } From bf724332d480d80578b1fd76e614c077a86fb831 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 10:53:30 -0500 Subject: [PATCH 14/31] implement tangential force history in mindlin/force and mindlin_rescale/force --- doc/src/pair_granular.rst | 134 ++++++++++++++++++++++++--------- src/GRANULAR/pair_granular.cpp | 83 ++++++++++++++------ 2 files changed, 158 insertions(+), 59 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index a8509a6019..ab735bdfde 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -140,7 +140,7 @@ where the force is computed as: \mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n} -Here, *a* is the radius of the contact zone, related to the overlap +Here, :math:`a` is the radius of the contact zone, related to the overlap :math:`\delta` according to: .. math:: @@ -167,7 +167,7 @@ following general form: \mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel} -Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n} \mathbf{n}` is the component of relative velocity along +Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n}\ \mathbf{n}` is the component of relative velocity along :math:`\mathbf{n}`. The optional *damping* keyword to the *pair_coeff* command followed by @@ -259,7 +259,9 @@ tangential model choices and their expected parameters are as follows: 1. *linear_nohistory* : :math:`x_{\gamma,t}`, :math:`\mu_s` 2. *linear_history* : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s` 3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` -4. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` +4. *mindlin/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` +5. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` +5. *mindlin_rescale/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal damping :math:`\eta_n` that determines the magnitude of the tangential @@ -272,7 +274,7 @@ gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by: .. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} + \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by: @@ -294,7 +296,7 @@ keyword also affects the tangential damping. The parameter literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall `, :ref:`Tsuji et al `, :ref:`Silbert et al `). The relative tangential velocity at the point of contact is given by -:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\mathbf{n}`, +:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\ \mathbf{n}`, :math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i` . The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` . @@ -314,22 +316,24 @@ form: .. math:: - F_{n0} = \|\mathbf{F}_ne + 2 F_{pulloff}\| + F_{n0} = \|\mathbf{F}_{ne} + 2 F_{pulloff}\| Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and :math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ . -The remaining tangential options all use accumulated tangential displacement, -or accumulated tangential force (i.e. contact history). This is discussed -in details below for accumulated tangential displacement in the context -of the *linear_history* option. The same treatment of the accumulated -displacement, or accumulated force, applies to the other options as well. +The remaining tangential options all use accumulated tangential +displacement (i.e. contact history), except for the options +*mindlin/force* and *mindlin_rescale/force*, that use accumulated +tangential force instead, and are discussed further below. +The accumulated tangential displacement is discussed in details below +in the context of the *linear_history* option. The same treatment of +the accumulated displacement applies to the other options as well. For *tangential linear_history*, the tangential force is given by: .. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} + \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated during the entire duration of the contact: @@ -390,27 +394,18 @@ overlap region) to induce a torque on each particle according to: For *tangential mindlin*\ , the :ref:`Mindlin ` no-slip solution is used which differs from the *linear_history* option by an additional factor -of *a*\ , the radius of the contact region. The tangential stiffness depends -on the radius of the contact region and the force must therefore be computed -incrementally. The accumulated tangential force characterizes the contact history. -The increment of the elastic component of the tangential force :math:`\mathbf{F}_{te}` -is given by: +of :math:`a`, the radius of the contact region. The tangential force is given by: .. math:: - \mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau - -The tangential force is given by: + \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} -.. math:: - \mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} - -Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}` +Here, :math:`a` is the radius of the contact region, given by :math:`a =\sqrt{R\delta}` for all normal contact models, except for *jkr*\ , where it is given implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see discussion above. To match the Mindlin solution, one should set -:math:`k_t = 8G_{eff}`, where :math:`G` is the effective shear modulus given by: +:math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by: .. math:: @@ -424,23 +419,80 @@ normal contact model that specifies material parameters :math:`E` and case, mixing of the shear modulus for different particle types *i* and *j* is done according to the formula above. +.. note:: + The radius of the contact region :math:`a` depends on the normal overlap. + As a result, the tangential force for *mindlin* can change due to + a variation in normal overlap, even with no change in tangential displacement. + +For *tangential mindlin/force*, the accumulated elastic tangential force +characterizes the contact history, instead of the accumulated tangential +displacement. This prevents the dependence of the tangential force on the +normal overlap as noted above. The tangential force is given by: + +.. math:: + + \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} + +The increment of the elastic component of the tangential force +:math:`\mathbf{F}_{te}` is given by: + +.. math:: + + \mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau + +The changes in frame of reference of the contacting pair of particles during +contact are accounted for by the same formula as above, replacing the +accumulated tangential displacement :math:`\xi`, by the accumulated tangential +elastic force :math:`F_{te}`. When the tangential force exceeds the critical +force, the tangential force is directly re-scaled to match the value for +the critical force: + +.. math:: + + \mathbf{F}_{te} = - \mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp} + +The same rules as those described for *mindlin* apply regarding the tangential +stiffness and mixing of the shear modulus for different particle types. The *mindlin_rescale* option uses the same form as *mindlin*\ , but the -magnitude of the tangential force is re-scaled as the contact +magnitude of the tangential displacement is re-scaled as the contact +unloads, i.e. if :math:`a < a_{t_{n-1}}`: + +.. math:: + + \mathbf{\xi} = \mathbf{\xi_{t_{n-1}}} \frac{a}{a_{t_{n-1}}} + +Here, :math:`t_{n-1}` indicates the value at the previous time +step. This rescaling accounts for the fact that a decrease in the +contact area upon unloading leads to the contact being unable to +support the previous tangential loading, and spurious energy is +created without the rescaling above (:ref:`Walton ` ). + +.. note:: + + For *mindlin*, a decrease in the tangential force already occurs as the + contact unloads, due to the dependence of the tangential force on the normal + force described above. By re-scaling :math:`\xi`, *mindlin_rescale* + effectively re-scales the tangential force twice, i.e., proportionally to + :math:`a^2`. This peculiar behavior results from use of the accumulated + tangential displacement to characterize the contact history. Although + *mindlin_rescale* remains available for historic reasons and backward + compatibility purposes, it should be avoided in favor of *mindlin_rescale/force*. + +The *mindlin_rescale/force* option uses the same form as *mindlin/force*, +but the magnitude of the tangential elastic force is re-scaled as the contact unloads, i.e. if :math:`a < a_{t_{n-1}}`: .. math:: \mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}} -Here, :math:`t_{n-1}` indicates the value at the previous time -step. This rescaling accounts for the fact that a decrease in the -contact area upon unloading leads to the contact being unable to -support the previous tangential loading, and spurious energy is -created without the rescaling above (:ref:`Walton ` ). See also -discussion in :ref:`Thornton et al, 2013 ` , particularly -equation 18(b) of that work and associated discussion. +This approach provides a better approximation of the :ref:`Mindlin-Deresiewicz ` +laws and is more consistent than *mindlin_rescale*. See discussions in +:ref:`Thornton et al, 2013 `, particularly equation 18(b) of that +work and associated discussion, and :ref:`Agnolin and Roux, 2007 `, +particularly Appendix A. ---------- @@ -477,7 +529,7 @@ exceeds a critical value: .. math:: - \mathbf{F}_{roll} = min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k} + \mathbf{F}_{roll} = \min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k} Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\|\mathbf{v}_{roll}\|` is the direction of the pseudo-force. As with tangential displacement, the rolling @@ -529,7 +581,7 @@ is then truncated according to: .. math:: - \tau_{twist} = min(\mu_{twist} F_{n,0}, \tau_{twist,0}) + \tau_{twist} = \min(\mu_{twist} F_{n,0}, \tau_{twist,0}) Similar to the sliding and rolling displacement, the angular displacement is rescaled so that it corresponds to the critical value @@ -794,3 +846,15 @@ Technology, 233, 30-46. .. _WaltonPC: **(Otis R. Walton)** Walton, O.R., Personal Communication + +** _Mindlin1953: + +**(Mindlin and Deresiewicz, 1953)** Mindlin, R.D., & Deresiewicz, H (1953). +Elastic Spheres in Contact under Varying Oblique Force. +J. Appl. Mech., ASME 20, 327-344. + +** _AgnolinRoux2007: + +**(Agnolin and Roux 2007)** Agnolin, I. & Roux, J-N. (2007). +Internal states of model isotropic granular packings. +I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302. \ No newline at end of file diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index ece1e9e211..e732453757 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -54,7 +54,8 @@ using namespace MathSpecial; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_SDS}; @@ -377,8 +378,11 @@ void PairGranular::compute(int eflag, int vflag) // tangential force, including history effects //**************************************** - // For linear forces, history = cumulative tangential displacement - // For nonlinear forces, history = cumulative elastic tangential force + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force // tangential component vt1 = vr1 - vn1; @@ -422,12 +426,15 @@ void PairGranular::compute(int eflag, int vflag) damp_normal_prefactor; if (tangential_history) { - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) { k_tangential *= a; } else if (tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE) { + TANGENTIAL_MINDLIN_RESCALE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) { k_tangential *= a; - // on unloading, rescale the shear force + // on unloading, rescale the shear displacements/force if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; @@ -435,11 +442,11 @@ void PairGranular::compute(int eflag, int vflag) history[2] *= factor; } } - // rotate and update displacements / force + // rotate and update displacements / force. // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON still valid when history is a force??? + if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON = 1e-10 can fail for small granular media: if (rsht > 0) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -454,22 +461,31 @@ void PairGranular::compute(int eflag, int vflag) history[1] *= scalefac; history[2] *= scalefac; } - // update tangential displacement / force history - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + // update history + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + // tangential displacement history[0] += vtr1*dt; history[1] += vtr2*dt; history[2] += vtr3*dt; - } else { // Eq. (18) Thornton et al., 2013 (http://dx.doi.org/10.1016/j.powtec.2012.08.012) + } else { + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 history[0] -= k_tangential*vtr1*dt; history[1] -= k_tangential*vtr2*dt; history[2] -= k_tangential*vtr3*dt; - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) - history[3] = a; } + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) + history[3] = a; } // tangential forces = history + tangential velocity damping - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { fs1 = -k_tangential*history[0] - damp_tangential*vtr1; fs2 = -k_tangential*history[1] - damp_tangential*vtr2; fs3 = -k_tangential*history[2] - damp_tangential*vtr3; @@ -485,7 +501,10 @@ void PairGranular::compute(int eflag, int vflag) shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); if (shrmag != 0.0) { - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE) { history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + @@ -760,7 +779,8 @@ void PairGranular::coeff(int narg, char **arg) //Defaults normal_model_one = tangential_model_one = -1; - roll_model_one = twist_model_one = 0; + roll_model_one = ROLL_NONE; + twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; int iarg = 2; @@ -846,7 +866,9 @@ void PairGranular::coeff(int narg, char **arg) iarg += 4; } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) || - (strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { + (strcmp(arg[iarg+1], "mindlin_rescale") == 0) || + (strcmp(arg[iarg+1], "mindlin/force") == 0) || + (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) { if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, " "not enough parameters provided for tangential model"); @@ -856,9 +878,15 @@ void PairGranular::coeff(int narg, char **arg) tangential_model_one = TANGENTIAL_MINDLIN; else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; + else if (strcmp(arg[iarg+1], "mindlin/force") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_FORCE; + else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE; tangential_history = 1; if ((tangential_model_one == TANGENTIAL_MINDLIN || - tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) && + tangential_model_one == TANGENTIAL_MINDLIN_RESCALE || + tangential_model_one == TANGENTIAL_MINDLIN_FORCE || + tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) && (strcmp(arg[iarg+2], "NULL") == 0)) { if (normal_model_one == HERTZ || normal_model_one == HOOKE) { error->all(FLERR, "NULL setting for Mindlin tangential " @@ -1040,7 +1068,8 @@ void PairGranular::init_style() } for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE) { + if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE || + tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) { size_history += 1; roll_history_index += 1; twist_history_index += 1; @@ -1510,6 +1539,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, // tangential force, including history effects //**************************************** + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + // tangential component vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -1545,9 +1580,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; if (tangential_history) { - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { - k_tangential *= a; - } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) { k_tangential *= a; } @@ -1555,7 +1588,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, history[2]*history[2]); // tangential forces = history + tangential velocity damping - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY) { + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { fs1 = -k_tangential*history[0] - damp_tangential*vtr1; fs2 = -k_tangential*history[1] - damp_tangential*vtr2; fs3 = -k_tangential*history[2] - damp_tangential*vtr3; @@ -1572,7 +1607,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs1 *= Fscrit/fs; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; - fs *= Fscrit/fs; + fs *= Fscrit/fs; } else fs1 = fs2 = fs3 = fs = 0.0; } From 2de98999c139c0694a263fb90556875add79e7e1 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 14:51:00 -0500 Subject: [PATCH 15/31] bug fix formula for frame of reference rotation for granular tangential history --- doc/src/pair_granular.rst | 2 +- src/GRANULAR/pair_granular.cpp | 57 ++++++++++++++++++---------------- 2 files changed, 31 insertions(+), 28 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index ab735bdfde..0b6e12ade7 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -361,7 +361,7 @@ work: .. math:: - \mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'}\| - \mathbf{n}\cdot\mathbf{\xi'}} + \mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'} - (\mathbf{n}\cdot\mathbf{\xi'})\mathbf{n}\|} Here, :math:`\mathbf{\xi'}` is the accumulated displacement prior to the current time step and :math:`\mathbf{\xi}` is the corrected diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index e732453757..9bde81b981 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -446,20 +446,20 @@ void PairGranular::compute(int eflag, int vflag) // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; // EPSILON = 1e-10 can fail for small granular media: - if (rsht > 0) { - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - // if rsht == shrmag, contacting pair has rotated 90 deg - // in one step, in which case you deserve a crash! - scalefac = shrmag/(shrmag - rsht); - history[0] -= rsht*nx; - history[1] -= rsht*ny; - history[2] -= rsht*nz; - // also rescale to preserve magnitude - history[0] *= scalefac; - history[1] *= scalefac; - history[2] *= scalefac; + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + // projection + history[0] -= rsht*nx; + history[1] -= rsht*ny; + history[2] -= rsht*nz; + + // also rescale to preserve magnitude + prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + scalefac = shrmag/prjmag; + history[0] *= scalefac; + history[1] *= scalefac; + history[2] *= scalefac; } // update history if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || @@ -559,19 +559,22 @@ void PairGranular::compute(int eflag, int vflag) rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; if (historyupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { // rotate into tangential plane - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - scalefac = rollmag/(rollmag - rolldotn); - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - // also rescale to preserve magnitude - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; + // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // projection + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + scalefac = rollmag/prjmag; + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; } history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; From 5ebac27fd5d45492b02ae1d46d1fe8f345593a9e Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Mon, 10 Aug 2020 15:15:47 -0500 Subject: [PATCH 16/31] safety for division by zero in scaling of the projection --- src/GRANULAR/pair_granular.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 9bde81b981..f59c00e90a 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -181,7 +181,7 @@ void PairGranular::compute(int eflag, int vflag) double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht,prjmag; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; @@ -456,11 +456,12 @@ void PairGranular::compute(int eflag, int vflag) // also rescale to preserve magnitude prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); - scalefac = shrmag/prjmag; + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; history[0] *= scalefac; history[1] *= scalefac; history[2] *= scalefac; - } + // update history if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || @@ -571,11 +572,12 @@ void PairGranular::compute(int eflag, int vflag) prjmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + history[rhist2]*history[rhist2]); - scalefac = rollmag/prjmag; + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; - } + history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; history[rhist2] += vrl3*dt; From f6d91b3b2c92d4f8c820265a3bcea5a6649c9daa Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 11 Aug 2020 15:02:37 -0400 Subject: [PATCH 17/31] move domain/comm commands --- src/reset_mol_ids.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index e6add7df9b..aa0fe37ce7 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -90,6 +90,18 @@ void ResetMolIDs::command(int narg, char **arg) lmp->init(); + // setup domain, communication + // exchange will clear map, borders will reset + // this is the map needed to lookup current global IDs for bond topology + + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + comm->setup(); + comm->exchange(); + comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + // reset molecule IDs reset(groupid); @@ -132,18 +144,6 @@ void ResetMolIDs::reset(char *groupid) modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", idchunk,groupid)); - // setup domain, communication - // exchange will clear map, borders will reset - // this is the map needed to lookup current global IDs for bond topology - - if (domain->triclinic) domain->x2lamda(atom->nlocal); - domain->pbc(); - domain->reset_box(); - comm->setup(); - comm->exchange(); - comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); - // invoke peratom method of compute fragment/atom // walks bond connectivity and assigns each atom a fragment ID // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0 From 3c69ebc6696f83f81f5b6f28bce35193d5bf93c6 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 11 Aug 2020 17:12:36 -0400 Subject: [PATCH 18/31] reset_mold_ids: add create_computes --- src/reset_mol_ids.cpp | 55 ++++++++++++++++++++++++++++++------------- src/reset_mol_ids.h | 10 +++++++- 2 files changed, 48 insertions(+), 17 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index aa0fe37ce7..efeb6a672a 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -17,7 +17,6 @@ #include "reset_mol_ids.h" #include -#include #include "atom.h" #include "domain.h" #include "comm.h" @@ -34,11 +33,24 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) { + cfa = NULL; + cca = NULL; + + // default settings + compressflag = 1; singleflag = 0; offset = -1; } +/* ---------------------------------------------------------------------- */ + +ResetMolIDs::~ResetMolIDs() +{ + modify->delete_compute(idfrag); + if (compressflag) modify->delete_compute(idchunk); +} + /* ---------------------------------------------------------------------- called as reset_mol_ids command in input script ------------------------------------------------------------------------- */ @@ -102,9 +114,13 @@ void ResetMolIDs::command(int narg, char **arg) comm->borders(); if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + // create computes + + create_computes(groupid); + // reset molecule IDs - reset(groupid); + reset(); // total time @@ -121,35 +137,49 @@ void ResetMolIDs::command(int narg, char **arg) } /* ---------------------------------------------------------------------- - called from command() and directly from fixes that update molecules + create computes used by reset_mol_ids ------------------------------------------------------------------------- */ -void ResetMolIDs::reset(char *groupid) +void ResetMolIDs::create_computes(char *groupid) { int igroup = group->find(groupid); if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); - int groupbit = group->bitmask[igroup]; + groupbit = group->bitmask[igroup]; // create instances of compute fragment/atom, compute reduce (if needed), // and compute chunk/atom. all use the group-ID for this command - const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM"; + idfrag = "reset_mol_ids_FRAGMENT_ATOM"; if (singleflag) modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid)); else modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid)); - const std::string idchunk = "reset_mol_ids_CHUNK_ATOM"; + idchunk = "reset_mol_ids_CHUNK_ATOM"; if (compressflag) modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", idchunk,groupid)); + int icompute = modify->find_compute(idfrag); + cfa = (ComputeFragmentAtom *) modify->compute[icompute]; + + + if (compressflag) { + icompute = modify->find_compute(idchunk); + cca = (ComputeChunkAtom *) modify->compute[icompute]; + } +} + +/* ---------------------------------------------------------------------- + called from command() and directly from fixes that update molecules +------------------------------------------------------------------------- */ + +void ResetMolIDs::reset() +{ // invoke peratom method of compute fragment/atom // walks bond connectivity and assigns each atom a fragment ID // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0 - int icompute = modify->find_compute(idfrag); - ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute]; cfa->compute_peratom(); double *fragIDs = cfa->vector_atom; @@ -173,8 +203,6 @@ void ResetMolIDs::reset(char *groupid) // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int if (compressflag) { - icompute = modify->find_compute(idchunk); - ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute]; cca->compute_peratom(); double *chunkIDs = cca->vector_atom; nchunk = cca->nchunk; @@ -221,9 +249,4 @@ void ResetMolIDs::reset(char *groupid) } } } - - // clean up - - modify->delete_compute(idfrag); - if (compressflag) modify->delete_compute(idchunk); } diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h index cb7c472724..62543f6eda 100644 --- a/src/reset_mol_ids.h +++ b/src/reset_mol_ids.h @@ -21,20 +21,28 @@ CommandStyle(reset_mol_ids,ResetMolIDs) #define LMP_RESET_MOL_IDS_H #include "pointers.h" +#include namespace LAMMPS_NS { class ResetMolIDs : protected Pointers { public: ResetMolIDs(class LAMMPS *); + ~ResetMolIDs(); void command(int, char **); - void reset(char *); + void create_computes(char *); + void reset(); private: + std::string idfrag, idchunk; int nchunk; + int groupbit; int compressflag; // 1 = contiguous values for new IDs int singleflag; // 0 = mol IDs of single atoms set to 0 tagint offset; // offset for contiguous mol ID values + + class ComputeFragmentAtom *cfa; + class ComputeChunkAtom *cca; }; } From 9d486d734b3365e835f3b28b17eaa115c2f741b1 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 11 Aug 2020 17:29:27 -0400 Subject: [PATCH 19/31] update bond/react for reset_mol_ids->create_computes --- src/USER-REACTION/fix_bond_react.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 16a9ca29f1..9b33630b59 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -171,6 +171,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg+1],"yes") == 0) { // default delete reset_mol_ids; reset_mol_ids = new ResetMolIDs(lmp); + reset_mol_ids->create_computes(group->names[igroup]); iarg += 2; } if (strcmp(arg[iarg+1],"no") == 0) { @@ -246,9 +247,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)"); strncpy(rxn_name[rxn],arg[iarg++],n); - int igroup = group->find(arg[iarg++]); - if (igroup == -1) error->all(FLERR,"Could not find fix group ID"); - groupbits[rxn] = group->bitmask[igroup]; + int groupid = group->find(arg[iarg++]); + if (groupid == -1) error->all(FLERR,"Could not find fix group ID"); + groupbits[rxn] = group->bitmask[groupid]; if (strncmp(arg[iarg],"v_",2) == 0) { n = strlen(&arg[iarg][2]) + 1; @@ -640,9 +641,9 @@ void FixBondReact::post_constructor() group->assign(cmd); if (stabilization_flag == 1) { - int igroup = group->find(exclude_group); + int groupid = group->find(exclude_group); // create exclude_group if not already existing, or use as parent group if static - if (igroup == -1 || group->dynamic[igroup] == 0) { + if (groupid == -1 || group->dynamic[groupid] == 0) { // create stabilization per-atom property cmd = std::string("bond_react_stabilization_internal"); id_fix3 = new char[cmd.size()+1]; @@ -672,7 +673,7 @@ void FixBondReact::post_constructor() strcat(exclude_group,"_REACT"); group->find_or_create(exclude_group); - if (igroup == -1) + if (groupid == -1) cmd = fmt::format("{} dynamic all property statted_tags",exclude_group); else cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group); @@ -3060,7 +3061,7 @@ void FixBondReact::update_everything() // done deleting atoms // reset mol ids - if (reset_mol_ids_flag) reset_mol_ids->reset(group->names[igroup]); + if (reset_mol_ids_flag) reset_mol_ids->reset(); // something to think about: this could done much more concisely if // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG From df497e4853ac26faf7d13d43efbcb7bfdc254aec Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 18 Aug 2020 17:53:07 -0400 Subject: [PATCH 20/31] bond/react: clarify bond-type-checking error --- src/USER-REACTION/fix_bond_react.cpp | 6 +++--- src/USER-REACTION/fix_bond_react.h | 4 +++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 9b33630b59..57bc42cffd 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -2079,7 +2079,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) for (int i = 0; i < twomol->natoms; i++) { if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) { char str[128]; - snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); + snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]); error->all(FLERR,str); } } @@ -2098,7 +2098,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) { if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) { char str[128]; - snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); + snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]); error->all(FLERR,str); } } @@ -2110,7 +2110,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn) if (onemol_batom == equivalences[i][1][myrxn]) { if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) { char str[128]; - snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); + snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]); error->all(FLERR,str); } } diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h index 1b98614064..61a1bd4213 100644 --- a/src/USER-REACTION/fix_bond_react.h +++ b/src/USER-REACTION/fix_bond_react.h @@ -242,8 +242,10 @@ E: Bond/react: Invalid template atom ID in map file Atom IDs in molecule templates range from 1 to the number of atoms in the template. E or W: Bond/react: Atom affected by reaction %s too close to template edge + Bond/react: Atom type affected by reaction %s too close to template edge + Bond/react: Bond type affected by reaction %s too close to template edge -This means an atom which changes type or connectivity during the +This means an atom (or bond) that changes type or connectivity during the reaction is too close to an 'edge' atom defined in the map file. This could cause incorrect assignment of bonds, angle, etc. Generally, this means you must include more atoms in your templates, such that there From e4ab49c2e5d0856b78047088f0442be2d801bb8f Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Tue, 18 Aug 2020 18:12:01 -0400 Subject: [PATCH 21/31] bond/react: bond-type-checking docs --- doc/src/Errors_messages.rst | 3 ++- doc/src/fix_bond_react.rst | 7 ++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index f3be94a239..96e21da681 100644 --- a/doc/src/Errors_messages.rst +++ b/doc/src/Errors_messages.rst @@ -502,7 +502,8 @@ Doc page with :doc:`WARNING messages ` *Bond/react: Unknown section in map file* Please ensure reaction map files are properly formatted. -*Bond/react: Atom affected by reaction too close to template edge* +*Bond/react: Atom type affected by reaction too close to template edge* +*Bond/react: Bond type affected by reaction too close to template edge* This means an atom which changes type or connectivity during the reaction is too close to an 'edge' atom defined in the map file. This could cause incorrect assignment of bonds, angle, etc. diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst index 82ab6beaa3..bbc6fbf864 100644 --- a/doc/src/fix_bond_react.rst +++ b/doc/src/fix_bond_react.rst @@ -213,9 +213,10 @@ surrounding topology. As described below, the bonding atom pairs of the pre-reacted template are specified by atom ID in the map file. The pre-reacted molecule template should contain as few atoms as possible while still completely describing the topology of all atoms affected -by the reaction (which includes all atoms that change atom type). For -example, if the force field contains dihedrals, the pre-reacted -template should contain any atom within three bonds of reacting atoms. +by the reaction (which includes all atoms that change atom type or +connectivity, and all bonds that change bond type). For example, if +the force field contains dihedrals, the pre-reacted template should +contain any atom within three bonds of reacting atoms. Some atoms in the pre-reacted template that are not reacting may have missing topology with respect to the simulation. For example, the From e742ae747532dd7e12558974216dab5fd45cc4d2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Aug 2020 00:31:55 -0400 Subject: [PATCH 22/31] fix RST syntax and spelling issues in granular pair style docs --- doc/src/pair_granular.rst | 8 ++++---- doc/utils/sphinx-config/false_positives.txt | 3 +++ 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 40f3888e8f..d846e79e35 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -261,7 +261,7 @@ tangential model choices and their expected parameters are as follows: 3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 4. *mindlin/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 5. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` -5. *mindlin_rescale/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` +6. *mindlin_rescale/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal damping :math:`\eta_n` that determines the magnitude of the tangential @@ -831,14 +831,14 @@ Technology, 233, 30-46. **(Otis R. Walton)** Walton, O.R., Personal Communication -** _Mindlin1953: +.. _Mindlin1953: **(Mindlin and Deresiewicz, 1953)** Mindlin, R.D., & Deresiewicz, H (1953). Elastic Spheres in Contact under Varying Oblique Force. J. Appl. Mech., ASME 20, 327-344. -** _AgnolinRoux2007: +.. _AgnolinRoux2007: **(Agnolin and Roux 2007)** Agnolin, I. & Roux, J-N. (2007). Internal states of model isotropic granular packings. -I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302. \ No newline at end of file +I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 2bf3e917cf..443478a9fc 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -43,6 +43,7 @@ Afshar agilio Agilio agni +Agnolin Ai Aidan aij @@ -599,6 +600,7 @@ Dequidt der dereference derekt +Deresiewicz Derjagin Derjaguin Derlet @@ -2219,6 +2221,7 @@ oxdna oxrna oxDNA oxRNA +packings padua Padua palegoldenrod From ee6ef98b9be42a4d20cce584b5acaa175666f98d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 21 Aug 2020 00:42:57 -0400 Subject: [PATCH 23/31] remove trailing whitespace --- doc/src/pair_granular.rst | 6 +++--- src/GRANULAR/pair_granular.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index d846e79e35..2fce2d24f1 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -410,7 +410,7 @@ discussion above. To match the Mindlin solution, one should set .. math:: G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1} - + where :math:`G` is the shear modulus, related to Young's modulus :math:`E` and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be achieved by specifying *NULL* for :math:`k_t`, in which case a @@ -476,8 +476,8 @@ created without the rescaling above (:ref:`Walton ` ). force described above. By re-scaling :math:`\xi`, *mindlin_rescale* effectively re-scales the tangential force twice, i.e., proportionally to :math:`a^2`. This peculiar behavior results from use of the accumulated - tangential displacement to characterize the contact history. Although - *mindlin_rescale* remains available for historic reasons and backward + tangential displacement to characterize the contact history. Although + *mindlin_rescale* remains available for historic reasons and backward compatibility purposes, it should be avoided in favor of *mindlin_rescale/force*. The *mindlin_rescale/force* option uses the same form as *mindlin/force*, diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index b16fd40236..1b92c2686d 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -447,7 +447,7 @@ void PairGranular::compute(int eflag, int vflag) history[0] -= rsht*nx; history[1] -= rsht*ny; history[2] -= rsht*nz; - + // also rescale to preserve magnitude prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -572,7 +572,7 @@ void PairGranular::compute(int eflag, int vflag) history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; - + history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; history[rhist2] += vrl3*dt; From 921b6d8135d4aa92c9050e7a33d6e05b77a94d6e Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Fri, 21 Aug 2020 13:20:53 -0500 Subject: [PATCH 24/31] relative threshold for contact frame update based on tangential critical force --- src/GRANULAR/pair_granular.cpp | 84 +++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 1b92c2686d..8bf2d1bf17 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -177,6 +177,7 @@ void PairGranular::compute(int eflag, int vflag) double tortwist1, tortwist2, tortwist3; double shrmag,rsht,prjmag; + bool frameupdate; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *history,*allhistory,**firsthistory; @@ -441,22 +442,29 @@ void PairGranular::compute(int eflag, int vflag) // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - // projection - history[0] -= rsht*nx; - history[1] -= rsht*ny; - history[2] -= rsht*nz; - - // also rescale to preserve magnitude - prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - if (prjmag > 0) scalefac = shrmag/prjmag; - else scalefac = 0; - history[0] *= scalefac; - history[1] *= scalefac; - history[2] *= scalefac; + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) + frameupdate = fabs(rsht) < EPSILON*Fscrit; + else + frameupdate = fabs(rsht)*k_tangential < EPSILON*Fscrit; + if (frameupdate) { + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + // projection + history[0] -= rsht*nx; + history[1] -= rsht*ny; + history[2] -= rsht*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; + history[0] *= scalefac; + history[1] *= scalefac; + history[2] *= scalefac; + } // update history if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || @@ -553,39 +561,41 @@ void PairGranular::compute(int eflag, int vflag) int rhist1 = rhist0 + 1; int rhist2 = rhist1 + 1; - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - if (historyupdate) { - // rotate into tangential plane - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - // projection - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - // also rescale to preserve magnitude - prjmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - if (prjmag > 0) scalefac = rollmag/prjmag; - else scalefac = 0; - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][1]; + Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; + if (historyupdate) { + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + frameupdate = fabs(rolldotn)*k_roll < EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // projection + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; + } history[rhist0] += vrl1*dt; history[rhist1] += vrl2*dt; history[rhist2] += vrl3*dt; } - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][1]; fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); if (fr > Frcrit) { From 6fc2ab07ef4797372063c349247877d0d9d02cd8 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Fri, 21 Aug 2020 14:52:39 -0400 Subject: [PATCH 25/31] reset_mol_ids: unique created computes --- src/USER-REACTION/fix_bond_react.cpp | 2 +- src/reset_mol_ids.cpp | 11 ++++++----- src/reset_mol_ids.h | 2 +- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 57bc42cffd..a1154c9221 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -171,7 +171,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg+1],"yes") == 0) { // default delete reset_mol_ids; reset_mol_ids = new ResetMolIDs(lmp); - reset_mol_ids->create_computes(group->names[igroup]); + reset_mol_ids->create_computes(id,group->names[igroup]); iarg += 2; } if (strcmp(arg[iarg+1],"no") == 0) { diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index efeb6a672a..1461a49f32 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg) // create computes - create_computes(groupid); + create_computes(NULL,groupid); // reset molecule IDs @@ -140,22 +140,23 @@ void ResetMolIDs::command(int narg, char **arg) create computes used by reset_mol_ids ------------------------------------------------------------------------- */ -void ResetMolIDs::create_computes(char *groupid) +void ResetMolIDs::create_computes(char *fixid, char *groupid) { int igroup = group->find(groupid); if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID"); groupbit = group->bitmask[igroup]; // create instances of compute fragment/atom, compute reduce (if needed), - // and compute chunk/atom. all use the group-ID for this command + // and compute chunk/atom. all use the group-ID for this command. + // 'fixid' allows for creating independent instances of the computes - idfrag = "reset_mol_ids_FRAGMENT_ATOM"; + idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid); if (singleflag) modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid)); else modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid)); - idchunk = "reset_mol_ids_CHUNK_ATOM"; + idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid); if (compressflag) modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes", idchunk,groupid)); diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h index 62543f6eda..9d2dd24bc8 100644 --- a/src/reset_mol_ids.h +++ b/src/reset_mol_ids.h @@ -30,7 +30,7 @@ class ResetMolIDs : protected Pointers { ResetMolIDs(class LAMMPS *); ~ResetMolIDs(); void command(int, char **); - void create_computes(char *); + void create_computes(char *, char *); void reset(); private: From 01dd80f35e7709bb15fcfeadbf6c6f9ceb960de6 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 23 Aug 2020 11:21:43 -0400 Subject: [PATCH 26/31] bond/react: actually make reset_mol_ids the default --- src/USER-REACTION/fix_bond_react.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index a1154c9221..29aa476148 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -168,21 +168,20 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: " "'reset_mol_ids' keyword has too few arguments"); - if (strcmp(arg[iarg+1],"yes") == 0) { // default - delete reset_mol_ids; - reset_mol_ids = new ResetMolIDs(lmp); - reset_mol_ids->create_computes(id,group->names[igroup]); - iarg += 2; - } - if (strcmp(arg[iarg+1],"no") == 0) { - reset_mol_ids_flag = 0; - iarg += 2; - } + if (strcmp(arg[iarg+1],"yes") == 0) ; // default + if (strcmp(arg[iarg+1],"no") == 0) reset_mol_ids_flag = 0; + iarg += 2; } else if (strcmp(arg[iarg],"react") == 0) { break; } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } + if (reset_mol_ids_flag) { + delete reset_mol_ids; + reset_mol_ids = new ResetMolIDs(lmp); + reset_mol_ids->create_computes(id,group->names[igroup]); + } + // set up common variables as vectors of length 'nreacts' // nevery, cutoff, onemol, twomol, superimpose file @@ -245,7 +244,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) : int n = strlen(arg[iarg]) + 1; if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)"); - strncpy(rxn_name[rxn],arg[iarg++],n); + strcpy(rxn_name[rxn],arg[iarg++]); int groupid = group->find(arg[iarg++]); if (groupid == -1) error->all(FLERR,"Could not find fix group ID"); From 63abb2dff919327995dc350b70f44800176d7eef Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 23 Aug 2020 11:32:54 -0400 Subject: [PATCH 27/31] fix broken reset_mol_ids command --- src/reset_mol_ids.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 1461a49f32..de76ef1732 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg) // create computes - create_computes(NULL,groupid); + create_computes("COMMAND",groupid); // reset molecule IDs From 10080079e3da86a39cc615a236242262710b91a9 Mon Sep 17 00:00:00 2001 From: Jacob Gissinger Date: Sun, 23 Aug 2020 11:44:48 -0400 Subject: [PATCH 28/31] ISO compliance --- src/reset_mol_ids.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index de76ef1732..7ce027a657 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg) // create computes - create_computes("COMMAND",groupid); + create_computes((char *) "COMMAND",groupid); // reset molecule IDs From 024e4c5f21c75189aa58db7e9edab01b0326d7f5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 24 Aug 2020 20:45:44 -0400 Subject: [PATCH 29/31] make formatting and doxygen decorations for utils functions consistent --- src/utils.cpp | 32 ++++++------- src/utils.h | 130 +++++++++++++++++++++++++------------------------- 2 files changed, 80 insertions(+), 82 deletions(-) diff --git a/src/utils.cpp b/src/utils.cpp index 8c19e889cd..f841086673 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -352,7 +352,7 @@ tagint utils::tnumeric(const char *file, int line, const char *str, Return string without leading or trailing whitespace ------------------------------------------------------------------------- */ -std::string utils::trim(const std::string & line) { +std::string utils::trim(const std::string &line) { int beg = re_match(line.c_str(),"\\S+"); int end = re_match(line.c_str(),"\\s+$"); if (beg < 0) beg = 0; @@ -365,7 +365,7 @@ std::string utils::trim(const std::string & line) { Return string without trailing # comment ------------------------------------------------------------------------- */ -std::string utils::trim_comment(const std::string & line) { +std::string utils::trim_comment(const std::string &line) { auto end = line.find_first_of("#"); if (end != std::string::npos) { return line.substr(0, end); @@ -377,7 +377,7 @@ std::string utils::trim_comment(const std::string & line) { return number of words ------------------------------------------------------------------------- */ -size_t utils::count_words(const char * text) { +size_t utils::count_words(const char *text) { size_t count = 0; const char * buf = text; char c = *buf; @@ -406,7 +406,7 @@ size_t utils::count_words(const char * text) { return number of words ------------------------------------------------------------------------- */ -size_t utils::count_words(const std::string & text) { +size_t utils::count_words(const std::string &text) { return utils::count_words(text.c_str()); } @@ -414,7 +414,7 @@ size_t utils::count_words(const std::string & text) { Return number of words ------------------------------------------------------------------------- */ -size_t utils::count_words(const std::string & text, const std::string & separators) { +size_t utils::count_words(const std::string &text, const std::string &separators) { size_t count = 0; size_t start = text.find_first_not_of(separators); @@ -435,7 +435,7 @@ size_t utils::count_words(const std::string & text, const std::string & separato Trim comment from string and return number of words ------------------------------------------------------------------------- */ -size_t utils::trim_and_count_words(const std::string & text, const std::string & separators) { +size_t utils::trim_and_count_words(const std::string &text, const std::string &separators) { return utils::count_words(utils::trim_comment(text), separators); } @@ -526,7 +526,7 @@ std::vector utils::split_words(const std::string &text) Return whether string is a valid integer number ------------------------------------------------------------------------- */ -bool utils::is_integer(const std::string & str) { +bool utils::is_integer(const std::string &str) { if (str.size() == 0) { return false; } @@ -542,7 +542,7 @@ bool utils::is_integer(const std::string & str) { Return whether string is a valid floating-point number ------------------------------------------------------------------------- */ -bool utils::is_double(const std::string & str) { +bool utils::is_double(const std::string &str) { if (str.size() == 0) { return false; } @@ -560,7 +560,7 @@ bool utils::is_double(const std::string & str) { strip off leading part of path, return just the filename ------------------------------------------------------------------------- */ -std::string utils::path_basename(const std::string & path) { +std::string utils::path_basename(const std::string &path) { #if defined(_WIN32) size_t start = path.find_last_of("/\\"); #else @@ -580,7 +580,7 @@ std::string utils::path_basename(const std::string & path) { join two paths ------------------------------------------------------------------------- */ -std::string utils::path_join(const std::string & a, const std::string & b) { +std::string utils::path_join(const std::string &a, const std::string &b) { #if defined(_WIN32) return fmt::format("{}\\{}", a, b); #else @@ -592,7 +592,7 @@ std::string utils::path_join(const std::string & a, const std::string & b) { try to open file for reading ------------------------------------------------------------------------- */ -bool utils::file_is_readable(const std::string & path) { +bool utils::file_is_readable(const std::string &path) { FILE * fp = fopen(path.c_str(), "r"); if(fp) { fclose(fp); @@ -607,7 +607,7 @@ bool utils::file_is_readable(const std::string & path) { specified ------------------------------------------------------------------------- */ -std::string utils::get_potential_file_path(const std::string& path) { +std::string utils::get_potential_file_path(const std::string &path) { std::string filepath = path; std::string filename = utils::path_basename(path); @@ -634,7 +634,7 @@ std::string utils::get_potential_file_path(const std::string& path) { if it has a DATE field, return the following word ------------------------------------------------------------------------- */ -std::string utils::get_potential_date(const std::string & path, const std::string & potential_name) { +std::string utils::get_potential_date(const std::string &path, const std::string &potential_name) { TextFileReader reader(path, potential_name); reader.ignore_comments = false; @@ -657,7 +657,7 @@ std::string utils::get_potential_date(const std::string & path, const std::strin if it has UNITS field, return following word ------------------------------------------------------------------------- */ -std::string utils::get_potential_units(const std::string & path, const std::string & potential_name) { +std::string utils::get_potential_units(const std::string &path, const std::string &potential_name) { TextFileReader reader(path, potential_name); reader.ignore_comments = false; @@ -710,7 +710,7 @@ double utils::get_conversion_factor(const int property, const int conversion) the strings "off" and "unlimited" result in -1.0; ------------------------------------------------------------------------- */ -double utils::timespec2seconds(const std::string & timespec) +double utils::timespec2seconds(const std::string ×pec) { double vals[3]; int i = 0; @@ -728,7 +728,7 @@ double utils::timespec2seconds(const std::string & timespec) if (!values.has_next()) break; vals[i] = values.next_int(); } - } catch (TokenizerException & e) { + } catch (TokenizerException &e) { return -1.0; } diff --git a/src/utils.h b/src/utils.h index 2751870d44..4b9d3f7e81 100644 --- a/src/utils.h +++ b/src/utils.h @@ -29,7 +29,7 @@ namespace LAMMPS_NS { namespace utils { - /** \brief Match text against a simplified regex pattern + /** Match text against a simplified regex pattern * * \param text the text to be matched against the pattern * \param pattern the search pattern, which may contain regexp markers @@ -37,14 +37,14 @@ namespace LAMMPS_NS { */ bool strmatch(const std::string &text, const std::string &pattern); - /** \brief Send message to screen and logfile, if available + /** Send message to screen and logfile, if available * * \param lmp pointer to LAMMPS class instance * \param mesg message to be printed */ void logmesg(LAMMPS *lmp, const std::string &mesg); - /** \brief return a string representing the current system error status + /** return a string representing the current system error status * * This is a wrapper around calling strerror(errno). * @@ -52,7 +52,7 @@ namespace LAMMPS_NS { */ std::string getsyserror(); - /** \brief safe wrapper around fgets() which aborts on errors + /** safe wrapper around fgets() which aborts on errors * or EOF and prints a suitable error message to help debugging * * \param srcname name of the calling source file (from FLERR macro) @@ -66,7 +66,7 @@ namespace LAMMPS_NS { void sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp, const char *filename, Error *error); - /** \brief safe wrapper around fread() which aborts on errors + /** safe wrapper around fread() which aborts on errors * or EOF and prints a suitable error message to help debugging * * \param srcname name of the calling source file (from FLERR macro) @@ -81,7 +81,7 @@ namespace LAMMPS_NS { void sfread(const char *srcname, int srcline, void *s, size_t size, size_t num, FILE *fp, const char *filename, Error *error); - /** \brief Report if a requested style is in a package or may have a typo + /** Report if a requested style is in a package or may have a typo * * \param style type of style that is to be checked for * \param name name of style that was not found @@ -91,8 +91,8 @@ namespace LAMMPS_NS { std::string check_packages_for_style(const std::string &style, const std::string &name, LAMMPS *lmp); - /** \brief Convert a string to a floating point number while checking - if it is a valid floating point or integer number + /** Convert a string to a floating point number while checking + * if it is a valid floating point or integer number * * \param file name of source file for error message * \param line in source file for error message @@ -104,8 +104,8 @@ namespace LAMMPS_NS { double numeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** \brief Convert a string to an integer number while checking - if it is a valid integer number (regular int) + /** Convert a string to an integer number while checking + * if it is a valid integer number (regular int) * * \param file name of source file for error message * \param line in source file for error message @@ -117,8 +117,8 @@ namespace LAMMPS_NS { int inumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** \brief Convert a string to an integer number while checking - if it is a valid integer number (bigint) + /** Convert a string to an integer number while checking + * if it is a valid integer number (bigint) * * \param file name of source file for error message * \param line in source file for error message @@ -130,8 +130,8 @@ namespace LAMMPS_NS { bigint bnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** \brief Convert a string to an integer number while checking - if it is a valid integer number (tagint) + /** Convert a string to an integer number while checking + * if it is a valid integer number (tagint) * * \param file name of source file for error message * \param line in source file for error message @@ -143,54 +143,51 @@ namespace LAMMPS_NS { tagint tnumeric(const char *file, int line, const char *str, bool do_abort, LAMMPS *lmp); - /** - * \brief Trim leading and trailing whitespace. Like TRIM() in Fortran. + /** Trim leading and trailing whitespace. Like TRIM() in Fortran. + * * \param line string that should be trimmed * \return new string without whitespace (string) */ std::string trim(const std::string &line); - /** - * \brief Trim anything from '#' onward + /** Trim anything from '#' onward + * * \param line string that should be trimmed * \return new string without comment (string) */ std::string trim_comment(const std::string &line); - /** - * \brief Count words in string + /** Count words in string + * * \param text string that should be searched * \param separators string containing characters that will be treated as whitespace * \return number of words found */ - size_t count_words(const std::string & text, const std::string & separators); + size_t count_words(const std::string &text, const std::string &separators); - /** - * \brief Count words in string, ignore any whitespace matching " \t\r\n\f" + /** Count words in string, ignore any whitespace matching " \t\r\n\f" + * * \param text string that should be searched - * \param separators string containing characters that will be treated as whitespace * \return number of words found */ - size_t count_words(const std::string & text); + size_t count_words(const std::string &text); - /** - * \brief Count words in C-string, ignore any whitespace matching " \t\r\n\f" + /** Count words in C-string, ignore any whitespace matching " \t\r\n\f" + * * \param text string that should be searched - * \param separators string containing characters that will be treated as whitespace * \return number of words found */ - size_t count_words(const char * text); + size_t count_words(const char *text); - /** - * \brief Count words in a single line, trim anything from '#' onward + /** Count words in a single line, trim anything from '#' onward + * * \param text string that should be trimmed and searched * \param separators string containing characters that will be treated as whitespace * \return number of words found */ - size_t trim_and_count_words(const std::string & text, const std::string & separators = " \t\r\n\f"); + size_t trim_and_count_words(const std::string &text, const std::string &separators = " \t\r\n\f"); - /** - * \brief Take text and split into non-whitespace words. + /** Take text and split into non-whitespace words. * * This can handle single and double quotes, escaped quotes, * and escaped codes within quotes, but due to using an STL @@ -203,21 +200,21 @@ namespace LAMMPS_NS { */ std::vector split_words(const std::string &text); - /** - * \brief Check if string can be converted to valid integer - * \param text string that should be checked + /** Check if string can be converted to valid integer + * + * \param str string that should be checked * \return true, if string contains valid integer, false otherwise */ - bool is_integer(const std::string & str); + bool is_integer(const std::string &str); - /** - * \brief Check if string can be converted to valid floating-point number - * \param text string that should be checked + /** Check if string can be converted to valid floating-point number + * + * \param str string that should be checked * \return true, if string contains valid floating-point number, false otherwise */ - bool is_double(const std::string & str); + bool is_double(const std::string &str); - /** \brief try to detect pathname from FILE pointer. + /** Try to detect pathname from FILE pointer. * * Currently only supported on Linux, otherwise will report "(unknown)". * @@ -228,12 +225,12 @@ namespace LAMMPS_NS { */ const char *guesspath(char *buf, int len, FILE *fp); - /** - * \brief Strip off leading part of path, return just the filename + /** Strip off leading part of path, return just the filename + * * \param path file path * \return file name */ - std::string path_basename(const std::string & path); + std::string path_basename(const std::string &path); /** * \brief Join two paths @@ -241,64 +238,65 @@ namespace LAMMPS_NS { * \param b second path * \return combined path */ - std::string path_join(const std::string & a, const std::string & b); + std::string path_join(const std::string &a, const std::string &b); /** * \brief Check if file exists and is readable * \param path file path * \return true if file exists and is readable */ - bool file_is_readable(const std::string & path); + bool file_is_readable(const std::string &path); - /** - * \brief Determine full path of potential file - * If file is not found in current directory, search LAMMPS_POTENTIALS folder + /** Determine full path of potential file. If file is not found in current directory, + * search LAMMPS_POTENTIALS folder + * * \param path file path * \return full path to potential file */ - std::string get_potential_file_path(const std::string& path); + std::string get_potential_file_path(const std::string &path); - /** - * \brief Read potential file and return DATE field if it is present + /** Read potential file and return DATE field if it is present + * * \param path file path * \param potential_name name of potential that is being read * \return DATE field if present */ - std::string get_potential_date(const std::string & path, const std::string & potential_name); + std::string get_potential_date(const std::string &path, const std::string &potential_name); - /** - * \brief Read potential file and return UNITS field if it is present + /** Read potential file and return UNITS field if it is present + * * \param path file path * \param potential_name name of potential that is being read * \return UNITS field if present */ - std::string get_potential_units(const std::string & path, const std::string & potential_name); + std::string get_potential_units(const std::string &path, const std::string &potential_name); enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1<<1 }; enum { UNKNOWN = 0, ENERGY }; - /** - * \brief Return bitmask of available conversion factors for a given propert + /** Return bitmask of available conversion factors for a given property + * * \param property property to be converted * \return bitmask indicating available conversions */ int get_supported_conversions(const int property); - /** - * \brief Return unit conversion factor for given property and selected from/to units + /** Return unit conversion factor for given property and selected from/to units + * * \param property property to be converted * \param conversion constant indicating the conversion * \return conversion factor */ double get_conversion_factor(const int property, const int conversion); - /** - * \brief Convert a time string to seconds - * The strings "off" and "unlimited" result in -1 + /** Convert a time string to seconds + * + * The strings "off" and "unlimited" result in -1 + * * \param timespec a string in the following format: ([[HH:]MM:]SS) * \return total in seconds */ - double timespec2seconds(const std::string & timespec); + double timespec2seconds(const std::string ×pec); } } From 14b66d1f84752addf96bfec1e500bc7cf62205e2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 25 Aug 2020 12:15:19 -0400 Subject: [PATCH 30/31] tweak test tolerance of reax/c tests for running on ubuntu 18.04 --- unittest/force-styles/tests/atomic-pair-reax_c.yaml | 2 +- unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/unittest/force-styles/tests/atomic-pair-reax_c.yaml b/unittest/force-styles/tests/atomic-pair-reax_c.yaml index a1fea19d79..15f7bf02ab 100644 --- a/unittest/force-styles/tests/atomic-pair-reax_c.yaml +++ b/unittest/force-styles/tests/atomic-pair-reax_c.yaml @@ -1,7 +1,7 @@ --- lammps_version: 21 Aug 2020 date_generated: Sun Aug 23 06:35:54 202 -epsilon: 5e-12 +epsilon: 7.5e-12 prerequisites: ! | pair reax/c fix qeq/reax diff --git a/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml b/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml index 1383339444..7d549efc18 100644 --- a/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml +++ b/unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml @@ -1,7 +1,7 @@ --- lammps_version: 21 Aug 2020 date_generated: Sun Aug 23 06:41:04 202 -epsilon: 5e-11 +epsilon: 1e-10 prerequisites: ! | pair reax/c fix qeq/reax From f89a0f9fe35a0a5567981f7a7de0f796ce61d38f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 26 Aug 2020 11:50:20 -0400 Subject: [PATCH 31/31] must not try to delete computes if they have not been created and their ids not yet set --- src/reset_mol_ids.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 7ce027a657..a6cfb69185 100644 --- a/src/reset_mol_ids.cpp +++ b/src/reset_mol_ids.cpp @@ -41,14 +41,17 @@ ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) { compressflag = 1; singleflag = 0; offset = -1; + + idfrag.clear(); + idchunk.clear(); } /* ---------------------------------------------------------------------- */ ResetMolIDs::~ResetMolIDs() { - modify->delete_compute(idfrag); - if (compressflag) modify->delete_compute(idchunk); + if (!idfrag.empty()) modify->delete_compute(idfrag); + if (compressflag && !idchunk.empty()) modify->delete_compute(idchunk); } /* ----------------------------------------------------------------------