\(\renewcommand{\AA}{\text{Å}}\)

pair_style sw/angle/table command

Syntax

pair_style style
  • style = sw/angle/table

Examples

pair_style sw/angle/table
pair_coeff * * spce.sw type

Used in example input script:

examples/PACKAGES/manybody_table/in.spce_sw

Description

Added in version 2Jun2022.

The sw/angle/table style is a modification of the original pair_style sw. It has been developed for coarse-grained simulations (of water) (Scherer1), but can be employed for all kinds of systems. It computes a modified 3-body Stillinger-Weber potential for the energy E of a system of atoms as

\[\begin{split} E & = \sum_i \sum_{j > i} \phi_2 (r_{ij}) + \sum_i \sum_{j \neq i} \sum_{k > j} \phi_3 (r_{ij}, r_{ik}, \theta_{ijk}) \\ \phi_2(r_{ij}) & = A_{ij} \epsilon_{ij} \left[ B_{ij} (\frac{\sigma_{ij}}{r_{ij}})^{p_{ij}} - (\frac{\sigma_{ij}}{r_{ij}})^{q_{ij}} \right] \exp \left( \frac{\sigma_{ij}}{r_{ij} - a_{ij} \sigma_{ij}} \right) \\ \phi_3(r_{ij},r_{ik},\theta_{ijk}) & = f^{\textrm{3b}}\left(\theta_{ijk}\right) \exp \left( \frac{\gamma_{ij} \sigma_{ij}}{r_{ij} - a_{ij} \sigma_{ij}} \right) \exp \left( \frac{\gamma_{ik} \sigma_{ik}}{r_{ik} - a_{ik} \sigma_{ik}} \right)\end{split}\]

where \(\phi_2\) is a two-body term and \(\phi_3\) is a three-body term. The summations in the formula are over all neighbors J and K of atom I within a cutoff distance \(a \sigma\). In contrast to the original sw style, sw/angle/table allows for a flexible three-body term \(f^{\textrm{3b}}\left(\theta_{ijk}\right)\) which is read in as a tabulated interaction. It can be parameterized with the csg_fmatch app of VOTCA as available at: https://gitlab.mpcdf.mpg.de/votca/votca.

Only a single pair_coeff command is used with the sw/angle/table style which specifies a modified Stillinger-Weber potential file with parameters for all needed elements. These are mapped to LAMMPS atom types by specifying N_el additional arguments after the “.sw” filename in the pair_coeff command, where N_el is the number of LAMMPS atom types:

  • “.sw” filename

  • N_el element names = mapping of SW elements to atom types

See the pair_coeff page for alternate ways to specify the path for the potential file.

As an example, imagine a file SiC.sw has Stillinger-Weber values for Si and C. If your LAMMPS simulation has 4 atoms types and you want the first 3 to be Si, and the fourth to be C, you would use the following pair_coeff command:

pair_coeff * * SiC.sw Si Si Si C

The first 2 arguments must be * * so as to span all LAMMPS atom types. The first three Si arguments map LAMMPS atom types 1,2,3 to the Si element in the SW file. The final C argument maps LAMMPS atom type 4 to the C element in the SW file. If a mapping value is specified as NULL, the mapping is not performed. This can be used when a sw/angle/table 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.

The (modified) Stillinger-Weber files have a “.sw” suffix. Lines that are not blank or comments (starting with #) define parameters for a triplet of elements. The parameters in a single entry correspond to the two-body and three-body coefficients in the formula above. Here, also the suffix “.sw” is used though the original Stillinger-Weber file format is supplemented with four additional lines per parameter block to specify the tabulated three-body interaction. A single entry then contains:

  • element 1 (the center atom in a 3-body interaction)

  • element 2

  • element 3

  • \(\epsilon\) (energy units)

  • \(\sigma\) (distance units)

  • a

  • \(\lambda\)

  • \(\gamma\)

  • \(\cos\theta_0\)

  • A

  • B

  • p

  • q

  • tol

  • filename

  • keyword

  • style

  • N

The A, B, p, and q parameters are used only for two-body interactions. The \(\lambda\) and \(\cos\theta_0\) parameters, only used for three-body interactions in the original Stillinger-Weber style, are read in but ignored in this modified pair style. The \(\epsilon\) parameter is only used for two-body interactions in this modified pair style and not for the three-body terms. The \(\sigma\) and a parameters are used for both two-body and three-body interactions. \(\gamma\) is used only in the three-body interactions, but is defined for pairs of atoms. The non-annotated parameters are unitless.

LAMMPS introduces an additional performance-optimization parameter tol that is used for both two-body and three-body interactions. In the Stillinger-Weber potential, the interaction energies become negligibly small at atomic separations substantially less than the theoretical cutoff distances. LAMMPS therefore defines a virtual cutoff distance based on a user defined tolerance tol. The use of the virtual cutoff distance in constructing atom neighbor lists can significantly reduce the neighbor list sizes and therefore the computational cost. LAMMPS provides a tol value for each of the three-body entries so that they can be separately controlled. If tol = 0.0, then the standard Stillinger-Weber cutoff is used.

The additional parameters filename, keyword, style, and N refer to the tabulated angular potential \(f^{\textrm{3b}}\left(\theta_{ijk}\right)\). The tabulated angular potential has to be of the format as used in the angle_style table command:

An interpolation tables of length N is created. The interpolation is done in one of 2 styles: linear or spline. For the linear style, the angle is used to find 2 surrounding table values from which an energy or its derivative is computed by linear interpolation. For the spline style, a cubic spline coefficients are computed and stored at each of the N values in the table. The angle is used to find the appropriate set of coefficients which are used to evaluate a cubic polynomial which computes the energy or derivative.

The filename specifies the file containing the tabulated energy and derivative values of \(f^{\textrm{3b}}\left(\theta_{ijk}\right)\). The keyword then specifies a section of the file. The format of this file is as follows (without the parenthesized comments):

# Angle potential for harmonic (one or more comment or blank lines)

HAM                           (keyword is the first text on line)
N 181 FP 0 0 EQ 90.0          (N, FP, EQ parameters)
                              (blank line)
1 0.0 200.5 2.5               (index, angle, energy, derivative)
2 1.0 198.0 2.5
...
181 180.0 0.0 0.0

A section begins with a non-blank line whose first character is not a “#”; blank lines or lines starting with “#” can be used as comments between sections. The first line begins with a keyword which identifies the section. The next line lists (in any order) one or more parameters for the table. Each parameter is a keyword followed by one or more numeric values.

The parameter “N” is required and its value is the number of table entries that follow. Note that this may be different than the N specified in the Stillinger-Weber potential file. Let Nsw = N in the “.sw” file, and Nfile = “N” in the tabulated angular file. What LAMMPS does is a preliminary interpolation by creating splines using the Nfile tabulated values as nodal points. It uses these to interpolate as needed to generate energy and derivative values at Ntable different points. The resulting tables of length Nsw are then used as described above, when computing energy and force for individual angles and their atoms. This means that if you want the interpolation tables of length Nsw to match exactly what is in the tabulated file (with effectively no preliminary interpolation), you should set Nsw = Nfile.

The “FP” parameter is optional. If used, it is followed by two values fplo and fphi, which are the second derivatives at the innermost and outermost angle settings. These values are needed by the spline construction routines. If not specified by the “FP” parameter, they are estimated (less accurately) by the first two and last two derivative values in the table.

The “EQ” parameter is also optional. If used, it is followed by a the equilibrium angle value, which is used, for example, by the fix shake command. If not used, the equilibrium angle is set to 180.0.

Following a blank line, the next N lines of the angular table file list the tabulated values. On each line, the first value is the index from 1 to N, the second value is the angle value (in degrees), the third value is the energy (in energy units), and the fourth is -dE/d(theta) (also in energy units). The third term is the energy of the 3-atom configuration for the specified angle. The last term is the derivative of the energy with respect to the angle (in degrees, not radians). Thus the units of the last term are still energy, not force. The angle values must increase from one line to the next. The angle values must also begin with 0.0 and end with 180.0, i.e. span the full range of possible angles.

Note that one angular potential file can contain many sections, each with a tabulated potential. LAMMPS reads the file section by section until it finds one that matches the specified keyword of appropriate section of the “.sw” file.

The Stillinger-Weber potential file must contain entries for all the elements listed in the pair_coeff command. It can also contain entries for additional elements not being used in a particular simulation; LAMMPS ignores those entries.

For a single-element simulation, only a single entry is required (e.g. SiSiSi). For a two-element simulation, the file must contain 8 entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that specify SW parameters for all permutations of the two elements interacting in three-body configurations. Thus for 3 elements, 27 entries would be required, etc.

As annotated above, the first element in the entry is the center atom in a three-body interaction. Thus an entry for SiCC means a Si atom with 2 C atoms as neighbors. The parameter values used for the two-body interaction come from the entry where the second and third elements are the same. Thus the two-body parameters for Si interacting with C, comes from the SiCC entry. The three-body angular potential \(f^{\textrm{3b}}\left(\theta_{ijk}\right)\) can in principle be specific to the three elements of the configuration. However, the user must ensure that it makes physically sense. Note also that the function \(\phi_3\) contains two exponential screening factors with parameter values from the ij pair and ik pairs. So \(\phi_3\) for a C atom bonded to a Si atom and a second C atom will depend on the three-body parameters for the CSiC entry, and also on the two-body parameters for the CCC and CSiSi entries. Since the order of the two neighbors is arbitrary, the three-body parameters and the tabulated angular potential for entries CSiC and CCSi should be the same. Similarly, the two-body parameters for entries SiCC and CSiSi should also be the same. The parameters used only for two-body interactions (A, B, p, and q) in entries whose second and third element are different (e.g. SiCSi) are not used for anything and can be set to 0.0 if desired. This is also true for the parameters in \(\phi_3\) that are taken from the ij and ik pairs (\(\sigma\), a, \(\gamma\))

Additional input files and reference data can be found at: https://gitlab.mpcdf.mpg.de/votca/votca/-/tree/master/csg-tutorials/spce/3body_sw


Mixing, shift, table, tail correction, restart, rRESPA info

For atom type pairs I,J and I != J, where types I and J correspond to two different element types, mixing is performed by LAMMPS as described above from values in the potential file, but not for the tabulated angular potential file.

This pair style does not support the pair_modify shift, table, and tail options.

This pair style does not write its information to binary restart files, since it is stored in potential files. Thus, you need to re-specify the pair_style and pair_coeff commands in an input script that reads a restart file.

This pair style can only be used via the pair keyword of the run_style respa command. It does not support the inner, middle, outer keywords.


Restrictions

This pair style is part of the MANYBODY package. It is only enabled if LAMMPS was built with that package. See the Build package page for more info.

This pair style requires the newton setting to be “on” for pair interactions.