fix heat and momentum conservation problem


I am experiencing an issue with momentum conservation when employing the “fix heat” routine to create a thermal flux across the interface of two nano-particles. Thermal energy is added to all atoms from one nano-particle and removed from all atoms in another particle:

Remove energy from particle 1

fix Qout 1 heat 1 -${flux}

Apply heat to particle 2

fix Qin 2 heat 1 ${flux}

The particles develop an increasing center-of-mass velocity in the periodic computational cell. Total energy is roughly conserved, but kinetic energy increases with time. I have ensured that there is no net initial momentum of the system, but the nature of the physical configuration necessitates some harmonic oscillations between the two particles, so each particle may have some nominal momentum individually. However, I would not expect a velocity scaling routine to systematically increase any initial momentum value, and the routine description mentions that aggregate momentum of the system should be conserved. This effect is only seen with the introduction of the “fix heat” routine, however. The particles remain in the center of the simulation cell when the “fix heat” routine is removed and all other conditions are kept constant.

Momentum fixes, like

Fix momentum of each particle individually

fix mv1 1 momentum {mvfix} linear 1 1 1 angular fix mv2 2 momentum {mvfix} linear 1 1 1 angular

do not work, since they drain energy from the system by effectively damping the harmonic oscillations of the two-particle system.

My current solution idea involves calculating the linear and angular velocities of each nano-particle, subtracting these values from each atom’s velocities, performing the energy flux routine (either by “fix heat” or an alternative) and adding back the linear and angular velocities calculated from the initial step to avoid energy damping. My problem here is that “fix heat” is a recurring or looped command and the “velocity” command is only performed a single time. In addition, the problem will become more difficult if I still find the “fix heat” routine is imparting some momentum to each particle even after explicitly zero-ing it before each scaling step.

Question #1: Can you explain/show the mathematics of the “fix heat” routine? Are all velocities scaled linearly? Do you have any idea why a net momentum is imparted or aggregate momentum is not conserved?
Question #2: Regarding my plan to explicitly subtract and then re-add momentum into each particle, is there a way to perform “velocity set” commands for every time step, before and after each instance of “fix heat” implementation?

Thanks for your time and efforts.
-Rich Salaway

The particles develop an increasing center-of-mass velocity in the periodic computational cell. Total energy is roughly conserved, but kinetic >energy increases with time.

This is consistent with the so-called "flying ice-cube" effect that
is sometimes seen in molecular systems where there are different
sets of degrees of freedom, e.g. bonds and translational.

Bascially anytime you are themostatting or otherwise adding/subtracting
energy to a system, you can aggravate this kind of problem.

The method used in fix_heat.cpp is simple; take a look at end_of_step().
It just scales the velocities relative to the V of the com.

Paul can comment on whether he's seen this kind of momentum
issue before when using 2 fix heat commands.


I haven't seen this kind of momentum issue before when using fix heat, but I also have not had a great deal of experience in using fix heat.

For the math behind fix heat, see section 6.6 of the following paper:

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

I'll send you the 11 MB pdf of the paper off list. Might be accessible to you from:

It is possible that we've implemented the algorithm incorrectly in LAMMPS, but I doubt it. Please let us know if you find otherwise. It should conserve momentum, so I'm not sure what the issue is with your simulation. May be something to do with your two heat groups not being truly independent of each other.

To answer your second question, yes I think it is possible to reset velocities on every timestep, but that seems drastic. I'd recommend resolving the problem in a more elegant manner --- figure out why you're getting the momentum gain and resolve the root cause.