# md.many_body¶

Overview

 Triplet Common three body potential documentation. RevCross Reversible crosslinker three-body potential to model bond swaps. SquareDensity Soft potential for simulating a van der Waals liquid. Tersoff Tersoff Potential.

Details

class hoomd.md.many_body.RevCross(nlist, r_cut=None)

Reversible crosslinker three-body potential to model bond swaps.

Parameters

RevCross specifies that the revcross three-body potential should be applied to every non-bonded particle pair in the simulation. 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{eqnarray*} V_{ij}(r) = 4 \varepsilon \left[ \left( \dfrac{ \sigma}{r_{ij}} \right)^{2n} - \left( \dfrac{ \sigma}{r_{ij}} \right)^{n} \right] \qquad r<r_{cut} \end{eqnarray*}

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

\begin{eqnarray*} 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)~,\ \end{eqnarray*}

where the two body potential is rewritten as:

\begin{eqnarray*} \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{eqnarray*}

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')]] = dict(sigma=0.0,n=0,epsilon=0,lambda3=0) and the only non-zero energy only between the different types rev_c.params[('A','B')] = dict(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 towards clusterization because the three-body term is not enough to compensate the energy of multiple bonds, so it may cause unphysical situations.

Use params dictionary to set potential coefficients. The coefficients must be set per unique pair of particle types.

params

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

Type

TypeParameter[tuple[particle_type, particle_type], dict]

Example:

nl = md.nlist.Cell()
bond_swap = md.many_body.RevCross(r_cut=1.3,nlist=nl)
bond_swap.params[(['A','B'],['A','B'])] = dict(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')] = dict(sigma=1,n=100,epsilon=10,lambda3=1)

class hoomd.md.many_body.SquareDensity(nlist, r_cut=None)

Soft potential for simulating a van der Waals liquid.

Parameters

SquareDensity specifies that the three-body potential should be applied to every non-bonded particle pair in the simulation, that is harmonic in the local density.

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)$

Use params dictionary to set potential coefficients. The coefficients must be set per unique pair of particle types.

params

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

• A (float, required) - $$A$$ - mean density (in units of volume^-1, default:0)

• B (float, required) - $$B$$ - coefficient of the harmonic density term (in units of energy*volumne^2)

Type

TypeParameter[tuple[particle_type, particle_type], dict]

Example:

nl = nlist.Cell()
sqd = md.many_body.SquareDensity(nl, 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, 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, r_cut=None)

Tersoff Potential.

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 smoothing function used in this work is exponential in nature as opposed to the sinusoid used by J. Tersoff 1988.

Tersoff specifies that the Tersoff three-body potential should be applied to every non-bonded particle pair in the simulation. 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.

$V_{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}
$$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}$$
$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}

The parameters of this potential are set via the params dictionary, they must be set for each unique pair of particle types.

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 (dimensionless, default: (1.0, 1.0))

• exp_factors (tuple[float, float]) - $$(\lambda_1, \lambda_2)$$ - exponential factors of the repulsive and attractive terms (in units of 1/length, default: 2.0)

• lambda3 (float) - $$\lambda_3$$ - exponential factor in $$\chi_{ij}$$ (in units of 1/length, default: 0.0)

• dimer_r (float) - $$r_D$$ - length shift in attractive and repulsive terms (in units of length, default: 1.5)

• cutoff_thickness (float) - $$r_{CT}$$ - distance which defines the different regions of the potential (in units of length, default: 0.2)

• alpha (float) - $$\alpha$$ - decay rate of the cutoff term $$f_C(r)$$ (dimensionless, default: 3.0)

• n (float) - $$n$$ - power in $$b_{ij}$$ (dimensionless, default: 0.0)

• gamma (float) - $$\gamma$$ - coefficient in $$b_{ij}$$ (dimensionless, default: 0.0)

• c (float) - $$c$$ - coefficient in $$g(\theta)$$ (dimensionless, default: 0.0)

• d (float) - $$d$$ - coefficient in $$g(\theta)$$ (dimensionless, default: 1.0)

• m (float) - $$m$$ - coefficient in $$g(\theta)$$ (dimensionless, default: 0.0)

Type

TypeParameter[tuple[particle_type, particle_type], dict]

Example:

nl = md.nlist.Cell()
tersoff = md.many_body.Tersoff(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, r_cut=None)

Common three body potential documentation.

Users should not invoke Triplet directly. It is a base class that provides common features to all standard triplet forces. Common documentation for all three-body potentials is documented here.

All triplet force commands specify that a given force be computed on all particles which have at least two other particles within a short range cutoff distance $$r_{\mathrm{cut}}$$.

The force $$\vec{F}$$ applied to each particle $$i$$ is:

\begin{eqnarray} \vec{F_i} = & -\nabla V(\vec r_{ij}, \vec r_{ik}) & r_{ij} < r_{\mathrm{cut}} \textrm{ and } r_{ik} < r_{\mathrm{cut}} \\ = & 0 & else \end{eqnarray}

Where $$\vec r_{\alpha \beta} = \vec r_{\alpha} - \vec r_{\beta}$$.

The following coefficients must be set per unique pair of particle types. See hoomd.md.many_body for information on how to set coefficients.

r_cut

r_cut (in distance units), optional: defaults to the value r_cut specified on construction

Type

TypeParameter [ tuple [particle_type, particle_type], float]

Warning

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