add pair_style nb3b/screened by Federica Lodesani

doc/src/Commands_pair.rst

@@ -221,6 +221,7 @@ OPT.
    * :doc:`multi/lucy <pair_multi_lucy>`
    * :doc:`multi/lucy/rx (k) <pair_multi_lucy_rx>`
    * :doc:`nb3b/harmonic <pair_nb3b_harmonic>`
+   * :doc:`nb3b/screened <pair_nb3b_screened>`
    * :doc:`nm/cut (o) <pair_nm>`
    * :doc:`nm/cut/coul/cut (o) <pair_nm>`
    * :doc:`nm/cut/coul/long (o) <pair_nm>`

doc/src/pair_nb3b_harmonic.rst

@@ -1,14 +1,20 @@
 .. index:: pair_style nb3b/harmonic
+.. index:: pair_style nb3b/screened
 
 pair_style nb3b/harmonic command
 ================================
 
+pair_style nb3b/screened command
+================================
+
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   pair_style nb3b/harmonic
+   pair_style style
+
+* style = *nb3b/harmonic* or *nb3b/screened*
 
 Examples
 """"""""
@@ -18,10 +24,14 @@ Examples
    pair_style nb3b/harmonic
    pair_coeff * * MgOH.nb3bharmonic Mg O H
 
+   pair_style nb3b/screened
+   pair_coeff * * PO.nb3b.screened P NULL O
+   pair_coeff * * SiOH.nb3b.screened Si O H
+
 Description
 """""""""""
 
-This pair style computes a non-bonded 3-body harmonic potential for the
+The pair style *nb3b/harmonic* computes a non-bonded 3-body harmonic potential for the
 energy E of a system of atoms as
 
 .. math::
@@ -33,7 +43,17 @@ prefactor. Note that the usual 1/2 factor is included in *K*\ . The form
 of the potential is identical to that used in angle_style *harmonic*,
 but in this case, the atoms do not need to be explicitly bonded.
 
-Only a single pair_coeff command is used with this style which
+Style *nb3b/screened* adds an additional exponentially decaying factor to
+the harmonic term, given by
+
+.. math::
+
+   E = K (\theta - \theta_0)^2 \exp \left(- \frac{r_{ij}}{\rho_{ij}} - \frac{r_{ik}}{\rho_{ik}} \right)
+
+where :math:`\rho_ij` and :math:`\rho_ik` are the screening factors along
+the two bonds. Note that the usual 1/2 factor is included in *K*.
+
+Only a single pair_coeff command is used with these styles which
 specifies a potential file with parameters for specified elements.
 These are mapped to LAMMPS atom types by specifying N additional
 arguments after the filename in the pair_coeff command, where N is the
@@ -61,8 +81,8 @@ type 4 to the C element in the potential file.  If a mapping value is
 specified as NULL, the mapping is not performed.  This can be used
 when the potential is used as part of the *hybrid* pair style.  The
 NULL values are placeholders for atom types that will be used with
-other potentials. An example of a pair_coeff command for use with the
+other potentials. Two examples of pair_coeff command for use with the
-*hybrid* pair style is:
+*hybrid* pair style are:
 
 .. code-block:: LAMMPS
 

doc/src/pair_style.rst

@@ -298,6 +298,7 @@ accelerated styles exist.
 * :doc:`multi/lucy <pair_multi_lucy>` - DPD potential with density-dependent force
 * :doc:`multi/lucy/rx <pair_multi_lucy_rx>` - reactive DPD potential with density-dependent force
 * :doc:`nb3b/harmonic <pair_nb3b_harmonic>` - non-bonded 3-body harmonic potential
+* :doc:`nb3b/screened <pair_nb3b_screened>` - non-bonded 3-body screened harmonic potential
 * :doc:`nm/cut <pair_nm>` - N-M potential
 * :doc:`nm/cut/coul/cut <pair_nm>` - N-M potential with cutoff Coulomb
 * :doc:`nm/cut/coul/long <pair_nm>` - N-M potential with long-range Coulomb

examples/nb3b/PSiO.nb3b.screened

@@ -0,0 +1 @@
+../../potentials/PSiO.nb3b.screened

examples/nb3b/README.nb3b.screened

@@ -0,0 +1,41 @@
+Additional information for pair style nb3b/screened example.
+
+This input uses the BMP potential (Bertani, M., Menziani, M. C.,
+& Pedone, A. (2021). Improved empirical force field for multicomponent
+oxide glasses and crystals. Physical Review Materials, 5(4), 045602).
+
+The interatomic potential is obtained by mixing coulombic interactions
+with other two-body functions (modified Morse and Buckingham), with a
+three-body interaction (screened harmonic function).
+
+The modified Morse is introduce with Tables and it is combined with
+Buckingham and coul/dsf to complete the two-body force field.
+The three-body interaction is computed with pair style nb3b/screened
+only for P-O-P angle. With the pair_coeff command the atom types
+involved in the three-body potential are specified:
+
+pair_coeff * * nb3b/screened PSiO.nb3b.screened  P NULL O
+
+And they appear in the same order as in the data file:
+type 1: P;
+type 2: Na, not involved in the three body so indicated as NULL;
+type 3: O.
+
+In the potential file PSiO.nb3b.screened there are the parameters of
+all permutations of the required atom types. In the example, also
+parameters for Si permutations appear but they are not employed in
+the current case.
+
+# i  j  k   K     theta0   rho  cutoff
+O   P   P    65.0  109.47   1.0   3.3
+O   P   Si  120.0  109.47   1.0   0.000
+O   P   O   0.000  0.000    1.0   0.000
+O   Si  P   120.0  109.47   1.0   0.000
+O   Si  Si   25.0  109.47   1.0   3.3
+(...)
+
+The rho value must be always higher than 0.0. Cutoff and rho are extracted
+only from symmetric interactions (e.g. O P P, O Si Si), because of that
+you do not need to specify those values for the asymmetric interactions
+(e.g. O P Si and O Si P).
+

examples/nb3b/in.nb3b.screened

@@ -0,0 +1,35 @@
+### LAMMPS input file
+
+units            metal
+dimension        3
+boundary         p p p
+atom_style       charge
+
+read_data  data.NaPO3
+
+pair_style hybrid/overlay coul/dsf 0.2 12.0 table spline 20000 buck 7.0 nb3b/screened
+
+pair_coeff  * *  coul/dsf
+pair_coeff 1  3 table Table_NP.dat  Pair_P-O 7.0
+pair_coeff 2  3 table Table_NP.dat  Pair_Na-O 7.0
+pair_coeff 3  3 table Table_NP.dat  Pair_O-O 7.0
+
+# pair Buckingham
+pair_coeff  1  1  buck  5.093669  0.905598  0.0
+
+# shrm
+pair_coeff * * nb3b/screened PSiO.nb3b.screened  P NULL O
+
+neighbor         1.0  bin
+neigh_modify     every 10 delay 10 check no
+timestep         0.002
+
+thermo 100
+thermo_style custom step density lx press pe evdwl ecoul temp
+
+#dump 2 all custom 10000 NaPO3-melt.lammpstrj  id type element x y z
+#dump_modify 2 element P  Na  O
+variable temp equal 2500
+
+fix 1 all npt temp ${temp} ${temp} $(100*dt) iso 1 1 $(1000*dt)
+run 1000

examples/nb3b/log.30Oct23.nb3b.screened.g++.1

@@ -0,0 +1,132 @@
+LAMMPS (3 Aug 2023 - Development - patch_2Aug2023-700-g901ed98d31)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
+  using 1 OpenMP thread(s) per MPI task
+### LAMMPS input file
+
+units            metal
+dimension        3
+boundary         p p p
+atom_style       charge
+
+read_data  data.NaPO3
+Reading data file ...
+  orthogonal box = (0 0 0) to (34.33782 34.33782 34.33782)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  3000 atoms
+  reading velocities ...
+  3000 velocities
+  read_data CPU = 0.019 seconds
+
+pair_style hybrid/overlay coul/dsf 0.2 12.0 table spline 20000 buck 7.0 nb3b/screened
+
+pair_coeff  * *  coul/dsf
+pair_coeff 1  3 table Table_NP.dat  Pair_P-O 7.0
+WARNING: 1 of 4995 force values in table Pair_P-O are inconsistent with -dE/dr.
+WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:466)
+pair_coeff 2  3 table Table_NP.dat  Pair_Na-O 7.0
+WARNING: 2 of 4995 force values in table Pair_Na-O are inconsistent with -dE/dr.
+WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:466)
+pair_coeff 3  3 table Table_NP.dat  Pair_O-O 7.0
+
+# pair Buckingham
+pair_coeff  1  1  buck  5.093669  0.905598  0.0
+
+# shrm
+pair_coeff * * nb3b/screened PSiO.nb3b.screened  P NULL O
+Reading nb3b/screened potential file PSiO.nb3b.screened with DATE: 2023-10-30
+
+neighbor         1.0  bin
+neigh_modify     every 10 delay 10 check no
+timestep         0.002
+
+thermo 100
+thermo_style custom step density lx press pe evdwl ecoul temp
+
+#dump 2 all custom 10000 NaPO3-melt.lammpstrj  id type element x y z
+#dump_modify 2 element P  Na  O
+variable temp equal 2500
+
+fix 1 all npt temp ${temp} ${temp} $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 ${temp} $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 0.2000000000000000111 iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 0.2000000000000000111 iso 1 1 2
+run 1000
+Neighbor list info ...
+  update: every = 10 steps, delay = 10 steps, check = no
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 13
+  ghost atom cutoff = 13
+  binsize = 6.5, bins = 6 6 6
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair coul/dsf, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+  (2) pair table, perpetual, skip from (5)
+      attributes: half, newton on, cut 8
+      pair build: skip
+      stencil: none
+      bin: none
+  (3) pair buck, perpetual, skip from (5)
+      attributes: half, newton on, cut 8
+      pair build: skip
+      stencil: none
+      bin: none
+  (4) pair nb3b/screened, perpetual, skip from (6)
+      attributes: full, newton on, cut 4.3
+      pair build: skip
+      stencil: none
+      bin: none
+  (5) neighbor class addition, perpetual, trim from (1)
+      attributes: half, newton on, cut 8
+      pair build: trim
+      stencil: none
+      bin: none
+  (6) neighbor class addition, perpetual
+      attributes: full, newton on, cut 4.3
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 12.11 | 12.11 | 12.11 Mbytes
+   Step        Density           Lx           Press          PotEng         E_vdwl         E_coul          Temp     
+         0   2.509108       34.33782       46464.417     -50581.248      2292.3853     -52873.634      2500         
+       100   2.4533672      34.595928      33156.229     -50182.221      2373.5614     -52555.782      1666.1883    
+       200   2.3649528      35.021791      29182.391     -50158.578      2367.6353     -52526.213      1885.4424    
+       300   2.2932707      35.382952      14450.387     -50100.892      2239.9823     -52340.875      2045.1557    
+       400   2.2282371      35.723887      16680.683     -50048.753      2281.6385     -52330.392      2257.1768    
+       500   2.1753698      36.010969      15871.062     -49981.163      2285.6757     -52266.839      2396.7925    
+       600   2.1285968      36.272824      15066.532     -49934.767      2282.4577     -52217.225      2540.7515    
+       700   2.0841139      36.529076      7572.3436     -49895.93       2222.7473     -52118.677      2577.5774    
+       800   2.0485057      36.739517      6642.7187     -49870.601      2230.3801     -52100.981      2556.7855    
+       900   2.0180051      36.923689      8318.7918     -49890.3        2255.1538     -52145.454      2538.111     
+      1000   1.990678       37.091879      7724.2804     -49928.89       2250.855      -52179.745      2492.3778    
+Loop time of 30.3284 on 1 procs for 1000 steps with 3000 atoms
+
+Performance: 5.698 ns/day, 4.212 hours/ns, 32.972 timesteps/s, 98.917 katom-step/s
+99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 27.347     | 27.347     | 27.347     |   0.0 | 90.17
+Neigh   | 2.8101     | 2.8101     | 2.8101     |   0.0 |  9.27
+Comm    | 0.067332   | 0.067332   | 0.067332   |   0.0 |  0.22
+Output  | 0.00047847 | 0.00047847 | 0.00047847 |   0.0 |  0.00
+Modify  | 0.087652   | 0.087652   | 0.087652   |   0.0 |  0.29
+Other   |            | 0.01585    |            |       |  0.05
+
+Nlocal:           3000 ave        3000 max        3000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:          11846 ave       11846 max       11846 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:         811471 ave      811471 max      811471 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 811471
+Ave neighs/atom = 270.49033
+Neighbor list builds = 100
+Dangerous builds not checked
+Total wall time: 0:00:30

examples/nb3b/log.30Oct23.nb3b.screened.g++.4

@@ -0,0 +1,132 @@
+LAMMPS (3 Aug 2023 - Development - patch_2Aug2023-700-g901ed98d31)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
+  using 1 OpenMP thread(s) per MPI task
+### LAMMPS input file
+
+units            metal
+dimension        3
+boundary         p p p
+atom_style       charge
+
+read_data  data.NaPO3
+Reading data file ...
+  orthogonal box = (0 0 0) to (34.33782 34.33782 34.33782)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  3000 atoms
+  reading velocities ...
+  3000 velocities
+  read_data CPU = 0.011 seconds
+
+pair_style hybrid/overlay coul/dsf 0.2 12.0 table spline 20000 buck 7.0 nb3b/screened
+
+pair_coeff  * *  coul/dsf
+pair_coeff 1  3 table Table_NP.dat  Pair_P-O 7.0
+WARNING: 1 of 4995 force values in table Pair_P-O are inconsistent with -dE/dr.
+WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:466)
+pair_coeff 2  3 table Table_NP.dat  Pair_Na-O 7.0
+WARNING: 2 of 4995 force values in table Pair_Na-O are inconsistent with -dE/dr.
+WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:466)
+pair_coeff 3  3 table Table_NP.dat  Pair_O-O 7.0
+
+# pair Buckingham
+pair_coeff  1  1  buck  5.093669  0.905598  0.0
+
+# shrm
+pair_coeff * * nb3b/screened PSiO.nb3b.screened  P NULL O
+Reading nb3b/screened potential file PSiO.nb3b.screened with DATE: 2023-10-30
+
+neighbor         1.0  bin
+neigh_modify     every 10 delay 10 check no
+timestep         0.002
+
+thermo 100
+thermo_style custom step density lx press pe evdwl ecoul temp
+
+#dump 2 all custom 10000 NaPO3-melt.lammpstrj  id type element x y z
+#dump_modify 2 element P  Na  O
+variable temp equal 2500
+
+fix 1 all npt temp ${temp} ${temp} $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 ${temp} $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 $(100*dt) iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 0.2000000000000000111 iso 1 1 $(1000*dt)
+fix 1 all npt temp 2500 2500 0.2000000000000000111 iso 1 1 2
+run 1000
+Neighbor list info ...
+  update: every = 10 steps, delay = 10 steps, check = no
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 13
+  ghost atom cutoff = 13
+  binsize = 6.5, bins = 6 6 6
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair coul/dsf, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+  (2) pair table, perpetual, skip from (5)
+      attributes: half, newton on, cut 8
+      pair build: skip
+      stencil: none
+      bin: none
+  (3) pair buck, perpetual, skip from (5)
+      attributes: half, newton on, cut 8
+      pair build: skip
+      stencil: none
+      bin: none
+  (4) pair nb3b/screened, perpetual, skip from (6)
+      attributes: full, newton on, cut 4.3
+      pair build: skip
+      stencil: none
+      bin: none
+  (5) neighbor class addition, perpetual, trim from (1)
+      attributes: half, newton on, cut 8
+      pair build: trim
+      stencil: none
+      bin: none
+  (6) neighbor class addition, perpetual
+      attributes: full, newton on, cut 4.3
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 7.188 | 7.188 | 7.189 Mbytes
+   Step        Density           Lx           Press          PotEng         E_vdwl         E_coul          Temp     
+         0   2.509108       34.33782       46464.417     -50581.248      2292.3853     -52873.634      2500         
+       100   2.4533672      34.595928      33156.229     -50182.221      2373.5614     -52555.782      1666.1883    
+       200   2.3649528      35.021791      29182.391     -50158.578      2367.6353     -52526.213      1885.4424    
+       300   2.2932707      35.382952      14450.387     -50100.892      2239.9823     -52340.875      2045.1557    
+       400   2.2282371      35.723887      16680.683     -50048.753      2281.6385     -52330.392      2257.1768    
+       500   2.1753698      36.010969      15871.062     -49981.163      2285.6757     -52266.839      2396.7925    
+       600   2.1285968      36.272824      15066.532     -49934.767      2282.4577     -52217.225      2540.7515    
+       700   2.0841139      36.529076      7572.3446     -49895.929      2222.7474     -52118.677      2577.5774    
+       800   2.0485057      36.739517      6642.7277     -49870.601      2230.3802     -52100.981      2556.7856    
+       900   2.0180051      36.923689      8318.7257     -49890.301      2255.153      -52145.454      2538.1115    
+      1000   1.990678       37.091879      7730.4071     -49928.905      2250.9166     -52179.822      2492.4197    
+Loop time of 9.72954 on 4 procs for 1000 steps with 3000 atoms
+
+Performance: 17.760 ns/day, 1.351 hours/ns, 102.780 timesteps/s, 308.339 katom-step/s
+99.3% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 7.1441     | 8.0123     | 8.622      |  19.6 | 82.35
+Neigh   | 0.7458     | 0.78431    | 0.80001    |   2.5 |  8.06
+Comm    | 0.26061    | 0.88607    | 1.7928     |  61.3 |  9.11
+Output  | 0.00024827 | 0.00027521 | 0.00035196 |   0.0 |  0.00
+Modify  | 0.03802    | 0.038395   | 0.038805   |   0.1 |  0.39
+Other   |            | 0.008186   |            |       |  0.08
+
+Nlocal:            750 ave         765 max         731 min
+Histogram: 1 0 0 0 1 0 0 1 0 1
+Nghost:         6607.5 ave        6671 max        6552 min
+Histogram: 2 0 0 0 0 0 0 0 1 1
+Neighs:         202867 ave      217674 max      182729 min
+Histogram: 1 0 0 0 1 0 0 0 1 1
+
+Total # of neighbors = 811469
+Ave neighs/atom = 270.48967
+Neighbor list builds = 100
+Dangerous builds not checked
+Total wall time: 0:00:09

src/.gitignore

@@ -1297,6 +1297,8 @@
 /pair_morse_soft.h
 /pair_nb3b_harmonic.cpp
 /pair_nb3b_harmonic.h
+/pair_nb3b_screened.cpp
+/pair_nb3b_screened.h
 /pair_nm_cut.cpp
 /pair_nm_cut.h
 /pair_nm_cut_coul_cut.cpp

src/MANYBODY/pair_nb3b_harmonic.cpp

@@ -37,18 +37,19 @@ using MathConst::MY_PI;
 #define DELTA 4
 #define SMALL 0.001
 
+static const char *substyle[] = { "nb3n/harmonic", "nb3b/screened" };
+
 /* ---------------------------------------------------------------------- */
 
-PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp)
+PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp), params(nullptr)
 {
+  variant = HARMONIC;
   single_enable = 0;
   restartinfo = 0;
   one_coeff = 1;
   manybody_flag = 1;
   centroidstressflag = CENTROID_NOTAVAIL;
   unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
-
-  params = nullptr;
 }
 
 /* ----------------------------------------------------------------------
@@ -191,9 +192,10 @@ void PairNb3bHarmonic::coeff(int narg, char **arg)
 
 void PairNb3bHarmonic::init_style()
 {
-  if (atom->tag_enable == 0) error->all(FLERR, "Pair style nb3b/harmonic requires atom IDs");
+  if (atom->tag_enable == 0)
+    error->all(FLERR, "Pair style {} requires atom IDs", substyle[variant]);
   if (force->newton_pair == 0)
-    error->all(FLERR, "Pair style nb3b/harmonic requires newton pair on");
+    error->all(FLERR, "Pair style {} requires newton pair on", substyle[variant]);
 
   // need a full neighbor list
 
@@ -222,14 +224,14 @@ void PairNb3bHarmonic::read_file(char *file)
   // open file on proc 0
 
   if (comm->me == 0) {
-    PotentialFileReader reader(lmp, file, "nb3b/harmonic", unit_convert_flag);
+    PotentialFileReader reader(lmp, file, substyle[variant], unit_convert_flag);
     char *line;
 
     // transparently convert units for supported conversions
 
     int unit_convert = reader.get_unit_convert();
     double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
-    while ((line = reader.next_line(NPARAMS_PER_LINE))) {
+    while ((line = reader.next_line(NPARAMS_PER_LINE + ((variant == SCREENED) ? 1 : 0)))) {
       try {
         ValueTokenizer values(line);
 
@@ -269,6 +271,10 @@ void PairNb3bHarmonic::read_file(char *file)
         params[nparams].kelement = kelement;
         params[nparams].k_theta = values.next_double();
         params[nparams].theta0 = values.next_double();
+        if (variant == SCREENED)
+          params[nparams].rho = values.next_double();
+        else
+          params[nparams].rho = 1.0; // dummy value
         params[nparams].cutoff = values.next_double();
 
         if (unit_convert) params[nparams].k_theta *= conversion_factor;
@@ -276,9 +282,9 @@ void PairNb3bHarmonic::read_file(char *file)
         error->one(FLERR, e.what());
       }
 
-      if (params[nparams].k_theta < 0.0 || params[nparams].theta0 < 0.0 ||
+      if ((params[nparams].k_theta < 0.0) || (params[nparams].theta0 < 0.0) ||
-          params[nparams].cutoff < 0.0)
+          (params[nparams].cutoff < 0.0) || (params[nparams].rho <= 0.0))
-        error->one(FLERR, "Illegal nb3b/harmonic parameter");
+        error->one(FLERR, "Illegal {} parameter", substyle[variant]);
 
       nparams++;
     }

src/MANYBODY/pair_nb3b_harmonic.h

@@ -35,23 +35,26 @@ class PairNb3bHarmonic : public Pair {
   void init_style() override;
 
   static constexpr int NPARAMS_PER_LINE = 6;
+  enum { HARMONIC = 0, SCREENED };
 
  protected:
   struct Param {
     double k_theta, theta0, cutoff;
+    double rho;    // added for screened harmonic style
     double cut, cutsq;
     int ielement, jelement, kelement;
   };
 
   double cutmax;    // max cutoff for all elements
   Param *params;    // parameter set for an I-J-K interaction
+  int variant;
 
   void allocate();
   void read_file(char *);
   void setup_params();
   void twobody(Param *, double, double &, int, double &);
-  void threebody(Param *, Param *, Param *, double, double, double *, double *, double *, double *,
+  virtual void threebody(Param *, Param *, Param *, double, double, double *, double *, double *,
-                 int, double &);
+                         double *, int, double &);
 };
 
 }    // namespace LAMMPS_NS

src/MANYBODY/pair_nb3b_screened.cpp

@@ -0,0 +1,89 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Federica Lodesani
+   (based on nb3b harmonic by Todd R. Zeitler and Stillinger-Weber pair style)
+------------------------------------------------------------------------- */
+
+#include "pair_nb3b_screened.h"
+
+#include <cmath>
+
+#define SMALL 0.001
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairNb3bScreened::PairNb3bScreened(LAMMPS *lmp) : PairNb3bHarmonic(lmp)
+{
+  variant = SCREENED;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairNb3bScreened::threebody(Param *paramij, Param *paramik, Param *paramijk, double rsq1,
+                                 double rsq2, double *delr1, double *delr2, double *fj, double *fk,
+                                 int eflag, double &eng)
+{
+  double dtheta, tk;
+  double r1, r2, c, s, a, a11, a12, a22;
+  double scr, t00, rho1inv, rho2inv;
+  double ratio1, ratio2;
+
+  // angle (cos and sin)
+
+  r1 = sqrt(rsq1);
+  r2 = sqrt(rsq2);
+
+  c = delr1[0] * delr2[0] + delr1[1] * delr2[1] + delr1[2] * delr2[2];
+  c /= r1 * r2;
+
+  if (c > 1.0) c = 1.0;
+  if (c < -1.0) c = -1.0;
+
+  s = sqrt(1.0 - c * c);
+  if (s < SMALL) s = SMALL;
+  s = 1.0 / s;
+
+  // force & energy
+
+  // Harmonic function multiplied by a screening function
+  //
+  // Uijk=k/2(theta-theta0)**2 * exp[-(rij/rhoij+rik/rhoik)]
+  //
+  rho1inv = paramij->rho;
+  rho2inv = paramik->rho;
+  scr = exp(-r1 * rho1inv - r2 * rho2inv);
+
+  dtheta = acos(c) - paramijk->theta0;
+  tk = paramijk->k_theta * dtheta * scr;
+  t00 = tk * dtheta;
+
+  if (eflag) eng = t00;
+
+  a = -2.0 * tk * s;
+  a11 = a * c / rsq1;
+  a12 = -a / (r1 * r2);
+  a22 = a * c / rsq2;
+  ratio1 = rho1inv / r1;
+  ratio2 = rho2inv / r2;
+
+  fj[0] = a11 * delr1[0] + a12 * delr2[0] + t00 * ratio1 * delr1[0];
+  fj[1] = a11 * delr1[1] + a12 * delr2[1] + t00 * ratio1 * delr1[1];
+  fj[2] = a11 * delr1[2] + a12 * delr2[2] + t00 * ratio1 * delr1[2];
+  fk[0] = a22 * delr2[0] + a12 * delr1[0] + t00 * ratio2 * delr2[0];
+  fk[1] = a22 * delr2[1] + a12 * delr1[1] + t00 * ratio2 * delr2[1];
+  fk[2] = a22 * delr2[2] + a12 * delr1[2] + t00 * ratio2 * delr2[2];
+}

src/MANYBODY/pair_nb3b_screened.h

@@ -0,0 +1,39 @@
+/* -*- c++ -*- ---------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(nb3b/screened,PairNb3bScreened);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_NB3B_SCREENED_H
+#define LMP_PAIR_NB3B_SCREENED_H
+
+#include "pair_nb3b_harmonic.h"
+
+namespace LAMMPS_NS {
+
+class PairNb3bScreened : public PairNb3bHarmonic {
+ public:
+  PairNb3bScreened(class LAMMPS *);
+
+ protected:
+  void threebody(Param *, Param *, Param *, double, double, double *, double *, double *, double *,
+                 int, double &) override;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif