I am calculating the dynamics of TIP4P/2005 flexible water, as
described in the input file below. The group of water is annealed from
1K to 300K with NVT thermostat.
However, the system equilibiates at 316.7 K, instead of the desired value 300K.
I used a ICE 1C structure with 624 H2O's. The averaged KE of 1E6 to
2E6 steps (1 ns to 2 ns) is 610.7731 Kcal/mol.
Hand calculation to this value seems to validate the temperature
calculation algorithm (DOF*0.5Kb*T=KE, where DOF=3*642-3 in this
case).
So I am guessing some incompatibility may occur between my TIP4P/2005f
with NVT algorithm?
Your help/comments/enlightenment is greatly appreciated.
Thank you so much for the enlightenment. I decreased the time step to 0.25 fs and ran for an annealing process from 1 K to 300 K for 500 ps, then thermalize it at 300K for 3.75 ns. Yet the resulting averaged temperature is still not 300 K but 316.6K…is there any mistake I overlooked?
Thank you so much for the enlightenment. I decreased the time step to 0.25
fs and ran for an annealing process from 1 K to 300 K for 500 ps, then
thermalize it at 300K for 3.75 ns. Yet the resulting averaged temperature is
still not 300 K but 316.6K....is there any mistake I overlooked?
perhaps. in the quoted input, you would define a new nvt fix
without unfixing the old one.
also, your formula for the nDOF was off, it was accounting
for the fact that waters have 3 atoms.
btw: what is the temperature that LAMMPS prints?
is it the same that you compute, or the desired 300K?
perhaps. in the quoted input, you would define a new nvt fix
without unfixing the old one.
Let me emphasize that this is a fatal mistake.
If you do not unfix the first thermostat, you are running
with 2 thermostats. LAMMPS should have printed
a warning that you are integrating atoms twice if
this is the case.
Also, if you let LAMMPS compute the temperature of the
TIP4P waters, it will get the number of DOF correct.
Thank you so much for the valuable inputs. I tried various alternatives in the setting but the resulting temperature still deviates from my desired value, namely, I want NVT at 300 K but the system hoovers at 316 K, all are shown in the log file.
Additionally, could you give me some hints in the source code so that I can look into it when running? fix_nvt.cpp does not seem to be an entry point to thermostatting…what function/subroutine does it call?
Great appreciations again for the valuable time/effort.
Thank you so much for the valuable inputs. I tried various alternatives in
the setting but the resulting temperature still deviates from my desired
value, namely, I want NVT at 300 K but the system hoovers at 316 K, all are
shown in the log file.
i set up a calculation using your potential
parameters starting from a pre-equilibrated
water structure using only one integrator/thermostat
and the temperature is fine.
all indications and experience hint at that
you are making some error in your input.
have a closer look at your output.
please provide a complete input deck that
can be easily and quickly run by somebody
else that demonstrates the aberration.
Additionally, could you give me some hints in the source code so that I can
look into it when running? fix_nvt.cpp does not seem to be an entry point to
thermostatting......what function/subroutine does it call?
??? what do you mean as "entry point"?
most certainly that class is the "entry".
it will then call methods in its base class
with the suitable settings to use the generalize
nose-hoover facilities in there as needed.
The temperature control works fine when I followed your suggestion by heating up a rigid model to the desired temperature. Then the flexible setting is turned on (i.e., enables the SHAKE algorithm).
However the averaged bond angle and length are not per literature…I am still sorting things out. Thanks again for the help!
The temperature control works fine when I followed your suggestion by
heating up a rigid model to the desired temperature. Then the flexible
setting is turned on (i.e., enables the SHAKE algorithm).
sorry, but enabling SHAKE makes your model *rigid*.
if your input is correct, your simulation should work correctly
regardless of whether you have a pre-equilibrated input data.
i only mentioned that, to indicated that i started from a different
input than you. incidentally, the input i used was equilibrated
with a different water model, so it needed equilibration as any
other input, only less.
However the averaged bond angle and length are not per literature....I am
still sorting things out. Thanks again for the help!
well, if you use shake, they cannot be like it. and if you don't
have shake, then you may have made a mistake converting
the force constants. there are many different conventions.
you have to compare the formula and overall settings very carefully.
great attention to detail is almost always required.