- 2018-04-25 16:24:25
The following is my notes about the GROMACS Soft-Core Potential for Free Energy Calculation
Theory
The soft-core potential between two states A and B is defined as
$\begin{split}
V_{sc}(r) &=(1-\lambda)V_A(r_A) + \lambda V_B(r_B)
&=(1-\lambda)V_A(r) + \lambda V_B(r) \;\;\;\; \text{if } \alpha \sigma=0
r_A &= \left(\alpha \sigma_A^6 \lambda^p+r^6\right)^{1/6}
r_B &= \left(\alpha \sigma_B^6 (1-\lambda)^p+r^6\right)^{1/6}
\end{split}$
When $\alpha=0$ or $\sigma=0$, the soft-core potential degenerates to the original potential. $\alpha$ is usually about 0.5 and $p=1$ or $p=2$. The recommended combination is $\alpha=0.5, p=1$.
For LJ potential,
$\alg
V &= { C_{12} \over r^{12} } - {C_6 \over r^6}
V_{sc} &=C_{12}\left[{1-\lambda \over (\alpha \sigma^6 \lambda^p+r^6)^2} + {\lambda \over (\alpha \sigma^6(1-\lambda)^6+r^6)^2}\right]
&-C_6 \left[{1-\lambda \over \alpha \sigma^6 \lambda^p+r^6} + {\lambda \over \alpha \sigma^6(1-\lambda)^6+r^6}\right]
\ealg$
Some special cases
$\alg
V_{sc}=0 &\Rightarrow V(r_A)=0 \; V(r_B)=0
r_A=r_B &\Rightarrow \lambda=0.5
\ealg$
If A state is totally decoupled, which means $V_A=0$, in this case
$\alg
V_{sc}(r) &= \lambda V_B(r_B)
&= \lambda \left\{ {C_{12} \over \left(\alpha \sigma^6 (1-\lambda)^p+r^6\right)^2} -{C_6 \over \alpha \sigma^6 (1-\lambda)^p+r^6} \right\}
\ealg$
Define
$\alg
s^6 &=\alpha\sigma^6
r_{sc} &= \left(r^6+ (1-\lambda)^p s^6\right)^{1/6}
\min r_{sc} &= r &\;\; \text{when } \lambda=1
\max r_{sc} &= (r^6+s^6)^{1/6} &\;\; \text{when } \lambda=0
\ealg$
Plot of $r_{sc}$ with $\l$
\[r_{sc}/r = \left(1+(1-\lambda)^p (s/r)^6\right)^{1/6}\]
Plot of $\max r_{sc}$ with $\l$
\[\max r_{sc}/r = (1+(s/r)^6)^{1/6}\]
Plot of $V_{sc}(r)$
\[\alg V_{sc}(r, s, \l) &= \lambda \left\\{ {C_{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{C_6 \over r^6+(1-\lambda)^p s^6} \right\\} \\ &=4\e \lambda \left\\{ {\s^{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{\s^6 \over r^6+(1-\lambda)^p s^6} \right\\} \ealg\]
GMX Options
In GMX, $\alpha, \sigma, p$ is corresponding to sc_alpha, sc_sigma, sc_power, respectively.
The manual said
-
sc-alpha: (0) the soft-core alpha parameter, a value of 0 results in linear interpolation of the LJ and Coulomb interactions -
sc-r-power: (6) the power of the radial term in the soft-core equation. Possible values are 6 and 48. 6 is more standard, and is the default. When 48 is used, then sc-alpha should generally be much lower (between 0.001 and 0.003). -
sc-power: (0) the power for lambda in the soft-core function, only the values 1 and 2 are supported -
sc-sigma: (0.3) [nm] the soft-core sigma for particles which have a C6 or C12 parameter equal to zero or a sigma smaller than sc-sigma
The description of sc-sigma is not correct, actually
$\sigma= \begin{cases}
& \sigma^*=(C_{12}/C_{6})^{1/6}
& \text{sc_sigma} \; \text{ if } C_6\times C_{12}=0
\end{cases}$
See the discussion
Testing
Script
Figures
Conclusion
-
GROMACS DOES support table potential for free energy calculations
-
GROMACS calculates the interaction only If
r_scless thanrvdw -
The
r_scfor any distance should be in the range of the table files, otherwise GROMACS will use a unkown value -
If
DispCorris used, be sureC6is the correct value -
If
sc-alphais not zero, be sure the $\s$ GROMACS used is the correct value -
If
DispCorrnot used, put total interaction in dispersion or repulsion column will force GROMACS to usesc-sigma -
If
DispCorris used andsc-alphais not zero, be sureC6andC12are both correct
