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

JR Schmidt schmidt at chem.wisc.edu
Tue Jul 14 22:30:32 CEST 2009

Thanks for the quick reply.  To summarize, it sounds like the last
(N+1st) particle is "special", and this is special case is incorporated
in several locations throughout the code.  Deletions are more
complicated, as one would ideally like to remove all molecules to
improve the statistics.

I believe one can workaround this by simply "swapping" the coordinates
of the molecule one wants to delete with the N+1st molecules.  Then the
do_force call will yield in interaction energy of the deleted molecule
with the rest of the system.  The \Delta U desired is just the negative
of this quantity.
> 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.
>>
>>
>> --
>> J.R. Schmidt
>> Assistant Professor of Chemistry
>> Room 8305D
>> Department of Chemistry
>> 1101 University Ave
>>
>> 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.
>>
>>
>
> _______________________________________________
> 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.
>

--
J.R. Schmidt
Assistant Professor of Chemistry
Room 8305D
Department of Chemistry