Merge pull request lammps#3960 from akohlmey/nb3b-screened

New pair_style nb3b/screened

.github/CODEOWNERS

@@ -61,6 +61,7 @@ src/GPU/pair_vashishta_gpu.*        @andeplane
 src/KOKKOS/pair_vashishta_kokkos.*  @andeplane
 src/MANYBODY/pair_vashishta_table.* @andeplane
 src/MANYBODY/pair_atm.*             @sergeylishchuk
+src/MANYBODY/pair_nb3b_screened.*   @flodesani
 src/REPLICA/*_grem.*                @dstelter92
 src/EXTRA-COMPUTE/compute_stress_mop*.* @RomainVermorel
 src/EXTRA-COMPUTE/compute_born_matrix.* @Bibobu @athomps

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
-*hybrid* pair style is:
+other potentials. Two examples of pair_coeff command for use with the
+*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    32.5  109.47   1.0   3.3
+O   P   Si   60.0  109.47   1.0   0.000
+O   P   O   0.000  0.000    1.0   0.000
+O   Si  P    60.0  109.47   1.0   0.000
+O   Si  Si   12.5  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).
+