glotzerlab / hoomd-blue

Molecular dynamics and Monte Carlo soft matter simulation on GPUs.
http://glotzerlab.engin.umich.edu/hoomd-blue
BSD 3-Clause "New" or "Revised" License
323 stars 127 forks source link

Failing to assign `shape` to `hoomd.md.pair.aniso.ALJ` for every particle type causes simulation to error incorrectly #1805

Closed josephburkhart closed 1 month ago

josephburkhart commented 1 month ago

Description

When using an ALJ force in a binary system (two types of particles), it is necessary to set the shape parameter for both types of particles. This requirement is not clearly communicated in the documentation. (Or at least, it wasn't clear to me.) When the user fails to set shape for both types, the simulation throws an error that is either incorrect or extremely unclear.

The MWE below demonstrates the error I'm referring to. Uncomment alj.shape["B"] = dict(vertices=vertices_a, faces=faces_a, rounding_radii=0.0) to prevent the error.

Script

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 ALJ force
vertices_a = [[-0.5, 0, 0], [0, 0.5, 0], [0.5, 0, 0], [0, -0.5, 0],
              [0, 0, -0.5], [0, 0, 0.5]]
faces_a = [[0, 1, 4], [1, 2, 4], [2, 3, 4], [3, 1, 4],
           [0, 1, 5], [1, 2, 5], [2, 3, 5], [3, 1, 5]]

cell = hoomd.md.nlist.Cell(0.5, exclusions=["body"])
alj = hoomd.md.pair.aniso.ALJ(nlist=cell, default_r_cut=2.5)

alj.shape["A"] = dict(vertices=vertices_a, faces=faces_a, rounding_radii=0.0)
# alj.shape["B"] = dict(vertices=vertices_a, faces=faces_a, rounding_radii=0.0) # <-- uncomment this line to prevent the error

alj.params[("A", "A")] = dict(epsilon=1.0, sigma_i=1, sigma_j=1, alpha=0)
alj.params[("B", ["A", "B"])] = dict(epsilon=0.0, sigma_i=1, sigma_j=1, alpha=0)

alj.r_cut[("A", "A")] = max(2**(1/6)*2*1, 1 + (2**(1/6)*0.15*2*1))

integrator.forces.append(alj)

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

Input files

N/A

Output

Traceback (most recent call last):
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operation.py", line 350, in _apply_typeparam_dict
    typeparam._attach(cpp_obj, simulation.state)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/typeparam.py", line 171, in _attach
    self.param_dict._attach(cpp_obj, self.name,
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/parameterdicts.py", line 518, in _attach
    self._single_setitem(key, parameters[key])
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/parameterdicts.py", line 461, in _single_setitem
    getattr(self._cpp_obj, self._setter)(key, item)
ValueError: dictionary update sequence element #0 has length 3; 2 is required

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operation.py", line 316, in _attach
    self._apply_typeparam_dict(self._cpp_obj, self._simulation)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operation.py", line 352, in _apply_typeparam_dict
    raise err.__class__(
ValueError: For <class 'hoomd.md.pair.aniso.ALJ'> in TypeParameter shape dictionary update sequence element #0 has length 3; 2 is required

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/joseph/GlotzerGroup/pentagonal-nanorods/mwe1.py", line 71, in <module>
    simulation.run(1000)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/simulation.py", line 555, in run
    self.operations._schedule()
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operations.py", line 211, in _schedule
    self.integrator._attach(sim)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operation.py", line 310, in _attach
    self._attach_hook()
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/md/integrate.py", line 307, in _attach_hook
    super()._attach_hook()
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/md/integrate.py", line 48, in _attach_hook
    self._forces._sync(self._simulation, self._cpp_obj.forces)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/syncedlist.py", line 246, in _sync
    raise err
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/syncedlist.py", line 242, in _sync
    self._register_item(item)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/data/syncedlist.py", line 190, in _register_item
    value._attach(self._simulation)
  File "/home/joseph/miniconda3/envs/glotzergroup/lib/python3.11/site-packages/hoomd/operation.py", line 318, in _attach
    raise type(err)(
ValueError: Error applying parameters for object of type <class 'hoomd.md.pair.aniso.ALJ'>.

Expected output

Instead of TypeParameterDict._single_setitem() throwing ValueError: dictionary update sequence element #0 has length 3; 2 is required, I would expect some class higher up to throw an error saying something like shape dictionary has not been provided for particle type 'B'.

Platform

CPU, Linux

Installation method

Conda-forge package

HOOMD-blue version

4.1.0

Python version

3.11.7

joaander commented 1 month ago

TypeParameter documents that parameters for all types (or type pairs) must be set in an important block: https://hoomd-blue.readthedocs.io/en/v4.7.0/module-hoomd-data.html#hoomd.data.typeparam.TypeParameter. ALJ.shape links directly to the TypeParameter documentation, as to all other attributes that are type parameters. How could this be documented more clearly?

And yes, there should be a user-friendly error message for this. I will investigate why it is not issued.

joaander commented 1 month ago

Also, I think you would want to set empty lists for vertices/faces and a r_cut=0 to implement a non-interacting type.

alj.shape["B"] = dict(vertices=[], faces=[], rounding_radii=0.0)
alj.r_cut[("B", ["A", "B"])] = 0