*4.2. Density Functional Theory Calculations*

In all of the *ab initio* calculations, we have used the solids-corrected Perdew-Burke-Ernzerhof (PBEsol) GGA exchange-correlation functional [29,30], and all structural optimisations were deemed converged when the atomic forces were less than 0.01 eV/Å.

A natural choice for the calculations on the double bubble clusters, due to its computational efficiency, is the DFT code FHI-aims [31] as it uses numeric atom-centred basis sets. These calculations were performed with the species defaults for the "tight" basis sets for accuracy (energies converged to 1 meV/atom) and with scalar ZORA relativistic treatment [32]. We have used the plane-wave DFT code VASP [33–36] to determined the equilibrium structures of the double bubble based framework (extended crystal systems—see Section 2), and, for comparison, wurtzite bulk ZnO and GaN. Within VASP, we employed the projector augmented wave (PAW) method [37] to describe the interactions between the cores (Zn:[Ar], Ga:[Ar], O:[He] and N:[He]) and the valence electrons.

To determine the equilibrium bulk structures avoiding the problem of Pulay stress, we have optimised the atomic coordinates at a series of different volumes, and fitted the resulting energy *versus* volume data to the Murnaghan equation of state.

We have found that for the framework systems, an energy cut-off of 500 eV, and Monkhorst-Pack *k*-point meshes of 8 × 8 × 6 and 1 × 1 × 1 for, respectively, the pure bulk wurtzite systems, and the (A)12@(B)48 systems, where A and B stand for either ZnO or GaN, provide convergence in total energy up to 10<sup>í</sup><sup>5</sup> eV for the framework systems, which is comparable with our double bubble cluster calculations.
