Recently, it was demonstrated that the conductivity of wurzite (wz) ZnO bicrystal samples containing ${0001}$ inversion domain boundaries (IDB) can be massively and reversibly tuned by mechanical loading. As a step towards a detailed microscopic understanding of this effect, we systematically investigate the atomic structure and chemical composition of such IDBs using density functional theory calculations. In total, 92 model geometries that differ in structure and/or chemical composition are constructed, optimized, and compared thermodynamically. The lack of higher symmetries in wz ZnO prohibits a straightforward calculation of individual grain boundary (GB) excess energies. However, we show that, in nonperiodic slab models of wz ${0001}$ IDBs, the additional surface contribution to the total energy may be approximated by that of corresponding zincblende (zb) surfaces; the latter can be obtained by a series of prism calculations. Subtracting these surface energies allows us to construct absolute GB energy diagrams for wz IDBs and compare their thermodynamic stability with other GBs known from the literature. We find that thermodynamically favored IDBs are characterized by fully (4-fold) coordinated atoms and possess relatively low excess energies that range from 45 to $95\phantom{\rule{4pt}{0ex}}\mathrm{meV}/{\AA{}}^{2}$, depending on the termination (Zn/Zn or O/O) of the IDB and the exchange-correlation functional used in the calculation (LDA, GGA, or $\text{GGA}+U$). The electronic properties of the GB deviate only weakly from those of the bulk and are rather insensitive towards compressive and tensile strains. Our results thus indicate that experimentally observed piezotronic properties of wz bicrystals are not an intrinsic property of the pristine GB itself, but originate, for example, from externally supplied trapped charges, defects, impurities, or dopants. Low-energy structure models identified here may also be transferable to other wz- or zb-type IDBs (e.g., GaN, AlN, SiC, etc.).