Problem of energy conservation using fix heat command

Dear lammps users,

My problem is that I set up a model of spc/e water boiling on Cu plate. In the simulation, I use “fix Shake all shake 0.0001 20 0 b 1 a 1” to make water rigid. I use “fix heat heat 1 0.0082” command to input heat flux to heat reservoir, then I account the energy change of water to calculate dE/dt( energy change divide time change), the value is around 0.13eV/ps, which is much larger than 0.0082ev/ps( around 16 times) that I input into the heat reservoir. There is only one heat flux input, so this result is against with energy conservation law.

Then I use “fix heat heat 1 0” command to simulate, but I get almost the same results as “fix heat heat 1 0.0082”, the energy change of water(dE/dt) is about 0.13eV/ps. I don’t know where the problem is.

During the simulation, the units is metal. It first in NVT to run 100000 timesteps, then in NVE to run 2500000 timesteps. The timestep is 0.001ps. I only set the heat reservoir, there is no cold reservoir.

I think for a long time but I can not solve this problem, so I need your help, could you please help me to solve this problem?


Fix nve doesn’t apply any “energy conservation law” to your simulations. It seems to me your model (structure, potential, and simulation settings) doesn’t conserve energy. Note that “fix heat” is just a diagnostic tool and does not affect dynamics at all. So if your model does not conserve energy, you should examine which part of your model needs to be improved.

Sorry – minor correction: fix heat does add heat (kinetic energy) to the system so it alters dynamics. Remove fix heat and see if you can conserve energy. If yes, add fix heat and examine your fix heat command.