[gmx-developers] Re: Implementation of Bennet estimator for chemical potential (JR Schmidt)

Michael Shirts mrshirts at gmail.com
Wed Jul 15 07:13:33 CEST 2009 

> 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)].

It seems like this could be implemented in postprocessing with a free
energy calculation run relatively easily. You do two runs, at lambda=0
and lambda=1, with one particle that is simulated with the free
energy.

For the lambda=1 run (particle disconnected), then the particle is
free to go anywhere in the system. So you collect configurations, and
do an md -rerun on those configurations, but now with lambda=0, and
collect the potential energy differences at these configurations.

The lambda=0 run doesn't necessarily have to have to be run with free
energies on, just a normal one. For this one, you reanalyze, each time
generating a new toplogy with a randomly generated particle being the
"free energy" particle, and rerun with lambda=1, and collect the
energies. This run will be a bit messier, since you do have to rerun
run grompp each time. But it could be very easily parallelized.

This does assume that you are turning off the solute-solvent terms,
and not the solute-solvent terms -- which I believe was being
implemented, but I can't remember!

Of course, if you have a very time consuming simulation, then this
might end up being a little too expensive to do the reruns. But
depending on the correlation time of the system, if you only need
every 1000 frame, it might not be so bad.

Note you don't necessarily want to accumulate the Fermi function; you
really want to be averaging the Fermi function of
(f\beta(\deltaU-\delta mu)) -- which you don't know to start out with.
If you collect the energy differences, you can postprocess and
determine the free energy self consistently.

Note also that if you don't want to solve self-consistently, then if
you have -guess- for the chemical potential, then accumulating
f\beta\delta(U-\delta(guess for mu)) will be a more efficient
estimator than accumulating f(\beta\deltaU). Or you can maintain
several averages, at different (guess for mu) values, and value
computed with the guess for mu closest to G will be the most precise.

Note also that if you collect different numbers of insertion or
deletion works, then instead of \delta U - \delta mu, \delta U-\delta
mu + T ln (N_del/N_insert) is the minimum variance constant. (it may
be N_insert/N_del, I can't quite remember at this point).

My understanding from David Kofke's work is that the particle
insertions contribute more to the information generated than the
deletions, so you would probably get better precision spending more
time on the insertions anyway.

Cheers,
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia

>
>
> 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
>
>
>
>
>
>
>
> ------------------------------
>
>
>
> Message: 3
>
> Date: Tue, 14 Jul 2009 22:05:32 +0200 (CEST)
>
> From: hess at sbc.su.se
>
> Subject: Re: [gmx-developers] Implementation of Bennet estimator for
>
>        chemical potential
>
> To: "Discussion list for GROMACS development"
>
>        gmx-developers at gromacs.org>
>
> Message-ID: 49323.81.234.234.16.1247601932.squirrel at mail.sbc.su.se>
>
> Content-Type: text/plain;charset=iso-8859-1
>
>
>
> 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.
>
> >
>
>
>
>
>
>
>
> ------------------------------
>
>
>
> Message: 4
>
> Date: Tue, 14 Jul 2009 15:30:32 -0500
>
> From: JR Schmidt schmidt at chem.wisc.edu>
>
> Subject: Re: [gmx-developers] Implementation of Bennet estimator for
>
>        chemical        potential
>
> To: Discussion list for GROMACS development
>
>        gmx-developers at gromacs.org>
>
> Message-ID: 4A5CEAE8.20000 at chem.wisc.edu>
>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> 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
>
>
>
>
>
>
>
> ------------------------------
>
>
>
> Message: 5
>
> Date: Tue, 14 Jul 2009 22:35:57 +0200 (CEST)
>
> From: hess at sbc.su.se
>
> Subject: Re: [gmx-developers] Implementation of Bennet estimator for
>
>        chemical potential
>
> To: "Discussion list for GROMACS development"
>
>        gmx-developers at gromacs.org>
>
> Message-ID: 49767.81.234.234.16.1247603757.squirrel at mail.sbc.su.se>
>
> Content-Type: text/plain;charset=iso-8859-1
>
>
>
> 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.
>
> >
>
>
>
>
>
>
>
> ------------------------------
>
>
>
> Message: 6
>
> Date: Wed, 15 Jul 2009 12:49:24 +1000
>
> From: Mark Abraham Mark.Abraham at anu.edu.au>
>
> Subject: Re: [gmx-developers] Single Point Energy Calculation
>
> To: Discussion list for GROMACS development
>
>        gmx-developers at gromacs.org>
>
> Message-ID: 4A5D43B4.9030709 at anu.edu.au>
>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> XAvier Periole wrote:
>
> > I think the best is to call mdrun using nstep = 0
>
>
>
> ... and an EM integrator.
>
>
>
> Mark
>
>
>
> >> Hi,
>
> >>
>
> >> I want to write an analysis script that calculates a single point
>
> >> energy without invoking mdrun. Since, the tpr file contains all the FF
>
> >> and topology information, is there a function that I can pass the TPR
>
> >> information to that will return the total PE of the system?
>
> >>
>
> >> Thanks,
>
> >>
>
> >> Ilya
>
> >>
>
> >>
>
> >> --
>
> >> Ilya Chorny Ph.D.
>
> >>
>
> >> _______________________________________________
>
> >> 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
>
> >&gt%

More information about the gromacs.org_gmx-developers mailing list