The thermodynamic binding networks (TBN) model was recently developed as a tool for studying engineered molecular systems. The TBN model allows one to reason about their behavior through a simplified abstraction that ignores details about molecular composition, focusing on two key determinants of a system's energetics common to any chemical substrate: how many molecular bonds are formed, and how many separate complexes exist in the system. We formulate as an integer program the NP-hard problem of computing stable configurations of a TBN (a.k.a., minimum energy: those that maximize the number of bonds and complexes). We provide open-source software that solves these formulations, and give empirical evidence that this approach enables dramatically faster computation of TBN stable configurations than previous approaches based on SAT solvers. Our setup can also reason about TBNs in which some molecules have unbounded counts. These improvements in turn allow us to efficiently automate verification of desired properties of practical TBNs. Finally, we show that the TBN's Graver basis (a kind of certificate of optimality in integer programming) has a natural interpretation as the "fundamental components" out of which locally minimal energy configurations are composed. This characterization helps verify correctness of not only stable configurations, but entire "kinetic pathways" in a TBN.