[EXTERNAL] Re: fix heat and momentum conservation problem

[Richard_Salaway]

Hi Steve and Paul,

Thanks so much for your prompt and thorough responses. Thank you for the Sandia Report as well. After looking it over and fix_heat.cpp, it does seem to me that any linear momentum in a group of atoms undergoing fix_heat should experience linear momentum conservation.

From fix_heat.cpp

v[i][0] = scale*v[i][0] - vsub[0]

is equivalent to

Vi’ = scale*(Vi - Vcm) + Vcm

which is essentially what I hoped to do with multiple LAMMPS command lines. My system consists of two carbon nanotubes forming a cross junction. Two separate fix_heat routines are applied to all atoms in each tube, individually. There typically exists minor periodic oscillations of the two-tube structure, with the separation distance between the two tubes oscillating much like a diatomic pair.

Unfortunately there is a net linear and angular momentum imparted on each tube individually, with the introduction of the fix heat routine. As the fix heat routine is applied to all atoms in each tube, I would expect any small linear momentum values to be conserved. However, the momentum of each tube (and therefore the kinetic energy) is affected by the fix_heat routine in this configuration.

While the scaling routine above should conserve linear momentum, I suspect angular momentum about a center of mass would be improperly scaled. My current idea is that the fix_heat routine imparts an improper angular velocity to the tube each time, and this may cause the tubes to effectively “roll along” each other (thus, the introduction of linear momentum).

I have had difficulty resolving this issue for several reasons:

1. The angular momentum adjustment must come at every time step to correspond to each incidence of fix_heat.
2. My atom_style (atomic) prohibits the use of straight forward commands like “set angmom”.
3. Any zeroing of angular momentum would artificial drain energy from the system. I believe I need an iterative procedure ie; a) removal of angular velocity, b) adjustment of energy through a scaling of velocities (fix heat) c) re-introduction of previous angular velocity.

I have tried to attached a movie file that demonstrates the non-zero linear momentum of individual tubes seen in the simulation. Though only a few sparse snapshots, the red tube is seen to transverse through periodic boundary conditions (right to left) and the blue tube goes through periodic boundaries into and out of the page. I’ve attached the LAMMPS script that produced this result. The tubes in the restart file did not have individual momentum, and I also zero total momentum by the “velocity all zero linear” command. Though I apply this command to all atoms at once in this file, I have seen the same behavior for all different types of momentum zeroing routines. For example, applying “fix momentum” to each tube individual initially suppresses the introduction of additional momentum, however the tubes develop an increasing momentum as soon as the “fix momentum” condition is relaxed. Unfortunately, “fix momentum” is not a valid solution b/c it drains kinetic energy from the system by quenching the natural oscillations between tubes.

I was hoping you could offer some advice or guidance:
Am I overlooking something?
Can you think of a way to resolve this issue using LAMMPS subroutines?
If your suggestion is altering the fix_heat routine, can you give some details about the relevant variables I need to adjust and their meaning?

Thanks again for your time and efforts.

-Rich Salaway

FLUXregppp.junction (1.95 KB)

[sjplimp]

I'll let Paul comment on whether it makes sense to try to conserve angular
momentum in what fix heat does. Angular momentum isn't conserved
generally in a periodic system. I don't understand how random ang momentum
changes (due to random fix heat perturbations) cause the 2 tubes
to roll along each other. The MPG movie was too jumpy for me
to see anything. Are you saying that applying fix heat to a single
nanotube causes it to spin in one direction around its long axis? I have a hard
time believing that. As opposed to just undergo small random oscillations
in an angular sense.

Steve

[Crozier_Paul_S]

I don't think that fix heat tries to conserve angular momentum. Isn't clear to me whether or not it might introduce some non-conservation of angular momentum, but it should conserve linear momentum.

My guess is that there's something else fundamentally problematic in the way you've got your simulation set up. Even though it appears that the blame lies on fix heat, it may be that the problem is in the way you're using the fix, or how it is interacting with other parts of the simulation you're trying to run.

I agree with Steve that it seems odd that applying fix heat to a single nanotube would cause it to spin. There's probably something else wrong in the way you've got your system set up.

Paul

[sjplimp]

I forgot, there is a fix momentum command that has
an angular option. You can use it every N steps to
zero the angular momentum of a group around its COM.
If you did this twice, once for each nanotube, you should
see no rotation. Unless, as Paul said, there is something
else in your model that is driving the rotation.

It is true that a nanotube has a pretty small moment
of inertia around its long axis, so it is maybe possible
that it will start spinning with small impulses. But I still
think the impulses from fix heat will be randomly in both
directions, so wouldn't think you'd get much net spin.

Steve

[_Russell_John_Thomas]

Dear Sir,
You might consider imposing harmonic springs to prevent rotation of
a rod or nanotube.

For example, consider a nanotube aligned with the x-axis with center of
mass at the origin. You can divide the nanotube into two goups where
atom type 1 has z > 0 and atom type 2 has z < 0. Then you can assign
harmonic springs for each group in the y-axis to prevent rotation.

Essentially you are tethering the center of mass of each half-cylinder
rather than tethering the center of mass of the full nanotube which can
rotate on its axis.

For example:
velocity all create 1.0 4928459 dist gaussian
fix 1 all nve
group up type 2
group down type 1
fix 2 up spring tether 50.0 NULL 0.0 NULL 0.0
fix 3 down spring tether 50.0 NULL 0.0 NULL 0.0
fix_modify 2 energy yes
fix_modify 3 energy yes
fix 4 all langevin 0.0 2100.0 1.0 12345 zero yes

This solved my problem with unwanted rotation.

Sincerely,
John