Temperature was not constant but increased in NVE simulation with TIP4P/Ice model

Dear LAMMPS users,

Hello,
I have a problem with no-constant temperature in NVE simulation with TIP4P/Ice model.

I’m trying to create an ice ball.
First, I gave the initial velocity so that the initial temperature was 10 K.
Using this initial condition (following data file), I performed NVT simulation with 30 K.
And in NVT, the temperature became 60 K within the first 500 fs.
After 3000 fs, the temperature oscillated around 30 K.

Next, I performed NVE simulation.
In this phase, the temperature gradually increased from 30 K to 34 K for 68 ps.
I expected that the temperature would become constant, but it wouldn’t.
The simulation is running, so I can show the temperature evolution until now in the following figure.

I’m not sure what makes the temperature increase in NVE.
I consider the possibility that timestep causes.
But in my simulation, I took timestep with 0.0078125 fs.
I think it is short enough.

Please give me any comments.

data.sep (8.3 MB)
in.water (1.2 KB)
temp

Hello,
Rigid water models are not good at conserving energy, so I am not too surprised by your results.
Usually one need to use a thermostat.
Best,
Simon

This is a ridiculously short timestep for a system with a water potential using fix shake.
You should be able to run with 1fs to 2fs. That should reduce the contribution to the temperature drift per timestep.

Another potential problem is using a cutoff Coulomb. Your forces and energies are still significant at the cutoff distance. LAMMPS offers long-range coulomb solvers that are tip4p compatible.

Thank you for your reply.

Rigid water models are not good at conserving energy, so I am not too surprised by your results.
Usually one need to use a thermostat.

Thank you for your information.
I didn’t know the properties.

Best regards,
Yuki

Here is another thought: with fix nve the conserved quantity is the total energy and not the kinetic energy. Since you are running at very low temperature, have you checked whether your system is fully equilibrated? It may not be and thus the increasing temperature simply an indication of the ongoing equilibration.

Thank you for your comments.

You should be able to run with 1fs to 2fs.

I know that some simulations used this timestep.
But the timestep with 1fs resulted in rapidly increasing temperature in my simulation, so I used a shorter timestep.

have you checked whether your system is fully equilibrated?

I did not have this idea, and I have not confirmed it. I will check. Thank you.

Another potential problem is using a cutoff Coulomb. Your forces and energies are still significant at the cutoff distance. LAMMPS offers long-range coulomb solvers that are tip4p compatible.

This comment is of great interest to me.
About this comment, I have another question.
Is it good to add a long-range Coulombic solver to TIP4P/Ice model?
Or should I use Ewald or PPPM model represented at the bottom of the LAMMPS documentation?
https://docs.lammps.org/Howto_tip4p.html

Best regards,
Yuki

Both Ewald and PPPM are long-ranged Coulombic solvers. :slight_smile:

I’m sorry for my messy questions.

TIP4P/Ice model uses the following parameters (in the LAMMPS document):
O mass = 15.9994, H mass = 1.008
O charge = -1.1794, H charge = 0.5897
r0 of OH bond = 0.9572, θ of HOH angle = 104.52∘
OM distance = 0.1577
LJ ϵ of O-O = 0.21084, LJ σ of O-O = 3.1668
LJ ϵ, σ of OH, HH = 0.0
Coulomb cutoff = 8.5

, and TIP4P with a long-range Coulombic solver:
O mass = 15.9994, H mass = 1.008
O charge = -1.0484, H charge = 0.5242
r0 of OH bond = 0.9572, θ of HOH angle = 104.52∘
OM distance = 0.1250
LJ ϵ of O-O = 0.16275, LJ σ of O-O = 3.16435
LJ ϵ, σ of OH, HH = 0.0

O, H charge, OM distance, and LJ parameters are different between the two models.
My question is
Can we use a long-range Coulombic solver with the TIP4P/Ice parameters?

Best,
Yuki

See kspace_style command — LAMMPS documentation (search it for “tip4p”).

Thank you for the information.
I check that we can customize the parameters.
But, when we use TIP4P/Ice model with a long Coulomb solver, can it represent the ice properties? (such as structure, density, …etc)

I am concerned that it may be able to reproduce liquid water but not be able to reproduce rigid water.

That is a concern that you can address yourself by surveying the published literature of studies that have simulated ice with the TIP4P model and what exact parameters were used and how well the desired properties were represented.

You have to keep in mind that a rigid water potential with simple LJ plus point charges is always going to be a compromise and will never represent all properties equally well. The differences between individual parameter sets tend to be small compared to the differences of either of those models to experimental properties of water at the same conditions. But also keep in mind that experiments can have errors. :wink:

Thank you for your kind comment!

I strongly recognize that I have to study some works.
I will keep in mind your comments.

Thank you!
Best regards