[gmx-developers] calculating poteintial energy into do_md()

Mark Abraham Mark.Abraham at anu.edu.au
Mon Jul 16 15:19:44 CEST 2012


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.

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
>
>             Since I dont need to change lambda along the simulation, I
>             prefer using a
>             modified input topology and running a normal REMD changing
>             something in
>             repl_ex.c and md.c.
>
>             Of course, any suggestions to improve performance are
>             wellcome :)
>
>
>
>
>             Cheers,
>
>                 Berk
>
>
>                   Thanks in advance,
>
>                   Francesco
>
>
>
>
>
>
>                 --
>                 gmx-developers mailing list
>                 gmx-developers at gromacs.org
>                 <mailto: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
>                 <mailto:gmx-developers-request at gromacs.org>.
>
>
>
>             -- 
>             Cordiali saluti, Dr.Oteri Francesco
>             -- 
>             gmx-developers mailing list
>             gmx-developers at gromacs.org <mailto: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
>             <mailto:gmx-developers-request at gromacs.org>.
>
>
>
>     -- 
>     gmx-developers mailing list
>     gmx-developers at gromacs.org <mailto: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
>     <mailto:gmx-developers-request at gromacs.org>.
>
>
>
>
> -- 
> Cordiali saluti, Dr.Oteri Francesco
>
>


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120716/794d6fef/attachment.html>


More information about the gromacs.org_gmx-developers mailing list