Brute force sampling would require DFT transition state optimizations and vibrational frequency analysis on O(10 5 )-O ) sites, which is not feasible. Therefore, most of the effort is devoted to qualitative studies. The first requirement is generating a reasonable model of functionalized amorphous silica, which is typically accomplished with classical molecular dynamics (MD) simulations. The simulations start from a periodic crystalline form of silica, which is heated up to a liquid form, then sliced and cooled down into a glassy amorphous structure. DFT simulations are then conducted with either periodic boundary conditions (PBC), or by carving out molecular cluster models. PBC simulations preserve the solid nature of the material, but they require large simulation cells to model disorder and they are computationally intensive. Cluster models are computationally cheaper, but one needs to minimize computational artifacts due to the finite cluster size. We have shown that large cluster models (about 300 atoms in size) can accurately reproduce a variety of experimental characterization techniques. However, small size models are really necessary to conduct efficient quantitative studies of reactivity. The goal of this work is to address a variety of issues related to simulations of reactivity on metal-doped amorphous silicates (although the same issues apply to most amorphous supports) that hinder our ability to evaluate effective reaction rates directly comparable to experiment. The first issue is whether small cluster models of the active site are really reliable representations of large models with extended silica framework. The second issue is whether taking a single transition state structure per active site is really representative of the activity of that site. The third issue is to determine what features (i.e., descriptors) of an active site correlate with its activity, and what makes one site more active than another. We use the epoxida-tion reaction of ethylene catalyzed by Nb-doped silica as a test case as it has been well characterized experimentally and theoretically, so the reaction mechanism is known. We compare the adsorption energy for H 2 O 2 and C 2 H 4 on the metal site and the activation energy for the reaction between large cluster models with approximately 300 atoms and minimalistic models with a single layer of Si atoms around the metal. Afterwards, we show that multiple TS structures can be found for each site that differ significantly in adsorption and activation energies. Finally, we use a machine learning (ML) model based on non-linear regression to correlate the activation energy with relevant descriptors. Although the data set is too small for quantitative prediction of k ef f , this analysis qualitatively reveals what makes a catalytic site more active than another for this test reaction.