Hello,
I’ve been running some simple simulations of two harmonically bonded molecules to see the effect of different thermostats on COM momentum conservation, the flying ice cube effect, and a few other things (as inspired by what I learned from the previous thread, “Atom drift (perhaps due to the flying ice cube effect”).
One thing I’m doing is putting two separate Nose-Hoover thermostats on each molecules. When I do this and print out the temperature, the temperature of the system quickly drops to 2/3 of the value I’m thermostatting to and stays there, e.g. if I put both thermostats to 300 K, the temperature goes to 200 K.
I’m assuming that I’m doing something very stupid that I can’t find. Any help? The run file is below, with system.data attached. I’m using the May 14, 2016 version of LAMMPS.
Thanks!
system.data (370 Bytes)
Hello,
I've been running some simple simulations of two harmonically bonded
molecules to see the effect of different thermostats on COM momentum
conservation, the flying ice cube effect, and a few other things (as
inspired by what I learned from the previous thread, "Atom drift (perhaps
due to the flying ice cube effect").
One thing I'm doing is putting two separate Nose-Hoover thermostats on each
molecules. When I do this and print out the temperature, the temperature of
the system quickly drops to 2/3 of the value I'm thermostatting to and stays
there, e.g. if I put both thermostats to 300 K, the temperature goes to 200
K.
I'm assuming that I'm doing something very stupid that I can't find. Any
help? The run file is below, with system.data attached. I'm using the May
14, 2016 version of LAMMPS.
this is a *tiny* system and "temperature" as such is not well defined
under these circumstances.
the discrepancy you see is due to the fact that temperature is
computed based on the number of degrees of freedom. for a bulk system
with periodic boundaries with periodic boundaries the temperature
compute will remove 3 DOFs due to translational invariance. normally,
that difference is quite small, but for your setup this becomes
massive, since with two individual thermostats, you remove a total of
6 DOFs (out of 12), while the global temperature compute will only
remove 3 DOFs (out of 12). to make those consistent, you have to make
certain, that the temperature compute is consistent. in your case, it
is best to turn off the translational invariance with:
compute_modify thermo_temp extra 0
compute_modify 1_temp extra 0
compute_modify 2_temp extra 0
a "cleaner" solution in this case would be to allow fractional DOFs
and distribute the globally removed 3 DOFs across the two subsystems
(1 1/2 each).
axel.
BTW: allowing fractional DOFs through compute_modify is extremely
simple to add to LAMMPS.
only two small modifications are needed, and i think this is a change
that should be included into the official LAMMPS version.
with the change below, you would only need to change your input like this:
compute_modify 1_temp extra 1.5
compute_modify 2_temp extra 1.5
and you should also get 300K as global "temperature".
axel.
diff --git a/src/compute.cpp b/src/compute.cpp
index 63a00da..662f60f 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -133,7 +133,7 @@ void Compute::modify_params(int narg, char **arg)
while (iarg < narg) {
if (strcmp(arg[iarg],"extra") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command");
- extra_dof = force->inumeric(FLERR,arg[iarg+1]);
+ extra_dof = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dynamic") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command");
diff --git a/src/compute.h b/src/compute.h
index 0c04286..ce47b49 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -149,7 +149,7 @@ class Compute : protected Pointers {
int instance_me; // which Compute class instantiation I am
double natoms_temp; // # of atoms used for temperature calculation
- int extra_dof; // extra DOF for temperature computes
+ double extra_dof; // extra DOF for temperature computes
int fix_dof; // DOF due to fixes
int dynamic; // recount atoms for temperature computes
int dynamic_user; // user request for temp compute to be dynamic
Efrem,
A different approach from Axel’s would be to thermostat each molecule relative to its center of mass.
compute 1 mol1 temp/com
compute 2 mol2 temp/com
fix_modify 1 temp 1
fix_modify 2 temp 2
This way, the molecule-molecule interaction will not be explicitly thermostatted, but upon equilibration the system will reach the target temperature on average with large fluctuations.
-David