[gmx-developers] Editing the GROMACS energy function

[Radhakrishna Bettadapura]

I want to add a term E to the MD potential energy function minimized by gromacs.  The function E operates on the set of input atoms in a complicated way and returns a real value.  I also have a closed-form representation of the first derivative of E with respect to each atom position.  

I've managed to dig up a few questions related to mine in the mailing list archives.  Specifically, 

http://lists.gromacs.org/pipermail/gmx-developers/2004-February/000767.html

says that this is as easy as modifying calc_f_el in sim_util.c.  However, the very next message

http://lists.gromacs.org/pipermail/gmx-developers/2004-February/000768.html

seems to contradict this, saying that force.c has to be modified. Further, these messages are from 2004, and I haven't been able to locate anything more recent that is also more concrete.  For instance, this message asks roughly the same question

http://lists.gromacs.org/pipermail/gmx-developers/2009-March/003118.html


but the thread peters out after a few messages with no obvious conclusion/solution.  

Anyway, I've explored the code a little bit, and I believe the following changes are necessary if I'm to add a term to the energy function.  

1) The function most relevant seems to be do_force_lowlevel (force.c). This function is called by do_force in sim_util.c, and appears to compute both the potential energy, stored in gmx_enerdata_t* enerd, and the force, stored in rvec f[].  It seems that the right way to proceed would be to do the following: 

a) Edit the gmx_enerdata_t data structure so that it contains an entry for the new term E. This entails editing F_NRE, the size of the array  term in that data structure (include/types/forcerec.h).  Presumably,  since F_NRE is the last entry in the enumeration data structure in include/types/idef.h, adding an entry like F_NEW to the enum data type just before F_NRE will automagically make (int)F_NRE return the right number of interactions. Correct?

b) Edit do_force_lowlevel such that  it adds the gradient due to my new energy function E to f[i], where i is each atom.  Also add the energy contribution due to E thus: E[F_NEW] = //funky new addition to energy function. 

2) One worry is that the force f[i]/position x[i] doesn't really point to the force on/position of the i-th atom.  I suspect this because of lines such as the following: 

//src/gmxlib/nonbonded.c
 nb_kernel_allvsallgb_sse2_double(fr,mdatoms,excl,x[0],f[0],egcoul,egnb,egpol,
                                                 &outeriter,&inneriter,&fr->AllvsAll_work);

Why are f[0] and x[0] special in the above function call?  

3) Finally, I need to edit the input processing subroutines in such a way that I can turn my interaction on and off at will.  I'm not sure how exactly to do this, but I suspect that the function parse_common_args in src/gmxlib/statutil.c is relevant.  I'd appreciate any "input" on this (pun unintended). 

4) I'd like to know if there's something I haven't covered in the above comments/ questions.  

Also, it'd be great if someone reading this were a) at the University of Texas at Austin, where I'm a graduate student, and b) willing to spare some of their time, not more than a half-hour, to go over the code with me.  I promise I'll make it worth their while; we could get a coffee at JP's Java!  


thanks for your time,

Radhakrishna Bettadapura
PhD candidate,
Dept. of Mechanical Engineering/Computational Visualization Center,
University of Texas at Austin
http://users.ices.utexas.edu/~brrk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120329/75f905ae/attachment.html>