docs/reference-manual/functions/free-energy-interactions.rst

   1 .. _feia:
   2
   3 Free energy interactions
   4 ------------------------
   5
   6 This section describes the :math:`\lambda`-dependence of the potentials
   7 used for free energy calculations (see sec. :ref:`fecalc`). All common
   8 types of potentials and constraints can be interpolated smoothly from
   9 state A (:math:`\lambda=0`) to state B (:math:`\lambda=1`) and vice
  10 versa. All bonded interactions are interpolated by linear interpolation
  11 of the interaction parameters. Non-bonded interactions can be
  12 interpolated linearly or via soft-core interactions.
  13
  14 Starting in |Gromacs| 4.6, :math:`\lambda` is a vector, allowing different
  15 components of the free energy transformation to be carried out at
  16 different rates. Coulomb, Lennard-Jones, bonded, and restraint terms can
  17 all be controlled independently, as described in the
  18 :ref:`mdp` options.
  19
  20 Harmonic potentials
  21 ~~~~~~~~~~~~~~~~~~~
  22
  23 The example given here is for the bond potential, which is harmonic in
  24 |Gromacs|. However, these equations apply to the angle potential and the
  25 improper dihedral potential as well.
  26
  27 .. math:: \begin{aligned}
  28           V_b     &=&{\frac{1}{2}}\left[{(1-{\lambda})}k_b^A +
  29                           {\lambda}k_b^B\right] \left[b - {(1-{\lambda})}b_0^A - {\lambda}b_0^B\right]^2  \\
  30           {\frac{\partial V_b}{\partial {\lambda}}}&=&{\frac{1}{2}}(k_b^B-k_b^A)
  31                           \left[b - {(1-{\lambda})}b_0^A + {\lambda}b_0^B\right]^2 +
  32                         \nonumber\\
  33                   & & \phantom{{\frac{1}{2}}}(b_0^A-b_0^B) \left[b - {(1-{\lambda})}b_0^A -{\lambda}b_0^B\right]
  34                         \left[{(1-{\lambda})}k_b^A + {\lambda}k_b^B \right]\end{aligned}
  35           :label: eqnfepharmpot
  36
  37 GROMOS-96 bonds and angles
  38 ~~~~~~~~~~~~~~~~~~~~~~~~~~
  39
  40 Fourth-power bond stretching and cosine-based angle potentials are
  41 interpolated by linear interpolation of the force constant and the
  42 equilibrium position. Formulas are not given here.
  43
  44 Proper dihedrals
  45 ~~~~~~~~~~~~~~~~
  46
  47 For the proper dihedrals, the equations are somewhat more complicated:
  48
  49 .. math:: \begin{aligned}
  50           V_d     &=&\left[{(1-{\lambda})}k_d^A + {\lambda}k_d^B \right]
  51                   \left( 1+ \cos\left[n_{\phi} \phi -
  52                             {(1-{\lambda})}\phi_s^A - {\lambda}\phi_s^B
  53                             \right]\right)\\
  54           {\frac{\partial V_d}{\partial {\lambda}}}&=&(k_d^B-k_d^A)
  55                    \left( 1+ \cos
  56                          \left[
  57                             n_{\phi} \phi- {(1-{\lambda})}\phi_s^A - {\lambda}\phi_s^B
  58                          \right]
  59                  \right) +
  60                  \nonumber\\
  61                   &&(\phi_s^B - \phi_s^A) \left[{(1-{\lambda})}k_d^A - {\lambda}k_d^B\right]
  62                   \sin\left[  n_{\phi}\phi - {(1-{\lambda})}\phi_s^A - {\lambda}\phi_s^B \right]\end{aligned}
  63           :label: eqnfeppropdihedral
  64
  65 **Note:** that the multiplicity :math:`n_{\phi}` can not be
  66 parameterized because the function should remain periodic on the
  67 interval :math:`[0,2\pi]`.
  68
  69 Tabulated bonded interactions
  70 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  71
  72 For tabulated bonded interactions only the force constant can
  73 interpolated:
  74
  75 .. math:: \begin{aligned}
  76                 V  &=& ({(1-{\lambda})}k^A + {\lambda}k^B) \, f \\
  77           {\frac{\partial V}{\partial {\lambda}}} &=& (k^B - k^A) \, f\end{aligned}
  78           :label: eqnfeptabbonded
  79
  80 Coulomb interaction
  81 ~~~~~~~~~~~~~~~~~~~
  82
  83 The Coulomb interaction between two particles of which the charge varies
  84 with :math:`{\lambda}` is:
  85
  86 .. math:: \begin{aligned}
  87           V_c &=& \frac{f}{{\varepsilon_{rf}}{r_{ij}}}\left[{(1-{\lambda})}q_i^A q_j^A + {\lambda}\, q_i^B q_j^B\right] \\
  88           {\frac{\partial V_c}{\partial {\lambda}}}&=& \frac{f}{{\varepsilon_{rf}}{r_{ij}}}\left[- q_i^A q_j^A + q_i^B q_j^B\right]\end{aligned}
  89           :label: eqnfepcoloumb
  90
  91 where :math:`f = \frac{1}{4\pi \varepsilon_0} = {138.935\,458}` (see
  92 chapter :ref:`defunits`).
  93
  94 Coulomb interaction with reaction field
  95 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  96
  97 The Coulomb interaction including a reaction field, between two
  98 particles of which the charge varies with :math:`{\lambda}` is:
  99
 100 .. math:: \begin{aligned}
 101           V_c     &=& f\left[\frac{1}{{r_{ij}}} + k_{rf}~ {r_{ij}}^2 -c_{rf}\right]
 102           \left[{(1-{\lambda})}q_i^A q_j^A + {\lambda}\, q_i^B q_j^B\right] \\
 103           {\frac{\partial V_c}{\partial {\lambda}}}&=& f\left[\frac{1}{{r_{ij}}} + k_{rf}~ {r_{ij}}^2 -c_{rf}\right]
 104           \left[- q_i^A q_j^A + q_i^B q_j^B\right]
 105           \end{aligned}
 106           :label: eqdVcoulombdlambda
 107
 108 **Note** that the constants :math:`k_{rf}` and :math:`c_{rf}` are
 109 defined using the dielectric constant :math:`{\varepsilon_{rf}}` of the
 110 medium (see sec. :ref:`coulrf`).
 111
 112 Lennard-Jones interaction
 113 ~~~~~~~~~~~~~~~~~~~~~~~~~
 114
 115 For the Lennard-Jones interaction between two particles of which the
 116 *atom type* varies with :math:`{\lambda}` we can write:
 117
 118 .. math:: \begin{aligned}
 119           V_{LJ}  &=&     \frac{{(1-{\lambda})}C_{12}^A + {\lambda}\, C_{12}^B}{{r_{ij}}^{12}} -
 120           \frac{{(1-{\lambda})}C_6^A + {\lambda}\, C_6^B}{{r_{ij}}^6}   \\
 121           {\frac{\partial V_{LJ}}{\partial {\lambda}}}&=&\frac{C_{12}^B - C_{12}^A}{{r_{ij}}^{12}} -
 122           \frac{C_6^B - C_6^A}{{r_{ij}}^6}
 123           \end{aligned}
 124           :label: eqdVljdlambda
 125
 126 It should be noted that it is also possible to express a pathway from
 127 state A to state B using :math:`\sigma` and :math:`\epsilon` (see
 128 :eq:`eqn. %s <eqnsigeps>`). It may seem to make sense physically to vary the
 129 force field parameters :math:`\sigma` and :math:`\epsilon` rather than
 130 the derived parameters :math:`C_{12}` and :math:`C_{6}`. However, the
 131 difference between the pathways in parameter space is not large, and the
 132 free energy itself does not depend on the pathway, so we use the simple
 133 formulation presented above.
 134
 135 Kinetic Energy
 136 ~~~~~~~~~~~~~~
 137
 138 When the mass of a particle changes, there is also a contribution of the
 139 kinetic energy to the free energy (note that we can not write the
 140 momentum :math:`\mathbf{p}` as
 141 m :math:`\mathbf{v}`, since that would result in the
 142 sign of :math:`{\frac{\partial E_k}{\partial {\lambda}}}` being
 143 incorrect \ :ref:`99 <refGunsteren98a>`):
 144
 145 .. math:: \begin{aligned}
 146           E_k      &=&     {\frac{1}{2}}\frac{\mathbf{p}^2}{{(1-{\lambda})}m^A + {\lambda}m^B}        \\
 147           {\frac{\partial E_k}{\partial {\lambda}}}&=&    -{\frac{1}{2}}\frac{\mathbf{p}^2(m^B-m^A)}{({(1-{\lambda})}m^A + {\lambda}m^B)^2}\end{aligned}
 148           :label: eqnfepekin
 149
 150 after taking the derivative, we *can* insert
 151 :math:`\mathbf{p}` = m :math:`\mathbf{v}`, such that:
 152
 153 .. math:: {\frac{\partial E_k}{\partial {\lambda}}}~=~    -{\frac{1}{2}}\mathbf{v}^2(m^B-m^A)
 154           :label: eqnfepekinderivative
 155
 156 Constraints
 157 ~~~~~~~~~~~
 158
 159 The constraints are formally part of the Hamiltonian, and therefore they
 160 give a contribution to the free energy. In |Gromacs| this can be
 161 calculated using the LINCS or the SHAKE algorithm. If we have
 162 :math:`k = 1 \ldots K` constraint equations :math:`g_k` for LINCS, then
 163
 164 .. math:: g_k     =       | \mathbf{r}_{k} | - d_{k}
 165           :label: eqnfepconstr
 166
 167 where :math:`\mathbf{r}_k` is the displacement vector
 168 between two particles and :math:`d_k` is the constraint distance between
 169 the two particles. We can express the fact that the constraint distance
 170 has a :math:`{\lambda}` dependency by
 171
 172 .. math:: d_k     =       {(1-{\lambda})}d_{k}^A + {\lambda}d_k^B
 173           :label: eqnfepconstrdistdep
 174
 175 Thus the :math:`{\lambda}`-dependent constraint equation is
 176
 177 .. math:: g_k     =       | \mathbf{r}_{k} | - \left({(1-{\lambda})}d_{k}^A + {\lambda}d_k^B\right).
 178           :label: eqnfepconstrlambda
 179
 180 The (zero) contribution :math:`G` to the Hamiltonian from the
 181 constraints (using Lagrange multipliers :math:`\lambda_k`, which are
 182 logically distinct from the free-energy :math:`{\lambda}`) is
 183
 184 .. math:: \begin{aligned}
 185           G           &=&     \sum^K_k \lambda_k g_k    \\
 186           {\frac{\partial G}{\partial {\lambda}}}    &=&     \frac{\partial G}{\partial d_k} {\frac{\partial d_k}{\partial {\lambda}}} \\
 187                       &=&     - \sum^K_k \lambda_k \left(d_k^B-d_k^A\right)\end{aligned}
 188           :label: eqnconstrfreeenergy
 189
 190 For SHAKE, the constraint equations are
 191
 192 .. math:: g_k     =       \mathbf{r}_{k}^2 - d_{k}^2
 193           :label: eqnfepshakeconstr
 194
 195 with :math:`d_k` as before, so
 196
 197 .. math:: \begin{aligned}
 198           {\frac{\partial G}{\partial {\lambda}}}    &=&     -2 \sum^K_k \lambda_k \left(d_k^B-d_k^A\right)\end{aligned}
 199           :label: eqnfepshakeconstr2
 200
 201 Soft-core interactions: Beutler *et al.*
 202 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 203
 204 .. _fig-softcore:
 205
 206 .. figure:: plots/softcore.*
 207    :height: 6.00000cm
 208
 209    Soft-core interactions at :math:`{\lambda}=0.5`, with :math:`p=2` and
 210    :math:`C_6^A=C_{12}^A=C_6^B=C_{12}^B=1`.
 211
 212 In a free-energy calculation where particles grow out of nothing, or
 213 particles disappear, using the simple linear interpolation of the
 214 Lennard-Jones and Coulomb potentials as described in
 215 :eq:`Equations %s <eqdVljdlambda>` and :eq:`%s <eqdVcoulombdlambda>` may lead to poor
 216 convergence. When the particles have nearly disappeared, or are close to
 217 appearing (at :math:`{\lambda}` close to 0 or 1), the interaction energy
 218 will be weak enough for particles to get very close to each other,
 219 leading to large fluctuations in the measured values of
 220 :math:`\partial V/\partial {\lambda}` (which, because of the simple
 221 linear interpolation, depends on the potentials at both the endpoints of
 222 :math:`{\lambda}`).
 223
 224 To circumvent these problems, the singularities in the potentials need
 225 to be removed. This can be done by modifying the regular Lennard-Jones
 226 and Coulomb potentials with “soft-core” potentials that limit the
 227 energies and forces involved at :math:`{\lambda}` values between 0 and
 228 1, but not *at* :math:`{\lambda}=0` or 1.
 229
 230 In |Gromacs| the soft-core potentials :math:`V_{sc}` are shifted versions
 231 of the regular potentials, so that the singularity in the potential and
 232 its derivatives at :math:`r=0` is never reached:
 233
 234 .. math:: \begin{aligned}
 235           V_{sc}(r) &=& {(1-{\lambda})}V^A(r_A) + {\lambda}V^B(r_B)
 236               \\
 237           r_A &=& \left(\alpha \sigma_A^6 {\lambda}^p + r^6 \right)^\frac{1}{6}
 238               \\
 239           r_B &=& \left(\alpha \sigma_B^6 {(1-{\lambda})}^p + r^6 \right)^\frac{1}{6}\end{aligned}
 240           :label: eqnfepsoftcore
 241
 242 where :math:`V^A` and :math:`V^B` are the normal “hard core” Van der
 243 Waals or electrostatic potentials in state A (:math:`{\lambda}=0`) and
 244 state B (:math:`{\lambda}=1`) respectively, :math:`\alpha` is the
 245 soft-core parameter (set with ``sc_alpha`` in the
 246 :ref:`mdp` file), :math:`p` is the soft-core :math:`{\lambda}`
 247 power (set with ``sc_power``), :math:`\sigma` is the radius
 248 of the interaction, which is :math:`(C_{12}/C_6)^{1/6}` or an input
 249 parameter (``sc_sigma``) when :math:`C_6` or :math:`C_{12}`
 250 is zero.
 251
 252 For intermediate :math:`{\lambda}`, :math:`r_A` and :math:`r_B` alter
 253 the interactions very little for :math:`r > \alpha^{1/6} \sigma` and
 254 quickly switch the soft-core interaction to an almost constant value for
 255 smaller :math:`r` (:numref:`Fig. %s <fig-softcore>`). The force is:
 256
 257 .. math:: F_{sc}(r) = -\frac{\partial V_{sc}(r)}{\partial r} =
 258            {(1-{\lambda})}F^A(r_A) \left(\frac{r}{r_A}\right)^5 +
 259           {\lambda}F^B(r_B) \left(\frac{r}{r_B}\right)^5
 260           :label: eqnfepsoftcoreforce
 261
 262 where :math:`F^A` and :math:`F^B` are the “hard core” forces. The
 263 contribution to the derivative of the free energy is:
 264
 265 .. math:: \begin{aligned}
 266           {\frac{\partial V_{sc}(r)}{\partial {\lambda}}} & = &
 267            V^B(r_B) -V^A(r_A)  +
 268                 {(1-{\lambda})}\frac{\partial V^A(r_A)}{\partial r_A}
 269                            \frac{\partial r_A}{\partial {\lambda}} +
 270                 {\lambda}\frac{\partial V^B(r_B)}{\partial r_B}
 271                            \frac{\partial r_B}{\partial {\lambda}}
 272           \nonumber\\
 273           &=&
 274            V^B(r_B) -V^A(r_A)  + \nonumber \\
 275            & &
 276            \frac{p \alpha}{6}
 277                  \left[ {\lambda}F^B(r_B) r^{-5}_B \sigma_B^6 {(1-{\lambda})}^{p-1} -
 278                        {(1-{\lambda})}F^A(r_A) r^{-5}_A \sigma_A^6 {\lambda}^{p-1} \right]\end{aligned}
 279           :label: eqnfepsoftcorederivative
 280
 281 The original GROMOS Lennard-Jones soft-core
 282 function\ :ref:`100 <refBeutler94>` uses :math:`p=2`, but :math:`p=1` gives a smoother
 283 :math:`\partial H/\partial{\lambda}` curve. Another issue that should be
 284 considered is the soft-core effect of hydrogens without Lennard-Jones
 285 interaction. Their soft-core :math:`\sigma` is set with
 286 ``sc_sigma`` in the :ref:`mdp` file. These
 287 hydrogens produce peaks in :math:`\partial H/\partial{\lambda}` at
 288 :math:`{\lambda}` is 0 and/or 1 for :math:`p=1` and close to 0 and/or 1
 289 with :math:`p=2`. Lowering ``sc_sigma``
 290 will decrease this effect, but it will also increase the interactions
 291 with hydrogens relative to the other interactions in the soft-core
 292 state.
 293
 294 When soft-core potentials are selected (by setting
 295 ``sc_alpha >0``), and the Coulomb and Lennard-Jones
 296 potentials are turned on or off sequentially, then the Coulombic
 297 interaction is turned off linearly, rather than using soft-core
 298 interactions, which should be less statistically noisy in most cases.
 299 This behavior can be overwritten by using the :ref:`mdp` option
 300 ``sc-coul`` to ``yes``. Note that the
 301 ``sc-coul`` is only taken into account when lambda states
 302 are used, not with ``couple-lambda0``  /
 303 ``couple-lambda1``, and you can still turn off soft-core
 304 interactions by setting ``sc-alpha=0``. Additionally, the
 305 soft-core interaction potential is only applied when either the A or B
 306 state has zero interaction potential. If both A and B states have
 307 nonzero interaction potential, default linear scaling described above is
 308 used. When both Coulombic and Lennard-Jones interactions are turned off
 309 simultaneously, a soft-core potential is used, and a hydrogen is being
 310 introduced or deleted, the sigma is set to ``sc-sigma-min``,
 311 which itself defaults to ``sc-sigma-default``.
 312
 313 Recently, a new formulation of the soft-core approach has been derived
 314 that in most cases gives lower and more even statistical variance than
 315 the standard soft-core path described above \ :ref:`101 <refPham2011>`,
 316 :ref:`102 <refPham2012>`. Specifically, we have:
 317
 318 .. math:: \begin{aligned}
 319           V_{sc}(r) &=& {(1-{\lambda})}V^A(r_A) + {\lambda}V^B(r_B)
 320               \\
 321           r_A &=& \left(\alpha \sigma_A^{48} {\lambda}^p + r^{48} \right)^\frac{1}{48}
 322               \\
 323           r_B &=& \left(\alpha \sigma_B^{48} {(1-{\lambda})}^p + r^{48} \right)^\frac{1}{48}\end{aligned}
 324           :label: eqnnewsoftcore
 325
 326 This “1-1-48” path is also implemented in |Gromacs|. Note that for this
 327 path the soft core :math:`\alpha` should satisfy
 328 :math:`0.001 < \alpha < 0.003`, rather than :math:`\alpha \approx
 329 0.5`.
 330
 331
 332 Soft-core interactions: Gapsys *et al.*
 333 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 334
 335 In this section we describe the functional form and parameters for
 336 the soft-cored non-bonded interactions using the formalism by Gapsys *et al.*\ :ref:`185 <refGapsys2012>`.
 337
 338 The Gapsys *et al.* soft-core is formulated to act on the level of van der Waals and electrostatic forces:
 339 the non-bonded interactions are linearized at a point defined as, :math:`r_{scLJ}` or :math:`r_{scQ}`, respectively.
 340 The linearization point depends on the state of the system as controlled by the :math:`\lambda` parameter and
 341 two parameters :math:`\alpha_Q` (set with ``sc-gapsys-scale-linpoint-q``) and :math:`\alpha_{LJ}` (set with ``sc-gapsys-scale-linpoint-lj``).
 342 The dependence on :math:`\lambda` guarantees that the end-states are properly represented by their hard-core potentials.
 343 :numref:`Fig. %s <fig-gapsyssc>` illustrates the behaviour of the linearization point, forces and integrated potential energies with respect
 344 to the parameters :math:`\alpha_Q` and :math:`\alpha_{LJ}`. The optimal choices of the parameter values have been systematically explored in :ref:`185 <refGapsys2012>`. These recommended values are set by default when ``sc-function=gapsys`` is selected: ``sc-gapsys-scale-linpoint-q=0.3`` and ``sc-gapsys-scale-linpoint-lj=0.85``.
 345
 346 .. _fig-gapsyssc:
 347
 348 .. figure:: plots/gapsys-sc.*
 349         :width: 15.0cm
 350
 351         Illustration of the soft-core parameter influence on the linearization point (top row),
 352         forces (middle row) and energies (bottom row)
 353         for van der Waals (left column) and electrostatic interactions (right column).
 354         The case of two interacting atoms is considered.
 355         In state A both atoms have charges of 0.5 and :math:`\sigma=0.3` nm, :math:`\epsilon=0.5` kJ/mol.
 356         In state B all the non-bonded interactions are set to zero.
 357         The parameter :math:`\lambda` is set to 0.5 and electrostatic interaction cutoff is 1 nm.
 358
 359 The parameter :math:`\alpha_{LJ}` is a unitless scaling factor in the range :math:`[0,1)`.
 360 It scales the position of the point from which the van der Waals force will be linearized.
 361 The linearization of the force is allowed in the range :math:`[0,F_{min}^{LJ})`,
 362 where setting :math:`\alpha_{LJ}=0` results in a standard hard-core van der Waals interaction.
 363 Setting :math:`\alpha_{LJ}` closer to 1 brings the force linearization point towards
 364 the minimum in the Lennard-Jones force curve (:math:`F_{min}^{LJ}`).
 365 This construct allows retaining the repulsion between two particles with non-zero C12 parameter at any :math:`\lambda` value.
 366
 367 The parameter :math:`\alpha_{Q}` has a unit of :math:`\frac{nm}{e^2}` and is defined in the range :math:`[0,\inf)`.
 368 It scales the position of the point from which the Coulombic force will be linearized.
 369 Even though in theory :math:`\alpha_{Q}` can be set to an arbitrarily large value,
 370 algorithmically the linearization point for the force is bound in the range :math:`[0,F_{rcoul}^{Q})`,
 371 where setting :math:`\alpha_{Q}=0` results in a standard hard-core Coulombic interaction.
 372 Setting :math:`\alpha_{Q}` to a larger value softens the Coulombic force.
 373
 374 In all the notations below, for simplicity, the distance between two atoms :math:`i` and :math:`j` is noted as :math:`r`, i.e. :math:`r=r_{ij}`.
 375
 376 Forces: van der Waals interactions
 377 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 378
 379 .. math:: \begin{aligned}
 380           \mathbf{F}_{ij}^{LJ}(\mathbf{r})=\begin{cases}
 381           (\frac{12C_{ij}^{(12)}}{r^{13}} - \frac{6C_{ij}^{(6)}}{r^7})\frac{\mathbf{r}}{r}, & \mbox{if } \mbox{ $r \geq r_{scLJ}$}
 382           \\
 383           \frac{d\mathbf{F}_{ij}^{LJ}}{dr}_{r=r_{scLJ}}r + \mathbf{F}_{ij}^{LJ}(r_{scLJ}), & \mbox{if } \mbox{ $r<r_{scLJ}$}
 384           \end{cases}\end{aligned}
 385           :label: eqvdwforces
 386
 387 where the switching point between the soft and hard-core Lennard-Jones forces
 388 :math:`r_{scLJ} = \alpha_{LJ}(\frac{26}{7}\sigma^6\lambda)^{\frac{1}{6}}` for state A, and
 389 :math:`r_{scLJ} = \alpha_{LJ}(\frac{26}{7}\sigma^6(1-\lambda))^{\frac{1}{6}}` for state B.
 390 In analogy to the Beutler *et al.* soft core version, :math:`\sigma` is the radius of the interaction, which is :math:`(C_{12}/C_6)^{1/6}`
 391 or an input parameter (set with ``sc-sigma-LJ-gapsys``) when C6 or C12 is zero. The default value for this parameter is ``sc-sigma-LJ-gapsys=0.3``.
 392
 393 Explicit expression:
 394
 395 .. math:: \begin{aligned}
 396           \mathbf{F}_{LJ}(\mathbf{r})=\begin{cases}
 397           \left(\frac{12C^{(12)}}{r^{13}} - \frac{6C^{(6)}}{r^7}\right)\frac{\mathbf{r}}{r}, & \mbox{if } \mbox{ $r \geq r_{scLJ}$}
 398           \\
 399           \left(-\frac{156C^{(12)}}{r_{scLJ}^{14}} + \frac{42C^{(6)}}{r_{scLJ}^{8}}\right)\mathbf{r} + \frac{168C^{(12)}}{r_{scLJ}^{13}} - \frac{48C^{(6)}}{r_{scLJ}^{7}}, & \mbox{if } \mbox{ $r<r_{scLJ}$}
 400           \end{cases}\end{aligned}
 401           :label: eqvdwforcesexpl
 402
 403 Forces: Coulomb interactions
 404 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 405
 406 .. math:: \begin{aligned}
 407           \mathbf{F}_{ij}^{Q}(\mathbf{r})=\begin{cases}
 408           \frac{q_{i}q_{j}}{4{\pi}{\varepsilon_0}{\varepsilon_r}r^{2}}\frac{\mathbf{r}}{r}, & \mbox{if } \mbox{ $r \geq r_{scQ} < r_{cutoffQ}$}
 409           \\
 410           \frac{d\mathbf{F}_{ij}^{Q}}{dr}_{r=r_{scQ}}r + \mathbf{F}_{ij}^{Q}(r_{scQ}), & \mbox{if } \mbox{ $r<r_{scQ} < r_{cutoffQ}$}
 411           \\
 412           \frac{d\mathbf{F}_{ij}^{Q}}{dr}_{r=r_{cutoffQ}}r + \mathbf{F}_{ij}^{Q}(r_{cutoffQ}), & \mbox{if } \mbox{ $r < r^{scQ} \geq r_{cutoffQ}$}
 413           \end{cases}\end{aligned}
 414           :label: eqqforces
 415
 416 where the switching point :math:`r^{sc}` between the soft and hard-core electrostatic forces is
 417 :math:`r_{scQ} = \alpha_Q(1+|q_iq_j|)\lambda^{\frac{1}{6}}` for state A, and
 418 :math:`r_{scQ} = \alpha_Q(1+|q_iq_j|)(1-\lambda)^{\frac{1}{6}}` for state B.
 419 The :math:`\lambda` dependence of the linearization point for both van der Waals and Coulombic interactions is of the same power :math:`1/6`.
 420
 421 Explicit expression:
 422
 423 .. math:: \begin{aligned}
 424           \mathbf{F}_{Q}(\mathbf{r})=\begin{cases}
 425           \frac{q_iq_j}{4{\pi}{\varepsilon_0}{\varepsilon_r}r^{2}}\frac{\mathbf{r}}{r}, & \mbox{if } \mbox{ $r \geq r_{scQ} < r_{cutoffQ}$}
 426           \\
 427           \frac{1}{4{\pi}{\varepsilon_0}{\varepsilon_r}}\big( -\frac{2q_{i}q_{j}}{r_{sc}^3}\mathbf{r} + \frac{3q_iq_j}{r_{sc}^2} \big), & \mbox{if } \mbox{ $r<r_{scQ} < r_{cutoffQ}$}
 428           \\
 429           \frac{1}{4{\pi}{\varepsilon_0}{\varepsilon_r}}\big( -\frac{2q_{i}q_{j}}{r_{cutoffQ}^3}\mathbf{r} + \frac{3q_iq_j}{r_{cutoffQ}^2} \big), & \mbox{if } \mbox{ $r < r_{scQ} \geq r_{cutoffQ}$}                        \end{cases}\end{aligned}
 430           :label: eqqforcesexpl
 431
 432 Energies: van der Waals interactions
 433 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 434
 435 Explicition definition of energies:
 436
 437 .. math:: \begin{aligned}
 438           V_{LJ}(r)=\begin{cases}
 439           \frac{C^{(12)}}{r^{12}} - \frac{C^{(6)}}{r^6}, & \mbox{if } \mbox{ $r \geq r_{scLJ}$}
 440           \\
 441           \left(\frac{78C^{(12)}}{r_{scLJ}^{14}} - \frac{21C^{(6)}}{r_{scLJ}^{8}}\right)r^2 - \left(\frac{168C^{(12)}}{r_{scLJ}^{13}} - \frac{48C^{(12)}}{r_{scLJ}^{7}}\right)r
 442           + \frac{91C^{(12)}}{r_{scLJ}^{12}} - \frac{28C^{(6)}}{r_{scLJ}^{6}}, & \mbox{if } \mbox{ $r<r_{scLJ}$}
 443           \end{cases}\end{aligned}
 444           :label: eqvdwener
 445
 446 Energies: Coulomb interactions
 447 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 448
 449 .. math:: \begin{aligned}
 450           V_{Q}(r)=\begin{cases}
 451           \frac{q_{i}q_{j}}{4{\pi}{\varepsilon_0}{\varepsilon_r}r}, & \mbox{if } \mbox{ $r \geq r_{scQ} < r_{cutoffQ}$}
 452           \\
 453           \frac{q_{i}q_{j}}{r_{scQ}^3}r^2 - \frac{3q_iq_j}{r_{scQ}^2}r + \frac{3q_iq_j}{r_{scQ}}, & \mbox{if } \mbox{ $r<r_{scQ} < r_{cutoffQ}$}
 454           \\
 455           \frac{q_{i}q_{j}}{r_{cutoffQ}^3}r^2 - \frac{3q_iq_j}{r_{cutoffQ}^2}r + \frac{3q_iq_j}{r_{cutoffQ}}, & \mbox{if } \mbox{ $r < r_{scQ} \geq r_{cutoffQ}$}
 456           \end{cases}\end{aligned}
 457           :label: eqqener
 458
 459 :math:`\partial H / \partial \lambda`: van der Waals interactions
 460 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 461
 462 Here we provide the explicit expressions of :math:`\partial H/ \partial \lambda` for Lennard-Jones potential, when :math:`r<r_{scLJ}`.
 463 For simplicity, in the expression below we use the notation :math:`r_{scLJ_A}=r_{scA}` and :math:`r_{scLJ_B}=r_{scB}`.
 464
 465 .. math:: \begin{aligned}
 466           \frac{\partial{H}}{\partial{\lambda}} &= V_{LJ}^B(r) - V_{LJ}^A(r) + (1-\lambda)\frac{\partial{V_{LJ}^A(r)}}{\partial{\lambda}} + \lambda\frac{\partial{V_{LJ}^B(r)}}{\partial{\lambda}} \\
 467           & =  \left(\frac{78C^{(12)}_B}{r_{scB}^{14}} - \frac{21C^{(6)}_B}{r_{scB}^{8}}\right)r^2 - \left(\frac{168C^{(12)}_B}{r_{scB}^{13}} - \frac{48C^{(12)}_B}{r_{scB}^{7}}\right)r
 468           + \frac{91C^{(12)}_B}{r_{scB}^{12}} - \frac{28C^{(6)}_B}{r_{scB}^{6}} \\
 469           & -  \left[\left(\frac{78C^{(12)}_A}{r_{scA}^{14}} - \frac{21C^{(6)}_A}{r_{scA}^{8}}\right)r^2 - \left(\frac{168C^{(12)}_A}{r_{scA}^{13}} - \frac{48C^{(12)}_A}{r_{scA}^{7}}\right)r
 470           + \frac{91C^{(12)}_A}{r_{scA}^{12}} - \frac{28C^{(6)}_A}{r_{scA}^{6}} \right]\\
 471           & +  \frac{14(\lambda-1)}{\lambda}\left[\left(\frac{13C^{(12)}_A}{r_{scA}^{14}} - \frac{2C^{(6)}_A}{r_{scA}^{8}}\right)r^2
 472           - \left(\frac{26C^{(12)}_A}{r_{scA}^{13}} - \frac{4C^{(6)}_A}{r_{scA}^{7}}\right)r
 473           + \frac{13C^{(12)}_A}{r_{scA}^{12}} - \frac{2C^{(6)}_A}{r_{scA}^{6}}\right] \\
 474           & +  \frac{14\lambda}{1-\lambda}\left[\left(\frac{13C^{(12)}_B}{r_{scB}^{14}} - \frac{2C^{(6)}_B}{r_{scB}^{8}}\right)r^2
 475           - \left(\frac{26C^{(12)}_B}{r_{scB}^{13}} - \frac{4C^{(6)}_B}{r_{scB}^{7}}\right)r
 476           + \frac{13C^{(12)}_B}{r_{scB}^{12}} - \frac{2C^{(6)}_B}{r_{scB}^{6}}\right] \end{aligned}
 477           :label: eqvdwdhdl
 478
 479 :math:`\partial H/ \partial \lambda` for Lennard-Jones potential, when :math:`r \geq r_{scLJ}` is calculated as a standard hard-core
 480 contribution to :math:`\partial H/ \partial \lambda`: :math:`\frac{\partial{H}}{\partial{\lambda}} = V_{LJ}^B(r) - V_{LJ}^A(r)`.
 481
 482 :math:`\partial H/ \partial \lambda` for Coulomb interactions
 483 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 484
 485 Here we provide the explicit expressions of :math:`\partial H/ \partial \lambda` for Coulomb potential, when :math:`r<r_{scQ}<r_{cutoffQ}`.
 486 For simplicity, in the expression below we use the notation :math:`r_{scQ_A}=r_{scA}` and :math:`r_{scQ_B}=r_{scB}`.
 487
 488 .. math:: \begin{aligned}
 489           \frac{\partial{H}}{\partial{\lambda}} &= V_Q^B(r) - V_Q^A(r) + (1-\lambda)\frac{\partial{V_Q^A(r)}}{\partial{\lambda}} + \lambda\frac{\partial{V_Q^B(r)}}{\partial{\lambda}} \\
 490           & =  \frac{q_{i}^Bq_{j}^B}{r_{scB}^3}r^2 - \frac{3q_i^Bq_j^B}{r_{scB}^2}r + \frac{3q^B_iq_j^B}{r_{scB}} \\
 491           & -  \left[\frac{q_{i}^Aq_{j}^A}{r_{scA}^3}r^2 - \frac{3q_i^Aq_j^A}{r_{scA}^2}r + \frac{3q^A_iq_j^A}{r_{scA}}\right] \\
 492           & +  \frac{\lambda-1}{2\lambda}\left[\frac{q_i^Aq_j^A}{r_{scA}^3}r^2 - \frac{2q_i^Aq_j^A}{r_{scA}^2}r + \frac{q_i^Aq_j^A}{r_{scA}}\right] \\
 493           & +  \frac{\lambda}{2(1-\lambda)}\left[\frac{q_i^Bq_j^B}{r_{scB}^3}r^2 - \frac{2q_i^Bq_j^B}{r_{scB}^2}r + \frac{q_i^Bq_j^B}{r_{scB}}\right] \end{aligned}
 494           :label: eqqdhdl
 495
 496 :math:`\partial H/ \partial \lambda` for Coulomb potential, when :math:`r < r_{scQ} \geq r_{cutoffQ}` is calculated using the same expression above
 497 by setting :math:`r_{scA}=r_{cutoffQ}` and :math:`r_{scB}=r_{cutoffQ}`.
 498
 499 :math:`\partial H/ \partial \lambda` for Coulomb potential, when :math:`r \geq r_{scQ} < r_{cutoffQ}` is calculated as a standard hard-core
 500 contribution to :math:`\partial H/ \partial \lambda`: :math:`\frac{\partial{H}}{\partial{\lambda}} = V_{Q}^B(r) - V_{Q}^A(r)`.
 501
 502
 503