From 047da30ec665f455a226be181a613b926c7e17d4 Mon Sep 17 00:00:00 2001 From: clpetix Date: Sat, 10 Feb 2024 12:11:29 -0600 Subject: [PATCH 01/10] Filter neighbor list. --- src/relentless/simulate/hoomd.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 276e587e..af268d53 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -70,6 +70,7 @@ def _call_v3(self, sim): # parse masses by type snap = sim["engine"]["_hoomd"].state.get_snapshot() sim.masses = self._get_masses_from_snapshot(sim, snap) + sim.bonds = self._get_bonds_from_snapshot(sim, snap) self._assert_dimension_safe(sim, snap) # create the potentials, defer attaching until later @@ -159,6 +160,9 @@ def _get_masses_from_snapshot(self, sim, snap): masses.update(masses_) return masses + def _get_bonds_from_snapshot(self, sim, snap): + return snap.bonds.group + def _assert_dimension_safe(self, sim, snap): if sim.dimension == 3: dim_safe = True @@ -1228,11 +1232,22 @@ def act(self, timestep): for i, j in self._rdf: query_args = dict(self._rdf[i, j].default_query_args) query_args.update(exclude_ii=(i == j)) + neighbors = ( + aabbs[j] + .query( + snap.particles.position[type_masks[i]], query_args + ) + .toNeighborList() + ) + bonds = numpy.vstack( + [snap.bonds.group, numpy.flip(snap.bonds.group, axis=1)] + ) + filter = ~(neighbors[:][:, None] == bonds).all(-1).any(1) # resetting when the samples are zero clears the RDF self._rdf[i, j].compute( aabbs[j], snap.particles.position[type_masks[i]], - neighbors=query_args, + neighbors=neighbors.filter(filter), reset=(self.num_samples == 0), ) From b1ba2523cd6bbec6f893e1de95d35384c98e1637 Mon Sep 17 00:00:00 2001 From: clpetix Date: Mon, 12 Feb 2024 13:44:06 -0600 Subject: [PATCH 02/10] Make bonds a private variable of the initializer. Ensure bond exclusions are only enforced if bonds are not None. --- src/relentless/simulate/hoomd.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index af268d53..a693d510 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -70,9 +70,8 @@ def _call_v3(self, sim): # parse masses by type snap = sim["engine"]["_hoomd"].state.get_snapshot() sim.masses = self._get_masses_from_snapshot(sim, snap) - sim.bonds = self._get_bonds_from_snapshot(sim, snap) + sim[self]["_bonds"] = self._get_bonds_from_snapshot(sim, snap) self._assert_dimension_safe(sim, snap) - # create the potentials, defer attaching until later neighbor_list = hoomd.md.nlist.Tree(buffer=sim.potentials.pair.neighbor_buffer) pair_potential = hoomd.md.pair.Table(nlist=neighbor_list) @@ -957,6 +956,7 @@ def _pre_run_v3(self, sim, sim_op): system=sim["engine"]["_hoomd"].state, rdf_params=self._get_rdf_params(sim), constraints=self._get_constrained_quantities(sim, sim_op), + bonds=sim[sim.initializer]["_bonds"], ) sim[self]["_hoomd_thermo_callback"] = hoomd.write.CustomWriter( trigger=self.every, @@ -1141,7 +1141,14 @@ class EnsembleAverageAction(Action): flags = [Action.Flags.PRESSURE_TENSOR] def __init__( - self, types, dimension, thermo, system, rdf_params=None, constraints=None + self, + types, + dimension, + thermo, + system, + rdf_params=None, + constraints=None, + bonds=None, ): if dimension not in (2, 3): raise ValueError("Only 2 or 3 dimensions supported") @@ -1152,6 +1159,7 @@ def __init__( self.system = system self.rdf_params = rdf_params self.constraints = constraints if constraints is not None else {} + self.bonds = bonds # this method handles all the initialization self.reset() @@ -1239,10 +1247,15 @@ def act(self, timestep): ) .toNeighborList() ) - bonds = numpy.vstack( - [snap.bonds.group, numpy.flip(snap.bonds.group, axis=1)] - ) - filter = ~(neighbors[:][:, None] == bonds).all(-1).any(1) + if snap.bonds is not None: + bonds = numpy.vstack( + [self.bonds, numpy.flip(self.bonds, axis=1)] + ) + # list intersect method from: + # https://stackoverflow.com/a/67113105 + filter = ~(neighbors[:, None] == bonds).all(-1).any(1) + else: + filter = numpy.full((len(neighbors[:])), True) # resetting when the samples are zero clears the RDF self._rdf[i, j].compute( aabbs[j], From 6e2694e6c7c800b5c921828dbbf0d8b64c15fd07 Mon Sep 17 00:00:00 2001 From: clpetix Date: Sat, 20 Apr 2024 12:01:30 -0500 Subject: [PATCH 03/10] Change neighborlist method to calculate from a full neighborlist. --- src/relentless/simulate/hoomd.py | 61 +++++++++++++++++++------------- 1 file changed, 37 insertions(+), 24 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index a693d510..b2d2a171 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1221,35 +1221,49 @@ def act(self, timestep): ) box = freud.box.Box.from_box(box_array, dimensions=dimensions) - # pre build aabbs per type and count particles by type - aabbs = {} - type_masks = {} - for i in self.types: - type_masks[ - i - ] = snap.particles.typeid == snap.particles.types.index(i) - if "N" not in self.constraints: - self._N[i] += numpy.sum(type_masks[i]) - - if compute_rdf: - aabbs[i] = freud.locality.AABBQuery( - box, snap.particles.position[type_masks[i]] - ) - # then do rdfs using the AABBs + # pre build aabbs for full system + aabbs = freud.locality.AABBQuery(box, snap.particles.position) + # then do rdfs using the aabbs if compute_rdf: for i, j in self._rdf: query_args = dict(self._rdf[i, j].default_query_args) query_args.update(exclude_ii=(i == j)) - neighbors = ( - aabbs[j] - .query( - snap.particles.position[type_masks[i]], query_args + neighbors = aabbs.query( + snap.particles.position, query_args + ).toNeighborList() + type_dict = dict( + zip(snap.particles.types, snap.particles.typeid) + ) + typeid_dict = dict( + zip( + numpy.arange(0, snap.particles.N), + snap.particles.typeid, ) - .toNeighborList() ) + neighbors_typeid = numpy.zeros( + numpy.shape(neighbors[:]), dtype=int + ) + row = 0 + for n1, n2 in neighbors: + neighbors_typeid[row] = [ + typeid_dict[n1], + typeid_dict[n2], + ] + row += 1 + type_mask = numpy.logical_or( + numpy.logical_and( + numpy.isin(neighbors_typeid[:, 0], type_dict[i]), + numpy.isin(neighbors_typeid[:, 1], type_dict[j]), + ), + numpy.logical_and( + numpy.isin(neighbors_typeid[:, 1], type_dict[i]), + numpy.isin(neighbors_typeid[:, 0], type_dict[j]), + ), + ) + neighbors = neighbors.filter(type_mask) if snap.bonds is not None: bonds = numpy.vstack( - [self.bonds, numpy.flip(self.bonds, axis=1)] + [self.bonds, numpy.flip(self.bonds, axis=1)], ) # list intersect method from: # https://stackoverflow.com/a/67113105 @@ -1258,12 +1272,11 @@ def act(self, timestep): filter = numpy.full((len(neighbors[:])), True) # resetting when the samples are zero clears the RDF self._rdf[i, j].compute( - aabbs[j], - snap.particles.position[type_masks[i]], + aabbs, + snap.particles.position, neighbors=neighbors.filter(filter), reset=(self.num_samples == 0), ) - self.num_samples += 1 if compute_rdf: self.num_rdf_samples += 1 From 1e9b4c95227f632dfd3cc970e5fe82db5841af65 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Mon, 22 Apr 2024 15:34:01 -0500 Subject: [PATCH 04/10] Simplify RDF calculations, including with bond exclusions --- src/relentless/simulate/hoomd.py | 148 +++++++++++++++++-------------- 1 file changed, 79 insertions(+), 69 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index b2d2a171..3920937f 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1197,84 +1197,94 @@ def act(self, timestep): raise NotImplementedError(f"HOOMD {_hoomd_version} not supported") if mpi.world.rank == 0: - if _hoomd_version.major == 3: - dimensions = snap.configuration.dimensions - box_array = numpy.array(snap.configuration.box) - elif _hoomd_version.major == 2: - dimensions = snap.box.dimensions - box_array = numpy.array( - [ - snap.box.Lx, - snap.box.Ly, - snap.box.Lz, - snap.box.xy, - snap.box.xz, - snap.box.yz, - ] - ) - if snap.box.dimensions == 2: - box_array[2] = 0.0 - box_array[-2:] = 0.0 - else: - raise NotImplementedError( - f"HOOMD {_hoomd_version} not supported" - ) - - box = freud.box.Box.from_box(box_array, dimensions=dimensions) - # pre build aabbs for full system - aabbs = freud.locality.AABBQuery(box, snap.particles.position) - # then do rdfs using the aabbs + # compute number of particles of each type + # save type masks for use in RDF calculations if needed + type_masks = {} + for i in self.types: + type_masks[ + i + ] = snap.particles.typeid == snap.particles.types.index(i) + if "N" not in self.constraints: + self._N[i] += numpy.sum(type_masks[i]) + + # then do rdf if requested if compute_rdf: - for i, j in self._rdf: - query_args = dict(self._rdf[i, j].default_query_args) - query_args.update(exclude_ii=(i == j)) - neighbors = aabbs.query( - snap.particles.position, query_args - ).toNeighborList() - type_dict = dict( - zip(snap.particles.types, snap.particles.typeid) + if _hoomd_version.major == 3: + dimensions = snap.configuration.dimensions + box_array = numpy.array(snap.configuration.box) + elif _hoomd_version.major == 2: + dimensions = snap.box.dimensions + box_array = numpy.array( + [ + snap.box.Lx, + snap.box.Ly, + snap.box.Lz, + snap.box.xy, + snap.box.xz, + snap.box.yz, + ] ) - typeid_dict = dict( - zip( - numpy.arange(0, snap.particles.N), - snap.particles.typeid, - ) + if snap.box.dimensions == 2: + box_array[2] = 0.0 + box_array[-2:] = 0.0 + else: + raise NotImplementedError( + f"HOOMD {_hoomd_version} not supported" ) - neighbors_typeid = numpy.zeros( - numpy.shape(neighbors[:]), dtype=int + box = freud.box.Box.from_box(box_array, dimensions=dimensions) + + # build aabb of all particles and generate a parent + # neighbor list with the RDF cutoff + aabb = freud.locality.AABBQuery(box, snap.particles.position) + neighbors = aabb.query( + snap.particles.position, + dict( + mode="ball", + r_max=self.rdf_params["stop"], + exclude_ii=True, + ), + ).toNeighborList() + + # filter bonds from the neighbor list if they are present + # bond exclusions apply regardless of order, so + # consider both (i,j) and (j,i) permutations + if snap.bonds is not None and len(neighbors[:]) > 0: + bonds = numpy.vstack( + [self.bonds, numpy.flip(self.bonds, axis=1)], ) - row = 0 - for n1, n2 in neighbors: - neighbors_typeid[row] = [ - typeid_dict[n1], - typeid_dict[n2], - ] - row += 1 - type_mask = numpy.logical_or( - numpy.logical_and( - numpy.isin(neighbors_typeid[:, 0], type_dict[i]), - numpy.isin(neighbors_typeid[:, 1], type_dict[j]), - ), - numpy.logical_and( - numpy.isin(neighbors_typeid[:, 1], type_dict[i]), - numpy.isin(neighbors_typeid[:, 0], type_dict[j]), - ), + # list intersect method from: + # https://stackoverflow.com/a/67113105 + bond_exclusion_filter = ( + ~(neighbors[:, None] == bonds).all(-1).any(1) ) - neighbors = neighbors.filter(type_mask) - if snap.bonds is not None: - bonds = numpy.vstack( - [self.bonds, numpy.flip(self.bonds, axis=1)], + neighbors.filter(bond_exclusion_filter) + + for i, j in self._rdf: + # make a neighbor list of all particles of type i + # whose neighbors are type j + neighbors_ij = neighbors.copy() + filter_ij = numpy.logical_and( + type_masks[i][neighbors[:, 0]], + type_masks[j][neighbors[:, 1]], + ) + # if the types are different, also consider + # neighbors of j that are type i + if j != i: + filter_ji = numpy.logical_and( + type_masks[j][neighbors[:, 0]], + type_masks[i][neighbors[:, 1]], ) - # list intersect method from: - # https://stackoverflow.com/a/67113105 - filter = ~(neighbors[:, None] == bonds).all(-1).any(1) - else: - filter = numpy.full((len(neighbors[:])), True) + filter_ij = numpy.logical_or(filter_ij, filter_ji) + neighbors_ij.filter(filter_ij) + + # it doesn't look like it but this calculation should + # only calculate g_ij because we prefiltered the neighbors + # # resetting when the samples are zero clears the RDF self._rdf[i, j].compute( - aabbs, + aabb, snap.particles.position, - neighbors=neighbors.filter(filter), + neighbors=neighbors_ij, reset=(self.num_samples == 0), ) self.num_samples += 1 From 7fbecd5893e2feb731370e93ec6de6e3c1f15193 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Sun, 28 Apr 2024 22:36:44 -0500 Subject: [PATCH 05/10] Histogram and normalize RDFs manually --- src/relentless/simulate/hoomd.py | 113 ++++++++++++++++++++----------- 1 file changed, 72 insertions(+), 41 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 3920937f..2b89e523 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1200,12 +1200,16 @@ def act(self, timestep): # compute number of particles of each type # save type masks for use in RDF calculations if needed type_masks = {} + N = {} for i in self.types: type_masks[ i ] = snap.particles.typeid == snap.particles.types.index(i) if "N" not in self.constraints: - self._N[i] += numpy.sum(type_masks[i]) + N[i] = numpy.sum(type_masks[i]) + self._N[i] += N[i] + else: + N[i] = self.constraints["N"][i] # then do rdf if requested if compute_rdf: @@ -1233,6 +1237,10 @@ def act(self, timestep): ) box = freud.box.Box.from_box(box_array, dimensions=dimensions) + # calculate average density per type + for i in self.types: + self._rdf_density[i] += N[i] / box.volume + # build aabb of all particles and generate a parent # neighbor list with the RDF cutoff aabb = freud.locality.AABBQuery(box, snap.particles.position) @@ -1259,34 +1267,19 @@ def act(self, timestep): ) neighbors.filter(bond_exclusion_filter) - for i, j in self._rdf: - # make a neighbor list of all particles of type i - # whose neighbors are type j - neighbors_ij = neighbors.copy() - filter_ij = numpy.logical_and( - type_masks[i][neighbors[:, 0]], - type_masks[j][neighbors[:, 1]], - ) - # if the types are different, also consider - # neighbors of j that are type i - if j != i: - filter_ji = numpy.logical_and( - type_masks[j][neighbors[:, 0]], - type_masks[i][neighbors[:, 1]], + for i in self.types: + for j in self.types: + filter_ij = numpy.logical_and( + type_masks[i][neighbors[:, 0]], + type_masks[j][neighbors[:, 1]], ) - filter_ij = numpy.logical_or(filter_ij, filter_ji) - neighbors_ij.filter(filter_ij) - - # it doesn't look like it but this calculation should - # only calculate g_ij because we prefiltered the neighbors - # - # resetting when the samples are zero clears the RDF - self._rdf[i, j].compute( - aabb, - snap.particles.position, - neighbors=neighbors_ij, - reset=(self.num_samples == 0), - ) + counts, _ = numpy.histogram( + neighbors.distances[filter_ij], + bins=self.rdf_params["bins"], + range=(0, self.rdf_params["stop"]), + ) + self._rdf_counts[i, j] += counts + self.num_samples += 1 if compute_rdf: self.num_rdf_samples += 1 @@ -1312,15 +1305,20 @@ def reset(self): self._N[i] = 0 if self.rdf_params is not None: - self._rdf = collections.PairMatrix(self.types) - for i, j in self._rdf: - self._rdf[i, j] = freud.density.RDF( - bins=self.rdf_params["bins"], - r_max=self.rdf_params["stop"], - normalize=(i == j), - ) + self._rdf_counts = {} + self._rdf_density = {} + self._rdf_num_origins = {} + for i in self.types: + self._rdf_density[i] = 0 + self._rdf_num_origins[i] = 0 + for j in self.types: + self._rdf_counts[i, j] = numpy.zeros( + self.rdf_params["bins"], dtype=int + ) else: - self._rdf = None + self._rdf_counts = None + self._rdf_density = None + self._rdf_num_origins = None self.num_samples = 0 self.num_rdf_samples = 0 @@ -1389,16 +1387,49 @@ def rdf(self): Radial distribution functions. """ if self.num_rdf_samples > 0: + bin_edges = numpy.linspace( + 0, self.rdf_params["stop"], self.rdf_params["bins"] + 1 + ) + bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + if self.dimension == 3: + bin_extents = (4 * numpy.pi / 3) * ( + bin_edges[1:] ** 3 - bin_edges[:-1] ** 3 + ) + elif self.dimension == 2: + bin_extents = numpy.pi * (bin_edges[1:] ** 2 - bin_edges[:-1] ** 2) + rdf = collections.PairMatrix(self.types) - for pair in rdf: + for i, j in rdf: if mpi.world.rank == 0: - gr = numpy.column_stack( - (self._rdf[pair].bin_centers, self._rdf[pair].rdf) - ) + density = { + k: self._rdf_density[k] / self.num_rdf_samples + for k in self.types + } + g = numpy.zeros_like(bin_centers) + if i == j: + if self._rdf_num_origins[i] > 0 and density[i] > 0: + g = self._rdf_counts[i, i] / ( + self._rdf_num_origins[i] * density[i] * bin_extents + ) + else: + # this takes the weighted average of g_ij and g_ji + num_ij_origins = ( + self._rdf_num_origins[i] + self._rdf_num_origins[j] + ) + if num_ij_origins > 0: + if density[j] > 0: + g += self._rdf_counts[i, j] / ( + num_ij_origins * density[j] * bin_extents + ) + if density[i] > 0: + g += self._rdf_counts[j, i] / ( + num_ij_origins * density[i] * bin_extents + ) + gr = numpy.column_stack((bin_centers, g)) else: gr = None gr = mpi.world.bcast_numpy(gr, root=0) - rdf[pair] = ensemble.RDF(gr[:, 0], gr[:, 1]) + rdf[i, j] = ensemble.RDF(gr[:, 0], gr[:, 1]) return rdf else: return None From ccfff2db9a7650478fb8863052da62f46852f621 Mon Sep 17 00:00:00 2001 From: clpetix Date: Mon, 29 Apr 2024 09:37:57 -0500 Subject: [PATCH 06/10] Set number of rdf origins. --- src/relentless/simulate/hoomd.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 2b89e523..84fdca3b 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1240,6 +1240,7 @@ def act(self, timestep): # calculate average density per type for i in self.types: self._rdf_density[i] += N[i] / box.volume + self._rdf_num_origins[i] += N[i] # build aabb of all particles and generate a parent # neighbor list with the RDF cutoff From 3d326dc07639d19c97faba723eb197db7702446a Mon Sep 17 00:00:00 2001 From: clpetix Date: Thu, 6 Jun 2024 16:11:16 -0500 Subject: [PATCH 07/10] Filter bond exclusions using Cantor pairing function. --- src/relentless/simulate/hoomd.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 393cc49f..6d1233a3 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1255,11 +1255,16 @@ def act(self, timestep): bonds = numpy.vstack( [self.bonds, numpy.flip(self.bonds, axis=1)], ) - # list intersect method from: - # https://stackoverflow.com/a/67113105 - bond_exclusion_filter = ( - ~(neighbors[:, None] == bonds).all(-1).any(1) - ) + # list intersect using Cantor Pairing Function (pi): + # https://en.wikipedia.org/wiki/Pairing_function + pi_bond = (bonds[:, 0] + bonds[:, 1]) * ( + bonds[:, 0] + bonds[:, 1] + 1 + ) / 2 + bonds[:, 1] + pi_neighbor = (neighbors[:, 0] + neighbors[:, 1]) * ( + neighbors[:, 0] + neighbors[:, 1] + 1 + ) / 2 + neighbors[:, 1] + bond_exclusion_filter = ~numpy.isin(pi_neighbor, pi_bond) + neighbors.filter(bond_exclusion_filter) for i in self.types: From 1d19b9314cc4dcb0bc5e7c5e5ce9ac9b2a173cb7 Mon Sep 17 00:00:00 2001 From: clpetix Date: Fri, 4 Oct 2024 15:16:26 -0500 Subject: [PATCH 08/10] Add bond potentials and exclusions. --- src/relentless/simulate/hoomd.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index f9c2df12..bb0c9e27 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -65,6 +65,7 @@ def _call_v3(self, sim): snap = sim["engine"]["_hoomd"].state.get_snapshot() sim.masses = self._get_masses_from_snapshot(sim, snap) sim[self]["_bonds"] = self._get_bonds_from_snapshot(sim, snap) + sim[self]["_angles"] = self._get_angles_from_snapshot(sim, snap) self._assert_dimension_safe(sim, snap) # create the potentials, defer attaching until later neighbor_list = hoomd.md.nlist.Tree(buffer=sim.potentials.pair.neighbor_buffer) @@ -156,6 +157,9 @@ def _get_masses_from_snapshot(self, sim, snap): def _get_bonds_from_snapshot(self, sim, snap): return snap.bonds.group + def _get_angles_from_snapshot(self, sim, snap): + return snap.angles.group + def _assert_dimension_safe(self, sim, snap): if sim.dimension == 3: dim_safe = True @@ -992,6 +996,7 @@ def _pre_run_v3(self, sim, sim_op): rdf_params=self._get_rdf_params(sim), constraints=self._get_constrained_quantities(sim, sim_op), bonds=sim[sim.initializer]["_bonds"], + angles=sim[sim.initializer]["_angles"], ) sim[self]["_hoomd_thermo_callback"] = hoomd.write.CustomWriter( trigger=self.every, @@ -1184,6 +1189,7 @@ def __init__( rdf_params=None, constraints=None, bonds=None, + angles=None, ): if dimension not in (2, 3): raise ValueError("Only 2 or 3 dimensions supported") @@ -1195,6 +1201,7 @@ def __init__( self.rdf_params = rdf_params self.constraints = constraints if constraints is not None else {} self.bonds = bonds + self.angles = angles # this method handles all the initialization self.reset() @@ -1308,6 +1315,21 @@ def act(self, timestep): neighbors.filter(bond_exclusion_filter) + if snap.angles is not None and len(neighbors[:]) > 0: + angles = numpy.vstack( + [self.angles, numpy.flip(self.angles, axis=1)], + ) + # list intersect using Cantor Pairing Function (pi): + # https://en.wikipedia.org/wiki/Pairing_function + pi_angle = (angles[:, 0] + angles[:, 2]) * ( + angles[:, 0] + angles[:, 2] + 1 + ) / 2 + angles[:, 2] + pi_neighbor = (neighbors[:, 0] + neighbors[:, 1]) * ( + neighbors[:, 0] + neighbors[:, 1] + 1 + ) / 2 + neighbors[:, 1] + bond_exclusion_filter = ~numpy.isin(pi_neighbor, pi_angle) + + neighbors.filter(bond_exclusion_filter) for i in self.types: for j in self.types: filter_ij = numpy.logical_and( From d05147b722303216b24ef5169a77405e51f84347 Mon Sep 17 00:00:00 2001 From: clpetix Date: Fri, 4 Oct 2024 15:21:51 -0500 Subject: [PATCH 09/10] Fix error in naming. --- src/relentless/simulate/hoomd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index bb0c9e27..730ca8db 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -1327,9 +1327,9 @@ def act(self, timestep): pi_neighbor = (neighbors[:, 0] + neighbors[:, 1]) * ( neighbors[:, 0] + neighbors[:, 1] + 1 ) / 2 + neighbors[:, 1] - bond_exclusion_filter = ~numpy.isin(pi_neighbor, pi_angle) + angle_exclusion_filter = ~numpy.isin(pi_neighbor, pi_angle) - neighbors.filter(bond_exclusion_filter) + neighbors.filter(angle_exclusion_filter) for i in self.types: for j in self.types: filter_ij = numpy.logical_and( From 0a9ad405003bc09caf475b72833b555ba9298a46 Mon Sep 17 00:00:00 2001 From: clpetix Date: Wed, 9 Oct 2024 12:19:52 -0500 Subject: [PATCH 10/10] Add angle exclusions. --- src/relentless/simulate/hoomd.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 730ca8db..46a04d54 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -68,7 +68,9 @@ def _call_v3(self, sim): sim[self]["_angles"] = self._get_angles_from_snapshot(sim, snap) self._assert_dimension_safe(sim, snap) # create the potentials, defer attaching until later - neighbor_list = hoomd.md.nlist.Tree(buffer=sim.potentials.pair.neighbor_buffer) + neighbor_list = hoomd.md.nlist.Tree( + buffer=sim.potentials.pair.neighbor_buffer, exclusions=["bond", "angle"] + ) pair_potential = hoomd.md.pair.Table(nlist=neighbor_list) r, u, f = sim.potentials.pair.pairwise_energy_and_force( sim.types, tight=True, minimum_num=2