In this quantum chemical study, we explore hydrogen bonding (H-bonding) and stacking interactions in different crystalline cellulose allomorphs; namely, cellulose I(β) and cellulose III(I). We consider a model system representing a cellulose crystalline core made from six cellobiose units arranged in three layers with two chains per layer. We calculate the contributions of intrasheet and intersheet interactions to the structure and stability in both cellulose I(β) and cellulose III(I) crystalline cores. Reference structures for this study were generated from molecular dynamics simulations of water-solvated cellulose I(β) and III(I) fibrils. A systematic analysis of various conformations describing different mutual orientations of cellobiose units is performed using the hybrid density functional theory with the M06-2X with 6-31+G(d,p) basis sets. We dissect the nature of the forces that stabilize the cellulose I(β) and cellulose III(I) crystalline cores and quantify the relative strength of H-bonding and stacking interactions. Our calculations demonstrate that individual H-bonding interactions are stronger in cellulose I(β) than in cellulose III(I); however, the total H-bonding contribution to stabilization is larger in cellulose III(I) because of the highly cooperative nature of the H-bonding network. In addition, we observe a significant contribution from cooperative stacking interactions to the stabilization of cellulose I(β). The theory of atoms-in-molecules (AIM) has been employed to characterize and quantify these intermolecular interactions. AIM analyses highlight the role of nonconventional CH···O H-bonding in the cellulose assemblies. Finally, we calculate molecular electrostatic potential maps for the cellulose allomorphs that capture the differences in chemical reactivity of the systems considered in our study.