Empirically determining thermostat damping factor

Question

Is there a way to empirically determine an appropriate thermostat damping factor given a timestep size and a numerical integration method? For example, I would like one for the surface and the molecule, and I will describe more below.

I am computing a set of simulations in LAMMPS, using ReaxFF, tracking the path of a single (~4 atom) gas molecule as it interacts with a (~600 atom) surface.

Two thermostats are imposed on the system, one for the surface, one for the molecule.

fix surfaceFix  surface  temp/csvr ${surfTemp} ${surfTemp} ${surfaceDampingFactor}
fix moleculeFix molecule temp/csvr ${molTemp}  ${molTemp}  ${moleculeDampingFactor}

I know intuitively that the surface temperature should be more tightly regulated than the temperature of the molecule, since the molecule trajectory is a dependent variable in the experiment. Yet I don't have a way of determining the value of each damping factor.

For instance with the molecule, errors are introduced in the limit of both small and large damping factors:

• In the case of weak thermostatting, the force calculations can diverge and we introduce an error due to the precision limits of numerical integration.
• In the case of strong thermostatting, we introduce an error in the energy transferred between the molecule and the surface.

Both of these errors could cause the effect of "hopping": where the molecule is provided with a non-physical velocity that causes it to jump into or out of locations that would otherwise be impossible.

There must be a connection between the timestep size + integration method, and the damping factor needed to correct for the deviation in forces that these produce. Naturally the damping required would also depend on the mechanics of the thermostat, and its effect on different degrees of freedom, but an answer addressing any thermostat is of interest.

Answer 1

To the best of my knowledge, there is no universal or even near-universal heuristic for what the time constant should be. The only rule of thumb I have encountered is that the time constant should be around 10x-1000x larger than the timestep of the simulation. Beyond that point, I would treat the damping constant as a hyperparameter, do a sensitivity analysis, and pick a value where changing the damping by e.g. 10x in either direction does not affect the results. In other words, if the simulation outcome is robust to the chosen value of damping, there is little for a picky reviewer to pick. In my own simulations, thermostatting between 10fs-1ps for simulation runs with a 1-fs timestep has never made a difference, unless there was energy getting added to the simulation mid-way through (e.g. from an applied electric field). But for those cases, NVE was more appropriate in the first place.

For the specific system in question, a gas molecule and a surface, heat exchange between the molecule and the surface physically happens via some density of vibrational states. For example, when the molecule sits on a surface, it could transfer energy back and forth with the surface’s phonon modes. As the surface has a lot of atoms, its density of vibrational states is large, and could very well be treated as continuous. In the simplest terms, this is how equipartitioning physically happens. If the phonon-phonon couplings or vibrational coherence times (dephrasing rates) for the molecule and the surface are known, those could offer physics-based starting points for the optimization of the damping hyperparameter. The molecular vibrations are typically more singular, but the molecule also has (typically low-energy and low-frequency, i.e. slow) rotations and translations.

The interplay between these slow degrees of freedom and the effective potential that keeps the molecule on the surface determines binding and unbinding. This is the effective potential in which the molecule moves as a whole. The correspondence of the simulation velocity and the accuracy with which the simulation describes the molecule’s positions (and the “hopping” in the question) depends on accuracy of thermalization on this relatively slow timescale. With translations and rotations being on the timescale of a terahertz (≈1 ps), and phonon thermalization typically faster, this should be the slowest timescale in the system.

Now, physically speaking, how would a hop (as described in the question) happen? The molecule could move far enough to unbind or switch position. Practically, this should not happen within one timestep of the simulation. Assuming that the timestep of the simulation is small enough that translation takes a large number of time steps, then the continued presence of extra energy in the translational degree of freedom is either physical, or will be taken out by the thermostat. Ultimately, judging whether a binding or unbinding event is physical or unphysical is going to be up to the human behind the simulation.

A second practical piece of advice I could offer is to verify whether exempting the center-of-mass motions of the molecule and/or the surface makes any difference. Sometimes that is sufficient to prevent unphysical accumulation of energy in the translational degrees of freedom, and I would consider it for simulations where one part of the system is expected to meaningfully translate past another, such as the molecule translating past the surface in question. E.g. for LAMMPS+Nosé:

fix molVibr molecule temp/com ${molTemp}  ${molTemp}  ${moleculeDampingFactor}