Download Notebook View in GitHub Open in Google Colab

Customizing Lennard-Jones Mixing

OpenMM’s standard NonbondedForce class for computing Lennard-Jones and electrostatic interactions implements the Lorentz-Berthelot mixing rules for Lennard-Jones parameters. As described in the user guide, after parameters \(\sigma_i\) and \(\epsilon_i\) have been specified for each particle \(i\) in an OpenMM System, the actual Lennard-Jones parameters used to compute the interaction between two particles \(i\) and \(j\) are

\[\begin{split}\begin{align}\sigma_{ij}&=\frac{\sigma_i+\sigma_j}{2}\\\epsilon_{ij}&=\sqrt{\epsilon_i\epsilon_j}.\end{align}\end{split}\]

You may want to use different mixing rules for Lennard-Jones interactions in a system, or fully customize the Lennard-Jones parameters employed for interactions between pairs of particles in different sets. NonbondedForce supports adding exceptions to specific pairs of particles, but exceptions should not be used to implement custom mixing rules, since exceptions are computed like bonded interactions that do not respect cutoffs.

Instead, you can use CustomNonbondedForce to implement the Lennard-Jones interactions. To handle the typical most general case, in which each particle is assigned one of several particle types, and every pair of particle types may have its own set of parameters, you can use a Discrete2DFunction to store tables of parameter values. To demonstrate this, we first create a System with some particles.

import numpy as np
import openmm

system = openmm.System()
for i in range(300):
    system.addParticle(1.0)

Now we create a CustomNonbondedForce representing a Lennard-Jones interaction, and add it to the System:

lj = openmm.CustomNonbondedForce("4*epsilon(type1, type2)*((sigma(type1, type2)/r)^12 - (sigma(type1, type2)/r)^6)")
system.addForce(lj);

Note that OpenMM has no built-in notion of “particle types”, but by adding a per-particle parameter called type using CustomNonbondedForce.addPerParticleParameter(), the names type1 and type2 in the custom expression will be automatically associated with the values of the type parameter for pairs of interacting particles.

lj.addPerParticleParameter("type");

To assign the actual type values for each particle, we can use CustomNonbondedForce.addParticle(). To make up some values for this example, we’ll set the type of the first 150 particles to 0, the next 100 to 1, and the last 50 to 2.

particle_types = [0] * 150 + [1] * 100 + [2] * 50
for particle_type in particle_types:
    lj.addParticle([particle_type])

Using CustomNonbondedForce.addTabulatedFunction(), we can then add tabulated functions for sigma and epsilon so that values for these parameters can be looked up by pairs of particle type indices.

lj.addTabulatedFunction("sigma", openmm.Discrete2DFunction(3, 3, np.array([
    [1.0, 1.3, 1.4],
    [1.3, 1.1, 1.5],
    [1.4, 1.5, 1.2]
]).flatten(order="F")))
lj.addTabulatedFunction("epsilon", openmm.Discrete2DFunction(3, 3, np.array([
    [1.0, 1.5, 2.0],
    [1.5, 2.0, 2.5],
    [2.0, 2.5, 3.0]
]).flatten(order="F")));

Note that:

• The matrices of parameters are symmetric, so that, e.g., sigma(0, 1) and sigma(1, 0) are equal. When CustomNonbondedForce’s energy expression is evaluated, pairs of particles may appear in any order, and there is no guarantee that they will be ordered consistently. Thus, the energy expression must be symmetric with respect to exchange of particles to obtain sensible results, and so in this example, these matrices must be symmetric.

• Discrete2DFunction takes a flat array of values, in column-major (i.e., “Fortran”) order. Here, therefore, we use 2D NumPy arrays, and flatten them with order="F". Since the arrays are symmetric in this example, this is not strictly needed, but it is a good practice to avoid surprises working with multidimensional tabulated functions in OpenMM.

• This example assumes the use of dimensionless units. Quantities with units passed to OpenMM’s custom force and tabulated function classes will be converted to OpenMM’s default units of nm, ps, amu, kJ/mol, etc. Results when using units should thus be valid as long as all custom expressions given to OpenMM are dimensionally consistent. It is the responsibility of the user to ensure that they are dimensionally consistent, as OpenMM does not verify this.

At this point, we could set any other options we desire on the CustomNonbondedForce, and set up other force field terms in the System as usual. CustomNonbondedForce supports several functions that the standard NonbondedForce supports, including setCutoffDistance(), setNonbondedMethod(), setSwitchingDistance(), setUseLongRangeCorrection(), and setUseSwitchingFunction(). Keep in mind that if your system includes both Lennard-Jones interactions with custom mixing rules and electrostatic interactions, you should use NonbondedForce for the electrostatics, setting all sigma and epsilon values to zero in the NonbondedForce, and then add a CustomNonbondedForce as described here.

Custom Expressions

You can populate the sigma and epsilon tables with arbitrary values, or even use some other mixing rule to calculate them. If, instead of using tabulated values for all particle type pairs, you want to have OpenMM compute values according to some mixing rule from sigma and epsilon values set for each particle, you can include your rule in the CustomNonbondedForce energy expression, and add per-particle parameters for these values instead of particle types, e.g.:

lj = openmm.CustomNonbondedForce(
    "4*epsilon*((sigma/r)^12 - (sigma/r)^6);"

    # Expressions can depend on sigma1, sigma2, epsilon1, and epsilon2
    # and must be symmetric with respect to exchange of 1 and 2
    "sigma = ...; epsilon = ..."
)
lj.addPerParticleParameter("sigma")
lj.addPerParticleParameter("epsilon")

for ...:
    # Specify sigma and epsilon for each particle
    lj.addParticle([..., ...])

Since you can use any energy expression in CustomNonbondedForce, both this approach and the tabular approach above can even be extended to pairwise interactions with arbitrary functional forms, not just Lennard-Jones interactions.

Exceptions and Exclusions

Unlike NonbondedForce, CustomNonbondedForce does not support “exceptions” that permit you to modify interaction parameters for particular pairs of particles, only “exclusions” that allow you to omit pairs entirely. If you only need to zero the Lennard-Jones parameters for certain pairs, you can use CustomNonbondedForce.addExclusion() or CustomNonbondedForce.createExclusionsFromBonds(). If you also want to include modified interactions, you can exclude them from the CustomNonbondedForce and add a CustomBondForce implementing the Lennard-Jones functional form for the pairs of interest. For instance, this could look like:

lj_exceptions = openmm.CustomBondForce("4*epsilon*((sigma/r)^12 - (sigma/r)^6)")
lj_exceptions.addPerBondParameter("sigma")
lj_exceptions.addPerBondParameter("epsilon")
system.addForce(lj_exceptions)

def add_exception(i, j, epsilon, sigma):
    lj.addExclusion(i, j)

    # For LJ interactions, if either sigma or epsilon is zero, the
    # exception force will evaluate to zero, so no bond is needed.
    if sigma and epsilon:
        lj_exceptions.addBond(i, j, [sigma, epsilon])