Associative bond swaps in molecular dynamics

The implementation of the method is useful.

Connections with experiments is not provided.

Authors implemented a three-body potential used for simulating vitrimers in HOOMD-blue software. As the authors mention, the method does not require using an MC step to alter the topology. The analyses presented in the paper are convincing, and I find the implementation in the software very useful. The manuscript is written clearly and in an acceptable format.

(1) One area that needs to be addressed is the connection with experiments. Specifically, how should one use this simulation toolkit and get consistent rheological data as obtained by Leibler and coworkers?
(2) As a follow-up to the previous comment, what is the MD time here? Is this the dimensionless time based on the Lennard-Jones potential parameters? It would be helpful if the authors showed the equations of motion and compared the time units with the real-time.
(3) The reference cited in Figure 3 and its caption does not match.

1- Well written, clear and to the point
2- This implementation will be very useful for simulations of network-forming systems

1- It does not fulfill any expectations listed in the journal's "acceptance criteria"
2- the "pressure" subsection seems a bit disconnected from the rest of the paper
3- the quality of some of the figures may be improved

The package developed by the authors seems a very useful addition to the toolset of a molecular simulator interested in network-forming systems, and as such I believe that the paper should be published. However, looking at the acceptance criteria listed on the SciPost Physics's website I cannot see this paper meeting any of the "expectations". Looking at the requirements I believe that a more appropriate journal for this paper would be SciPost Physics Codebases, and this is why I'm ticking the "Accept in alternative journal" box.
Now about the paper content. The paper is clear and well-written, and I have only a few comments that the authors may find useful:
- The "pressure" subsection seems a bit out of context. I understand the importance of computing the thermodynamics right, but it's not clear whether the relation they present is already included in their HOOMD package or if it's there just for reference in case someone needs to compute the pressure in their simulations.
- Fig. 5 is a bit weird, as it looks like the authors simulated very small systems and used periodic-boundary conditions to draw larger images. If this is the case then it should be explicitly stated or (even better) the simulations should be re-run with a larger number of particles. Moreover, looking at fig. 4 it looks like even with the swap the system has not equilibrated yet at the end of the simulation. I think that the authors should pick a system which can get to the ground state (i.e. ~100% bonds formed) with the swap mechanism.
- Fig. 3 has a similar issue in that the three curves are a bit noisy. If it is not too much work I'd ask the authors to run longer simulations to clean up the data.
- Even though the authors explicitly write that the 3-body algorithm is not their own, there are a few places scattered in the paper where a distracted reader may get the wrong impression (see for instance the first sentence of IV.B or in V., "the most important feature of our method")
- I think Fig. 7 may very well be included in the main text, since it is already introduced therein.
-It would be nice to understand what is the impact of the 3-body part on the performance of the code. The authors state that using fewer beads interacting with this 3-body potential makes the code go faster, but some hard numbers would help an interested reader.

1- Add details about pressure calculation in systems with 3-body swaps in HOOMD
2- Perform simulations of larger systems that get to the equilibrium state with the swap for Figg. 4 and 5
3- Make sure that it is clear everywhere in the text that the authors are presenting an implementation of an already-existing algorithm
4- (optional) Improve the quality of the curves in Fig. 3
5- (optional) Add Fig. 7 to the main text
6- (optional) add some details/ballpark numbers about the impact of using a certain fraction of "swapping beads" in a simulation

- Versatile algorithm
- Highly reproducible

Algorithms for bond swapping have been recently introduced in both Monte Carlo (MC) and Molecular Dynamics (MD) simulations: they allow, e.g., for simulating self-healing and restructuring properties of network forming materials as well as for equilibrating associative systems faster. While in MC simulations a bond swap was implemented as a stochastic event, a possible -- quite efficient --implementation for MD simulations relies on the introduction of a three-body potential, which compensates the energy gain associated to the formation of a double bond.

This contribution presents (and provides in an official repository) an implementation of the three-body potential for swappable bonds into the framework of the HOOMD-blue toolkit.Two examples are provided as bed tests: 1. the same binary mixture of mutually attractive spheres used in the reference publication by F. Sciortino (reference 19 of this paper): the reproduction of the reference behavior is fully achieved; 2. a binary mixture of di- and tri-functional patchy colloids: the swap algorithm is proven to achieve an effective annealing towards to the ground state, thus guaranteeing a faster equilibration, compared to self-assembly based on simple pair potentials.

The contribution is significantly useful for many purposes within the area of soft matter simulations and I support its publication. Here below I drop some observations to possibly improve the quality of the manuscript.

- As the aim of part 2 is to demonstrate numerical efficiency, estimates for three-dimensional systems would be interesting.
- How is a bond between two patches defined? The gaussian attraction between them has tails...
- Even though the proof of principle is in place, it would be interesting to see the full equilibration curve of the patchy colloid system --until the fully bonded state is reached (it should be doable as it is a system of 120 particles in two-dimensions)
- Why is the box shape rhombic in part 1?
- From the brown and purple curves in Figure 4, it seems the systems freeze almost instantaneously into the network state, is that so?
- I naively assumed the 2body-interactions between patches would be described by (1) plus (3), while the 3body-interaction would be described by (2). But I am not completely sure about it because of the dangling sentence "Each particle is a WCA repulsive disk with two or three attractive patches represented by either a three-body potential like eqns (1,2) if swaps are used", where "either" suggests a missing "or". What is the two-body potential?
- The reference to the use of the RevCross is only for insiders...
- Error in the text (while it is correct in the caption of figure 6): lighter (green) and heavier (orange) A particles.
- The difference between the networks shown in figure 5 is a bit qualitative and relies on a trained eye.