[gmx-developers] Implementing a new type of VdW interaction

[Mark Abraham]

On 5/05/2011 8:26 PM, Lee-Ping Wang wrote:
> Dear Mark,
>
> Thanks for your help.  Last night, I did an early implementation that
> doesn't use the tables.  I put the combination rule into "gmx_mtop_t" as a
> temporary solution.
>
> I do need a full implementation as the testing of my method will include an
> automatic force-matching parameterization, which tunes the parameters
> automatically (see my article on the method,
> http://dx.doi.org/10.1063/1.3519043).
>
> I think I understand your explanation, but I can't find the combination rule
> stored in idef.  The interaction type is indeed passed, but I need more than
> that.

You need "more than that" only if you want more than one VDW interaction 
functional form to be able to be used in the simulation system (e.g. 
normal Buckingham and your form).  And even if you do have multiple VDW 
functional forms (and it will not be easy to use any existing machinery 
to do this), it's already the job of grompp to do all the combining for 
all possible atom-type pairings, so the path of least resistance should 
be to make grompp do that for both functional forms, and mdrun only 
needs machinery to match the pre-combined parameters to the functional form.

If you've only got your modified form, then what's the problem with 
letting grompp do the combining? With or without a new VDW ifunc type?

>    I wish to pass the combination rule in addition to the interaction
> type, because it would be very convenient to use the existing "BHAM" type
> rather than defining an entirely new one (Remark: I have implemented a new
> interaction type for fluctuating charge force fields; please let me know if
> you're interested in adding this to the main version.)

Defining a new interaction type is easy. Crafting combination rule 
machinery into mdrun when it's already found in grompp is hard. See the 
files I pointed out.

> The reason is because my proposed interaction is Buckingham-like; it will
> never be used at the same time as the normal Buckingham potential.  The
> functional form is similar to Buckingham and the number of parameters is the
> same, but there is no singularity at r=0.  This should allow us to model
> somewhat softer repulsions without worrying about the singularity in the
> normal Buckingham potential.  Also, the parameters are given using sigma,
> epsilon, and gamma (a measure of the repulsive strength), which is a little
> more intuitive than the combination of A, b, and c.
>
> Since the functional form is different, I need to add a new case statement
> in the nb_generic.c kernel, and I don't know how to do that without (1)
> defining a new interaction type, which seems quite cumbersome or (2) passing
> the combination rule from grompp into the forcerecord.

By the time code gets into the kernels, the neighbour lists have been 
constructed pre-filtered for interaction type. See init_neighbor_list in 
src/mdlib/ns.c. You either have to hijack the existing Buckingham ivdw 
value, or craft a new one. Either way, I see no value in combining 
parameters in mdrun.

> The head of my parameter file looks like this:
>
> [ defaults ]
> 2    2    no        0.5  0.5
>
> [ atomtypes ]
> H    1    1.0080    0.0000    A    2.655e-01  4.602e-01  1.000e+01
> O    8   15.9950    0.0000    A    3.405e-01  5.648e-02  9.000e+00
>
> Note the second "two" in the defaults section, I want to use this as the
> trigger for the new code.  Also note the order of magnitude of the
> parameters, they are different than for the normal Buckingham interaction.

OK, so I infer you plan to have only this VDW interaction type available.

> The relevant part of gmxdump looks like this:
>
>     ffparams:
>        atnr=2
>        ntypes=6
>           functype[0]=BHAM, a= 3.40500000e-01, b= 5.64800000e-02, c=
> 9.00000000e+00
>           functype[1]=BHAM, a= 3.03000000e-01, b= 1.61220644e-01, c=
> 9.50000000e+00
>           functype[2]=BHAM, a= 3.03000000e-01, b= 1.61220644e-01, c=
> 9.50000000e+00
>           functype[3]=BHAM, a= 2.65500000e-01, b= 4.60200000e-01, c=
> 1.00000000e+01
>           functype[4]=MORSE, b0= 9.71554431e-02, cb= 2.76414222e+02, beta=
> 2.85560502e+01
>           functype[5]=UREY_BRADLEY, theta= 1.15666551e+02, ktheta=
> 3.64462646e+02, r13= 6.69829626e-02, kUB= 9.21828916e+03
>        fudgeQQ              = 0.5
>
> I don't see the combination rule.  Am I looking in the wrong place?

You're looking for the wrong thing. I said "grompp records [the 
parameters] together with their type". So above you have BHAM type 
(which later provokes the right ivdw type) with pre-combined parameters 
for O-O, O-H, H-O and H-H interactions (deduced from the values).

I'm not sure why the hetero interactions are duplicated.

Mark

> Thanks!
>
> - Lee-Ping
>
> From: gmx-developers-bounces at gromacs.org
> [mailto:gmx-developers-bounces at gromacs.org] On Behalf Of Mark Abraham
> Sent: Thursday, May 05, 2011 2:02 AM
> To: Discussion list for GROMACS development
> Subject: Re: [gmx-developers] Implementing a new type of VdW interaction
>
>
>
> On 05/05/11, Lee-Ping<leeping at MIT.EDU>  wrote:
> Dear all,
>
> I would like to implement a new type of vdW interaction.  It's similar
> to the Buckingham interaction, so I was thinking of adding the
> interaction in nb_generic.c with ivdw set to 4. 
>
> First, test the idea with the nonbonded interaction tables - see manual and
> webpage. Don't get hands dirty unnecessarily
>
>
>
> ivdw will be set to four when the Buckingham interaction is turned on
> with the combination rule set to 2 (currently this is unused).
> Unfortunately, it seems like the integer corresponding to the
> combination rule (comb) is used only in grompp, and when we get to mdrun
> it is lost. 
>
> For LJ parameters, I understand why this makes sense; because grompp
> does all the parameter processing, and the nonbonded kernels always see
> c6 and c12 and then apply the same formula; they don't care about the
> user input because grompp takes care of all that.
>
> Right, but grompp records them together with their type. Check out gmxdump
> on a .tpr
> and compare with stuff in src/gmxlib/ifunc.c
>
>
>
> Since my interaction function actually uses a different functional form,
> I can't simply combine the parameters differently in grompp and pass
> them into the Buckingham kernel.  I need a switch in the forcerecord
> that decides between whether to call the normal Buckingham interaction
> or my custom one.
>
> Use of Buckingham is triggered by mdrun noticing the types recorded in the
> .tpr, so you can do a similar thing.
>
> Mark
>
>
> Is my analysis correct?  How do I pass the combination rule into "mdrun"
> so I can trigger my custom VdW interaction?
>
> (My current idea is to add the variable into the main gmx_mtop_t
> topology.  Please let me know if there's a better idea.)
>
> Thank you,
>
> - Lee-Ping
>