diff --git a/doc/src/Howto_tip4p.rst b/doc/src/Howto_tip4p.rst index d1400ec479..47a1b9b578 100644 --- a/doc/src/Howto_tip4p.rst +++ b/doc/src/Howto_tip4p.rst @@ -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: - * :doc:`pair_style tip4p/cut ` + * :doc:`pair_style tip4p/cut ` * :doc:`pair_style lj/cut/tip4p/cut ` or these commands for a long-range Coulomb treatment: diff --git a/doc/src/pair_atm.rst b/doc/src/pair_atm.rst index 3e89bba243..c0b06504a0 100644 --- a/doc/src/pair_atm.rst +++ b/doc/src/pair_atm.rst @@ -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:`\nu` = prefactor (energy/distance\^9 units) -:math:`K` can be specified in one of two ways. An explicit numeric value can -be used, as in the second example above. :math:`J \leq K` is required. LAMMPS +:math:`K` can be specified in one of two ways. An explicit numeric value +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 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}`, diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 3043aa34ac..e8f09515b0 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -93,12 +93,19 @@ Syntax pair_style coul/long cutoff pair_style coul/wolf alpha cutoff 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/long otype htype btype atype qdist cutoff -* cutoff = global cutoff for Coulombic interactions -* kappa = Debye length (inverse distance units) -* alpha = damping parameter (inverse distance units) + * otype,htype = atom types (numeric or type label) for TIP4P O and H + * btype,atype = bond and angle types (numeric or type label) for TIP4P waters + * qdist = distance from O atom to massless charge (distance units) Examples """""""" @@ -138,6 +145,12 @@ Examples pair_style tip4p/long 1 2 7 8 0.15 10.0 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 """"""""""" @@ -307,6 +320,11 @@ Coulombic solver (Ewald or PPPM). 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. +.. note:: + + If using type labels, the type labels must be defined before calling + the :doc:`pair_coeff ` command. + See the :doc:`Howto tip4p ` page for more information on how to use the TIP4P pair styles and lists of parameters to set. 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 -part of the KSPACE package. The *coul/cut/global* and *coul/exclude* are -part of the EXTRA-PAIR package. A pair style is only enabled if LAMMPS was -built with its corresponding package. See the :doc:`Build package ` +part of the KSPACE package. The *coul/cut/global*, *coul/exclude* styles are +part of the EXTRA-PAIR package. The *tip4p/cut* style is part of the MOLECULE +package. A pair style is only enabled if LAMMPS was built with its +corresponding package. See the :doc:`Build package ` doc page for more info. Related commands diff --git a/doc/src/pair_e3b.rst b/doc/src/pair_e3b.rst index 8ae6d10b82..f88a5331b9 100644 --- a/doc/src/pair_e3b.rst +++ b/doc/src/pair_e3b.rst @@ -10,7 +10,7 @@ Syntax pair_style e3b Otype -* Otype = atom type for oxygen +* Otype = atom type (numeric or type label) for oxygen .. 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_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: .. 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. 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 ` command. + The *pair_coeff* command must have at least one keyword/value pair, as described above. The *preset* keyword sets the potential parameters to the values used in :ref:`(Tainter 2011) ` or diff --git a/doc/src/pair_fep_soft.rst b/doc/src/pair_fep_soft.rst index af8d0668c0..3a9b6801b2 100644 --- a/doc/src/pair_fep_soft.rst +++ b/doc/src/pair_fep_soft.rst @@ -97,8 +97,8 @@ Syntax cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (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) - otype,htype = atom types for TIP4P O and H - btype,atype = bond and angle types for TIP4P waters + otype,htype = atom types (numeric or type label) for TIP4P O and H + btype,atype = bond and angle types (numeric or type label) for TIP4P waters qdist = distance from O atom to massless charge (distance units) n, alpha_LJ, alpha_C = parameters of the soft-core potential 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 cutoff = global cutoff for Coulomb interactions (distance units) *tip4p/long/soft* args = otype htype btype atype qdist n alpha_C cutoff - otype,htype = atom types for TIP4P O and H - btype,atype = bond and angle types for TIP4P waters + otype,htype = atom types (numeric or type label) for TIP4P O and H + btype,atype = bond and angle types (numeric or type label) for TIP4P waters qdist = distance from O atom to massless charge (distance units) n, alpha_C = parameters of the soft-core potential cutoff = global cutoff for Coulomb interactions (distance units) @@ -161,6 +161,13 @@ Examples pair_coeff * * 0.155 3.1536 1.0 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 9.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 ` command, after 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 ` command. + 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 lj/charmm/coul/long ` style. In the soft version the parameters diff --git a/doc/src/pair_hbond_dreiding.rst b/doc/src/pair_hbond_dreiding.rst index f9b2a010fd..7e73f23b08 100644 --- a/doc/src/pair_hbond_dreiding.rst +++ b/doc/src/pair_hbond_dreiding.rst @@ -38,6 +38,9 @@ Examples 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 + 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 """"""""""" @@ -144,7 +147,7 @@ in the examples above. For the *hbond/dreiding/lj* style the list of coefficients is as follows: -* K = hydrogen atom type = 1 to Ntypes +* K = hydrogen atom type = 1 to Ntypes, or type label * donor flag = *i* or *j* * :math:`\epsilon` (energy units) * :math:`\sigma` (distance units) @@ -156,7 +159,7 @@ follows: For the *hbond/dreiding/morse* style the list of coefficients is as follows: -* K = hydrogen atom type = 1 to Ntypes +* K = hydrogen atom type = 1 to Ntypes, or type label * donor flag = *i* or *j* * :math:`D_0` (energy units) * :math:`\alpha` (1/distance units) @@ -240,7 +243,10 @@ heading) the following commands could be included in an input script: Restrictions """""""""""" - none + +This pair style can only be used if LAMMPS was built with the +MOLECULE package. See the :doc:`Build package ` doc page +for more info. Related commands """""""""""""""" diff --git a/doc/src/pair_lj_cut_tip4p.rst b/doc/src/pair_lj_cut_tip4p.rst index f30e474f9a..7889e12bb5 100644 --- a/doc/src/pair_lj_cut_tip4p.rst +++ b/doc/src/pair_lj_cut_tip4p.rst @@ -28,14 +28,14 @@ Syntax .. parsed-literal:: *lj/cut/tip4p/cut* args = otype htype btype atype qdist cutoff (cutoff2) - otype,htype = atom types for TIP4P O and H - btype,atype = bond and angle types for TIP4P waters + otype,htype = atom types (numeric or type label) for TIP4P O and H + btype,atype = bond and angle types (numeric or type label) for TIP4P waters qdist = distance from O atom to massless charge (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) *lj/cut/tip4p/long* args = otype htype btype atype qdist cutoff (cutoff2) - otype,htype = atom types for TIP4P O and H - btype,atype = bond and angle types for TIP4P waters + otype,htype = atom types (numeric or type label) for TIP4P O and H + btype,atype = bond and angle types (numeric or type label) for TIP4P waters qdist = distance from O atom to massless charge (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) @@ -55,6 +55,13 @@ Examples pair_coeff * * 100.0 3.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 """"""""""" @@ -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 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 ` command. + See the :doc:`Howto tip4p ` page for more information on how to use the TIP4P pair styles and lists of parameters to set. Note that the neighbor list cutoff for Coulomb interactions is diff --git a/doc/src/pair_lj_long.rst b/doc/src/pair_lj_long.rst index 050d1c0602..ba824fcfd9 100644 --- a/doc/src/pair_lj_long.rst +++ b/doc/src/pair_lj_long.rst @@ -44,8 +44,8 @@ Syntax flag_coul = *long* or *off* *long* = use Kspace long-range summation for Coulombic 1/r term *off* = omit Coulombic term - otype,htype = atom types for TIP4P O and H - btype,atype = bond and angle types for TIP4P waters + otype,htype = atom types (numeric or type label) for TIP4P O and H + btype,atype = bond and angle types (numeric or type label) for TIP4P waters qdist = distance from O atom to massless charge (distance units) cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) @@ -66,6 +66,13 @@ Examples pair_coeff * * 100.0 3.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 """"""""""" @@ -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 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 ` command. + See the :doc:`Howto tip4p ` page for more information on how to use the TIP4P pair style. Note that the neighbor list cutoff for Coulomb interactions is effectively extended diff --git a/doc/src/pair_srp.rst b/doc/src/pair_srp.rst index 7122c40c1e..b1c3c140d3 100644 --- a/doc/src/pair_srp.rst +++ b/doc/src/pair_srp.rst @@ -15,7 +15,7 @@ Syntax pair_style srp/react cutoff btype dist react-id keyword value ... * 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* * react-id = id of either fix bond/break or fix bond/create * zero or more keyword/value pairs may be appended @@ -23,7 +23,7 @@ Syntax .. 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* Examples @@ -52,6 +52,10 @@ Examples pair_coeff 1 2 none pair_coeff 2 2 srp 100.0 0.8 + labelmap bond 1 C-C + pair_style hybrid srp 0.8 C-C mid + + Description @@ -117,6 +121,11 @@ is used. between atom type *bptype* and all other types of atoms. An error will be flagged if :doc:`pair_style hybrid ` is not used. +.. note:: + + If using type labels, the type labels must be defined before calling + the :doc:`pair_coeff ` command. + The optional *exclude* keyword determines if forces are computed between first neighbor (directly connected) bonds. For a setting of *no*, first neighbor forces are computed; for *yes* they are not @@ -207,4 +216,3 @@ Chem Phys, 136 (13) 134903, 2012. .. _Palkar: **(Palkar)** Palkar V, Kuksenok O, J. Phys. Chem. B, 126 (1), 336-346, 2022 - diff --git a/doc/src/pair_write.rst b/doc/src/pair_write.rst index b62383cb7b..151d3dee5e 100644 --- a/doc/src/pair_write.rst +++ b/doc/src/pair_write.rst @@ -10,7 +10,7 @@ Syntax 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 * style = *r* or *rsq* or *bitmap* * 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 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 """"""""""" @@ -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. 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 ` coefficients. If the style 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. diff --git a/src/EXTRA-PAIR/pair_e3b.cpp b/src/EXTRA-PAIR/pair_e3b.cpp index efd8c5e9af..cc81cee03f 100644 --- a/src/EXTRA-PAIR/pair_e3b.cpp +++ b/src/EXTRA-PAIR/pair_e3b.cpp @@ -373,7 +373,7 @@ void PairE3B::allocateE3B() void PairE3B::settings(int narg, char **arg) { 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 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 = * * int n = atom->ntypes; for (int i = 1; i <= n; i++) diff --git a/src/EXTRA-PAIR/pair_e3b.h b/src/EXTRA-PAIR/pair_e3b.h index c61bb10798..e5a3293dbb 100644 --- a/src/EXTRA-PAIR/pair_e3b.h +++ b/src/EXTRA-PAIR/pair_e3b.h @@ -36,6 +36,7 @@ class PairE3B : public Pair { protected: //potential parameters + std::string typeO_str; int typeO; double ea, eb, ec; //three body energies double k3; //three body exponential decay (units inverse length) diff --git a/src/FEP/pair_lj_cut_tip4p_long_soft.cpp b/src/FEP/pair_lj_cut_tip4p_long_soft.cpp index d0a0d846a9..3cd444d685 100644 --- a/src/FEP/pair_lj_cut_tip4p_long_soft.cpp +++ b/src/FEP/pair_lj_cut_tip4p_long_soft.cpp @@ -408,18 +408,18 @@ void PairLJCutTIP4PLongSoft::settings(int narg, char **arg) { if (narg < 9 || narg > 10) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); - nlambda = utils::numeric(FLERR,arg[5],false,lmp); - alphalj = utils::numeric(FLERR,arg[6],false,lmp); - alphac = utils::numeric(FLERR,arg[7],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + qdist = utils::numeric(FLERR, arg[4], false, lmp); + nlambda = utils::numeric(FLERR, arg[5], false, lmp); + alphalj = utils::numeric(FLERR, arg[6], 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; - 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 @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/FEP/pair_lj_cut_tip4p_long_soft.h b/src/FEP/pair_lj_cut_tip4p_long_soft.h index 3d54376a20..7bbf2e36d3 100644 --- a/src/FEP/pair_lj_cut_tip4p_long_soft.h +++ b/src/FEP/pair_lj_cut_tip4p_long_soft.h @@ -30,6 +30,7 @@ class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft { ~PairLJCutTIP4PLongSoft() override; void compute(int, int) override; void settings(int, char **) override; + void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart_settings(FILE *fp) override; @@ -38,6 +39,7 @@ class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft { double memory_usage() override; protected: + std::string typeH_str, typeO_str, typeA_str, typeB_str; int typeH, typeO; // atom types of TIP4P water H and O atoms int typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/FEP/pair_tip4p_long_soft.cpp b/src/FEP/pair_tip4p_long_soft.cpp index 09351f9e05..525112fc89 100644 --- a/src/FEP/pair_tip4p_long_soft.cpp +++ b/src/FEP/pair_tip4p_long_soft.cpp @@ -376,16 +376,35 @@ void PairTIP4PLongSoft::settings(int narg, char **arg) { if (narg != 8) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + qdist = utils::numeric(FLERR, arg[4], false, lmp); - nlambda = utils::numeric(FLERR,arg[5],false,lmp); - alphac = utils::numeric(FLERR,arg[6],false,lmp); + nlambda = utils::numeric(FLERR, arg[5], 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); } /* ---------------------------------------------------------------------- diff --git a/src/FEP/pair_tip4p_long_soft.h b/src/FEP/pair_tip4p_long_soft.h index af30d77b05..b503f610fb 100644 --- a/src/FEP/pair_tip4p_long_soft.h +++ b/src/FEP/pair_tip4p_long_soft.h @@ -30,6 +30,7 @@ class PairTIP4PLongSoft : public PairCoulLongSoft { ~PairTIP4PLongSoft() override; void compute(int, int) override; void settings(int, char **) override; + void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart_settings(FILE *fp) override; @@ -38,6 +39,7 @@ class PairTIP4PLongSoft : public PairCoulLongSoft { double memory_usage() override; protected: + std::string typeH_str, typeO_str, typeA_str, typeB_str; int typeH, typeO; // atom types of TIP4P water H and O atoms int typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/KSPACE/pair_lj_cut_tip4p_long.cpp b/src/KSPACE/pair_lj_cut_tip4p_long.cpp index 32a04e2761..2420ab9c3d 100644 --- a/src/KSPACE/pair_lj_cut_tip4p_long.cpp +++ b/src/KSPACE/pair_lj_cut_tip4p_long.cpp @@ -425,15 +425,15 @@ void PairLJCutTIP4PLong::settings(int narg, char **arg) { if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + 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; - 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 @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_cut_tip4p_long.h b/src/KSPACE/pair_lj_cut_tip4p_long.h index 3a07ce31d3..4624375b64 100644 --- a/src/KSPACE/pair_lj_cut_tip4p_long.h +++ b/src/KSPACE/pair_lj_cut_tip4p_long.h @@ -30,6 +30,7 @@ class PairLJCutTIP4PLong : public PairLJCutCoulLong { ~PairLJCutTIP4PLong() override; void compute(int, int) override; void settings(int, char **) override; + void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart_settings(FILE *fp) override; @@ -38,6 +39,7 @@ class PairLJCutTIP4PLong : public PairLJCutCoulLong { double memory_usage() override; protected: + std::string typeH_str, typeO_str, typeA_str, typeB_str; int typeH, typeO; // atom types of TIP4P water H and O atoms int typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/KSPACE/pair_lj_long_tip4p_long.cpp b/src/KSPACE/pair_lj_long_tip4p_long.cpp index 187b22a78a..de561156a3 100644 --- a/src/KSPACE/pair_lj_long_tip4p_long.cpp +++ b/src/KSPACE/pair_lj_long_tip4p_long.cpp @@ -1440,16 +1440,16 @@ void PairLJLongTIP4PLong::settings(int narg, char **arg) if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR, "Coulomb cut not supported in pair_style lj/long/tip4p/long"); - typeO = utils::inumeric(FLERR,arg[1],false,lmp); - typeH = utils::inumeric(FLERR,arg[2],false,lmp); - typeB = utils::inumeric(FLERR,arg[3],false,lmp); - typeA = utils::inumeric(FLERR,arg[4],false,lmp); - qdist = utils::numeric(FLERR,arg[5],false,lmp); + typeO_str = arg[1]; + typeH_str = arg[2]; + typeB_str = arg[3]; + typeA_str = arg[4]; + 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; - 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 @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_long_tip4p_long.h b/src/KSPACE/pair_lj_long_tip4p_long.h index ba87117302..549580bc62 100644 --- a/src/KSPACE/pair_lj_long_tip4p_long.h +++ b/src/KSPACE/pair_lj_long_tip4p_long.h @@ -33,6 +33,7 @@ class PairLJLongTIP4PLong : public PairLJLongCoulLong { void compute_middle() override; void compute_outer(int, int) override; void settings(int, char **) override; + void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart_settings(FILE *fp) override; @@ -41,6 +42,7 @@ class PairLJLongTIP4PLong : public PairLJLongCoulLong { double memory_usage() override; protected: + std::string typeH_str, typeO_str, typeA_str, typeB_str; int typeH, typeO; // atom types of TIP4P water H and O atoms int typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/KSPACE/pair_tip4p_long.cpp b/src/KSPACE/pair_tip4p_long.cpp index 637a272e49..0608968fdf 100644 --- a/src/KSPACE/pair_tip4p_long.cpp +++ b/src/KSPACE/pair_tip4p_long.cpp @@ -392,13 +392,32 @@ void PairTIP4PLong::settings(int narg, char **arg) { if (narg != 6) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + 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); } /* ---------------------------------------------------------------------- diff --git a/src/KSPACE/pair_tip4p_long.h b/src/KSPACE/pair_tip4p_long.h index a0f2915743..890c70c8a6 100644 --- a/src/KSPACE/pair_tip4p_long.h +++ b/src/KSPACE/pair_tip4p_long.h @@ -30,6 +30,7 @@ class PairTIP4PLong : public PairCoulLong { ~PairTIP4PLong() override; void compute(int, int) override; void settings(int, char **) override; + void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart_settings(FILE *fp) override; @@ -38,6 +39,7 @@ class PairTIP4PLong : public PairCoulLong { double memory_usage() override; protected: + std::string typeH_str, typeO_str, typeA_str, typeB_str; int typeH, typeO; // atom types of TIP4P water H and O atoms int typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/MANYBODY/pair_atm.cpp b/src/MANYBODY/pair_atm.cpp index a37da6c10a..671b77b206 100644 --- a/src/MANYBODY/pair_atm.cpp +++ b/src/MANYBODY/pair_atm.cpp @@ -226,11 +226,27 @@ void PairATM::coeff(int narg, char **arg) if (!allocated) allocate(); int ilo,ihi,jlo,jhi,klo,khi; - 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[2],1,atom->ntypes,klo,khi,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_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; for (int i = ilo; i <= ihi; i++) { diff --git a/src/MISC/pair_srp.cpp b/src/MISC/pair_srp.cpp index 31f5b85760..8decae5aca 100644 --- a/src/MISC/pair_srp.cpp +++ b/src/MISC/pair_srp.cpp @@ -348,14 +348,12 @@ void PairSRP::settings(int narg, char **arg) if (atom->tag_enable == 0) 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 if (strcmp(arg[1],"*") == 0) { btype = 0; } else { - btype = utils::inumeric(FLERR,arg[1],false,lmp); - if ((btype > atom->nbondtypes) || (btype <= 0)) - error->all(FLERR,"Illegal pair_style command"); + btype_str = arg[1]; } // settings @@ -383,20 +381,10 @@ void PairSRP::settings(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"bptype") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command"); - bptype = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if ((bptype < 1) || (bptype > atom->ntypes)) - error->all(FLERR,"Illegal bond particle type for srp"); + bptype_str = arg[iarg+1]; iarg += 2; } 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"); 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 int ilo, ihi, jlo, jhi; utils::bounds(FLERR,arg[0], 1, bptype, ilo, ihi, error); diff --git a/src/MISC/pair_srp.h b/src/MISC/pair_srp.h index 3fb1563728..159e3fbbfe 100644 --- a/src/MISC/pair_srp.h +++ b/src/MISC/pair_srp.h @@ -51,6 +51,7 @@ class PairSRP : public Pair { double **a0; double **srp; double cut_global; + std::string bptype_str, btype_str; int bptype; int btype; class Fix *f_srp; diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index dbd7db7780..274f8bc2a3 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -319,35 +319,35 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) if (!allocated) allocate(); int ilo,ihi,jlo,jhi,klo,khi; - 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[2],1,atom->ntypes,klo,khi,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_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); int donor_flag; if (strcmp(arg[3],"i") == 0) donor_flag = 0; else if (strcmp(arg[3],"j") == 0) donor_flag = 1; else error->all(FLERR,"Incorrect args for pair coefficients"); - double epsilon_one = utils::numeric(FLERR,arg[4],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[5],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[4], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[5], false, lmp); 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_outer_one = cut_outer_global; if (narg > 8) { - cut_inner_one = utils::numeric(FLERR,arg[7],false,lmp); - cut_outer_one = utils::numeric(FLERR,arg[8],false,lmp); + cut_inner_one = utils::numeric(FLERR, arg[7], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[8], false, lmp); } if (cut_inner_one>cut_outer_one) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); 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 if (nparams == maxparam) { maxparam += CHUNK; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), "pair:params"); // make certain all addional allocated storage is initialized diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 5cc45ea234..c8bc0a627d 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -243,37 +243,37 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) if (!allocated) allocate(); int ilo,ihi,jlo,jhi,klo,khi; - 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[2],1,atom->ntypes,klo,khi,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_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); int donor_flag; if (strcmp(arg[3],"i") == 0) donor_flag = 0; else if (strcmp(arg[3],"j") == 0) donor_flag = 1; else error->all(FLERR,"Incorrect args for pair coefficients"); - double d0_one = utils::numeric(FLERR,arg[4],false,lmp); - double alpha_one = utils::numeric(FLERR,arg[5],false,lmp); - double r0_one = utils::numeric(FLERR,arg[6],false,lmp); + double d0_one = utils::numeric(FLERR, arg[4], false, lmp); + double alpha_one = utils::numeric(FLERR, arg[5], false, lmp); + double r0_one = utils::numeric(FLERR, arg[6], false, lmp); 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_outer_one = cut_outer_global; if (narg > 9) { - cut_inner_one = utils::numeric(FLERR,arg[8],false,lmp); - cut_outer_one = utils::numeric(FLERR,arg[9],false,lmp); + cut_inner_one = utils::numeric(FLERR, arg[8], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[9], false, lmp); } if (cut_inner_one>cut_outer_one) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); 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 if (nparams == maxparam) { maxparam += CHUNK; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), "pair:params"); // make certain all addional allocated storage is initialized diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp index 124e7e4ed8..b8946b67f4 100644 --- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp +++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp @@ -430,15 +430,15 @@ void PairLJCutTIP4PCut::settings(int narg, char **arg) { if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + 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; - 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_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"); 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; utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.h b/src/MOLECULE/pair_lj_cut_tip4p_cut.h index 75769274a8..a30845274a 100644 --- a/src/MOLECULE/pair_lj_cut_tip4p_cut.h +++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.h @@ -50,6 +50,7 @@ class PairLJCutTIP4PCut : public Pair { double **epsilon, **sigma; 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 typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/MOLECULE/pair_tip4p_cut.cpp b/src/MOLECULE/pair_tip4p_cut.cpp index 73a5651e6b..c914cb0ec0 100644 --- a/src/MOLECULE/pair_tip4p_cut.cpp +++ b/src/MOLECULE/pair_tip4p_cut.cpp @@ -375,12 +375,12 @@ void PairTIP4PCut::settings(int narg, char **arg) { if (narg != 6) error->all(FLERR,"Illegal pair_style command"); - typeO = utils::inumeric(FLERR,arg[0],false,lmp); - typeH = utils::inumeric(FLERR,arg[1],false,lmp); - typeB = utils::inumeric(FLERR,arg[2],false,lmp); - typeA = utils::inumeric(FLERR,arg[3],false,lmp); - qdist = utils::numeric(FLERR,arg[4],false,lmp); - cut_coul = utils::numeric(FLERR,arg[5],false,lmp); + typeO_str = arg[0]; + typeH_str = arg[1]; + typeB_str = arg[2]; + typeA_str = arg[3]; + qdist = utils::numeric(FLERR, arg[4], false, lmp); + cut_coul = utils::numeric(FLERR, arg[5], false, lmp); cut_coulsq = cut_coul * cut_coul; 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"); 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; utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); diff --git a/src/MOLECULE/pair_tip4p_cut.h b/src/MOLECULE/pair_tip4p_cut.h index 73dced7432..e7baeeaa2f 100644 --- a/src/MOLECULE/pair_tip4p_cut.h +++ b/src/MOLECULE/pair_tip4p_cut.h @@ -45,6 +45,7 @@ class PairTIP4PCut : public Pair { double cut_coul, 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 typeA, typeB; // angle and bond types of TIP4P water double alpha; // geometric constraint parameter for TIP4P diff --git a/src/pair.cpp b/src/pair.cpp index 036e21ec28..5421108eba 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1795,8 +1795,8 @@ void Pair::write_file(int narg, char **arg) // parse arguments - int itype = utils::inumeric(FLERR,arg[0],false,lmp); - int jtype = utils::inumeric(FLERR,arg[1],false,lmp); + int itype = utils::expand_type_int(FLERR, arg[0], Atom::ATOM, lmp); + int jtype = utils::expand_type_int(FLERR, arg[1], Atom::ATOM, lmp); if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes) 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"); - double inner = utils::numeric(FLERR,arg[4],false,lmp); - double outer = utils::numeric(FLERR,arg[5],false,lmp); + double inner = utils::numeric(FLERR, arg[4], false, lmp); + double outer = utils::numeric(FLERR, arg[5], false, lmp); if (inner <= 0.0 || inner >= outer) 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 // while printing the date and check that units are consistent. 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)) { error->one(FLERR,"Trying to append to a table file " "with UNITS: {} while units are {}", 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); fp = fopen(table_file.c_str(),"a"); } else { @@ -1847,12 +1847,12 @@ void Pair::write_file(int narg, char **arg) } if (fp == nullptr) 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", - force->pair_style,itype,jtype); + fprintf(fp, "# Pair potential %s for atom types %d %d: i,r,energy,force\n", + force->pair_style, itype, jtype); 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) - 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 @@ -1870,7 +1870,7 @@ void Pair::write_file(int narg, char **arg) double *eamfp_hold; 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)) 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]; q[0] = q[1] = 1.0; if (narg == 10) { - q[0] = utils::numeric(FLERR,arg[8],false,lmp); - q[1] = utils::numeric(FLERR,arg[9],false,lmp); + q[0] = utils::numeric(FLERR, arg[8], false, lmp); + q[1] = utils::numeric(FLERR, arg[9], false, lmp); } double *q_hold; @@ -1893,10 +1893,10 @@ void Pair::write_file(int narg, char **arg) int masklo,maskhi,nmask,nshiftbits; 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; 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; } @@ -1922,7 +1922,7 @@ void Pair::write_file(int narg, char **arg) } 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; } else e = f = 0.0; 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 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 (comm->me == 0) fclose(fp); @@ -2010,4 +2010,3 @@ double Pair::memory_usage() bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double); return bytes; } - diff --git a/src/utils.cpp b/src/utils.cpp index 27d89f1bf5..dc6ee459a7 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -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. - Not guaranteed to return a valid type. - For example, type <= 0 or type > Ntypes is checked in calling routine. + Not guaranteed to return a valid type if param verify = 0 (default) + 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, - LAMMPS *lmp) + LAMMPS *lmp, bool verify) { 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; - return out; + return type; } /* ---------------------------------------------------------------------- diff --git a/src/utils.h b/src/utils.h index 4ba5f0abe2..ebf1709505 100644 --- a/src/utils.h +++ b/src/utils.h @@ -409,7 +409,7 @@ This functions adds the following case to :cpp:func:`utils::bounds()