[gmx-developers] modifying mdrun for non-pairwise-additive potentials

[Erik Lindahl]

Hi Graham,

Sorry for the delay - I've been home with our daughter for two weeks.

>
>I hoped at first that f_coul, f_lj6 and f_lj12 would be accumulated
>separately for each particle, then summed, so I could get at them there,
>but instead it seems that they are added up in the inner loops. I don't
>
You're right - we have to sum the there for efficiency reasons.

>
>have the godlike programming powers needed to alter these, so I'm thinking
>instead of changing function force in force.c to call do_fnbf twice :
>
>  do_fnbf(...forcerec fr, ... rvec f[], .. 
>  do_fnbf(...forcerec fr0, ... rvec f0[], .. 
>
>where fr0 has all the c12_ab parameters set to 0, so will get zero for
>all the hydrophobic forces. Thus f[i] - f0[i] is the hydrophobic force,
>which I can fix and add back to f0. 
>
>Is there any reason why this won't work? I don't mind the big performance
>hit. 
>
No, it should work just splendid, but with a performance hit.

To improve, you should have a look at the loop type encoded in each 
neighborlist. There is a switch statement in fnbf.c that determines 
which innerloop we actually call (e.g. coulomb+lj or lj only). The best 
way to do things might actually be to edit stuff in fnbf.c instead, and 
send different force arrays (coulomb or LJ) to the inner loops. The case 
where we call a LJ+coulomb inner loop should then simply be replaced 
with two calls to LJ-only and Coulomb-only loops, respectively.

This will be significantly faster than using zero-value parameters.

If you use water in your system you will notice that there is no LJ-only 
water (or water-water) loop. This is because only the first particle 
(OW) has any LJ parameters, so you can simply use the non-water loop if 
you only want the LJ part.

The main advantage with this approach is that you won't have to study 
the details of the forcerec, just add two temporary force arrays and do 
some bookkeeping.

>
>
>Also, could somebody help me with what exactly to I have to do in forcerec
>fr0? I guess should be an almost exact copy of fr, (so I'll copy the ints,
>reals etc, and set internal pointers to point to the corresponding parts
>of fr). The exception is c12 parameters; are these fr->vdwtab or
>fr->coulvdwtab or something else?  
>
This is an optimization hack; the fr->coulvdwtab contains both the 
coulomb and LJ parameters (i.e. 12 floats per table point), while the 
vdwtab contains only LJ (8 floats) and coultab only coulomb (4 floats). 
It improves the cache access pattern significantly!

If you call any loop that does tabulated interactions for both coulomb 
and LJ you must use the combined table, and if you call a LJ- or 
coulomb-only loop you should use the corresponding non-combined table.

Cheers,

Erik