diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a5529d9..33ea182 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,7 @@ set(COMPONENT_NAME azplugins)
 # TODO: List all host C++ source code files in _${COMPONENT_NAME}_sources.
 set(_${COMPONENT_NAME}_sources
     ConstantFlow.cc
+    DynamicBondUpdater.cc
     export_ImagePotentialBondHarmonic.cc
     module.cc
     ParabolicFlow.cc
@@ -11,6 +12,7 @@ set(_${COMPONENT_NAME}_sources
     )
 if (ENABLE_HIP)
     list(APPEND _${COMPONENT_NAME}_sources
+    DynamicBondUpdaterGPU.cc
     export_ImagePotentialBondHarmonicGPU.cc
     VelocityComputeGPU.cc
     )
@@ -18,6 +20,7 @@ endif()
 
 # TODO: List all GPU C++ source code files in _${COMPONENT_NAME}_cu_sources.
 set(_${COMPONENT_NAME}_cu_sources
+    DynamicBondUpdaterGPU.cu
     ImagePotentialBondGPUKernel.cu
     VelocityComputeGPU.cu
     VelocityFieldComputeGPU.cu
@@ -31,6 +34,7 @@ set(python_files
     conftest.py
     external.py
     flow.py
+    update.py
     pair.py
     wall.py
     )
diff --git a/src/DynamicBondUpdater.cc b/src/DynamicBondUpdater.cc
new file mode 100644
index 0000000..0f463c9
--- /dev/null
+++ b/src/DynamicBondUpdater.cc
@@ -0,0 +1,872 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// This file is part of the azplugins project, released under the Modified BSD License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdater.cc
+ * \brief Definition of DynamicBondUpdater
+ */
+
+#include "DynamicBondUpdater.h"
+#include "hoomd/RandomNumbers.h"
+#include "RNGIdentifiers.h"
+#include "hoomd/md/NeighborListTree.h"
+
+#include <set>
+#include <algorithm>
+
+namespace hoomd
+{
+namespace azplugins
+{
+
+/*!
+ * \param sysdef System definition
+ *
+ * The system is initialized in a configuration that will be not forming any bonds.
+ * This constructor requires that the user properly initialize the system via setters.
+ */
+DynamicBondUpdater::DynamicBondUpdater(std::shared_ptr<SystemDefinition> sysdef,
+                                       std::shared_ptr<Trigger> trigger,
+                                       std::shared_ptr<md::NeighborList> pair_nlist,
+                                       std::shared_ptr<ParticleGroup> group_1,
+                                       std::shared_ptr<ParticleGroup> group_2,
+                                       uint16_t seed)
+          : Updater(sysdef, trigger),
+           m_group_1(group_1),
+           m_group_2(group_2),
+           m_groups_identical(false),
+           m_r_cut(0),
+           m_probability(0.0),
+           m_bond_type(0),
+           m_max_bonds_group_1(0),
+           m_max_bonds_group_2(0),
+           m_seed(seed),
+           m_pair_nlist(pair_nlist),
+           m_pair_nlist_exclusions_set(true),
+           m_box_changed(true),
+           m_max_N_changed(true)
+    {
+      m_exec_conf->msg->notice(5) << "Constructing DynamicBondUpdater" << std::endl;
+
+      m_pdata->getBoxChangeSignal().connect<DynamicBondUpdater, &DynamicBondUpdater::slotBoxChanged>(this);
+      m_pdata->getGlobalParticleNumberChangeSignal().connect<DynamicBondUpdater, &DynamicBondUpdater::slotNumParticlesChanged>(this);
+
+      m_bond_data = m_sysdef->getBondData();
+
+      m_pair_internal_nlist = std::shared_ptr<hoomd::md::NeighborList>(
+        new hoomd::md::NeighborListTree(sysdef, 0.0));
+      m_pair_internal_nlist->setStorageMode(hoomd::md::NeighborList::full);
+
+      setGroupOverlap();
+
+      setCutoffs();
+
+      m_max_bonds = 4;
+      m_max_bonds_overflow = 0;
+      m_num_all_possible_bonds = 0;
+
+
+    }
+
+
+DynamicBondUpdater::DynamicBondUpdater(std::shared_ptr<SystemDefinition> sysdef,
+                         std::shared_ptr<Trigger> trigger,
+                         std::shared_ptr<md::NeighborList> pair_nlist,
+                         std::shared_ptr<ParticleGroup> group_1,
+                         std::shared_ptr<ParticleGroup> group_2,
+                         uint16_t seed,
+                         const Scalar r_cut,
+                         const Scalar probability,
+                         unsigned int max_bonds_group_1,
+                         unsigned int max_bonds_group_2,
+                         unsigned int bond_type
+                         )
+        : Updater(sysdef, trigger),
+         m_group_1(group_1),
+         m_group_2(group_2),
+         m_groups_identical(false),
+         m_r_cut(r_cut),
+         m_probability(probability),
+         m_bond_type(bond_type),
+         m_max_bonds_group_1(max_bonds_group_1),
+         m_max_bonds_group_2(max_bonds_group_2),
+         m_seed(seed),
+         m_pair_nlist(pair_nlist),
+         m_pair_nlist_exclusions_set(true),
+         m_box_changed(true),
+         m_max_N_changed(true)
+    {
+
+    m_exec_conf->msg->notice(5) << "Constructing DynamicBondUpdater" << std::endl;
+
+    m_pdata->getBoxChangeSignal().connect<DynamicBondUpdater, &DynamicBondUpdater::slotBoxChanged>(this);
+    m_pdata->getGlobalParticleNumberChangeSignal().connect<DynamicBondUpdater, &DynamicBondUpdater::slotNumParticlesChanged>(this);
+
+    m_bond_data = m_sysdef->getBondData();
+
+    m_pair_internal_nlist = std::shared_ptr<hoomd::md::NeighborList>(
+        new hoomd::md::NeighborListTree(sysdef, 0.0));
+
+    setGroupOverlap();
+    m_pair_internal_nlist->setStorageMode(hoomd::md::NeighborList::full);
+    setCutoffs();
+
+    m_max_bonds = 4;
+    m_max_bonds_overflow = 0;
+    m_num_all_possible_bonds = 0;
+
+
+    }
+
+DynamicBondUpdater::~DynamicBondUpdater()
+    {
+    m_exec_conf->msg->notice(5) << "Destroying DynamicBondUpdater" << std::endl;
+
+    m_pdata->getBoxChangeSignal().disconnect<DynamicBondUpdater, &DynamicBondUpdater::slotBoxChanged>(this);
+    m_pdata->getGlobalParticleNumberChangeSignal().disconnect<DynamicBondUpdater, &DynamicBondUpdater::slotNumParticlesChanged>(this);
+
+    }
+
+/*!
+* \param timestep Timestep update is called
+*/
+void DynamicBondUpdater::update(uint64_t timestep)
+    {
+      //Scalar test = m_pair_internal_nlist->getRCut(pybind11::make_tuple('A','A'));
+
+      // don't do anything if either one of the groups is  empty
+      if (m_group_1->getNumMembers() == 0 || m_group_2->getNumMembers() == 0)
+      return;
+      // don't do anything if maximum number of bonds is zero
+      if (m_max_bonds_group_1 == 0 || m_max_bonds_group_2 == 0)
+      return;
+
+      // update properties that depend on the box
+      if (m_box_changed)
+      {
+        checkBoxSize();
+        m_box_changed = false;
+      }
+
+      // update properties that depend on the number of particles
+      if (m_max_N_changed)
+      {
+        allocateParticleArrays();
+        m_max_N_changed = false;
+      }
+
+      // rebuild the list of possible bonds until there is no overflow
+      bool overflowed = false;
+
+      do
+      {
+        m_pair_internal_nlist->compute(timestep);
+
+        ArrayHandle<unsigned int> h_n_neigh(m_pair_internal_nlist->getNNeighArray(), access_location::host, access_mode::read);
+        for (unsigned int group_idx = 0; group_idx < m_group_1->getNumMembers(); group_idx++)
+          {
+          unsigned int i = m_group_1->getMemberIndex(group_idx);
+          const unsigned int n_neigh = h_n_neigh.data[i];
+          if (n_neigh > m_max_bonds)
+            m_max_bonds = n_neigh;
+          }
+
+        overflowed = m_max_bonds < m_max_bonds_overflow;
+        // if we overflowed, need to reallocate memory and re-traverse the neighbor list
+        if (overflowed)
+        {
+          resizePossibleBondlists();
+        }
+      } while (overflowed);
+
+      filterPossibleBonds();
+      // this function is not easily implemented on the GPU, uses addBondedGroup()
+      makeBonds(timestep);
+
+    }
+
+// todo: should go into helper class/separate file?
+// bonds need to be sorted such that dublicates end up next to each other, otherwise
+// unique will not work properly. If the  bond length of the potential bond is different, we can
+// sort according to that, but there might be the case where multiple possible bond lengths are exactly identical,
+// e.g. particles on a lattice.
+// This is hiracical sorting: first according to possible bond distance r_ab_sq, then after first tag_a, last after second tag_b.
+// Should work given that the tags are oredered within each pair, (tag_a,tag_b,d_ab_sq) with tag_a < tag_b.
+
+// todo: because it's possible uniquely map a pair to a single int (and back!) with a pairing function,
+// the possible bond array could be restructured into a different data structure?
+// if we don't keep the possible bond length, a unsigned int array could hold all information needed
+// when would that lead to problems with overflow from 0.5*(tag_a+tag_b)*(tag_a+tag_b+1)+tag_b being too large?
+
+bool SortBonds(Scalar3 i, Scalar3 j)
+    {
+      const Scalar r_sq_1 = i.z;
+      const Scalar r_sq_2 = j.z;
+      if (r_sq_1==r_sq_2)
+      {
+        const unsigned int tag_11 = __scalar_as_int(i.x);
+        const unsigned int tag_21 = __scalar_as_int(j.x);
+        if (tag_11==tag_21)
+        {
+          const unsigned int tag_12 = __scalar_as_int(i.y);
+          const unsigned int tag_22 = __scalar_as_int(j.y);
+          return tag_22>tag_12;
+        }
+        else
+        {
+          return tag_21>tag_11;
+        }
+      }
+      else
+      {
+        return r_sq_2>r_sq_1;
+      }
+    }
+
+// todo: migrate to separate file/class?
+// Cantor paring function can also be used for comparison
+bool CompareBonds(Scalar3 i, Scalar3 j)
+    {
+      const unsigned int tag_11 = __scalar_as_int(i.x);
+      const unsigned int tag_12 = __scalar_as_int(i.y);
+      const unsigned int tag_21 = __scalar_as_int(j.x);
+      const unsigned int tag_22 = __scalar_as_int(j.y);
+
+      if ((tag_11==tag_21 && tag_12==tag_22))
+      {
+        return true;
+      }
+      else
+      {
+        return false;
+      }
+    }
+
+
+bool DynamicBondUpdater::CheckisExistingLegalBond(Scalar3 i)
+    {
+      const unsigned int tag_1 = __scalar_as_int(i.x);
+      const unsigned int tag_2 = __scalar_as_int(i.y);
+
+      // (0,0,0.0) is the default "empty" value - todo: have a "invalid" unsigned int ? how to fill/reset with memset()?
+      if (tag_1==0 && tag_2==0 )
+      {
+        return true;
+      }
+      else
+      {
+        return isExistingBond(tag_1,tag_2);
+      }
+    }
+
+void DynamicBondUpdater::calculateExistingBonds()
+    {
+      {
+      // reset exisitng bond list
+      ArrayHandle<unsigned int> h_rtag(m_pdata->getRTags(), access_location::host, access_mode::read);
+      ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::overwrite);
+      memset((void*)h_n_existing_bonds.data,0,sizeof(unsigned int)*m_pdata->getMaxN());
+
+      ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::overwrite);
+      memset((void*)h_existing_bonds_list.data,0,sizeof(unsigned int)*m_group_1->getNumMembers()*m_existing_bonds_list_indexer.getH());
+      }
+      {
+      ArrayHandle<typename BondData::members_t> h_bonds(m_bond_data->getMembersArray(), access_location::host, access_mode::read);
+
+      // for each of the bonds in the system - regardless of their type
+      const unsigned int size = (unsigned int)m_bond_data->getN();
+      for (unsigned int i = 0; i < size; i++)
+      {
+        // lookup the tag of each of the particles participating in the bond
+        const typename BondData::members_t& bond = h_bonds.data[i];
+        unsigned int tag1 = bond.tag[0];
+        unsigned int tag2 = bond.tag[1];
+
+       // @mphoward: The next section is the one that needed to be dublicated to accomodate
+       // the changes in hoomd (GPUArray and GPUVector)
+
+        //AddtoExistingBonds(tag1,tag2);
+
+        // BEGIN DynamicBondUpdater::AddtoExistingBonds
+        assert(tag1 <= m_pdata->getMaximumTag());
+        assert(tag2 <= m_pdata->getMaximumTag());
+
+        bool overflowed = false;
+
+        ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::readwrite);
+
+
+        // resize the list if necessary
+        if (h_n_existing_bonds.data[tag1] == m_existing_bonds_list_indexer.getH())
+        overflowed = true;
+        if (h_n_existing_bonds.data[tag2] == m_existing_bonds_list_indexer.getH())
+        overflowed = true;
+
+        if (overflowed) resizeExistingBondList();
+
+        ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::readwrite);
+
+        // add tag_b to tag_a's existing bonds list
+        unsigned int pos_a = h_n_existing_bonds.data[tag1];
+        assert(pos_a < m_existing_bonds_list_indexer.getH());
+        h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag1,pos_a)] = tag2;
+        h_n_existing_bonds.data[tag1]++;
+
+        // add tag_a to tag_b's existing bonds list
+        unsigned int pos_b = h_n_existing_bonds.data[tag2];
+        assert(pos_b < m_existing_bonds_list_indexer.getH());
+        h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag2,pos_b)] = tag1;
+        h_n_existing_bonds.data[tag2]++;
+
+        // END DynamicBondUpdater::AddtoExistingBonds
+
+      }
+      }
+
+    }
+
+/*! \param tag1 First particle tag in the pair
+    \param tag2 Second particle tag in the pair
+    \return true if the particles \a tag1 and \a tag2 are bonded
+*/
+bool DynamicBondUpdater::isExistingBond(unsigned int tag1, unsigned int tag2)
+    {
+      {
+      ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::read);
+      ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::read);
+      unsigned int n_existing_bonds = h_n_existing_bonds.data[tag1];
+
+      for (unsigned int i = 0; i < n_existing_bonds; i++)
+      {
+        if (h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag1,i)] == tag2)
+        return true;
+      }
+      return false;
+      }
+    }
+
+
+// @mphoward: This is the function that I originally had that doesn't work anymore due to
+// the changes in hoomd (GPUArray and GPUVector).
+
+/*! \param tag1 First particle tag in the pair
+    \param tag2 Second particle tag in the pair
+    adds a bond between the tag1 and tag2 to the existing bonds list
+*/
+void DynamicBondUpdater::AddtoExistingBonds(unsigned int tag_a,unsigned int tag_b)
+    {
+
+      // BEGIN DynamicBondUpdater::AddtoExistingBonds
+      assert(tag_a <= m_pdata->getMaximumTag());
+      assert(tag_b <= m_pdata->getMaximumTag());
+
+      bool overflowed = false;
+
+      ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::readwrite);
+
+
+      // resize the list if necessary
+      if (h_n_existing_bonds.data[tag_a] == m_existing_bonds_list_indexer.getH())
+      overflowed = true;
+      if (h_n_existing_bonds.data[tag_b] == m_existing_bonds_list_indexer.getH())
+      overflowed = true;
+
+      if (overflowed) resizeExistingBondList();
+
+      ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::readwrite);
+
+      // add tag_b to tag_a's existing bonds list
+      unsigned int pos_a = h_n_existing_bonds.data[tag_a];
+      assert(pos_a < m_existing_bonds_list_indexer.getH());
+      h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag_a,pos_a)] = tag_b;
+      h_n_existing_bonds.data[tag_a]++;
+
+      // add tag_a to tag_b's existing bonds list
+      unsigned int pos_b = h_n_existing_bonds.data[tag_b];
+      assert(pos_b < m_existing_bonds_list_indexer.getH());
+      h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag_b,pos_b)] = tag_a;
+      h_n_existing_bonds.data[tag_b]++;
+
+      // END DynamicBondUpdater::AddtoExistingBonds
+
+    }
+
+// grows the existing bonds list and its indexer when needed
+void DynamicBondUpdater::resizeExistingBondList()
+    {
+      {
+      unsigned int new_height = m_existing_bonds_list_indexer.getH() + 1;
+      m_existing_bonds_list.resize(m_pdata->getRTags().size(), new_height);
+      // update the indexer
+      m_existing_bonds_list_indexer = Index2D((unsigned int)m_existing_bonds_list.getPitch(), new_height);
+      m_exec_conf->msg->notice(6) << "DynamicBondUpdater: (Re-)size existing bond list, new size " << new_height << " bonds per particle " << std::endl;
+      }
+    }
+
+// grows the all possible bonds list when needed in increments of 4, inspired by the neighbor list
+void DynamicBondUpdater::resizePossibleBondlists()
+    {
+      {
+      // round up to nearest multiple of 4
+      m_max_bonds_overflow = (m_max_bonds_overflow > 4) ? (m_max_bonds_overflow + 3) & ~3 : 4;
+      m_max_bonds = m_max_bonds_overflow;
+      m_max_bonds_overflow = 0;
+      unsigned int size = m_group_1->getNumMembers()*m_max_bonds;
+      m_all_possible_bonds.resize(size);
+      m_num_all_possible_bonds=0;
+
+      m_exec_conf->msg->notice(6) << "DynamicBondUpdater: (Re-)size possible bond list, new size " << m_max_bonds << " bonds per particle " << std::endl;
+      }
+    }
+
+
+
+// allocates all arrays depending on the particles and groups
+void DynamicBondUpdater::allocateParticleArrays()
+    {
+      {
+      GPUArray<Scalar3> all_possible_bonds(m_group_1->getNumMembers() *m_max_bonds, m_exec_conf);
+      m_all_possible_bonds.swap(all_possible_bonds);
+
+      GPUArray<unsigned int> n_existing_bonds(m_pdata->getRTags().size(), m_exec_conf);
+      m_n_existing_bonds.swap(n_existing_bonds);
+
+      GPUArray<unsigned int> existing_bonds_list(m_pdata->getRTags().size(),1, m_exec_conf);
+      m_existing_bonds_list.swap(existing_bonds_list);
+      m_existing_bonds_list_indexer = Index2D((unsigned int)m_existing_bonds_list.getPitch(), (unsigned int)m_existing_bonds_list.getHeight());
+
+      ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::overwrite);
+      ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::overwrite);
+
+      memset(h_n_existing_bonds.data, 0, sizeof(unsigned int)*m_n_existing_bonds.getNumElements());
+      memset(h_existing_bonds_list.data, 0, sizeof(unsigned int)*m_existing_bonds_list.getNumElements());
+      }
+      calculateExistingBonds();
+    }
+
+
+
+/*! This function takes the information about neighbors between group_2 and group_1 saved in m_nlist and
+* m_n_neigh and copies pairs within the m_r_cut cutoff distance into the m_all_possible_bonds array.
+* Then, all invalid (0,0,0), dublicated, and existing bonds  are removed from m_all_possible_bonds. It
+* is sorted by distance  (shortest to longest) between the two particles in the possible bond.
+*/
+void DynamicBondUpdater::filterPossibleBonds()
+    {
+      //copy data from h_n_list to h_all_possible_bonds
+      {
+
+      ArrayHandle<unsigned int> h_nlist(m_pair_internal_nlist->getNListArray(), access_location::host, access_mode::read);
+      ArrayHandle<unsigned int> h_n_neigh(m_pair_internal_nlist->getNNeighArray(), access_location::host, access_mode::read);
+      ArrayHandle<size_t> h_n_head_list(m_pair_internal_nlist->getHeadList(), access_location::host, access_mode::read);
+
+
+      ArrayHandle<Scalar4> h_postype(m_pdata->getPositions(), access_location::host, access_mode::read);
+      ArrayHandle<unsigned int> h_tag(m_pdata->getTags(), access_location::host, access_mode::read);
+
+      ArrayHandle<Scalar3> h_all_possible_bonds(m_all_possible_bonds, access_location::host, access_mode::readwrite);
+      unsigned int size = m_group_1->getNumMembers()*m_max_bonds;
+      memset((void*) h_all_possible_bonds.data, 0, sizeof(Scalar3)*size);
+
+      const BoxDim& box = m_pdata->getBox();
+
+      // Loop over all particles in group 1
+      for (unsigned int group_idx = 0; group_idx < m_group_1->getNumMembers(); group_idx++)
+      {
+
+        unsigned int i = m_group_1->getMemberIndex(group_idx);
+        const unsigned int tag_i = h_tag.data[i];
+        const Scalar4 postype_i = h_postype.data[i];
+
+
+        unsigned int n_curr_bond = 0;
+        const Scalar r_cutsq = m_r_cut*m_r_cut;
+
+        const unsigned int n_neigh = h_n_neigh.data[i];
+        const size_t head = h_n_head_list.data[i];
+
+        // loop over all neighbors of this particle
+        for (unsigned int l=0; l<n_neigh; ++l)
+        {
+
+          // get index of neighbor from neigh_list
+          const unsigned int j = h_nlist.data[head + l];
+
+          if (m_group_2->isMember(j))
+          {
+            Scalar4 postype_j = h_postype.data[j];
+            const unsigned int tag_j = h_tag.data[j];
+
+            Scalar3 drij = make_scalar3(postype_j.x,postype_j.y,postype_j.z)
+            - make_scalar3(postype_i.x,postype_i.y,postype_i.z);
+
+            // apply periodic boundary conditions
+            drij = box.minImage(drij);
+            Scalar dr_sq = dot(drij,drij);
+
+            if (dr_sq < r_cutsq)
+            {
+              if (n_curr_bond < m_max_bonds)
+              {
+                Scalar3 d;
+                if(m_groups_identical)
+                {
+                  // sort the two tags in this possible bond pair if groups identical
+                  const unsigned int tag_a = tag_j > tag_i ? tag_i : tag_j;
+                  const unsigned int tag_b = tag_j > tag_i ? tag_j : tag_i;
+                  d = make_scalar3(__int_as_scalar(tag_a),__int_as_scalar(tag_b),dr_sq);
+                }
+                else
+                {
+                  d = make_scalar3(__int_as_scalar(tag_i),__int_as_scalar(tag_j),dr_sq);
+                }
+                h_all_possible_bonds.data[group_idx*m_max_bonds + n_curr_bond] = d;
+              }
+              ++n_curr_bond;
+            }
+          }
+        }
+      }
+
+
+      //now sort and select down
+      m_num_all_possible_bonds = 0;
+
+      // remove a possible bond if it already exists. It also removes zeros, e.g.
+      // (0,0,0), which fill the unused spots in the array.
+      auto last2 = std::remove_if(h_all_possible_bonds.data,
+                   h_all_possible_bonds.data + size,
+                   [this](Scalar3 i) {return CheckisExistingLegalBond(i); });
+
+      m_num_all_possible_bonds = (unsigned int) std::distance(h_all_possible_bonds.data,last2);
+
+      // then sort array by distance between particles in the found possible bond pairs
+      // performance is better if remove_if happens before sort
+      std::sort(h_all_possible_bonds.data, h_all_possible_bonds.data + m_num_all_possible_bonds, SortBonds);
+
+      // now make sure each possible bond is in the array only once by comparing tags
+      auto last = std::unique(h_all_possible_bonds.data, h_all_possible_bonds.data + m_num_all_possible_bonds, CompareBonds);
+      m_num_all_possible_bonds = (unsigned int)std::distance(h_all_possible_bonds.data,last);
+
+
+      // at this point, the sub-array: h_all_possible_bonds[0,m_num_all_possible_bonds]
+      // should contain only unique entries of possible bonds which are not yet formed.
+      }
+      }
+
+/*! This function actually creates the bonds by looping over the entries in m_all_possible_bonds
+*  and adding them to the system (m_bond_data->addBondedGroup), to the existing bonds, as well as
+*  to the neighbor list used by the rest of the simulation if the exclusions of that neighbor list
+*  should be updated.
+*
+* Note: this function is very hard to parallelize on the GPU since we need to go through the bonds sequentially
+* to prevent forming too many bonds in one step. Have not found a good way of doing this on the GPU.
+*/
+void DynamicBondUpdater::makeBonds(uint64_t timestep)
+  {
+
+    ArrayHandle<Scalar3> h_all_possible_bonds(m_all_possible_bonds, access_location::host, access_mode::read);
+
+    ArrayHandle<unsigned int> h_n_existing_bonds(m_n_existing_bonds, access_location::host, access_mode::readwrite);
+
+    // we need to count how many bonds are in the h_all_possible_bonds array for a given tag
+    // so that we don't end up forming too many bonds in one step. "AddtoExistingBonds" increases the count in
+    // h_n_existing_bonds in the for loop below as we go, so no extra bookkeeping should be needed.
+    // This also makes it very difficult to do on the GPU.
+
+    ArrayHandle<unsigned int> h_rtag(m_pdata->getRTags(), access_location::host, access_mode::read);
+
+    //todo: can this for loop be simplified/parallelized?
+    for (unsigned int i = 0; i < m_num_all_possible_bonds; i++)
+    {
+      Scalar3 d = h_all_possible_bonds.data[i];
+
+      unsigned int tag_i = __scalar_as_int(d.x);
+      unsigned int tag_j = __scalar_as_int(d.y);
+
+     //todo: put in other external criteria here, e.g. max number of bonds possible in one step, etc.
+     //todo: randomize which bonds are formed or keep them ordered by their distances ?
+     //todo: would it be faster/better to create the rng outside of the loop?
+    hoomd::RandomGenerator rng(
+                hoomd::Seed(hoomd::azplugins::detail::RNGIdentifier::DynamicBondUpdater,
+                            timestep,
+                            m_seed),
+                hoomd::Counter());
+
+     hoomd::UniformDistribution<Scalar> uniform(0, 1);
+     const Scalar random = uniform(rng);
+
+      if ((m_max_bonds_group_1 > h_n_existing_bonds.data[tag_i]) &&
+          (m_max_bonds_group_2 > h_n_existing_bonds.data[tag_j]) &&
+          (random < m_probability))
+        {
+          m_bond_data->addBondedGroup(Bond(m_bond_type,tag_i,tag_j));
+          //AddtoExistingBonds(tag_i,tag_j);
+
+         // @mphoward: The next section is the one that needed to be dublicated to accomodate
+         // the changes in hoomd (GPUArray and GPUVector)
+
+          // BEGIN DynamicBondUpdater::AddtoExistingBonds
+          assert(tag_i <= m_pdata->getMaximumTag());
+          assert(tag_j <= m_pdata->getMaximumTag());
+
+          bool overflowed = false;
+
+          // resize the list if necessary
+          if (h_n_existing_bonds.data[tag_i] == m_existing_bonds_list_indexer.getH())
+          overflowed = true;
+          if (h_n_existing_bonds.data[tag_j] == m_existing_bonds_list_indexer.getH())
+          overflowed = true;
+
+          if (overflowed) resizeExistingBondList();
+
+          ArrayHandle<unsigned int> h_existing_bonds_list(m_existing_bonds_list, access_location::host, access_mode::readwrite);
+
+          // add tag_b to tag_a's existing bonds list
+          unsigned int pos_a = h_n_existing_bonds.data[tag_i];
+          assert(pos_a < m_existing_bonds_list_indexer.getH());
+          h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag_i,pos_a)] = tag_j;
+          h_n_existing_bonds.data[tag_i]++;
+
+          // add tag_a to tag_b's existing bonds list
+          unsigned int pos_b = h_n_existing_bonds.data[tag_j];
+          assert(pos_b < m_existing_bonds_list_indexer.getH());
+          h_existing_bonds_list.data[m_existing_bonds_list_indexer(tag_j,pos_b)] = tag_i;
+          h_n_existing_bonds.data[tag_j]++;
+
+          // END DynamicBondUpdater::AddtoExistingBonds
+
+
+          if (m_pair_nlist_exclusions_set)
+            {
+            m_pair_nlist -> addExclusion(tag_i,tag_j);
+            m_pair_internal_nlist -> addExclusion(tag_i,tag_j);
+            }
+        }
+      }
+
+    }
+
+
+
+/*!
+* Check that the largest neighbor search radius is not bigger than twice the shortest box size.
+* Raises an error if this condition is not met.
+*/
+void DynamicBondUpdater::checkBoxSize()
+    {
+      const BoxDim& box = m_pdata->getBox();
+      const uchar3 periodic = box.getPeriodic();
+
+      // check that rcut fits in the box
+      Scalar3 nearest_plane_distance = box.getNearestPlaneDistance();
+      Scalar rmax = m_r_cut;
+
+      if ((periodic.x && nearest_plane_distance.x <= rmax * 2.0) ||
+      (periodic.y && nearest_plane_distance.y <= rmax * 2.0) ||
+      (m_sysdef->getNDimensions() == 3 && periodic.z && nearest_plane_distance.z <= rmax * 2.0))
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: Simulation box is too small! Particles would be interacting with themselves." << std::endl;
+        throw std::runtime_error("Error in DynamicBondUpdater, Simulation box too small.");
+      }
+    }
+
+/*! Calculates if the two groups have overlap or not. Returns an error if  partial
+* overlap is detected.
+*/
+void DynamicBondUpdater::setGroupOverlap()
+    {
+      if(m_group_1->getNumMembers()==0)
+      {
+        m_exec_conf->msg->warning() << "DynamicBondUpdater: First group group_1 appears to be empty. No bonds will be formed. " << std::endl;
+      }
+
+      if(m_group_2->getNumMembers()==0)
+      {
+        m_exec_conf->msg->warning() << "DynamicBondUpdater: Second group group_2 appears to be empty. No bonds will be formed. " << std::endl;
+      }
+      {
+      //check if the two groups are either identical or have no overlap
+      ArrayHandle<unsigned int> h_index_group_1(m_group_1->getIndexArray(), access_location::host, access_mode::read);
+
+      // count particles which are in both groups. Should be either zero of them or all of them.
+      unsigned int overlap = 0;
+      for (unsigned int i=0; i<m_group_1->getNumMembers(); ++i)
+      {
+        unsigned int idx = h_index_group_1.data[i];
+        if (m_group_2->isMember(idx))
+        overlap++;
+      }
+
+      if( overlap>0 && overlap != m_group_1->getNumMembers())
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: group 1 and group 2 have " << overlap  << " overlaps. Partially overlapping groups are not implemented." << std::endl;
+        throw std::runtime_error("Partial overlapping groups in DynamicBondUpdater");
+      }
+
+      if(overlap==m_group_1->getNumMembers())
+      {
+        m_groups_identical=true;
+      }
+      }
+    }
+
+/*! Sets cutoffs based on types present in the two groups to save some performance from
+* the neighbor list.
+*/
+void DynamicBondUpdater::setCutoffs()
+    {
+    {
+
+    // set all rcuts to zero first, then only set the one between group1 and group 2 particle types to be m_r_cut
+     unsigned int NTypes = m_pdata -> getNTypes();
+      for (unsigned int i=0; i< NTypes; ++i)
+      {
+        for (unsigned int j=0; j< NTypes; ++j)
+          {
+            m_pair_internal_nlist->setRcut(i,j,0);
+          }
+      }
+
+    ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(),
+                               access_location::host,
+                               access_mode::read);
+
+
+    ArrayHandle<unsigned int> h_index_group_1(m_group_1->getIndexArray(), access_location::host, access_mode::read);
+
+    // finding all types in group 1
+    std::set<unsigned int> types_group_1;
+    for (unsigned int i=0; i<m_group_1->getNumMembers(); ++i)
+      {
+        unsigned int idx = h_index_group_1.data[i];
+        Scalar4 type = h_pos.data[idx];
+        types_group_1.insert(__scalar_as_int(type.w));
+      }
+
+    std::set<unsigned int> types_group_2;
+    if (m_groups_identical == false)
+    {
+    ArrayHandle<unsigned int> h_index_group_2(m_group_2->getIndexArray(), access_location::host, access_mode::read);
+
+    // finding all types in group 2
+    for (unsigned int i=0; i<m_group_2->getNumMembers(); ++i)
+      {
+        unsigned int idx = h_index_group_2.data[i];
+        Scalar4 type =  h_pos.data[idx];
+        types_group_2.insert(__scalar_as_int(type.w));
+      }
+    }
+    else
+    {
+    types_group_2 = types_group_1;
+    }
+
+      // looping over types in group 1 and 2 to set the cutoff to m_r_cut
+      for (const auto& element_1 : types_group_1)
+      {
+        for (const auto& element_2 : types_group_2)
+        {
+          m_pair_internal_nlist->setRcut(element_1,element_2,m_r_cut);
+          m_pair_internal_nlist->setRcut(element_2,element_1,m_r_cut);
+        }
+      }
+
+    }
+    }
+
+
+// Check that the given cutoff value is valid
+void DynamicBondUpdater::checkRcut()
+    {
+      if (m_r_cut <= 0.0)
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: Requested cutoff distance is less than or equal to zero" << std::endl;
+        throw std::runtime_error("Error initializing DynamicBondUpdater");
+      }
+      checkBoxSize();
+    }
+
+void DynamicBondUpdater::checkProbability()
+{
+      if (m_probability < 0.0)
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: Requested probability  is less than zero" << std::endl;
+        throw std::runtime_error("Error initializing DynamicBondUpdater");
+      }
+      else if (m_probability > 1.0)
+      {
+         m_exec_conf->msg->error() << "DynamicBondUpdater: Requested probability is larger than one" << std::endl;
+        throw std::runtime_error("Error initializing DynamicBondUpdater");
+      }
+
+}
+
+void DynamicBondUpdater::checkMaxBondsGroup()
+{
+      if (m_max_bonds_group_1 < 0.0 or m_max_bonds_group_2 < 0.0 )
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: Max number of bonds that groups can form is negative. Check parameters." << std::endl;
+        throw std::runtime_error("Error initializing DynamicBondUpdater");
+      }
+      else if (m_max_bonds_group_1 > 10 or m_max_bonds_group_2 > 10  )
+      {
+         m_exec_conf->msg->warning() << "DynamicBondUpdater: Requested number of bonds that can form is very large. This can lead to performance issues." << std::endl;
+      }
+
+}
+
+// Check that the given bond type is valid
+void DynamicBondUpdater::checkBondType()
+    {
+      if (m_bond_type >= m_bond_data -> getNTypes())
+      {
+        m_exec_conf->msg->error() << "DynamicBondUpdater: bond type id " << m_bond_type  << " is not a valid bond type." << std::endl;
+        throw std::runtime_error("Invalid bond type for DynamicBondUpdater");
+      }
+    }
+
+
+namespace detail
+{
+/*!
+* \param m Python module to export to
+*/
+
+void export_DynamicBondUpdater(pybind11::module& m)
+    {
+
+    pybind11::class_< DynamicBondUpdater, Updater,std::shared_ptr<DynamicBondUpdater>>(
+      m,
+      "DynamicBondUpdater")
+    .def(pybind11::init<std::shared_ptr<SystemDefinition>,
+      std::shared_ptr<Trigger>,
+      std::shared_ptr<md::NeighborList>,
+      std::shared_ptr<ParticleGroup>,
+      std::shared_ptr<ParticleGroup>,
+      uint16_t>())
+    .def(pybind11::init<std::shared_ptr<SystemDefinition>,
+      std::shared_ptr<Trigger>,
+      std::shared_ptr<md::NeighborList>,
+      std::shared_ptr<ParticleGroup>,
+      std::shared_ptr<ParticleGroup>,
+      uint16_t,
+      Scalar,
+      Scalar,
+      unsigned int,
+      unsigned int,
+      unsigned int>())
+      .def_property("r_cut", &DynamicBondUpdater::getRcut, &DynamicBondUpdater::setRcut)
+      .def_property("probability", &DynamicBondUpdater::getProbability, &DynamicBondUpdater::setProbability)
+      .def_property("bond_type",&DynamicBondUpdater::getBondType, &DynamicBondUpdater::setBondType)
+      .def_property("max_bonds_group_1", &DynamicBondUpdater::getMaxBondsGroup1, &DynamicBondUpdater::setMaxBondsGroup1)
+      .def_property("max_bonds_group_2", &DynamicBondUpdater::getMaxBondsGroup2, &DynamicBondUpdater::setMaxBondsGroup2);
+    }
+} // end namespace detail
+
+} // end namespace azplugins
+
+} // end namespace hoomd
\ No newline at end of file
diff --git a/src/DynamicBondUpdater.h b/src/DynamicBondUpdater.h
new file mode 100644
index 0000000..caa376d
--- /dev/null
+++ b/src/DynamicBondUpdater.h
@@ -0,0 +1,239 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdater.h
+ * \brief Declaration of DynamicBondUpdater
+ */
+
+#ifndef AZPLUGINS_DYNAMIC_BOND_UPDATER_H_
+#define AZPLUGINS_DYNAMIC_BOND_UPDATER_H_
+
+#ifdef NVCC
+#error This header cannot be compiled by nvcc
+#endif
+
+#ifndef __HIPCC__
+#include <pybind11/pybind11.h>
+#endif
+
+#include "hoomd/AABBTree.h"
+#include "hoomd/ParticleGroup.h"
+#include "hoomd/Updater.h"
+#include "hoomd/md/NeighborList.h"
+
+#ifndef __HIPCC__
+#define HOSTDEVICE __host__ __device__
+#else
+#define HOSTDEVICE
+#endif // __HIPCC__
+
+#ifndef PYBIND11_EXPORT
+#define PYBIND11_EXPORT __attribute__((visibility("default")))
+#endif
+
+namespace hoomd
+    {
+
+namespace azplugins
+    {
+
+class PYBIND11_EXPORT DynamicBondUpdater : public Updater
+    {
+    public:
+    //! Simple constructor
+    DynamicBondUpdater(std::shared_ptr<SystemDefinition> sysdef,
+                       std::shared_ptr<Trigger> trigger,
+                       std::shared_ptr<md::NeighborList> pair_nlist,
+                       std::shared_ptr<ParticleGroup> group_1,
+                       std::shared_ptr<ParticleGroup> group_2,
+                       uint16_t seed);
+
+    //! Constructor with parameters
+    DynamicBondUpdater(std::shared_ptr<SystemDefinition> sysdef,
+                       std::shared_ptr<Trigger> trigger,
+                       std::shared_ptr<md::NeighborList> pair_nlist,
+                       std::shared_ptr<ParticleGroup> group_1,
+                       std::shared_ptr<ParticleGroup> group_2,
+                       uint16_t seed,
+                       const Scalar r_cut,
+                       const Scalar probability,
+                       unsigned int max_bonds_group_1,
+                       unsigned int max_bonds_group_2,
+                       unsigned int bond_type);
+
+    //! Destructor
+    virtual ~DynamicBondUpdater();
+
+    //! update
+    virtual void update(uint64_t timestep);
+
+    //! Set the cutoff distance for finding bonds
+    /*!
+     * \param r_cut cutoff distance between particles for finding potential bonds
+     */
+    void setRcut(Scalar r_cut)
+        {
+        m_r_cut = r_cut;
+        checkRcut();
+        setCutoffs();
+        }
+    //! Get the cutoff distance between particles for finding bonds
+    Scalar getRcut()
+        {
+        return m_r_cut;
+        }
+    //! Set the bond type of the dynamically formed bonds
+    /*!
+     * \param bond_type type of bonds to be formed
+     */
+    void setBondType(unsigned int bond_type)
+        {
+        // m_bond_data = m_sysdef->getBondData();
+        // m_bond_type = m_bond_data->getTypeByName(bond_type);
+        m_bond_type = bond_type;
+        checkBondType();
+        }
+    //! Get the bond type of the dynamically formed bonds
+    unsigned int getBondType()
+        {
+        return m_bond_type;
+        }
+
+    //! Set the maximum number of bonds on particles in group_1
+    /*!
+     * \param max_bonds_group_1 max number of bonds formed by particles in group 1
+     */
+    void setMaxBondsGroup1(unsigned int max_bonds_group_1)
+        {
+        m_max_bonds_group_1 = max_bonds_group_1;
+        checkMaxBondsGroup();
+        }
+    //! Get the maximum number of bonds on particles in group_1
+    unsigned int getMaxBondsGroup1()
+        {
+        return m_max_bonds_group_1;
+        }
+    //! Set the maximum number of bonds on particles in group_2
+    /*!
+     * \param max_bonds_group_2 max number of bonds formed by particles in group_2
+     */
+    void setMaxBondsGroup2(unsigned int max_bonds_group_2)
+        {
+        m_max_bonds_group_2 = max_bonds_group_2;
+        checkMaxBondsGroup();
+        }
+    //! Get the maximum number of bonds on particles in group_2
+    unsigned int getMaxBondsGroup2()
+        {
+        return m_max_bonds_group_2;
+        }
+    //! set probablilty
+    void setProbability(Scalar probability)
+        {
+        m_probability = probability;
+        checkProbability();
+        }
+    //! Get the probability
+    Scalar getProbability()
+        {
+        return m_probability;
+        }
+
+    protected:
+    std::shared_ptr<BondData> m_bond_data; //!< Bond data
+
+    std::shared_ptr<ParticleGroup> m_group_1; //!< First particle group to form bonds with
+    std::shared_ptr<ParticleGroup> m_group_2; //!< Second particle group to form bonds with
+    bool m_groups_identical; //!< whether or not the two groups are identical. Set true if they are
+                             //!< identical.
+
+    Scalar m_r_cut; //!< cutoff for the bond forming criterion
+    Scalar
+        m_probability; //!< probability of bond formation if bond can be formed (i.e. within cutoff)
+    unsigned int m_bond_type; //!< Type id of the bond to form
+    unsigned int
+        m_max_bonds_group_1; //!< maximum number of bonds which can be formed by the first group
+    unsigned int
+        m_max_bonds_group_2; //!< maximum number of bonds which can be formed by the second group
+    uint16_t m_seed;         //!< seed for random number generator for bond probability
+
+    unsigned int m_max_bonds; //!<  maximum number of possible bonds (or neighbors) which can be
+                              //!<  found, is resized if overflow is triggered
+    unsigned int m_max_bonds_overflow; //!< registers if there is an overflow in  maximum number of
+                                       //!< possible bonds
+
+    GPUArray<Scalar3>
+        m_all_possible_bonds; //!< list of possible bonds, size: NumMembers(group_1)*m_max_bonds
+    unsigned int m_num_all_possible_bonds; //!< number of valid possible bonds at the beginning of
+                                           //!< m_all_possible_bonds
+
+    GPUArray<unsigned int>
+        m_existing_bonds_list; //!< List of existing bonded particles referenced by tag
+    GPUArray<unsigned int>
+        m_n_existing_bonds;                 //!< Number of existing bonds for a given particle tag
+    unsigned int m_max_existing_bonds_list; //!< maximum number of  bonds in list of existing bonded
+                                            //!< particles
+    Index2D
+        m_existing_bonds_list_indexer; //!< Indexer for accessing the by-tag bonded particle list
+
+    std::shared_ptr<md::NeighborList>
+        m_pair_nlist; //!< The hoomd neighborlist, only used if exclusions of the newly formed bonds
+                      //!< need to be set
+    std::shared_ptr<md::NeighborList> m_pair_internal_nlist;
+    bool m_pair_nlist_exclusions_set; //!< whether or not the bonds are set as exclusions in the
+                                      //!< hoomd particle neighborlist. Set to true when
+                                      //!< m_pair_nlist is set
+
+    //! filter out existing and doublicate bonds from all found possible bonds
+    virtual void filterPossibleBonds();
+
+    bool CheckisExistingLegalBond(Scalar3 i); // this acesses info in m_existing_bonds_list_tag.
+                                              // todo: rename to something sensible
+    void calculateExistingBonds();
+    void makeBonds(uint64_t timestep);
+
+    void AddtoExistingBonds(unsigned int tag1, unsigned int tag2);
+    bool isExistingBond(unsigned int tag1,
+                        unsigned int tag2); // this acesses info in m_existing_bonds_list_tag
+    void checkBoxSize();
+    void checkRcut();
+    void checkBondType();
+    void checkProbability();
+    void setGroupOverlap();
+    void setCutoffs();
+    void checkMaxBondsGroup();
+    void resizePossibleBondlists();
+    void resizeExistingBondList();
+    void allocateParticleArrays();
+
+    //! Notification of a box size change
+    void slotBoxChanged()
+        {
+        m_box_changed = true;
+        }
+
+    //! Notification of total particle number change
+    void slotNumParticlesChanged()
+        {
+        m_max_N_changed = true;
+        }
+
+    bool m_box_changed;   //!< Flag if box dimensions changed
+    bool m_max_N_changed; //!< Flag if total number of particles changed
+    };
+
+namespace detail
+    {
+//! Export the Evaporator to python
+void export_DynamicBondUpdater(pybind11::module& m);
+    } // end namespace detail
+
+    } // end namespace azplugins
+
+    } // end namespace hoomd
+
+#endif // AZPLUGINS_TYPE_UPDATER_H_
diff --git a/src/DynamicBondUpdaterGPU.cc b/src/DynamicBondUpdaterGPU.cc
new file mode 100644
index 0000000..846c895
--- /dev/null
+++ b/src/DynamicBondUpdaterGPU.cc
@@ -0,0 +1,244 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdaterGPU.cc
+ * \brief Definition of DynamicBondUpdaterGPU
+ */
+
+#include "DynamicBondUpdaterGPU.h"
+#include "DynamicBondUpdaterGPU.cuh"
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+
+/*!
+ * \param sysdef System definition
+ *
+ * The system is initialized in a configuration that will be invalid on the
+ * first check of the types and region. This constructor requires that the user
+ * properly initialize the system via setters.
+ */
+DynamicBondUpdaterGPU::DynamicBondUpdaterGPU(std::shared_ptr<SystemDefinition> sysdef,
+                                             std::shared_ptr<Trigger> trigger,
+                                             std::shared_ptr<md::NeighborList> pair_nlist,
+                                             std::shared_ptr<ParticleGroup> group_1,
+                                             std::shared_ptr<ParticleGroup> group_2,
+                                             uint16_t seed)
+    : DynamicBondUpdater(sysdef, trigger, pair_nlist, group_1, group_2, seed),
+      m_num_nonzero_bonds_flag(m_exec_conf), m_max_bonds_overflow_flag(m_exec_conf)
+    {
+    // only one GPU is supported
+    if (!m_exec_conf->isCUDAEnabled())
+        {
+        throw std::runtime_error("Cannot initialize DynamicBondUpdaterGPU on a CPU device.");
+        }
+
+    m_tuner_copy_nlist.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                              m_exec_conf,
+                                              "dynamic_bonding_copy_nlist"));
+    m_tuner_filter_bonds.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                                m_exec_conf,
+                                                "dynamic_bonding_filter_bonds"));
+
+    m_autotuners.insert(m_autotuners.end(), {m_tuner_copy_nlist, m_tuner_filter_bonds});
+
+    m_exec_conf->msg->notice(5) << "Constructing DynamicBondUpdaterGPU" << std::endl;
+    }
+
+DynamicBondUpdaterGPU::DynamicBondUpdaterGPU(std::shared_ptr<SystemDefinition> sysdef,
+                                             std::shared_ptr<Trigger> trigger,
+                                             std::shared_ptr<md::NeighborList> pair_nlist,
+                                             std::shared_ptr<ParticleGroup> group_1,
+                                             std::shared_ptr<ParticleGroup> group_2,
+                                             uint16_t seed,
+                                             const Scalar r_cut,
+                                             const Scalar probability,
+                                             unsigned int max_bonds_group_1,
+                                             unsigned int max_bonds_group_2,
+                                             unsigned int bond_type)
+    : DynamicBondUpdater(sysdef,
+                         trigger,
+                         pair_nlist,
+                         group_1,
+                         group_2,
+                         seed,
+                         r_cut,
+                         probability,
+                         max_bonds_group_1,
+                         max_bonds_group_2,
+                         bond_type),
+      m_num_nonzero_bonds_flag(m_exec_conf), m_max_bonds_overflow_flag(m_exec_conf)
+    {
+    // only one GPU is supported
+    if (!m_exec_conf->isCUDAEnabled())
+        {
+        throw std::runtime_error("Cannot initialize DynamicBondUpdaterGPU on a CPU device.");
+        }
+
+    m_tuner_copy_nlist.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                              m_exec_conf,
+                                              "dynamic_bonding_copy_nlist"));
+    m_tuner_filter_bonds.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)},
+                                                m_exec_conf,
+                                                "dynamic_bonding_filter_bonds"));
+
+    m_autotuners.insert(m_autotuners.end(), {m_tuner_copy_nlist, m_tuner_filter_bonds});
+
+    m_exec_conf->msg->notice(5) << "Constructing DynamicBondUpdaterGPU" << std::endl;
+    }
+
+DynamicBondUpdaterGPU::~DynamicBondUpdaterGPU()
+    {
+    m_exec_conf->msg->notice(5) << "Destroying DynamicBondUpdaterGPU" << std::endl;
+    }
+
+void DynamicBondUpdaterGPU::filterPossibleBonds()
+    {
+    // copy data from m_n_list to d_all_possible_bonds. nlist saves indices and the existing bonds
+    // have to be stored by tags so copy data first then sort out existing bonds
+    const unsigned int size = m_group_1->getNumMembers() * m_max_bonds;
+    ArrayHandle<unsigned int> d_n_existing_bonds(m_n_existing_bonds,
+                                                 access_location::device,
+                                                 access_mode::read);
+    ArrayHandle<unsigned int> d_existing_bonds_list(m_existing_bonds_list,
+                                                    access_location::device,
+                                                    access_mode::read);
+
+    ArrayHandle<unsigned int> d_nlist(m_pair_internal_nlist->getNListArray(),
+                                      access_location::device,
+                                      access_mode::read);
+    ArrayHandle<unsigned int> d_n_neigh(m_pair_internal_nlist->getNNeighArray(),
+                                        access_location::device,
+                                        access_mode::read);
+    ArrayHandle<size_t> d_n_head_list(m_pair_internal_nlist->getHeadList(),
+                                      access_location::device,
+                                      access_mode::read);
+
+    ArrayHandle<unsigned int> d_tag(m_pdata->getTags(), access_location::device, access_mode::read);
+    ArrayHandle<Scalar4> d_pos(m_pdata->getPositions(), access_location::device, access_mode::read);
+    ArrayHandle<unsigned int> d_index_group_1(m_group_1->getIndexArray(),
+                                              access_location::device,
+                                              access_mode::read);
+
+    ArrayHandle<Scalar3> d_all_possible_bonds(m_all_possible_bonds,
+                                              access_location::device,
+                                              access_mode::readwrite);
+    cudaMemset((void*)d_all_possible_bonds.data,
+               0,
+               sizeof(Scalar3) * m_all_possible_bonds.getNumElements());
+
+    const BoxDim& box = m_pdata->getBox();
+
+    if (m_groups_identical)
+        {
+        m_tuner_copy_nlist->begin();
+        gpu::copy_possible_bonds(d_all_possible_bonds.data,
+                                 d_pos.data,
+                                 d_tag.data,
+                                 d_index_group_1.data,
+                                 d_index_group_1.data,
+                                 d_n_neigh.data,
+                                 d_nlist.data,
+                                 d_n_head_list.data,
+                                 box,
+                                 m_max_bonds,
+                                 m_r_cut,
+                                 m_groups_identical,
+                                 m_group_1->getNumMembers(),
+                                 m_group_1->getNumMembers(),
+                                 m_tuner_copy_nlist->getParam()[0]);
+        if (m_exec_conf->isCUDAErrorCheckingEnabled())
+            CHECK_CUDA_ERROR();
+        m_tuner_copy_nlist->end();
+        }
+    else
+        {
+        ArrayHandle<unsigned int> d_index_group_2(m_group_2->getIndexArray(),
+                                                  access_location::device,
+                                                  access_mode::read);
+        m_tuner_copy_nlist->begin();
+        gpu::copy_possible_bonds(d_all_possible_bonds.data,
+                                 d_pos.data,
+                                 d_tag.data,
+                                 d_index_group_1.data,
+                                 d_index_group_2.data,
+                                 d_n_neigh.data,
+                                 d_nlist.data,
+                                 d_n_head_list.data,
+                                 box,
+                                 m_max_bonds,
+                                 m_r_cut,
+                                 m_groups_identical,
+                                 m_group_2->getNumMembers(),
+                                 m_group_1->getNumMembers(),
+                                 m_tuner_copy_nlist->getParam()[0]);
+        if (m_exec_conf->isCUDAErrorCheckingEnabled())
+            CHECK_CUDA_ERROR();
+        m_tuner_copy_nlist->end();
+        }
+
+    // filter out the existing bonds - based on neighbor list exclusion handling
+    m_tuner_filter_bonds->begin();
+    gpu::filter_existing_bonds(d_all_possible_bonds.data,
+                               d_n_existing_bonds.data,
+                               d_existing_bonds_list.data,
+                               m_existing_bonds_list_indexer,
+                               size,
+                               m_tuner_filter_bonds->getParam()[0]);
+    if (m_exec_conf->isCUDAErrorCheckingEnabled())
+        CHECK_CUDA_ERROR();
+    m_tuner_filter_bonds->end();
+
+    m_num_all_possible_bonds = 0;
+
+    gpu::remove_zeros_and_sort_possible_bond_array(d_all_possible_bonds.data,
+                                                   size,
+                                                   m_num_nonzero_bonds_flag.getDeviceFlags());
+
+    m_num_all_possible_bonds = m_num_nonzero_bonds_flag.readFlags();
+
+    // at this point, the sub-array: d_all_possible_bonds[0,m_num_all_possible_bonds]
+    // should contain only unique entries of possible bonds which are not yet formed.
+    }
+
+namespace detail
+    {
+/*!
+ * \param m Python module to export to
+ */
+void export_DynamicBondUpdaterGPU(pybind11::module& m)
+    {
+    namespace py = pybind11;
+
+    py::class_<DynamicBondUpdaterGPU, DynamicBondUpdater, std::shared_ptr<DynamicBondUpdaterGPU>>(
+        m,
+        "DynamicBondUpdaterGPU")
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>,
+                            std::shared_ptr<Trigger>,
+                            std::shared_ptr<md::NeighborList>,
+                            std::shared_ptr<ParticleGroup>,
+                            std::shared_ptr<ParticleGroup>,
+                            uint16_t>())
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>,
+                            std::shared_ptr<Trigger>,
+                            std::shared_ptr<md::NeighborList>,
+                            std::shared_ptr<ParticleGroup>,
+                            std::shared_ptr<ParticleGroup>,
+                            uint16_t,
+                            Scalar,
+                            Scalar,
+                            unsigned int,
+                            unsigned int,
+                            unsigned int>());
+    }
+
+    } // end namespace detail
+    } // end namespace azplugins
+
+    } // end namespace hoomd
diff --git a/src/DynamicBondUpdaterGPU.cu b/src/DynamicBondUpdaterGPU.cu
new file mode 100644
index 0000000..3e3eb1c
--- /dev/null
+++ b/src/DynamicBondUpdaterGPU.cu
@@ -0,0 +1,387 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdaterGPU.cu
+ * \brief Definition of kernel drivers and kernels for DynamicBondUpdaterGPU
+ */
+
+#include "DynamicBondUpdaterGPU.cuh"
+#include "hoomd/HOOMDMath.h"
+#include "hoomd/md/NeighborListGPUTree.cuh"
+
+#include <thrust/device_vector.h>
+#include <thrust/reduce.h>
+#include <thrust/remove.h>
+#include <thrust/sort.h>
+#include <thrust/unique.h>
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+
+// todo: migrate to separate file/class.
+// this sorts according distance, then first tag, then second tag
+struct SortBondsGPU
+    {
+    __host__ __device__ bool operator()(const Scalar3& i, const Scalar3& j)
+        {
+        const Scalar r_sq_1 = i.z;
+        const Scalar r_sq_2 = j.z;
+        if (r_sq_1 == r_sq_2)
+            {
+            const unsigned int tag_11 = __scalar_as_int(i.x);
+            const unsigned int tag_21 = __scalar_as_int(j.x);
+            if (tag_11 == tag_21)
+                {
+                const unsigned int tag_12 = __scalar_as_int(i.y);
+                const unsigned int tag_22 = __scalar_as_int(j.y);
+                return tag_22 > tag_12;
+                }
+            else
+                {
+                return tag_21 > tag_11;
+                }
+            }
+        else
+            {
+            return r_sq_2 > r_sq_1;
+            }
+        }
+    };
+
+// returns true if given possible bond is zero, e.g. (0,0,0.0)
+// possible bonds are ordered, such that tag_a < tag_b in (tag_a,tag_b,rsq)
+// meaning we only need to check tag_b == 0
+struct isZeroBondGPU
+    {
+    __host__ __device__ bool operator()(const Scalar3& i)
+        {
+        const unsigned int tag_1 = __scalar_as_int(i.y);
+        return !(bool)tag_1;
+        }
+    };
+
+struct CompareBondsGPU
+    {
+    __host__ __device__ bool operator()(const Scalar3& i, const Scalar3& j)
+        {
+        const unsigned int tag_11 = __scalar_as_int(i.x);
+        const unsigned int tag_12 = __scalar_as_int(i.y);
+        const unsigned int tag_21 = __scalar_as_int(j.x);
+        const unsigned int tag_22 = __scalar_as_int(j.y);
+
+        if ((tag_11 == tag_21 && tag_12 == tag_22)) // should work because pairs are ordered
+            {
+            return true;
+            }
+        else
+            {
+            return false;
+            }
+        }
+    };
+
+namespace gpu
+    {
+
+//! Number of elements of the exclusion list to process in each batch
+const unsigned int FILTER_BATCH_SIZE = 4;
+
+namespace kernel
+    {
+
+__global__ void copy_nlist_possible_bonds(Scalar3* d_all_possible_bonds,
+                                          const Scalar4* d_postype,
+                                          const unsigned int* d_tag,
+                                          const unsigned int* d_sorted_indexes,
+                                          const unsigned int* d_sorted_indexes_group_2,
+                                          const unsigned int* d_n_neigh,
+                                          const unsigned int* d_nlist,
+                                          const size_t* d_n_head_list,
+                                          const BoxDim box,
+                                          const unsigned int max_bonds,
+                                          const Scalar r_cut,
+                                          const bool groups_identical,
+                                          const unsigned int size_group_2,
+                                          const unsigned int N)
+    {
+    // one thread per particle in group_1
+    const unsigned int idx = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (idx >= N)
+        return;
+
+    // idx = group index , pidx = actual particle index
+    const unsigned int pidx_i = d_sorted_indexes[idx];
+    unsigned int n_curr_bond = 0;
+    const Scalar r_cutsq = r_cut * r_cut;
+
+    // get all information for this particle
+    Scalar4 postype_i = d_postype[pidx_i];
+    const unsigned int tag_i = d_tag[pidx_i];
+    const unsigned int n_neigh = d_n_neigh[pidx_i];
+    const size_t head = d_n_head_list[pidx_i];
+
+    // loop over all neighbors of this particle
+    for (unsigned int j = 0; j < n_neigh; ++j)
+        {
+        // get index of neighbor from neigh_list
+        const unsigned int pidx_j = d_nlist[head + j];
+
+        // test if j is in group 2
+
+        // wrapper for pointer needed for thrust
+        thrust::device_ptr<const unsigned int> d_sorted_indexes_group_2_wrap(
+            d_sorted_indexes_group_2);
+
+        auto iter = thrust::find(thrust::device,
+                                 d_sorted_indexes_group_2_wrap,
+                                 d_sorted_indexes_group_2_wrap + size_group_2,
+                                 pidx_j);
+
+        if (iter != d_sorted_indexes_group_2_wrap + size_group_2)
+            {
+            Scalar4 postype_j = d_postype[pidx_j];
+            const unsigned int tag_j = d_tag[pidx_j];
+
+            Scalar3 drij = make_scalar3(postype_j.x, postype_j.y, postype_j.z)
+                           - make_scalar3(postype_i.x, postype_i.y, postype_i.z);
+
+            // apply periodic boundary conditions (FLOPS: 12)
+            drij = box.minImage(drij);
+
+            // same as on the cpu, just not during the tree traversal
+            Scalar dr_sq = dot(drij, drij);
+
+            if (dr_sq < r_cutsq)
+                {
+                if (n_curr_bond < max_bonds)
+                    {
+                    Scalar3 d;
+                    if (groups_identical)
+                        {
+                        // sort the two tags in this possible bond pair if groups are identical
+                        const unsigned int tag_a = tag_j > tag_i ? tag_i : tag_j;
+                        const unsigned int tag_b = tag_j > tag_i ? tag_j : tag_i;
+                        d = make_scalar3(__int_as_scalar(tag_a), __int_as_scalar(tag_b), dr_sq);
+                        }
+                    else
+                        {
+                        d = make_scalar3(__int_as_scalar(tag_i), __int_as_scalar(tag_j), dr_sq);
+                        }
+                    d_all_possible_bonds[idx * max_bonds + n_curr_bond] = d;
+                    }
+                ++n_curr_bond;
+                }
+            }
+        }
+    }
+
+/*! \param d_all_possible_bonds all possible bonds list
+    \param d_n_existing_bonds Number of existing for each particle
+    \param d_existing_bonds_list List of exitsting for each particle
+    \param exli Indexer for indexing into d_existing_bonds_list
+    \param size Length of d_all_possible_bonds
+    \param ex_start Start filtering  d_all_possible_bonds from existing bond number \a ex_start
+
+    the kernel filter_existing_bonds() processes the all possible bonds list \a d_nlist and removes
+   any entries that already exist. To allow for an arbitrary large number of existing bonds, these
+   are processed in batch sizes of FILTER_BATCH_SIZE. The kernel must be called multiple times in
+   order to fully remove all exclusions from the nlist.
+
+    \note The driver filter_existing_bonds properly makes as many calls as are necessary, it only
+   needs to be called once.
+
+    \b Implementation
+
+    One thread is run for each particle. Existing bonds \a ex_start, \a ex_start + 1, ... are loaded
+   in for that particle (or the thread returns if there are no exclusions past that point). The
+   thread then loops over the neighbor list, comparing each entry to the list of exclusions. If the
+   entry is not excluded, it is written back out. \a d_n_neigh is updated to reflect the current
+   number of particles in the list at the end of the kernel call.
+*/
+__global__ void filter_existing_bonds(Scalar3* d_all_possible_bonds,
+                                      const unsigned int* d_n_existing_bonds,
+                                      const unsigned int* d_existing_bonds_list,
+                                      const Index2D exli,
+                                      const unsigned int size,
+                                      const unsigned int ex_start)
+    {
+    // compute the bond index this thread operates on
+    const unsigned int idx = blockDim.x * blockIdx.x + threadIdx.x;
+
+    // quit now if this thread is processing past the end of the list of all possible bonds
+    if (idx >= size)
+        return;
+
+    Scalar3 current_bond = d_all_possible_bonds[idx];
+    unsigned int tag_1 = __scalar_as_int(current_bond.x);
+    unsigned int tag_2 = __scalar_as_int(current_bond.y);
+
+    if (tag_1 == 0 && tag_2 == 0)
+        return;
+
+    // const unsigned int n_neigh = d_n_neigh[idx];
+    const unsigned int n_ex = d_n_existing_bonds[tag_1];
+
+    // quit now if the ex_start flag is past the end of n_ex
+    if (ex_start >= n_ex)
+        return;
+
+    // count the number of existing bonds to process in this thread
+    const unsigned int n_ex_process = n_ex - ex_start;
+
+    // load the existing bond list into "local" memory - fully unrolled loops should dump this into
+    // registers
+    unsigned int l_existing_bonds_list[FILTER_BATCH_SIZE];
+#pragma unroll
+    for (unsigned int cur_ex_idx = 0; cur_ex_idx < FILTER_BATCH_SIZE; cur_ex_idx++)
+        {
+        if (cur_ex_idx < n_ex_process)
+            {
+            l_existing_bonds_list[cur_ex_idx]
+                = d_existing_bonds_list[exli(tag_1, cur_ex_idx + ex_start)];
+            }
+        else
+            {
+            l_existing_bonds_list[cur_ex_idx] = 0xffffffff;
+            }
+        }
+
+    // test if excluded
+    bool excluded = false;
+#pragma unroll
+    for (unsigned int cur_ex_idx = 0; cur_ex_idx < FILTER_BATCH_SIZE; cur_ex_idx++)
+        {
+        if (tag_2 == l_existing_bonds_list[cur_ex_idx])
+            excluded = true;
+        }
+    // if it is excluded, overwrite that entry with (0,0,0).
+    if (excluded)
+        {
+        d_all_possible_bonds[idx] = make_scalar3(__int_as_scalar(0), __int_as_scalar(0), 0.0);
+        }
+    }
+
+    } // end namespace kernel
+
+cudaError_t remove_zeros_and_sort_possible_bond_array(Scalar3* d_all_possible_bonds,
+                                                      const unsigned int size,
+                                                      int* d_max_non_zero_bonds)
+    {
+    if (size == 0)
+        return cudaSuccess;
+    // wrapper for pointer needed for thrust
+    thrust::device_ptr<Scalar3> d_all_possible_bonds_wrap(d_all_possible_bonds);
+
+    isZeroBondGPU zero;
+    thrust::device_ptr<Scalar3> last0
+        = thrust::remove_if(d_all_possible_bonds_wrap, d_all_possible_bonds_wrap + size, zero);
+    unsigned int l0 = thrust::distance(d_all_possible_bonds_wrap, last0);
+
+    // sort remainder by distance, should make all identical bonds consequtive
+    SortBondsGPU sort;
+    thrust::sort(d_all_possible_bonds_wrap, d_all_possible_bonds_wrap + l0, sort);
+
+    // thrust::unique only removes identical consequtive elements, so sort above is needed.
+    CompareBondsGPU comp;
+    thrust::device_ptr<Scalar3> last1
+        = thrust::unique(d_all_possible_bonds_wrap, d_all_possible_bonds_wrap + l0, comp);
+    unsigned int l1 = thrust::distance(d_all_possible_bonds_wrap, last1);
+
+    *d_max_non_zero_bonds = l1;
+
+    return cudaSuccess;
+    }
+
+cudaError_t filter_existing_bonds(Scalar3* d_all_possible_bonds,
+                                  unsigned int* d_n_existing_bonds,
+                                  const unsigned int* d_existing_bonds_list,
+                                  const Index2D& exli,
+                                  const unsigned int size,
+                                  const unsigned int block_size)
+    {
+    static unsigned int max_block_size = UINT_MAX;
+    if (max_block_size == UINT_MAX)
+        {
+        cudaFuncAttributes attr;
+        cudaFuncGetAttributes(&attr, (const void*)kernel::filter_existing_bonds);
+        max_block_size = attr.maxThreadsPerBlock;
+        }
+
+    unsigned int run_block_size = min(block_size, max_block_size);
+
+    // determine parameters for kernel launch
+    int n_blocks = size / run_block_size + 1;
+
+    // split the processing of the full exclusion list up into a number of batches
+    unsigned int n_batches = (unsigned int)ceil(double(exli.getH()) / double(FILTER_BATCH_SIZE));
+    unsigned int ex_start = 0;
+    for (unsigned int batch = 0; batch < n_batches; batch++)
+        {
+        kernel::filter_existing_bonds<<<n_blocks, run_block_size>>>(d_all_possible_bonds,
+                                                                    d_n_existing_bonds,
+                                                                    d_existing_bonds_list,
+                                                                    exli,
+                                                                    size,
+                                                                    ex_start);
+
+        ex_start += FILTER_BATCH_SIZE;
+        }
+
+    return cudaSuccess;
+    }
+
+cudaError_t copy_possible_bonds(Scalar3* d_all_possible_bonds,
+                                const Scalar4* d_postype,
+                                const unsigned int* d_tag,
+                                const unsigned int* d_sorted_indexes,
+                                const unsigned int* d_sorted_indexes_group_2,
+                                const unsigned int* d_n_neigh,
+                                const unsigned int* d_nlist,
+                                const size_t* d_n_head_list,
+                                const BoxDim box,
+                                const unsigned int max_bonds,
+                                const Scalar r_cut,
+                                const bool groups_identical,
+                                const unsigned int size_group_2,
+                                const unsigned int N,
+                                const unsigned int block_size)
+    {
+    static unsigned int max_block_size = UINT_MAX;
+
+    if (max_block_size == UINT_MAX)
+        {
+        cudaFuncAttributes attr;
+        cudaFuncGetAttributes(&attr, (const void*)kernel::copy_nlist_possible_bonds);
+        max_block_size = attr.maxThreadsPerBlock;
+        }
+    unsigned int run_block_size = min(block_size, max_block_size);
+
+    kernel::copy_nlist_possible_bonds<<<N / run_block_size + 1, run_block_size>>>(
+        d_all_possible_bonds,
+        d_postype,
+        d_tag,
+        d_sorted_indexes,
+        d_sorted_indexes_group_2,
+        d_n_neigh,
+        d_nlist,
+        d_n_head_list,
+        box,
+        max_bonds,
+        r_cut,
+        groups_identical,
+        size_group_2,
+        N);
+    return cudaSuccess;
+    }
+
+    } // end namespace gpu
+    } // end namespace azplugins
+
+    } // end namespace hoomd
diff --git a/src/DynamicBondUpdaterGPU.cuh b/src/DynamicBondUpdaterGPU.cuh
new file mode 100644
index 0000000..74d13bf
--- /dev/null
+++ b/src/DynamicBondUpdaterGPU.cuh
@@ -0,0 +1,78 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdaterGPU.cuh
+ * \brief Declaration of kernel drivers for DynamicBondUpdaterGPU
+ */
+
+#ifndef AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_CUH_
+#define AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_CUH_
+
+#include "hoomd/BoxDim.h"
+#include "hoomd/HOOMDMath.h"
+#include "hoomd/Index1D.h"
+#include "hoomd/ParticleData.cuh"
+#include <cuda_runtime.h>
+#include <iostream>
+
+#include "hoomd/md/NeighborListGPUTree.cuh"
+
+#ifdef NVCC
+#define DEVICE __device__ __forceinline__
+#define HOSTDEVICE __host__ __device__ __forceinline__
+#else
+#define DEVICE
+#define HOSTDEVICE
+#endif
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+namespace gpu
+    {
+
+cudaError_t sort_possible_bond_array(Scalar3* d_all_possible_bonds, const unsigned int size);
+
+cudaError_t filter_existing_bonds(Scalar3* d_all_possible_bonds,
+                                  unsigned int* d_n_existing_bonds,
+                                  const unsigned int* d_existing_bonds_list,
+                                  const Index2D& exli,
+                                  const unsigned int size,
+                                  const unsigned int block_size);
+
+cudaError_t copy_possible_bonds(Scalar3* d_all_possible_bonds,
+                                const Scalar4* d_postype,
+                                const unsigned int* d_tag,
+                                const unsigned int* d_sorted_indexes,
+                                const unsigned int* d_sorted_indexes_group_2,
+                                const unsigned int* d_n_neigh,
+                                const unsigned int* d_nlist,
+                                const size_t* d_n_head_list,
+                                const BoxDim box,
+                                const unsigned int max_bonds,
+                                const Scalar r_cut,
+                                const bool groups_identical,
+                                const unsigned int size_group_2,
+                                const unsigned int N,
+                                const unsigned int block_size);
+
+cudaError_t remove_zeros_and_sort_possible_bond_array(Scalar3* d_all_possible_bonds,
+                                                      const unsigned int size,
+                                                      int* d_max_non_zero_bonds);
+
+//! Sentinel for an invalid particle (e.g., ghost)
+const unsigned int NeighborListTypeSentinel = 0xffffffff;
+
+    } // end namespace gpu
+    } // end namespace azplugins
+    } // end namespace hoomd
+
+#undef DEVICE
+#undef HOSTDEVICE
+
+#endif // AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_CUH_
diff --git a/src/DynamicBondUpdaterGPU.h b/src/DynamicBondUpdaterGPU.h
new file mode 100644
index 0000000..ebecf88
--- /dev/null
+++ b/src/DynamicBondUpdaterGPU.h
@@ -0,0 +1,84 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2025, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+// Maintainer: astatt
+
+/*!
+ * \file DynamicBondUpdaterGPU.h
+ * \brief Declaration of DynamicBondUpdaterGPU.h
+ */
+
+#ifndef AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_H_
+#define AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_H_
+
+#ifdef NVCC
+#error This header cannot be compiled by nvcc
+#endif
+
+#include "DynamicBondUpdater.h"
+#include "DynamicBondUpdaterGPU.cuh"
+#include "hoomd/Autotuner.h"
+
+#include "hip/hip_runtime.h"
+#include "hoomd/md/NeighborListGPUTree.cuh"
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+
+//! Particle type updater on the GPU
+/*!
+ * See DynamicBondUpdater for details. This class inherits and minimally implements
+ * the CPU methods from DynamicBondUpdater on the GPU.
+ */
+class PYBIND11_EXPORT DynamicBondUpdaterGPU : public DynamicBondUpdater
+    {
+    public:
+    //! Simple constructor
+    DynamicBondUpdaterGPU(std::shared_ptr<SystemDefinition> sysdef,
+                          std::shared_ptr<Trigger> trigger,
+                          std::shared_ptr<md::NeighborList> pair_nlist,
+                          std::shared_ptr<ParticleGroup> group_1,
+                          std::shared_ptr<ParticleGroup> group_2,
+                          uint16_t seed);
+
+    //! Constructor with parameters
+    DynamicBondUpdaterGPU(std::shared_ptr<SystemDefinition> sysdef,
+                          std::shared_ptr<Trigger> trigger,
+                          std::shared_ptr<md::NeighborList> pair_nlist,
+                          std::shared_ptr<ParticleGroup> group_1,
+                          std::shared_ptr<ParticleGroup> group_2,
+                          uint16_t seed,
+                          const Scalar r_cut,
+                          const Scalar probability,
+                          unsigned int max_bonds_group_1,
+                          unsigned int max_bonds_group_2,
+                          unsigned int bond_type);
+
+    //! Destructor
+    virtual ~DynamicBondUpdaterGPU();
+
+    protected:
+    //! filter out existing and doublicate bonds from all found possible bonds
+    virtual void filterPossibleBonds();
+
+    private:
+    std::shared_ptr<Autotuner<1>> m_tuner_copy_nlist;   //!< Tuner for the primitive-copy kernel
+    std::shared_ptr<Autotuner<1>> m_tuner_filter_bonds; //!< Tuner for existing bond filter
+
+    GPUFlags<int> m_num_nonzero_bonds_flag;           //!< GPU flag for the number of valid bonds
+    GPUFlags<unsigned int> m_max_bonds_overflow_flag; //!< GPU flag for overflow
+    };
+
+namespace detail
+    {
+//! Export DynamicBondUpdaterGPU to python
+void export_DynamicBondUpdaterGPU(pybind11::module& m);
+    } // end namespace detail
+
+    } // end namespace azplugins
+    } // end namespace hoomd
+
+#endif // AZPLUGINS_DYNAMIC_BOND_UPDATER_GPU_H_
diff --git a/src/RNGIdentifiers.h b/src/RNGIdentifiers.h
index 283f5cb..d7df23f 100644
--- a/src/RNGIdentifiers.h
+++ b/src/RNGIdentifiers.h
@@ -24,6 +24,7 @@ struct RNGIdentifier
     static const uint8_t TwoStepBrownianFlow = 201;
     static const uint8_t TwoStepLangevinFlow = 202;
     static const uint8_t ParticleEvaporator = 203;
+    static const uint8_t DynamicBondUpdater = 204;
     };
 
     } // end namespace detail
diff --git a/src/__init__.py b/src/__init__.py
index 3e71abf..70b9f6e 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -4,6 +4,6 @@
 
 """azplugins."""
 
-from hoomd.azplugins import bond, compute, external, flow, pair, wall
+from hoomd.azplugins import bond, compute, external, flow, update, pair, wall
 
 __version__ = "1.2.0"
diff --git a/src/module.cc b/src/module.cc
index 1ab81c7..6c42349 100644
--- a/src/module.cc
+++ b/src/module.cc
@@ -77,6 +77,9 @@ void export_PotentialPairPerturbedLennardJones(pybind11::module&);
 // dpd
 void export_PotentialPairDPDThermoGeneralWeight(pybind11::module&);
 
+// update
+void export_DynamicBondUpdater(pybind11::module&);
+
 // wall
 void export_PotentialExternalWallColloid(pybind11::module&);
 void export_PotentialExternalWallLJ93(pybind11::module&);
@@ -106,6 +109,9 @@ void export_PotentialPairPerturbedLennardJonesGPU(pybind11::module&);
 // dpd
 void export_PotentialPairDPDThermoGeneralWeightGPU(pybind11::module&);
 
+// update
+void export_DynamicBondUpdaterGPU(pybind11::module&);
+
 // wall
 void export_PotentialExternalWallColloidGPU(pybind11::module&);
 void export_PotentialExternalWallLJ93GPU(pybind11::module&);
@@ -149,6 +155,9 @@ PYBIND11_MODULE(_azplugins, m)
     // dpd pair
     export_PotentialPairDPDThermoGeneralWeight(m);
 
+    // updater
+    export_DynamicBondUpdater(m);
+
     // wall
     export_PotentialExternalWallColloid(m);
     export_PotentialExternalWallLJ93(m);
@@ -178,6 +187,9 @@ PYBIND11_MODULE(_azplugins, m)
     // dpd pair
     export_PotentialPairDPDThermoGeneralWeightGPU(m);
 
+    // updater
+    export_DynamicBondUpdaterGPU(m);
+
     // wall
     export_PotentialExternalWallColloidGPU(m);
     export_PotentialExternalWallLJ93GPU(m);
diff --git a/src/pytest/CMakeLists.txt b/src/pytest/CMakeLists.txt
index 157b326..62b99a4 100644
--- a/src/pytest/CMakeLists.txt
+++ b/src/pytest/CMakeLists.txt
@@ -3,6 +3,7 @@ set(test_files
     __init__.py
     test_bond.py
     test_compute.py
+    test_dynamic_bond.py
     test_external.py
     test_flow.py
     test_pair.py
diff --git a/src/pytest/test_dynamic_bond.py b/src/pytest/test_dynamic_bond.py
new file mode 100644
index 0000000..009bfe9
--- /dev/null
+++ b/src/pytest/test_dynamic_bond.py
@@ -0,0 +1,201 @@
+# Copyright (c) 2018-2020, Michael P. Howard
+# Copyright (c) 2021-2025, Auburn University
+# Part of azplugins, released under the BSD 3-Clause License.
+
+import hoomd
+import numpy
+
+
+def test_setters_getters(simulation_factory, one_particle_snapshot_factory):
+    # make one particle test configuration
+    simulation_factory(one_particle_snapshot_factory(position=[0, 0, 0], L=20))
+
+    u = hoomd.azplugins.update.DynamicBond(
+        nlist=hoomd.md.nlist.Cell(buffer=0.4),
+        r_cut=1,
+        trigger=hoomd.trigger.Periodic(period=10),
+        bond_type=0,
+        group_1=hoomd.filter.All(),
+        group_2=hoomd.filter.All(),
+        max_bonds_group_1=0,
+        max_bonds_group_2=0,
+    )
+
+    assert numpy.equal(u.r_cut, 1)
+    u.r_cut = 1.5
+    assert numpy.equal(u.r_cut, 1.5)
+
+    u.max_bonds_group_1 = 3
+    assert numpy.equal(u.max_bonds_group_1, 3)
+
+    u.max_bonds_group_2 = 7
+    assert numpy.equal(u.max_bonds_group_2, 7)
+
+
+def test_form_bonds_same_group(simulation_factory):
+    snap = hoomd.Snapshot()
+    if snap.communicator.rank == 0:
+        snap.configuration.box = [20, 20, 20, 0, 0, 0]
+        snap.particles.N = 6
+        snap.particles.types = ["A", "B"]
+        snap.particles.position[:, 0] = (0, 1, 2, 3, 4, 6)
+        snap.particles.position[:, 1] = (0, 0, 0, 0, 0, 0)
+        snap.particles.position[:, 2] = (0, 0, 0, 0, 0, 0)
+        # dynamic bond operates on groups, so typeids should not matter at all
+        snap.particles.typeid[:] = [0, 0, 1, 1, 0, 0]
+        snap.bonds.types = ["bond"]
+
+    sim = simulation_factory(snap)
+    nl = hoomd.md.nlist.Cell(buffer=0.4)
+
+    group = hoomd.filter.Tags([0, 1, 2, 3, 4, 5])
+    form_bonds = hoomd.azplugins.update.DynamicBond(
+        nlist=nl,
+        r_cut=1.1,
+        trigger=hoomd.trigger.Periodic(period=1),
+        bond_type=0,
+        group_1=group,
+        group_2=group,
+        max_bonds_group_1=2,
+        max_bonds_group_2=2,
+    )
+
+    # test bond formation. All are in the same group, so bonds should be formed
+    # between 0-1, 1-2, 2-3, and 3-4 (but not 5, too far away)
+    integrator = hoomd.md.Integrator(dt=0)
+    nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
+    integrator.methods.append(nve)
+    sim.operations.updaters.append(form_bonds)
+    sim.run(1)
+    s = sim.state.get_snapshot()
+
+    numpy.testing.assert_almost_equal(4, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([0, 1], s.bonds.group[0])
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[1])
+    numpy.testing.assert_array_almost_equal([2, 3], s.bonds.group[2])
+    numpy.testing.assert_array_almost_equal([3, 4], s.bonds.group[3])
+
+    # same as before, we will have formed bonds 0-1-2-3-4 in the first step, so nothing
+    # new should be formed
+    sim.run(1)
+    s = sim.state.get_snapshot()
+
+    numpy.testing.assert_almost_equal(4, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([0, 1], s.bonds.group[0])
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[1])
+    numpy.testing.assert_array_almost_equal([2, 3], s.bonds.group[2])
+    numpy.testing.assert_array_almost_equal([3, 4], s.bonds.group[3])
+
+    # now put particle 5 in range of partice 2 and 3, but no new bonds should
+    # be formed since max_bonds_group_1=max_bonds_group_2=2 and that would be the third
+    # bond on those particles
+    s.particles.position[5] = (2.5, 0, 0)
+    sim.state.set_snapshot(s)
+    sim.run(1)
+    s = sim.state.get_snapshot()
+
+    numpy.testing.assert_almost_equal(4, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([0, 1], s.bonds.group[0])
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[1])
+    numpy.testing.assert_array_almost_equal([2, 3], s.bonds.group[2])
+    numpy.testing.assert_array_almost_equal([3, 4], s.bonds.group[3])
+
+
+def test_form_bonds_same_group_priority(simulation_factory):
+    snap = hoomd.Snapshot()
+    if snap.communicator.rank == 0:
+        snap.configuration.box = [20, 20, 20, 0, 0, 0]
+        snap.particles.N = 6
+        snap.particles.types = ["A", "B"]
+        snap.particles.position[:, 0] = (0, 1, 2, 3, 4, 2.5)
+        snap.particles.position[:, 1] = (0, 0, 0, 0, 0, 0)
+        snap.particles.position[:, 2] = (0, 0, 0, 0, 0, 0)
+        # dynamic bond operates on groups, so typeids should not matter at all
+        snap.particles.typeid[:] = [0, 0, 1, 1, 0, 0]
+        snap.bonds.types = ["bond"]
+
+    sim = simulation_factory(snap)
+    nl = hoomd.md.nlist.Cell(buffer=0.4)
+
+    group = hoomd.filter.Tags([0, 1, 2, 3, 4, 5])
+    form_bonds = hoomd.azplugins.update.DynamicBond(
+        nlist=nl,
+        r_cut=1.1,
+        trigger=hoomd.trigger.Periodic(period=1),
+        bond_type=0,
+        group_1=group,
+        group_2=group,
+        max_bonds_group_1=2,
+        max_bonds_group_2=2,
+    )
+
+    # test bond formation. All are in the same group, so bonds should be formed
+    # between 0-1, 1-2, 2-3, and 3-4 (but not 5, too far away)
+    integrator = hoomd.md.Integrator(dt=0)
+    nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
+    integrator.methods.append(nve)
+    sim.operations.updaters.append(form_bonds)
+    sim.run(1)
+    s = sim.state.get_snapshot()
+
+    # here, particle 5 is in range of partice 2 and 3, so in principle the bonds
+    # 0-1, 1-2 ,2-3, 3-4 (all length 1), 2-5 and 3-5 (length 0.5) can be formed.
+    # The all_possible_bonds array  is sorted by bond distance, so bonds 2-5
+    # and 3-5 are formed first. Then in order of index 0-1, and 1-2.
+    # Particle 2 now has two bonds, so 2-3 can't be formed, but 3-4 can.
+    numpy.testing.assert_almost_equal(5, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([2, 5], s.bonds.group[0])
+    numpy.testing.assert_array_almost_equal([3, 5], s.bonds.group[1])
+    numpy.testing.assert_array_almost_equal([0, 1], s.bonds.group[2])
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[3])
+    numpy.testing.assert_array_almost_equal([3, 4], s.bonds.group[4])
+
+
+def test_update_bond_two_groups(simulation_factory):
+    snap = hoomd.Snapshot()
+    if snap.communicator.rank == 0:
+        snap.configuration.box = [20, 20, 20, 0, 0, 0]
+        snap.particles.N = 4
+        snap.particles.types = ["A", "B"]
+        snap.particles.position[:, 0] = (0, 0.9, 1.1, 2)
+        snap.particles.position[:, 1] = (0, 0, 0, 0)
+        snap.particles.position[:, 2] = (0, 0, 0, 0)
+        snap.particles.typeid[:] = [1, 0, 1, 0]
+        snap.bonds.types = ["bond"]
+
+    sim = simulation_factory(snap)
+    nl = hoomd.md.nlist.Cell(buffer=0.4)
+
+    group_1 = hoomd.filter.Tags([0, 1])
+    group_2 = hoomd.filter.Tags([2, 3])
+    form_bonds = hoomd.azplugins.update.DynamicBond(
+        nlist=nl,
+        r_cut=1.0,
+        trigger=hoomd.trigger.Periodic(period=1),
+        bond_type=0,
+        group_1=group_1,
+        group_2=group_2,
+        max_bonds_group_1=1,
+        max_bonds_group_2=2,
+    )
+
+    # test bond formation between particle 1-2
+    # particle 0,1 are in the same group, so even if their distance is
+    # below r_cut, no bond should be formed, same for goes for particle 2,3
+    integrator = hoomd.md.Integrator(dt=0)
+    nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
+    integrator.methods.append(nve)
+    sim.operations.updaters.append(form_bonds)
+    sim.run(1)
+    s = sim.state.get_snapshot()
+
+    numpy.testing.assert_almost_equal(1, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[0])
+
+    # after this bond 1-2 is formed, there shouldn't be a second one formed
+    # e.g. a dublicate, even when the simulation continues to run
+    sim.run(10)
+    s = sim.state.get_snapshot()
+
+    numpy.testing.assert_almost_equal(1, s.bonds.N)
+    numpy.testing.assert_array_almost_equal([1, 2], s.bonds.group[0])
diff --git a/src/update.py b/src/update.py
new file mode 100644
index 0000000..a304c51
--- /dev/null
+++ b/src/update.py
@@ -0,0 +1,190 @@
+# Copyright (c) 2018-2020, Michael P. Howard
+# Copyright (c) 2021-2025, Auburn University
+# Part of azplugins, released under the BSD 3-Clause License.
+
+"""Updaters."""
+
+import hoomd
+
+from hoomd.azplugins import _azplugins
+from hoomd.data.parameterdicts import ParameterDict
+from hoomd.data.typeconverter import OnlyTypes
+
+
+class DynamicBond(hoomd.operation.Updater):
+    R"""Update bonds dynamically during simulation.
+
+    Args:
+        r_cut (float): Distance cutoff for making bonds between particles
+        probability (float): Probability of bond formation, between 0 and 1, default = 1
+        bond_type (str): Type of bond to be formed
+        group_1 (:py:mod: `hoomd.filter.ParticleFilter`): First particle group
+        group_2 (:py:mod:`hoomd.filter.ParticleFilter`): Second particle group
+        max_bonds_1 (int): Maximum number of bonds a particle in group_1 can have
+        max_bonds_2 (int): Maximum number of bonds a particle in group_2 can have
+        seed (int): Seed to the pseudo-random number generator
+        nlist (:py:mod:`hoomd.md.nlist`): NeighborList (optional) to update exclusions
+        period (int): Particle types will be updated every *period* time steps
+        phase (int): When -1, start on the current time step. Otherwise, execute
+                     on steps where *(step + phase) % period* is 0.
+
+        Forms bonds of type bond_type between particles in group_1 and group_2 during
+        the simulation, if particle distances are shorter than r_cut. If the
+        neighborlist used for the pair potential in the simulation is given as a
+        parameter nlist, the neighbor list exclusions will be updated to include the
+        newly formed bonds. Each particle has a number of maximum bonds which it can
+        form, given by max_bonds_1 for particles in group_1 and max_bonds_2 for group_2.
+
+        The particles in the two groups group_1 and group_2 should be completely
+        separate with no common elements, e.g. two different types, or the two
+        groups should be identical, where now max_bonds_1 needs to be equal to
+        max_bonds_2.
+
+
+        .. warning::
+            If the groups group_1 and group_2 are modified during the simulation, this
+            Updater will not be updated to reflect the changes. It is the user's
+            responsibility to ensure that the groups do not change as long as this
+            updater is active.
+
+        Examples::
+
+            azplugins.update.dynamic_bond(
+                nlist=nl,
+                r_cut=1.0,
+                bond_type=0,
+                group_1=hoomd.group.type(type="A"),
+                group_2=hoomd.group.type(type="B"),
+                max_bonds_1=3,
+                max_bonds_2=2,
+            )
+    """
+
+    _ext_module = _azplugins
+    _cpp_class_name = "DynamicBondUpdater"
+
+    def __init__(
+        self,
+        trigger,
+        nlist,
+        group_1,
+        group_2,
+        bond_type=None,
+        max_bonds_group_1=None,
+        max_bonds_group_2=None,
+        r_cut=None,
+        seed=0,
+        probability=1,
+    ):
+        super().__init__(trigger)
+
+        params = ParameterDict(
+            r_cut=OnlyTypes(float, allow_none=True),
+            nlist=OnlyTypes(hoomd.md.nlist.NeighborList, strict=True, allow_none=False),
+            group_1=OnlyTypes(hoomd.filter.ParticleFilter, allow_none=False),
+            group_2=OnlyTypes(hoomd.filter.ParticleFilter, allow_none=False),
+            max_bonds_group_1=OnlyTypes(int, allow_none=True),
+            max_bonds_group_2=OnlyTypes(int, allow_none=True),
+            bond_type=OnlyTypes(int, strict=True, allow_none=True),
+            seed=OnlyTypes(int, strict=True, allow_none=False),
+            probability=float(probability),
+        )
+
+        params.update(
+            dict(
+                r_cut=r_cut,
+                nlist=nlist,
+                group_1=group_1,
+                group_2=group_2,
+                max_bonds_group_1=max_bonds_group_1,
+                max_bonds_group_2=max_bonds_group_2,
+                bond_type=bond_type,
+                seed=seed,
+                probability=probability,
+            )
+        )
+
+        self._param_dict.update(params)
+        # self.set_params(r_cut,probability,bond_type,max_bonds_1,max_bonds_2,nlist)
+
+    @property
+    def bond_type(self):
+        """bond_type (int): Type of the bonds that are made by this updater."""
+        return self._cpp_obj.bond_type
+
+    @bond_type.setter
+    def bond_type(self, value):
+        if value is not None:
+            self._param_dict["bond_type"] = value
+            self._cpp_obj.setBondType(value)
+
+    @property
+    def probability(self):
+        """Probability (float): Probability of forming bonds. Value between 0 and 1."""
+        return self._cpp_obj.probability
+
+    @probability.setter
+    def probability(self, value):
+        self._param_dict["probability"] = value
+        self._cpp_obj.probability = value
+
+    @property
+    def r_cut(self):
+        """r_cut (float): Distance cutoff for making bonds between particles."""
+        return self._cpp_obj.r_cut
+
+    @r_cut.setter
+    def r_cut(self, value):
+        self._cpp_obj.r_cut = value
+        self._param_dict["r_cut"] = value
+
+    @property
+    def max_bonds_group_1(self):
+        """max_bonds_1 (int): Max bonds on group 1."""
+        return self._cpp_obj.max_bonds_group_1
+
+    @max_bonds_group_1.setter
+    def max_bonds_group_1(self, value):
+        self._cpp_obj.max_bonds_group_1 = value
+        self._param_dict["max_bonds_group_1"] = value
+
+    @property
+    def max_bonds_group_2(self):
+        """max_bonds_2 (int): Max bonds on group 2."""
+        return self._cpp_obj.max_bonds_group_2
+
+    @max_bonds_group_2.setter
+    def max_bonds_group_2(self, value):
+        self._cpp_obj.max_bonds_group_2 = value
+        self._param_dict["max_bonds_group_2"] = value
+
+    def _attach_hook(self):
+        """Create the c++ mirror class."""
+        sim = self._simulation
+        if isinstance(self._simulation.device, hoomd.device.CPU):
+            cpp_class = getattr(self._ext_module, self._cpp_class_name)
+        else:
+            cpp_class = getattr(self._ext_module, self._cpp_class_name + "GPU")
+
+        if self.group_1 is not None:
+            group_1 = sim.state._get_group(self.group_1)
+        else:
+            group_1 = None
+
+        if self.group_2 is not None:
+            group_2 = sim.state._get_group(self.group_2)
+        else:
+            group_2 = None
+
+        self.nlist._attach(sim)
+
+        self._cpp_obj = cpp_class(
+            sim.state._cpp_sys_def,
+            self.trigger,
+            self.nlist._cpp_obj,
+            group_1,
+            group_2,
+            self.seed,
+        )
+
+        super()._attach_hook()