…
but I didn’t find that any of the suggestions actually solved the problem.
I followed several suggestions and played around with the settings, but my conclusion is that there is indeed a problem with the combination of fix heat & shake, even for high precision of shake and kspace and a low time step. For example, in an NVE simulation of 2916 SPC/E water molecules at T=400 K, the following settings conserve the energy very well over a time of 2 ns: pppm 1.e-6, dt 0.5, shake 1.e-10 (see first output bellow). Likewise with Ewald summation. I know that I could use a higher time step, but it could then be too high for the non-equilibrium simulation in the presence of a temperature gradient. PLEASE NOTE, that the outputs correspond to two independent runs, meaning that the gradient simulation (second output bellow) starts at a slightly different energy. The reason for that is simply that I equilibrated the input configuration which I actually used for the gradient simulation with a time step of 1 fs, and I wanted to present a comparison for the same time step.
The second output shows that, with the same settings, I lose about 3.5 kcal/mole every 0.5 ns when I use fix heat. The maximum temperature (ave/spatial) in the system for my settings is ~500 K. To rule out a time step issue, I also ran a separate NVE simulation at 500 K without fix heat and the energy was well conserved again. Please find my input script for the gradient simulation bellow. That simulation started with a binary input file (nve.res), which was the result of an NVT simulation (which you can also find bellow) and a subsequent short NVE run. I omitted the script for the NVE run, as it is basically the same as the gradient one but without fix heat and the regions.
I would appreciate further suggestions and comments on how to get rid of this problem.
Are you trying to run NEMD simulations? If yes, try using two thermostats instead for creating the gradient as they are simple and give you right away what the temp of the hot/cold zones would be. If you need to use fit heat for another reason, I pass the baton to the next user as I have no experience dealing with that fix.
Carlos
Thanks for the comment, Carlos. My first attempt was actually to use simple velocity rescaling for the hot/cold regions.
However, the energy fluctuations are quite large, if you do it this way. The fix heat algorithm rescales in a more clever way and handles this very well, whilst conserving the momentum. I cannot really see how you could achieve that with two independent thermostats.
Do you have something particular in mind? If yes, I’d very much appreciate a sketch.
Basically I am interested in creating a nonequilibrium stationary state, so yes NEMD.
Thanks Niall!. That is definitely true for the total momentum. But fix heat conserves the momentum of the individual region during the rescaling step. For fix momentum you cannot specify a region, right?
Let’s assume you can. If you independently thermostat that region and fix its momentum I could imagine that this leads to some competing mechanisms, as fixing the momentum after rescaling, for example, would probably change the temperature again.
However, maybe this is actually not a problem. I just like fix heat because it does everything in one go.
I thought I should maybe highlight a passage from a previous post: http://lammps.sandia.gov/threads/msg30671.html, which suggests that it is not a problem of “fix heat + fix shake” but rather “fix heat” alone:
>> 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.
I guess it all depends on what you are after. I have no expertise measuring the effects of this lack of local momentum conservation in the thermal conductivity. The few times I have worked on the topic the thermostats just worked fine for my systems as I always cross-validated the numbers with RNEMD simulations. Yet, there might be cases where such effect could lead to problems.
You seem to be a very careful tester (a physicist?) and that certainly has helped you to identified other issues with the code. So, I let you keep pushing the developers when raising the topic about a potential bug in fix heat.
Thanks, Carlos! Eventually, I might be interested in using two thermostats instead of “fix_heat” depending on the outcome of this discussion. When you say, use two independent thermostats, you suggest that I should apply the thermostats to different regions in the box, right? I tried something similar to this (for water), but as mentioned before I found that “fix heat” was an improvement with respect to energy conservation - apart from the slight drift. Can you think of improvements to the code snippet bellow?
create regions of low and high temperature and apply thermostats
region Thi_region block 0 INF 0 INF {zlo_Thi} {zhi_Thi}
region Tlo_region block 0 INF 0 INF {zlo_Tlo} {zhi_Tlo}
this takes 3 degrees of freedom for every atom -> 9N_atom/2 dof instead of 6N_atom/2
compute cTlo water temp/region Tlo_region
compute cThi water temp/region Thi_region
rescale temperature to correct for the 3N_atom dof
Thanks, Carlos! Eventually, I might be interested in using two thermostats
instead of "fix_heat" depending on the outcome of this discussion. When you
say, use two independent thermostats, you suggest that I should apply the
thermostats to different regions in the box, right?
Yes.
I tried something similar to this (for water), but as mentioned before I
found that "fix heat" was an improvement with respect to energy
conservation - apart from the slight drift. Can you think of improvements
to the code snippet bellow?
I have never used SHAKE thus have no feeling for how such integrator works.
Have you tried using langevin thermostats instead of temp/rescale?
For computing the heat flux I usually do it in a region that doesn't share
overlap boundaries with the thermostated regions (temperature gradient is
global).
Carlos
PS: I am far from being considered an expert in this type of simulations.
And I mean it, which is obvious once you check the number of papers (ONE) I
have on this topic. I just have no fear to share out loud what I may have
done/encountered
This sounds like a problem similar to swapping momentum of molecules that use shake (like in fix thermal/conductivity). In that case, a solution is to apply the same change in velocity to all the atoms in the frozen molecule. The idea is to change the center of mass velocity but preserve the relative velocities for shake. E.g, see D. Bedrov and G. Smith, J. Chem. Phys. 113, 8080 (2000).
This (the RNEMD behaviour) I have seen as well dealing with the fix/thermal conductivity with bonded potentials. Kind of funny because Muller Plathe had in one of the early papers an explanation on how to do it for molecular systems (swap COM velocities but retain relative ones). I believe the reason why the fix in lammps acts on atoms and not molecules is historical. Mainly because the original developer of the fix
worked on atomic fluids not molecular ones. He may still be around reading theses posts to comment on it.
My solution to mitigate the energy drift was to lower the time step.
I guess I should start looking at the code and calculate the total momentum of the two regions before and after “fix heat” was applied. Maybe there is a tolerance for the error, which can be increased. With regard to the time step, I am trying 0.3 fs at the moment. At this low time step I am losing less energy, about -1.5 kcal/mole as opposed to -3.5 kcal/mole in 0.5 ns. I will look into that.
I checked the fix-heat function and in my opinion it actually does what it should. I tried various things including adding the heat with respect to the velocity of the centre of mass of the whole system, such that both reservoirs use the same frame of reference. It did not help, I always got a drift.
However, I now think that it might be an intrinsic problem related to the fact that the whole system starts to rotate due to the two opposite fix heat commands. I think adding/removing rotational kinetic energy at different locations creates an average net torque causing the box to rotate. For example, if you use “fix momentum” not only to fix the linear momentum but also the angular, you will find a constant rapid decrease in the total energy of the system.
I do not really understand, why you would lose energy this way. Maybe this theory does not even make much sense, but it is as far as I could get with my analysis so far. Has anyone experience with this sort of problem?
I have never dealt with fix heat in my simulations, never mind looking into its source code. The only suggestions that come to mind are: Does the energy drift changes by increasing the system size of your liquid water? Have you tried to see if the same thing happens with an atomic (LJ) fluid?
But I think intrinsically that fix heat and fix shake
are not fully compatible. When you just add random
energy to atoms in a water molecule that is being constrained
by SHAKE, some of the velocity change you induce
will be directed along the bonds, which means the
next iteration of fix shake will take that motion
out via the constraint forces. So I don’t think
the total energy of the system will be conserved, i.e.
fix shake will remove energy. Which isn’t something
it normally does if you just run NVE with conservative
Other people reported loss of energy for LJ too, but I have not looked at it myself. However, I will check this out and also how the rate of energy loss depends on the system size. That would be interesting to see. Thanks for the suggestions (please also read comment to “all”).
@Steve:
Thanks for your comment. Well, in an ideal simulation, the Lagrange multipliers should ensure that there is no velocity component along the connection lines - if you think of a pendulum for example. However, in the simulation there might be a non-vanishing contribution in this direction and I don’t know, whether it can actually explain the effect. To be honest, I doubt it because the problem was also reported for a simple LJ fluid in a previous post, one year ago. Furthermore, the test I did (please see comment bellow addressed to “all”), seems to rule this out as the (main) source of this particular problem.
@All:
The energy remains conserved when I flip the reservoirs with a certain frequency, e.g. every 10th time step. The idea was that this should prevent the box from taking up angular momentum, if that were to be the problem. At the same time rounding errors and errors in combination with SHAKE should keep showing up. What I did was running the same input file + binary input configuration with lammps-25Nov13 and my adapted version of lammps-25Nov13, with just this little modification in fix_heat.cpp with the same number of cores:
void FixHeat::end_of_step() {
if (update->ntimestep % 10 == 0) heat_input *= -1;
I don’t know of any conceptual reason why 2 equal/opposite fix
heat commands should lose energy for an atomic system.
Paul - are you aware of this issue?
However, if I understand your input script
you are doing it for SHAKE molecules and using
regions. That means that there will be water moleclues
which straddle the region boundary. This means
the velocity of one (or two) atoms in the water
will be rescaled, but not the other(s). I believe
that means you have added a component of motion
to the molecule that is along a bond and will be immediately
removed by SHAKE. Hence you didn’t really add as much
energy as you thought. And it could lead to a downward energy
drift.
Assuming a homogeneous density throughout the box, I think that the effect you were mentioning actually cancels out on average. Some molecules lose energy, some gain energy. However, in my setup an inhomogeneous density profile establishes and the density is considerably higher close to the cold reservoir. This might explain the drift in the energy, because there are more molecules effected in the cold region and downscaling might be more likely to remove energy combined with SHAKE.
That also explains perfectly well, why the flipping simulation did not show a drift: because the density remained homogeneous and the effect canceled out again. If this is true, a simple adaption of the fix_heat function would solve the problem: I will just take into account those atoms, where the whole molecule is inside the region. Then there should be no drift! I will get back to you about this tomorrow.
