[gmx-developers] calculating poteintial energy into do_md()
Shirts, Michael (mrs5pt)
mrs5pt at eservices.virginia.edu
Mon Jul 16 18:23:07 CEST 2012
OK, let's go back to the basics:
For REST 2, you need to scale the solute-solute interactions by T_0/T_M, and
the solute-solvent interactions by sqrt(T_0/T_M), with the solvent-solvent
interactions unchanged.
So if the protein charges are scaled by sqrt of the temperature ratio at
each new topology, then their interactions with themselves will be changed
appropriately, and their interactions with the water will be changed
appropriately. This can be done with lambda currently (but could be done
in other ways as well) by setting the correct initial and final charges
correctly, and choosing the right lambdas. So actually, no extra PME calls
required.
And this should be the same with the LJ epsilons -- if they are scaled by
sqrt of the temperature ratio, then they will have the proper interactions.
Again, can be done by lambda, in the same way as the charges.
Solute nonbonded should be scaled with the full temperature ratio.
Do I have this right?
This should be trivial to set up with the new Hamiltonian exchange
formalism, at the cost of some computational speed. I can update the manual
to describe explicitly how to set it up.
With some tweaking, it should be possible to make nearly costless in 5.0,
especially with some other ideas I have for free energy calculations, but
can't be done in 4.6, since the release is imminent.
In terms of how to do it best now -- I can't really say. I'd argue you'd be
better of taking the performance hit than having to recode and debug
something, but that's your choice. I'm happy to look at any configuration
and offer help on how to set it up using lambda and possibly speed things
up.
But as Berk mentioned, better to wait until after the accelerated nonbonded
interactions are added back in.
Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821
> From: Berk Hess <hess at kth.se>
> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Date: Mon, 16 Jul 2012 17:56:38 +0200
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] calculating poteintial energy into do_md()
>
> 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
>>
>
> --
> 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