md.many_body

Overview

Triplet

Base class triplet force.

RevCross

Reversible crosslinker three-body force.

SquareDensity

Soft force for simulating a van der Waals liquid.

Tersoff

Tersoff force.

Details

Implement many body potentials.

Triplet force classes apply a force and virial on every particle in the simulation state commensurate with the potential energy:

\[U_\mathrm{many-body} = \frac{1}{2} \sum_{i=0}^\mathrm{N_particles-1} \sum_{j \ne i} \sum_{j \ne k} U(\vec{r}_{ij}, \vec{r}_{ik})\]

where \(\vec{r}_{ij} = \mathrm{minimum\_image}(\vec{r}_j - \vec{r}_i)\). Triplet applies a short range cutoff for performance and assumes that both \(U(\vec{r}_{ij}, \vec{r}_{ik})\) and its derivatives are 0 when \(r_{ij} > r_\mathrm{cut}\) or \(r_{ik} > r_\mathrm{cut}\).

Specifically, the force \(\vec{F}\) applied to each particle \(i\) is:

\[\begin{split}\vec{F_i} = \begin{cases} -\nabla V(\vec r_{ij}, \vec r_{ik}) & r_{ij} < r_{\mathrm{cut}} \land r_{ik} < r_{\mathrm{cut}} \\ 0 & \mathrm{otherwise} \end{cases}\end{split}\]

The per particle energy terms are:

\[U_i = \frac{1}{2} \sum_{j \ne i} \sum_{j \ne k} U(\vec{r}_{ij}, \vec{r}_{ik}) [r_{ij} < r_{\mathrm{cut}} \land r_{ik} < r_{\mathrm{cut}}]\]
class hoomd.md.many_body.RevCross(nlist, default_r_cut=None)

Bases: Triplet

Reversible crosslinker three-body force.

Parameters
  • nlist (hoomd.md.nlist) – Neighbor list

  • default_r_cut (float) – Default cutoff radius \([\mathrm{length}]\).

RevCross computes the revcross three-body force on every particle in the simulation state. Despite the fact that the revcross potential accounts for the effects of third bodies, it is actually just a combination of two body potential terms. It can thus use type-pair parameters similar to those of the pair potentials.

The RevCross potential has been described in detail in S. Ciarella and W.G. Ellenbroek 2019. It is based on a generalized Lennard-Jones pairwise attraction to form bonds between interacting particles:

\[\begin{split}U_{ij}(r) = \begin{cases} 4 \varepsilon \left[ \left( \dfrac{ \sigma}{r_{ij}} \right)^{2n} - \left( \dfrac{ \sigma}{r_{ij}} \right)^{n} \right] & r < r_\mathrm{cut} \\ 0 & r \ge r_\mathrm{cut} \end{cases}\end{split}\]

Then an additional three-body repulsion is evaluated to compensate the bond energies imposing single bond per particle condition:

\[v^{\left( 3b \right)}_{ijk} = \lambda \epsilon \hat{v}^{ \left( 2b \right)}_{ij} \left(\vec{r}_{ij}\right) \cdot \hat{v}^{ \left( 2b \right)}_{ik} \left(\vec{r}_{ik}\right)~,\]

where the two body potential is rewritten as:

\[\begin{split}\hat{v}^{ \left( 2b \right)}_{ij}\left(\vec{r}_{ij}\right) = \begin{cases} 1 & r \le r_{min} \\ - \dfrac{v_{ij}\left(\vec{r}_{ij}\right)}{\epsilon} & r > r_{min} \\ \end{cases}\end{split}\]

Attention

The RevCross potential models an asymmetric interaction between two different chemical moieties that can form a reversible bond. This requires the definition of (at least) two different types of particles. A reversible bond is only possible between two different species, otherwise \(v^{\left( 3b \right)}_{ijk}\), would prevent any bond. In our example we then set the interactions for types A and B with rev_c.params[[('A','B'),('A','B')]] to {"sigma": 0.0, "n": 0, "epsilon": 0, "lambda3": 0} and the only non-zero energy only between the different types with setting rev_c.params[('A','B')] to {"sigma":1, "n": 100, "epsilon": 100, "lambda3": 1}. Notice that the number of the minority species corresponds to the maximum number of bonds.

This three-body term also tunes the energy required for a bond swap through the coefficient:- \(\lambda\) - lambda3 (unitless) in S. Ciarella and W.G. Ellenbroek 2019 is explained that setting \(\lambda=1\) corresponds to no energy requirement to initiate bond swap, while this energy barrier scales roughly as \(\beta \Delta E_\text{sw} =\beta \varepsilon(\lambda-1)\).

Note

Choosing \(\lambda=1\) pushes the system to cluster because the three-body term is not enough to compensate the energy of multiple bonds, so it may cause nonphysical situations.

params

The revcross potential parameters. The dictionary has the following keys:

  • epsilon (float, required) - \(\varepsilon\) \([\mathrm{energy}]\)

  • sigma (float, required) - \(\sigma\) \([\mathrm{length}]\)

  • n (float, required) - \(n\) \([\mathrm{dimensionless}]\)

  • lambda3 (float, required) - \(\lambda_3\) \([\mathrm{dimensionless}]\)

Type: TypeParameter [tuple [particle_type, particle_type], dict]

Example:

nl = md.nlist.Cell()
bond_swap = md.many_body.RevCross(default_r_cut=1.3,nlist=nl)
bond_swap.params[('A', 'A'), ('B', 'B')] = {
    "sigma":0,"n": 0, "epsilon": 0, "lambda3": 0}
# a bond can be made only between A-B and not A-A or B-B
bond_swap.params[('A','B')] = {
    "sigma": 1, "n": 100, "epsilon": 10, "lambda3": 1}
class hoomd.md.many_body.SquareDensity(nlist, default_r_cut=None)

Bases: Triplet

Soft force for simulating a van der Waals liquid.

Parameters
  • nlist (hoomd.md.nlist) – Neighbor list

  • default_r_cut (float) – Default cutoff radius \([\mathrm{length}]\).

SquareDensity that the square density three-body force should on every particle in the simulation state.

The self energy per particle takes the form:

\[\Psi^{ex} = B (\rho - A)^2\]

which gives a pair-wise additive, three-body force:

\[\vec{f}_{ij} = \left( B (n_i - A) + B (n_j - A) \right) w'_{ij} \vec{e}_{ij}\]

Here, \(w_{ij}\) is a quadratic, normalized weighting function,

\[w(x) = \frac{15}{2 \pi r_{c,\mathrm{weight}}^3} (1-r/r_{c,\mathrm{weight}})^2\]

The local density at the location of particle i is defined as

\[n_i = \sum\limits_{j\neq i} w_{ij} \left(\big| \vec r_i - \vec r_j \big|\right)\]
params

The SquareDensity potential parameters. The dictionary has the following keys:

  • A (float, required) - \(A\) - mean density (default:0) \([\mathrm{length}^{-2}]\) in 2D and \([\mathrm{length}^{-3}]\) in 3D

  • B (float, required) - \(B\) - coefficient of the harmonic density term \([\mathrm{energy} \cdot \mathrm{length}^4]\) in 2D and \([\mathrm{energy} \cdot \mathrm{length}^6]\) in 3D

Type: TypeParameter [tuple [particle_type, particle_type], dict]

Example:

nl = nlist.Cell()
sqd = md.many_body.SquareDensity(nl, default_r_cut=3.0)
sqd.params[('A', 'B')] = dict(A=1.0, B=2.0)
sqd.params[('B', 'B')] = dict(A=2.0, B=2.0, default_r_on=1.0)

For further details regarding this multibody potential, see

[1] P. B. Warren, “Vapor-liquid coexistence in many-body dissipative particle dynamics” Phys. Rev. E. Stat. Nonlin. Soft Matter Phys., vol. 68, no. 6 Pt 2, p. 066702, 2003.

class hoomd.md.many_body.Tersoff(nlist, default_r_cut=None)

Bases: Triplet

Tersoff force.

Parameters

The Tersoff potential is a bond-order potential based on the Morse potential that accounts for the weakening of individual bonds with increasing coordination number. It does this by computing a modifier to the attractive term of the potential. The modifier contains the effects of third-bodies on the bond energies. The potential also includes a smoothing function around the cutoff. The implemented smoothing function is exponential in nature as opposed to the sinusoid used by J. Tersoff 1988.

Tersoff computes the Tersoff three-body force on every particle in the simulation state. Despite the fact that the Tersoff potential accounts for the effects of third bodies, the species of the third body is irrelevant. It can thus use type-pair parameters similar to those of the pair potentials:

\[U_{ij}(r) = \frac{1}{2} f_C(r_{ij}) \left[f_R(r_{ij}) + b_{ij}f_A(r_{ij})\right]\]

where

\[ \begin{align}\begin{aligned}f_R(r) = A_1e^{\lambda_1(r_D-r)}\\f_A(r) = A_2e^{\lambda_2(r_D-r)}\end{aligned}\end{align} \]
\[\begin{split}f_C(r) = \begin{cases} 1 & r < r_{\mathrm{cut}} - r_{CT} \\ \exp \left[-\alpha\frac{x(r)^3}{x(r)^3 - 1} \right] & r_{\mathrm{cut}} - r_{CT} < r < r_{\mathrm{cut}} \\ 0 & r > r_{\mathrm{cut}} \end{cases}\end{split}\]
\[b_{ij} = (1 + \gamma^n\chi_{ij}^n)^{\frac{-1}{2n}}\]

In the definition of \(f_C(r)\), there is a quantity \(x(r)\), which is defined as

\[x(r) = \frac{r - (r_{\mathrm{cut}} - r_{CT})}{r_{CT}}\]

which ensures continuity between the different regions of the potential. In the definition of \(b_{ij}\), there is a quantity \(\chi_{ij}\) which is defined as

\[ \begin{align}\begin{aligned}\chi_{ij} = \sum_{k \neq i,j} f_C(r_{ik}) \cdot e^{\lambda_3^3 |r_{ij} - r_{ik}|^3} \cdot g(\theta_{ijk})\\g(\theta_{ijk}) = 1 + \frac{c^2}{d^2} - \frac{c^2}{d^2 + |m - \cos(\theta_{ijk})|^2}\end{aligned}\end{align} \]
params

The Tersoff potential parameters. The dictionary has the following keys:

  • magnitudes (tuple[float, float]) - \((A_1, A_2)\) - Magnitudes of the repulsive and attractive terms (default: (1.0, 1.0)) \([\mathrm{energy}]\)

  • exp_factors (tuple[float, float]) - \((\lambda_1, \lambda_2)\) - exponential factors of the repulsive and attractive terms (default: 2.0) \([\mathrm{length}^{-1}]\)

  • lambda3 (float) - \(\lambda_3\) - exponential factor in \(\chi_{ij}\) (default: 0.0) \([\mathrm{length}^{-1}]\)

  • dimer_r (float) - \(r_D\) - length shift in attractive and repulsive terms (default: 1.5) \([\mathrm{length}]\)

  • cutoff_thickness (float) - \(r_{CT}\) - distance which defines the different regions of the potential (default: 0.2) \([\mathrm{length}]\)

  • alpha (float) - \(\alpha\) - decay rate of the cutoff term \(f_C(r)\) (default: 3.0) \([\mathrm{dimensionless}]\)

  • n (float) - \(n\) - power in \(b_{ij}\) (default: 0.0) \([\mathrm{dimensionless}]\)

  • gamma (float) - \(\gamma\) - coefficient in \(b_{ij}\) (default: 0.0) \([\mathrm{dimensionless}]\)

  • c (float) - \(c\) - coefficient in \(g(\theta)\) (default: 0.0) \([\mathrm{dimensionless}]\)

  • d (float) - \(d\) - coefficient in \(g(\theta)\) (default: 1.0) \([\mathrm{dimensionless}]\)

  • m (float) - \(m\) - coefficient in \(g(\theta)\) (default: 0.0) \([\mathrm{dimensionless}]\)

Type: TypeParameter [tuple [particle_type, particle_type], dict]

Example:

nl = md.nlist.Cell()
tersoff = md.many_body.Tersoff(default_r_cut=1.3, nlist=nl)
tersoff.params[('A', 'B')] = dict(magnitudes=(2.0, 1.0), lambda3=5.0)
class hoomd.md.many_body.Triplet(nlist, default_r_cut=None)

Bases: Force

Base class triplet force.

Note

Triplet is the base class for many-body triplet forces. Users should not instantiate this class directly.

r_cut

r_cut \([\mathrm{length}]\), optional: defaults to the value default_r_cut specified on construction.

Type: TypeParameter [tuple [particle_type, particle_type], float])

nlist

Neighbor list used to compute the triplet potential.

Type: hoomd.md.nlist.NeighborList

Warning

Currently HOOMD-blue does not support reverse force communication between MPI domains on the GPU. Since reverse force communication is required for the calculation of three-body forces, attempting to use this potential on the GPU with MPI will result in an error.