How should I modify my molecular dynamics simulation so that known chemical reactions can be modelled?
Since you've written your own code, you will probably be able to judge this best, but I'll give you some pointers as to where I'd start.

As far as I've gathered, you simulate a NVT ensemble of a pair of chemicals.

First of all you have to decide whether you want explicit representation of the solvent or not. Hint: The latter is computationally far more expensive and I doubt that you have implemented all the optimizations that you'd find in larger MD software.

This is because this choice can drastically influence the diffusion properties of your chemicals. In a coarse-grained force field they might still be able to diffuse around pretty much unhindered, while when using explicit atomistic representation of the solvent, the properties may change to a larger degree.

The way I'd do it would be to check whether the two areas of the molecules that react with each other are close to one another, meaning that some distance criteria $(d&lt;d_\text{rxn.})$ has to apply. When this criteria is met, you can make the reaction happen by creating new bonds and changing the atom types that are bonded on-the-fly. It wouldn't be so tricky to implement this part, I guess. Caveat: This will only work if you have a reaction whose mechanism is close to the proposed scheme. You might introduce nontrivial errors otherwise.

Then you can model the "real" concentration curves of all the substances, and simultaneously change your system to reflect the change in number of molecules. I wouldn't be surprised that you'd find that the pressure would slowly decrease as the reaction takes place (since you keep the volume constant).