Associative bond swaps in molecular dynamics

Thank you for sending us the reports on our manuscript "Associative bond swaps in molecular dynamics". We appreciate your efforts to provide an in-depth review of our work.
All the three Reviewers expressed positive comments on the validity, significance, originality, clarity, formatting and grammar, as well as insightful comments to further improve the manuscript.
Following the Reviewers’ suggestions, we have performed a series of additional calculations and have revised the manuscript accordingly.
Below we address all of the Reviewer’s comments and questions, including the corresponding changes in the manuscript. We believe these changes have significantly improved the quality of our manuscript and we hope that with this revision our work will be suitable for publication in Scipost Physics.

*Report 1:

(1.1) As the aim of part 2 is to demonstrate numerical efficiency, estimates for three-dimensional systems would be interesting.
- While the part about network self assembly is in 2d, the rest of the paper is about 3d dumbbells, so Figs. 3,6-8 already refer to 3d systems.

(1.2) How is a bond between two patches defined? The gaussian attraction between them has tails...
- That is true. We now explicitly say in page 4 that: We define two patches bonded if they are closer than d=2 sigma_p=0.2 sigma_WCA.

(1.3) 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)
- Following the suggestion of Ref.2, we decided to rerun a larger self-assembly system (from N=120 to N=1000). We stop the simulations when the swapping system forms more than 90\% of the bonds, which corresponds to ~5days of cpu time. The new figure 4 gives a much more complete picture of the equilibration process.

(1.3) Why is the box shape rhombic in part 1?
- There was no need for the rhombic shape and it was there for historic reasons. We decided to re-run new simulations in square boxes.

(1.4) From the brown and purple curves in Figure 4, it seems the systems freeze almost instantaneously into the network state, is that so?
- Yes the referee is correct. While the small size of the system in the old version of Fig.4 exacerbated the difference between strong and weaker bonds, it is still evident that most of the bonds of the non-swapping mixtures appear almost instantaneously on our timescale.
This strengthen our conclusion that a mechanism to rearrange bonds is extremely valuable in many circumstances.

(1.5) 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?
- We edited the sentence to clarify that the gaussian attraction we mention is the interaction that we use when modelling non-swappable systems.

(1.6) The reference to the use of the RevCross is only for insiders...
- The three-body potential discussed in this paper can be directly used from HOOMD-blue with the name 'RevCross', so we decided to mention this fact to make the paper more pedagogical for possible HOOMD-blue users.

(1.7) Error in the text (while it is correct in the caption of figure 6): lighter (green) and heavier (orange) A particles.
- Thank you for pointing this out. We fixed this.

(1.8) The difference between the networks shown in figure 5 is a bit qualitative and relies on a trained eye.
- We now highlight the differences a bit more in the caption where we say that: [for strong bonds] There are in fact long sequences of bonded red beads that should not be present at equilibrium.

*Report 2:

(2.1) Add details about pressure calculation in systems with 3-body swaps in HOOMD
- We now explicitly say that the 3-body pressure is automatically calculated by HOOMD, in the same way it can be done for any other potential.
For more details about the calculation we say that: " In the Supplemental Material of Ref. [20] we show how to derive them from the standard virial approach"

(2.2) Perform simulations of larger systems that get to the equilibrium state with the swap for Figg. 4 and 5
- Following the suggestion of the Referee we run simulations of N=1000 units to update Fig.4-5.

(2.3) Make sure that it is clear everywhere in the text that the authors are presenting an implementation of an already-existing algorithm
- We reworked several statements to make this point more clear.

(2.4) (optional) Improve the quality of the curves in Fig. 3
- We think the improvement in the figure that might be achievable would not be worth the effort required to realise it. Since the referee marked this suggestions as optional, we decided not to pursue it.

(2.5) (optional) Add Fig. 7 to the main text
- We prefer to keep this figure in the supplement.

(2.6) (optional) add some details/ballpark numbers about the impact of using a certain fraction of "swapping beads" in a simulation
- Thank you for this very nice suggestion. We added a sec.IV C to specifically discuss this.

*Report 3:

(3.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?
- We agree with this point since studying vitrimers was also our original motivation for developing the method. We added sec.IV D to discuss more details about how vitrimers can be modelled.

(3.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.
- In sec.IV D we also discuss about the MD timescale. In particular, we mention that the timescale does not depend on the three-body potential, but rather on the coarse-grain parameters that would also need to be chosen in a simulation that does not use the three-body approach.

(3.3) The reference cited in Figure 3 and its caption does not match.
- Fixed. Thank you.