Add thermostat/integrator compatibility checks

[jngrad]

Fixes #3850

Description of changes:

• add missing thermostat combination decorators
• document policy on allowed thermostat combinations
• prevent usage of incompatible thermostat/integrator combinations

[KaiSzuttor]

is no combination of SD with any other thermostat allowed?

[RudolfWeeber]

is no combination of SD with any other thermostat allowed?

I think so.

[KaiSzuttor]

is no combination of SD with any other thermostat allowed?

I think so.

and why?

[mkuron]

That's how SD works. It puts Brownian fluctuations on all particles because it's not only a hydrodynamics solver, but also a momentum-conserving thermostat. Just like you wouldn't put a Langevin thermostat on particles coupled to a thermalized LB, you shouldn't combine SD with other thermostats.

[KaiSzuttor]

yeah, but you can for example combine LB + DPD because both conserve momentum and add fluctuations on the particles

[RudolfWeeber]

The SD integrator choses velocities randomly from the Maxwell distribution. I don’t see how this interacts safely with anything else which is based on an equation of motion including inertia. That aside: safety first. If someone convinces us of a safe combination, it can be enabled then.

[KaiSzuttor]

well, it should be documented then. we have to tell the user why we do not allow the combination (even if the reason is that we don't know if works).

[biermanncarl]

Thanks JN. This fixes that both thermostats could be enabled, which is not known to correctly work (yet). However, I still can combine the SD integrator and the Langevin thermostat (and possibly other thermostats, too). Should this be allowed?

I think that maybe, the SD integrator should only be allowed to work with the SD thermostat (and the other way round is obvious, because the SD thermostat depends on the SD intergator).

[RudolfWeeber]

Thanks JN. This fixes that both thermostats could be enabled, which is not known to correctly work (yet). However, I still can combine the SD integrator and the >Langevin thermostat (and possibly other thermostats, too). Should this be allowed?

No, it shoudn’t. That should probably be checked in on_integration_start(). The underlying issue is that in Espresso, thermostats and integrators are independent, whereas that’s not so in physics…

[mkuron]

I still can combine the SD integrator and the Langevin thermostat

It would be an odd combination, but probably not outright wrong. It's somewhat like combining an unthermalized LB with a Langevin thermostat -- do we currently disallow that?

[jngrad]

AFAIK the integrator and thermostat classes don't know about each other in the python interface.

[jngrad]

well, it should be documented then. we have to tell the user why we do not allow the combination (even if the reason is that we don't know if works).

Currently none of the combination restrictions are documented. This is essentially an allowlist of combinations that are known to work together, which gets expanded whenever we have hard data confirming that two thermostats work together, in opposition to a denylist that would prevent specific combinations for which we can justify they produce wrong physics. I'll write something in the docs to explain this policy.

[RudolfWeeber]

AFAIK the integrator and thermostat classes don't know about each other in the python interface.

Yes. It needs to be done in the core, e.g., in on_integration_start(). Once that check is in, this PR is ready, imo.

[jngrad]

There are more issues:

• we can initialize thermostats without a seed to crash the script (the error message from the un-initialized boost::optional is displayed as a char pointer)

  import espressomd
  system = espressomd.System(box_l=[1., 1., 1.])
  system.time_step = 0.01
  system.cell_system.skin = 0.4
  system.part.add(pos=(0, 0, 0), ext_force=[1,0,0])
  system.integrator.set_brownian_dynamics()
  system.integrator.run(1)

  output:

  terminate called after throwing an instance of 'char const*'
  

• we can use a thermostat after it has been "destructed" in the core with turn_off(), leading to divisions by zero:

  import espressomd
  system = espressomd.System(box_l=[1., 1., 1.])
  system.time_step = 0.01
  system.cell_system.skin = 0.4
  system.part.add(pos=(0, 0, 0), ext_force=[1,0,0])
  system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=42)
  system.thermostat.turn_off()
  system.integrator.set_brownian_dynamics()
  system.integrator.run(1)
  print(system.part[:].pos)

  output:

  [[inf nan nan]]