# 7. Non-bonded interactions¶

In *ESPResSo*, interactions are set up and investigated by the `espressomd.interactions`

module. There are
mainly two types of interactions: non-bonded and bonded interactions.

Non-bonded interactions only depend on the *type* of the two particles
involved. This also applies to the electrostatic interaction; however,
due to its long-ranged nature, it requires special care and *ESPResSo* handles it
separately with a number of state-of-the-art algorithms. To specify particle
type and charge see Setting up particles.

A bonded interaction defines an interaction between a number of specific particles; it only applies to the set of particles for which it has been explicitly set. A bonded interaction between a set of particles has to be specified explicitly by the command, while the command is used to define the interaction parameters.

## 7.1. Isotropic non-bonded interactions¶

Non-bonded interaction are configured via the `espressomd.interactions.NonBondedInteraction`

class, which is a member of `espressomd.system.System`

:

```
system.non_bonded_inter[type1, type2]
```

This command defines an interaction between all particles of type `type1`

and
`type2`

. Possible interaction types and their parameters are
listed below.

For many non-bonded interactions, it is possible to artificially cap the forces, which often allows to equilibrate the system much faster. See the subsection Capping the force during warmup for more details.

### 7.1.1. Tabulated interaction¶

Note

Feature `TABULATED`

required.

The interface for tabulated interactions are implemented in the
`TabulatedNonBonded`

class. They can be configured
via the following syntax:

```
system.non_bonded_inter[type1, type2].tabulated.set_params(
min='min', max='max', energy='energy', force='force')
```

This defines an interaction between particles of the types `type1`

and
`type2`

according to an arbitrary tabulated pair potential by linear interpolation.
`force`

specifies the tabulated forces and `energy`

the energies as a function of the
separation distance. `force`

and `energy`

have to have the same length \(N_\mathrm{points}\).
Take care when choosing the number of points, since a copy of each lookup
table is kept on each node and must be referenced very frequently.
The maximal tabulated separation distance also acts as the effective cutoff
value for the potential.

The values of \(r\) are assumed to be equally distributed between \(r_\mathrm{min}\) and \(r_\mathrm{max}\) with a fixed distance of \((r_\mathrm{max}-r_\mathrm{min})/(N_\mathrm{points}-1)\).

### 7.1.2. Lennard-Jones interaction¶

Note

Feature `LENNARD_JONES`

required.

The interface for the Lennard-Jones interaction is implemented in
`LennardJonesInteraction`

. The Lennard-Jones parameters
can be set via:

```
system.non_bonded_inter[type1, type2].lennard_jones.set_params(**kwargs)
```

This command defines the traditional (12-6)-Lennard-Jones interaction
between particles of the types `type1`

and `type2`

. For a description of the input arguments
see `LennardJonesInteraction`

. The potential is defined by

The traditional Lennard-Jones potential is the “work-horse” potential of particle–particle interactions in coarse-grained simulations. It is a simple model for the van-der-Waals interaction, and is attractive at large distance, but strongly repulsive at short distances. \(r_\mathrm{off} + \sigma\) corresponds to the sum of the radii of the interaction particles. At this distance, the potential is \(V_\mathrm{LJ}(r_\mathrm{off} + \sigma) = 4 \epsilon c_\mathrm{shift}\). The minimum of the potential is at \(V_\mathrm{LJ}(r_\mathrm{off} + 2^\frac{1}{6}\sigma) = -\epsilon + 4 \epsilon c_\mathrm{shift}\). Beyond this value the interaction is attractive. Beyond the distance \(r_\mathrm{cut}\) the potential is cut off and the interaction force is zero.

If \(c_\mathrm{shift}\) is set to the string `'auto'`

, the shift will be
automatically computed such that the potential is continuous at the
cutoff radius.

The net force on a particle can be capped by using force capping, see section Capping the force during warmup

An optional additional parameter can be used to restrict the interaction
from a *minimal* distance \(r_\mathrm{min}\). This is an
optional parameter, set to \(0\) by default.

A special case of the Lennard-Jones potential is the
Weeks–Chandler–Andersen (WCA) potential, which one obtains by putting
the cutoff into the minimum and shifting the potential to be continuous, choosing
\(r_\mathrm{cut}=2^\frac{1}{6}\sigma\) and \(c_\mathrm{shift}=\) `'auto'`

. The WCA
potential is purely repulsive, and is often used to mimic hard sphere repulsion.

### 7.1.3. Generic Lennard-Jones interaction¶

Note

Feature `LENNARD_JONES_GENERIC`

required.

The interface for the generic Lennard-Jones interactions is implemented in
`espressomd.interactions.GenericLennardJonesInteraction`

. They
are configured via the syntax:

```
system.non_bonded_inter[type1, type2].generic_lennard_jones.set_params(**kwargs)
```

This command defines a generalized version of the Lennard-Jones
interaction (see Lennard-Jones interaction) between particles of the
types `type1`

and `type2`

. The potential is defined by

Note that the prefactor 4 of the standard LJ potential is missing, so the normal LJ potential is recovered for \(b_1=b_2=4\), \(e_1=12\) and \(e_2=6\).

The optional `LJGEN_SOFTCORE`

feature activates a softcore version of
the potential, where the following transformations apply:
\(\epsilon \rightarrow \lambda \epsilon\) and
\(r-r_\mathrm{off} \rightarrow \sqrt{(r-r_\mathrm{off})^2 +
(1-\lambda) \delta \sigma^2}\). \(\lambda\) allows to tune the strength of the
interaction, while \(\delta\) varies how smoothly the potential goes to zero as
\(\lambda\rightarrow 0\). Such a feature allows one to perform
alchemical transformations, where a group of atoms can be slowly turned
on/off during a simulation.

### 7.1.4. Weeks–Chandler–Andersen interaction¶

Note

Feature `WCA`

required.

The interface for the Weeks–Chandler–Andersen interactions is implemented in
`espressomd.interactions.WCAInteraction`

. They
are configured via the syntax:

```
system.non_bonded_inter[type1, type2].wca.set_params(**kwargs)
```

This command defines a Weeks-Chandler-Andersen interaction between particles of the
types `type1`

and `type2`

. The potential is defined by

### 7.1.5. Lennard-Jones cosine interaction¶

Note

Feature `LJCOS`

and/or `LJCOS2`

required.

```
system.non_bonded_inter[type1, type2].lennard_jones_cos.set_params(**kwargs)
system.non_bonded_inter[type1, type2].lennard_jones_cos2.set_params(**kwargs)
```

`espressomd.interactions.LennardJonesCosInteraction`

and
`espressomd.interactions.LennardJonesCos2Interaction`

specifies
a Lennard-Jones interaction with cosine tail [Soddemann *et al.*, 2001]
between particles of the types `type1`

and `type2`

. The first variant
behaves as follows: Until the minimum of the Lennard-Jones potential
at \(r_\mathrm{min} = r_\mathrm{off} + 2^{\frac{1}{6}}\sigma\), it
behaves identical to the unshifted Lennard-Jones potential
(\(c_\mathrm{shift}=0\)). Between \(r_\mathrm{min}\) and \(r_\mathrm{cut}\), a cosine is used to
smoothly connect the potential to 0, i.e.,

where \(\alpha = \pi\left[(r_\mathrm{cut} - r_\mathrm{off})^2-(r_\mathrm{min} - r_\mathrm{off})^2\right]^{-1}\) and \(\beta = \pi - \left(r_\mathrm{min} - r_\mathrm{off}\right)^2\alpha\).

In the second variant, the cutoff radius is \(r_\mathrm{cut}=r_\mathrm{min} + \omega\), where \(r_\mathrm{min} = r_\mathrm{off} + 2^{\frac{1}{6}}\sigma\) as in the first variant. The potential between \(r_\mathrm{min}\) and \(r_\mathrm{cut}\) is given by

For \(r < r_\mathrm{min}\), \(V(r)\) is implemented as normal Lennard-Jones interaction with \(c_\mathrm{shift} = 0\).

The net force on a particle can be capped by using force capping, see section Capping the force during warmup

### 7.1.6. Smooth step interaction¶

Note

Feature `SMOOTH_STEP`

required.

The interface for the smooth-step interaction is implemented in
`espressomd.interactions.SmoothStepInteraction`

. The smooth-step parameters
can be set via:

```
system.non_bonded_inter[type1, type2].smooth_step.set_params(**kwargs)
```

This defines a smooth step interaction between particles of the types `type1`

and `type2`

, for which the potential is

for \(r<r_\mathrm{cut}\), and \(V(r)=0\) elsewhere. With \(n\) around 10, the first term creates a short range repulsion similar to the Lennard-Jones potential, while the second term provides a much softer repulsion. This potential therefore introduces two length scales, the range of the first term, \(d\), and the range of the second one, \(\sigma\), where in general \(d<\sigma\).

### 7.1.7. BMHTF potential¶

Note

Feature `BMHTF_NACL`

required.

The interface for the smooth-step interaction is implemented in
`espressomd.interactions.BMHTFInteraction`

. The parameters of the BMHTF potential
can be set via:

```
system.non_bonded_inter[type1, type2].bmhtf.set_params(**kwargs)
```

This defines an interaction with the *short-ranged part* of the
Born–Meyer–Huggins–Tosi–Fumi potential between particles of the types `type1`

and `type2`

, which is often used to simulate NaCl crystals. The potential is
defined by:

where \(\epsilon_\mathrm{shift}\) is automatically chosen such that \(V(r_\mathrm{cut})=0\). For \(r\ge r_\mathrm{cut}\), the \(V(r)=0\).

For NaCl, the parameters should be chosen as follows:

types |
\(A\) \(\left(\mathrm{kJ}/\mathrm{mol}\right)\) |
\(B\) \(\left(\mathring{\mathrm{A}}^{-1}\right)\) |
\(C\) \(\left(\mathring{\mathrm{A}}^6 \mathrm{kJ}/\mathrm{mol})\right)\) |
\(D\) \(\left(\mathring{\mathrm{A}}^8 \mathrm{kJ}/\mathrm{mol}\right)\) |
\(\sigma\) \(\left(\mathring{\mathrm{A}}\right)\) |
---|---|---|---|---|---|

Na-Na |
25.4435 |
3.1546 |
101.1719 |
48.1771 |
2.34 |

Na-Cl |
20.3548 |
3.1546 |
674.4793 |
837.0770 |
2.755 |

Cl-Cl |
15.2661 |
3.1546 |
6985.6786 |
14031.5785 |
3.170 |

The cutoff can be chosen relatively freely because the potential decays fast; a value around 10 seems reasonable.

In addition to this short ranged interaction, one needs to add a Coulombic, long-ranged part. If one uses elementary charges, a charge of \(q=+1\) for the Na-particles, and \(q=-1\) for the Cl-particles, the corresponding prefactor of the Coulomb interaction is \(\approx 1389.3549\,\mathrm{kJ}/\mathrm{mol}\).

### 7.1.8. Morse interaction¶

Note

Feature `MORSE`

required.

The interface for the Morse interaction is implemented in
`espressomd.interactions.MorseInteraction`

. The Morse interaction parameters
can be set via:

```
system.non_bonded_inter[type1, type2].morse.set_params(**kwargs)
```

This defines an interaction using the Morse potential between particles
of the types `type1`

and `type2`

. It serves similar purposes as the Lennard-Jones
potential, but has a deeper minimum, around which it is harmonic. This
models the potential energy in a diatomic molecule.

For \(r < r_\mathrm{cut}\), this potential is given by

where is again chosen such that \(V(r_\mathrm{cut})=0\). For \(r\ge r_\mathrm{cut}\), the \(V(r)=0\).

### 7.1.9. Buckingham interaction¶

Note

Feature `BUCKINGHAM`

required.

The interface for the Buckingham interaction is implemented in
`espressomd.interactions.BuckinghamInteraction`

. The Buckingham interaction parameters
can be set via:

```
system.non_bonded_inter[type1, type2].morse.set_params(**kwargs)
```

This defines a Buckingham interaction between particles of the types *type1* and *type2*,
for which the potential is given by

for \(r_\mathrm{discont} < r < r_\mathrm{cut}\). Below \(r_\mathrm{discont}\), the potential is linearly continued towards \(r=0\), similarly to force capping, see below. Above \(r=r_\mathrm{cut}\), the potential is \(0\).

### 7.1.10. Soft-sphere interaction¶

Note

Feature `SOFT_SPHERE`

required.

The interface for the soft-sphere interaction is implemented in
`espressomd.interactions.SoftSphereInteraction`

. The Soft-sphere parameters
can be set via:

```
system.non_bonded_inter[type1, type2].soft_sphere.set_params(**kwargs)
```

This defines a soft-sphere interaction between particles of the types `type1`

and `type2`

, which is defined by a single power law:

for \(r<r_\mathrm{cut}\), and \(V(r)=0\) above. There is no shift implemented currently, which means that the potential is discontinuous at \(r=r_\mathrm{cut}\). Therefore energy calculations should be used with great caution.

### 7.1.11. Hat interaction¶

Note

Feature `HAT`

required.

The interface for the hat interaction is implemented in
`espressomd.interactions.HatInteraction`

. The hat parameters
can be set via:

```
system.non_bonded_inter[type1, type2].hat.set_params(**kwargs)
```

This defines a simple force ramp between particles of two types. The maximal force acts at zero distance and zero force is applied at distances \(r_c\) and bigger. For distances smaller than \(r_c\), the force is given by

for distances exceeding \(r_c\), the force is zero.

The potential energy is given by

which is zero for distances bigger than \(r_c\) and continuous at distance \(0\).

This is the standard conservative DPD potential and can be used in combination with Dissipative Particle Dynamics (DPD).

### 7.1.12. Hertzian interaction¶

Note

Feature `HERTZIAN`

required.

The interface for the Hertzian interaction is implemented in
`espressomd.interactions.HertzianInteraction`

. The Hertzian interaction parameters
can be set via:

```
system.non_bonded_inter[type1, type2].hertzian.set_params(**kwargs)
```

This defines an interaction according to the Hertzian potential between
particles of the types `type1`

and `type2`

. The Hertzian potential is defined by

The potential has no singularity and is defined everywhere; the potential has a non-differentiable maximum at \(r=0\), where the force is undefined.

### 7.1.13. Gaussian interaction¶

Note

Feature `GAUSSIAN`

required.

The interface for the Gaussian interaction is implemented in
`espressomd.interactions.GaussianInteraction`

. The Gaussian interaction parameters
can be set via:

```
system.non_bonded_inter[type1, type2].gaussian.set_params(**kwargs)
```

This defines an interaction according to the Gaussian potential between
particles of the types `type1`

and `type2`

. The Gaussian potential is defined by

The Gaussian potential is smooth except at the cutoff, and has a finite overlap energy of \(\epsilon\). It can be used to model overlapping polymer coils.

Currently, there is no shift implemented, which means that the potential is discontinuous at \(r=r_\mathrm{cut}\). Therefore use caution when performing energy calculations. However, you can often choose the cutoff such that the energy difference at the cutoff is less than a desired accuracy, since the potential decays very rapidly.

### 7.1.14. DPD interaction¶

Note

Feature `DPD`

required.

This is a special interaction that is to be used in conjunction with the Dissipative Particle Dynamics (DPD) thermostat. The parameters can be set via:

```
system.non_bonded_inter[type1, type2].dpd.set_params(**kwargs)
```

This command defines an interaction between particles of the types `type1`

and `type2`

that contains velocity-dependent friction and noise according
to a temperature set by `espressomd.thermostat.Thermostat.set_dpd()`

.
The parameters for the interaction are

`gamma`

`weight_function`

`k`

`r_cut`

`trans_gamma`

`trans_weight_function`

`trans_r_cut`

which will be explained below. The interaction only has an effect if the DPD thermostat is activated.

The interaction consists of a friction force \(\vec{F}_{ij}^{D}\) and a random force \(\vec{F}_{ij}^R\) added to the other interactions between particles \(i\) and \(j\). It is decomposed into a component parallel and perpendicular to the distance vector \(\vec{r}_{ij}\) of the particle pair. The friction force contributions of the parallel part are

for the dissipative force and

for the random force. This introduces the friction coefficient \(\gamma_\parallel\) (parameter `gamma`

) and the weight function
\(w_\parallel\). The thermal energy \(k_B T\) is not set by the interaction,
but by the DPD thermostat (`espressomd.thermostat.Thermostat.set_dpd()`

)
to be equal for all particles. The weight function can be specified via the `weight_function`

switch.
The possible values for `weight_function`

are 0 and 1, corresponding to the
order of \(w_\parallel\):

Both weight functions are set to zero for \(r_{ij}>r^\text{cut}_\parallel\) (parameter `r_cut`

).

In case the `weight_function`

1 is selected the parameter `k`

can be chosen. \(k = 1\) is the
default and recovers the linear ramp. \(k > 1\) enhances the dissipative nature of the interaction
and thus yields higher Schmidt numbers [Yaghoubi *et al.*, 2015].

The random force has the properties

and is numerically discretized to a random number \(\overline{\eta}\) for each spatial component for each particle pair drawn from a uniform distribution with properties

For the perpendicular part, the dissipative and random force are calculated analogously

and introduce the second set of parameters prefixed with `trans_`

.
For \(w_\bot (r_{ij})\) (parameter `trans_weight_function`

)
the same options are available as for \(w_\parallel (r_{ij})\).

Note: This interaction does *not* conserve angular momentum.

A more detailed description of the interaction can be found in [Soddemann *et al.*, 2003].

### 7.1.15. Thole correction¶

Note

Requires features `THOLE`

and `ELECTROSTATICS`

.

Note

`THOLE`

is only implemented for the P3M electrostatics solver.

The Thole correction is closely related to simulations involving Particle polarizability with thermalized cold Drude oscillators. In this context, it is used to correct for overestimation of induced dipoles at short distances. Ultimately, it alters the short-range electrostatics of P3M to result in a damped Coulomb interaction potential \(V(r) = \frac{q_1 q_2}{r} \cdot (1- e^{-s r} (1 + \frac{s r}{2}) )\). The Thole scaling coefficient \(s\) is related to the polarizabilities \(\alpha\) and Thole damping parameters \(a\) of the interacting species via \(s = \frac{ (a_i + a_j) / 2 }{ (\alpha_i \alpha_j)^{1/6} }\). Note that for the Drude oscillators, the Thole correction should be applied only for the dipole part \(\pm q_d\) added by the Drude charge and not on the total core charge, which can be different for polarizable ions. Also note that the Thole correction acts between all dipoles, intra- and intermolecular. Again, the accuracy is related to the P3M accuracy and the split between short-range and long-range electrostatics interaction. It is configured by:

```
system = espressomd.System(box_l=[1, 1, 1])
system.non_bonded_inter[type_1,type_2].thole.set_params(scaling_coeff=<float>, q1q2=<float>)
```

- with parameters:
`scaling_coeff`

: The scaling coefficient \(s\).`q1q2`

: The charge factor of the involved charges.

Because the scaling coefficient depends on the *mixed* polarizabilities and the
nonbonded interaction is controlled by particle types, each Drude charge with a
unique polarizability has to have a unique type. Each Drude charge type has
a Thole correction interaction with all other Drude charges and all Drude
cores, except the one it’s connected to. This exception is handled internally
by disabling Thole interaction between particles connected via Drude bonds.
Also, each Drude core has a Thole correction interaction with all other Drude
cores and Drude charges. To assist with the bookkeeping of mixed scaling
coefficients, the helper method `add_drude_particle_to_core()`

(see
Particle polarizability with thermalized cold Drude oscillators)
collects all core types, Drude types and relevant parameters when a Drude
particle is created. The user already provided all the information when
setting up the Drude particles, so the simple call:

```
add_all_thole(<system>, <verbose>)
```

given the `espressomd.System()`

object, uses this information to create all
necessary Thole interactions. The method calculates the mixed scaling
coefficient `s`

and creates the non-bonded Thole interactions between the
collected types to cover all the Drude-Drude, Drude-core and core-core
combinations. No further calls of `add_drude_particle_to_core()`

should
follow. Set `verbose`

to `True`

to print out the coefficients, charge factors
and involved types.

The samples folder contains the script `/samples/drude_bmimpf6.py`

with a
fully polarizable, coarse grained ionic liquid where this approach is applied.
To use the script, compile espresso with the following features:

```
#define EXTERNAL_FORCES
#define MASS
#define THERMOSTAT_PER_PARTICLE
#define ROTATION
#define ROTATIONAL_INERTIA
#define ELECTROSTATICS
#define VIRTUAL_SITES_RELATIVE
#define LENNARD_JONES
#define THOLE
```

## 7.2. Anisotropic non-bonded interactions¶

### 7.2.1. Gay–Berne interaction¶

Note

Feature `GAY_BERNE`

required.

The interface for a Gay–Berne interaction is provided by the
`espressomd.interactions.GayBerneInteraction`

class. Interaction
parameters can be set via:

```
system.non_bonded_inter[type1, type2].gay_berne.set_params(**kwargs)
```

This defines a Gay–Berne potential for prolate and oblate particles
between particles types `type1`

and `type2`

. The Gay–Berne potential is an
anisotropic version of the classic Lennard-Jones potential, with
orientational dependence of the range \(\sigma_0\) and the well-depth
\(\varepsilon_0\).

Assume two particles with orientations given by the unit vectors \(\mathbf{\hat{u}}_i\) and \(\mathbf{\hat{u}}_j\) and intermolecular vector \(\mathbf{r} = r\mathbf{\hat{r}}\). If \(r<r_\mathrm{cut}\), then the interaction between these two particles is given by

otherwise \(V(r)=0\). The reduced radius is

where

and

Note that contrary to the original paper [Gay and Berne, 1981], here \(\varepsilon_0\) is not raised to the power of \(\nu\), in agreement with the convention used in the Gay–Berne literature.

The parameters \(\chi = \left(k_1^{2} - 1\right)/\left(k_1^{2} + 1\right)\) and \(\chi' = \left(k_2^{1/\mu} - 1\right)/\left(k_2^{1/\mu} + 1\right)\) are responsible for the degree of anisotropy of the molecular properties. \(k_1\) is the molecular elongation, and \(k_2\) is the ratio of the potential well depths for the side-by-side and end-to-end configurations. The exponents and are adjustable parameters of the potential. Several Gay–Berne parametrizations exist, the original one being \(k_1 = 3\), \(k_2 = 5\), \(\mu = 2\) and \(\nu = 1\).