Density functional theory (DFT) based Born-Oppenheimer molecular dynamics (BOMD) simulations were carried out using the CP2K 33 software suite. We have specifically chosen to use this software package since the implementation of linear-scaling Kohn-Sham approaches in CP2K allows robust and efficient electronic structure calculations for large systems. The Perdew-Burke-Ernzerhof (PBE) was used for the DFT calculations with Grimmes's D3 method to account for dispersion forces. To obtain reasonable accuracy, we utilized a DZVP (double zeta valence polarized) basis set for Al in the DFT calculations and the TZV2P basis set for C, O, H, and P atoms with the Goedecker, Teter, and Hutter (GTH) pseudopotentials for atomic core electrons. A similar basis set for Al atoms was also used to accurately calculate methane activation and alkane dehydrogenation 39 on the γ-Al 2 O 3 surface. The orbital transformation method with an electronic gradient tolerance value of 1 × 10 -5 atomic units was adopted as the convergence criteria for the SCF cycle. A kinetic energy cut-off of 400 Ry for the auxiliary plane-wave basis with a 0.5 fs timestep was employed to integrate the equations of motion. The initial guess was furnished by the stable predictor-corrector extrapolation method at each molecular dynamics step. We have carried out several molecular dynamics simulations in the 200 -1000 • C temperature range in increments of 100 • C (i.e., a total of nine independent trajectories were considered).