Confusion regarding hoomd.md.nlist.Cell

[mgaimann]

Hello, using HOOMD version 4.7, I would like to compute a custom force by using hoomd.md.nlist.Cell instead of a naive implementation that I have been using for rapid prototyping purposes:

class MyCustomForce(hoomd.md.force.Custom):
    def __init__(
            self,
            nlist: hoomd.md.nlist.NeighborList,
            strength: float,
            r_cut: float,
            box_size: float,
            device: str = "cpu",
    ):
        self.nlist = nlist
        self.strength = strength
        self.r_cut = r_cut
        self.box_size = box_size
        self.device = device
        super().__init__(aniso=False)  # no torque

          with self.cpu_local_force_arrays as force_arrays:
              with self._state.cpu_local_snapshot as snap:
                  with self.nlist.cpu_local_nlist_arrays as nlist_arrays:
                      # agent_mask = snap.particles.typeid == 0
                      positions = snap.particles.position._coerce_to_ndarray()[:-1]

                      pair_list = self.nlist.pair_list
                      # pair_list = nlist_arrays._nlist_obj.local_pair_list
                      pair_list_naive = compute_pair_list(self.box_size, positions, 1.0)

                      # driver_id = np.where(snap.particles.typeid == 1)[0]
                      driver_id = len(positions)  # driver is last
                      pair_list_driver_filtered = pair_list[~np.any(pair_list == driver_id, axis=1)]
                      assert compare_arrays_fully_sorted(pair_list_driver_filtered, pair_list_naive)

                      force_pl = get_custom_force_using_pair_list(
                          pair_list_driver_filtered,
                          positions,
                          self.strength,
                          self.r_cut,
                          self.box_size,
                      )

                      force_naive = get_custom_force_naive(
                          positions,
                          self.strength,
                          self.r_cut,
                          self.box_size,
                      )
                      force_arrays.force[:-1] = force_pl

I have a special particle (the "driver" particle with type 1) that is externally controlled and should be located at the end of the positions array, and N agent particles with type 0.

I am using a hoomd.md.nlist.Cell object as follows:

    nlist = hoomd.md.nlist.Cell(
        buffer=0.0,
        check_dist=False,
        # deterministic=True,
        rebuild_check_delay=1,
        default_r_cut=1.0,
    )

and added it to the operations array: simulation.operations.computes.append(nlist).

Essentially, I am looking at the difference between the HOOMD computed pair_list and the pair list that I am computing myself for unit testing purposes here (pair_list_naive).
I am removing the driver particle in the force computation here since the custom force should only compute interactions between agents (of type 0) (-> pair_list_driver_filtered).

In my unit tests, I can confirm that in the first step, pair_list_driver_filtered and pair_list_naive are the same (checked using the assertion assert compare_arrays_fully_sorted(pair_list_driver_filtered, pair_list_naive). In the second step, however, this assertion fails.

What confuses me is that

• nlist_arrays._nlist_obj.shortest_rebuild is set to 100 (I would expect 1)
• nlist_arrays._nlist_obj.num_builds is set to 0, even after simulation.run() has been called for a second time.

Am I using the access via self.nlist.pair_list correctly here, and I don't actually need with self.nlist.cpu_local_nlist_arrays as nlist_arrays? Or should I make use of the latter and proceed as outlined in this discussion?

Thank you for your help!

[joaander]

You do not need the with self.nlist.cpu_local_nlist_arrays as nlist_arrays to use the .pair_list attribute.

The neighbor list is not building because you disabled distance checks:

hoomd-blue/hoomd/md/NeighborList.cc

Lines 1252 to 1256 in 6df3e7d

	if (!m_force_update && !shouldCheckDistance(timestep))
	{
	m_last_check_result = false;
	return false;
	}

Set check_dist=True (the default). With buffer == 0.0, it will not actually check the distance:

hoomd-blue/hoomd/md/NeighborList.cc

Lines 1284 to 1314 in 6df3e7d

	// not a forced update, perform the distance check to determine
	// if the list needs to be updated - no dist check needed if r_buff is tiny
	// it also needs to be updated if m_rebuild_check_delay is 0, or the check period is hit
	// when distance checks are disabled
	if (m_r_buff < 1e-6
	|| (!m_dist_check
	&& (m_rebuild_check_delay == 0
	|| (m_rebuild_check_delay > 1
	&& timestep == (m_last_updated_tstep + m_rebuild_check_delay)))))
	{
	result = true;
	}
	else
	{
	result = distanceCheck(timestep);
	}
	if (result)
	{
	// record update histogram - but only if the period is positive
	if (timestep > m_last_updated_tstep)
	{
	uint64_t period = timestep - m_last_updated_tstep;
	if (period >= m_update_periods.size())
	period = m_update_periods.size() - 1;
	m_update_periods[period]++;
	}
	m_last_updated_tstep = timestep;
	m_updates += 1;
	}

[joaander]

I don't know how you compute your forces, so I can only guess. pair_list should list each neighboring pair only once. Do you account for both i->j and j->i interactions resulting from each pair?

The internal neighbor list is stored by index, not tag. On the GPU, when deterministic=False, those indices can be in any order based on the order of thread execution. Floating point math is not associative, so the same initial condition will lead to a different trajectory when forces are summed in a different order. Setting deterministic=True is only one step to making binary reproducible simulations on the GPU. You also need to set all the autotuned kernel parameters to fixed values and possibly take other steps as well.

pair_list is a mapping from those indices to tags coupled with a filter that removes duplicates. Depending on the current particle index order, a given pair might appear in (i,j) or (j,i) ordering.

[mgaimann]

The deterministic setting does also affect single-threaded CPU computations, right?
When using a single thread, the particle order should not change, i.e. indices would be equal to tags, right?

The force implementation is quite simple, and should be reciprocal for a pair (i,j):

@numba.njit(fastmath=True)
def get_lymburn_repulsion_force_using_pair_list(
        pair_list: np.ndarray,
        positions: np.ndarray,
        strength: float,
        r_cut: float,
        box_size: float
) -> np.ndarray:
    """Calculate repulsive forces between neighboring agents using a pair list.

    Computes forces that repel agents from each other based on their positions,
    using a precomputed neighbor list for efficiency. The force magnitude follows
    a 1/r potential, scaled by the specified strength.

    This is more efficient than the naive O(n²) implementation when using a
    neighbor list that has already been computed.

    Args:
        pair_list (np.ndarray): List of neighboring particle pairs, shape (n_pairs, 2).
        positions (np.ndarray): Positions of all agents, shape (n_agents, n_dim).
        strength (float): Strength of the repulsion force.
        r_cut (float): Cutoff radius for interactions.
        box_size (float): Size of the simulation box.

    Returns:
        np.ndarray: Repulsion forces for all agents, shape (n_agents, n_dim).
    """
    forces = np.zeros(positions.shape, dtype=np.float64)
    r_cut_sq = r_cut * r_cut

    for i in range(pair_list.shape[0]):
        idx_a = pair_list[i, 0]
        idx_b = pair_list[i, 1]

        # Calculate displacement vector from agent i to agent j
        distance_vector = positions[idx_b] - positions[idx_a]
        distance_vector = apply_pbcs_to_vector(distance_vector, box_size)
        distance_sq = np.sum(distance_vector * distance_vector)

        # Only apply forces if agents are within cutoff radius
        if distance_sq < r_cut_sq:
            distance = np.sqrt(distance_sq)
            distance_unit_vector = distance_vector / distance

            # Calculate force using 1/r potential
            this_force = strength / distance * distance_unit_vector

            # Apply equal and opposite forces to both agents
            forces[idx_a] -= this_force
            forces[idx_b] += this_force

    return forces

[mgaimann]

print(snap.particles.tag) during the simulation (in the set_forces method) yields the following:

HOOMDArray([  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107
 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
 198 199 200])
HOOMDArray([ 83 199 131  17  63   8 107 102 197 161 166  99 195 134  41   7 156 148
  73  12  77  46 153  66  13 124  75  30  76  90 165 112  96 180  32  86
 117  27  23  54  98 137 187 130  52 115 196 168 162  15 160 111 178 119
 133  38 164  47 140  39 100  24  40  45 108 151 136  88 142 198 125  20
 145 170 118 101 179  95  79 152 158   2  62 176  71   3 200  55  16 189
 171  11  29 144 147 121 116 173  85  14  43  93 177 104 181   0 128 123
  42   6  67  25  68  44  91 174  48 169 109  53   1  18 141  65 194 110
 139 175  69  28  31  64 106  22 132 157  50 103 163 113  80  21  56 185
  89  74 122  81 183  60 154  84  94  51  37  97 186 126 155   9 167   5
  87 188 172  72  78  57 193 192  33 146  49  34 190 159 120  58  59 182
  70 127 184 150  92 138  19 114 191  61 129   4 149  82  35  36 135  10
 143  26 105])
HOOMDArray([ 83 199 131  17  63   8 107 102 197 161 166  99 195 134  41   7 156 148
  73  12  77  46 153  66  13 124  75  30  76  90 165 112  96 180  32  86
 117  27  23  54  98 137 187 130  52 115 196 168 162  15 160 111 178 119
 133  38 164  47 140  39 100  24  40  45 108 151 136  88 142 198 125  20
 145 170 118 101 179  95  79 152 158   2  62 176  71   3 200  55  16 189
 171  11  29 144 147 121 116 173  85  14  43  93 177 104 181   0 128 123
  42   6  67  25  68  44  91 174  48 169 109  53   1  18 141  65 194 110
 139 175  69  28  31  64 106  22 132 157  50 103 163 113  80  21  56 185
  89  74 122  81 183  60 154  84  94  51  37  97 186 126 155   9 167   5
  87 188 172  72  78  57 193 192  33 146  49  34 190 159 120  58  59 182
  70 127 184 150  92 138  19 114 191  61 129   4 149  82  35  36 135  10
 143  26 105])
...

and the tags then stay like that during the unit test.

Just to double-check: Does this mean the position of the particle with tag 83 is stored at index 0 now?
And consequently, a pair (83,199) in the pair_list corresponds to the indices in the tag array, i.e. index 0 and 1. Meaning that to compute the force between these two particles, I have to retrieve the position of particles 83 and 199 at indices 0 and 1, and then write the results into the force array at indices 0 and 1.

Which could explain why my code works on the first step (when tags happen to equal indices) but breaks afterwards?

[joaander]

pair_list stores pairs of tags. You are using cpu_local_force_arrays and cpu_local_snapshot which are accessed by index. To get the current indices, use snap.particles.rtag[tag] instead of the tags directly on these lines:

        idx_a = pair_list[i, 0]
        idx_b = pair_list[i, 1]

It doesn't help you right now, but I have designed hoomd-rs to avoid all of this needless complexity. In hoomd-rs, you only need to write your custom pair force routine as a function of two site's properties --- all the index, tag, ghost, and spatial data structure complexities are already implemented elsewhere. Your custom potential is also compiled to machine code and will therefore run 100's to 1000's of times faster than interpreted Python code. I hope to release hoomd-rs open source before the end of the year. The Monte Carlo implementation is nearly feature-complete, and I'll be sharing benchmarks soon. We have a working prototype for MD, but there is still a lot of work to make it production-ready.

[mgaimann]

Okay, thank you very much! Also very exciting to see that hooms-rs is coming :)

One last thing is now not working: when comparing the pair list that I compute naively and the pair list that I retrieve from self.nlist.pair_list, I noticed that at the beginning of the simulation, the self.nlist.num_builds counter remains at zero, until timestep = 3 (passed in set_forces). There is exactly one mismatch between the two pair lists, which always happens for timestep = 1.

TIME:  0
NUM BUILDS:  0
Arrays are equal (within tolerance)

TIME:  1
NUM BUILDS:  0
Rows in arr1 but not in arr2: 11
[[ 29  42]
 [ 31 123]
 [ 33  37]
 [ 39 101]
 [ 46  90]
 [ 56 122]
 [ 60  69]
 [ 66  98]
 [110 174]
 [139 174]
 [174 175]]

Rows in arr2 but not in arr1: 12
[[  3  88]
 [  4 143]
 [ 18 177]
 [ 35 195]
 [ 36  58]
 [ 42 106]
 [ 81 132]
 [ 97 122]
 [120 155]
 [128 177]
 [135 197]
 [148 180]]
Shape mismatch: (636, 2) vs (637, 2)

TIME:  2
NUM BUILDS:  0
Arrays are equal (within tolerance)

TIME:  3
NUM BUILDS:  1
Arrays are equal (within tolerance)

TIME:  4
NUM BUILDS:  2
Arrays are equal (within tolerance)

...

Can you help me understand, please, why the num_builds counter remains at 0 until the third time step? I would expect the num_builds increment to match the timestep increment. Is there something special happening with the pair list at the beginning of the simulation that I am missing, maybe?

[joaander]

I do not have an explanation for that behavior. You are welcome to debug the needsUpdating method I linked to above and submit a pull request that fixes it.

buildNlist is only called when needsUpdating returns true and the only code path that might return true without incrementing m_updates is this:

hoomd-blue/hoomd/md/NeighborList.cc

Lines 1239 to 1248 in dc2b62a

	if (m_last_checked_tstep == timestep)
	{
	if (m_force_update)
	{
	// force update is counted only once per time step
	m_force_update = false;
	return true;
	}
	return m_last_check_result;
	}

It also appears that num_builds might overcount the force builds:

hoomd-blue/hoomd/md/NeighborList.h

Lines 409 to 413 in dc2b62a

	//! Get the number of updates
	virtual uint64_t getNumUpdates()
	{
	return m_updates + m_forced_updates;
	}

[mphoward]

The particles will be resorted periodically throughout the simulation, even on a single CPU:

https://hoomd-blue.readthedocs.io/en/v5.4.0/hoomd/tune/particlesorter.html#hoomd.tune.ParticleSorter