# [gmx-developers] Calculating shifted LJ potentials

Michael Lerner mglerner at gmail.com
Thu Jul 3 23:28:10 CEST 2014

```Hi all,

I was wondering if you could point me at an explicit formula for shifted LJ
potentials in GROMACS. I've emailed gmx-users at gromacs.org a couple of times
with no response, and was hoping that you might be able to help. When I try
to use the formulas in the GROMACS manual to calculate the LJ potential by
hand, I get a different answer than that reported by GROMACS. My guess is
that I'm just using the wrong formula for the shifted potential (I applied
formula 4.29 from the GROMACS manual to the 6 term and the 12 term
individually, and added the results together), so hopefully just pointing
me at an explicit formula will be enough.

I've included a full description below (calculations by hand, a full set of
GROMACS input files and commands, and Python code that calculates the
potential analytically) in case it helps.

Thank you,
-Michael

P.S. I can, of course, send the GROMACS files along as attachments if you
want, but I assume that attachments are poor form on mailing lists.

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.

--
Michael Lerner
Department of Physics and Astronomy
Earlham College - Drawer 111