Dear all,
I’m having a problem with keeping the total linear and angular momentum of a set of atoms at 0. I have stripped down my simulation to the simplest components, and the problem persists.
I find that, either in NVT or NVE, running a minimisation between runs causes the total angular momentum to change. This might be proper behavior; I can’t see much information in the documentation on how minimization is intended to handle velocities.
A similar problem occurs when using some fixes, such as fix viscous. I use the following steps:
fix NVTfix all nvt temp 2500.0 2500.0 0.1
minimize 1.0e-5 0 1000 10000
velocity all create 2500.0 23543 mom yes rot yes
fix equilibrate_viscous all viscous 0.1
run 40000
unfix equilibrate_viscous
run 100000
At the start of the second run command, immediately after the unfix command, the system appears to get a significant spinning ‘kick’, violating both linear and angular momentum. It looks like a football being kicked.
I have attached my complete input scripts for the last case. I would like to know whether this is expected behavior. If it is, is there an efficient way to relax a system which does not violate momentum conservation?
Regards,
Gil Rutter
LJ_cluster.lmp (855 Bytes)
LJ_cluster_test.in (1.16 KB)
Dear all,
I'm having a problem with keeping the total linear and angular momentum of a
set of atoms at 0. I have stripped down my simulation to the simplest
components, and the problem persists.
I find that, either in NVT or NVE, running a minimisation between runs
causes the total angular momentum to change. This might be proper behavior;
I can't see much information in the documentation on how minimization is
intended to handle velocities.
that is the point, the minimization does not use
or change velocities. so a change to the total
angular moment should be expected when the
minimization changes the coordinates.
A similar problem occurs when using some fixes, such as fix viscous. I use
the following steps:
fix NVTfix all nvt temp 2500.0 2500.0 0.1
minimize 1.0e-5 0 1000 10000
velocity all create 2500.0 23543 mom yes rot yes
fix equilibrate_viscous all viscous 0.1
run 40000
unfix equilibrate_viscous
you'd probably be better off to use fix nve in
combination with fix langevin here for the
same effect. i am not 100% certain, but i
could imagine some unwanted interaction
between fix nvt and fix viscous.
run 100000
At the start of the second run command, immediately after the unfix command,
the system appears to get a significant spinning 'kick', violating both
linear and angular momentum. It looks like a football being kicked.
in general, is is a good idea to remove total translational
and rotational momentum when you make drastic changes
changes to a system, that is sensitive to accumulating
those, like yours seems to be.
I have attached my complete input scripts for the last case. I would like to
know whether this is expected behavior. If it is, is there an efficient way
to relax a system which does not violate momentum conservation?
see above my suggestion about using fix langevin.
but still i would recommend to remove total momenta
when switching from equilibration to production,
and perhaps even occasionally remove them during
the run, if your system is sensitive to accumulating
those for technical reasons.
cheers,
axel.
If you want to be really servere about the momentum conservation you should use a different compute for the temperature of the thermostats, the default compute acts in all degrees of freedom, you should always remove the center of mass degree of freedom to calculate the temperature of the system, do something like:
fix f1 all nvt temp 1000.0 1000.0 0.1
compute Tcm all temp/com
fix_modify f1 temp Tcm.
I think that the ‘kick’ is to be expected: you are suddenly removing all viscosity of your system.
A personal question: is it really necessary to use minimize + fix viscous + create velocity to equilibrate the system? I’ve never used all this things together, usually (in my simulations…) just the create velocity + thermostat are enough to a (really) rapid equilibration.
Regards,
Rodrigo Freitas
Thanks for your reply, Axel.
The fact that making alterations to the simulation causes major violations of angular and linear momentum conservation is quite worrying. Could you please explain where the kick which I witness at that unfix command originates? Is LAMMPS programmed in such a way that some commands are intended to be suffixed with a command to restore these conservations? My problem is that it’s not clear that simply zeroing the angular and linear momentum would preserve the ensemble from which I was originally sampling.
Regards,
Gil
Thanks for your reply, Axel.
The fact that making alterations to the simulation causes major violations
of angular and linear momentum conservation is quite worrying. Could you
please explain where the kick which I witness at that unfix command
originates? Is LAMMPS programmed in such a way that some commands are
intended to be suffixed with a command to restore these conservations? My
problem is that it's not clear that simply zeroing the angular and linear
momentum would preserve the ensemble from which I was originally sampling.
but in both cases (minimization or fix viscous)
you already "change" the ensemble. there is
nothing to preserve.
axel
Rodrigo,
If you want to be really servere about the momentum conservation you should use a different compute for the temperature of the thermostats, the default compute acts in all degrees of freedom, you should always remove the center of mass degree of freedom to calculate the temperature of the system, do something like:
fix f1 all nvt temp 1000.0 1000.0 0.1
compute Tcm all temp/com
fix_modify f1 temp Tcm.
Thanks for this. Unfortunately, I would really like my system to conserve its momentum properly so that this isn’t needed!
I think that the ‘kick’ is to be expected: you are suddenly removing all viscosity of your system.
The particles’ velocities go to 0 in the viscosity regime after relaxation, so I wouldn’t expect any significant forces or velocities at all when I remove the fix - especially not ones that violate momentum conservation.
A personal question: is it really necessary to use minimize + fix viscous + create velocity to equilibrate the system? I’ve never used all this things together, usually (in my simulations…) just the create velocity + thermostat are enough to a (really) rapid equilibration.
The reason I have a hodgepodge of relaxation techniques is that I am experimenting with this momentum conservation problem - I don’t actually need them all for my simulation.
Gil
The reason I have a hodgepodge of relaxation techniques is that I am
experimenting with this momentum conservation problem - I don't actually
need them all for my simulation.
equilibrium is a state function, i.e. it doesn't matter
how you get there. i simply think it is not reasonable
to expect perfect momentum conservation while you
manipulate your system. in re fix viscous, keep in mind
that fix nvt manipulates accelerations, not velocities!
thus i think it is completely valid approach to use
whatever tricks you can find to coax your system
(close!) to equilibrium and then just clear the
angular momentum with
velocity zero angular
and then equilibrate a bit more with your final
simulation setup and *then* start production.
that is how it should be anyway.
axel.
Axel,
I’ve done a bit more research on what exactly is the problem.
but in both cases (minimization or fix viscous)
you already “change” the ensemble. there is
nothing to preserve.
You are correct. However, my original simulation is intended to run in a multi-canonical ensemble, using fix external to modify forces in a way that does not intrinsically violate momentum conservation. I find that momentum conservation is not obeyed here too, and in this case it is a real problem.
I have altered the code of fix_external.cpp as follows:
void FixExternal::post_force(int vflag)
{
// invoke the callback in driver program
// it will multiply the forces by STMD_factor
double STMD_factor = 1.0;
(this->callback)(ptr_caller,update->ntimestep,atom->nlocal,atom->tag,atom->x,&STMD_factor);
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (STMD_factor > 1.0){
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
f[i][0] *= STMD_factor;
f[i][1] *= STMD_factor;
f[i][2] *= STMD_factor;
}
}
}
}
I’ve also adjusted various other bits of the source code, so that the callback function wants the correct parameters. The behavior I observe with this is that the system behaves normally for large time periods of the simulation, then suddenly gets linear and angular momentum and starts moving as a whole.
This may be an unwanted interaction with the NVT thermostat, in any case, is the problem obvious?
Gil
Axel,
I've done a bit more research on what exactly is the problem.
but in both cases (minimization or fix viscous)
you already "change" the ensemble. there is
nothing to preserve.
You are correct. However, my original simulation is intended to run in a
multi-canonical ensemble, using fix external to modify forces in a way that
does not intrinsically violate momentum conservation. I find that momentum
conservation is not obeyed here too, and in this case it is a real problem.
I have altered the code of fix_external.cpp as follows:
void FixExternal::post_force(int vflag)
{
// invoke the callback in driver program
// it will multiply the forces by STMD_factor
double STMD_factor = 1.0;
(this->callback)(ptr_caller,update->ntimestep,atom->nlocal,atom->tag,atom->x,&STMD_factor);
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (STMD_factor > 1.0){
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit) {
f[i][0] *= STMD_factor;
f[i][1] *= STMD_factor;
f[i][2] *= STMD_factor;
}
}
}
}
I've also adjusted various other bits of the source code, so that the
callback function wants the correct parameters. The behavior I observe with
this is that the system behaves normally for large time periods of the
simulation, then suddenly gets linear and angular momentum and starts moving
as a whole.
This may be an unwanted interaction with the NVT thermostat, in any case, is
the problem obvious?
hmm... you are scaling all forces by the same factor.
do all your atoms have the exact same mass?
axel.
In this case, they do. But that is not necessary for STMD to work, including conservation of linear and rotational momentum.
Gil
In this case, they do. But that is not necessary for STMD to work, including
conservation of linear and rotational momentum.
well i would suggest to do an experiment:
comment out the callback and do a test run
with scaling factors 1.0 2.0 5.0 10.0 0.1 0.01
and see how this turns out. and of course
compare the 1.0 case to not using fix external
at all.
thanks,
axel.
Gil,
If I understand correctly, your original complaint is that angular
momentum is not conserved when you run fix nvt with your fix external.
I don't find this behavior surprising. In general, any deviation
from Hamiltonian mechanics (fix nve) will violate all conservation
laws. Nose-hoover dynamics (fix nvt) is a special case. It is
non-Hamiltonian but still conserves linear and angular momentum (I
think). Also, angular momentum is never truly conserved in a periodic
system. Each time an atom i crosses a periodic boundary, there is a
kick of order L x p_i. However, once you modify the Nose-Hoover
dynamics, all bets are off, unless you can prove that the new
equations of motion do indeed conserve angular momentum (by taking the
time derivative and showing that it is identical to zero).
Aidan
Aidan,
Thanks for your response. The system receives a ‘kick’, just like with fix viscous. These impulses occur at random-seeming times, and either give the system linear and angular momenta, or zero the linear and angular momenta. Given that using fixed boundaries doesn’t help, one would expect more continuous changes with your explanation, I think.
Gil
Hi
I’ve found that scaling factors below 1.0 tend to cause problems, and other scaling factors don’t. In my code, as you can see, values of less than 1.0 aren’t allowed and there are still problems, which suggests to me that the bug might arise when the value is lowered. I’ll let you know if I find anything else interesting.
Gil
Okay, I was barking up the wrong tree - I think fix external is fine. I have now made the problem occur using an extremely simple set-up of 5 Lennard-Jones atoms, with no fixes except NVT. Also replacing fix NVT with fix temp/rescale and NVE produces similar problems.
My scripts are attached again. In this case, the problem occurs after about 8.5 *10^6 time-steps and is very obvious if you visualize the simulation. It should be very easy for anyone to run these scripts and see what they get. Could someone please run them so I know if this issue is on my side only, or not?
Gil.
LJ_cluster_test.in (1.1 KB)
LJ_cluster.lmp (345 Bytes)
Dear all,
Sorry if any messages have been ignored, it seems that were was a problem with my inbox.
Gil