Using “fix momentum” and “fix nve”: How the total energy is conserved?

Jafary_Zadeh_Mehdi · February 25, 2013, 11:00am

Dear LAMMPS developers/users,
I am using LAMMPS to study surface diffusion of a cluster on a substrate. My system is consisted of two groups of atoms: 1000 atoms in the “substrate" group and 100 atoms is the “cluster” group.
I used the “fix nve” command to run my simulations. To set the temperature of the system I used separate “velocity” commands for each group of atoms with no initial rotation or translation motion. I also used “fix momentum” with “angular” keyword to freeze the rotational motion of the “cluster” group during the simulation.
According to the LAMMPS documentation, when we use "fix momentum " with “angular” keyword on a group of atoms (which is the “cluster” group in my case), the angular momentum is zeroed by subtracting a rotational component from each atom of that group . So, my question is:
Is there any possibility of energy leakage in the system when I use “fix nve” and “fix momentum with “angular” keyword on the cluster group of my system?
Here are my commands in my LAMMPS script:
group substrate id <= 1000

1000 atoms in group substrate

group cluster id > 1000

100 atoms in group cluster

fix cluster_fixMom cluster momentum 1 angular
velocity substrate create {nveEnergyDoublingTemp} 6456222 rot yes dist gaussian velocity cluster create {nveEnergyDoublingTemp} 6456222 rot yes dist gaussian
fix 1 all nve
run 10000000
Regards,
Mehdi

akohlmey · February 25, 2013, 11:04am

Dear LAMMPS developers/users,
I am using LAMMPS to study surface diffusion of a cluster on a substrate. My
system is consisted of two groups of atoms: 1000 atoms in the “substrate"
group and 100 atoms is the “cluster” group.
I used the "fix nve" command to run my simulations. To set the temperature
of the system I used separate "velocity" commands for each group of atoms
with no initial rotation or translation motion. I also used “fix
momentum” with “angular” keyword to freeze the rotational motion of the
“cluster” group during the simulation.
According to the LAMMPS documentation, when we use "fix momentum " with
“angular” keyword on a group of atoms (which is the “cluster” group in my
case), the angular momentum is zeroed by subtracting a rotational component
from each atom of that group . So, my question is:
Is there any possibility of energy leakage in the system when I use "fix
nve" and “fix momentum with "angular" keyword on the cluster group of my
system?

why ask this question when you can answer it yourself easily by just
running the simulation?

axel.

Jafary_Zadeh_Mehdi · February 25, 2013, 11:22am

Thanks for the fast reply Axel.
Actually, I have run this simulation in different temperatures and I have observed that in all cases the total energy of the system is conserved (which means that the NVE ensemble works!).
But, I cannot understand why? While the “fix momentum” with angular keyword subtracts a rotational component from a subgroup of atoms ("cluster"group in my system), it looks like some sort of kinetic energy leakage in the system. However, “fortunately” it does not happened in practice, and I wonder why?

Thanks a lot,
Mehdi

akohlmey · February 25, 2013, 1:21pm

Thanks for the fast reply Axel.
Actually, I have run this simulation in different temperatures and I have
observed that in all cases the total energy of the system is conserved
(which means that the NVE ensemble works!).
But, I cannot understand why? While the "fix momentum" with angular keyword

this is a different question. why don't you ask what you *really* want
to know right away?

subtracts a rotational component from a subgroup of atoms ("cluster"group in
my system), it looks like some sort of kinetic energy leakage in the system.

only if the "cluster" part of your system will actually pick up
rotational kinetic energy. have you measured how much of that it will
pick up, if at all?

again, this is something that can be easily assessed through simple
experimentation with your test system.

However, "fortunately" it does not happened in practice, and I wonder why?

there are always three explanations for things like this: 1) you do
have an error and it is being compensated by a second error that goes
into the other direction, 2) what you consider conserved energy is not
that well conserved after all, 3) your system doesn't behave the way
you expect or the contributions that you expect to be significant are
not.

these are all issues that can be easily tested by yourself. it is just
a matter to building inputs piece by piece and observing the impact of
each feature added. it is very good practice to do this.

axel.

Jafary_Zadeh_Mehdi · February 27, 2013, 2:42pm

Dear Axel,
Thanks a lot for your consideration and reply. As a brief reminder of my problem: I have a system consisted of two groups of atoms: an admolecule and a substrate. I run my simulation using fix nve while I ONLY froze the rotational motion of the admolecule using “fix momentum 1 angular” command. My time step was 1 fs. Although people might expect there must be some “energy leakage” in my simulations, but exact monitoring of energies indicates that my total energy (as well as PE and KE) is conserved during long time run (10,000,000 and beyond). So, my question was: why the energy is conserved in my simulation (which of course is a good event for me)?!

You replied me:
“there are always three explanations for things like this:”
So, let’s analyze your answer from bottom to top:

“3) your system doesn’t behave the way you expect or the contributions that you expect to be significant are not.”
No! I have precisely checked the movies, and I found that the system behaves as I expected: The molecule has a translational diffusive motion on the substrate without any rotation. (For your information, without applying fix momentum angular command, the natural diffusive motion of my admolecule also consisted of rotational motion around its center of mass).

what you consider conserved energy is not that well conserved after all,
I have precisely monitored the thermodynamics outputs in my log files. The total PE, KE and total energy of my system remains constant with an acceptable fluctuation around the equilibrium average. (Though I will back to this fluctuations in the next section).

1)you do have an error and it is being compensated by a second error that goes
into the other direction,
Maybe! So, please see whether my explanation is reasonable for you:
I have a system consisted of about 1000 atoms. Using fix momentum 1 angular, I slightly manipulated the velocities of about 50 atoms in my “admolecule” group at every time step to freeze its rotational motion. Since I eliminate 3 dimensional rotational motion of the admolecue, I “kill” 3 degrees of freedom in my system. Consequently, based on the equipartition theorem, the “turbulence” I impose to my system by adjusting the velocities of admolecule is in the order of “31/2k_BT”. On the other hand, in my simulations, I observed that the standard deviation of the total kinetic energy (and potential energy) of my system, which is due to the existence of energy fluctuations raised form Verlet integration technique, is more than one order of magnitude higher than 31/2k_BT. So, since the imposed turbulence in the kinetic energy (according to applying fix momentum angular) is much smaller than the numerical uncertainties of the energy calculations in fix nve, no accumulative leakage is observable during my simulations. (It is also noteworthy that in practice, the actual turbulence due to “fix momentum angular” can be much smaller than 31/2k_BT, because the rotational motion of the admolecule may not have enough time to be thermalized).

So, is my description above reasonable/convincing for you to describe “why there is no energy leakage” in my simulations?

Thanks a lot for your all your passionate considerations and help,
Mehdi

akohlmey · March 1, 2013, 8:11pm

Dear Axel,
Thanks a lot for your consideration and reply. As a brief reminder of my
problem: I have a system consisted of two groups of atoms: an admolecule and
a substrate. I run my simulation using fix nve while I ONLY froze the
rotational motion of the admolecule using “fix momentum 1 angular” command.
My time step was 1 fs. Although people might expect there must be some
“energy leakage” in my simulations, but exact monitoring of energies
indicates that my total energy (as well as PE and KE) is conserved during
long time run (10,000,000 and beyond). So, my question was: why the energy
is conserved in my simulation (which of course is a good event for me)?!

You replied me:
“there are always three explanations for things like this:”
So, let’s analyze your answer from bottom to top:

“3) your system doesn't behave the way you expect or the contributions that
you expect to be significant are not.”
No! I have precisely checked the movies, and I found that the system behaves
as I expected: The molecule has a translational diffusive motion on the
substrate without any rotation. (For your information, without applying fix
momentum angular command, the natural diffusive motion of my admolecule also
consisted of rotational motion around its center of mass).

2) what you consider conserved energy is not that well conserved after all,
I have precisely monitored the thermodynamics outputs in my log files. The
total PE, KE and total energy of my system remains constant with an
acceptable fluctuation around the equilibrium average. (Though I will back
to this fluctuations in the next section).

1)you do have an error and it is being compensated by a second error that
goes
into the other direction,
Maybe! So, please see whether my explanation is reasonable for you:
I have a system consisted of about 1000 atoms. Using fix momentum 1 angular,
I slightly manipulated the velocities of about 50 atoms in my “admolecule”
group at every time step to freeze its rotational motion. Since I eliminate
3 dimensional rotational motion of the admolecue, I “kill” 3 degrees of
freedom in my system. Consequently, based on the equipartition theorem, the
“turbulence” I impose to my system by adjusting the velocities of admolecule
is in the order of “3*1/2*kBT”. On the other hand, in my simulations, I
observed that the standard deviation of the total kinetic energy (and
potential energy) of my system, which is due to the existence of energy
fluctuations raised form Verlet integration technique, is more than one
order of magnitude higher than 3*1/2*kBT. So, since the imposed turbulence
in the kinetic energy (according to applying fix momentum angular) is much
smaller than the numerical uncertainties of the energy calculations in fix
nve, no accumulative leakage is observable during my simulations. (It is
also noteworthy that in practice, the actual turbulence due to “fix momentum
angular” can be much smaller than 3*1/2*kBT, because the rotational motion
of the admolecule may not have enough time to be thermalized).

So, is my description above reasonable/convincing for you to describe
“why there is no energy leakage” in my simulations?

you have missed the most important statement in my e-mail. let me
rephrase my statement as a question:

why do you expect your cluster to accumulate significant rotational
kinetic in the first place?

if your cluster doesn't have significant rotational kinetic energy in
the first place then fix momentum does not need to remove it. in other
words. re-run your calculation *without* fix momentum and observe (and
compute!) the rotational kinetic energy of that cluster.

axel.

Jafary_Zadeh_Mehdi · March 6, 2013, 9:12am

Hi Axel,

Thanks a lot for you’re following up and sorry for my a bit late reply.
“ re-run your calculation without fix momentum and observe (and compute!) the rotational kinetic energy of that cluster”.
I have performed this simulation, and without fix momentum 1 angular command, my cluster rotates. It takes some time that its rotational degrees of freedom (DOF) become thermalized and acquire average kinetic energy in the order of k_BT.

More specifically, my cluster is a C60 molecule on a graphene substrate. Our group has published a paper about surface diffusion of C60 on graphene in which we used just fix nve command to perform the simulations in the microcanonical ensemble to get rid of any possible artificial effects due to thermostats. At the beginning of simulations, we distributed the velocities of our cluster (C60 molecule) and substrate (graphene sheet) using separate velocity commands at desired temperature, so that there is no initial rotational/translational momentum in the C60 molecule. However, after a certain amount of time, the rotational/translational degrees of freedom of C60 thermalized and assured a thermal energy in the order of k_BT . (As u can see in the Fig. 3 of our attached paper, we have reported that there is an interplay (in the form of anti-correlation) between rotational and translational kinetic energies of the C60 molecule.

The purpose of our new set of simulations is to investigate the effect of ONLY C60 rotational Degreed of freedom (DOF) on its diffusive motion. So, in contrast to our previous simulations, we JUST added a fix fix_id C60group momentum 1 angular command to only freeze the rotational motion of cluster. As we have discussed before, the result is: the C60 molecule does not rotate, and meanwhile the energy of the system is conserved (with acceptable fluctuations).

“if your cluster doesn’t have significant rotational kinetic energy in the first place then fix momentum does not need to remove it.”
And your suggestion might be also right, as I explained in previous email,
Since the C60 does not have any initial rotational motion, and because at each time step, we freeze its rotation using fix momentum, its rotational DOF does not have enough time to acquire energy and be thermalized. So, it might be possible that the amount of rotational energy which must be removed at each step by fix momentum command is negligible in this case. However, the fix momentum command at the current LAMMPS code does not report how much kinetic energy it subtract from the system, so we cannot have an answer for this question. A possible practice is that someone rewrites the fix momentum module for this purpose.

Best regards,
Mehdi

Kinetic_nanofriction_a_mechanism_transition_2012.pdf (1.13 MB)