[gmx-developers] calculating poteintial energy into do_md()
Berk Hess
hess at kth.se
Mon Jul 16 18:25:24 CEST 2012
On 07/16/2012 06:08 PM, francesco oteri wrote:
>
>
> 2012/7/16 Berk Hess <hess at kth.se <mailto:hess at kth.se>>
>
> Using lambda scaling for non-bonded interactions is going to be
> very slow.
>
>
> If it is done with the actual implementation, you are right. But if
> mdrun recognize rest2=yes thzn it rescales the
> parameters at teh beginning. Then a normal MD run goes.
>
>
> But this would scale the protein-protein non-bonded interactions
> with lambda^2.
>
>
>
> Why? computing charge-charge interaction is:
>
> S*qi*S*qj = (S^2 )*qi*qj
> but S=sqrt(Bi/B)
> so what I obtain is:
> (Bi/B)* qi*qj
>
> That is exactly what I want: rescaling by (Bi/B) protein-protein
> interactions
Ah, sorry, that works indeed.
Cheers,
Berk
>
> 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
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120716/3356925e/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list