salt water energy instability

Hello lammps-users,
I am running two analogous simulations of a salt water solution in a nanochannel. The channel is periodic in x and y, and bounded in z with a Lennard-Jones type wall (personally modified 9-3 wall to a 10-4-3 wall). The simulations include + and - ions. The interatomic potentials are all lj/cut/coul/long, where long-range electrostatics are computed using PPPM in slab mode. An additional background electric field is applied (fix efield) to replicate a charged parallel plate capacitor.
The difference between the two simulations is the solvent. In one case, the solvent is a TIP3P water model with explicit H2O and partial charges. The shake algorithm is used for the bond vibrations. The dielectric constant is 1. In the second case, the solvent is a Lennard-Jones particle with equivalent mass to the H2O molecule and the same Lennard-Jones parameters as the TIP3P oxygen. To account for dielectric effects of the polar solvent, a dielectric constant of 80 is used. I have made an effort to keep everything else the same/analogous.
Both cases equilibrate to 300K using NVT. After the systems are equilibrated for a long time (at least 2 million timesteps), I switch the NVT to NVE. In the non-polar solvent case, the temperature remains steady at 300K. In the polar solvent case, the temperature (and total energy) rises.
I suspect the shake algorithm is the culprit for the lack of energy conservation. Can anyone confirm or dispute this? Does anyone have experience running shake without a thermostat? Thank you in advance.

Jonathan

You could try 2 things, independently. Crank up the SHAKE tolerance.
Run with flexible TIP3P water, which doesn't use SHAKE.

Note that if you are running non-equilibrium MD, e.g. flow thru the channel,
then the fix efield adds energy to the system.

Steve

Hello lammps-users,
I am running two analogous simulations of a salt water solution in a
nanochannel. The channel is periodic in x and y, and bounded in z with a
Lennard-Jones type wall (personally modified 9-3 wall to a 10-4-3 wall).
The simulations include + and - ions. The interatomic potentials are all
lj/cut/coul/long, where long-range electrostatics are computed using PPPM in
slab mode. An additional background electric field is applied (fix efield)
to replicate a charged parallel plate capacitor.
The difference between the two simulations is the solvent. In one case,
the solvent is a TIP3P water model with explicit H2O and partial charges.
The shake algorithm is used for the bond vibrations. The dielectric
constant is 1. In the second case, the solvent is a Lennard-Jones particle
with equivalent mass to the H2O molecule and the same Lennard-Jones
parameters as the TIP3P oxygen. To account for dielectric effects of the
polar solvent, a dielectric constant of 80 is used. I have made an effort
to keep everything else the same/analogous.
Both cases equilibrate to 300K using NVT. After the systems are
equilibrated for a long time (at least 2 million timesteps), I switch the
NVT to NVE. In the non-polar solvent case, the temperature remains steady
at 300K. In the polar solvent case, the temperature (and total energy)
rises.
I suspect the shake algorithm is the culprit for the lack of energy
conservation. Can anyone confirm or dispute this? Does anyone have
experience running shake without a thermostat? Thank you in advance.

no. the reason for the increase in energy is the use of fix efield.
this is adding a force to all charged particles and thus energy.

axel.

Thanks for the suggestions. Actually my efield is applied in the wall-normal direction, so the ions serve to shield the background efield (i.e., no net force added at steady state). There is no flow in the steady state, and furthermore, in the nonpolar solvent case, there is no change in energy/temperature. I will try playing with the shake tolerance. Thanks again!

Jonathan