diff --git a/doc/src/pair_sw.rst b/doc/src/pair_sw.rst index a83f446a07..b2305dab5e 100644 --- a/doc/src/pair_sw.rst +++ b/doc/src/pair_sw.rst @@ -24,17 +24,16 @@ Syntax pair_style style keyword values * style = *sw* or *sw/mod* -* keyword = *maxdelcs* or *skip_threebody* +* keyword = *maxdelcs* or *threebody* .. parsed-literal:: - *maxdelcs* value = delta1 delta2 (optional) + *maxdelcs* value = delta1 delta2 (optional, sw/mod only) delta1 = The minimum thershold for the variation of cosine of three-body angle delta2 = The maximum threshold for the variation of cosine of three-body angle - *skip_threebody* value = *on* or *off* (optional) - off (default) = Compute both the three-body and two-body terms of the potential - on = Compute only the two-body term of the potential - + *threebody* value = *on* or *off* (optional, sw only) + on (default) = Compute both the three-body and two-body terms of the potential + off = Compute only the two-body term of the potential Examples """""""" @@ -48,7 +47,7 @@ Examples pair_style sw/mod maxdelcs 0.25 0.35 pair_coeff * * tmd.sw.mod Mo S S - pair_style hybrid sw sw skip_threebody on + pair_style hybrid sw threebody on sw threebody off pair_coeff * * sw 1 mW_xL.sw mW NULL pair_coeff 1 2 sw 2 mW_xL.sw mW xL pair_coeff 2 2 sw 2 mW_xL.sw mW xL @@ -77,22 +76,25 @@ three-body term. The summations in the formula are over all neighbors J and K of atom I within a cutoff distance :math:`a `\sigma`. The *sw/mod* style is designed for simulations of materials when -distinguishing three-body angles are necessary, such as borophene -and transition metal dichalcogenides, which cannot be described -by the original code for the Stillinger-Weber potential. -For instance, there are several types of angles around each Mo atom in `MoS_2`, -and some unnecessary angle types should be excluded in the three-body interaction. -Such exclusion may be realized by selecting proper angle types directly. -The exclusion of unnecessary angles is achieved here by the cut-off function (`f_C(\delta)`), -which induces only minimum modifications for LAMMPS. +distinguishing three-body angles are necessary, such as borophene and +transition metal dichalcogenides, which cannot be described by the +original code for the Stillinger-Weber potential. For instance, there +are several types of angles around each Mo atom in `MoS_2`, and some +unnecessary angle types should be excluded in the three-body +interaction. Such exclusion may be realized by selecting proper angle +types directly. The exclusion of unnecessary angles is achieved here by +the cut-off function (`f_C(\delta)`), which induces only minimum +modifications for LAMMPS. Validation, benchmark tests, and applications of the *sw/mod* style can be found in :ref:`(Jiang2) ` and :ref:`(Jiang3) `. -The *sw/mod* style computes the energy E of a system of atoms, whose potential -function is mostly the same as the Stillinger-Weber potential. The only modification -is in the three-body term, where the value of :math:`\delta = \cos \theta_{ijk} - \cos \theta_{0ijk}` -used in the original energy and force expression is scaled by a switching factor :math:`f_C(\delta)`: +The *sw/mod* style computes the energy E of a system of atoms, whose +potential function is mostly the same as the Stillinger-Weber +potential. The only modification is in the three-body term, where the +value of :math:`\delta = \cos \theta_{ijk} - \cos \theta_{0ijk}` used in +the original energy and force expression is scaled by a switching factor +:math:`f_C(\delta)`: .. math:: @@ -103,28 +105,38 @@ used in the original energy and force expression is scaled by a switching factor 0 & \left| \delta \right| > \delta_2 \end{array} \right. \\ -This cut-off function decreases smoothly from 1 to 0 over the range :math:`[\delta_1, \delta_2]`. -This smoothly turns off the energy and force contributions for :math:`\left| \delta \right| > \delta_2`. -It is suggested that :math:`\delta 1` and :math:`\delta_2` to be the value around -:math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with -:math:`\theta_1` and :math:`\theta_2` as the different types of angles around an atom. -For borophene and transition metal dichalcogenides, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`. -This value enables the cut-off function to exclude unnecessary angles in the three-body SW terms. +This cut-off function decreases smoothly from 1 to 0 over the range +:math:`[\delta_1, \delta_2]`. This smoothly turns off the energy and +force contributions for :math:`\left| \delta \right| > \delta_2`. It is +suggested that :math:`\delta 1` and :math:`\delta_2` to be the value +around :math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with +:math:`\theta_1` and :math:`\theta_2` as the different types of angles +around an atom. For borophene and transition metal dichalcogenides, +:math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`. This value enables +the cut-off function to exclude unnecessary angles in the three-body SW +terms. .. note:: - The cut-off function is just to be used as a technique to exclude some unnecessary angles, - and it has no physical meaning. It should be noted that the force and potential are inconsistent - with each other in the decaying range of the cut-off function, as the angle dependence for the - cut-off function is not implemented in the force (first derivation of potential). - However, the angle variation is much smaller than the given threshold value for actual simulations, - so the inconsistency between potential and force can be neglected in actual simulations. + The cut-off function is just to be used as a technique to exclude + some unnecessary angles, and it has no physical meaning. It should be + noted that the force and potential are inconsistent with each other + in the decaying range of the cut-off function, as the angle + dependence for the cut-off function is not implemented in the force + (first derivation of potential). However, the angle variation is + much smaller than the given threshold value for actual simulations, + so the inconsistency between potential and force can be neglected in + actual simulations. -The *skip_threebody* keyword determines whether or not the three-body term of the potential is calculated. -Skipping this significantly increases the speed of the calculation, with the tradeoff that :math:\lambda_{ijk} -is forcibly set to 0. If the keyword is used with the pair styles, sw/gpu, sw/intel, or sw/kokkos, -:math:\lambda_{ijk} will still be forcibly set to 0, but no additional benefits will be gained. This keyword -cannot be used for variants of the sw/mod or sw/angle/table pair styles. +The *threebody* keyword is optional and determines whether or not the +three-body term of the potential is calculated. The default value is +"on" and it is only available for the plain *sw* pair style variants, +but not available for the *sw/mod* and :doc:`sw/angle/table +` pair style variants. To turn off the threebody +contributions all :math:`\lambda_{ijk}` parameters from the potential +file are forcibly set to 0. In addition the pair style implementations +may employ code optimizations for the *threebody off* setting that can +result in significant speedups versus the default. Only a single pair_coeff command is used with the *sw* and *sw/mod* styles which specifies a Stillinger-Weber potential file with parameters for all @@ -297,8 +309,9 @@ Related commands Default """"""" -The default values for the *maxdelcs* setting of the *sw/mod* pair -style are *delta1* = 0.25 and *delta2* = 0.35`. +The default value for the *threebody* setting of the "sw" pair style is +"on", the default values for the "*maxdelcs* setting of the *sw/mod* +pair style are *delta1* = 0.25 and *delta2* = 0.35`. ---------- diff --git a/doc/src/pair_sw_angle_table.rst b/doc/src/pair_sw_angle_table.rst index ff88711d7a..65e73d4445 100644 --- a/doc/src/pair_sw_angle_table.rst +++ b/doc/src/pair_sw_angle_table.rst @@ -296,7 +296,8 @@ for pair interactions. Related commands """""""""""""""" -:doc:`pair_coeff `, :doc:`pair_style sw `, :doc:`pair_style threebody/table ` +:doc:`pair_coeff `, :doc:`pair_style sw `, +:doc:`pair_style threebody/table ` ---------- diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index 3cc402cb5f..b615480754 100644 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -44,6 +44,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp) manybody_flag = 1; centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); + skip_threebody_flag = false; params = nullptr; @@ -172,24 +173,24 @@ void PairSW::compute(int eflag, int vflag) delr1[1] = x[j][1] - ytmp; delr1[2] = x[j][2] - ztmp; rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - + double fjxtmp,fjytmp,fjztmp; fjxtmp = fjytmp = fjztmp = 0.0; - + for (kk = jj+1; kk < numshort; kk++) { k = neighshort[kk]; ktype = map[type[k]]; ikparam = elem3param[itype][ktype][ktype]; ijkparam = elem3param[itype][jtype][ktype]; - + delr2[0] = x[k][0] - xtmp; delr2[1] = x[k][1] - ytmp; delr2[2] = x[k][2] - ztmp; rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - + threebody(¶ms[ijparam],¶ms[ikparam],¶ms[ijkparam], rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); - + fxtmp -= fj[0] + fk[0]; fytmp -= fj[1] + fk[1]; fztmp -= fj[2] + fk[2]; @@ -199,7 +200,7 @@ void PairSW::compute(int eflag, int vflag) f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; - + if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2); } f[j][0] += fjxtmp; @@ -233,18 +234,15 @@ void PairSW::allocate() void PairSW::settings(int narg, char ** arg) { - // Default - skip_threebody_flag = false; - + // process optional keywords int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"skip_threebody") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style command"); - skip_threebody_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg],"threebody") == 0) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style sw", error); + skip_threebody_flag = !utils::logical(FLERR,arg[iarg+1],false,lmp); one_coeff = !skip_threebody_flag; iarg += 2; - } else error->all(FLERR,"Illegal pair_style command"); + } else error->all(FLERR, "Illegal pair_style sw keyword: {}", arg[iarg]); } } @@ -305,6 +303,8 @@ void PairSW::read_file(char *file) PotentialFileReader reader(lmp, file, "sw", unit_convert_flag); char *line; + if (skip_threebody_flag) utils::logmesg(lmp, " disabling sw potential three-body terms\n"); + // transparently convert units for supported conversions int unit_convert = reader.get_unit_convert(); @@ -369,6 +369,7 @@ void PairSW::read_file(char *file) params[nparams].epsilon *= conversion_factor; } + // turn off three-body term if (skip_threebody_flag) params[nparams].lambda = 0; if (params[nparams].epsilon < 0.0 || params[nparams].sigma < 0.0 || diff --git a/src/MANYBODY/pair_sw_angle_table.cpp b/src/MANYBODY/pair_sw_angle_table.cpp index 7667f77491..5b9339780f 100644 --- a/src/MANYBODY/pair_sw_angle_table.cpp +++ b/src/MANYBODY/pair_sw_angle_table.cpp @@ -759,5 +759,5 @@ void PairSWAngleTable::uf_lookup(ParamTable *pm, double x, double &u, double &f) void PairSWAngleTable::settings(int narg, char **/*arg*/) { - if (narg != 0) error->all(FLERR,"Illegal pair_style command"); -} \ No newline at end of file + if (narg != 0) error->all(FLERR,"Illegal pair_style sw/angle/table command"); +} diff --git a/src/MANYBODY/pair_sw_mod.cpp b/src/MANYBODY/pair_sw_mod.cpp index ce24952fc7..e6d17b0733 100644 --- a/src/MANYBODY/pair_sw_mod.cpp +++ b/src/MANYBODY/pair_sw_mod.cpp @@ -41,21 +41,18 @@ PairSWMOD::PairSWMOD(LAMMPS *lmp) : PairSW(lmp) void PairSWMOD::settings(int narg, char **arg) { - // process optional keywords - + // process optional keywords (and not (yet) optional keywords from parent class). int iarg = 0; - while (iarg < narg) { if (strcmp(arg[iarg],"maxdelcs") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style sw/mod", error); delta1 = utils::numeric(FLERR,arg[iarg+1],false,lmp); delta2 = utils::numeric(FLERR,arg[iarg+2],false,lmp); iarg += 3; if ((delta1 < 0.0) || (delta1 > 1.0) || (delta2 < 0.0) || (delta2 > 1.0) || (delta1 > delta2)) - error->all(FLERR,"Illegal values for maxdelcs keyword"); - } else error->all(FLERR,"Illegal pair_style command"); + error->all(FLERR, "Out of range value(s) for pair style sw/mod maxdelcs keyword"); + } else error->all(FLERR, "Illegal pair_style sw/mod keyword: {}", arg[iarg]); } - PairSW::settings(narg-iarg,arg+iarg); } /* ---------------------------------------------------------------------- */