hoomd.md should warn or error on submission of certain invalid parameters for pair forces

[josephburkhart]

Description

When assigning parameters for pairwise forces, it is currently possible to assign values that will cause critical errors such as divide-by-zero in the simulation's execution. It would be beneficial to detect these cases prior to simulation execution and either warn the user about them or throw an error. An example is given below.

In the following MWE, I create a binary ("A" and "B") system of particles that interact with Gaussian pairwise forces. "A" particles are cores, "B" particles are constituents. I only want interactions to happen for B particles, so I set the parameters for A-A and A-B pairs to zero. This should cause an error, since σ is in the denominator in the force's definition.

Indeed, the simulation does throw an error: RuntimeError: Particle with unique tag 17 has NaN for its position.. Within this MWE, the cause of the problem is pretty clear, but when working with systems that are more complex, this problem can be much more difficult to diagnose - this is evidenced by the fact that it took me, @tcmoore3 and @joaander multiple meetings and many hours to figure out that this was causing fatal errors in my simulations.

import hoomd
import itertools
import numpy
import gsd.hoomd

# Initialize Frame
frame = gsd.hoomd.Frame()
frame.particles.types = ["A", "B"]

positions = list(itertools.product(numpy.linspace(-2, 2, 2, endpoint=False), repeat=3))
frame.particles.N = len(positions)
frame.particles.position = positions
frame.particles.typeid = [0] * frame.particles.N
frame.configuration.box = [8, 8, 8, 0, 0, 0]
frame.particles.mass = [1] * frame.particles.N
frame.particles.moment_inertia = [1, 1, 1] * frame.particles.N
frame.particles.orientation = [(1, 0, 0, 0)] * frame.particles.N

# Initialize Simulation
simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)
simulation.create_state_from_snapshot(frame)

# Create rigid bodies
rigid = hoomd.md.constrain.Rigid()
rigid.body["A"] = {
    "constituent_types": ["B", "B"],
    "positions": [[-0.5, 0, 0], [0.5, 0, 0]],
    "orientations": [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)]
}
rigid.create_bodies(simulation.state)

# Create integrator
integrator = hoomd.md.Integrator(dt=0.0001, integrate_rotational_dof=True)
simulation.operations.integrator = integrator
integrator.rigid = rigid

method = hoomd.md.methods.ConstantVolume(
    filter=hoomd.filter.Rigid(("center", "free")),
    thermostat=hoomd.md.methods.thermostats.Bussi(kT=1.0)
)

integrator.methods.append(method)

# Create, configure, and append the Gaussian force
cell = hoomd.md.nlist.Cell(0.5, exclusions=["body"])
gaussian = hoomd.md.pair.Gaussian(nlist=cell, default_r_cut=2.5)

gaussian.params[("B", "B")] = dict(epsilon=1.0, sigma=1.0)
gaussian.params[("A", ["A", "B"])] = dict(epsilon=0, sigma=0) # <-- change sigma to 1.0 to prevent the error 

integrator.forces.append(gaussian)

# Thermalize and run
simulation.state.thermalize_particle_momenta(
    filter=hoomd.filter.Rigid(("center", "free")),
    kT=1.0
)
simulation.run(1000)

Proposed solution

I am not yet familiar enough with Hoomd's architecture to directly propose a solution, but I would expect that the value check should happen in the same place where types are checked at simulation execution. With some assistance/direction from @joaander, this issue may offer a good opportunity for me to make my first pull request in this repo.

Additional context

No response

[joaander]

For future reference, there is make_example_simulation which can make these MWE's easier to write.

HOOMD has rich Python abstractions for validating parameters. Here is an example:

hoomd-blue/hoomd/md/minimize/fire.py

Lines 220 to 224 in 04c6c71

	pdict = ParameterDict(
	dt=float(dt),
	integrate_rotational_dof=bool(integrate_rotational_dof),
	min_steps_adapt=OnlyTypes(int, preprocess=positive_real),
	finc_dt=float(finc_dt),

hoomd-blue/hoomd/data/typeconverter.py

Lines 82 to 91 in 04c6c71

	def positive_real(number):
	"""Ensure that a value is positive."""
	try:
	float_number = float(number)
	except Exception as err:
	raise TypeConversionError(
	f"{number} not convertible to float.") from err
	if float_number <= 0:
	raise TypeConversionError("Expected a number greater than zero.")
	return float_number

You can apply positive_real in the same way to sigma in the Gaussian and ExpandedGaussian potentials here:

hoomd-blue/hoomd/md/pair/pair.py

Lines 274 to 279 in 04c6c71

	def __init__(self, nlist, default_r_cut=None, default_r_on=0., mode='none'):
	super().__init__(nlist, default_r_cut, default_r_on, mode)
	params = TypeParameter(
	'params', 'particle_types',
	TypeParameterDict(epsilon=float, sigma=float, len_keys=2))
	self._add_typeparam(params)

hoomd-blue/hoomd/md/pair/pair.py

Lines 331 to 337 in 04c6c71

	params = TypeParameter(
	'params', 'particle_types',
	TypeParameterDict(epsilon=float,
	sigma=float,
	delta=float,
	len_keys=2))
	self._add_typeparam(params)