<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 <<a
moz-do-not-send="true" href="mailto:b.reuter@uni-kassel.de">b.reuter@uni-kassel.de</a>>
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 <<a moz-do-not-send="true"
href="mailto:b.reuter@uni-kassel.de"
target="_blank">b.reuter@uni-kassel.de</a>>
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>
> I mean a drift of the total energy in NVE -
while with Nose-Hoover the<br>
> drift is in the Conserved-Energy quantity of
g_energy (the total<br>
> energy shows no drift with Noose-Hoover...).<br>
><br>
> Am 15/07/15 um 17:38 schrieb Bernhard:<br>
>> I also did a NVE simulation with the same
parameters, system and<br>
>> starting conditions but with manually set
rlist=1.012nm (since<br>
>> verlet-buffer-drift doesnt work in NVE):<br>
>> There I also got a linear drift (but
smaller) of 0.78 kJ/mol/ps<br>
>> (3.436*10^-5 kJ/mol/ps per atom).<br>
>> For comparison reasons I also did a NVT
Nose-Hoover Simulation with<br>
>> manually set rlist=1.012nm:<br>
>> There I got a comparable linear drift of
0.67 kJ/mol/ps (2.94*10^-5<br>
>> kJ/mol/ps per atom).<br>
>> So no differences between NVE and NVT so
far in my opinion...<br>
>><br>
>> Best,<br>
>> Bernhard<br>
>><br>
>> Am 15/07/15 um 17:10 schrieb Shirts,
Michael R. (mrs5pt):<br>
>>> The conserved quantity in nose-hoover
is not quite as good as the<br>
>>> conserved energy, which should have
no drift at all. For NH, the<br>
>>> conserved quantity should drift as a
random Gaussian process with mean<br>
>>> zero (i.e. go with sqrt(N)). It
shouldn't be drifting linearly.<br>
>>><br>
>>> I would check to see if your system
conserved energy when run with NVE<br>
>>> (use the endpoint of the NPT
simulation). It's easier to diagnose any<br>
>>> problems with an NVE simulation,
which should have virtually no<br>
>>> drift, vs<br>
>>> a NVT simulation, which has random
noise drift. Odds are, if there<br>
>>> is a<br>
>>> problem with the NVT simulation, it
will also show up in the NVE<br>
>>> simulation if only the thermostat is
removed.<br>
>>><br>
>>> 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>
>>> for tests of whether the ensemble
generated is correct.<br>
>>><br>
>>> Best,<br>
>>> ~~~~~~~~~~~~<br>
>>> Michael Shirts<br>
>>> Associate 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>
>>><br>
>>> On 7/15/15, 10:58 AM, "Bernhard" <<a
moz-do-not-send="true"
href="mailto:b.reuter@uni-kassel.de"
target="_blank">b.reuter@uni-kassel.de</a>>
wrote:<br>
>>><br>
>>>> Dear Gromacs Users and
Developers,<br>
>>>><br>
>>>> I have a problem regarding energy
conservation in my 10ns NVT<br>
>>>> protein+water+ions (22765 atoms)
production (minimization and<br>
>>>> equilibration for more than 15ns
was carried out in NPT before)<br>
>>>> simulations using a Nose-Hoover
thermostat (tau=2.5ps).<br>
>>>> On first glance everything looks
fine - the potential, kinetic and<br>
>>>> total<br>
>>>> energy are nearly perfectly
constant (with normal fluctuations) - but<br>
>>>> when I checked the
"Conserved-Energy" quantity that g_energy outputs
I<br>
>>>> had to recognize a significant
(nearly perfectly) linear downward<br>
>>>> drift<br>
>>>> of this "to-be-conserved"
quantity of around 1.7 kJ/mol/ps (7.48*10^-5<br>
>>>> kJ/mol/ps per atom).<br>
>>>> This appears somehow disturbing
to me since I would expect that this<br>
>>>> Conserved-Energy is the conserved
energy of the extended Nose-Hoover<br>
>>>> Hamiltonian - which should by
definition be conserved.<br>
>>>><br>
>>>> If it would be a drift caused by
normal round-off error due to single<br>
>>>> precision I would expect it to
grow with Sqrt(N) and not with N<br>
>>>> (linear)<br>
>>>> (N=number of steps).<br>
>>>> So I would like to know if this
is a normal behaviour and also what<br>
>>>> could cause this (buffer size,
precision, constraints etc)?<br>
>>>> Also I would like to know, if I
am correct with my guess that the<br>
>>>> "Conserved-Energy" quantity is in
this case the energy of the extended<br>
>>>> Nose-Hoover Hamiltonian?<br>
>>>> The .mdp file is atatched (don't
be confused about rlist=1 - since Im<br>
>>>> using the Verlet-scheme the
verlet-buffer-drift option should be by<br>
>>>> default active and determine the
rlist value (Verlet buffer-size)<br>
>>>> automatically).<br>
>>>><br>
>>>> Best regards,<br>
>>>> Bernhard<br>
>>>><br>
>>>> ; Run parameters<br>
>>>> integrator = md ;
leap-frog integrator<br>
>>>> nsteps = 10000000 ;
10000 ps = 10 ns<br>
>>>> dt = 0.001 ; 1 fs<br>
>>>> ; Output control<br>
>>>> nstxout = 5000 ;
save coordinates every ps<br>
>>>> nstvout = 5000 ;
save velocities every ps<br>
>>>> nstxtcout = 1000 ; xtc
compressed trajectory output every ps<br>
>>>> nstenergy = 1000 ; save
energies every ps<br>
>>>> nstlog = 5000 ;
update log file every ps<br>
>>>> ; Bond parameters<br>
>>>> continuation = yes ;
continue from NPT<br>
>>>> constraint_algorithm = lincs ;
holonomic constraints<br>
>>>> constraints = h-bonds ; all
bonds (even heavy atom-H bonds)<br>
>>>> constrained<br>
>>>> lincs_iter = 1 ;
accuracy of LINCS<br>
>>>> lincs_order = 4 ; also
related to accuracy<br>
>>>> ; Neighborsearching<br>
>>>> cutoff-scheme = Verlet ;
Verlet cutoff-scheme instead of<br>
>>>> group-scheme (no charge-groups
used)<br>
>>>> ns_type = grid ;
search neighboring grid cells<br>
>>>> nstlist = 10 ; 10
fs<br>
>>>> rlist = 1.0 ;
short-range neighborlist cutoff (in nm)<br>
>>>> rcoulomb = 1.0 ;
short-range electrostatic cutoff (in nm)<br>
>>>> rvdw = 1.0 ;
short-range van der Waals cutoff (in nm)<br>
>>>> ; Electrostatics<br>
>>>> coulombtype = PME ;
Particle Mesh Ewald for long-range<br>
>>>> electrostatics<br>
>>>> pme_order = 4 ; cubic
interpolation<br>
>>>> fourierspacing = 0.12 ;
grid spacing for FFT<br>
>>>> ; Temperature coupling is on<br>
>>>> tcoupl = nose-hoover ;
modified Berendsen thermostat<br>
>>>> tc-grps = Protein
Non-Protein ; two coupling groups - more<br>
>>>> accurate<br>
>>>> tau_t = 2.5 2.5 ;
time constant, in ps<br>
>>>> ref_t = 300 300 ;
reference temperature, one for each<br>
>>>> group, in K<br>
>>>> ; Pressure coupling is off<br>
>>>> pcoupl = no ; no
pressure coupling in NVT<br>
>>>> ; Periodic boundary conditions<br>
>>>> pbc = xyz ; 3-D PBC<br>
>>>> ; Dispersion correction<br>
>>>> DispCorr = EnerPres ;
account for cut-off vdW scheme<br>
>>>> ; Velocity generation<br>
>>>> gen_vel = no ;
don¹t assign velocities from Maxwell<br>
>>>> distribution<br>
>>>><br>
>>>> --<br>
>>>> Gromacs Developers mailing list<br>
>>>><br>
>>>> * Please search the archive at<br>
>>>> <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>
>>>> before<br>
>>>> 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><br>
>>>><br>
>>>> 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>
>><br>
><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>