[EXTERNAL] Re: Energy conservation problem with fix heat + fix shake (SPC/E water)

You also might try using a smaller heat flux, every timestep (Nevery=1), and/or a larger system size in order to minimize the magnitude of the velocity rescaling that fix heat does.

Paul

Paul, thanks for all the suggestion.

With regard to the settings, I did about 12 preceding simulations over 2 ns for various temperatures, time steps, kspace solvers and shake settings (see initial posts). I am very confident that the settings I use (pppm 1.e-6) conserve the energy perfectly well (within 2 ns) at the highest temperature that builds up during the simulations. (I do know that pppm is not supposed to conserve the energy exactly, but the energy you lose because of that is negligible to the problem I am facing.) The spce/e configurations were generated from a lattice structure, followed by a 100 ps nvt run and a subsequent 200 ps nve equilibration. Not very long, but should be sufficient. I will follow all your suggestions, but I decided to start to look into what Steve was saying about SHAKE removing some energy, followed by looking at LJ.

I found that Steve was right to some extent. I modified various methods (group->vcm, group->ke and group->mass and FixHeat->end_of_step) to take into account an additional mask, so that only atoms are considered for rescaling when the whole molecule is in the region. That lowers the energy drift, but does not completely solve the problem (see “output water_mod”, “output water_orig”). Next I performed LJ simulations to find the remaining problem. Again there is a drift, but I have to admit that I am not a 100% confident about my settings, as I use LAMMPS with LJ for the first time. I took the “melt” example, added the regions and decreased the time step to 0.002. Please have a look at the output with fix heat (output lj_grad) and without fix heat (output lj_simple) and the attached input file lj.lmp. As you can see there is again a problem, which slowly builds up during the simulations.

I think there is actually an easy explanation to all this:

For the simple LJ case as well as for my modified water + shake, each fix works perfectly well! Momentum is conserved and +heat is added to that region. However, both fixes (for two different regions) use a different frame of reference. Kinetic energy depends on the frame of reference. The regions might move with different speeds with respect to a common frame of reference, e.g. that one where the v_com of the total system is 0, and therefore you will, overall, never add/remove the exactly same amount of energy, unless they move with the same speed. In my opinion this is an intrinsic problem: region 1 needs to know how fast the centre of mass of region 2 moves in order to resolve the problem. Alternatively, one could use the v_com of the entire box as common frame of reference. However, my attempts to do that have not been successful yet. What do you think about this theory?

Peter

lj.lmp (1.8 KB)

Maybe I am wrong about the necessity of a common frame of reference. I played around with the code and could not get any improvement.

However, another thing that I found which might be relevant is that the box picks up angular momentum. For example, if I apply “fix 1 all momentum 1 linear 1 1 1 angular” this constantly takes out energy from the system. If I am not mistaken this means that the box starts rotating due to the setup and the way fix heat rescales the velocities. This is not really surprising, but could you imagine this leading to a problem with energy conservation? Maybe in combination with periodic boundary conditions? I have no experience with rotating boxes.

I rarely use fix momentum and I would generally advise against it.

I think it could mess with energy conservation, independent of fix heat.
It’s more designed for extreme cases where the system is going to

acquire a large drift. In your system, if it is periodic, I suppose it could
drift a bit due to the 2 energy changes, however I imagine it would
be minor, and you could monitor energy conservation by correcting
for it with compute temp/com on the entire system.

Steve

Steve, thanks for the hint.

Unfortunately, omitting the “fix momentum 1 linear 1 1 1” does not improve the situation, but it is good to be aware of that.
I cannot further look into this problem anymore, but I hope that the list will be kept updated, in case someone finds a solution to/explanation for this problem.

Thanks a lot for all your help!
Peter

Dear LAMMPS users,

I carried out some energy conservation studies for a LJ system containing two regions that are thermostated with fix_region. As reported in previous emails, the energy is not conserved for SPC/E water, but not even for LJ.

I carried out several more runs for a LJ system of 1728 atoms with a temperature gradient along the z-direction. To be sure that fix_momentum does not interfere with energy conservation, I removed the command for one set of runs. Also I carried out NVE simulations without a gradient at the highest kinetic energy that appears in any of the runs to show that the energy is well conserved with my settings.

It really seems that there is a bug in the way fix_heat works. Maybe this is a naive question, but going back to my “flipping experiment”, where I exchanged the hot/cold reservoirs every other timestep and did not see a loss of energy: Would it be possible that by alternating the order of fix_heat (flipping hot<->cold), I also changed the numerical order of the splitting procedure to a second order method (e.g. Strang-Splitting at a multiple of the small dt)? I mean, maybe the integrator (2-nd order in dt, velocity) is coupled with fix_heat (1st-order in dt) with a first-order splitting procedure and the overall scheme is unstable.

Please find attached the plots of a comparison with “fix_momentum” and one without “fix_momentum” for various fluxes, as well as a sample input script for generating the data (lj.lmp). The full archive with all the results, plot scripts, etc… can be downloaded from here:

https://dl.dropboxusercontent.com/u/3692937/lj_1728.tar

I would be happy to carry out more simulations, if someone of the developer team is willing to help me look into this.

Best regards,

Peter

w_fixmom.png

wo_fixmom.png

lj.lmp (2.12 KB)

Thanks for the careful analysis. I think Paul is the person to look into this, as he did
the initial development of fix heat.

Steve

Thanks, Steve!

Yes, thanks Peter for your careful analysis of this, and for sharing your results. Your results make a compelling argument that imposing a couple of equal and opposite fix heat commands can lead to more energy drift than seen in a vanilla NVE simulation. This is counter to what I’d thought before, but perhaps shouldn’t be all that surprising since fix heat does involve velocity rescaling, thereby slightly messing with the time integration and energy conservation. However, I do think that the 0.1% worst case drift in energy you show after a million timesteps is still relatively small, and probably acceptable for many users’ applications. I don’t have a good feel for how aggressive your thermal gradient and flux are relative to what the typical user would use, but you’ve clearly shown that the more aggressive you pull/push energy into the system, the more energy drift you see.

I think that fix heat is implemented correctly in LAMMPS in that it follows the recipe outlined in section 6.6 here:

SANDIA REPORT

SAND2004-4778

Unlimited Release

Printed September 2004

A Robust, Coupled Approach for Atomistic-Continuum Simulation

S. Aubry, D.J. Bammann, J.J. Hoyt, R.E. Jones, C.J. Kimmer, P.A. Klein, G.J. Wagner, E.B. Webb III, J.A. Zimmerman

(Let me know if you need me to send you the report.)

Paul

Thanks a lot for your message, Paul.

I had a look at the relevant section in the Appendix of the report, which explains the algorithm, and from what I can see and from previously testing the code, I completely agree that the FixHeat::end_of_step() procedure is correct! However, do they also talk about how the procedure needs to be coupled to the integrator? I have the feeling, that although your implementation is correct, the coupling scheme might not be. Velocity rescaling is closely related to the Gauss-Thermostat which imposes dT/dt = 0. Therefore I think of fix_heat more like solving a first order differential equation that is coupled to a second-order differential equation (Newton’s law). Now talking about a different notion of order, problems can arise if a) the order of the algorithms solving the individual equations do not match b) they do match, but the coupling procedure has the wrong order (e.g. exp(dt (A+B)) = exp(dt A) exp(dt B) + O(dt) instead of exp(dt (A+B)) = exp(dt/2 A) exp(dt B) exp(dt/2 A)+ O(dt^2), where operators A and B evolve the relevant quantities in time, respectively.)

I therefore think it might be worth (and quite cheap) trying to call the same fix_heat procedure twice (once before the integration step and once afterwards) with dt/2, as this would be a higher order coupling of the two individual differential equations.

With regard to the magnitude of the thermal gradient and the flux, for the LJ-system I used somewhat arbitrary values. However in the simulation of SPC/E water, you lose an energy equivalent of about 0.5 K/ns for a system of 4096 molecules with a timestep of 0.5 fs. This is for a gradient of about 5 K/Amstrong, which is quite strong. I admit that this is a rather extreme situation. However, even if you lose only half the energy simulating a moderate gradient, it would still lead to a temperature drift of 10 K in 40 ns which I wouldn’t dare to neglect.

A few weeks ago, Steve suggested that rescaling the velocities of molecules that partly reach out of the region introduces some velocity components along the bonds, which will be subsequently removed by SHAKE. That turned out to be true and I tested that by restricting the rescaling-algorithm to molecules that are fully contained within the region. With this you can considerably reduce the energy loss. The rest is based on a different problem - the “lj-problem”. I don’t have an exact estimate of how much the molecules region out of the box contribute to the total energy loss, but I think it is probably between 20-40%.

Do you think it would be worth trying a different splitting scheme? Would LAMMPS easily allow for that?

Best, Peter

LAMMPS does allow for a different splitting scheme as you suggest, and I think it wouldn’t be all that hard to do. If you’re willing to give it a try, I would indeed be interested in what you find.

Paul

Hi Paul/Steve,

I am afraid the splitting scheme does not seem to be the problem. Especially, because the individual operators for the rescaling steps (warm/cold) seem to commute as the regions are non-overlapping and they are also linear in dt, respectively. Do you think a simulation of ideal gas (if possibly in LAMMPS) might help?

There must be some bug, as the algorithm should just conserve the energy sharply (w.r.t. floating point operations). Is it likely that you or somebody else from the developer team will look into this or is it up to me?

Obviously, I am very interesting in knowing what is going on, but I think finding the problem now could also avoid future posts of this kind and avoid troubles.

Thanks for your help up to this point,

Peter

My suspicion is that it is inherently more difficult to perfectly conserve energy when doing velocity rescaling as we’re trying to do with equal and opposite fix heat commands and a fix nve command. I think that fix heat has been correctly implemented as per the formula in the SAND report I referenced, but would welcome you or others verifying that that is true. Your ideal gas idea should be doable with LAMMPS and might help shed some light, but I think I’d better leave these verification exercises up to you or others who suspect there may be a bug in the fix heat implementation.

Paul

Hi Paul/Steve,

sorry for my late reply. I have been working on the problem for a few days
now and I think I am getting closer to fixing it. For test purposes I would
need to update the forces within the "fix_heat" method as I am correcting
the coordinates. If I don't do that, the initial_integration with Verlet
uses wrong_forces.

I had a look at verlet.cpp and checked some of the flags to determine, what
force-methods I need to call for the update. I think for my LJ-System
something like:

  //force_clear();
  int nall = atom->nlocal + atom->nghost;
  for (int i = atom->nlocal; i < nall; i++) {
    f[i][0] = 0.0;
    f[i][1] = 0.0;
    f[i][2] = 0.0;
  }
  force->pair->compute(1,1);
  timer->stamp(TIME_PAIR);
  comm->reverse_comm();
  timer->stamp(TIME_COMM);

should do it. However, it doesn't seem to work like that. Could you please
advice me on how to update the forces at the end of fix_heat, such that the
next Verlet update can take into account modifications of the coordinates?

Many thanks,

Peter

I’m not sure what you’re asking. Fix heat only
rescales velocities. Those are done at the end of the
timestep, so fix heat provides an end_of_step() method.

If you also want to modify forces, after pairwise contributions,
then provide a post_force() method, like several other
fixes do.

Paul can comment as to whether that makes sense. I don’t
see why fix heat should need to do that.

Steve

Dear LAMMPS users and developers,

We recently proposed an extension to the heat exchange (HEX) algorithm (Ikeshoji and Hafskjold 1994), which I would like to make publicly available to the LAMMPS community. Our new enhanced HEX (eHEX) algorithm resolves the energy conservation problem of the original algorithm, which we discussed in this thread quite some time ago (http://sourceforge.net/p/lammps/mailman/message/31705840/). Our paper explains in detail the cause of the energy drift and is available on arXiv (http://arxiv.org/pdf/1507.07081.pdf) prior to publication in J. Chem. Phys.

I implemented the eHEX algorithm in LAMMPS (9Dec14) as “fix ehex” and tested it for a LJ system and SPC/E water. In both cases, we found at least a hundredfold decrease in the energy drift. The HEX algorithm can still be used with the new fix by providing the keyword “hex”. In order to make the additional coordinate integration of the eHEX algorithm compatible with (small) rigid bodies, I implemented some additional methods in FixShake and FixRattle to allow the new fix to update constraining forces as a “friend class”.

If the developers consider the eHEX algorithm a useful feature for LAMMPS, I would be happy to integrate it in the latest version of LAMMPS and update the doc pages.

Best wishes,
Peter

Very cool that your original frustrations with the heat algorithm turned out into a new one with superior
performance. This certainly illustrates nicely the whole discovery cycle, i.e., a case is found where an old/established method fails, tests are done to ensure the problem is not on the method implementation, and finally, after (probably) some long hours of thinking you come up with a finer solution that contains the old method as a limiting case.
Congrats on the article,
Carlos

yes, sounds great. You can send the files to me or Axel.

thanks Peter,
Steve

Dear Steve and Carlos

Thank you very much for your comments!

I will send you the source files in a separate email, Steve.

Thanks,
Peter