Dear Lammps users,
I have two queries. They both relate to my attempts to perform non-equilibrium simulations of spc/e water under a temperature gradient via the “fix heat” command.
The first question is slightly more general and is an observation that the heat fix appears to progressively remove energy from the system. I have performed various simulations using the “fix heat” command (For a Lennard-Jones, Dumbbell and Spc/e water) and in all cases the total energy of the system appears to decrease over time. The extent to which it decreases from preliminary tests appears to be related to the rate of energy addition/removal.
Obviously in all cases I am supplementing each addition of heat, with the same magnitude of removal at a separate location. I was fairly convinced that I have been implementing the fixes correctly and was wondering if this could be an effect of variable accuracy/rounding error upon the distribution of the added/removed momenta to the system. This is just an extreme rookie’s guess so apologies if it is a misguided one.
Secondly, and of more importance to myself I have recently noticed that my spc/e simulations are showing zero for the intra-molecular energy component in the log file. I was under the impression my data file and input script for spc/e water was correct and it was infact only my investigation of the aforementioned non-equilibrium simulations energy drift which lead me to notice the zero molecular contribution.
I have attached the simplest input script and data file which reproduces my issue with the spc/e intra-molecular energy.
Thank you in advance, Jeff.
in.spce (618 Bytes)
data.spce (414 KB)
I can’t really comment on “fix heat.” However, the reason why you don’t see an internal energy for SPC/E water is a result of the SHAKE algorithm. Since the algorithm forces the bond lengths and angles to take fixed values, the values given for the bond and angle strengths are immaterial. You get the same results for the geometry of the system, independent of those values. As a result, it doesn’t make sense to report a value for those energies, because there is no “scale” to measure them on.
A system with a mix of SHAKEn and unSHAKEn bonds would have an intramolecular contribution to the internal energy. SPC/E water, where all bonds are already SHAKEn, would not.
I can't really comment on "fix heat." However, the reason why you don't see
an internal energy for SPC/E water is a result of the SHAKE algorithm. Since
the algorithm forces the bond lengths and angles to take fixed values, the
values given for the bond and angle strengths are immaterial. You get the
even more so. bonds and angles that are "shaken" are removed from
their respective bond/angle lists and thus never computed.
Paul can likely comment on using fix heat with SHAKE-constrained
I haven't tried using fix heat with fix shake. Fix heat scales velocities in order to add/remove energy, and fix shake imposes bond/angle constraints by adding constraint forces. As you might imagine, there could be some interplay between the two, with possibly unintended consequences. But I think that they might work OK together if the rate of heat add/removal is gradual enough. I'd recommend testing to make sure you get your specified fix heat energy addition/removal rate by tracking your system's energy. You should probably start with a non-shake system with only one fix heat and then add layers of complexity to make sure the simulations are performing as you expect. If you find that fix heat and fix shake don't work together, please let us know so that we can add a warning message.
I noticed from your example that you're using fix nvt. I realize that your example is not a fix heat example. But I wanted to make sure and stress the fix heat documentation statement: "This fix should not normally be used on atoms that have their temperature controlled by another fix - e.g. fix nvt or fix langevin fix."