Berendsen thermostat

LAMMPS LAMMPS Beginners
1

Hello,

I am new to lammps, I have been studying the stability of adamantane following simulated annealing. I performed several tests to determine the optimal time step and damping parameter (tdamp) using the Berendsen thermostat. In my procedure, I first heated the adamantane structure up to 3000 K, equilibrated it at this temperature, and then cooled it down to 0 K.

Surprisingly, the adamantane remained stable at 3000 K, which was not expected. After some reading, I found that this might be due to the flying ice cube problem, where energy equipartition is violated in simulations using certain thermostats, leading to non-physical results.

To address this, some sources suggest applying the command
fix all momentum 1 linear 1 1 1 angular rescale
to remove accumulated center-of-mass motion. After adding this command, the results changed adamantane was no longer stable.

My question is: Is this the correct way to resolve the issue? Should this momentum fix be applied whenever using the Berendsen thermostat, or only during the relaxation (constant-temperature, canonical ensemble) phase? Since this command can significantly alter the results, I want to be sure about the proper usage

2

You are right, so a safer way sometimes is to only use this command during equilibration, not production, and also to not invoke that command every timestep, like that :

fix all momentum 1000 linear 1 1 1 angular rescale

Best,
Simon

3

Just saying, if you worry about the equipartition of energy, then the first thing to consider is probably switching to a different thermostat, e.g. stochastic velocity rescaling (fix temp/csvr). If you found difficulty in dealing with 0K (where some thermostats do not work properly) then you can change the initial/final temperature to e.g. 1K or 10K (and give initial velocities with the velocity command), which should be essentially same for most purposes. You will likely still need fix momentum due to the numerical noise.

Another thing to notice is that by default the temperature calculation in LAMMPS only removes translational DOF (i.e. 3N-3), which could be a problem for single molecules whose global rotation is also removed (so the “real” DOF is 3N-6). See the documentation of compute temp and compute_modify for how to change this behavior.