The Empirical Valence Bond (EVB) method offers a suitable framework to obtain reactive potentials through the coupling of non-reactive force fields. However, most of the implemented functional forms for the coupling terms depend on complex spatial coordinates, which precludes the computation of the stress tensor for condensed phase systems and prevents the possibility to carry out EVB molecular dynamics in the isothermal-isobaric (NPT) ensemble. In this work, we make use of coupling terms that depend on the energy gaps, defined as the energy differences between the participating non-reactive force fields, and derive an expression for the EVB stress tensor suitable for computations. Implementation of this new methodology is tested for a model of a single reactive malonaldehyde solvated in non-reactive water. Computed densities and classical probability distributions in the NPT ensemble reveals a negligible role of the reactive potential in the limit of low concentrated solutions, thus corroborating the validity of standard approximations customarily adopted for EVB simulations.