[gmx-users] Calculating shifted VDW energies

[Michael Lerner]

Hi all,

The manual lists an explicit formula for an electrostatic shift with
r1=0.0. Could someone point me at a similar explicit formula for a VDW
shift including both r1 and rc? If so, I'm sure I could double-check my
numbers myself.

Thanks,
-Michael


On Wed, Jun 25, 2014 at 11:40 AM, Michael Lerner <mglerner at gmail.com> wrote:

> 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
>



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