In this study, we built up a comprehensive computational dataset for C 20 -C 60 fullerenes with 5770 structures and computed 12 fundamental properties at DFT-level. We conducted statistical analysis on stability, electronic properties, and solubility, and assessed the Pearson correlation coefficients among them. The observed weak correlations between E g and both E b (r = 0.13) and logP (r = -0.24) suggest that favorable electronic properties can be independently adjusted for particular applications without compromising the stability and solubility. To gain a deeper understanding of the structure-property relationships in fullerenes, we introduced various multi-dimensional topological features, as well as geometric measures to further explore the nature and origin of these properties expanding beyond the commonly used IPR. Among them, atom, bond, and hexagon features prove Table : The binding energies (E b ) of four large-size fullerenes (Cn) calculated with B3LYP-D3/6-311G* compared to the values predicted by linear fitting models with atom, bond, pentagon, and hexagon features respectively across C 20 -C 60 fullerenes, together with various statistical measures, i.e., the mean deviation (MD), the overall mean absolute deviation (MAD), the standard deviation (SD), and the maximal absolute error (MAX). The units for all the values are given in eV/atom. The molecular structures of these large-size fullerenes are shown in Figure .