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

\[\begin{split}\label{eq:lj} V_\mathrm{LJ}(r) = \begin{cases} 4 \epsilon \left[ \left(\frac{\sigma}{r-r_\mathrm{off}}\right)^{12} - \left(\frac{\sigma}{r-r_\mathrm{off}}\right)^6+c_\mathrm{shift}\right] & \mathrm{if~} r_\mathrm{min}+r_\mathrm{off} < r < r_\mathrm{cut}+r_\mathrm{off}\\ 0 & \mathrm{otherwise} \end{cases}.\end{split}\]

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

\[\begin{split}\label{eq:lj-generic} V_\mathrm{LJ}(r) = \begin{cases} \epsilon\left[b_1\left(\frac{\sigma}{r-r_\mathrm{off}}\right)^{e_1} -b_2\left(\frac{\sigma}{r-r_\mathrm{off}}\right)^{e_2}+c_\mathrm{shift}\right] & \mathrm{if~} r_\mathrm{min}+r_\mathrm{off} < r < r_\mathrm{cut}+r_\mathrm{off}\\ 0 & \mathrm{otherwise} \end{cases}\ .\end{split}\]

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

\[\begin{split}\label{eq:wca} V_\mathrm{WCA}(r) = \begin{cases} 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 + \frac{1}{4} \right] & \mathrm{if~} r < \sigma 2^{\frac{1}{6}}\\ 0 & \mathrm{otherwise} \end{cases}.\end{split}\]

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.,

\[V(r)=\frac{1}{2}\epsilon\left(\cos\left[\alpha(r - r_\mathrm{off})^2 + \beta\right]-1\right),\]

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

\[V(r)=-\epsilon\cos^2\left[\frac{\pi}{2\omega}(r - r_\mathrm{min})\right].\]

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

\[V(r)= \left(d/r\right)^n + \epsilon/(1 + \exp\left[2k_0 (r - \sigma)\right])\]

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:

\[V(r)= A\exp\left[B(\sigma - r)\right] - C r^{-6} - D r^{-8} + \epsilon_\mathrm{shift},\]

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

\[V(r)=\epsilon\left(\exp\left[-2 \alpha \left(r - r_\mathrm{min}\right)\right] - 2\exp\left[-\alpha\left(r - r_\mathrm{min}\right)\right]\right) - \epsilon_\mathrm{shift},\]

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

\[V(r)= A \exp(-B r) - C r^{-6} - D r^{-4} + \epsilon_\mathrm{shift}\]

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:

\[V(r)=a\left(r-r_\mathrm{offset}\right)^{-n}\]

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

\[F(r)=F_{\text{max}} \cdot \left( 1 - \frac{r}{r_c} \right),\]

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

The potential energy is given by

\[V(r)=F_{\text{max}} \cdot (r-r_c) \cdot \left( \frac{r+r_c}{2r_c} - 1 \right),\]

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

\[\begin{split}V(r)= \begin{cases} \epsilon\left(1-\frac{r}{\sigma}\right)^{5/2} & r < \sigma\\ 0 & r \ge \sigma. \end{cases}\end{split}\]

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

\[\begin{split}V(r) = \begin{cases} \epsilon\,e^{-\frac{1}{2}\left(\frac{r}{\sigma}\right)^{2}} & r < r_\mathrm{cut}\\ 0 & r \ge r_\mathrm{cut} \end{cases}\end{split}\]

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

\[\vec{F}_{ij}^{D} = -\gamma_\parallel w_\parallel (r_{ij}) (\hat{r}_{ij} \cdot \vec{v}_{ij}) \hat{r}_{ij}\]

for the dissipative force and

\[\vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\parallel w_\parallel (r_{ij}) } \eta_{ij}(t) \hat{r}_{ij}\]

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\):

\[\begin{split}w_\parallel (r_{ij}) = \left\{ \begin{array}{clcr} 1 & , \; \text{weight_function} = 0 \\ {( 1 - (\frac{r_{ij}}{r^\text{cut}_\parallel})^k} )^2 & , \; \text{weight_function} = 1 \end{array} \right.\end{split}\]

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

\[<\eta_{ij}(t)> = 0 , <\eta_{ij}^\alpha(t)\eta_{kl}^\beta(t')> = \delta_{\alpha\beta} \delta_{ik}\delta_{jl}\delta(t-t')\]

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

\[<\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt\]

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

\[\vec{F}_{ij}^{D} = -\gamma_\bot w^D (r_{ij}) (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{v}_{ij}\]

\[\vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\bot w_\bot (r_{ij})} (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{\eta}_{ij}\]

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

\[V(\mathbf{r}_{ij}, \mathbf{\hat{u}}_i, \mathbf{\hat{u}}_j) = 4 \varepsilon(\mathbf{\hat{r}}_{ij}, \mathbf{\hat{u}}_i, \mathbf{\hat{u}}_j) \left( \tilde{r}_{ij}^{-12}-\tilde{r}_{ij}^{-6} \right),\]

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

\[\tilde{r}=\frac{r - \sigma(\mathbf{\hat{r}}, \mathbf{\hat{u}}_i, \mathbf{\hat{u}}_j)+\sigma_0}{\sigma_0},\]

where

\[\sigma( \mathbf{\hat{r}}, \mathbf{\hat{u}}_i, \mathbf{\hat{u}}_j) = \sigma_{0} \left\{ 1 - \frac{1}{2} \chi \left[ \frac{ \left( \mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_i + \mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_j \right)^{2} } {1 + \chi \mathbf{\hat{u}}_i \cdot \mathbf{\hat{u}}_j } + \frac{ \left( \mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_i - \mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_j \right)^{2} } {1 - \chi \mathbf{\hat{u}}_i \cdot \mathbf{\hat{u}}_j} \right] \right\}^{-\frac{1}{2}}\]

and

\[\begin{split}\begin{gathered} \varepsilon(\mathbf{\hat{r}}, \mathbf{\hat{u}}_i, \mathbf{\hat{u}}_j) = \\ \varepsilon_0 \left( 1- \chi^{2}(\mathbf{\hat{u}}_i \cdot \mathbf{\hat{u}}_j)^{2} \right)^{-\frac {\nu}{2}} \left[1-\frac {\chi'}{2} \left( \frac { (\mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_i+ \mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_j)^{2}} {1+\chi' \, \mathbf{\hat{u}}_i \cdot \mathbf{\hat{u}}_j }+ \frac {(\mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_i-\mathbf{\hat{r}} \cdot \mathbf{\hat{u}}_j)^{2}} {1-\chi' \, \mathbf{\hat{u}}_i \cdot \mathbf{\hat{u}}_j } \right) \right]^{\mu}.\end{gathered}\end{split}\]

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\).