[gmx-developers] Implementation of Bennet estimator for chemical potential

[JR Schmidt]

I am currently doing some chemical potential in GROMACS using the 
particle insertion method.  (For a variety of technical reasons related 
to the eventual system of interest, I would like to avoid thermodynamic 
integration.)  I would like to implement the Bennet estimator of 
chemical potential, which is the minimum variance estimator [see, for 
example, JCP, 123, 054105 (2005)].

Fundamentally, this requires doing BOTH trial insertions into an N 
particular ensemble, and particle DELETIONS from a N+1 particle 
ensemble, and accumulating the average of f(\beta*\Delta U), where f is 
a Fermi function.  (In practice, as long as N is large enough, both can 
be done from the same N particle trajectory.)  The former is already 
implemented in GROMACS via the TPI method.  I would like to implement 
the latter.

I located the do_tpi subroutine in minimize.c, and I follow the basic 
logic of it.  However, what I do NOT understand is how the call to 
"do_force" in minimize.c somehow returns in  enerd->term[F_EPOT] the 
component of the energy corresponding to the Nth particle, rather than 
the total energy of the system.  Obviously this is due to one of the 
arguments, but through the maze of subsequent function calls I can't 
track it down.

I would appreciate any guidance from the more experienceddevelopers.  I 
am an seasoned coder, but not familiar at all with the innards of 
GROMACS.  Yet I think this should be a rather simple addition if I can 
simply figure out this particular issue.

Thank you in advance.

-- 
J.R. Schmidt
Assistant Professor of Chemistry
Room 8305D
Department of Chemistry
University of Wisconsin-Madison
1101 University Ave
Madison, WI 53706

Phone: (608) 262-2996
Fax: (608) 262-9918
E-mail: schmidt at chem.wisc.edu
http://www.chem.wisc.edu/users/schmidt