## 3 Computational Methods  
We used two computational packages in this study: VASP [26] and Orthogonal Linear Combination of Atomic Orbitals (OLCAO) [27]. Both are based on the density functional theory. For structural optimization and elastic properties calculations using VASP, we used the PAW-PBE potential [28] with generalized gradient approximation (GGA) for the exchange correlation potential. We used a relatively high energy cutoff of 600 eV. The electronic and ionic force convergence criteria are set at 10\({}^{-9}\) eV and 10\({}^{-7}\) eV/A respectively. A 6 \(\times\) 6 \(\times\) 6 k-point mesh was used for C1, C2, and C3 crystals and a less dense mesh of 4 \(\times\) 4 \(\times\) 4 for the two larger crystals C4 and C5. No discernable difference in the relaxed structure was observed using other exchange correlation functions such as hybrid functional PBE0, HSE03, and applying van der Waals correction in the case of K\({}_{2}\)Mg(H\({}_{2}\)P\({}_{2}\)O\({}_{7}\))\({}_{2}\).2H\({}_{2}\)O. The calculation of mechanical properties using the optimized structure has been described before [29; 30]. Specifically, we apply a small strain \(\epsilon\) (\(\pm\) 0.25%) to the crystal. The elastic tensor elements are obtained using the stress (\(\sigma_{\rm j}\)) vs strain (\(\epsilon_{\rm i}\)) response analysis scheme to the relaxed structure and then obtained the elastic coefficients \(C_{\rm ij}\) (i, j = 1, 2, 3, 4, 5, 6 (using Hooks law) and compliance coefficient \(S_{\rm ij}\) by solving the set of linear equations From the calculated \(C_{\rm ij}\) and \(S_{\rm ij}\), other mechanical properties such as bulk modulus (K), shear modulus (G), [PERSON]'s ratio (\(\eta\)), and Young's modulus (E) are obtained using Voight-Reuss-Hill (VRH) polycrystals approximation [31; 32]. For electronic structure and interatomic bonding, we use the in-house developed OLCAO with the VASP-relaxed structure as input. A more localized minimal basis (MB) is used for the calculation of effective charge \(Q_{\alpha}^{*}\) and bond order (BO) values. The BO is the overlap population \(\rho_{\alpha\beta}\) between any pair of atoms (\(\alpha\), \(\beta\)) based on Mulliken population analysis [33; 34]. Mulliken analysis is only meaningful for comparisons between different materials if the same basis set in the same computational package is used as in the present case. Mulliken analysis is more effective if a localized basis is used, and we always use the minimal basis as implemented in OLCAO method with the same atomic basis for all atomic species [27]. The same approach has been successfully applied to many different materials systems including complex biomolecular systems [35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45]. \[Q_{\alpha}^{*}=\sum_{i}\sum_{m,occ}\sum_{j,\beta}C_{i\alpha}^{*m}C_{\beta \beta}^{m}S_{i\alpha,j\beta} \tag{1}\]  
\[\rho_{\alpha\beta}=\sum_{m,occ}\sum_{i,j}C_{i\alpha}^{*m}C_{\beta\beta}^{m}S_{ i\alpha,j\beta} \tag{2}\]  
In the above equation, \(S_{i\alpha,j\beta}\) are the overlap integrals between the \(i^{th}\) orbital in \(\alpha^{th}\) atom in the \(j^{th}\) orbital in \(\beta^{th}\) atom. \(C^{m}_{\beta\beta}\) is the eigenvector coefficients of the \(m^{th}\) occupied band. The BO from equation (2) defines the relative strength of the bond. Total bond order density (TBOD) is obtained by normalizing the total BO (TBO) with cell volume. TBOD is a single metric to access the internal cohesion in the crystal and can be decomposed partial components (PBOD) for any structural units or groups of bonded atomic pairs.