The isomerisation energies and thermal stabilities of the entire set of C40 fullerene isomers are investigated using molecular dynamics simulations with the recently developed machine-learning-based Gaussian Approximation Potential (GAP-20) force field. Our simulations predict high statistical correlation between the relative isomerisation energy distribution of all C40 isomers and their thermal fragmentation temperatures. The most energetically and thermally stable isomer is the symmetric C40-D2 isomer and its fragmentation temperature is 3875 ± 25 K. In contrast, the least stable isomer is the D5d nanorod-shaped isomer which is 146.8 kcal mol−1 higher in energy and decomposes at 2975 ± 25 K, i.e., around 900 K below the most stable isomer. In addition, we show that most isomers lie in the range between 20 and 60 kcal mol−1 above the most stable isomer, and nearly 60% of isomers decompose in the temperature range of 3425–3675 K. These computational results provide valuable insights into the synthesis of C40 fullerenes at high-temperatures.