[lammps-users] water temperature rising

Hello all,
I have a two-fold question. I am simulating thermal transport in a water-based system. For my control simulation, basically I have a water box (~ 30 Angs x 30 Angs x 100 Angs) and after I equilibrate it to some target temperature (using fix nvt to say, 300K), I use fix heat on both ends in the long direction with fix nve to get a temperature gradient. The boundaries are p p f, with the z-direction bounded by LJ-93 walls. I am using the TIP3 model, with pppm in slab mode for the long range electrostatics. Here are my questions/problems:

  1. I have noticed that I need a timestep of 0.05 femtoseconds or smaller, otherwise the temperature will rise when the thermostat is turned off. I guess this has to do with integration error in the Verlet algorithm. Can anybody confirm this, particularly anybody who has simulated water systems?

  2. So now, I equilibrate my system with a timestep of 0.5 femtoseconds, but when I switch to nve, I decrease the timestep to 0.05 femtoseconds. This seemingly works ok at first. I run the simulation in nve for 10 million timesteps, and the temperature remains around the target temperature. But I intend to disregard these steps as transient, and run the simulation for as long as possible in “steady state.” So following the 10 million steps, I take the restart file and continue it for another 10 million steps. But the temperature immediately begins to rise after the restart. Within 1 million time steps, the temperature rises 2 orders of magnitude. Eventually the simulation crashes. I checked that I am still using the same conditions as before the restart, so I’m not sure where the problem is. Has anybody experienced anything like this before?

Thanks!

Jonathan

Sounds to me like you’ve got an energy conservation problem. You shouldn’t have to use such a small timestep size. If you are using shake for your waters, you should be able to get away with up to 2 fs timesteps. The problem isn’t the Verlet algorithm.

With such a fast rise in temperature, you’ve probably got an error somewhere in your force-field definition, or perhaps with how you’ve got your walls set up. I’d recommend trying to track down the culprit by eliminating the fancy stuff (fix heat, the walls, slab mode, PPPM, etc), then running vanilla NVE with on a well equilibrated system and a reasonable timestep size and see if you get energy conservation. Then you could put the fancy stuff back in a piece at a time and figure out what is causing the lack of energy conservation.

Paul

Hi Paul,
Thanks for your reply. Regarding the timestep issue, I’ve been using shake. I had been playing around with the timestep for a while, and ultimately I found that the small timestep was the solution.
I am trying your recommendation, but one thing I can’t remove is the PPPM in slab mode. I have observed that without the long range electrostatics, the water box will “freeze”, i.e. the atoms seem to get stuck and vibrate in place. But as long range electrostatics are included, this problem goes away.

Jonathan

You should be able to remove both PPPM and slab mode for testing purposes. You could try a long cutoff with smoothing and it should work. If not, there’s probably a bigger problem with your force field. But if you can get it working (good energy conservation with a reasonable timestep size) with PPPM and slab mode, you may not need to perform test runs without them.

A good exercise would be to try running a simple 3D periodic box of water with shake and a 2 fs timestep, NVE, and no long-range electrostatics. You should see quite good energy conservation.

A really small timestep is not a very good solution — it tends to just mask the real underlying problem that you are still trying to find. Probably a problem with your force field setup on perhaps your initial configuration.

Paul

I see. Ok sounds good, I will try it. Thanks!

Jonathan

Hi Paul,
I ran a series of tests with a 2 femtosecond timestep, but it seems consistent with what I saw before. Here is a list of things I varied.

Boundary: fully perioidic vs. fixed in z-direction (with LJ walls)
Coulombics: “lj/charmm/coul/charmm 10.0 13.0” vs. “lj/charmm/coul/long 10.0 13.0” with pppm
Heat flux: off vs. on

All the cases I ran with pair style “lj/charmm/coul/charmm 10.0 13.0” gave me the “freezing” effect that I saw before. All the cases with pppm gave me an exponential increase in temperature. It’s uncertain whether the boundary and heat flux options played a role.

I wasn’t able to re-create the restart problem (where it is stable before the restart but unstable after). I believe the only differences are domain size and timestep.

Jonathan

Hi Jonathan. I’d recommend going back to running a simple water box and trying to achieve good energy conservation without freezing. Please see the attached for an example input and data file.

Paul

water_equil.data (421 KB)

water.in (667 Bytes)

Hi all,
An update on this issue. I am using a development version of LAMMPS, so it appears I didn’t have all the latest bug fixes in place. When I used the released version, it seems all the problems go away. I will update later if I notice anything.

Jonathan