[gmx-developers] Calculating shifted LJ potentials
Michael Lerner
mglerner at gmail.com
Tue Jul 8 20:40:56 CEST 2014
Hi all,
After poking at this for a while, I decided to graph things. The formulas
in the manual work for alpha=1 (i.e. electrostatics). For higher alpha, the
potential goes to zero at rc, but the force does not. You can see this for
yourself by plugging r=rc into equation (4.28) in the manual. After some
simplification, you end up with
F_s = alpha/rc**(alpha+1) - 1/rc**(alpha+1)
which zeros out for alpha=1, but not for higher alpha. I'm pretty sure
there's a typo in the manual, and that both A and B in (4.27) should have a
leading factor of alpha. When I do that, the force is smooth, and matches
up with the results I get from GROMACS. So, I think there's a typo in the
manual (and in at least one paper), but not a bug in the code.
After figuring that out, I did find the formulas with a leading alpha in at
least one paper (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).
I made some pictures along the way to convince myself I was doing the right
thing. If anyone else is curious, you can see the IPython Notebook at
http://nbviewer.ipython.org/gist/mglerner/750cada28c651c0a3507
Cheers,
-Michael
On Thu, Jul 3, 2014 at 5:27 PM, Michael Lerner <mglerner at gmail.com> wrote:
> 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
>
--
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/20140708/e45c9505/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list