TIP4P energy conservation problem with NVE

Dear lammps users,

My system consists in 280 TIP4P2005 particles with 120 Lennard-Jones particles that are frozen. In the beginning, I obtain the desired temperature (300K) for the TIP4P particles via the Berendsen thermostat, while the LJ particles are frozen by the fix “setforce”. However, when I restart the MD calculation using the NVE time-integrator, the temperature and the energy of the water particles decrease. I was reading in previous discussions that the “fix shake” option can be the crucial point for the energy conservation. However, I haven’t obtained any good result changing the parameters of this fix. My input file is:

TIP4P/2005 water model

echo both
units real # m=g/mol, d=Ang, t=fs, E=Kcal/mol, T=K, P=atm, Dynamic viscoty=Poise, charge=e (negative=-)

Dear lammps users,

My system consists in 280 TIP4P2005 particles with 120 Lennard-Jones particles that are frozen. In the beginning, I obtain the desired temperature (300K) for the TIP4P particles via the Berendsen thermostat, while the LJ particles are frozen by the fix “setforce”.

It is not clear how you can obtain any temperature with the particles frozen.

However, when I restart the MD calculation using the NVE time-integrator, the temperature and the energy of the water particles decrease.

By equipartition temperature is going to drop because some of the kinetic energy converts to potential energy. You should monitor the total energy.

Ray

Dear Ray,

What I do is to fix the Berendsen thermostat just to the water molecules. I use this command: “fix 5 temp/berendsen $T $T 10000” and 10000 time steps are equal to 5 ps. The total energy is decreasing too. What I think is happening is that the fixed particles are absorving energy to the water molecules, and then because of setting all the forces to 0 for the frozen particles, the absorbed energy is lost. I realized that if I decrease the time step to 0.1 fs for example, the energy decreases in a slower way. This would happen if the absorption hypothesis is true. Indeed, if the system is well equilibrated and some of the kinetic energy of the water molecules is being given to the frozen particles as potential energy, then if the energy of the system is conserved, this energy should be given back to the water molecules in a sufficiently long simulation and the mean value of temperature should stay constant, not decreasing.
I have tried to remove the set force command, increase the mass of the frozen particles and integrate all the particles with NVE. However, the “frozen” particles move.

Do you think that the system might be badly equilibrated, i.e., that the “10000” value is too large? Do you think that the energy is being lost via the frozen particles? If yes, do you have any suggestion to fix that?

Thank you!

When you specify a compute to calc temperature, e.g. compute temp,
you specify a group. If you exclude frozen atoms from the group,
you will only get the temp of the mobile atoms.

You can output the value from any compute with the thermo_style
command, or make it the default “temp” via the thermo_modify command.

Steve

Hi Steve,
I do not see how the output of the temperature of the entire system would help me. Indeed, I have done it, and it works as I expected. The temperature of the frozen atoms is 0, and the temperature of the water molecules while doing Berendsen is the correct one (300 K). When I print the temperature of the entire system, the output is less than 300 K, which makes sense because I have some particles that are frozen. However, the problem comes when I do NVE, because the temperature of the water molecules drops down (the temperature of the frozen atoms is still 0, because they are fixed). And sorry, when I was adding before the command for the Berensen I forgot to add the group. The correct command that I use is: “fix 5 water temp/berendsen $T $T 10000"

Jon

I’m confused by the thread now. Can you post a simple
input script (as small a system as possible) that just does
SHAKE with fix temp/berendsen for some waters and
no other contsraints, where the temperature does not
thermostat to what you expect?

Steve

Jon,

The “problem” you are describing is not a problem but an expected outcome. A system with constraints like yours is definitely going to lose energy so that the total energy will descrease in an NVE ensemble.

Ray

Dear Steve and Ray,

I deduce from Steve’s answer that I did not explain the problem well enough, I will try to do it better.

The problem is not related with the thermostat. My system consists on water molecules between two fixed walls made by Lennard-Jones particles. The simulation is divided in 2 parts. In the first part, I equilibrate the system using the Berendsen thermostat just for the water molecules (“fix 5 water temp/berendsen $T $T 10000”) and I set the force to 0 for the LJ particles, thus they are frozen. This way, I get the desired temperature for the water molecules (300 K) and in this part of the simulation everything works fine.

After the equilibration, what I do is to unfix the Berendsen thermostat and apply the NVE time-integrator for the water molecules, but still maintain the set force equal to 0 for the LJ particles. The problem that I find here, or at least is unexpected for me, is that the temperature of the water molecules starts decreasing. The system is clearly loosing energy, and my analysis suggests that not only the energy is passing from kinetic energy of the water molecules to potential energy, but also is being lost from the point of view of the entire system. My guess is that the kinetic energy that the LJ particles get in each time step is removed by the “set force” command, and therefore, this energy is being lost. However, I do not have a clear solution for this problem. I though to increase the mass of the LJ particles and add hem in the NVE integrator, but this did not work.

Ray, this seems to be expectable for you. However, I was thinking on a simple example of two particles interacting between them where the position of one of them is fixed. I would expect in this case the energy to be constant, why shouldn’t it be in my case? Do you think that the problem comes from the fact that i’m not adding the LJ particles in the NVE integrator?

I was thinking on adding a very strong harmonic potential to each LJ particle to maintain their position within the walls and add them in the NVE time-integrator. When I get the answer I will let you know.

Thanks in advance,

Jon

If the wall LJ particles are simply frozen (or not integrated), then
they are like an external field on the water molecules. They should

not affect the energy of the remainder of the system.

Once you start NVE, does the T of the water

decrease indefinitely, or does it come to a new equilibrium.

Steve