<html>
  <head>
    <meta content="text/html; charset=windows-1252"
      http-equiv="Content-Type">
  </head>
  <body bgcolor="#FFFFFF" text="#000000">
    VDW is shifted (vdw-modifier=Potential-shift). The simulations are
    run in single precision.<br>
    Currently Im re-runing them in DP.<br>
    <br>
    Best,<br>
    Bernhard<br>
    <br>
    <div class="moz-cite-prefix">Am 15/07/15 um 21:09 schrieb Shirts,
      Michael R. (mrs5pt):<br>
    </div>
    <blockquote
      cite="mid:D1CC2797.6ED7C%25mrs5pt@eservices.virginia.edu"
      type="cite">
      <meta http-equiv="Content-Type" content="text/html;
        charset=windows-1252">
      <div>
        <div>Just wondering, Is this with single or double precision?</div>
        <div>
          <div><br>
          </div>
          <div>Also, is the VDW interaction tapered or shifted?  (check
            the log file or mdp out) If not, that can easily cause lack
            of energy conservation.</div>
          <div><br>
          </div>
          <div>Best,</div>
          ~~~~~~~~~~~~<br>
          Michael Shirts<br>
          Associate Professor<br>
          Department of Chemical Engineering<br>
          University of Virginia<br>
          <a moz-do-not-send="true" href="michael.shirts@virginia.edu">michael.shirts@virginia.edu</a><br>
          (434) 243-1821</div>
      </div>
      <div><br>
      </div>
      <span id="OLK_SRC_BODY_SECTION">
        <div style="font-family:Calibri; font-size:14pt;
          text-align:left; color:black; BORDER-BOTTOM: medium none;
          BORDER-LEFT: medium none; PADDING-BOTTOM: 0in; PADDING-LEFT:
          0in; PADDING-RIGHT: 0in; BORDER-TOP: #b5c4df 1pt solid;
          BORDER-RIGHT: medium none; PADDING-TOP: 3pt">
          <span style="font-weight:bold">From: </span>Bernhard &lt;<a
            moz-do-not-send="true" href="mailto:b.reuter@uni-kassel.de">b.reuter@uni-kassel.de</a>&gt;<br>
          <span style="font-weight:bold">Reply-To: </span>"<a
            moz-do-not-send="true"
            href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>"
          &lt;<a moz-do-not-send="true"
            href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>&gt;<br>
          <span style="font-weight:bold">Date: </span>Wednesday, July
          15, 2015 at 11:09 AM<br>
          <span style="font-weight:bold">To: </span>"<a
            moz-do-not-send="true"
            href="mailto:gromacs.org_gmx-developers@maillist.sys.kth.se">gromacs.org_gmx-developers@maillist.sys.kth.se</a>"
          &lt;<a moz-do-not-send="true"
            href="mailto:gromacs.org_gmx-developers@maillist.sys.kth.se">gromacs.org_gmx-developers@maillist.sys.kth.se</a>&gt;<br>
          <span style="font-weight:bold">Subject: </span>Re:
          [gmx-developers] Drift in Conserved-Energy with Nose-Hoover
          thermostat<br>
        </div>
        <div><br>
        </div>
        <div>
          <div bgcolor="#FFFFFF" text="#000000">Finally I also tried a
            NVE simulation (rlist=1.012nm) with the same system,
            parameters and starting conditions but
            <u>without</u><u> any constraints</u> (constraints=none):
            The result is an equal linear drift (-0.78 kJ/mol/ps or
            -3.43*10^-5 kJ/mol/ps per atom) of the total energy as for
            NVE  with constraints=h-bonds or of the conserved energy for
            Nose-Hoover with constraints=h-bonds (rlist=1.012nm in all
            cases; lincs-iter=2 for NVE and lincs-iter=1 for NVT).<br>
            So constraints should be ruled out as cause of the drift...<br>
            <br>
            Best,<br>
            Bernhard<br>
            <br>
            <div class="moz-cite-prefix">Am 15/07/15 um 18:14 schrieb
              Mark Abraham:<br>
            </div>
            <blockquote
cite="mid:CAMNuMARv=T2G2hKBtn1hY_WNVXsaiUh-u9o0B31DVebzBDO8vQ@mail.gmail.com"
              type="cite">
              <div dir="ltr">Hi,
                <div><br>
                </div>
                <div>This looks approximately normal for constrained
                  dynamics of systems with water. See e.g. fig 8 of <a
                    moz-do-not-send="true"
                    href="http://dx.doi.org/10.1016/j.cpc.2013.06.003">http://dx.doi.org/10.1016/j.cpc.2013.06.003</a><br>
                  <br>
                  Mark</div>
                <div><br>
                  <div class="gmail_quote">
                    <div dir="ltr">On Wed, Jul 15, 2015 at 5:57 PM
                      Bernhard &lt;<a moz-do-not-send="true"
                        href="mailto:b.reuter@uni-kassel.de">b.reuter@uni-kassel.de</a>&gt;
                      wrote:<br>
                    </div>
                    <blockquote class="gmail_quote" style="margin:0 0 0
                      .8ex;border-left:1px #ccc solid;padding-left:1ex">
                      I now also tried with bigger manual rlist values
                      of 1.02 and 1.05: This<br>
                      doesnt help - the drift of the conserved energy
                      with Nose-Hoover gets<br>
                      even bigger: -1,8 kJ/mol/ps (7.9*10^-5 kJ/mol/ps
                      per atom) and -2.62<br>
                      kJ/mol/ps (1.15*10^-4 kJ/mol/ps per atom).<br>
                      <br>
                      Am 15/07/15 um 17:41 schrieb Bernhard:<br>
                      &gt; I mean a drift of the total energy in NVE -
                      while with Nose-Hoover the<br>
                      &gt; drift is in the Conserved-Energy quantity of
                      g_energy (the total<br>
                      &gt; energy shows no drift with Noose-Hoover...).<br>
                      &gt;<br>
                      &gt; Am 15/07/15 um 17:38 schrieb Bernhard:<br>
                      &gt;&gt; I also did a NVE simulation with the same
                      parameters, system and<br>
                      &gt;&gt; starting conditions but with manually set
                      rlist=1.012nm (since<br>
                      &gt;&gt; verlet-buffer-drift doesnt work in NVE):<br>
                      &gt;&gt; There I also got a linear drift (but
                      smaller) of 0.78 kJ/mol/ps<br>
                      &gt;&gt; (3.436*10^-5 kJ/mol/ps per atom).<br>
                      &gt;&gt; For comparison reasons I also did a NVT
                      Nose-Hoover Simulation with<br>
                      &gt;&gt; manually set rlist=1.012nm:<br>
                      &gt;&gt; There I got a comparable linear drift of
                      0.67 kJ/mol/ps (2.94*10^-5<br>
                      &gt;&gt; kJ/mol/ps per atom).<br>
                      &gt;&gt; So no differences between NVE and NVT so
                      far in my opinion...<br>
                      &gt;&gt;<br>
                      &gt;&gt; Best,<br>
                      &gt;&gt; Bernhard<br>
                      &gt;&gt;<br>
                      &gt;&gt; Am 15/07/15 um 17:10 schrieb Shirts,
                      Michael R. (mrs5pt):<br>
                      &gt;&gt;&gt; The conserved quantity in nose-hoover
                      is not quite as good as the<br>
                      &gt;&gt;&gt; conserved energy, which should have
                      no drift at all.  For NH, the<br>
                      &gt;&gt;&gt; conserved quantity should drift as a
                      random Gaussian process with mean<br>
                      &gt;&gt;&gt; zero (i.e. go with sqrt(N)). It
                      shouldn't be drifting linearly.<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt; I would check to see if your system
                      conserved energy when run with NVE<br>
                      &gt;&gt;&gt; (use the endpoint of the NPT
                      simulation).  It's easier to diagnose any<br>
                      &gt;&gt;&gt; problems with an NVE simulation,
                      which should have virtually no<br>
                      &gt;&gt;&gt; drift, vs<br>
                      &gt;&gt;&gt; a NVT simulation, which has random
                      noise drift.  Odds are, if there<br>
                      &gt;&gt;&gt; is a<br>
                      &gt;&gt;&gt; problem with the NVT simulation, it
                      will also show up in the NVE<br>
                      &gt;&gt;&gt; simulation if only the thermostat is
                      removed.<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt; Also, consider looking at <a
                        moz-do-not-send="true"
                        href="http://pubs.acs.org/doi/abs/10.1021/ct300688p"
                        rel="noreferrer" target="_blank">
                        http://pubs.acs.org/doi/abs/10.1021/ct300688p</a><br>
                      &gt;&gt;&gt; for tests of whether the ensemble
                      generated is correct.<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt; Best,<br>
                      &gt;&gt;&gt; ~~~~~~~~~~~~<br>
                      &gt;&gt;&gt; Michael Shirts<br>
                      &gt;&gt;&gt; Associate Professor<br>
                      &gt;&gt;&gt; Department of Chemical Engineering<br>
                      &gt;&gt;&gt; University of Virginia<br>
                      &gt;&gt;&gt; <a moz-do-not-send="true"
                        href="mailto:michael.shirts@virginia.edu"
                        target="_blank">
                        michael.shirts@virginia.edu</a><br>
                      &gt;&gt;&gt; (434) 243-1821<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt; On 7/15/15, 10:58 AM, "Bernhard" &lt;<a
                        moz-do-not-send="true"
                        href="mailto:b.reuter@uni-kassel.de"
                        target="_blank">b.reuter@uni-kassel.de</a>&gt;
                      wrote:<br>
                      &gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; Dear Gromacs Users and
                      Developers,<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; I have a problem regarding energy
                      conservation in my 10ns NVT<br>
                      &gt;&gt;&gt;&gt; protein+water+ions (22765 atoms)
                      production (minimization and<br>
                      &gt;&gt;&gt;&gt; equilibration for more than 15ns
                      was carried out in NPT before)<br>
                      &gt;&gt;&gt;&gt; simulations using a Nose-Hoover
                      thermostat (tau=2.5ps).<br>
                      &gt;&gt;&gt;&gt; On first glance everything looks
                      fine - the potential, kinetic and<br>
                      &gt;&gt;&gt;&gt; total<br>
                      &gt;&gt;&gt;&gt; energy are nearly perfectly
                      constant (with normal fluctuations) - but<br>
                      &gt;&gt;&gt;&gt; when I checked the
                      "Conserved-Energy" quantity that g_energy outputs
                      I<br>
                      &gt;&gt;&gt;&gt; had to recognize a significant
                      (nearly perfectly) linear downward<br>
                      &gt;&gt;&gt;&gt; drift<br>
                      &gt;&gt;&gt;&gt; of this "to-be-conserved"
                      quantity of around 1.7 kJ/mol/ps (7.48*10^-5<br>
                      &gt;&gt;&gt;&gt; kJ/mol/ps per atom).<br>
                      &gt;&gt;&gt;&gt; This appears somehow disturbing
                      to me since I would expect that this<br>
                      &gt;&gt;&gt;&gt; Conserved-Energy is the conserved
                      energy of the extended Nose-Hoover<br>
                      &gt;&gt;&gt;&gt; Hamiltonian - which should by
                      definition be conserved.<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; If it would be a drift caused by
                      normal round-off error due to single<br>
                      &gt;&gt;&gt;&gt; precision I would expect it to
                      grow with Sqrt(N) and not with N<br>
                      &gt;&gt;&gt;&gt; (linear)<br>
                      &gt;&gt;&gt;&gt; (N=number of steps).<br>
                      &gt;&gt;&gt;&gt; So I would like to know if this
                      is a normal behaviour and also what<br>
                      &gt;&gt;&gt;&gt; could cause this (buffer size,
                      precision, constraints etc)?<br>
                      &gt;&gt;&gt;&gt; Also I would like to know, if I
                      am correct with my guess that the<br>
                      &gt;&gt;&gt;&gt; "Conserved-Energy" quantity is in
                      this case the energy of the extended<br>
                      &gt;&gt;&gt;&gt; Nose-Hoover Hamiltonian?<br>
                      &gt;&gt;&gt;&gt; The .mdp file is atatched (don't
                      be confused about rlist=1 - since Im<br>
                      &gt;&gt;&gt;&gt; using the Verlet-scheme the
                      verlet-buffer-drift option should be by<br>
                      &gt;&gt;&gt;&gt; default active and determine the
                      rlist value (Verlet buffer-size)<br>
                      &gt;&gt;&gt;&gt; automatically).<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; Best regards,<br>
                      &gt;&gt;&gt;&gt; Bernhard<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; ; Run parameters<br>
                      &gt;&gt;&gt;&gt; integrator    = md        ;
                      leap-frog integrator<br>
                      &gt;&gt;&gt;&gt; nsteps        = 10000000    ;
                      10000 ps = 10 ns<br>
                      &gt;&gt;&gt;&gt; dt        = 0.001        ; 1 fs<br>
                      &gt;&gt;&gt;&gt; ; Output control<br>
                      &gt;&gt;&gt;&gt; nstxout        = 5000        ;
                      save coordinates every ps<br>
                      &gt;&gt;&gt;&gt; nstvout        = 5000        ;
                      save velocities every ps<br>
                      &gt;&gt;&gt;&gt; nstxtcout    = 1000        ; xtc
                      compressed trajectory output every ps<br>
                      &gt;&gt;&gt;&gt; nstenergy    = 1000        ; save
                      energies every ps<br>
                      &gt;&gt;&gt;&gt; nstlog        = 5000        ;
                      update log file every ps<br>
                      &gt;&gt;&gt;&gt; ; Bond parameters<br>
                      &gt;&gt;&gt;&gt; continuation    = yes        ;
                      continue from NPT<br>
                      &gt;&gt;&gt;&gt; constraint_algorithm = lincs    ;
                      holonomic constraints<br>
                      &gt;&gt;&gt;&gt; constraints    = h-bonds    ; all
                      bonds (even heavy atom-H bonds)<br>
                      &gt;&gt;&gt;&gt; constrained<br>
                      &gt;&gt;&gt;&gt; lincs_iter    = 1        ;
                      accuracy of LINCS<br>
                      &gt;&gt;&gt;&gt; lincs_order    = 4        ; also
                      related to accuracy<br>
                      &gt;&gt;&gt;&gt; ; Neighborsearching<br>
                      &gt;&gt;&gt;&gt; cutoff-scheme    = Verlet    ;
                      Verlet cutoff-scheme instead of<br>
                      &gt;&gt;&gt;&gt; group-scheme (no charge-groups
                      used)<br>
                      &gt;&gt;&gt;&gt; ns_type        = grid        ;
                      search neighboring grid cells<br>
                      &gt;&gt;&gt;&gt; nstlist        = 10        ; 10
                      fs<br>
                      &gt;&gt;&gt;&gt; rlist        = 1.0        ;
                      short-range neighborlist cutoff (in nm)<br>
                      &gt;&gt;&gt;&gt; rcoulomb    = 1.0        ;
                      short-range electrostatic cutoff (in nm)<br>
                      &gt;&gt;&gt;&gt; rvdw        = 1.0        ;
                      short-range van der Waals cutoff (in nm)<br>
                      &gt;&gt;&gt;&gt; ; Electrostatics<br>
                      &gt;&gt;&gt;&gt; coulombtype    = PME        ;
                      Particle Mesh Ewald for long-range<br>
                      &gt;&gt;&gt;&gt; electrostatics<br>
                      &gt;&gt;&gt;&gt; pme_order    = 4        ; cubic
                      interpolation<br>
                      &gt;&gt;&gt;&gt; fourierspacing    = 0.12        ;
                      grid spacing for FFT<br>
                      &gt;&gt;&gt;&gt; ; Temperature coupling is on<br>
                      &gt;&gt;&gt;&gt; tcoupl        = nose-hoover    ;
                      modified Berendsen thermostat<br>
                      &gt;&gt;&gt;&gt; tc-grps        = Protein
                      Non-Protein    ; two coupling groups - more<br>
                      &gt;&gt;&gt;&gt; accurate<br>
                      &gt;&gt;&gt;&gt; tau_t        = 2.5    2.5    ;
                      time constant, in ps<br>
                      &gt;&gt;&gt;&gt; ref_t        = 300     300    ;
                      reference temperature, one for each<br>
                      &gt;&gt;&gt;&gt; group, in K<br>
                      &gt;&gt;&gt;&gt; ; Pressure coupling is off<br>
                      &gt;&gt;&gt;&gt; pcoupl        = no         ; no
                      pressure coupling in NVT<br>
                      &gt;&gt;&gt;&gt; ; Periodic boundary conditions<br>
                      &gt;&gt;&gt;&gt; pbc        = xyz        ; 3-D PBC<br>
                      &gt;&gt;&gt;&gt; ; Dispersion correction<br>
                      &gt;&gt;&gt;&gt; DispCorr    = EnerPres    ;
                      account for cut-off vdW scheme<br>
                      &gt;&gt;&gt;&gt; ; Velocity generation<br>
                      &gt;&gt;&gt;&gt; gen_vel        = no        ;
                      don¹t assign velocities from Maxwell<br>
                      &gt;&gt;&gt;&gt; distribution<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; --<br>
                      &gt;&gt;&gt;&gt; Gromacs Developers mailing list<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; * Please search the archive at<br>
                      &gt;&gt;&gt;&gt; <a moz-do-not-send="true"
                        href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List"
                        rel="noreferrer" target="_blank">
http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</a><br>
                      &gt;&gt;&gt;&gt; before<br>
                      &gt;&gt;&gt;&gt; posting!<br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; * Can't post? Read <a
                        moz-do-not-send="true"
                        href="http://www.gromacs.org/Support/Mailing_Lists"
                        rel="noreferrer" target="_blank">
                        http://www.gromacs.org/Support/Mailing_Lists</a><br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; * For (un)subscribe requests
                      visit<br>
                      &gt;&gt;&gt;&gt; <a moz-do-not-send="true"
href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers"
                        rel="noreferrer" target="_blank">
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</a><br>
                      &gt;&gt;&gt;&gt;<br>
                      &gt;&gt;&gt;&gt; or send a mail to <a
                        moz-do-not-send="true"
                        href="mailto:gmx-developers-request@gromacs.org"
                        target="_blank">
                        gmx-developers-request@gromacs.org</a>.<br>
                      &gt;&gt;<br>
                      &gt;<br>
                      <br>
                      --<br>
                      Gromacs Developers mailing list<br>
                      <br>
                      * Please search the archive at <a
                        moz-do-not-send="true"
                        href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List"
                        rel="noreferrer" target="_blank">
http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</a>
                      before posting!<br>
                      <br>
                      * Can't post? Read <a moz-do-not-send="true"
                        href="http://www.gromacs.org/Support/Mailing_Lists"
                        rel="noreferrer" target="_blank">
                        http://www.gromacs.org/Support/Mailing_Lists</a><br>
                      <br>
                      * For (un)subscribe requests visit<br>
                      <a moz-do-not-send="true"
href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers"
                        rel="noreferrer" target="_blank">https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</a>
                      or send a mail to
                      <a moz-do-not-send="true"
                        href="mailto:gmx-developers-request@gromacs.org"
                        target="_blank">
                        gmx-developers-request@gromacs.org</a>.<br>
                    </blockquote>
                  </div>
                </div>
              </div>
              <br>
              <fieldset class="mimeAttachmentHeader"></fieldset>
              <br>
            </blockquote>
            <br>
          </div>
        </div>
      </span>
      <br>
      <fieldset class="mimeAttachmentHeader"></fieldset>
      <br>
    </blockquote>
    <br>
  </body>
</html>