[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
801 National Road West
Richmond, IN 47374-4095
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20140703/a80153cb/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list