<div dir="ltr">[Disclaimer: tl; dr the entire discussion.]<div class="gmail_extra"><br></div><div class="gmail_extra">The drift profile plot (Fig 8) and related discussion in the paper Mark linked should cover most if not all your questions. Both constraints effect and cancellation of errors is explained there.</div><div class="gmail_extra"><br clear="all"><div><div class="gmail_signature">--<br>Szilárd</div></div>
<br><div class="gmail_quote">On Thu, Jul 16, 2015 at 11:25 AM, Bernhard <span dir="ltr">&lt;<a href="mailto:b.reuter@uni-kassel.de" target="_blank">b.reuter@uni-kassel.de</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
  
    
  
  <div bgcolor="#FFFFFF" text="#000000">
    I see... I didn&#39;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&#39;t resulting from a round-off error (since DP doesn&#39;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&#39;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&#39;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>Am 16/07/15 um 08:10 schrieb Mark
      Abraham:<br>
    </div><div><div class="h5">
    <blockquote 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 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">
          <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 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 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&#39;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&#39;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 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 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, &quot;Bernhard&quot; &lt;<a 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
                      &quot;Conserved-Energy&quot; 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 &quot;to-be-conserved&quot;
                      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; &quot;Conserved-Energy&quot; 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&#39;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 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&#39;t post? Read <a 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 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 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 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&#39;t post? Read <a 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 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 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 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&#39;t post? Read <a 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 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 href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@gromacs.org</a>.</blockquote>
      </div>
      <br>
      <fieldset></fieldset>
      <br>
    </blockquote>
    <br>
  </div></div></div>

<br>--<br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a 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&#39;t post? Read <a 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 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 href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br></blockquote></div><br></div></div>