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

hess at sbc.su.se hess at sbc.su.se
Wed Jul 15 09:17:28 CEST 2009 

Hi,

You can indeed do it with lambda-coupling of one particle.
But it is of course far more efficient to use all particles.
This can be done without changes to the code.
You can put each molecule in a separate energy group.
That you can simulate directly, or do a rerun if simulation
with so many energy groups is too slow.
You then need to write a script to sum the g_energy output
of the whole matrix of LJ and Coulomb energy terms into a vector
for all molecules.
The maximum number of energy groups is, I think, 255.
So if you have more than 255 molecules, you can make 254
groups of one molecule and one group of the rest.

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

More information about the gromacs.org_gmx-developers mailing list