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