# Capped Lennard-Jones potential

## What it does

Driving the Lennard-Jones scaling parameter of a slice to zero with the usual potential leaves a singularity
at $r = 0$ that makes free-energy calculations hard to converge. A scaling parameter added with
`includeLJDecapping=True` multiplies the *decapping* term of its slice, which splits each Lennard-Jones pair
of the slice into a capped term and a decapping term {cite}`Correa_2022`, in the spirit of the capping of
{cite}`Buelens_2012` and the linear basis functions of {cite}`Naden_2014`:

$$
U(r) = \lambda_{LJ}\, u_{cap}(r) + \lambda_D \left[ u_{LJ}(r) - u_{cap}(r) \right],
$$

where $u_{LJ}$ is the usual potential, $\lambda_{LJ}$ is the scaling parameter that multiplies the
Lennard-Jones contributions of the slice (1 if there is none), $\lambda_D$ is the decapping parameter, and

$$
u_{cap}(r) = \begin{cases}
4 \epsilon\, c(r/\sigma), & r < \sigma \\
u_{LJ}(r), & r \ge \sigma
\end{cases}
$$

is a polynomial core. It is defined by two orders, $m_0$ and $m_1$: $c(x)$ is the polynomial of degree
$m_0 + m_1$ whose first $m_0$ derivatives vanish at $x = 0$ and whose value and first $m_1$ derivatives
match those of $x^{-12} - x^{-6}$ at $x = 1$. The first condition leaves the terms $a_0$ and $a_k x^k$ for
$k = m_0 + 1, \dots, m_0 + m_1$, and the second fixes their $m_1 + 1$ coefficients. The core thus joins the
usual potential smoothly, up to the $m_1$-th derivative, at $r = \sigma$ and takes the finite value
$4 \epsilon\, a_0$ at $r = 0$. The default orders, $m_0 = m_1 = 3$, give

$$
c(x) = 29.8 - 360 x^4 + 547.2 x^5 - 217 x^6,
$$

with a core of $119.2\,\epsilon$. The decapping term vanishes for $r \ge \sigma$, and
$\lambda_D = \lambda_{LJ}$ gives the usual potential.

The energy is linear in both parameters. The derivative with respect to $\lambda_{LJ}$ is the capped van der
Waals energy of the slice, the one with respect to $\lambda_D$ is its decapping energy, and their sum is the
usual van der Waals energy ({doc}`derivatives`). One evaluation therefore gives the energy at any pair of
values, which is what reweighting estimators such as MBAR need.

## Setting it up

```python
force.addGlobalParameter("lambdaQ", 1.0)
force.addGlobalParameter("lambdaLJ", 1.0)
force.addGlobalParameter("lambdaD", 1.0)
force.addScalingParameter("lambdaQ", 0, 1, True, False)
force.addScalingParameter("lambdaLJ", 0, 1, False, True)
force.addScalingParameter("lambdaD", 0, 1, False, False, False, True)    # includeLJDecapping
```

The orders are set for the whole force with `setCappingOrders`, from 1 to 6 each. Lowering $m_0$ or raising
$m_1$ raises the core, $4 \epsilon\, a_0$, given here in units of $\epsilon$ (rows: $m_0$; columns: $m_1$):

| | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| **1** | 12.0 | 92.0 | 440.0 | 1640.0 | 5184.0 | 14523.4 |
| **2** | 8.0 | 50.0 | 206.0 | 684.0 | 1969.7 | 5108.1 |
| **3** | 6.0 | 32.4 | 119.2 | 361.3 | 964.4 | 2345.7 |
| **4** | 4.8 | 23.2 | 78.1 | 219.7 | 550.7 | 1268.7 |
| **5** | 4.0 | 17.7 | 55.4 | 146.7 | 348.7 | 766.8 |
| **6** | 3.4 | 14.1 | 41.7 | 104.7 | 237.8 | 502.3 |

```python
force.setCappingOrders(2, 4)    # a core of 684 epsilon
```

## Coulomb decapping

The capped core is finite while the Coulomb attraction of opposite charges is not, so with the decapping term off
they can collapse onto it wherever the Coulomb term is on. One answer is to keep the Coulomb term off until the
decapping term is on; the other is to cap the Coulomb interaction too, with a scaling parameter added with
`includeCoulombDecapping=True`. One parameter can multiply both decapping terms:

```python
force.addScalingParameter("lambdaD", 0, 1, False, False, True, True)    # both decapping terms
```

or each can have its own:

```python
force.addGlobalParameter("lambdaDQ", 1.0)
force.addScalingParameter("lambdaD", 0, 1, False, False, False, True)   # Lennard-Jones decapping
force.addScalingParameter("lambdaDQ", 0, 1, False, False, True)         # Coulomb decapping
```

Then the Coulomb term of every pair in the slice, whatever its $\epsilon$, is capped below
$r = \sigma$ in the same way as the Lennard-Jones one. The $1/r$ part of its direct-space potential is replaced by
$c_Q(r/\sigma)/\sigma$, where $c_Q(x)$ is the polynomial of the same orders $m_0$ and $m_1$ whose first $m_0$
derivatives vanish at $x = 0$ and whose value and first $m_1$ derivatives match those of $1/x$ at $x = 1$. The
default orders give

$$
c_Q(x) = 1.75 - 5.25 x^4 + 7 x^5 - 2.5 x^6,
$$

a core 1.75 times the bare potential at $r = \sigma$. The cap is a statement about $1/r$ alone, as the Lennard-Jones
one is about $x^{-12} - x^{-6}$: the rest of the direct-space potential, which is $-\mathrm{erf}(\alpha r)/r$ with
Ewald, PME and LJPME or the reaction-field terms with a cutoff, and the reciprocal space are unchanged, and the
Coulomb decapping term of the slice is

$$
\frac{q_i q_j}{4 \pi \epsilon_0} \left[ \frac{1}{r} - \frac{c_Q(r/\sigma)}{\sigma} \right]
$$

for $r < \sigma$, which is the same for every nonbonded method and for exceptions. With both interactions decapped,
the potential is bounded below everywhere in the space of the four parameters, as long as the Coulomb decapping term
is not switched on while the Lennard-Jones core is weak, which would let opposite charges collapse. The energy stays
linear in every parameter: the derivative with respect to $\lambda_Q$ is the capped Coulomb energy of the slice (its
capped direct-space part plus the reciprocal one), and the one with respect to the Coulomb decapping parameter is its
Coulomb decapping energy. A slice can also have its Coulomb interaction decapped and its Lennard-Jones one left as
it is, or the reverse.

Since every pair of a Coulomb-decapped slice is capped, $\sigma$ matters for particles with $\epsilon = 0$ too, and the
sigma restriction below covers them. A TIP3P hydrogen, for instance, carries a placeholder $\sigma = 1$ nm that
makes $\sigma_{ij}$ exceed a typical cutoff; set it to the Coulomb cap radius you want (its $\epsilon = 0$ keeps it out of every
Lennard-Jones term), such as zero, which makes $\sigma_{ij}$ half the oxygen's for its pairs. For opposite charges
the capped Coulomb potential is a finite well, deepest at $r = 0$, so decapping should still not lag too far behind
the Coulomb term; raising $m_0$ lowers the core (1.33 at $m_0 = 5$, $m_1 = 2$).

A coupling path for Lennard-Jones decapping alone switches on the capped term first ($\lambda_{LJ}$ from 0 to 1 with $\lambda_D = 0$), then
the decapping term, and finally the Coulomb term. The stages can overlap: {cite}`Correa_2022` couples all
three in one sweep of a master parameter $\lambda$ with smoothstep functions $h_k(\lambda)$ that start and
finish at $(0, 0.7)$, $(0.4, 0.8)$ and $(0.6, 1.0)$ for the capped, decapping and Coulomb terms. OpenMM
cannot make one global parameter a function of another, so set `lambdaLJ`, `lambdaD` and `lambdaQ` from
$\lambda$ at each state; for $\lambda$-dynamics, apply the chain rule to the three derivatives.

The capped term is the usual potential beyond $\sigma$, so the cutoff, the switching function, LJPME and the
dispersion correction apply to it exactly and are multiplied by $\lambda_{LJ}$. The decapping term is not
switched, so the $\sigma$ of every Lennard-Jones pair of a decapped slice (pairs with $\epsilon = 0$
excluded) cannot exceed the cutoff distance, nor the switching distance with a switching function; this is
checked when a Context is created and by `updateParametersInContext`, but not for parameter offsets. Decapping is available on all
platforms and with every nonbonded method, for direct-space pairs and 1-4 exceptions. Like the other
properties of a scaling parameter, the decapping roles and the capping orders are fixed when a
{OpenMM}`Context` is created, while $\lambda_{LJ}$ and $\lambda_D$ are its global parameters.

## Restrictions

* Only one scaling parameter per slice can multiply each of its decapping terms.
* The capped core is finite, so charges of opposite sign can collapse onto it. Turn the Coulomb contribution
  of a slice on only while $\lambda_D > 0$, or cap the Coulomb interaction too with `includeCoulombDecapping`.
* With Coulomb decapping, the $\sigma$ of every pair of charged particles of a decapped slice and of every
  exception with a nonzero charge product, those with $\epsilon = 0$ included, cannot exceed the cutoff
  distance. Coulomb interactions are not switched, so the switching distance does not restrict it.
* With the default orders, the height of the core, $119.2\,\epsilon$, is only a few $k_B T$ for atoms with
  small $\epsilon$, which then overlap often at $\lambda_D = 0$. The decapping energy of such overlaps is
  large, so raise the core with `setCappingOrders` or place enough states near $\lambda_D = 0$ to keep the
  estimates converged. Higher orders give larger coefficients that partly cancel near $r = \sigma$, which
  costs some single-precision accuracy there.
* With an inner force group ({doc}`mts`), the inner group takes all of the capped and decapping terms below
  $\sigma$. The $\sigma$ of every Lennard-Jones pair of a Lennard-Jones-decapped slice (pairs with
  $\epsilon = 0$ excluded) and of every charged pair of a Coulomb-decapped slice cannot exceed the inner
  switching distance. This is checked when a Context is created and by `updateParametersInContext`, but
  not for parameter offsets.