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

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

glotzerlab / hoomd-blue