Hi Aidan,
Thanks for picking up on the mixup.
Your diagnosis of the two thermostats messing each other up sounds promising, so I tried to implement the code you provided, but at the first line of “fix_modify 1 temp tempB” I get the error of, “ERROR: Illegal fix_modify command (…/fix_rigid_nh.cpp:1302).” I’m using the June 28, 2014 version of LAMMPS, and this is what I see in the source code:
int FixRigidNH::modify_param(int narg, char **arg)
I believe this error message was due to a bug that was recently fixed in fix_rigid_nh.cpp
Hi Aidan,
Thanks for picking up on the mixup.
Your diagnosis of the two thermostats messing each other up sounds
promising, so I tried to implement the code you provided, but at the first
line of "fix_modify 1 temp tempB" I get the error of, "ERROR: Illegal
fix_modify command (../fix_rigid_nh.cpp:1302)." I'm using the June 28, 2014
version of LAMMPS, and this is what I see in the source code:
int FixRigidNH::modify_param(int narg, char **arg){
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
<This is where I get the errorDo you not get this peculiar error? I'm not sure why 'fix_modify temp' would
require a fix that modifies the pressure...I see this same line of code in
the most current version of LAMMPS too.
the *current* version of LAMMPS does not have this line as it is a bug
to have this test there. it was removed on june 6th 2016. this was
part of a consolidation of multiple nose-hoover thermostat fixes which
had, under rare circumstances not properly initialized data and also
had inconsistent and slightly confusing naming conventions for
internal variables.
let me repeat, with the current version of LAMMPS, i cannot reproduce
the original issue you are reporting, even if only using a single
velocity command. all i see is a little creep of the framework.
i would definitely recommend to update to the latest patchlevel of
LAMMPS and re-check your simulations, there have been several bugfixes
and improvements since for the rigid body fixes, especially the ones
with nose-hoover thermostat. it is quite possible, that the effect you
are seeing, is due to one of those bugs.
axel.
I reran the simulation using the current development version of LAMMPS (June 28, 2016), which did indeed get rid of the error message. Unfortunately, it didn’t get rid of the system drift.
Axel, would you be willing to send me a log file in which you print out the COM of the total system? When I do this by adding the lines:
#define computes
compute comT all com
compute comB Benzene com
compute comM MOF com
compute tempB Benzene temp/com
compute tempM MOF temp/com
thermo_style custom step temp pe etotal c_comT[1] c_comT[2] c_comT[3] c_comB[1] c_comB[2] c_comB[3] c_comM[1] c_comM[2] c_comM[3]
to the simulation file I had sent out in my first email, I see that comT drifts. If I understand you correctly, you’re saying that in your simulation comT doesn’t drift, since the drift in comB and comM cancel each other out. If the comT not staying in place is due to the machine I’m running on, it’d be good to know. If so, and the comT stays in place on another machine, I’d move to another machine and be fine with tethering comM with a spring. Since I’m seeing comT move though, that tells me that something is funny with the physics. I’m attaching my runfile and log file (with ‘-axel’ suffixes) showing the drift to this email.
Aidan, I ran the input file you sent that uses ‘compute temp/com’ and zeroes out velocities after equilibration, but unfortunately I still get system drift during production. I’m attaching the runfile and log file (with ‘-aidan’ suffixes) to this email as well so you can see that comT continues to drift even after the equilibration period.
Does all this maybe signal that I need to give up on combining ‘fix NVT’ with ‘fix rigid NVT,’ and that I should instead put some stiff bonds/angles into my benzene model so I can use a single ‘fix NVT’ to time integrate everything?
Thanks again for the extensive help.
run-aidan (1.84 KB)
run-axel (1.11 KB)
log.lammps-aidan.gz (163 KB)
log.lammps-axel.gz (104 KB)
I still get system drift during production.
Drift happens. I would not necessarily expect such a complex simulation to be completely free of drift. You have two subsystems interacting with each other, each controlled by a separate thermostat, and one of the systems is imposing rigid constraints. It is not surprising that some drift occurs, although I think using temp/com and rezeroing the COM velocities after equilibration probably helped quite a bit.
Looking at log.lammps-aidan, I see that after rezeroing, comT starts out small, but then seems to slowing sample a somewhat random distribution centered on zero, with standard deviation about 5e-05 A/timestep (~10 m/s, much less than thermal speed of atoms). I think I have seen this behavior before. As long as the drift velocity remains small and bounded, there should be no consequences for the integrity of the simulation, but it may caused you difficulty analyzing the output. So:
If the drift is not causing you any trouble, I would leave it alone.
If it is causing you trouble, there are ways to suppress it, for example by gently tethering some or all of your framework atoms. Regardless of what method you use, you should verify that the underlying structural and dynamical properties are not altered
Drift happens.
Aidan
Sorry to bother again, but…
If comT were performing a random walk centered on its initial coordinate, I wouldn’t be so worried. I might just chalk it up to integration precision errors or some other noise source and be fine. However, I don’t see this when I look at log.lammps-aiden. See the attached graph that shows comT (in the x-direction) during the production timesteps. It grows massively in a single direction.
What is especially odd to me is that when I print out the forces (in the x-direction) via:
I inferred my random walk observation by taking the derivative of the graph that you show above. I observed a distribution of velocities centered on zero, with standard deviation 5e-05 A/timestep, and a long correlation time. Your graph is consistent with random walk behavior, albeit with a long correlation time. You say it grows massively in a single direction. I don’t know what that means. A random walk does indeed grow without bound, but not in a single direction (for t >> correlation time).
The overall behavior is probably due to a weak interaction between the two independent thermostats. Your statements about total force appear correct, but are only relevant to pure NVE dynamics. The thermostats do not modify the forces directly, and so it is not surprising that the average properties of the forces behave the same as NVE.
Bottom lines:
Aidan
Just adding this to the thread for the future reference of others having a similar problem. No need to respond.
One reason I had suspected that an error was present is because in the two other MD codes I’m working with (Tinker and RASPA2), no COM drift occurs for a system very similar to the one I was working with here. This led me to suspect that something wrong was going on in the one simulation package in which COM drift was occurring (LAMMPS).
I looked into the source code of Tinker and RASPA2, and I see that actually those other packages would exhibit COM drift, but their source codes automatically hide it in ways that LAMMPS allows the user to decide for himself. In Tinker, the linear momentum of the COM is zeroed every timestep (similar to LAMMPS’s ‘fix 1 all momentum 1 linear 1 1 1’) and in RASPA2 the coordinates are moved so that the COM of the framework maintains its position (similar to LAMMPS’ 'fix 1 MOF recenter INIT INIT INIT shift all).
In fact, and as Aidan pointed out, the linear momentum of the COM § should not be conserved in the canonical ensemble; instead Pe^{\eta} is conserved (see doi.org/10.1063/1.1378321). If the linear momentum is set to 0 initially, it should remain at 0 (because 0e^{\eta}=P*e^{\eta} only when P is 0), but if minor numerical errors bring it away from 0, it will oscillate along with the thermostat’s \eta parameter. I looked at the data from some longer runs, and Aidan is correct that it appears to be a random walk.
I believe that when I ran the simulations with just the MOF or just the benzene, these small numerical errors that brought the linear momentum away from 0 didn’t occur because the ‘fix’ was able to control for them, but when two ‘fixes’ were in place the small errors crept in; in fact, when I ran a simulation with just the MOF present (so just 1 ‘fix’) but with an initial non-zero linear momentum, I saw the linear momentum oscillating rather than being kept constant. This is further evidence that the linear momentum is not conserved with the Nose-Hoover thermostat but should in fact oscillate, and that there’s no problem with LAMMPS.
Given this, I think that it’s completely proper way to put ‘fix 1 all momentum 1 linear 1 1 1’ at the beginning of the run, since it will just be making adjustments to numerical error, and the linear momentum should stay at 0 if it starts out at 0. This makes inconsequential the decision of whether to modify the calculated temperature to remove the bias of the COM motion, since it will be negligibly small.
Just adding this to the thread for the future reference of others having a
similar problem. No need to respond.One reason I had suspected that an error was present is because in the two
other MD codes I'm working with (Tinker and RASPA2), no COM drift occurs for
a system very similar to the one I was working with here. This led me to
suspect that something wrong was going on in the one simulation package in
which COM drift was occurring (LAMMPS).I looked into the source code of Tinker and RASPA2, and I see that actually
those other packages would exhibit COM drift, but their source codes
automatically hide it in ways that LAMMPS allows the user to decide for
himself. In Tinker, the linear momentum of the COM is zeroed every timestep
(similar to LAMMPS's 'fix 1 all momentum 1 linear 1 1 1') and in RASPA2 the
coordinates are moved so that the COM of the framework maintains its
position (similar to LAMMPS' 'fix 1 MOF recenter INIT INIT INIT shift all).In fact, and as Aidan pointed out, the linear momentum of the COM (P) should
not be conserved in the canonical ensemble; instead P*e^{\eta} is conserved
(see doi.org/10.1063/1.1378321). If the linear momentum is set to 0
initially, it should remain at 0 (because 0*e^{\eta}=P*e^{\eta} only when P
is 0), but if minor numerical errors bring it away from 0, it will oscillate
along with the thermostat's \eta parameter. I looked at the data from some
longer runs, and Aidan is correct that it appears to be a random walk.I believe that when I ran the simulations with just the MOF or just the
benzene, these small numerical errors that brought the linear momentum away
from 0 didn't occur because the 'fix' was able to control for them, but when
two 'fixes' were in place the small errors crept in; in fact, when I ran a
simulation with just the MOF present (so just 1 'fix') but with an initial
non-zero linear momentum, I saw the linear momentum oscillating rather than
being kept constant. This is further evidence that the linear momentum is
not conserved with the Nose-Hoover thermostat but should in fact oscillate,
and that there's no problem with LAMMPS.Given this, I think that it's completely proper way to put 'fix 1 all
momentum 1 linear 1 1 1' at the beginning of the run, since it will just be
making adjustments to numerical error, and the linear momentum should stay
at 0 if it starts out at 0. This makes inconsequential the decision of
whether to modify the calculated temperature to remove the bias of the COM
motion, since it will be negligibly small.
you are overlooking one detail here: the command 'fix 1 all momentum 1
linear 1 1 1' will remove energy systematically.
since a finite size system will have finite size oscillation of the
COM momentum there is an optimum period for removing the momentum with
minimal impact. but in any case, you'll have a systematic energy
removal unless you also use the "rescale". many years back, i've been
involved in a project about pseudorotations in CH5+ with
car-parrinello MD and adding removed energy back had a noticeable
effect. it may be less in your case, however, if you prefer to use fix
momentum, i would advise against calling it at every step and for
using the rescale flag.
still, since part of your system is a stationary framework, i am still
very much in favor of simply adding a weak tether to that framework
group with fix spring for two reasons: a) you are less interested in
the framework dynamics and since it would macroscopically be part of a
large stationary object, that tether is actually a quite physically
meaningful representation of the environment for your framework. and
b) the diffusing liquid should have its momentum removed automatically
through the tether of the framework, without having to modify its
dynamics directly. as i assume, that this is the dynamics you care
about, you should have the minimum taint this way, in my personal
opinion.
axel.
I agree with Axel’s suggestion to use a weak tether instead of fix momentum, for the reasons he stated.
Also, I think it is important to note that the drift you are observing is probably not due to some kind of inaccuracy or error in the LAMMPS time integration scheme, but rather it is due to the composite structure of your simulation. Nose-Hoover-style thermostats, such as the one implemented as fix nvt in LAMMPS do indeed have the nice mathematical property of conserving zero total momentum in a system with zero initial total momentum, as elegantly described in the paper you cited, with which I am very familiar. However, that is only true for a dynamical systems composed of an isolated Hamiltonian sub-system and a suitable set of extended non-Hamiltonian variables interacting with the Hamiltonian subsystem. In your case, you have built a composite dynamical system consisting of two Hamiltonian subsystems, each with its own set of extended non-Hamiltonian variables. The two Hamiltonian subsystems are interacting with each other. In addition, one of the Hamiltonian subsystems is imposing rigidity constraints, the properties of which I don’t know much about. It is not surprising that such a composite system does not conserve zero total momentum. The best that you can hope for is that the total momentum remains bounded and small (no flying ice cubes), which is what you are getting.
Aidan