Similar to our previous work on the metal doped amorphous silica library (METASIL), the clusters were generated from amorphous silica slabs obtained with classical molecular dynamics (MD) simulations. However, in this work the slabs were taken from a recent paper by Thompson and co-workers, specifically, those with 1K/ps cooling rate, 2 OH/nm 2 density, and the elongated-bond functionalization procedure. The detailed procedure is described in Ref. 25, but here we summarize the key points. The MD simulations were carried out with the LAMMPS software and the BKS force field was used to describe the interactions within the silica. The simulations started from a crystalline β-cristobalite simulation cell, which is composed of 576 atoms. The crystal was heated to 8000 K for 1 ns to create molten silica; the simulation cell was then expanded in the z-dimension to form the liquid silica slab. After equilibration for 3 ns at high temperature, the system was cooled down at a 1K/ps cooling rate to 298 K. At this point, the amorphous silica slab was allowed to equilibrate for 1.5 ns. Both sides of a slab (orthogonal to the z-axis) were functionalized with silanol groups by breaking elongated siloxane bonds at the surface, and capping the other broken bond with a H atom (i.e., a water molecule was added to the system for each new silanol group, as in the hydrolysis reaction). The functionalized slabs were relaxed and equilibrated for 1 ns at 298K. This protocol was not meant to precisely mimic a specific synthetic route for the production and functionalization of silica surfaces, as the experimental procedures occur at time scales that are not amenable to computer simulations. Nevertheless, simulations of characterization techniques using the functionalized slabs were in agreement with available experimental data, which indicates that the slabs are reliable. The clusters were carved out from the MD structures as described in Ref. 31. However, in this paper we formed one cluster from each slab. The radius of each cluster bowl was approximately 12 Å and the dangling Si-O bonds were capped with hydrogen, so that each model contains approximately 275 non-H atoms on average. These clusters are larger than those in METASIL, and they were relaxed using the procedure we developed in Ref. 35 to minimize the so-called procedural disorder. The latter is an artifact of cluster simulations, where the boundary between the MD and QM region may create an artificial straining of the active site that can affect reaction energies. We showed that procedural disorder can be minimized by pre-optimizing the entire cluster (except for the capping atoms) at the target QM level of theory before doping the cluster with a metal center. After doping, only a smaller part of the cluster needs to be relaxed to adjust to the presence of the metal. Doping was performed by replacing a central Si -OH group with a Nb --O (Nboxo) group in a distorted-tetrahedral structure (Nb V ), as suggested by experimental data. After doping with a Nb site, the region within 9 Å from the metal was relaxed again. This relaxed zone is larger than the 8 Å radius recommended in Ref. 35, therefore, we do not expect any sizable procedural disorder effects in this work. We built fifteen large clusters, which will be called "Si-12Å" in the rest of the paper, see an example in Figure . The full geometry of the Si-12Å clusters is reported in the supporting information (SI).