6. 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.
Todo
IMPLEMENT: print interaction list
6.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.
Todo
Implement this functionality: If the interaction is omitted, the command returns the currently defined interaction between the two types using the syntax to define the interaction
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.
6.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)\).
6.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 not set or it is set to the string auto, the shift will be automatically computed such that the potential is continuous at the cutoff radius. If is not set, it is set to \(0\).
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, choosing \(r_\mathrm{cut}=2^\frac{1}{6}\sigma\). The WCA potential is purely repulsive, and is often used to mimic hard sphere repulsion.
When coupling particles to a Shan-Chen fluid, if the affinity interaction is set, the Lennard-Jones potential is multiplied by the function
where \(\alpha_i\) is the affinity to the \(i\)-th fluid component (see Affinity interaction), and the order parameter \(\phi\) is calculated from the fluid component local density as \(\phi=\frac{\rho_1 - \rho_2}{\rho_1+\rho_2}\). For example, if the affinities are chosen so that the first component is a good solvent (\(\alpha_1=1\)) and the second one is a bad solvent (\(\alpha_2=0\)), then, if the two particles are both in a region rich in the first component, then \(\phi\simeq1\), and \(A(r)\simeq0\) for \(r>r_\mathrm{cut}+2^{\frac{1}{6}}\sigma\). Therefore, the interaction potential will be very close to the WCA one. Conversely, if both particles are in a region rich in the second component, then \(\phi\simeq-1\), and \(A(r)\simeq 1\), so that the potential will be very close to the full LJ one. If the cutoff has been set large enough, the particle will experience the attractive part of the potential, mimicking the effective attraction induced by the bad solvent.
6.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 net force on a particle can be capped by using force capping system.non_bonded_inter.set_force_cap(max)
, see
section Capping the force during warmup
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.
6.1.4. 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 [SDunwegK01]
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
6.1.5. 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\).
6.1.6. 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}\).
6.1.7. 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\).
6.1.8. 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\).
6.1.9. 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.
6.1.10. Membrane-collision interaction¶
Note
Feature MEMBRANE_COLLISION
required.
This defines a membrane collision interaction between particles of the
types type1
and type2
, where particle of type1
belongs to one OIF or OIF-like object and
particle of type2
belongs to another such object.
It is very similar to soft-sphere interaction, but it takes into account the local outward normal vectors on the surfaces of the two objects to determine the direction for repulsion of objects (i.e. determine whether the two membranes are intersected). It is inversely proportional to the distance of nodes of membranes that are not crossed and saturating with growing distance of nodes of crossed membranes.
In order to work with the OIF objects, both OIF objects need to be created
using OifCellType class with keyword normal=1
, because this implicitly sets up the
bonded out-direction interaction, which computes the outward normal
vector.
The membrane-collision interaction for non-intersected membranes is then defined by:
for \(d<d_\mathrm{cut}\) and \(V(d)=0\) above. For intersected membranes, it is defined as \(V(-d)\). There is no shift implemented currently, which means that the potential is discontinuous at \(d=d_\mathrm{cut}\). Therefore energy calculations should be used with great caution.
6.1.11. Hat interaction¶
Note
Feature HAT
required.
The interface for the Lennard-Jones 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).
6.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 nondifferentiable maximum at \(r=0\), where the force is undefined.
6.1.13. Gaussian¶
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.
6.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, for a general description of the algorithm see there. The parameters can be set via:
system.non_bonded_inter[type1, type2].dpd.set_params(**kwargs)
This command defines a velocity-dependent interaction
between particles of the types type1
and type2
. For a description of the input arguments
see espressomd.interactions.DPDInteraction
. The interaction
only has an effect if the DPD thermostat activated, and acts according to the
temperature of the thermostat.
The interaction consists of a dissipative force \(\vec{F}_{ij}^{D}\) and a random force \(\vec{F}_{ij}^R\), and is decomposed into a component parallel and perpendicular to the distance vector of the particle pair \(\vec{F}_{ij}\). The parameters for the two parts can be chosen independently. The force contributions of the parallel part are
for the dissipative force and
for the random force. Here \(w^D\) and \(w^R\) are weight functions that can be specified via the weight_function parameter of the interaction. The dissipative and random weight function are related by the dissipation-fluctuation theorem:
The possible values for weight_function are 0 and 1, corresponding to the order of \(w^D\):
For the perpendicular part, the dissipative force is calculated by
The random force by
The parameters define the strength of the friction and the cutoff in the same way as above. Note: This interaction does not conserve angular momentum.
6.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()
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 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 LANGEVIN_PER_PARTICLE
#define ROTATION
#define ROTATIONAL_INERTIA
#define ELECTROSTATICS
#define VIRTUAL_SITES_RELATIVE
#define LENNARD_JONES
#define THOLE
#define GHOSTS_HAVE_BONDS
6.2. Anisotropic non-bonded interactions¶
6.2.1. Gay-Berne interaction¶
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 \(\epsilon_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
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\).
6.2.2. Affinity interaction¶
Todo
Not implemented yet.
inter affinity
Instead of defining a new interaction, this command acts as a modifier for existing interactions, so that the conditions of good/bad solvent associated to the two components of a Shan-Chen fluid. The two types must match those of the interaction that one wants to modify, and the two affinity values and are values between 0 and 1. A value of 1 (of 0) indicates that the component acts as a good (bad) solvent. The specific functional form depends on the interaction type and is listed in the interaction section. So far, only the standard Lennard-Jones interaction is modified by the interaction.