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

[hess at sbc.su.se]

Hi,

There are a few small "tricks" that make this work.
In the do_force call in do_tpi the bonded flag is not given,
such that you get only the non-bonded interactions.
In ns.c the fr->n_tpi parameter sets that only the non-bonded
interactions with the test molecule are calculated.
Furthermore, there in several places in the code, mainly in force.c,
the fr->n_tpi parameter is used to get only the test molecule
contribution, for instance for the long-range dispersion correction.

I would assume particle deletion is more complicated,
since you probably want to try to remove all molecules,
not just the last one as in insertition.

I don't know if it is of use to you, but I implemented the BAR
equivalent of TPI, including g_bar, in the git master branch.

Berk

> 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
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>