[gmx-developers] calculating poteintial energy into do_md()
Berk Hess
hess at kth.se
Mon Jul 16 17:56:38 CEST 2012
Using lambda scaling for non-bonded interactions is going to be very slow.
But this would scale the protein-protein non-bonded interactions with
lambda^2.
Is that what you want?
Cheers,
Berk
On 07/16/2012 05:06 PM, francesco oteri wrote:
> In the paper they just say that non bonded rescaling is obtained
> rescaling
> LG epsilon by Bi/B and the protein charges by sqrt(Bi/B).
>
> I did simulation using the lambda approach and I qualitatively reproduced
> the results.
>
> Just to give my contribution, I think the most efficient way to
> implement REST2
> is using the lambda approach as follows:
>
> 1) in .mdp introducing a new keyword like rest2 = yes/no
> 2) using the init_lambda value to rescale the parameters at the
> bootstrap of md
> (es. in do_md or init_replica_exchange).
> 3) running the code as normal MD and swapping the states
>
> Since point 1 implies modification of grompp and point 2 requires the
> knowledge of
> the data structure, I am just implementing point 3 while point 2 is
> demanded to external
> scripts.
>
> Francesco
>
> 2012/7/16 Berk Hess <hess at kth.se <mailto:hess at kth.se>>
>
> On 07/16/2012 04:42 PM, Shirts, Michael (mrs5pt) wrote:
>
> One question, since better REST2 is one item on the todo list
> for 5.0, and I
> want to be thinking about the possibilities.
>
> For solute-solvent energy rescaling, how does one calculate
> the changes in
> solute-solvent energy without doing multiple PME calls to
> decompose the
> electrostatic energy?
>
> I'm guessing that by since your are load multiple topologies
> of the system,
> for each of those topologies, you can explicitly write new
> nonbonded
> pairwise parameter terms between the water and protein. This
> is a pain to
> do (need new ij terms for every parameter pair) but
> straightforward in the
> end to automate. Is this how your are planning to do it?
>
> Without PME, then it's straightforward to decompose into
> solute-solvent
> energy / solute-solute / solvent-solvent electrostatics,
> though for adding
> long term to Gromacs, we'd want to support PME.
>
> No, even without PME this is not straightforward.
> Using a plain cut-off's for electrostatics is horrible, so I'd say
> you want to use reaction-field.
> But with reaction-filed the correction terms are non pair-wise, as
> with PME.
>
> Cheers,
>
> Berk
>
>
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu <mailto:michael.shirts at virginia.edu>
> (434)-243-1821
>
>
> From: francesco oteri <francesco.oteri at gmail.com
> <mailto:francesco.oteri at gmail.com>>
> Reply-To: Discussion list for GROMACS development
> <gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>>
> Date: Mon, 16 Jul 2012 16:13:16 +0200
> To: Discussion list for GROMACS development
> <gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>>
> Subject: Re: [gmx-developers] calculating poteintial
> energy into do_md()
>
> Hi
>
> 2012/7/16 Mark Abraham <Mark.Abraham at anu.edu.au
> <mailto:Mark.Abraham at anu.edu.au>>
>
> On 16/07/2012 10:59 PM, francesco oteri wrote:
>
> Hi,
> what I am trying to do is implementing the REST2
> tecnique [1]
> The goal is enanching the conforlational sampling
> splitting the system in
> two parts: solvent and solute.
> The two groups have a different APPARENT temperature
> because the bonded
> interaction within the solute atoms are rescaled by a
> factor Bi/B (
> Bi=1/Ti*Kb, Ti = apparent temperature for the i-th
> replica, B=Boltzmann
> constant) and the protein-water interactions are
> rescaled by sqrt(Bi/B).
> So the system temperature, regulated by the
> thermostat, is equal along the
> different replicas but the dynamic of the different
> replicas
> changes becaus of the different Boltzmann factor.
> To implement the tecnique, the naive approach (adopted
> in [2]) is:
> 1) Generating N-different topology (N= number of
> replicas).
> 2) Removing same check at the beginning of the Replica
> Exchange code, in
> order to run N replicas with equal run parameters
> (temperature, pression,
> ecc. ecc.)
> 3) Changing the acceptance ratio formula.
>
> This approach has been developped for gmx4.0.3 but I
> need to use
> gmx4.5.5 because it has a native implementation of the
> CHARMM force-field.
> Since the code between the two versions is different I
> am trying to figure
> out how to reinvent the wheel.
>
> As shown in [3] and adopted in [2] the same result
> can be obtained using
> the lamda dynamics in gromacs.In this case, only one
> topology has to be
> generated. The state A correspond to the replica at
> the lowest apparent
> temperature, while the state B represents the highest
> temperature replica.
> Intermediate replicas are generated changing the
> init_lambda value.
>
> What has been found, this approach results in
> slowest run.
>
>
> So, right now I got stuck at the point 3. In fact
> to find:
>
> delta = B( Va(Xb) + Vb(Xa) - Va(Xa) - Vb(Xb))
>
> I need to calculate Va(Xb) and Vb(Xa)
>
>
> Va(Xb) is determined by the potential of replica a and
> the coordinates of
> replica b, so it is most effectively evaluated on the
> processor that has
> the coordinates of replica b. If the topology is the
> same at the replicas a
> and b, then can Va(Xb) be computed by replica b by
> rescaling suitably the
> quantities that contributed to Vb(Xb), based on the
> knowledge that replica
> a is the exchange partner? Then the quantities that
> are exchanged before
> the test are Va(Xb) and Vb(Xa). Should be as fast as
> regular REMD.
>
>
> This is my idea, indeed
>
>
> Mark
>
>
>
> [1] http://pubs.acs.org/doi/abs/10.1021/jp204407d
> [2]http://pubs.acs.org/doi/abs/10.1021/ct100493v
> [3]
> http://onlinelibrary.wiley.com/doi/10.1002/jcc.21703/abstract
>
>
>
>
>
> 2012/7/16 Berk Hess <hess at kth.se <mailto:hess at kth.se>>
>
> On 7/16/12 14:18 , Shirts, Michael (mrs5pt) wrote:
>
> Actually I did it but I the performance
> decrease a lot (2x).
>
> In particular I run the same run with
> free_energy = yes and
> free_energy =
> no and in the second case the run is
> 2times faster.
> I found in the mailing list that this
> performance drop is known.
>
> What are the particular parameters you are
> trying to change? I wrote the
> new Hamiltonian Exchange code, so I'm
> interested in making sure it's as
> useful as possible. There might be some cases
> where this performance
> drop
> can be avoided, and It's likely the PME
> parameters, since doing
> repeated
> that's the only thing that really makes free
> energy slow.
>
> UNLESS of course you are changing ALL of the
> atoms, which which case all
> of
> them will go through the free energy code, and
> it certainly will be
> slower.
> There are some possible ways to fix this, but
> that will not be until 5.0.
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> <mailto:michael.shirts at virginia.edu>
> (434)-243-1821
>
>
> From: francesco oteri
> <francesco.oteri at gmail.com
> <mailto:francesco.oteri at gmail.com>>
>
> Reply-To: Discussion list for GROMACS
> development <
> gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>>
> Date: Mon, 16 Jul 2012 14:08:49 +0200
> To: Discussion list for GROMACS
> development <gmx-developers at gromacs.org
> <mailto:gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] calculating
> poteintial energy into do_md()
>
> Dear Berk
>
> 2012/7/16 Berk Hess <hess at kth.se
> <mailto:hess at kth.se>>
>
> On 7/16/12 13:05 , francesco oteri wrote:
>
> Dear gromacs developers,
> I am trying to implement in gromacs
> 4.5.5 a particular Hemiltonian
> Replica
> Exchange tecnique.
> Right now, at every exchange attempt
> in do_md, I figured out how to
> access at the potential
> energy of replica A (=configuration A
> at temperature A) and B
> (=configuration B at temperature B)
> and so on.
>
> In my case, each replica has the same
> temperature, but there is a
> different Hemiltonian equation for
> every replica.
> The different Hemiltonian are obtained
> simply changing the force field
> parameters in the input topology so I dont
> need to modify anything in gromacs.
>
> But at every exchange attempt I
> have to test if the configuration B
> can
> exist in the state A so I need to
> calculate
> its potential energy using the force
> field data of replica A.
>
> I found that function sum_epot
> calculates the potential energy but I
> suspect that it uses values calculated in
> do_force since sum_epot is called by
> do_force_lowlevel in turn called
> by
> do_force.
>
>
>
> So my question is, should I call
> do_force in replica A with
> coordinates
> from replica B reach my goal?
>
> Yes.
> Are you changing bonded and/or
> non-bonded parameters?
> Some non-bonded parameters might be
> preprocessed, so you might need
> reprocess those
> and them reprocesses again to get back
> to the original state.
>
>
> Which parameters need such a
> preprocessing?
>
>
>
> Note that 4.6 will have proper
> Hamiltonian replica exchange, but that
>
> will
> use the lambda coupling
> parameter approach. If you need to do
> something similar, it might be
> much
> simpler to use this code.
>
>
> Actually I did it but I the
> performance decrease a lot (2x).
>
> In particular I run the same run with
> free_energy = yes and
> free_energy =
> no and in the second case the run is
> 2times faster.
> I found in the mailing list that this
> performance drop is known.
>
> I think the performance drop is due to the,
> temporarily, missing sse
>
> non-bonded kernels.
> They should be back somewhere this week.
>
> Cheers,
>
> Berk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120716/f4f6305a/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list