%matplotlib inline
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
(Note: The Gromacs Manual gives an excellent description of all of this background material)
The GROMACS force field is a typical molecular dynamics force field, including bonds, angles, torsions and non-bonded interactions. There are two types of non-bonded interactions: electrostatic and van der Waals. Electrostatics are modeled with a Coulomb potential:
$V_c(r_{ij}) = f\frac{q_iq_j}{\epsilon_rr_{ij}}$ with $f = \frac{1}{4\pi\epsilon_0} = 138.935485$
and resulting force
$\vec{F}_j(\vec{r}_{ij}) = f\frac{q_iq_j}{\epsilon_rr^2_{ij}}\frac{\vec{r}_{ij}}{r_{ij}} $
van der Waals forces are modeled with a 6-12 Leonard-Jones potential
$ V_{LJ}(r_{ij}) = \frac{C^{(12)}_{ij}} {r^{12}_{ij}} - \frac{C^{(6)}_{ij}} {r^{6}_{ij}}$ where the superscripted numbers in parenthesis denote the name of the constant (i.e. "C twelve")
and resulting force
$ \vec{F}_{i}(\vec{r}_{ij}) = \left(12\frac{C^{(12)}_{ij}} {r^{13}_{ij}} - 6\frac{C^{(6)}_{ij}} {r^{7}_{ij}}\right)\frac{\vec{r}_{ij}}{r_{ij}}$
In order to avoid calculating nonbonded interactions between all $N^2$ possible atom pairs, we want some way to say that, beyond a certain distance, the interaction forces are close enough to zero that we'll just call them zero. The obvious first choice is to pick a cutoff distance and set the potentials to zero at (and after) that distance. The obvious first problem with this is that the derivative of the potential (i.e. the force) would be infinite at the cutoff distance. You get all sorts of nasty effects, including heating. So, instead, you want to smoothly transition the potential to zero. While you're at it, you also want to smoothly transition the force to zero. GROMACS and CHARMM differ in what exactly to call the transition functions, but I'll stick with GROMACS nomenclature here and call this a "shift" function. Section 4 of the GROMACS manual works this out. In the end, we get the following (equations 4.27-4.30):
For a pure Coulomb or Leonard-Jones interaction $F(r) = F_\alpha(r) = r^{-(\alpha + 1)}$ we can write the shifted force and shifted potential as
$ \begin{align*} F_s(r) & = & F_\alpha(r) \quad&r < r_1 \\ F_s(r) & = & \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3 \quad&r_1 \le r < r_c \\ F_s(r) & = & 0 \quad&rc \le r \end{align*} $
and (only middle one is in the manual)
$ \begin{align*} \Phi_s(r) & = & \Phi_\alpha(r) - C \quad&r < r_1 \\ \Phi_s(r) & = & \frac{1}{r^{\alpha}} - \frac{A}{3}(r-r_1)^3 + \frac{B}{4}(r-r_1)^4 - C \quad&r_1 \le r < r_c \\ \Phi_s(r) & = & 0 \quad&rc \le r \end{align*} $
with
$ \begin{align*} A & = & -\frac{(\alpha+4)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^2} \\ B & = & \frac{(\alpha+3)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^3} \\ C & = & \frac{1}{r^{\alpha}_c} - \frac{A}{3}(r_c-r_1)^3 - \frac{B}{4}(r_c-r_1)^4 \end{align*} $
You can see that $C$ is chosen to zero out the potential at $r_c$.
We have gently shifted the force to zero between a start value $r_1$ and a cutoff value $r_c$. This whole thing is also known as a "force shift" since we've applied a shift function
$ S(r) = A(r-r_1)^2 + B(r-r_1)^3 $
to the force in between $r_1$ and $r_c$
The above formulas for $A$ and $B$ come from all of the versions of the GROMACS manual that I checked, as well as [JCTC2006]. I didn't find this until after doing the poking below, but the correct versions of the formulas are in [PCCP2010]
[JCTC2006]: David van der Spoel and Paul J. van Maaren, "The Origin of Layer Structure Artifacts in Simulations of Liquid Water", J. Chem Theory Comput. 2006, 2, 1-11
[PCCP2010]: Siewert J Marrink, Xavier Periole, D. Peter Tieleman and Alex H de Vries, "Comment on 'On using a too large integration time step in molecular dynamics simulations of coarse-grained molecular models' by M. Wigner, D. Trzesniak, R. Barton and W. F. van Gunsteren, Phys. Chem. Chem. Phys., 2009, 11, 1934", Phys. Chem. Chem. Phys, 2010, 12, 2254-2256
We can implement that quickly in Python and see what the results look like. Let's try to reproduce Figure 4.4 from the GROMACS manual, which shows the Coulomb force, Shifted Force, and Shift Function.
We can compute the force two ways: we can use the explicit formula above, or we can take the numerical derivative of the potential. We'll do both.
def get_abc(alpha,r1,rc):
""" NOTE: THIS IS WRONG, AS EXPLAINED LATER!"""
A = -((alpha+4)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**2)
B = ((alpha+3)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**3)
C = (1/rc**alpha) - (A/3)*(rc-r1)**3 - (B/4)*(rc-r1)**4
return A,B,C
def get_s(r,alpha,r1,rc):
A,B,C = get_abc(alpha,r1,rc)
S = A*(r-r1)**2 + B*(r-r1)**3
return S
def get_switched_force(r,alpha,r1,rc):
# Here, alpha should be 1 for electrostatics.
# The definition is F(r) = F_\alpha(r) = r**(-(alpha+1))
unswitched = alpha/r**(alpha+1)
switched = unswitched + get_s(r,alpha,r1,rc)
switched[r<r1] = unswitched[r<r1]
switched[rc<=r] = 0.0
return switched
def get_switched_potential(r,alpha,r1,rc):
unswitched = 1/r**alpha
A,B,C = get_abc(alpha,r1,rc)
result = (1/r**alpha) - (A/3)*(r-r1)**3 - (B/4)*(r-r1)**4 - C
result[r>rc] = 0.0
result[r<r1] = unswitched[r<r1] - C
return result
# And reproducing Fig 4.4
rmin,rmax,rstep=0.8,4.2,0.01
r = np.arange(rmin,rmax,rstep)
r1,rc = 2.0,4.0
analytic_noswitch = 1/r**2
analytic_switch = get_switched_force(r=r,alpha=1.0,r1=r1,rc=rc)
computed_switch = -np.gradient(get_switched_potential(r=r,alpha=1.0,r1=r1,rc=rc),rstep)
fig = plt.figure()
ax = plt.subplot(1,1,1)
plt.plot(r,analytic_noswitch,'k-',label='Normal Force')
plt.plot(r,analytic_switch,'g-',linewidth=1,label='Shifted Force')
plt.plot(r,computed_switch,'r--',linewidth=2,label='Shifted Force from Potential')
plt.plot(r,analytic_switch-analytic_noswitch,'k--',label='Shift function')
ax.set_ylabel('f(r)')
ax.set_xlabel('r')
leg = ax.legend(fancybox=True,numpoints=1,loc='best')
leg.get_frame().set_alpha(0.5)
plt.grid(True)
ax.set_ylim(-0.5,1.5)
ax.set_xlim(0,rmax+1)
plt.title('Comparing Coulomb Force, $r_1={r1}$, $r_c={rc}$'.format(r1=r1,rc=rc))
<matplotlib.text.Text at 0x7fc263b1a0d0>
They picked their limits well; the shift function (especially the shift function you'd get when computing it directly from the potential) can do some crazy things at lower values of r. Those kinds of numerical errors are handled by the MD program, and are a topic for another time.
The electrostatics look good. What about the van der Waals?
Here, I'll make two plots, one near the minimum, and one near $r_c$. This came up in the context of MARTINI for me, so I'll use a MARTINI-related example, two NC3 particles. Those are type Q0, which means the LJ parameters are $C^{(6)} = 0.15091$, $C^{(12)} = 0.16267\times 10^{-2}$. I also ran this example through GROMACS, and stored the results in the file gromacs-results-shift-0.15091-0.0016267.dat
, which I load up below for comparison. I've worked with CHARMM long enough that I tend to keep measurements in Angstroms instead of nanometers, so you'll see a couple of factors of 10 to convert things below.
canned_grom_s = np.array([[ 1.00000000e+00, 1.62655066e+09, 7.97663190e+01], [ 1.10000000e+00, 5.18231616e+08, 7.13482740e+01], [ 1.20000000e+00, 1.82394544e+08, 6.43340610e+01], [ 1.30000000e+00, 6.97899280e+07, 5.83998490e+01], [ 1.40000000e+00, 2.86725980e+07, 5.33143350e+01], [ 1.50000000e+00, 1.25242920e+07, 4.89079400e+01], [ 1.60000000e+00, 5.77020450e+06, 4.50534320e+01], [ 1.70000000e+00, 2.78578025e+06, 4.16535640e+01], [ 1.80000000e+00, 1.40173312e+06, 3.86326940e+01], [ 1.90000000e+00, 7.31753312e+05, 3.59311030e+01], [ 2.00000000e+00, 3.94786031e+05, 3.35010300e+01], [ 2.10000000e+00, 2.19384953e+05, 3.13038060e+01], [ 2.20000000e+00, 1.25211258e+05, 2.93078100e+01], [ 2.30000000e+00, 7.32095391e+04, 2.74869120e+01], [ 2.40000000e+00, 4.37526797e+04, 2.58193510e+01], [ 2.50000000e+00, 2.66734902e+04, 2.42868460e+01], [ 2.60000000e+00, 1.65577715e+04, 2.28739190e+01], [ 2.70000000e+00, 1.04484170e+04, 2.15674110e+01], [ 2.80000000e+00, 6.69197510e+03, 2.03560180e+01], [ 2.90000000e+00, 4.34399707e+03, 1.92300300e+01], [ 3.00000000e+00, 2.85401245e+03, 1.81810050e+01], [ 3.10000000e+00, 1.89528650e+03, 1.72016090e+01], [ 3.20000000e+00, 1.27049036e+03, 1.62854080e+01], [ 3.30000000e+00, 8.58550842e+02, 1.54267680e+01], [ 3.40000000e+00, 5.84058655e+02, 1.46207150e+01], [ 3.50000000e+00, 3.99388000e+02, 1.38628320e+01], [ 3.60000000e+00, 2.74075012e+02, 1.31492150e+01], [ 3.70000000e+00, 1.88389023e+02, 1.24763520e+01], [ 3.80000000e+00, 1.29412094e+02, 1.18411350e+01], [ 3.90000000e+00, 8.85925290e+01, 1.12407570e+01], [ 4.00000000e+00, 6.02145960e+01, 1.06726840e+01], [ 4.10000000e+00, 4.04234580e+01, 1.01346480e+01], [ 4.20000000e+00, 2.65962160e+01, 9.62458600e+00], [ 4.30000000e+00, 1.69344880e+01, 9.14063300e+00], [ 4.40000000e+00, 1.01959190e+01, 8.68108300e+00], [ 4.50000000e+00, 5.51683200e+00, 8.24439600e+00], [ 4.60000000e+00, 2.29282800e+00, 7.82914600e+00], [ 4.70000000e+00, 9.89690000e-02, 7.43405300e+00], [ 4.80000000e+00, -1.36525900e+00, 7.05792000e+00], [ 4.90000000e+00, -2.31307700e+00, 6.69966200e+00], [ 5.00000000e+00, -2.89638000e+00, 6.35827300e+00], [ 5.10000000e+00, -3.22364500e+00, 6.03282300e+00], [ 5.20000000e+00, -3.37248500e+00, 5.72246100e+00], [ 5.30000000e+00, -3.39849400e+00, 5.42639400e+00], [ 5.40000000e+00, -3.34148000e+00, 5.14388800e+00], [ 5.50000000e+00, -3.22990500e+00, 4.87425900e+00], [ 5.60000000e+00, -3.08405200e+00, 4.61687000e+00], [ 5.70000000e+00, -2.91831100e+00, 4.37113900e+00], [ 5.80000000e+00, -2.74279300e+00, 4.13650600e+00], [ 5.90000000e+00, -2.56453300e+00, 3.91245800e+00], [ 6.00000000e+00, -2.38833400e+00, 3.69851400e+00], [ 6.10000000e+00, -2.21739400e+00, 3.49422000e+00], [ 6.20000000e+00, -2.05375600e+00, 3.29915300e+00], [ 6.30000000e+00, -1.89863500e+00, 3.11291100e+00], [ 6.40000000e+00, -1.75266700e+00, 2.93512100e+00], [ 6.50000000e+00, -1.61607600e+00, 2.76543000e+00], [ 6.60000000e+00, -1.48879800e+00, 2.60350000e+00], [ 6.70000000e+00, -1.37058900e+00, 2.44901900e+00], [ 6.80000000e+00, -1.26107300e+00, 2.30168600e+00], [ 6.90000000e+00, -1.15980400e+00, 2.16121700e+00], [ 7.00000000e+00, -1.06629300e+00, 2.02734300e+00], [ 7.10000000e+00, -9.80035000e-01, 1.89981000e+00], [ 7.20000000e+00, -9.00527000e-01, 1.77837500e+00], [ 7.30000000e+00, -8.27273000e-01, 1.66280100e+00], [ 7.40000000e+00, -7.59801000e-01, 1.55287200e+00], [ 7.50000000e+00, -6.97663000e-01, 1.44837600e+00], [ 7.60000000e+00, -6.40432000e-01, 1.34910700e+00], [ 7.70000000e+00, -5.87716000e-01, 1.25487400e+00], [ 7.80000000e+00, -5.39147000e-01, 1.16549200e+00], [ 7.90000000e+00, -4.94382000e-01, 1.08077900e+00], [ 8.00000000e+00, -4.53109000e-01, 1.00056500e+00], [ 8.10000000e+00, -4.15040000e-01, 9.24684000e-01], [ 8.20000000e+00, -3.79907000e-01, 8.52977000e-01], [ 8.30000000e+00, -3.47469000e-01, 7.85290000e-01], [ 8.40000000e+00, -3.17502000e-01, 7.21472000e-01], [ 8.50000000e+00, -2.89802000e-01, 6.61381000e-01], [ 8.60000000e+00, -2.64182000e-01, 6.04877000e-01], [ 8.70000000e+00, -2.40472000e-01, 5.51823000e-01], [ 8.80000000e+00, -2.18516000e-01, 5.02087000e-01], [ 8.90000000e+00, -1.98172000e-01, 4.55542000e-01], [ 9.00000000e+00, -1.79309000e-01, 4.12063000e-01], [ 9.10000000e+00, -1.61813000e-01, 3.71527000e-01], [ 9.20000000e+00, -1.45595000e-01, 3.33814000e-01], [ 9.30000000e+00, -1.30580000e-01, 2.98810000e-01], [ 9.40000000e+00, -1.16698000e-01, 2.66400000e-01], [ 9.50000000e+00, -1.03885000e-01, 2.36473000e-01], [ 9.60000000e+00, -9.20810000e-02, 2.08918000e-01], [ 9.70000000e+00, -8.12330000e-02, 1.83629000e-01], [ 9.80000000e+00, -7.12880000e-02, 1.60499000e-01], [ 9.90000000e+00, -6.22000000e-02, 1.39426000e-01], [ 1.00000000e+01, -5.39230000e-02, 1.20306000e-01], [ 1.01000000e+01, -4.64140000e-02, 1.03039000e-01], [ 1.02000000e+01, -3.96340000e-02, 8.75240000e-02], [ 1.03000000e+01, -3.35420000e-02, 7.36650000e-02], [ 1.04000000e+01, -2.81010000e-02, 6.13620000e-02], [ 1.05000000e+01, -2.32750000e-02, 5.05210000e-02], [ 1.06000000e+01, -1.90260000e-02, 4.10450000e-02], [ 1.07000000e+01, -1.53200000e-02, 3.28410000e-02], [ 1.08000000e+01, -1.21200000e-02, 2.58150000e-02], [ 1.09000000e+01, -9.39300000e-03, 1.98730000e-02], [ 1.10000000e+01, -7.10100000e-03, 1.49230000e-02], [ 1.11000000e+01, -5.21000000e-03, 1.08740000e-02], [ 1.12000000e+01, -3.68400000e-03, 7.63400000e-03], [ 1.13000000e+01, -2.48500000e-03, 5.11300000e-03], [ 1.14000000e+01, -1.57500000e-03, 3.21900000e-03], [ 1.15000000e+01, -9.18000000e-04, 1.86200000e-03], [ 1.16000000e+01, -4.73000000e-04, 9.53000000e-04], [ 1.17000000e+01, -2.01000000e-04, 4.02000000e-04], [ 1.18000000e+01, -6.00000000e-05, 1.19000000e-04], [ 1.19000000e+01, -8.00000000e-06, 1.50000000e-05], [ 1.20000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.21000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.22000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.23000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.24000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.25000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.26000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.27000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.28000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.29000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.30000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.31000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.32000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.33000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.34000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.35000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.36000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.37000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.38000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.39000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.40000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.41000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.42000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.43000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.44000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.45000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.46000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.47000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.48000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.49000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.50000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.51000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.52000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.53000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.54000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.55000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.56000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.57000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.58000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.59000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.60000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.61000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.62000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.63000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.64000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.65000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.66000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.67000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.68000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.69000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.70000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.71000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.72000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.73000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.74000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.75000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.76000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.77000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.78000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.79000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.80000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.81000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.82000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.83000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.84000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.85000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.86000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.87000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.88000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.89000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.90000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.91000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.92000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.93000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.94000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.95000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.96000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.97000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.98000000e+01, 0.00000000e+00, 0.00000000e+00], [ 1.99000000e+01, 0.00000000e+00, 0.00000000e+00]])
print canned_grom_s.shape
(190, 3)
def v_vdw(r,c6,c12):
v_vdw = c12/r**12 - c6/r**6
return v_vdw
def v_vdw_switch(r,r1,rc,c12,c6):
vs12 = get_switched_potential(r,alpha=12,r1=r1,rc=rc)
vs6 = get_switched_potential(r,alpha=6,r1=r1,rc=rc)
v_vdw_switch = c12*vs12 - c6*vs6
return v_vdw_switch
def plot_vdw_potentials():
c6 = 0.15091E-00
c12 = 0.16267E-02
try:
grom_s = np.loadtxt('gromacs-comparisons/dat/gromacs-results-shift-0.15091-0.0016267.dat')
except IOError:
grom_s = canned_grom_s
r = np.arange(1,20,0.1)
analytic_noswitch_potential = v_vdw(r=r*.1,c12=c12,c6=c6)
analytic_switch_potential = v_vdw_switch(r=r*.1,c12=c12,c6=c6,r1=0.9,rc=1.2)
fig = plt.figure()
ax = plt.subplot(111)
plt.plot(r,analytic_noswitch_potential,'k.-',label='Normal potential')
plt.plot(r,analytic_switch_potential,'g--',linewidth=3,label='Our shifted potential')
plt.plot(grom_s[:,0],grom_s[:,1],'r-',label='GROMACS shifted potential')
leg = ax.legend(fancybox=True)
leg.get_frame().set_alpha(0.5)
ax.set_ylim(-4,.1)
ax.set_ylim(-3.5,-3)
ax.set_xlim(5,5.5)
ax.set_ylabel('$V_{LJ}$')
ax.set_xlabel('r (Angstroms)')
plt.grid(True)
fig = plt.figure()
ax = plt.subplot(111)
plt.plot(r,analytic_noswitch_potential,'k-',label='Normal potential')
plt.plot(r,analytic_switch_potential,'g--',linewidth=3,label='Our shifted potential')
plt.plot(grom_s[:,0],grom_s[:,1],'r-',label='GROMACS shifted potential')
leg = ax.legend(fancybox=True)
leg.get_frame().set_alpha(0.5)
ax.set_ylim(-.1,.1)
ax.set_xlim(10,13)
ax.set_ylabel('$V_{LJ}$')
ax.set_xlabel('r (Angstroms)')
plt.grid(True)
plot_vdw_potentials()
There are a few things to notice here
That last part is pretty worrying. Let's see if it's real, or just a visual artifact.
def plot_vdw_forces():
c6 = 0.15091E-00
c12 = 0.16267E-02
try:
grom_s = np.loadtxt('gromacs-comparisons/dat/gromacs-results-shift-0.15091-0.0016267.dat')
except IOError:
grom_s = canned_grom_s
r = np.arange(1,20,0.1)
analytic_switch_potential = v_vdw_switch(r=r*.1,c12=c12,c6=c6,r1=0.9,rc=1.2)
computed_switched_force = -np.gradient(analytic_switch_potential,0.1)
gromacs_computed_switched_force = -np.gradient(grom_s[:,1],grom_s[:,0][1] - grom_s[:,0][0])
fig = plt.figure()
ax = plt.subplot(111)
plt.plot(r,computed_switched_force,'g--',linewidth=3,label='Shifted force from our potential')
plt.plot(grom_s[:,0],gromacs_computed_switched_force,'k-',label='Shifted force from GROMACS potential')
leg = ax.legend(fancybox=True,loc=4)
leg.get_frame().set_alpha(0.5)
ax.set_ylim(-.1,.001)
ax.set_xlim(10,13)
ax.set_ylabel('$F_{LJ}$')
ax.set_xlabel('r (Angstroms)')
plt.grid(True)
plot_vdw_forces()
Eek. So, as we guessed, our forces aren't matching up properly (that diagonal line straddling $r=12.0$ would become vertical if we made the $r$ step size increasingly small; right now, it's set to 0.1). Meanwhile, GROMACS (the code, not the manual) is doing just fine.
So, something is wrong with the formula. Let's take a look at the shifted force in the shifted region:
$ F_s(r) = \frac{\alpha}{r^{\alpha+1}} + A(r-r_1)^2 + B(r-r_1)^3 $
with
$ \begin{align*} A & = & -\frac{(\alpha+4)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^2} \\ B & = & \frac{(\alpha+3)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^3} \\ C & = & \frac{1}{r^{\alpha}_c} - \frac{A}{3}(r_c-r_1)^3 + \frac{B}{4}(r_c-r_1)^4 \end{align*} $
This should be zero at $r=r_c$. However, if we plug $r=r_c$ into $F_s(r)$ and simplify a bit, we get
$ \begin{align*} F_s(r_c) &=& \frac{\alpha}{r_c^{\alpha+1}}& + A(r_c-r_1)^2 + B(r_c-r_1)^3 \\ &=& \frac{\alpha}{r_c^{\alpha+1}}& -\frac{(\alpha+4)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^2}(r_c-r_1)^2 + \frac{(\alpha+3)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^3}(r_c-r_1)^3 \\ &=& \frac{\alpha}{r_c^{\alpha+1}}& -\frac{(\alpha+4)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}} + \frac{(\alpha+3)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}} \\ &=& \frac{\alpha}{r_c^{\alpha+1}}& -\frac{(\alpha+4 - \alpha - 3)r_c - (\alpha+1 - \alpha - 1)r_1}{r_c^{\alpha + 2}} \\ F_s(r_c) &=& \frac{\alpha}{r_c^{\alpha+1}}& - \frac{1}{r_c^{\alpha+1}} \end{align*} $
which explains everything above. If we have $\alpha=1$ (i.e. the electrostatics case), this goes to zero and matches up with what we expect (and with the explicit Coulomb functions given in the manual and papers). If we don't, however, we get a discontinuity in the force. Now that we know what's going on, we see that things will cancel properly if we stick a leading $\alpha$ onto the last two terms in $F_s$. I.e. we replace $A$ with $\alpha A$ and $B$ with $\alpha B$:
$ \begin{align*} A & = & -\alpha\frac{(\alpha+4)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^2} \\ B & = & \alpha\frac{(\alpha+3)r_c - (\alpha+1)r_1}{r_c^{\alpha + 2}(r_c - r_1)^3} \\ C & = & \frac{1}{r^{\alpha}_c} - \frac{A}{3}(r_c-r_1)^3 - \frac{B}{4}(r_c-r_1)^4 \end{align*} $
Let's try it:
def get_abc(alpha,r1,rc):
""" This is now correct """
A = -((alpha+4)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**2)
A = alpha*A
B = ((alpha+3)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**3)
B = alpha*B
C = (1/rc**alpha) - (A/3)*(rc-r1)**3 - (B/4)*(rc-r1)**4
return A,B,C
plot_vdw_potentials()
plot_vdw_forces()
And all was good in the world.