tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

restore user facing keyword back to "threebody" defaulting to "on"

Browse Source 
also some minor updates:
- streamline description in the documentation, add links/references
- print message when disabling threebody terms
- improve error messages, simplify argument processing
This commit is contained in:
Axel Kohlmeyer
2022-07-02 17:08:22 -04:00
parent f826d175f0
commit f38a417a32
5 changed files with 75 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

91
doc/src/pair_sw.rst
View File

@ -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) <Jiang2>` and :ref:`(Jiang3) <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_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`.
----------

3
doc/src/pair_sw_angle_table.rst
View File

@ -296,7 +296,8 @@ for pair interactions.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style sw <pair_sw>`, :doc:`pair_style threebody/table <pair_threebody_table>`
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style sw <pair_sw>`,
:doc:`pair_style threebody/table <pair_threebody_table>`
----------

29
src/MANYBODY/pair_sw.cpp
View File

@ -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(&params[ijparam],&params[ikparam],&params[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 ||

4
src/MANYBODY/pair_sw_angle_table.cpp
View File

@ -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");
}
if (narg != 0) error->all(FLERR,"Illegal pair_style sw/angle/table command");
}

11
src/MANYBODY/pair_sw_mod.cpp
View File

@ -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);
}
/* ---------------------------------------------------------------------- */