Because non-reactive molecular force fields cannot account for bond breaking and forming, the simulation of the curing reaction usually involves an iteration of three steps (embedded into classical molecular dynamics): the first being a reactive-groups search, followed by a topology change, and the structural relaxation of the reacted system. The identification of reactive pairs evokes a change in connectivity and in the molecular structure of the system, which in turn requires an update of the molecular mechanics parameters defining the molecular structure after the reaction. This is necessary to continue with the time propagation of the system after (local) structural relaxation. While the formulation of these steps is in principle universal, details of their implementation may vary significantly. For example, the change of connectivity is typically based on finding pre-defined binding partners within a cutoff distance. The latter is often increased over time to facilitate further searching of reacting atoms in a progressively less diffusive environment and curing. While typical search radii of 4 -10 Å are applied, values up to 14 Å were used to achieve higher curing rates. Another challenge with non-reactive MD approaches is the parametrization of the reactants and molecular entities appearing in the system along the curing process. The drastic changes in the connectivity and the transformed ("reacted") chemical groups are associated with completely different electrostatic interactions, which need to be accounted for. An accurate updating of the partial atomic charges during the curing process, and hence a correct calculation of electrostatic interactions, plays a key role in the bonding dynamics of the polymerizing system. It is thus surprising that many computational studies investigating the polymerization reaction of epoxy resins show a lack of methodological detail regarding the derivation procedure of partial atomic charges, as noted recently by Demir and Walsh. Indeed, partial charges are omitted, not mentioned or insufficient information is presented to allow for reproducibility. Besides determining appropriate partial charges, maintaining charge neutrality of the system is an important constraint that should apply to all procedures modeling the curing reaction. This demands partial atomic charges be updated after each reaction event to keep a net sum of zero, while also reflecting the post-reaction electrostatic environment. To this end, self-consistent schemes based on equalization of chemical potentials (QEq) , and on electronegativity equalization (EEM) are available. These methods show accuracy and transferability, and have been applied to polymerization, but remain computationally demanding, especially in the case of large systems simulating numerous rounds of cross-linking. Therefore, streamlining the curing process would be ideal, but this has not been achieved so far. The result of the curing simulation is a reliable molecular scale model of an amorphous epoxy resin which allows to relate thermodynamic and mechanical properties to the molecular structure of the resin. However, all these properties depend on the thermodynamic conditions, and their extraction from the simulation data often involves approximate procedures. Further difficulty is the required sampling and large system sizes that capture the complex topology of the macroscopic network. Therefore, obtaining meaningful mesoscopic data comparable to experiments is another demanding task, which, however, needs to be completed to advance the field.