Is there a way to determine an appropriate thermostat damping factor given a timestep size and a numerical integration method?

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}

Question: How do I empirically determine a damping factor for the thermostats of the surface and the molecule?

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.

  • 3
    $\begingroup$ For a start, I'd strongly suggest not using Berendsen thermostats and barostats for production runs, because they don't sample rigourously from the Boltzmann distribution. I'd say it's best to use Langevin dynamics for that, and if you really need continuous/"deterministic" trajectories, then I'd use Nosé-Hoover chains instead. $\endgroup$
    – Godzilla
    Sep 8 '20 at 23:22
  • $\begingroup$ @Godzilla Thanks. I hadn't updated the question, but for the molecule, I have since planned to substitute the berendsen thermostat out for temp/csvr instead, after reading this doi.org/10.1063/1.2408420. As I understand it, temp/csvr addresses canonical resampling, but should also provide deterministic trajectories since it won't change c.o.m motion in the molecule. I'll compare this against the Nosé-Hoover chain, and consider using a different thermostat for the surface. $\endgroup$
    – geo
    Sep 9 '20 at 8:21
  • $\begingroup$ @Godzilla how are we doing with this question? Has the user provided sufficient information for the question to be answered by somebody? Geo: have you done the comparison against the Nosé-Hoover chain as you said in the last comment? I'd really like to get this question out of the unanswered queue if possible, since if that queue gets too big, no one will see the other unanswered questions. Are you aware of the LAMMPS chat room? chat.stackexchange.com/rooms/109805/lammps Enter and say "hello" at least once so that we remember you if anything useful for you comes up! $\endgroup$ Oct 18 '20 at 18:32
  • $\begingroup$ @NikeDattani Well I think this is one of those questions that we need a follow-up on, I am getting a déjà vu here :) $\endgroup$
    – Godzilla
    Oct 20 '20 at 16:03
  • 1
    $\begingroup$ Related: mattermodeling.stackexchange.com/q/4245/5 $\endgroup$ Mar 10 at 16:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.