Merge pull request #4262 from jrgissing/type-labels-for-pair_style_commands

Type label support for pair_style commands
This commit is contained in:
Axel Kohlmeyer
2024-08-19 15:02:23 -04:00
committed by GitHub
34 changed files with 404 additions and 144 deletions

View File

@ -25,7 +25,7 @@ There are two ways to implement TIP4P water in LAMMPS:
This can be done with the following pair styles for Coulomb with a cutoff: This can be done with the following pair styles for Coulomb with a cutoff:
* :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>` * :doc:`pair_style tip4p/cut <pair_coul>`
* :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>` * :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
or these commands for a long-range Coulomb treatment: or these commands for a long-range Coulomb treatment:

View File

@ -70,8 +70,8 @@ above, or in the restart files read by the
* :math:`K` = atom type of the third atom (1 to :math:`N_{\text{types}}`) * :math:`K` = atom type of the third atom (1 to :math:`N_{\text{types}}`)
* :math:`\nu` = prefactor (energy/distance\^9 units) * :math:`\nu` = prefactor (energy/distance\^9 units)
:math:`K` can be specified in one of two ways. An explicit numeric value can :math:`K` can be specified in one of two ways. An explicit numeric value
be used, as in the second example above. :math:`J \leq K` is required. LAMMPS or type label can be used, as in the second example above. LAMMPS
sets the coefficients for the other 5 symmetric interactions to the same sets the coefficients for the other 5 symmetric interactions to the same
values. E.g. if :math:`I = 1`, :math:`J = 2`, :math:`K = 3`, then these 6 values. E.g. if :math:`I = 1`, :math:`J = 2`, :math:`K = 3`, then these 6
values are set to the specified :math:`\nu`: :math:`\nu_{123}`, values are set to the specified :math:`\nu`: :math:`\nu_{123}`,

View File

@ -93,12 +93,19 @@ Syntax
pair_style coul/long cutoff pair_style coul/long cutoff
pair_style coul/wolf alpha cutoff pair_style coul/wolf alpha cutoff
pair_style coul/streitz cutoff keyword alpha pair_style coul/streitz cutoff keyword alpha
* cutoff = global cutoff for Coulombic interactions
* kappa = Debye length (inverse distance units)
* alpha = damping parameter (inverse distance units)
.. code-block:: LAMMPS
pair_style tip4p/cut otype htype btype atype qdist cutoff pair_style tip4p/cut otype htype btype atype qdist cutoff
pair_style tip4p/long otype htype btype atype qdist cutoff pair_style tip4p/long otype htype btype atype qdist cutoff
* cutoff = global cutoff for Coulombic interactions * otype,htype = atom types (numeric or type label) for TIP4P O and H
* kappa = Debye length (inverse distance units) * btype,atype = bond and angle types (numeric or type label) for TIP4P waters
* alpha = damping parameter (inverse distance units) * qdist = distance from O atom to massless charge (distance units)
Examples Examples
"""""""" """"""""
@ -138,6 +145,12 @@ Examples
pair_style tip4p/long 1 2 7 8 0.15 10.0 pair_style tip4p/long 1 2 7 8 0.15 10.0
pair_coeff * * pair_coeff * *
pair_style tip4p/cut OW HW HW-OW HW-OW-HW 0.15 12.0
labelmap atom 1 OW 2 HW
labelmap bond 1 HW-OW
labelmap angle 1 HW-OW-HW
pair_coeff * *
Description Description
""""""""""" """""""""""
@ -307,6 +320,11 @@ Coulombic solver (Ewald or PPPM).
atom. For example, if the atom ID of an O atom in a TIP4P water atom. For example, if the atom ID of an O atom in a TIP4P water
molecule is 500, then its 2 H atoms must have IDs 501 and 502. molecule is 500, then its 2 H atoms must have IDs 501 and 502.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
See the :doc:`Howto tip4p <Howto_tip4p>` page for more information See the :doc:`Howto tip4p <Howto_tip4p>` page for more information
on how to use the TIP4P pair styles and lists of parameters to set. on how to use the TIP4P pair styles and lists of parameters to set.
Note that the neighbor list cutoff for Coulomb interactions is Note that the neighbor list cutoff for Coulomb interactions is
@ -381,9 +399,10 @@ Restrictions
"""""""""""" """"""""""""
The *coul/long*, *coul/msm*, *coul/streitz*, and *tip4p/long* styles are The *coul/long*, *coul/msm*, *coul/streitz*, and *tip4p/long* styles are
part of the KSPACE package. The *coul/cut/global* and *coul/exclude* are part of the KSPACE package. The *coul/cut/global*, *coul/exclude* styles are
part of the EXTRA-PAIR package. A pair style is only enabled if LAMMPS was part of the EXTRA-PAIR package. The *tip4p/cut* style is part of the MOLECULE
built with its corresponding package. See the :doc:`Build package <Build_package>` package. A pair style is only enabled if LAMMPS was built with its
corresponding package. See the :doc:`Build package <Build_package>`
doc page for more info. doc page for more info.
Related commands Related commands

View File

@ -10,7 +10,7 @@ Syntax
pair_style e3b Otype pair_style e3b Otype
* Otype = atom type for oxygen * Otype = atom type (numeric or type label) for oxygen
.. code-block:: LAMMPS .. code-block:: LAMMPS
@ -50,6 +50,10 @@ Examples
pair_style hybrid/overlay e3b 1 lj/cut/tip4p/long 1 2 1 1 0.15 8.5 pair_style hybrid/overlay e3b 1 lj/cut/tip4p/long 1 2 1 1 0.15 8.5
pair_coeff * * e3b preset 2011 pair_coeff * * e3b preset 2011
pair_style e3b OW
labelmap atom 1 C 2 H 3 O 4 N 5 OW 6 HW
pair_coeff * * Ea 35.85 Eb -240.2 Ec 449.3 E2 108269.9 K3 1.907 K2 4.872 Rc3 5.2 Rc2 5.2 Rs 5.0 bondL 0.9572
Used in example input script: Used in example input script:
.. parsed-literal:: .. parsed-literal::
@ -102,6 +106,11 @@ atom type is inferred from the ordering of the atoms.
Each water molecule must have consecutive IDs with the oxygen first. Each water molecule must have consecutive IDs with the oxygen first.
This pair style does not test that this criteria is met. This pair style does not test that this criteria is met.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
The *pair_coeff* command must have at least one keyword/value pair, as The *pair_coeff* command must have at least one keyword/value pair, as
described above. The *preset* keyword sets the potential parameters to described above. The *preset* keyword sets the potential parameters to
the values used in :ref:`(Tainter 2011) <Tainter2011>` or the values used in :ref:`(Tainter 2011) <Tainter2011>` or

View File

@ -97,8 +97,8 @@ Syntax
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
*lj/cut/tip4p/long/soft* args = otype htype btype atype qdist n alpha_LJ alpha_C cutoff (cutoff2) *lj/cut/tip4p/long/soft* args = otype htype btype atype qdist n alpha_LJ alpha_C cutoff (cutoff2)
otype,htype = atom types for TIP4P O and H otype,htype = atom types (numeric or type label) for TIP4P O and H
btype,atype = bond and angle types for TIP4P waters btype,atype = bond and angle types (numeric or type label) for TIP4P waters
qdist = distance from O atom to massless charge (distance units) qdist = distance from O atom to massless charge (distance units)
n, alpha_LJ, alpha_C = parameters of the soft-core potential n, alpha_LJ, alpha_C = parameters of the soft-core potential
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
@ -125,8 +125,8 @@ Syntax
n, alpha_C = parameters of the soft-core potential n, alpha_C = parameters of the soft-core potential
cutoff = global cutoff for Coulomb interactions (distance units) cutoff = global cutoff for Coulomb interactions (distance units)
*tip4p/long/soft* args = otype htype btype atype qdist n alpha_C cutoff *tip4p/long/soft* args = otype htype btype atype qdist n alpha_C cutoff
otype,htype = atom types for TIP4P O and H otype,htype = atom types (numeric or type label) for TIP4P O and H
btype,atype = bond and angle types for TIP4P waters btype,atype = bond and angle types (numeric or type label) for TIP4P waters
qdist = distance from O atom to massless charge (distance units) qdist = distance from O atom to massless charge (distance units)
n, alpha_C = parameters of the soft-core potential n, alpha_C = parameters of the soft-core potential
cutoff = global cutoff for Coulomb interactions (distance units) cutoff = global cutoff for Coulomb interactions (distance units)
@ -161,6 +161,13 @@ Examples
pair_coeff * * 0.155 3.1536 1.0 pair_coeff * * 0.155 3.1536 1.0
pair_coeff 1 1 0.155 3.1536 1.0 9.5 pair_coeff 1 1 0.155 3.1536 1.0 9.5
pair_style lj/cut/tip4p/long/soft OW HW HW-OW HW-OW-HW 0.15 2.0 0.5 10.0 9.8
labelmap atom 1 OW 2 HW
labelmap bond 1 HW-OW
labelmap angle 1 HW-OW-HW
pair_coeff * * 0.155 3.1536 1.0
pair_coeff OW OW 0.155 3.1536 1.0 9.5
pair_style lj/charmm/coul/long 2.0 0.5 10.0 8.0 10.0 pair_style lj/charmm/coul/long 2.0 0.5 10.0 8.0 10.0
pair_style lj/charmm/coul/long 2.0 0.5 10.0 8.0 10.0 9.0 pair_style lj/charmm/coul/long 2.0 0.5 10.0 8.0 10.0 9.0
pair_coeff * * 0.28 3.1 1.0 pair_coeff * * 0.28 3.1 1.0
@ -275,6 +282,11 @@ TIP4P water model and before the cutoffs. The activation parameter lambda is
supplied as an argument of the :doc:`pair_coeff <pair_coeff>` command, after supplied as an argument of the :doc:`pair_coeff <pair_coeff>` command, after
epsilon and sigma and before the optional cutoffs. epsilon and sigma and before the optional cutoffs.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
Style *lj/charmm/coul/long/soft* implements a soft-core version of the modified Style *lj/charmm/coul/long/soft* implements a soft-core version of the modified
12-6 LJ potential used in CHARMM and documented in the :doc:`pair_style 12-6 LJ potential used in CHARMM and documented in the :doc:`pair_style
lj/charmm/coul/long <pair_charmm>` style. In the soft version the parameters lj/charmm/coul/long <pair_charmm>` style. In the soft version the parameters

View File

@ -38,6 +38,9 @@ Examples
pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse 2 9.0 11.0 90 pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse 2 9.0 11.0 90
pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9 11 90 pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9 11 90
labelmap atom 1 C 2 O 3 H
pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9 11 90
Description Description
""""""""""" """""""""""
@ -144,7 +147,7 @@ in the examples above.
For the *hbond/dreiding/lj* style the list of coefficients is as For the *hbond/dreiding/lj* style the list of coefficients is as
follows: follows:
* K = hydrogen atom type = 1 to Ntypes * K = hydrogen atom type = 1 to Ntypes, or type label
* donor flag = *i* or *j* * donor flag = *i* or *j*
* :math:`\epsilon` (energy units) * :math:`\epsilon` (energy units)
* :math:`\sigma` (distance units) * :math:`\sigma` (distance units)
@ -156,7 +159,7 @@ follows:
For the *hbond/dreiding/morse* style the list of coefficients is as For the *hbond/dreiding/morse* style the list of coefficients is as
follows: follows:
* K = hydrogen atom type = 1 to Ntypes * K = hydrogen atom type = 1 to Ntypes, or type label
* donor flag = *i* or *j* * donor flag = *i* or *j*
* :math:`D_0` (energy units) * :math:`D_0` (energy units)
* :math:`\alpha` (1/distance units) * :math:`\alpha` (1/distance units)
@ -240,7 +243,10 @@ heading) the following commands could be included in an input script:
Restrictions Restrictions
"""""""""""" """"""""""""
none
This pair style can only be used if LAMMPS was built with the
MOLECULE package. See the :doc:`Build package <Build_package>` doc page
for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""

View File

@ -28,14 +28,14 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*lj/cut/tip4p/cut* args = otype htype btype atype qdist cutoff (cutoff2) *lj/cut/tip4p/cut* args = otype htype btype atype qdist cutoff (cutoff2)
otype,htype = atom types for TIP4P O and H otype,htype = atom types (numeric or type label) for TIP4P O and H
btype,atype = bond and angle types for TIP4P waters btype,atype = bond and angle types (numeric or type label) for TIP4P waters
qdist = distance from O atom to massless charge (distance units) qdist = distance from O atom to massless charge (distance units)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
*lj/cut/tip4p/long* args = otype htype btype atype qdist cutoff (cutoff2) *lj/cut/tip4p/long* args = otype htype btype atype qdist cutoff (cutoff2)
otype,htype = atom types for TIP4P O and H otype,htype = atom types (numeric or type label) for TIP4P O and H
btype,atype = bond and angle types for TIP4P waters btype,atype = bond and angle types (numeric or type label) for TIP4P waters
qdist = distance from O atom to massless charge (distance units) qdist = distance from O atom to massless charge (distance units)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
@ -55,6 +55,13 @@ Examples
pair_coeff * * 100.0 3.0 pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0 pair_coeff 1 1 100.0 3.5 9.0
pair_style lj/cut/tip4p/long OW HW HW-OW HW-OW-HW 0.15 12.0
labelmap atom 1 OW 2 HW
labelmap bond 1 HW-OW
labelmap angle 1 HW-OW-HW
pair_coeff * * 100.0 3.0
pair_coeff OW OW 100.0 3.5 9.0
Description Description
""""""""""" """""""""""
@ -81,6 +88,11 @@ with a long-range Coulombic solver (Ewald or PPPM).
atom. For example, if the atom ID of an O atom in a TIP4P water atom. For example, if the atom ID of an O atom in a TIP4P water
molecule is 500, then its 2 H atoms must have IDs 501 and 502. molecule is 500, then its 2 H atoms must have IDs 501 and 502.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
See the :doc:`Howto tip4p <Howto_tip4p>` page for more information See the :doc:`Howto tip4p <Howto_tip4p>` page for more information
on how to use the TIP4P pair styles and lists of parameters to set. on how to use the TIP4P pair styles and lists of parameters to set.
Note that the neighbor list cutoff for Coulomb interactions is Note that the neighbor list cutoff for Coulomb interactions is

View File

@ -44,8 +44,8 @@ Syntax
flag_coul = *long* or *off* flag_coul = *long* or *off*
*long* = use Kspace long-range summation for Coulombic 1/r term *long* = use Kspace long-range summation for Coulombic 1/r term
*off* = omit Coulombic term *off* = omit Coulombic term
otype,htype = atom types for TIP4P O and H otype,htype = atom types (numeric or type label) for TIP4P O and H
btype,atype = bond and angle types for TIP4P waters btype,atype = bond and angle types (numeric or type label) for TIP4P waters
qdist = distance from O atom to massless charge (distance units) qdist = distance from O atom to massless charge (distance units)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
@ -66,6 +66,13 @@ Examples
pair_coeff * * 100.0 3.0 pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0 pair_coeff 1 1 100.0 3.5 9.0
pair_style lj/long/tip4p/long long long OW HW HW-OW HW-OW-HW 0.15 12.0
labelmap atom 1 OW 2 HW
labelmap bond 1 HW-OW
labelmap angle 1 HW-OW-HW
pair_coeff * * 100.0 3.0
pair_coeff OW OW 100.0 3.5 9.0
Description Description
""""""""""" """""""""""
@ -113,6 +120,11 @@ massless charge site are specified as pair_style arguments.
atom. For example, if the atom ID of an O atom in a TIP4P water atom. For example, if the atom ID of an O atom in a TIP4P water
molecule is 500, then its 2 H atoms must have IDs 501 and 502. molecule is 500, then its 2 H atoms must have IDs 501 and 502.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
See the :doc:`Howto tip4p <Howto_tip4p>` page for more See the :doc:`Howto tip4p <Howto_tip4p>` page for more
information on how to use the TIP4P pair style. Note that the information on how to use the TIP4P pair style. Note that the
neighbor list cutoff for Coulomb interactions is effectively extended neighbor list cutoff for Coulomb interactions is effectively extended

View File

@ -15,7 +15,7 @@ Syntax
pair_style srp/react cutoff btype dist react-id keyword value ... pair_style srp/react cutoff btype dist react-id keyword value ...
* cutoff = global cutoff for SRP interactions (distance units) * cutoff = global cutoff for SRP interactions (distance units)
* btype = bond type to apply SRP interactions to (can be wildcard, see below) * btype = bond type (numeric, type label, or wildcard) to apply SRP interactions to
* distance = *min* or *mid* * distance = *min* or *mid*
* react-id = id of either fix bond/break or fix bond/create * react-id = id of either fix bond/break or fix bond/create
* zero or more keyword/value pairs may be appended * zero or more keyword/value pairs may be appended
@ -23,7 +23,7 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*bptype* value = atom type for bond particles *bptype* value = atom type (numeric or type label) for bond particles
*exclude* value = *yes* or *no* *exclude* value = *yes* or *no*
Examples Examples
@ -52,6 +52,10 @@ Examples
pair_coeff 1 2 none pair_coeff 1 2 none
pair_coeff 2 2 srp 100.0 0.8 pair_coeff 2 2 srp 100.0 0.8
labelmap bond 1 C-C
pair_style hybrid srp 0.8 C-C mid
Description Description
@ -117,6 +121,11 @@ is used.
between atom type *bptype* and all other types of atoms. An error between atom type *bptype* and all other types of atoms. An error
will be flagged if :doc:`pair_style hybrid <pair_hybrid>` is not used. will be flagged if :doc:`pair_style hybrid <pair_hybrid>` is not used.
.. note::
If using type labels, the type labels must be defined before calling
the :doc:`pair_coeff <pair_coeff>` command.
The optional *exclude* keyword determines if forces are computed The optional *exclude* keyword determines if forces are computed
between first neighbor (directly connected) bonds. For a setting of between first neighbor (directly connected) bonds. For a setting of
*no*, first neighbor forces are computed; for *yes* they are not *no*, first neighbor forces are computed; for *yes* they are not
@ -207,4 +216,3 @@ Chem Phys, 136 (13) 134903, 2012.
.. _Palkar: .. _Palkar:
**(Palkar)** Palkar V, Kuksenok O, J. Phys. Chem. B, 126 (1), 336-346, 2022 **(Palkar)** Palkar V, Kuksenok O, J. Phys. Chem. B, 126 (1), 336-346, 2022

View File

@ -10,7 +10,7 @@ Syntax
pair_write itype jtype N style inner outer file keyword Qi Qj pair_write itype jtype N style inner outer file keyword Qi Qj
* itype,jtype = 2 atom types * itype,jtype = 2 atom types (numeric or type label)
* N = # of values * N = # of values
* style = *r* or *rsq* or *bitmap* * style = *r* or *rsq* or *bitmap*
* inner,outer = inner and outer cutoff (distance units) * inner,outer = inner and outer cutoff (distance units)
@ -26,6 +26,9 @@ Examples
pair_write 1 3 500 r 1.0 10.0 table.txt LJ pair_write 1 3 500 r 1.0 10.0 table.txt LJ
pair_write 1 1 1000 rsq 2.0 8.0 table.txt Yukawa_1_1 -0.5 0.5 pair_write 1 1 1000 rsq 2.0 8.0 table.txt Yukawa_1_1 -0.5 0.5
labelmap atom 1 C 2 H
pair_write C H 500 r 1.0 10.0 table.txt LJ
Description Description
""""""""""" """""""""""
@ -42,7 +45,7 @@ compared against the entry in the file, if present, and pair_write
will refuse to add a table if the units are not the same. will refuse to add a table if the units are not the same.
The energy and force values are computed at distances from inner to The energy and force values are computed at distances from inner to
outer for 2 interacting atoms of type itype and jtype, using the outer for 2 interacting atoms of type *itype* and *jtype*, using the
appropriate :doc:`pair_coeff <pair_coeff>` coefficients. If the style appropriate :doc:`pair_coeff <pair_coeff>` coefficients. If the style
is *r*, then N distances are used, evenly spaced in r; if the style is is *r*, then N distances are used, evenly spaced in r; if the style is
*rsq*, N distances are used, evenly spaced in r\^2. *rsq*, N distances are used, evenly spaced in r\^2.

View File

@ -373,7 +373,7 @@ void PairE3B::allocateE3B()
void PairE3B::settings(int narg, char **arg) void PairE3B::settings(int narg, char **arg)
{ {
if (narg != 1) error->all(FLERR, "Illegal pair_style command"); if (narg != 1) error->all(FLERR, "Illegal pair_style command");
typeO = utils::inumeric(FLERR, arg[0], false, lmp); typeO_str = arg[0];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -387,6 +387,9 @@ void PairE3B::coeff(int narg, char **arg)
//1=* 2=* 3/4=1st keyword/value //1=* 2=* 3/4=1st keyword/value
if (narg < 4) error->all(FLERR, "There must be at least one keyword given to pair_coeff"); if (narg < 4) error->all(FLERR, "There must be at least one keyword given to pair_coeff");
if (typeO_str.size() > 0)
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
// clear setflag since coeff() called once with I,J = * * // clear setflag since coeff() called once with I,J = * *
int n = atom->ntypes; int n = atom->ntypes;
for (int i = 1; i <= n; i++) for (int i = 1; i <= n; i++)

View File

@ -36,6 +36,7 @@ class PairE3B : public Pair {
protected: protected:
//potential parameters //potential parameters
std::string typeO_str;
int typeO; int typeO;
double ea, eb, ec; //three body energies double ea, eb, ec; //three body energies
double k3; //three body exponential decay (units inverse length) double k3; //three body exponential decay (units inverse length)

View File

@ -408,18 +408,18 @@ void PairLJCutTIP4PLongSoft::settings(int narg, char **arg)
{ {
if (narg < 9 || narg > 10) error->all(FLERR,"Illegal pair_style command"); if (narg < 9 || narg > 10) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
nlambda = utils::numeric(FLERR,arg[5],false,lmp); nlambda = utils::numeric(FLERR, arg[5], false, lmp);
alphalj = utils::numeric(FLERR,arg[6],false,lmp); alphalj = utils::numeric(FLERR, arg[6], false, lmp);
alphac = utils::numeric(FLERR,arg[7],false,lmp); alphac = utils::numeric(FLERR, arg[7], false, lmp);
cut_lj_global = utils::numeric(FLERR,arg[8],false,lmp); cut_lj_global = utils::numeric(FLERR, arg[8], false, lmp);
if (narg == 9) cut_coul = cut_lj_global; if (narg == 9) cut_coul = cut_lj_global;
else cut_coul = utils::numeric(FLERR,arg[9],false,lmp); else cut_coul = utils::numeric(FLERR, arg[9], false, lmp);
// reset cutoffs that have been explicitly set // reset cutoffs that have been explicitly set
@ -431,6 +431,25 @@ void PairLJCutTIP4PLongSoft::settings(int narg, char **arg)
} }
} }
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCutTIP4PLongSoft::coeff(int narg, char **arg)
{
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
PairLJCutCoulLongSoft::coeff(narg, arg);
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
init specific to this pair style init specific to this pair style
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -30,6 +30,7 @@ class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft {
~PairLJCutTIP4PLongSoft() override; ~PairLJCutTIP4PLongSoft() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override; void init_style() override;
double init_one(int, int) override; double init_one(int, int) override;
void write_restart_settings(FILE *fp) override; void write_restart_settings(FILE *fp) override;
@ -38,6 +39,7 @@ class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft {
double memory_usage() override; double memory_usage() override;
protected: protected:
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -376,16 +376,35 @@ void PairTIP4PLongSoft::settings(int narg, char **arg)
{ {
if (narg != 8) error->all(FLERR,"Illegal pair_style command"); if (narg != 8) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
nlambda = utils::numeric(FLERR,arg[5],false,lmp); nlambda = utils::numeric(FLERR, arg[5], false, lmp);
alphac = utils::numeric(FLERR,arg[6],false,lmp); alphac = utils::numeric(FLERR, arg[6], false, lmp);
cut_coul = utils::numeric(FLERR,arg[7],false,lmp); cut_coul = utils::numeric(FLERR, arg[7], false, lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairTIP4PLongSoft::coeff(int narg, char **arg)
{
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
PairCoulLongSoft::coeff(narg, arg);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -30,6 +30,7 @@ class PairTIP4PLongSoft : public PairCoulLongSoft {
~PairTIP4PLongSoft() override; ~PairTIP4PLongSoft() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override; void init_style() override;
double init_one(int, int) override; double init_one(int, int) override;
void write_restart_settings(FILE *fp) override; void write_restart_settings(FILE *fp) override;
@ -38,6 +39,7 @@ class PairTIP4PLongSoft : public PairCoulLongSoft {
double memory_usage() override; double memory_usage() override;
protected: protected:
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -425,15 +425,15 @@ void PairLJCutTIP4PLong::settings(int narg, char **arg)
{ {
if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command"); if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
cut_lj_global = utils::numeric(FLERR,arg[5],false,lmp); cut_lj_global = utils::numeric(FLERR, arg[5], false, lmp);
if (narg == 6) cut_coul = cut_lj_global; if (narg == 6) cut_coul = cut_lj_global;
else cut_coul = utils::numeric(FLERR,arg[6],false,lmp); else cut_coul = utils::numeric(FLERR, arg[6], false, lmp);
// reset cutoffs that have been explicitly set // reset cutoffs that have been explicitly set
@ -445,6 +445,25 @@ void PairLJCutTIP4PLong::settings(int narg, char **arg)
} }
} }
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCutTIP4PLong::coeff(int narg, char **arg)
{
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
PairLJCutCoulLong::coeff(narg, arg);
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
init specific to this pair style init specific to this pair style
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -30,6 +30,7 @@ class PairLJCutTIP4PLong : public PairLJCutCoulLong {
~PairLJCutTIP4PLong() override; ~PairLJCutTIP4PLong() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override; void init_style() override;
double init_one(int, int) override; double init_one(int, int) override;
void write_restart_settings(FILE *fp) override; void write_restart_settings(FILE *fp) override;
@ -38,6 +39,7 @@ class PairLJCutTIP4PLong : public PairLJCutCoulLong {
double memory_usage() override; double memory_usage() override;
protected: protected:
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -1440,16 +1440,16 @@ void PairLJLongTIP4PLong::settings(int narg, char **arg)
if (!((ewald_order^ewald_off)&(1<<1))) if (!((ewald_order^ewald_off)&(1<<1)))
error->all(FLERR, error->all(FLERR,
"Coulomb cut not supported in pair_style lj/long/tip4p/long"); "Coulomb cut not supported in pair_style lj/long/tip4p/long");
typeO = utils::inumeric(FLERR,arg[1],false,lmp); typeO_str = arg[1];
typeH = utils::inumeric(FLERR,arg[2],false,lmp); typeH_str = arg[2];
typeB = utils::inumeric(FLERR,arg[3],false,lmp); typeB_str = arg[3];
typeA = utils::inumeric(FLERR,arg[4],false,lmp); typeA_str = arg[4];
qdist = utils::numeric(FLERR,arg[5],false,lmp); qdist = utils::numeric(FLERR, arg[5], false, lmp);
cut_lj_global = utils::numeric(FLERR,arg[6],false,lmp); cut_lj_global = utils::numeric(FLERR, arg[6], false, lmp);
if (narg == 8) cut_coul = cut_lj_global; if (narg == 8) cut_coul = cut_lj_global;
else cut_coul = utils::numeric(FLERR,arg[7],false,lmp); else cut_coul = utils::numeric(FLERR, arg[7], false, lmp);
// reset cutoffs that have been explicitly set // reset cutoffs that have been explicitly set
@ -1462,6 +1462,25 @@ void PairLJLongTIP4PLong::settings(int narg, char **arg)
} }
} }
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJLongTIP4PLong::coeff(int narg, char **arg)
{
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
PairLJLongCoulLong::coeff(narg, arg);
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
init specific to this pair style init specific to this pair style
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -33,6 +33,7 @@ class PairLJLongTIP4PLong : public PairLJLongCoulLong {
void compute_middle() override; void compute_middle() override;
void compute_outer(int, int) override; void compute_outer(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override; void init_style() override;
double init_one(int, int) override; double init_one(int, int) override;
void write_restart_settings(FILE *fp) override; void write_restart_settings(FILE *fp) override;
@ -41,6 +42,7 @@ class PairLJLongTIP4PLong : public PairLJLongCoulLong {
double memory_usage() override; double memory_usage() override;
protected: protected:
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -392,13 +392,32 @@ void PairTIP4PLong::settings(int narg, char **arg)
{ {
if (narg != 6) error->all(FLERR,"Illegal pair_style command"); if (narg != 6) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
cut_coul = utils::numeric(FLERR,arg[5],false,lmp); cut_coul = utils::numeric(FLERR, arg[5], false, lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairTIP4PLong::coeff(int narg, char **arg)
{
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
PairCoulLong::coeff(narg, arg);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -30,6 +30,7 @@ class PairTIP4PLong : public PairCoulLong {
~PairTIP4PLong() override; ~PairTIP4PLong() override;
void compute(int, int) override; void compute(int, int) override;
void settings(int, char **) override; void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override; void init_style() override;
double init_one(int, int) override; double init_one(int, int) override;
void write_restart_settings(FILE *fp) override; void write_restart_settings(FILE *fp) override;
@ -38,6 +39,7 @@ class PairTIP4PLong : public PairCoulLong {
double memory_usage() override; double memory_usage() override;
protected: protected:
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -226,11 +226,27 @@ void PairATM::coeff(int narg, char **arg)
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi,jlo,jhi,klo,khi; int ilo,ihi,jlo,jhi,klo,khi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
utils::bounds(FLERR,arg[2],1,atom->ntypes,klo,khi,error); utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM);
double nu_one = utils::numeric(FLERR,arg[3],false,lmp); // if explicit numeric value or type label used, reorder such at I < J < K
// note that, in the above case, I is already guaranteed to be less than J
if (jlo == jhi && klo == khi && klo < jlo) {
klo = jlo;
jlo = khi;
khi = klo;
jhi = jlo;
if (ilo == ihi && jlo == jhi && jlo < ilo) {
jlo = ilo;
ilo = jhi;
jhi = jlo;
ihi = ilo;
}
}
double nu_one = utils::numeric(FLERR, arg[3], false, lmp);
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {

View File

@ -348,14 +348,12 @@ void PairSRP::settings(int narg, char **arg)
if (atom->tag_enable == 0) if (atom->tag_enable == 0)
error->all(FLERR,"Pair_style srp requires atom IDs"); error->all(FLERR,"Pair_style srp requires atom IDs");
cut_global = utils::numeric(FLERR,arg[0],false,lmp); cut_global = utils::numeric(FLERR, arg[0], false, lmp);
// wildcard // wildcard
if (strcmp(arg[1],"*") == 0) { if (strcmp(arg[1],"*") == 0) {
btype = 0; btype = 0;
} else { } else {
btype = utils::inumeric(FLERR,arg[1],false,lmp); btype_str = arg[1];
if ((btype > atom->nbondtypes) || (btype <= 0))
error->all(FLERR,"Illegal pair_style command");
} }
// settings // settings
@ -383,20 +381,10 @@ void PairSRP::settings(int narg, char **arg)
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"bptype") == 0) { } else if (strcmp(arg[iarg],"bptype") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command"); if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command");
bptype = utils::inumeric(FLERR,arg[iarg+1],false,lmp); bptype_str = arg[iarg+1];
if ((bptype < 1) || (bptype > atom->ntypes))
error->all(FLERR,"Illegal bond particle type for srp");
iarg += 2; iarg += 2;
} else error->all(FLERR,"Illegal pair srp command"); } else error->all(FLERR,"Illegal pair srp command");
} }
// reset cutoffs if explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= bptype; i++)
for (j = i; j <= bptype; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -409,6 +397,24 @@ void PairSRP::coeff(int narg, char **arg)
error->all(FLERR,"PairSRP: Incorrect args for pair coeff"); error->all(FLERR,"PairSRP: Incorrect args for pair coeff");
if (!allocated) allocate(); if (!allocated) allocate();
if (btype_str.size() > 0)
btype = utils::expand_type_int(FLERR, btype_str, Atom::BOND, lmp);
if ((btype > atom->nbondtypes) || (btype <= 0))
error->all(FLERR,"Invalid bond type {} for pair style srp", btype);
if (bptype_str.size() > 0)
bptype = utils::expand_type_int(FLERR, bptype_str, Atom::ATOM, lmp);
if ((bptype < 1) || (bptype > atom->ntypes))
error->all(FLERR,"Invalid bond particle type {} for pair style srp", bptype);
// reset cutoffs if explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= bptype; i++)
for (j = i; j <= bptype; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
// set ij bond-bond cutoffs // set ij bond-bond cutoffs
int ilo, ihi, jlo, jhi; int ilo, ihi, jlo, jhi;
utils::bounds(FLERR,arg[0], 1, bptype, ilo, ihi, error); utils::bounds(FLERR,arg[0], 1, bptype, ilo, ihi, error);

View File

@ -51,6 +51,7 @@ class PairSRP : public Pair {
double **a0; double **a0;
double **srp; double **srp;
double cut_global; double cut_global;
std::string bptype_str, btype_str;
int bptype; int bptype;
int btype; int btype;
class Fix *f_srp; class Fix *f_srp;

View File

@ -319,35 +319,35 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg)
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi,jlo,jhi,klo,khi; int ilo,ihi,jlo,jhi,klo,khi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
utils::bounds(FLERR,arg[2],1,atom->ntypes,klo,khi,error); utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM);
int donor_flag; int donor_flag;
if (strcmp(arg[3],"i") == 0) donor_flag = 0; if (strcmp(arg[3],"i") == 0) donor_flag = 0;
else if (strcmp(arg[3],"j") == 0) donor_flag = 1; else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
else error->all(FLERR,"Incorrect args for pair coefficients"); else error->all(FLERR,"Incorrect args for pair coefficients");
double epsilon_one = utils::numeric(FLERR,arg[4],false,lmp); double epsilon_one = utils::numeric(FLERR, arg[4], false, lmp);
double sigma_one = utils::numeric(FLERR,arg[5],false,lmp); double sigma_one = utils::numeric(FLERR, arg[5], false, lmp);
int ap_one = ap_global; int ap_one = ap_global;
if (narg > 6) ap_one = utils::inumeric(FLERR,arg[6],false,lmp); if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, lmp);
double cut_inner_one = cut_inner_global; double cut_inner_one = cut_inner_global;
double cut_outer_one = cut_outer_global; double cut_outer_one = cut_outer_global;
if (narg > 8) { if (narg > 8) {
cut_inner_one = utils::numeric(FLERR,arg[7],false,lmp); cut_inner_one = utils::numeric(FLERR, arg[7], false, lmp);
cut_outer_one = utils::numeric(FLERR,arg[8],false,lmp); cut_outer_one = utils::numeric(FLERR, arg[8], false, lmp);
} }
if (cut_inner_one>cut_outer_one) if (cut_inner_one>cut_outer_one)
error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
double cut_angle_one = cut_angle_global; double cut_angle_one = cut_angle_global;
if (narg == 10) cut_angle_one = utils::numeric(FLERR,arg[9],false,lmp) * MY_PI/180.0; if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0;
// grow params array if necessary // grow params array if necessary
if (nparams == maxparam) { if (nparams == maxparam) {
maxparam += CHUNK; maxparam += CHUNK;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),
"pair:params"); "pair:params");
// make certain all addional allocated storage is initialized // make certain all addional allocated storage is initialized

View File

@ -243,37 +243,37 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg)
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi,jlo,jhi,klo,khi; int ilo,ihi,jlo,jhi,klo,khi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
utils::bounds(FLERR,arg[2],1,atom->ntypes,klo,khi,error); utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM);
int donor_flag; int donor_flag;
if (strcmp(arg[3],"i") == 0) donor_flag = 0; if (strcmp(arg[3],"i") == 0) donor_flag = 0;
else if (strcmp(arg[3],"j") == 0) donor_flag = 1; else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
else error->all(FLERR,"Incorrect args for pair coefficients"); else error->all(FLERR,"Incorrect args for pair coefficients");
double d0_one = utils::numeric(FLERR,arg[4],false,lmp); double d0_one = utils::numeric(FLERR, arg[4], false, lmp);
double alpha_one = utils::numeric(FLERR,arg[5],false,lmp); double alpha_one = utils::numeric(FLERR, arg[5], false, lmp);
double r0_one = utils::numeric(FLERR,arg[6],false,lmp); double r0_one = utils::numeric(FLERR, arg[6], false, lmp);
int ap_one = ap_global; int ap_one = ap_global;
if (narg > 7) ap_one = utils::inumeric(FLERR,arg[7],false,lmp); if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, lmp);
double cut_inner_one = cut_inner_global; double cut_inner_one = cut_inner_global;
double cut_outer_one = cut_outer_global; double cut_outer_one = cut_outer_global;
if (narg > 9) { if (narg > 9) {
cut_inner_one = utils::numeric(FLERR,arg[8],false,lmp); cut_inner_one = utils::numeric(FLERR, arg[8], false, lmp);
cut_outer_one = utils::numeric(FLERR,arg[9],false,lmp); cut_outer_one = utils::numeric(FLERR, arg[9], false, lmp);
} }
if (cut_inner_one>cut_outer_one) if (cut_inner_one>cut_outer_one)
error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
double cut_angle_one = cut_angle_global; double cut_angle_one = cut_angle_global;
if (narg > 10) cut_angle_one = utils::numeric(FLERR,arg[10],false,lmp) * MY_PI/180.0; if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0;
// grow params array if necessary // grow params array if necessary
if (nparams == maxparam) { if (nparams == maxparam) {
maxparam += CHUNK; maxparam += CHUNK;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),
"pair:params"); "pair:params");
// make certain all addional allocated storage is initialized // make certain all addional allocated storage is initialized

View File

@ -430,15 +430,15 @@ void PairLJCutTIP4PCut::settings(int narg, char **arg)
{ {
if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command"); if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
cut_lj_global = utils::numeric(FLERR,arg[5],false,lmp); cut_lj_global = utils::numeric(FLERR, arg[5], false, lmp);
if (narg == 6) cut_coul = cut_lj_global; if (narg == 6) cut_coul = cut_lj_global;
else cut_coul = utils::numeric(FLERR,arg[6],false,lmp); else cut_coul = utils::numeric(FLERR, arg[6], false, lmp);
cut_coulsq = cut_coul * cut_coul; cut_coulsq = cut_coul * cut_coul;
cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist); cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist);
@ -461,6 +461,16 @@ void PairLJCutTIP4PCut::coeff(int narg, char **arg)
error->all(FLERR,"Incorrect args for pair coefficients"); error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
int ilo,ihi,jlo,jhi; int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);

View File

@ -50,6 +50,7 @@ class PairLJCutTIP4PCut : public Pair {
double **epsilon, **sigma; double **epsilon, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset; double **lj1, **lj2, **lj3, **lj4, **offset;
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -375,12 +375,12 @@ void PairTIP4PCut::settings(int narg, char **arg)
{ {
if (narg != 6) error->all(FLERR,"Illegal pair_style command"); if (narg != 6) error->all(FLERR,"Illegal pair_style command");
typeO = utils::inumeric(FLERR,arg[0],false,lmp); typeO_str = arg[0];
typeH = utils::inumeric(FLERR,arg[1],false,lmp); typeH_str = arg[1];
typeB = utils::inumeric(FLERR,arg[2],false,lmp); typeB_str = arg[2];
typeA = utils::inumeric(FLERR,arg[3],false,lmp); typeA_str = arg[3];
qdist = utils::numeric(FLERR,arg[4],false,lmp); qdist = utils::numeric(FLERR, arg[4], false, lmp);
cut_coul = utils::numeric(FLERR,arg[5],false,lmp); cut_coul = utils::numeric(FLERR, arg[5], false, lmp);
cut_coulsq = cut_coul * cut_coul; cut_coulsq = cut_coul * cut_coul;
cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist); cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist);
@ -396,6 +396,16 @@ void PairTIP4PCut::coeff(int narg, char **arg)
error->all(FLERR,"Incorrect args for pair coefficients"); error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
// set atom types from pair_style command unless we were restarted
// and the types are already set and the strings are empty.
if (typeO_str.size() > 0) {
typeO = utils::expand_type_int(FLERR, typeO_str, Atom::ATOM, lmp, true);
typeH = utils::expand_type_int(FLERR, typeH_str, Atom::ATOM, lmp, true);
typeB = utils::expand_type_int(FLERR, typeB_str, Atom::BOND, lmp, true);
typeA = utils::expand_type_int(FLERR, typeA_str, Atom::ANGLE, lmp, true);
}
int ilo,ihi,jlo,jhi; int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);

View File

@ -45,6 +45,7 @@ class PairTIP4PCut : public Pair {
double cut_coul, cut_coulsq; double cut_coul, cut_coulsq;
double cut_coulsqplus; // extended value for cut_coulsq double cut_coulsqplus; // extended value for cut_coulsq
std::string typeH_str, typeO_str, typeA_str, typeB_str;
int typeH, typeO; // atom types of TIP4P water H and O atoms int typeH, typeO; // atom types of TIP4P water H and O atoms
int typeA, typeB; // angle and bond types of TIP4P water int typeA, typeB; // angle and bond types of TIP4P water
double alpha; // geometric constraint parameter for TIP4P double alpha; // geometric constraint parameter for TIP4P

View File

@ -1795,8 +1795,8 @@ void Pair::write_file(int narg, char **arg)
// parse arguments // parse arguments
int itype = utils::inumeric(FLERR,arg[0],false,lmp); int itype = utils::expand_type_int(FLERR, arg[0], Atom::ATOM, lmp);
int jtype = utils::inumeric(FLERR,arg[1],false,lmp); int jtype = utils::expand_type_int(FLERR, arg[1], Atom::ATOM, lmp);
if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes) if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes)
error->all(FLERR,"Invalid atom types in pair_write command"); error->all(FLERR,"Invalid atom types in pair_write command");
@ -1810,8 +1810,8 @@ void Pair::write_file(int narg, char **arg)
if (n < 2) error->all(FLERR, "Must have at least 2 table values"); if (n < 2) error->all(FLERR, "Must have at least 2 table values");
double inner = utils::numeric(FLERR,arg[4],false,lmp); double inner = utils::numeric(FLERR, arg[4], false, lmp);
double outer = utils::numeric(FLERR,arg[5],false,lmp); double outer = utils::numeric(FLERR, arg[5], false, lmp);
if (inner <= 0.0 || inner >= outer) if (inner <= 0.0 || inner >= outer)
error->all(FLERR,"Invalid cutoffs in pair_write command"); error->all(FLERR,"Invalid cutoffs in pair_write command");
@ -1829,13 +1829,13 @@ void Pair::write_file(int narg, char **arg)
// - if the file already exists, print a message about appending // - if the file already exists, print a message about appending
// while printing the date and check that units are consistent. // while printing the date and check that units are consistent.
if (platform::file_is_readable(table_file)) { if (platform::file_is_readable(table_file)) {
std::string units = utils::get_potential_units(table_file,"table"); std::string units = utils::get_potential_units(table_file, "table");
if (!units.empty() && (units != update->unit_style)) { if (!units.empty() && (units != update->unit_style)) {
error->one(FLERR,"Trying to append to a table file " error->one(FLERR,"Trying to append to a table file "
"with UNITS: {} while units are {}", "with UNITS: {} while units are {}",
units, update->unit_style); units, update->unit_style);
} }
std::string date = utils::get_potential_date(table_file,"table"); std::string date = utils::get_potential_date(table_file, "table");
utils::logmesg(lmp,"Appending to table file {} with DATE: {}\n", table_file, date); utils::logmesg(lmp,"Appending to table file {} with DATE: {}\n", table_file, date);
fp = fopen(table_file.c_str(),"a"); fp = fopen(table_file.c_str(),"a");
} else { } else {
@ -1847,12 +1847,12 @@ void Pair::write_file(int narg, char **arg)
} }
if (fp == nullptr) if (fp == nullptr)
error->one(FLERR,"Cannot open pair_write file {}: {}",table_file, utils::getsyserror()); error->one(FLERR,"Cannot open pair_write file {}: {}",table_file, utils::getsyserror());
fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n", fprintf(fp, "# Pair potential %s for atom types %d %d: i,r,energy,force\n",
force->pair_style,itype,jtype); force->pair_style, itype, jtype);
if (style == RLINEAR) if (style == RLINEAR)
fprintf(fp,"\n%s\nN %d R %.15g %.15g\n\n",arg[7],n,inner,outer); fprintf(fp, "\n%s\nN %d R %.15g %.15g\n\n", arg[7], n, inner, outer);
if (style == RSQ) if (style == RSQ)
fprintf(fp,"\n%s\nN %d RSQ %.15g %.15g\n\n",arg[7],n,inner,outer); fprintf(fp, "\n%s\nN %d RSQ %.15g %.15g\n\n", arg[7], n, inner, outer);
} }
// initialize potentials before evaluating pair potential // initialize potentials before evaluating pair potential
@ -1870,7 +1870,7 @@ void Pair::write_file(int narg, char **arg)
double *eamfp_hold; double *eamfp_hold;
Pair *epair = force->pair_match("^eam",0); Pair *epair = force->pair_match("^eam",0);
if (epair) epair->swap_eam(eamfp,&eamfp_hold); if (epair) epair->swap_eam(eamfp, &eamfp_hold);
if ((comm->me == 0) && (epair)) if ((comm->me == 0) && (epair))
error->warning(FLERR,"EAM pair style. Table will not include embedding term"); error->warning(FLERR,"EAM pair style. Table will not include embedding term");
@ -1879,8 +1879,8 @@ void Pair::write_file(int narg, char **arg)
double q[2]; double q[2];
q[0] = q[1] = 1.0; q[0] = q[1] = 1.0;
if (narg == 10) { if (narg == 10) {
q[0] = utils::numeric(FLERR,arg[8],false,lmp); q[0] = utils::numeric(FLERR, arg[8], false, lmp);
q[1] = utils::numeric(FLERR,arg[9],false,lmp); q[1] = utils::numeric(FLERR, arg[9], false, lmp);
} }
double *q_hold; double *q_hold;
@ -1893,10 +1893,10 @@ void Pair::write_file(int narg, char **arg)
int masklo,maskhi,nmask,nshiftbits; int masklo,maskhi,nmask,nshiftbits;
if (style == BMP) { if (style == BMP) {
init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits); init_bitmap(inner, outer, n, masklo, maskhi, nmask, nshiftbits);
int ntable = 1 << n; int ntable = 1 << n;
if (comm->me == 0) if (comm->me == 0)
fprintf(fp,"\n%s\nN %d BITMAP %.15g %.15g\n\n",arg[7],ntable,inner,outer); fprintf(fp, "\n%s\nN %d BITMAP %.15g %.15g\n\n", arg[7], ntable, inner, outer);
n = ntable; n = ntable;
} }
@ -1922,7 +1922,7 @@ void Pair::write_file(int narg, char **arg)
} }
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
e = single(0,1,itype,jtype,rsq,1.0,1.0,f); e = single(0, 1, itype, jtype, rsq, 1.0, 1.0, f);
f *= r; f *= r;
} else e = f = 0.0; } else e = f = 0.0;
if (comm->me == 0) fprintf(fp,"%8d %- 22.15g %- 22.15g %- 22.15g\n",i+1,r,e,f); if (comm->me == 0) fprintf(fp,"%8d %- 22.15g %- 22.15g %- 22.15g\n",i+1,r,e,f);
@ -1931,7 +1931,7 @@ void Pair::write_file(int narg, char **arg)
// restore original vecs that were swapped in for // restore original vecs that were swapped in for
double *tmp; double *tmp;
if (epair) epair->swap_eam(eamfp_hold,&tmp); if (epair) epair->swap_eam(eamfp_hold, &tmp);
if (atom->q) atom->q = q_hold; if (atom->q) atom->q = q_hold;
if (comm->me == 0) fclose(fp); if (comm->me == 0) fclose(fp);
@ -2010,4 +2010,3 @@ double Pair::memory_usage()
bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double); bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double);
return bytes; return bytes;
} }

View File

@ -995,17 +995,42 @@ char *utils::expand_type(const char *file, int line, const std::string &str, int
/* ------------------------------------------------------------------------- /* -------------------------------------------------------------------------
Expand type string to integer-valued numeric type from labelmap. Expand type string to integer-valued numeric type from labelmap.
Not guaranteed to return a valid type. Not guaranteed to return a valid type if param verify = 0 (default)
For example, type <= 0 or type > Ntypes is checked in calling routine. In this case, type <= 0 or type > Ntypes should be checked in calling routine.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
int utils::expand_type_int(const char *file, int line, const std::string &str, int mode, int utils::expand_type_int(const char *file, int line, const std::string &str, int mode,
LAMMPS *lmp) LAMMPS *lmp, bool verify)
{ {
char *typestr = expand_type(file, line, str, mode, lmp); char *typestr = expand_type(file, line, str, mode, lmp);
int out = inumeric(file, line, typestr ? typestr : str, false, lmp); int type = inumeric(file, line, typestr ? typestr : str, false, lmp);
if (verify) {
int errorflag = 0;
if (type <= 0) errorflag = 1;
int nmax;
switch (mode) {
case Atom::ATOM:
nmax = lmp->atom->ntypes;
break;
case Atom::BOND:
nmax = lmp->atom->nbondtypes;
break;
case Atom::ANGLE:
nmax = lmp->atom->nangletypes;
break;
case Atom::DIHEDRAL:
nmax = lmp->atom->ndihedraltypes;
break;
case Atom::IMPROPER:
nmax = lmp->atom->nimpropertypes;
break;
}
if (type > nmax) errorflag = 1;
if (errorflag) lmp->error->all(file, line, "{} type {} is out of bounds ({}-{})",
labeltypes[mode], type, 1, nmax);
}
delete[] typestr; delete[] typestr;
return out; return type;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -409,7 +409,7 @@ This functions adds the following case to :cpp:func:`utils::bounds() <LAMMPS_NS:
* *
* This function has the same arguments as expand_type() but returns an integer value */ * This function has the same arguments as expand_type() but returns an integer value */
int expand_type_int(const char *file, int line, const std::string &str, int mode, LAMMPS *lmp); int expand_type_int(const char *file, int line, const std::string &str, int mode, LAMMPS *lmp, bool verify = false);
/*! Check grid reference for valid Compute or Fix which produces per-grid data /*! Check grid reference for valid Compute or Fix which produces per-grid data
* *
@ -425,6 +425,7 @@ This functions adds the following case to :cpp:func:`utils::bounds() <LAMMPS_NS:
* \param nevery frequency at which caller will access fix for per-grid info, * \param nevery frequency at which caller will access fix for per-grid info,
* ignored when reference is to a compute * ignored when reference is to a compute
* \param lmp pointer to top-level LAMMPS class instance * \param lmp pointer to top-level LAMMPS class instance
* \param verify check bounds for interaction type
* \return id ID of Compute or Fix * \return id ID of Compute or Fix
* \return igrid which grid is referenced (0 to N-1) * \return igrid which grid is referenced (0 to N-1)
* \return idata which data on grid is referenced (0 to N-1) * \return idata which data on grid is referenced (0 to N-1)