Neighbor Lists Access for Custom Pair Force

[mgaimann]

Hi everyone, I am very new to hoomd-blue and have the following question:

When trying to implement a custom pair force using hoomd.md.force.Custom and hoomd.md.nlist.Cell I am receiving the error DataAccessError: The property cpu_local_nlist_arrays is unavailable until the simulation runs for 0 or more steps..
I found the problem to be related to this issue and this PR.

Minimal working example:

import hoomd
import gsd.hoomd
import numpy as np

# integrator
integrator = hoomd.md.Integrator(
    dt=0.001, integrate_rotational_dof=False
)
nvt = hoomd.md.methods.ConstantVolume(  # implements Newton's law
    filter=hoomd.filter.All()
)
integrator.methods.append(nvt)

# custom force
class CustomPairForce(hoomd.md.force.Custom):
    def __init__(self, nlist, strength):
        self.nlist = nlist
        self.strength = strength
        super().__init__(aniso=False)  # no torque

    def set_forces(self, timestep):
        with self.cpu_local_force_arrays as force_arrays:
            with self.nlist.cpu_local_nlist_arrays as nlist_arrays:
                pass
                # use neighbor information to update forces

cell = hoomd.md.nlist.Cell(buffer=0.4)
custom_pair_force = CustomPairForce(
    nlist=cell, strength=1.0
)
integrator.forces.append(custom_pair_force)

# initial frame
config = {"n_agents": 100, "box_size": 16.0}
frame = gsd.hoomd.Frame()
frame.configuration.box = [config['box_size'], config['box_size'], config['box_size'], 0, 0, 0]

# add agents
frame.particles.N = config['n_agents']
frame.particles.position = np.random.uniform(-config['box_size'] / 2, config['box_size'] / 2, (config['n_agents'], 3))
frame.particles.typeid = [0] * config['n_agents']

# simulation
cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)
simulation.create_state_from_gsd(filename='initial_frame.gsd')
simulation.operations.integrator = integrator
simulation.run(0)

I found that Cell is a subclass of NeighborList, which is a Compute subclass. Could the order of execution of neighbor computation, force computation, and integration matter here somehow?
I was expecting that run(0) would enable the neighbor computation without changing the system state following this tutorial page.

Thank you for your help!

[joaander]

Add:

simulation.operations.computes.append(cell)

before calling simulation.run.

#1820 changes the error message to make this requirement more clear.