[gmx-developers] Implementation of Bennet estimator for chemical potential
hess at sbc.su.se
hess at sbc.su.se
Tue Jul 14 22:35:57 CEST 2009
Swapping coordinates is probably the simplest solution code-wise.
But that would mean the neighbor search would have to be repeated
for all molecules, which is computationally very expensive.
> 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.
>>> 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
>>> gmx-developers mailing list
>>> gmx-developers at gromacs.org
>>> 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
>> 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
> 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
> gmx-developers mailing list
> gmx-developers at gromacs.org
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
More information about the gromacs.org_gmx-developers