[gmx-users] Calculating shifted VDW energies

[Michael Lerner]

Hi all,

I'm having trouble replicating the VDW energy from a shifted potential.
When I calculate unshifted potentials by hand, I get an answer that agrees
with GROMACS. When I calculate shifted potentials, I don't. I've included
the intermediate numbers I get as well as a simple GROMACS test case, and a
Python function that I was using to calculate the shifted potential. I
would really appreciate you help, as I'm probably making a simple mistake.
Thanks!

Here's what I think should be a simple example, two particles of the same
type, 1.0 nm apart

sigma = 1
epsilon = 0.25

V = 4*epsilon*sigma**6 = 1     [kJ mol^{-1}nm^6]
W = 4*epsilon*sigma**12 = 1  [kJ mol^{-1}nm^{12}]

So, V_LJ should be be W/r**12 - V/r**6 = 1/1 - 1/1 = 0. I get that when I
run a simple simulation.

Meanwhile, I think I'm doing wrong with shift functions. I'm trying to
replicate

vdwtype = Shift
rvdw_switch = 0.9
rvdw = 1.2

If I understand the manual correctly, I want to make repeated use of
equation Phi from 4.29 (and thus also need equations 4.27 and 4.30). Since
V and W are both 1 here, I want

Phi(r=1, alpha=12, r1=0.9, rc=1.2) - Phi(r=1, alpha=6, r1=0.9, rc=1.2)

If I call the first set of constants A12, B12 and C12, and the second A6,
B6, C6, I get

A12 = -( (12 + 4)*1.2 - (12+1)*0.9 ) / (1.2**(12+2) * (1.2-0.9)**2)
= -6.490547151887453
B12 = ((12 + 3)*1.2 - (12+1)*0.9) / (1.2**(12+2) * (1.2-0.9)**3) =
18.17353202528487
C12 = (1/1.2**12) - (A12/3)*(1.2-0.9)**3 - (B12/4)*(1.2-0.9)**4
= 0.13377017680040035
PHI12 = 1/1**12 - (A12/3)*(1-0.9)**3 - (B12/4)*(1-0.9)**4 - C12
= 0.8679390006162634

and

A6 = -( (6 + 4)*1.2 - (6+1)*0.9 ) / (1.2**(6+2) * (1.2-0.9)**2)
= -14.729309159553942
B6 = ((6 + 3)*1.2 - (6+1)*0.9) / (1.2**(6+2) * (1.2-0.9)**3)
= 38.761339893563004
C6 = (1/1.2**6) - (A6/3)*(1.2-0.9)**3 - (B6/4)*(1.2-0.9)**4
= 0.38897004583190453
PHI6 = 1/1**6 - (A6/3)*(1-0.9)**3 - (B6/4)*(1-0.9)**4 - C6
= 0.6149706903906076

So I think the shifted potential here should be PHI12 - PHI6
= 0.2529683102256558

Meanwhile, when I use GROMACS, I get 0.284678

I'm sure I'm making some simple mistake in my calculations, but I can't
find it. Perhaps I've misunderstood how to handle a shift function for VDW
interactions.

In order to test with GROMACS, I use the following commands:

grompp -f simple-energy.mdp -c martini-twoparticle-1.00.gro -p
martini-twoparticle-tworesi.top -o martini-twoparticle-1.00.tpr
mdrun -v -deffnm martini-twoparticle-1.00
echo "1 2" | g_energy -f martini-twoparticle-1.00.edr -s
martini-twoparticle-1.00.tpr -o martini-twoparticle-1.00.xvg

with the following files:

--- begin simple-energy.mdp ---
title                    = Simple Energy
integrator               = md
unconstrained_start      = yes
nsteps                   = 0
nstlist                  = 1
pbc                      = no
rlist                    = 100.0
coulombtype              = Shift
rcoulomb_switch          = 0.0
rcoulomb                 = 1.2
epsilon_r                = 15
vdw_type                 = Shift
rvdw_switch              = 0.9
rvdw                     = 1.2
--- end simple-energy.mdp ---

--- begin martini-two-particle-two-resi.top ---
#include "martini-test-short.itp"

[ system ]
; name
TWO PARTICLES

[ molecules ]
; name number
one 2
--- end martini-two-particle-two-resi.top ---

--- begin martini-test-short.itp ---
; Topology for test particles

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
  1             1               no              1.0     1.0

[ atomtypes ]
;name at.num mass charge ptype V(c6) W(c12)
D     1      72.0 0.000  A     1.0 1.0

[ nonbond_params ]
; i j funda c6 c12
   D              D     1       1.0000000000000000 1.0000000000000000

[ moleculetype ]
; molname nrexcl
one 0

[ atoms]
;id     type   resnr residu atom cgnr    charge
   1    D        1    TWO    D1   1       -1.00

[ bonds ]

--- end martini-test-short.itp ---

--- begin martini-twoparticle-1.00.gro ---
TWO PARTICLES
    2
    1ONE  D1       1   0.000   0.000   0.000
    2ONE  D1       1   1.000   0.000   0.000
 100.00000 100.00000 100.00000

--- end martini-twoparticle-1.00.gro ---

Finally, here's the Python code I was using to try to generate shifted
potentials:

def get_switched(r,alpha,r1,rc):
    unswitched = 1/r**alpha
    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 (1/r) - (5/(3*rc)) + (5*r**3)/(3*rc**4) - (r**4)/(rc**5)
    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]
    return result

def v_vdw_switch(r,r1,rc,c12,c6):
    vs12 = get_switched(r,alpha=12,r1=r1,rc=rc)
    vs6 = get_switched(r,alpha=6,r1=r1,rc=rc)
    v_vdw_switch = c12*vs12 - c6*vs6
    return v_vdw_switch

as a side note, the above function *does* replicate the standard
electrostatics shift with r1=0.0.

Thanks,
-Michael

-- 
Michael Lerner
Department of Physics and Astronomy
Earlham College - Drawer 111
801 National Road West
Richmond, IN   47374-4095