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

[hess at sbc.su.se]

Hi,

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.

Berk

> 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.
>>>
>>> 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.
>>>
>>>
>>
>> _______________________________________________
>> 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
> 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.
>