<html>
  <head>
    <meta content="text/html; charset=windows-1252"
      http-equiv="Content-Type">
  </head>
  <body bgcolor="#FFFFFF" text="#000000">
    I see... I didn't consider the implicit constraints in rigid water
    models like TIP3P (which I used). My bad! Thx for pointing that out,
    Mark!<br>
    <br>
    I rerun the simulations (NVT, Nose-Hoover) in double precision and
    still get perfect linear drifts of the Conserved Energy - but
    upwards and with different reaction to variations of rlist.<br>
    To compare (rlist= auto means verlet-buffer-drift=0.005; all
    energies in kJ/mol/ps per atom):<br>
    rlist        Conserved-En (SP)    Conserved-En (DP)<br>
    auto        -7.14*10^-5            4.31*10^-5<br>
    1.012      -2.94*10^-5            8.82*10^-5<br>
    1.02        -7.94*10^-5            4.27*10^-5<br>
    1.05        -1.15*10^-4            4.93*10^-7<br>
    <br>
    Interesting is the reversion of the slope for DP compared to SP.<br>
    Even more interesting is the two orders of magnitude smaller drift
    of 4.93*10^-7 for rlist=1.05 for DP compared to the average (even
    three orders smaller compared to the SP value of -1.15*10^-4 for
    rlist=1.05).<br>
    Since the drift of the conserved energy for DP and rlist=1.05 is
    only 5*10^-4 % of the average conserved energy in 100ps simulation
    time, one could even say that the conserved energy is constant in
    this case.<br>
    While in the SP case the drift for rlist=1.05 is the biggest and
    quite significant.<br>
    Seems like the nearly perfect conservation for DP and rlist=1.05 is
    something like a lucky cancellation of errors?<br>
    <br>
    My conclusions of this are at least:<br>
    - The drift isn't resulting from a round-off error (since DP doesn't
    improve the behaviour significantly).<br>
    - It may be an algorithmic error (since it goes linear with the
    number of steps N) that is at least influenced by the    size of the
    verlet buffer (rlist).<br>
    - Its not related to the thermostat since NVE doesn't change
    anything.<br>
    <br>
    But were is the reversion of the slope for DP compared to SP coming
    from? <br>
    Can this really be related to the constraints (mainly of the water
    since switching of the constraints for the protein doesn't change
    anything)?<br>
    <br>
    At last a simple technical question: Were can I find the value of
    rlist if verlet-buffer-drift is switched on and rlist is calculated
    automatically?<br>
    <br>
    Best,<br>
    Bernhard<br>
    <br>
    <div class="moz-cite-prefix">Am 16/07/15 um 08:10 schrieb Mark
      Abraham:<br>
    </div>
    <blockquote
cite="mid:CAMNuMAT-nK+BorLv1KMdT55C7H3Yw8ywMnJD_589p0Jq2v+P6Q@mail.gmail.com"
      type="cite">
      <p dir="ltr">Hi,</p>
      <p dir="ltr">Likely your water still used SETTLE, as the mdp
        option merely disables the automatic conversion of normal bonds
        to constraints. </p>
      <p dir="ltr">Mark </p>
      <br>
      <div class="gmail_quote">
        <div dir="ltr">On Wed, 15 Jul 2015 19:12 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">
          <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>Am 15/07/15 um 18:14 schrieb Mark Abraham:<br>
            </div>
          </div>
          <div bgcolor="#FFFFFF" text="#000000">
            <blockquote 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"
                    target="_blank">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"
                        target="_blank">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></fieldset>
              <br>
            </blockquote>
            <br>
          </div>
          --<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>.</blockquote>
      </div>
      <br>
      <fieldset class="mimeAttachmentHeader"></fieldset>
      <br>
    </blockquote>
    <br>
  </body>
</html>