[gmx-users] Modifying gromacs

[Mark Abraham]

On 29/10/2010 10:14 AM, Sai Pooja wrote:
> Hi,
> I want to modify the contribution to Interaction energy of different 
> groups - say I have groups A and B and I want the energy to be scaled 
> as E_AA + 0.1E_AB + 0.5*E_BB. Interaction parameters of each of these 
> groups are set by a forcefield
> Now after multiple correspondence on the gromacs list I have concluded 
> that there are 3 ways of doing this:
> 1. Using tables - for this I would have to list non-bonded parameters 
> for all atoms such that the combination rule and the table-potential 
> is used. For the table for BB interaction, scale the Coloumb and VdW 
> interactions in the tables by a factor of 0.5 and so on...

Sure. I think that for the above example, you'd need only 2 (maybe 3) 
table files.

> However, since tables would have to be supplied for pairs too 
> (tablep), it may not be accurate to supply 6-12 tables with coulomb 
> potential for these pairs. I am using CHARMM and the 1998 paper on 
> CHARMM says that in some specific cases the 1-4 interactions many be 
> scaled which makes me doubt this approach.

Yeah, something extra will be required here. Obviously, test and develop 
on a small toy system that you can compute by hand in a spreadsheet.
> 2. Forcefield parameters - By defining scaled [nonbonded_params] for 
> all relevant atoms. This will change the VdW interactions, but not 
> sure about the Coulomb interactions.

The Coulomb interactions are based on atomic charges, and there's no 
ready way to scale them differently for different interactions.

> 3. Modifying gromacs - by passing a parameter lambda to gromacs which 
> scales the force/potential by a factor lambda when gromacs calculates 
> force/potential.
> For implementing option 3, which programs in the gromacs package would 
> be the bes tstarting points for editing the energy contributions of 
> different groups/atoms?

I can think of two ways of approaching this. Efficiency requires that 
you use the energy group mechanism for your respective groups. Then you 
need to identify the generated non-bonded lists and do the necessary 
book-keeping to allocate scale factors to them. Then, either

a) pass the scale factor all the way into the non-bonded kernels and use 
them suitably there, or
b) create a copy of the force and energy arrays for each required scale 
factor, pass matching arrays into the appropriate kernels, and do the 
scaling on the respective arrays after all kernel contributions have 
been made, and then add the corresponding array elements.

Either way, you'll need a sound understanding of the code in 
src/gmxlib/nonbonded.c (and the kernels it calls, and how its data 
structures get created in the neighbour-searching). Make a test system, 
like a dipeptide in a small box of water, with energy groups, and use a 
debugger to see how things work.

a) is the easiest to develop. The kernels already have a "user data" 
pointer you can use, and it would be a fairly simple matter to adapt the 
generic C kernels to do your scaling. You will give up a lot of 
performance, however, unless you are prepared to adapt the 
hardware-optimized kernels similarly (not advised).

b) will use somewhat more memory, have a smaller code-change footprint, 
keep essentially the same performance

> As a more general question, how does one run a generalized Hamiltonian 
> REM on gromacs?

You can't. GROMACS REMD requires a normal MD Hamiltonian supplied in a .tpr.

Mark