As mentioned that the classical molecular dynamics (MD) simulations could analyze the thermal stability and structural transformation of qHPC 60 , qHPC 58 B 1 N 1 and qHPC 54 B 3 N 3 , we plot the total-energy at various temperatures in Figure (a). The curve for the qHPC 60 (black line), is lower in energy than qHPC 58 B 1 N 1 and qHPC 54 B 3 N 3 indicating its greater stability. The higher stability of pristine qHPC 60 is also supported by our first-principles simulation results predicting its notably higher negative formation energy value. However, the qHPC 60 , qHPC 58 B 1 N 1 and qHPC 54 B 3 N 3 show similar trend for the drop in their energies. Each curve in Figure (a) can be divided into three regions based on the slope: before the jump in the total energy (I), the jump itself (II) and after the jump (III). Region I corresponds to the phase where the qHPC 60 cages are a bit structurally deformed as the temperature rises, as shown in Figure ), S2(a-b) and S3(ab). In this wide range, it contains large numbers of local configurations of pentagons, hexagons and heptagons, etc., which keeps the overall number of C-C bonds unchanged. As a result, the system's total energy does not increase sharply with temperature. The region where larger rings were lost and linear atomic chains (LACs) formed before the lattice atomization started occurring in Region (II). The qHPC 60 system begins to lose structural integrity, signalling the onset of the melting process at around 3900K. This causes the cage to rupture and unravel, as shown in Figures ), S2(c) and S3(c). At the final stage, Region III with T > 3900 K [See Figs. S1(d)-S2(d)], the phase transition occurs corresponding to the gaseous phase of the system where only atoms/short atomic chains are observed. The system looks uniform, as it should be in gas, leading to a sharp energy increase.