<div dir="ltr">Hi all,<div><br></div><div>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</div>
<div><br></div><div>F_s = alpha/rc**(alpha+1) - 1/rc**(alpha+1)</div><div><br></div><div>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.</div>
<div><br></div><div>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).</div>
<div><br></div><div>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 <a href="http://nbviewer.ipython.org/gist/mglerner/750cada28c651c0a3507">http://nbviewer.ipython.org/gist/mglerner/750cada28c651c0a3507</a></div>
<div><br></div><div>Cheers,</div><div>-Michael</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Jul 3, 2014 at 5:27 PM, Michael Lerner <span dir="ltr"><<a href="mailto:mglerner@gmail.com" target="_blank">mglerner@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi all,<div><br></div><div><div>I was wondering if you could point me at an explicit formula for shifted LJ potentials in GROMACS. I've emailed <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a> 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.</div>
<div><br></div><div>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.</div><div>
<br></div><div>Thank you,</div><div>-Michael</div><div><br></div><div>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.</div><div>
<br></div><div><br></div><div><br></div><div><br></div><div>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!</div>
<div><br></div><div>Here's what I think should be a simple example, two particles of the same type, 1.0 nm apart</div><div><br></div><div>sigma = 1</div><div>epsilon = 0.25</div><div><br></div><div>V = 4*epsilon*sigma**6 = 1 [kJ mol^{-1}nm^6]</div>
<div>W = 4*epsilon*sigma**12 = 1 [kJ mol^{-1}nm^{12}]</div><div><br></div><div>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.</div><div><br></div><div>Meanwhile, I think I'm doing wrong with shift functions. I'm trying to replicate</div>
<div><br></div><div>vdwtype = Shift</div><div>rvdw_switch = 0.9</div><div>rvdw = 1.2</div><div><br></div><div>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</div>
<div><br></div><div>Phi(r=1, alpha=12, r1=0.9, rc=1.2) - Phi(r=1, alpha=6, r1=0.9, rc=1.2)</div><div><br></div><div>If I call the first set of constants A12, B12 and C12, and the second A6, B6, C6, I get</div><div><br></div>
<div>A12 = -( (12 + 4)*1.2 - (12+1)*0.9 ) / (1.2**(12+2) * (1.2-0.9)**2) = -6.490547151887453</div><div>B12 = ((12 + 3)*1.2 - (12+1)*0.9) / (1.2**(12+2) * (1.2-0.9)**3) = 18.17353202528487</div><div>C12 = (1/1.2**12) - (A12/3)*(1.2-0.9)**3 - (B12/4)*(1.2-0.9)**4 = 0.13377017680040035</div>
<div>PHI12 = 1/1**12 - (A12/3)*(1-0.9)**3 - (B12/4)*(1-0.9)**4 - C12 = 0.8679390006162634</div><div><br></div><div>and</div><div><br></div><div>A6 = -( (6 + 4)*1.2 - (6+1)*0.9 ) / (1.2**(6+2) * (1.2-0.9)**2) = -14.729309159553942</div>
<div>B6 = ((6 + 3)*1.2 - (6+1)*0.9) / (1.2**(6+2) * (1.2-0.9)**3) = 38.761339893563004</div><div>C6 = (1/1.2**6) - (A6/3)*(1.2-0.9)**3 - (B6/4)*(1.2-0.9)**4 = 0.38897004583190453</div><div>PHI6 = 1/1**6 - (A6/3)*(1-0.9)**3 - (B6/4)*(1-0.9)**4 - C6 = 0.6149706903906076</div>
<div><br></div><div>So I think the shifted potential here should be PHI12 - PHI6 = 0.2529683102256558</div><div><br></div><div>Meanwhile, when I use GROMACS, I get 0.284678</div><div><br></div><div>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. </div>
<div><br></div><div>In order to test with GROMACS, I use the following commands:</div><div><br></div><div>grompp -f simple-energy.mdp -c martini-twoparticle-1.00.gro -p martini-twoparticle-tworesi.top -o martini-twoparticle-1.00.tpr</div>
<div>mdrun -v -deffnm martini-twoparticle-1.00</div><div>echo "1 2" | g_energy -f martini-twoparticle-1.00.edr -s martini-twoparticle-1.00.tpr -o martini-twoparticle-1.00.xvg</div><div><br></div><div>with the following files:</div>
<div><br></div><div>--- begin simple-energy.mdp ---</div><div>title = Simple Energy</div><div>integrator = md</div><div>unconstrained_start = yes</div><div>nsteps = 0</div>
<div>nstlist = 1</div><div>pbc = no</div><div>rlist = 100.0</div><div>coulombtype = Shift </div><div>rcoulomb_switch = 0.0</div><div>rcoulomb = 1.2</div>
<div>epsilon_r = 15</div><div>vdw_type = Shift </div><div>rvdw_switch = 0.9</div><div>rvdw = 1.2</div><div>--- end simple-energy.mdp ---</div><div><br></div>
<div>--- begin martini-two-particle-two-resi.top ---</div><div>#include "martini-test-short.itp"</div><div><br></div><div>[ system ]</div><div>; name</div><div>TWO PARTICLES</div><div><br></div><div>[ molecules ]</div>
<div>; name number</div><div>one 2</div><div>--- end martini-two-particle-two-resi.top ---</div><div><br></div><div>--- begin martini-test-short.itp ---</div><div>; Topology for test particles</div><div><br></div><div>[ defaults ]</div>
<div>; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ</div><div> 1 1 no 1.0 1.0</div><div><br></div><div>[ atomtypes ]</div><div>;name at.num mass charge ptype V(c6) W(c12)</div>
<div>D 1 72.0 0.000 A 1.0 1.0</div><div><br></div><div>[ nonbond_params ]</div><div>; i j funda c6 c12</div><div> D D 1 1.0000000000000000 1.0000000000000000</div><div><br></div><div>
[ moleculetype ]</div><div>; molname nrexcl</div><div>one 0</div><div><br></div><div>[ atoms]</div><div>;id type resnr residu atom cgnr charge</div><div> 1 D 1 TWO D1 1 -1.00</div><div>
<br>
</div><div>[ bonds ]</div><div><br></div><div>--- end martini-test-short.itp ---</div><div><br></div><div>--- begin martini-twoparticle-1.00.gro ---</div><div>TWO PARTICLES</div><div> 2</div><div> 1ONE D1 1 0.000 0.000 0.000</div>
<div> 2ONE D1 1 1.000 0.000 0.000</div><div> 100.00000 100.00000 100.00000</div><div><br></div><div>--- end martini-twoparticle-1.00.gro ---</div><div><br></div><div>Finally, here's the Python code I was using to try to generate shifted potentials:</div>
<div><br></div><div>def get_switched(r,alpha,r1,rc):</div><div> unswitched = 1/r**alpha</div><div> A = -((alpha+4)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**2)</div><div> B = ((alpha+3)*rc-(alpha+1)*r1)/(rc**(alpha+2)*(rc-r1)**3)</div>
<div> C = (1/rc**alpha) - (A/3)*(rc-r1)**3 - (B/4)*(rc-r1)**4</div><div> #return (1/r) - (5/(3*rc)) + (5*r**3)/(3*rc**4) - (r**4)/(rc**5)</div><div> result = (1/r**alpha) - (A/3)*(r-r1)**3 - (B/4)*(r-r1)**4 - C</div>
<div> result[r>rc] = 0.0</div><div> result[r<r1] = unswitched[r<r1]</div><div> return result</div><div><br></div><div>def v_vdw_switch(r,r1,rc,c12,c6):</div><div> vs12 = get_switched(r,alpha=12,r1=r1,rc=rc)</div>
<div> vs6 = get_switched(r,alpha=6,r1=r1,rc=rc)</div><div> v_vdw_switch = c12*vs12 - c6*vs6</div><div> return v_vdw_switch</div><div><br></div><div>as a side note, the above function *does* replicate the standard electrostatics shift with r1=0.0.</div>
<span class="HOEnZb"><font color="#888888">
<div><br></div>-- <br>Michael Lerner<br><div><div>Department of Physics and Astronomy</div><div>Earlham College - Drawer 111</div><div>801 National Road West</div><div>Richmond, IN 47374-4095</div></div>
</font></span></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>Michael Lerner<br><div><div>Department of Physics and Astronomy</div><div>Earlham College - Drawer 111</div><div>801 National Road West</div><div>Richmond, IN 47374-4095</div>
</div>
</div>