This is the accepted version of the journal article: González, Danilo; Camino, Bruno; Heras-Domingo, Javier; [et al.]. «BCN-M : a free computational tool for generating Wulff-like nanoparticle models with controlled stoichiometry». Journal of physical chemistry C, Vol. 124, issue 1 (Jan. 2020), p. 1227-1237. DOI 10.1021/acs.jpcc.9b10506 This version is available at https://ddd.uab.cat/record/279250 under the terms of the license BCN-M: A Free Computational tool for Generating Wulff-like Nanoparticle Models with Controlled Stoichiometry Danilo González, Bruno Camino, Javier Heras-Domingo, Albert Rimola, Luis Rodríguez-Santiago,* Xavier Solans-Monfort,* Mariona Sodupe Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain 1 Abstract A thorough knowledge of the atomic structure of nanomaterials is of high importance to understand their properties. This requires developing nanoparticle models, which is not always straightforward, particularly in the case of non-pure metallic systems. The Bulk Cut Nanoparticle Models (BCN-M) computational tool generates Wulff-like models for binary materials with controlled stoichiometry automatically with none or little need for further manipulation from the user. The models are obtained exclusively by introducing the structure of the bulk material, its symmetry, the surface energies of the most representative surfaces and information about surface termination as input data. The algorithm produces different structural model sets and the quality of these models is evaluated using different criteria: i) the deviation from ideal Wulff shape; ii) the global coordination of surface metal atoms and iii) the polarity of the model. BCN-M has been applied to fifteen different materials leading to a variety of models that cover the most relevant binary ionic structures and symmetries (cubic, tetragonal, hexagonal and monoclinic). The resulting models can be used for structure analysis of ideal systems as well as their simulation. BCN-M is available as a free web platform (https://bcnm.qf.uab.cat) or as a downloadable utility and it is expected to be an important tool for the design of future nanomaterials. 2 Introduction The interest of the scientific community in nanochemistry has exponentially grown over the last decade.1–4 One of the key features of nanomaterials is exhibiting properties that significantly differ from those of the bulk material. A variety of nanostructured materials such as graphene-based sheets (2D), carbon nanotubes (1D), metal and semiconductor nanoparticles (0D), within others, have been synthesized and engineered in order to have potential applications in several fields as wide as optoelectronics,5–7 biomedicine8– 12 and catalysis.13–21 In the last years, great progress has been done in terms of characterization of nanostructured materials by means of experimental techniques, even at an atomistic level.22,23 This characterization allows establishing relationships between the structure of the nanomaterials and their chemical properties. As a powerful complementary tool, computational chemistry is usually used for stablishing these relationships. Indeed, combination of both computational modeling and experimental data has been demonstrated to be a successful strategy to obtain detailed and accurate information of nanoparticles (NPs) and their properties.24–26 Unfortunately, computational modeling of NPs is not trivial, mainly because of the large variety of NP sizes and shapes that can be synthesized for any particular material. There are two main strategies for designing realistic structural NP models: i) the bottom-up27– 30 and ii) the top-down approaches.29,31–34 The bottom-up procedure builds up NPs by joining atomic/molecular/cluster sized entities and it is useful when one aims to work in the lowest limits of the nanoworld. In this approach the most stable atomic structure is not constrained to be bulk-like and accordingly many combinatorial possibilities can be 3 low-energy structures. Thus, a confident theoretical study must unravel the rich energy landscape including all the NP isomers. This is commonly performed with global optimization algorithms.35–39 In contrast, the top-down approach starts from the bulk structure. It does not explore all possible isomers but it is suitable for exploring larger nanosized materials.29,31–34 The NP model is obtained by reducing the size from the bulk material to the nanoscale. For monocrystalline systems, the size reduction results in NP morphologies that are linked to the thermodynamic equilibrium crystal shape of the bulk system, which depends on the energy cost of forming the different surfaces, edges and vertices as well as strain effects. In this context, in absence of strain effects, the nanoparticle is large enough and monocrystalline, the edge and vertices effects are negligible and the usual way to define the crystal morphology is by using the Wulff construction, which is based on the Gibbs- Wulff theorem.34,40–45 The Gibbs-Wulff theorem postulates that under thermodynamic equilibrium, a crystal acquires the structure that satisfies the requirement of minimum total surface energy. Thus, the surface contribution can be estimated by a vector perpendicular to the surface, (i) from the center of a crystallite, whose length (l) is proportional to its surface free energy (γi): 𝑙! = 𝑐! 𝛾! Accordingly, the surfaces presenting the shortest vectors are predominant. The Wulff construction results in a polyhedron that depends only on the ratios between surface energies and symmetry point group.46,47 Therefore, the higher the number of crystallographic facets used to construct the Wulff shape, the higher the difficulty in building it. The range of validity of the NPs based on the Wulff construction method has 4 been largely discussed and it is usually accepted for NPs sizes larger than 2-3 nm depending on the material.29,36 There are several available codes capable of building Wulff construction plots from the surface energy ratios and the symmetry group of the material.48–50 Some tools allow obtaining atomistic models for single component materials. As an example, the suite of Python tools implemented in the Atomic Simulation Environment (ASE) allows generating Wulff shaped nanoparticle models of metals.51 Moreover, multicomponent nanoparticle models for some materials where stoichiometry is obtained without further manipulation can be generated with some other pieces of software.52–54 However, the procedure is not general for all binary systems, since in most cases, post- processing is mandatory to obtain stoichiometric and charge electroneutral models. The most usual post-processing operation is the hand cleaning, i.e., the removal of the excess atoms or the addition of missing ones to reach the desired structure. As a consequence, the success on obtaining a realistic model strongly depends on the user abilities and subjectivity and accordingly, the models can easily suffer from human bias and mistakes are likely to happen. At this point, it is worth mentioning that there are several other aspects to take into consideration. Two of them are: i) the presence of singly coordinated atoms in the final nanoparticle model, and ii) the contribution of polar surfaces. Singly coordinated atoms are sometimes mandatory to reach stoichiometric and electroneutral models. However, they are usually highly unstable and thus, at least for some particular cases, nanoparticle models derived from metal reduction or oxidation, by removal of singly coordinated species, are probably more realistic. This redox process maintains electroneutrality, but 5 the resulting model is not stoichiometric. On the other hand, polar surfaces are also highly unstable, and their formation is generally accompanied by an important surface reconstruction to reduce the polarity.55–58 Three main kinds of reconstruction are possible, the occurrence of which depends on the material and the synthetic conditions: i) generation of some metal and/or non-metal vacancies in a few of the outermost layers, which maintains the stoichiometry of the material, ii) redistribution of the electron density at the outermost layers due to redox processes, which leads to neutral non-stoichiometric surfaces (anion- or metal-rich), and iii) incorporation of external groups present in the media such as hydrides or hydroxyls. While the former is too specific of the material under study (hence requiring further manipulation by the user), the other two can be partially included in the automatic nanoparticle generation process. In this contribution, we present the computational user-friendly tool entitled Bulk Cut Nanoparticle Models (BCN-M), which is supported in the form of a free web platform (https://bcnm.qf.uab.cat) and as a downloadable utility. BCN-M is able to generate atomistic Wulff-like NPs with controlled stoichiometry in a systematic and reproducible way for a large variety of ionic binary compounds, including the most representative metal oxides and metal chalcogenides. In particular, we tested the performance of the BCN-M computational tool by constructing NP models of fifteen systems as test cases that cover the most common cubic (MgO (rock salt structure), CeO2 (fluorite structure), Li2O (antifluorite structure) and In2O3 (Bixbyite structure)), tetragonal (TiO2 anatase and MO2 rutile) and hexagonal (Wurtzite like metal chalcogenides) ionic systems, as well as materials involving either polar and non-polar surfaces in the Wulff construction. Several different kinds of models can be obtained depending on the material under study and 6 user specifications: i) stoichiometric and electroneutral nanoparticles; ii) nanoparticle models without singly coordinated atoms and an excess of the less abundant ion and iii) models involving polar surfaces that can either be metal or non-metal rich or saturated with H or OH groups. It is noteworthy that automation is intended to save time for the user and remove the bias introduced by the manual modification. However, we are aware that this cannot substitute the scientific rigor of the user who should carefully validate the models and have a thorough knowledge of the material under consideration. Methodology Computational details. The Gibbs surface energies required to build the Wulff construction are approximated to the total electronic surface energy arising from density functional theory (DFT) calculations, already reported in the literature.59–77 Moreover, in order to analyze the differences and similarities between the isomers of a given nanoparticle model, DFT calculations are carried out on the (RuO2)33 isomers. As in our previous works with RuO2,62,63 we use the PBE78 functional as implemented in the VASP code.79,80 Grimme (D2) empirical correction is added to account for dispersion interactions.81 Noticeably, PBE-D2 leads to lattice parameters that are in better agreement with experimental data than when using PBE-D3.62,82 The core electrons are described by ultrasoft pseudopotential and the external ones by plane wave basis set, with a kinetic energy cutoff of 500 eV. With the aim of minimizing replica interactions, a cubic cell of 25 Å is employed. The k point grid includes the G point only. Computational platform algorithm. BCN-M is a computational tool written in Python and uses standard and scientific Python libraries (NumPy,83 SciPy,84 ASE,51 pymatgen85 7 and PyYAML86). Its algorithm can be divided in four sections: i) initial nanoparticle generation, ii) nanoparticle symmetry and chemical evaluation, iii) nanoparticle refinement and iv) quality indices. The nanoparticle refinement includes the generation of the different types of nanoparticles (stoichiometric, without singly coordinated atoms and with contribution of polar surfaces). A description of each section is provided below. The flow chart of BCN-M is shown in Figure 1. Initial nanoparticle generation. The input data required for the generation of the final nanoparticle models (NPF) is: i) the irreducible atoms, types and positions in fractional coordinates; ii) the conventional cell parameters; iii) the space group number of the crystal (1-230); iv) the Miller indices of the planes constituting the NP and the related surface energies; v) the size of the desired NP in Å, defined with a central value, a range around this size and a step size; and vi) the type of surface termination that the final model should have. According to our experience, we suggest starting by exploring the NPs between 1.0 and 2.0 nm with a step size of 1.0 Å. This will result in at least one NP model. With this initial model, the user will be able to estimate the size of the subsequent realistic models and identify the most relevant center for the material. This would reduce the computational time specially when trying to construct large nanoparticle models. In relation to the nanoparticle center, the user can optionally set the center of the nanoparticle. However, by default, the tool explores three different types of nanoparticle centers: i) the positions of the irreducible atoms, ii) the origin of coordinates and iii) several special positions that present the lower symmetry multiplicity in the unit cell (See table S2 of the Supporting Information for the centers considered for each symmetry group). Within all tested examples, these centers guarantee obtaining at least one set of nanoparticle models. 8 Figure 1. BCN-M flow chart. AnBm material formula. NP0 initial nanoparticle. NPF Final nanoparticle model. The flow chart of the nanoparticle without dangling atoms and the nanoparticle with polar surfaces generation can be found in Figure S1 and S2. With this input data, BCN-M analyzes if any of the surfaces included in the input data is polar and then it constructs n initial NP models, which are called NPs0. The polarity of the surface is evaluated by the pymatgen library.87 If one polar surface is found, but the 9 user has not specified that there are polar surfaces, BCN-M warns the user. The NPs0 are built by replicating the unit cell along the planes selected by the user. The number of cells in this bulk expansion is estimated by the size of the nanoparticle and the surface energy in each direction. The final cut is performed by the proportional distance in each direction, adopting a shape as close as possible to the Wulff morphology. In order to achieve this goal, ASE is used to apply the symmetry operators to the irreducible atoms and to identify the directions that are equivalent by symmetry.88 Symmetry and chemical nanoparticle evaluation. The characterization and evaluation module is called to analyze the suitability of the generated NP0 through the following criteria: i) contribution of equivalent faces, ii) Wulff-like structure, and iii) absence of undercoordinated most abundant ions. If the NP0 does not fulfill one of these criteria the model is discarded, and other centers and sizes are explored. Next, a detailed description of each criterion is reported. Contribution of equivalent faces. Surfaces that are equivalent by symmetry have the same surface energy. Therefore, they should all present the same surface area in the final Wulff-like model. This is evaluated by determining the surface area from the length of the normal vector of each facet. The length of these vectors is computed as the largest distances between the center of the nanoparticle and the outermost atoms of the surfaces. This allows obtaining the ratio of each facet in the total surface area and evaluate the equivalence of the facets related by symmetry. BCN-M considers that equivalent facets contribute equally if the area difference is smaller than 5%. Wulff-like structure. The different families of surfaces in a nanoparticle model should contribute to the global shape with the same (or very similar) area ratio than the one 10 predicted by the Wulff construction. This is verified by computing the distance between the nanoparticle center and the equimolar Gibbs surface41,44 in each direction.89 With these distances, BCN-M computes the contribution of each family of surfaces to the nanoparticle model applying the Wulff-Gibbs theorem and compares it to that obtained with the usual Wulff construction procedure. If the percentage of presence of each individual facet does not follow the same order with the two methods, the model is discarded. Absence of undercoordinated less abundant centers. The presence of undercoordinated atoms is usually associated with less stable structures and thus they should be avoided in the construction of the NP models. As a consequence, the less abundant ions (e. g. Ti in TiO2 and O in Cu2O) are classified as “reasonably coordinated” and “undercoordinated” centers. The default threshold criterion chosen to identify the so- called undercoordinated atoms is two vacancies of the maximum coordination (max. coord – 2). This means, for example, that in a system where the maximum coordination for the less abundant atom is six (e.g., the case of TiO2), each center whose coordination is three or lower will be considered undercoordinated. The presence of a single undercoordinated center is sufficient for the tool to discard the model. It is worth mentioning that, for some materials, the default threshold is too strict and thus, no nanoparticle models are obtained. This is, for example, the case of MgO, where corners are tricoordinated and Cu2O, where the maximum coordination of oxygen is 2. For these and similar cases, the user can modify the undercoordinated atom threshold to half of the maximum coordination. 11 Nanoparticle refinement. The accepted NP0 models are Wulff-like and do not present undercoordinated more abundant centers. However, surface termination has not yet been analyzed. The nanoparticle refinement section is the responsible for achieving more realistic structures. Four potential terminations are possible: i) stoichiometric nanoparticles, ii) nanoparticles without singly coordinated atoms, iii) models involving polar surfaces that are either metal- or non-metal-rich, and iv) models involving polar surfaces whose polarity is compensated by the addition of H or OH groups. The suitability of each termination depends on the material and synthesis procedure, thus the user should have thorough knowledge on the material. These different possibilities are achieved through three different procedures, which are described below. Stoichiometric nanoparticles. The key issue in this procedure is to ensure that the material stoichiometry is respected in the final model. Therefore, for each generated NP0 the stoichiometry of the model is verified. Models presenting an excess of the less abundant ion in the general formula are discarded. For instance, all TiO2 models where the O/Ti ratio is smaller than 2 are discarded, while for Cu2O, the discarded models present less Cu+ cations than twice the number of O2- atoms. If the model is non- stoichiometric because of an excess of the most abundant element, the extra atoms are removed in the “Remove ion excess loop” (see supplementary information). In this loop, only dangling ions are eliminated. For systems with a 1:1 stoichiometry, atoms in excess in the NP0 are the ones that are removed, regardless of their nature. In this cleaning process, two scenarios can happen: either the number of ions to delete matches the number of dangling ions of the NP0, or only a fraction of the dangling ions must be removed in order to achieve stoichiometry. The first case will only lead to one 12 final NP, whereas the second one leads to many different models because several combinations exist for the removal of atoms. With the aim of obtaining models with larger global coordination numbers, the cleaning process starts eliminating those dangling atoms linked to the most coordinated sites. This is done iteratively until stoichiometry is achieved. If in any of these steps the number of sites exceeds the number of dangling atoms to remove, the program assumes a random procedure that is repeated 1000 times (see Supporting Information for further details). While this should be enough to obtain a set of models containing the most representative structures, the user can modify this parameter for a more exhaustive sampling. It should be noted that some of the models obtained may be equivalent by symmetry and thus, only one structure should be kept. For models below 2 nm, BCN Model recognizes these equivalent structures by using the Social Permutation Invariant Coordinates (SPRINT), which are a kind of topological coordinates that reflect the connectivity between atoms (see below). These coordinates are used in several other computational tools/simulations for the same purpose90,91 and are based on the principal eigenvalue and eigenvector of the adjacency matrix of chemical bonds.91 The major advantage of the SPRINT coordinates is that they are rotationally and translationally invariant. Overall, two models presenting the same SPRINT coordinates are equally connected and are considered equivalent, and therefore one of the two models is discarded. Nanoparticle models without singly coordinated ions. As already mentioned above, the generation of stoichiometric models frequently leads to singly coordinated atoms on the surface. However, these structures are likely unstable and thus, at least for some 13 particular cases, metal reduction or oxidation to avoid the presence of singly coordinated species are energetically more favorable. In this context, BCN-M is able to generate models without dangling atoms, if this is specified by the user. This is done by following the same described procedure for the generation of stoichiometric nanoparticles except the “Remove Ion Excess Loop”, since, in this case, BCN-M removes all singly coordinated ions, regardless of whether stoichiometry is preserved or not. Nanoparticles involving polar surfaces. Polar surfaces are usually highly energetic and only contribute to the Wulff construction in some particular cases. Due to the instability of polar surfaces, surface reconstruction usually takes place.55–58 BCN-M is able to include some common terminations of polar surfaces, i.e., metal-terminated, anion- terminated and addition of H or OH groups. The nanoparticle refinement for models including polar surfaces is called with the appropriate keyword. A NP0 can present one or more non-equivalent facets, in which at least one of these families of facets is polar. BCN-M computes the stoichiometry of the outermost plane and determines if it is metal- or anion-rich. If the ions in excess match with the required termination, BCN-M assumes that the NP0 polar surfaces do not require further refinements. In the case that the ions in excess do not match with the required termination, the NP0 is increased in several steps in this polar direction until the required termination is achieved, the maximum displacement being the largest interatomic distance between the different irreducible atoms. Once, the polar surface satisfies this requirement, BCN-M removes all dangling atoms of the remaining non-polar surfaces. Finally, if the user specifies having an OH terminated surface, BCN-M constructs first the metal rich nanoparticle and then adds the H or OH groups to the polar surface. 14 Quality indices. The characterization and evaluation module is again called to compute the Wulff-like index, the distribution coefficient and the global coordination quality indices to guide the users in selecting the most relevant model for their particular application. The distribution coefficient index is associated with the dipole moment of the generated NPs, while the global coordination is used to identify the most saturated models. These quality indices are briefly described below. Wulff-like quality index. This index is computed as the sum of the deviations of each surface family from the ideal Wulff polyhedron (see above), in terms of absolute differences between percentages. Distribution coefficient. The distribution coefficient (DC) index is used to roughly evaluate the dipole moment of the nanoparticle model and is calculated as the standard deviation of the center of mass of the metals (CMM) and that of the anions (CMA). If the distribution coefficient is 0, the anions are distributed homogeneously around cations of the NP. The rule of thumb is: the larger the DC, the more heterogeneous the ions distribution and thus, the larger the polarity of the model. At this point, it is worth mentioning that two models can have the same DC but can differ in their structure. Global coordination. Similarly to the absence of undercoordinated cations criterion, the global coordination index is used to distinguish between models with a different percentage of less coordinated ions. The reported value is the sum of all cation and anion coordination numbers. The final outputs are XYZ format files so that they can be opened with a wide variety of visualization programs. The comment line of the XYZ output includes the distribution 15 coefficient, the center of the nanoparticle model, the global coordination and the Wulff like criterion. Results and discussion In this section we present a series of NP models obtained with the BCN-M tool. The systems used as a test case cover the most representative binary ionic compounds with different symmetries. These include cubic, tetragonal, hexagonal and monoclinic crystallographic systems (Table 1). Moreover, they have been selected to explore the suitability of BCN-M in generating stoichiometric nanoparticles with and without dangling atoms, non-stoichiometric models without dangling atoms, and models involving polar surfaces. It is worth noting that the generated models arise from cutting the bulk without any further manipulations. Consequently, surface reconstruction is not taken into account explicitly. In practice, the resulting models are thought to be used as a starting point for further simulations or to analyze ideal structures. Indeed, similarly to what have been done with fully metallic systems,92 BCN-M models can be used to determine the number of different active sites (tips, corners, edges, surfaces, etc.) as a function of the nanoparticle size and thus, deduce the optimal ideal size for one particular catalytic process. The section is organized in different subsections: i) presentation of the obtained models, ii) Influence of the optional parameters in the generated models and iii) electronic structure analysis of the isomers obtained for (RuO2)33. Models. Fifteen different test examples have been considered to ensure the performance of the here described algorithm. Table S3 of the Supporting Information reports the parameters (fractional coordinates, cell parameters, space group, Miller indices, surface energies) to be introduced in BCN-M as input data. 16 Table 1. Test cases explored in this contribution Non-polar surfaces Material Space group MgO (rock salt) Cubic 𝐹𝑚3$𝑚 (225) CeO2 (Fluorite) Cubic 𝐹𝑚3$𝑚 (225) Li2O (anti-fluorite) Cubic 𝐹𝑚3$𝑚 (225) Cu2O Cubic 𝑃𝑚3$𝑛 (224) In2O3 (Bixbyite) Cubic 𝐼𝑎3$ (206) TiO2 (Anatase) Tetragonal 𝐼4! /𝑎𝑚𝑑 (141) MO2 (M=Ti, Ru, Ir) (Rutile) Tetragonal 𝑃4" /𝑚𝑛𝑚 (136) CuO Monoclinic 𝐶2/𝑐 (16) Polar Surfaces Material Space group ReO3 Cubic $ 𝑚 (221) 𝑃𝑚3 PdO Tetragonal 𝑃4" 𝑚𝑚𝑐 (131) ZnO, ZnS, CdSe (Wurzite) Hexagonal 𝑃6# 𝑚𝑐 (186) Tables 2 and 3 summarize the main features of the obtained models. Table 2 is focused on materials whose polar surfaces do not contribute to the Wulff structure. Table 3 is devoted to models with at least one polar surface. In each case, the nanoparticle center, stoichiometry, size and number of isomers are reported. For each center we only report two models. Models up to 5nm can be found in the Supporting Information (Table S4 and S5). Figure 2 shows two models per example without polar surfaces (stoichiometric nanoparticle and nanoparticle without dangling atoms) and Figure 3 reports two models for the materials including polar surfaces (metal or non-metal rich and with additional H or OH groups). It should be noted that the size difference between two contiguous models with the same center corresponds approximately to the twice of the interlayer distance in the direction defining the size of the model. 17 The models reported here range between 1 and 5 nm and aim to illustrate the output arising from BCN-M. The 1 to 5 nm range corresponds to the lower limit of the nanoworld and, for most of the cases, the provided nanoparticle models fall below the limit for which Wulff-like models are the thermodynamic preferred structure. However, we believe that several of the reported models can be used in further simulations as models for larger nanoparticles, in a similar manner to what is done in surface modelling, in which slab thickness rarely exceeds 1.5 nm. Table 2. Selected NP models without polar surfaces obtained with BCN-M. Two possible centers (1 and 2) are reported. For each one, the size of the NP model and total number of isomers after removing equivalent structures are shown. Material Center Formulaa Formulab Sizec Isomers MgO (0.25,0.25,0.25) (MgO)32,(MgO)108 N/A 1.0,1.7 1, 1 CeO2 Cation (CeO2)19,(CeO2)85 Ce19O32,Ce85O160 1.1, 2.2 1, 1000 (0.5,0.0,0.0) (CeO2)44, (CeO2)146 Ce44O80,Ce146O280 1.6, 2.7 2,1 Li2O Anion (Li2O)19,(Li2O)85 Li32O19, Li160O85 1.2, 2.1 1,999 (0.5,0.0,0.0) (Li2O)44, (Li2O)146 Li80O44, Li280O146 1.6,2.5 2,999 Cu2O Cation (Cu2O)40,(Cu2O)112 Cu55O40, Cu175O112 1.7, 2.6 12,1000 In2O3 (0.5, 0.0, 0.0) (In2O3)22, (In2O3)73 In44O60, In146O210 1.5, 2.6 1, 29 N/A N/A N/A N/A N/A e RuO2 (0.5, 0.5, 0.0) (RuO2)24,(RuO2)90, Ru24O42,Ru90O170 0.9, 1.5 19, 1000 Cation (RuO2)33,(RuO2)115 Ru33O60,Ru115O220 1.2, 1.9 7, 7 e IrO2 (0.5, 0.5, 0.0) (IrO2)24,(IrO2)90 Ir24O42,Ir90O170 0.9, 1.5 19, 986 Cation (IrO2)33,(IrO2)115 Ir33O60,Ir115O220 1.2, 1.8 7, 644 TiO2e (0.5, 0.5, 0.0) (TiO2)60,d (TiO2)240 Ti60O114,Ti240O470 2.1, 3.3 19, 1000 d Cation (TiO2)69, (TiO2)215 Ti69O132,Ti215O420 2.4, 3.0 7,1000 TiO2d (0.0, 0.0, 0.5) (TiO2)10, (TiO2)84 N/A 1.0, 2.9 1, 1 Cation (TiO2)35, (TiO2)165 N/A 2.0, 3.9 1, 1 CuO Anion (CuO)40,(CuO)110 Cu40O36, Cu110O105 1.6, 2.4 429,1000 a Stoichiometric nanoparticle Nanoparticle model without dangling atoms c In nm, measured b as the largest distance between cations, d Anatase, e Rutile, f two smaller models ((TiO2)42 and (TiO2)51) are also obtained and their Cartesian Coordinates reported in the ESI. 18 Figure 2. Selected NP models of the test cases that do not include polar surfaces Table 3. Selected NP models containing polar surfaces obtained with BCN-M. Two possible centers (1 and 2) are reported. Material Center Formulaa Formulab Sizec ReO3 Cation Re27O108,Re125O450 H54Re27O108,H150Re125O450 1.3, 2.6 (0.5, 0.5, 0.5) Re64O240,Re216O756 H96Re64O240,H216Re216O756 1.9, 3.2 PdO anion Pd28O35,Pd40O49 N/A 1.4, 1.9 d ZnO (0.0, 0.0, 0.0) Zn60O72, Zn216O243 H24Zn60O72,H54Zn216O243 1.4, 2.5 ZnSd (0.0, 0.0, 0.0) Zn108S96,Zn324S297 2.8, 3.9 d CdSe (0.0, 0.0, 0.0) Cd108Se120,Cd210Se231 H24Cd108Se120,H42Cd210Se231 3.1, 3.7 a Formula of models with either metal or non-metal reach termination, b formula of models with additional OH groups, c In nm, measured as the largest distance between cations, d Wurtzite 19 Figure 3. Selected NP models of the test cases that include polar surfaces. Optional parameters. Nanoparticle center. The position of the center of the nanoparticle for the Wulff construction is a key issue in the nanoparticle model generation and it allows obtaining different families of models (See Table 2). By default, BCN-M considers three types of centers: i) the atomic positions of the irreducible atoms, ii) the origin of coordinates of the irreducible unit cell, and iii) several special positions with low symmetry multiplicity, with the aim of exploring the center of the motifs that are being repeated in the crystal structure. Noticeably, the first type of center usually leads to models with an odd number of AXBY units, while the other two usually leads to an even number of AXBY units. However, for some particular cases only one set of centers generates stoichiometric nanoparticles. 20 Figure 4a shows two RuO2 nanoparticle models as an illustrative example of the importance of nanoparticle center. When the center of the Wulff construction is set at a cation position, models with an odd number of MO2 units are obtained. In contrast, when the center is located between two cations of the central row of the nanoparticle model (i.e., fractional coordinates of 0.5, 0.5, 0.0), another set of Wulff-like NP models, which are characterized by the presence of an even number of MO2 units, is obtained. The models with an even number of MO2 units present a larger portion of the {110} family of facets than the models with an odd number of units, the Wulff-like criterion suggesting that models the former are closer to the Wulff morphology for IrO2 and RuO2. Overall, our experience suggests that exploring these three center types for the Wulff construction is sufficient to ensure the generation of realistic and stoichiometric models. However, deciding which one of them is more appropriate in each case is not straightforward. Polar surface termination. Polar surfaces are generally highly unstable and thus, they are less common than non-polar ones. However, for some particular nanoparticles their presence is mandatory to achieve the Wulff morphology. The formation of polar surfaces is usually accompanied by an important surface reconstruction that reduces the polarity. As mentioned, three main reconstructions can occur:55–58 i) generation of metal and/or non-metal vacancies in the outermost layers, ii) electron density redistribution to achieve neutral non-stoichiometric surfaces (anion or metal reach), and iii) incorporation of new surface groups such as hydrides or hydroxyls. The computed surface energies for polar facets are usually associated with either a metal or an anion termination and consequently, the final model obtained with these surface energies cannot be stoichiometric. 21 Figure 4. a) Selected examples of RuO2 nanoparticle obtained with different centers. b) Selected examples of ZnO nanoparticles showing the three possible terminations in polar surfaces metal rich, anion rich or OH-terminated. c) Three models of (IrO2)33 that are considered equivalent by the SPRINT Coordinates. Atoms in blue correspond to the dangling oxygen atoms. Wurzite-like zinc oxide is a paradigmatic example. Experimental evidences suggest the presence of terraces and/or OH groups when synthetizing films involving the (001) surface depending on the conditions.56,57,93,94 Calculations for ZnO predict that the most stable termination of the polar basal (001) surface implies a shell of zinc atoms at the outermost layer in oxygen-poor environments and an oxygen termination in oxygen-rich conditions.65 As already mentioned, BCN-M is able to distinguish the presence of polar surfaces and treat them differently. Three different terminations are possible: i) metal- rich; ii) non-metal-rich, and iii) OH terminated. Generating models with BCN-M that include terraces and other large reconstructions described in the literature is at present 22 not possible. This kind of reconstructions highly depend on the material and synthetic conditions and thus, their modeling can hardly be automatized. Experienced users can construct this kind of models by using any of the three resulting models as starting points by removing/adding atoms. Figure 4b reports three models associated with a relatively small nanoparticle. The nanoparticle morphology depends on the ratio between the surface energy of the polar (001) surface and that of the stoichiometric (100) one. Since the surface energy for the (001) facet depends on the termination, while that of the (100) does not, the O-terminated and Zn-terminated nanoparticles present different morphologies. Moreover, the OH groups are added perpendicularly to the (001) facet. This leads to M-O-H bond angles of 180º, which are unrealistic. However, the most stable structure is difficult to determine a priori due the large number of H-bonding possibilities and thus, the OH terminated structures should be taken as preliminary models. Redundant structures. The “Remove ion excess loop” removes dangling ions until stoichiometry is achieved. When the number of dangling ions is larger than the number of atoms to be removed, several isomers can be generated, and equivalent models have to be identified. Figure 4c shows the representative example of (MO2)33 by showing three illustrative species that the SPRINT coordinates identify as identical structures. For this system the program generates 28 models, which are reduced to only 7 with different connectivity after using the SPRINT coordinates. However, it is worth mentioning that for nanoparticles larger than 2 nm, the number of potential isomers is so large (following a binomial coefficient of the number of dangling atoms in the last step of the “remove ion excess loop” and the number of atoms to be removed) that few equivalent models are usually generated from the statistical distribution. Thus, and since the computation 23 of the SPRINT coordinates would become very expensive for models larger than 2 nm, BCN-M turns off the SPRINT command Electronic structure of the different isomers. As a final test, the similarities/differences among the isomers generated by the BCN-M algorithm are analyzed by performing DFT calculations on the seven (RuO2)33 NPf isomers. (RuO2)33 has been taken as representative example due to the limited number of isomers and its reduced size. We expect that conclusions can be extrapolated to the other systems. The seven (RuO2)33 isomers have the same global coordination and Wulff-like quality criterion value and only differ on the distribution coefficient. Figure 6 shows the optimized geometries for the seven models and the relative energies of the optimized structures with respect to the most stable isomer (in kJ mol-1). Figure 7 displays: a) the distribution of the Ru-O distance and O-Ru-O angle corresponding to atoms at the surface, and b) the projected Density of States (PDOS) associated with the d orbitals of surface Ru atoms and the p orbitals of the surface oxygen atoms of two selected isomers. The PDOS of all other (RuO2)33 species can be found in Figure S3 of the Supporting Information. The initial structures obtained with the BCN-M algorithm shows six dangling O atoms. Four of these dangling atoms are located at the four tetracoordinated Ru centers of the intersection between two surfaces of the {011} family and one surface of the {110} family. The other two are located at pentacoordinated Ru centers of the intersection of two surfaces of the {110} family and one of the {011}. Thus, the different isomers arise from the different distribution of the two oxygens in these eight sites. 24 Figure 6. Optimized structures of the seven (RuO2)33 isomers, distribution coefficient (DC) in pm and relative energies with respect to the most stable isomer in kJ mol-1. The values in parenthesis correspond to the relative energies per RuO2 unit. Dangling oxygen atoms are represented in blue. During the DFT optimizations, the seven isomers behave similarly. No surface reconstruction is observed, but a larger reorganization than that observed on surfaces takes place.63 In general, the Ru-O distances are shorter than those of the bulk, which are 1.942 Å (apical) and 1.998 Å (equatorial), and the effect is more pronounced for the Ru-O bond trans to a vacancy (at least 0.06 Å). The ruthenium center exposed on the (110) surface also moves inside the nanoparticle enlarging the Oax-Ru-Oeq angle up to around 1000 for most cases, as compared to the 900 of the bulk. The major reorganization is observed on ruthenium centers presenting dangling atoms. The Ru- Odangling distance is even shorter (1.67-1.68 Å in the tetracoordinated centers and 1.70 Å in the pentacoordinated ones), which suggests the formation of a Ru=Oxo species to 25 maintain ruthenium formal oxidation state to +4. Moreover, the coordination environment of these ruthenium centers becomes a distorted trigonal bipyramid for pentacoordinated ones and a distorted tetrahedron for the tetracoordinated ones. This implies several angles between 1500 and 1050 as can be seen in Figure 8a. This reorganization occurs for all the isomers considered, regardless of the position of the dangling oxygen atoms in the NP. As a consequence, the largest energy difference between two isomers is rather small (27.3 kJ mol-1 or 0.8 kJ mol-1 per RuO2 unit), the relative stabilities between isomers arising mainly from subtle reconstruction differences between the seven nanoparticles. Figure 7. a) Ru-O distance and O-Ru-O angle distribution of the atoms on the nanoparticle surface. b) Projected Density of States of the most stable nanoparticle isomer (4) and the highest in energy one (6) The Projected Density Of States (PDOS) of the valence Ru d orbitals and O p orbitals of surface atoms are also very close for all isomers. All of them present two different bands. 26 The first one appears between -9.0/-8.8 eV and -2.1/-1.9 eV and the second one spreads between -2.1/-1.9 and +4.9/+5.1 eV. Main contribution to the former comes from ruthenium d orbitals (74-75%, see Table S6 of the Supporting Information) whereas the major contribution to the second band comes from the p orbitals of oxygen (65 – 67 %). For the two bands, several distinctive peaks can be distinguished, the highest ones appearing at very similar values and presenting similar shapes. Overall, the seven structures generated by BCN-M have similar geometric and electronic structure, pointing out they all could behave similarly. Conclusions Construction of atomistic nanoparticle models of binary ionic compounds is not always straightforward because it usually involves post-processing to obtain stoichiometric and charge electroneutral systems. This implies manual modification of the NP structures, which is under user subjectivity and bias. Bulk Cut Nanoparticle Models (BCN-M) is a computational tool delivered as a free web platform (http://bcnm.qf.uab.cat) and downloadable utility able to generate Wulff like nanoparticle models with controlled stoichiometry and surface termination. Models are constructed from the crystallographic data of the selected material (cell parameters of the conventional unit cell, space group symmetry and fractional coordinates of the irreducible atoms), the surface energies of the main crystallographic facets and information about the desired surface termination. BCN-M explores several different nanoparticle centers to ensure different sets of nanoparticles usually covering even and odd number of AXBY units. Moreover, the quality of the resulting models is analysed with: i) the Wulff-like criterion, 27 which measures the deviation from ideal Wulff shape; ii) the global coordination of the atoms of the model and iii) the Distribution Coefficient, which evaluates the polarity of the model. BCN-M has been tested and applied to fifteen cases covering the most representative cubic, tetragonal, hexagonal and monoclinic crystallographic systems as well as systems involving only non-polar surfaces or a combination of both polar and non-polar ones. The resulting models are listed, and their Cartesian Coordinates reported for further uses. BCN-M models are generated from cutting the bulk structure, and accordingly they do not take into account the most likely surface reconstruction events. In this context, these models can be used as starting points for further simulations or structure analysis of ideal systems. Because of its user-friendly handling, BCN-M can facilitate the developments of structure – property relationships in nanoscience i.e. atomistic models of nanoparticles can provide information on the number and nature of surface active sites, thereby relating their size and shape with its surface properties. In addition, BCN- M tool can help evaluating how changes in the relative surface energies, which can be modified through capping agents, can influence the nanoparticle shape and surface properties, assisting in the design of new catalysts. Acknowledgements. The authors gratefully acknowledge financial support from MINECO (CTQ2017-89132-P) and the Generalitat de Catalunya (2017SGR1323) as well as the computational time in the Red Española de Supercomputación (Projects QCM- 2018-2-0009, QCM-2018-3-0021, QCM-2019-1-0043 and QS-2019-2-0031). We are indebted to the anonymous reviewers for their constructive suggestions and fruitful 28 remarks. XSM is grateful for the Professor Agregat Serra Húnter position. AR is indebted to “Ramón y Cajal” program. Supporting Information. The supporting information including the description of the process followed in the “remove ions excess loop”, the nanoparticle centers considered as function of the Bravais lattices, input data for obtaining the models for the test systems, details on larger nanoparticle models for each test example, the Projected Density of States of the Ru d and O p orbitals of the surface atoms of the seven (RuO2)33 isomers, the Flow Chart for obtaining: i) nanoparticle models without dangling atoms and ii) nanoparticles involving polar surfaces and the Cartesian Coordinates of all models reported in a .xyz file can be found free of charge at the ACS Publication website. 29 References (1) Sapsford, K. E.; Algar, W. R.; Berti, L.; Gemmill, K. B.; Casey, B. J.; Oh, E.; Stewart, M. H.; Medintz, I. L. Functionalizing Nanoparticles with Biological Molecules: Developing Chemistries That Facilitate Nanotechnology. Chem. Rev. 2013, 113, 1904–2074. (2) Polshettiwar, V.; Luque, R.; Fihri, A.; Zhu, H.; Bouhrara, M.; Basset, J.-M. Magnetically Recoverable Nanocatalysts. Chem. Rev. 2011, 111, 3036–3075. (3) Polshettiwar, V.; Varma, R. S. Green Chemistry by Nano-Catalysis. Green Chem. 2010, 12, 743. (4) Somorjai, G. A.; Frei, H.; Park, J. Y. Advancing the Frontiers in Nanocatalysis, Biointerfaces, and Renewable Energy Conversion by Innovations of Surface Techniques. J. Am. Chem. Soc. 2009, 131, 16589–16605. (5) Liu, D.; Kelly, T. L. Perovskite Solar Cells with a Planar Heterojunction Structure Prepared Using Room-Temperature Solution Processing Techniques. Nat. Photonics 2014, 8, 133–138. (6) Kim, H.-S.; Lee, C.-R.; Im, J.-H.; Lee, K.-B.; Moehl, T.; Marchioro, A.; Moon, S.-J.; Humphry-Baker, R.; Yum, J.-H.; Moser, J. E.; Grätzel, M.; Park, N.-G. Lead Iodide Perovskite Sensitized All-Solid-State Submicron Thin Film Mesoscopic Solar Cell with Efficiency Exceeding 9%. Sci. Rep. 2012, 2, 591. (7) Etgar, L.; Gao, P.; Xue, Z.; Peng, Q.; Chandiran, A. K.; Liu, B.; Nazeeruddin, M. K.; Grätzel, M. Mesoscopic CH3 NH3 PbI3/TiO2 Heterojunction Solar Cells. J. Am. Chem. Soc. 2012, 134, 17396–17399. (8) Nordmeyer, D.; Stumpf, P.; Gröger, D.; Hofmann, A.; Enders, S.; Riese, S. B.; Dernedde, J.; Taupitz, M.; Rauch, U.; Haag, R.; Rühl, E.; Graf, C. Iron Oxide 30 Nanoparticles Stabilized with Dendritic Polyglycerols as Selective MRI Contrast Agents. Nanoscale 2014, 6, 9646–9654. (9) Wu, K. C.-W.; Yamauchi, Y.; Hong, C.-Y.; Yang, Y.-H.; Liang, Y.-H.; Funatsu, T.; Tsunoda, M. Biocompatible, Surface Functionalized Mesoporous Titania Nanoparticles for Intracellular Imaging and Anticancer Drug Delivery. Chem. Commun. 2011, 47, 5232. (10) Lagopati, N.; Kitsiou, P. V.; Kontos, A. I.; Venieratos, P.; Kotsopoulou, E.; Kontos, A. G.; Dionysiou, D. D.; Pispas, S.; Tsilibary, E. C.; Falaras, P. Photo-Induced Treatment of Breast Epithelial Cancer Cells Using Nanostructured Titanium Dioxide Solution. J. Photochem. Photobiol. A Chem. 2010, 214, 215–223. (11) McNamara, K.; Tofail, S. A. M. Nanoparticles in Biomedical Applications. Adv. Phys. X 2017, 2, 54–88. (12) McNamara, K.; Tofail, S. A. M. Nanosystems: The Use of Nanoalloys, Metallic, Bimetallic, and Magnetic Nanoparticles in Biomedical Applications. Phys. Chem. Chem. Phys. 2015, 17, 27981–27995. (13) Xie, X.; Li, Y.; Liu, Z.-Q.; Haruta, M.; Shen, W. Low-Temperature Oxidation of CO Catalysed by Co3O4 Nanorods. Nature 2009, 458, 746–749. (14) Campbell, C. T. PHYSICS: The Active Site in Nanoparticle Gold Catalysis. Science 2004, 306, 234–235. (15) Grass, M. E.; Zhang, Y.; Butcher, D. R.; Park, J. Y.; Li, Y.; Bluhm, H.; Bratlie, K. M.; Zhang, T.; Somorjai, G. A. A Reactive Oxide Overlayer on Rhodium Nanoparticles during CO Oxidation and Its Size Dependence Studied by In Situ Ambient- Pressure X-Ray Photoelectron Spectroscopy. Angew. Chemie Int. Ed. 2008, 47, 8893–8896. 31 (16) Kim, I. H.; Seo, H. O.; Park, E. J.; Han, S. W.; Kim, Y. D. Low Temperature CO Oxidation over Iron Oxide Nanoparticles Decorating Internal Structures of a Mesoporous Alumina. Sci. Rep. 2017, 7, 40497. (17) Yu, X.; Shavel, A.; An, X.; Luo, Z.; Ibáñez, M.; Cabot, A. Cu2ZnSnS4-Pt and Cu2ZnSnS4-Au Heterostructured Nanoparticles for Photocatalytic Water Splitting and Pollutant Degradation. J. Am. Chem. Soc. 2014, 136, 9236–9239. (18) Mangrulkar, P. A.; Polshettiwar, V.; Labhsetwar, N. K.; Varma, R. S.; Rayalu, S. S. Nano-Ferrites for Water Splitting: Unprecedented High Photocatalytic Hydrogen Production under Visible Light. Nanoscale 2012, 4, 5202–5209. (19) Jaramillo, T. F.; Jørgensen, K. P.; Bonde, J.; Nielsen, J. H.; Horch, S.; Chorkendorff, I. Identification of Active Edge Sites for Electrochemical H2 Evolution from MoS2 Nanocatalysts. Science 2007, 317, 100–102. (20) Astruc, D.; Lu, F.; Aranzaes, J. R. Nanoparticles as Recyclable Catalysts: The Frontier between Homogeneous and Heterogeneous Catalysis. Angew. Chemie Int. Ed. 2005, 44, 7852–7872. (21) Polshettiwar, V.; Asefa, T. Introduction to Nanocatalysis. In Nanocatalysis Synthesis and Applications; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2013; pp 1–9. (22) Titus, D.; James Jebaseelan Samuel, E.; Roopan, S. M. Nanoparticle Characterization Techniques. In Green Synthesis, Characterization and Applications of Nanoparticles; Elsevier, 2019; pp 303–319. (23) Mourdikoudis, S.; Pallares, R. M.; Thanh, N. T. K. Characterization Techniques for Nanoparticles: Comparison and Complementarity upon Studying Nanoparticle Properties. Nanoscale 2018, 10, 12871–12934. 32 (24) Billinge, S. J. L.; Levin, I. The Problem with Determining Atomic Structure at the Nanoscale. Science 2007, 316, 561–565. (25) Gong, X.-Q.; Selloni, A.; Batzill, M.; Diebold, U. Steps on Anatase TiO2(101). Nat. Mater. 2006, 5, 665–670. (26) Weissenrieder, J.; Kaya, S.; Lu, J.-L.; Gao, H.-J.; Shaikhutdinov, S.; Freund, H.-J.; Sierka, M.; Todorova, T. K.; Sauer, J. Atomic Structure of a Thin Silica Film on a Mo(112) Substrate: A Two-Dimensional Network of SiO4 Tetrahedra. Phys. Rev. Lett. 2005, 95, 076103. (27) Chen, M.; Felmy, A. R.; Dixon, D. A. Structures and Stabilities of (MgO)n Nanoclusters. J. Phys. Chem. A 2014, 118, 3136–3146. (28) Bromley, S. T.; Flikkema, E. Columnar-to-Disk Structural Transition in Nanoscale (SiO2)N Clusters. Phys. Rev. Lett. 2005, 95, 185505. (29) Lamiel-Garcia, O.; Cuko, A.; Calatayud, M.; Illas, F.; Bromley, S. T. Predicting Size- Dependent Emergence of Crystallinity in Nanomaterials: Titania Nanoclusters versus Nanocrystals. Nanoscale 2017, 9, 1049–1058. (30) Chen, M.; Dixon, D. A. Machine-Learning Approach for the Development of Structure–Energy Relationships of ZnO Nanoparticles. J. Phys. Chem. C 2018, 122, 18621–18639. (31) Guo, H.; Barnard, A. S. Thermodynamic Modelling of Nanomorphologies of Hematite and Goethite. J. Mater. Chem. 2011, 21, 11566. (32) Lamiel-Garcia, O.; Ko, K. C.; Lee, J. Y.; Bromley, S. T.; Illas, F. When Anatase Nanoparticles Become Bulklike: Properties of Realistic TiO2 Nanoparticles in the 1–6 Nm Size Range from All Electron Relativistic Density Functional Theory Based Calculations. J. Chem. Theory Comput. 2017, 13, 1785–1793. 33 (33) Bruno, M.; Massaro, F. R.; Prencipe, M.; Demichelis, R.; De La Pierre, M.; Nestola, F. Ab Initio Calculations of the Main Crystal Surfaces of Forsterite (Mg2SiO4): A Preliminary Study to Understand the Nature of Geochemical Processes at the Olivine Interface. J. Phys. Chem. C 2014, 118, 2498–2506. (34) Ringe, E.; Van Duyne, R. P.; Marks, L. D. Wulff Construction for Alloy Nanoparticles. Nano Lett. 2011, 11, 3399–3403. (35) Woodley, S. M.; Lazauskas, T.; Illingworth, M.; Carter, A. C.; Sokol, A. A. What Is the Best or Most Relevant Global Minimum for Nanoclusters? Predicting, Comparing and Recycling Cluster Structures with WASP@N. Faraday Discuss. 2018, 211, 593–611. (36) Viñes, F.; Lamiel-Garcia, O.; Illas, F.; Bromley, S. T. Size Dependent Structural and Polymorphic Transitions in ZnO: From Nanocluster to Bulk. Nanoscale 2017, 9, 10067–10074. (37) Loschen, C.; Migani, A.; Bromley, S. T.; Illas, F.; Neyman, K. M. Density Functional Studies of Model Cerium Oxide Nanoparticles. Phys. Chem. Chem. Phys. 2008, 10, 5730. (38) Migani, A.; Neyman, K. M.; Illas, F.; Bromley, S. T. Exploring Ce3+/Ce4+ Cation Ordering in Reduced Ceria Nanoparticles Using Interionic-Potential and Density- Functional Calculations. J. Chem. Phys. 2009, 131, 064701. (39) Wales, D. J.; Doye, J. P. K. Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms. J. Phys. Chem. A 1997, 101, 5111–5116. (40) Rahm, J. M.; Erhart, P. Beyond Magic Numbers: Atomic Scale Equilibrium Nanoparticle Shapes for Any Size. Nano Lett. 2017, 17, 5775–5781. 34 (41) Cleveland, C. L.; Landman, U. The Energetics and Structure of Nickel Clusters: Size Dependence. J. Chem. Phys. 1991, 94, 7376–7396. (42) Barnard, A. S. A Thermodynamic Model for the Shape and Stability of Twinned Nanostructures. J. Phys. Chem. B 2006, 110, 24498–24504. (43) Marks, L. D. Modified Wulff Constructions for Twinned Particles. J. Cryst. Growth 1983, 61, 556–566. (44) Marks, L. D.; Peng, L. Nanoparticle Shape, Thermodynamics and Kinetics. J. Phys. Condens. Matter 2016, 28, 053001. (45) Wang, Y.; He, J.; Liu, C.; Chong, W. H.; Chen, H. Thermodynamics versus Kinetics in Nanosynthesis. Angew. Chemie Int. Ed. 2015, 54, 2022–2051. (46) Barmparis, G. D.; Lodziana, Z.; Lopez, N.; Remediakis, I. N. Nanoparticle Shapes by Using Wulff Constructions and First-Principles Calculations. Beilstein J. Nanotechnol. 2015, 6, 361–368. (47) Novell-Leruth, G.; Carchini, G.; López, N. On the Properties of Binary Rutile MO2 Compounds, M = Ir, Ru, Sn, and Ti: A DFT Study. J. Chem. Phys. 2013, 138, 194706. (48) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; Noël, Y.; Causà, M.; Rérat, M.; Kirtman, B. CRYSTAL14: A Program for the Ab Initio Investigation of Crystalline Solids. Int. J. Quantum Chem. 2014, 114, 1287–1317. (49) Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272–1276. (50) Zucker, R. V.; Chatain, D.; Dahmen, U.; Hagège, S.; Carter, W. C. New Software Tools for the Calculation and Display of Isolated and Attached Interfacial-Energy 35 Minimizing Particle Shapes. J. Mater. Sci. 2012, 47, 8290–8302. (51) Hjorth Larsen, A.; Jørgen Mortensen, J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Bjerre Jensen, P.; Kermode, J.; Kitchin, J. R.; Leonhard Kolsbjerg, E.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; et al. The Atomic Simulation Environment—a Python Library for Working with Atoms. J. Phys. Condens. Matter 2017, 29, 273002. (52) Scopece, D. SOWOS : An Open-Source Program for the Three-Dimensional Wulff Construction. J. Appl. Crystallogr. 2013, 46, 811–816. (53) Mathew, K.; Singh, A. K.; Gabriel, J. J.; Choudhary, K.; Sinnott, S. B.; Davydov, A. V.; Tavazza, F.; Hennig, R. G. MPInterfaces: A Materials Project Based Python Tool for High-Throughput Computational Screening of Interfacial Systems. Comput. Mater. Sci. 2016, 122, 183–190. (54) Chatzigoulas, A.; Karathanou, K.; Dellis, D.; Cournia, Z. NanoCrystal: A Web- Based Crystallographic Tool for the Construction of Nanoparticles Based on Their Crystal Habit. J. Chem. Inf. Model. 2018, 58, 2380–2386. (55) Diebold, U.; Li, S.-C.; Schmid, M. Oxide Surface Science. Annu. Rev. Phys. Chem. 2010, 61, 129–148. (56) Noguera, C. Polar Oxide Surfaces. J. Phys. Condens. Matter 2000, 12, R367– R410. (57) Goniakowski, J.; Finocchi, F.; Noguera, C. Polarity of Oxide Surfaces and Nanostructures. Reports Prog. Phys. 2008, 71, 016501. (58) Andersen, T. K.; Fong, D. D.; Marks, L. D. Pauling’s Rules for Oxide Surfaces. Surf. Sci. Rep. 2018, 73, 213–232. 36 (59) Rogal, J.; Reuter, K.; Scheffler, M. Thermodynamic Stability of PdO Surfaces. Phys. Rev. B 2004, 69, 075421. (60) Seriani, N. Ab Initio Thermodynamics of Lithium Oxides: From Bulk Phases to Nanoparticles. Nanotechnology 2009, 20, 445703. (61) Pantaleone, S.; Rimola, A.; Sodupe, M. Canonical, Deprotonated, or Zwitterionic? A Computational Study on Amino Acid Interaction with the TiO2 (101) Anatase Surface. J. Phys. Chem. C 2017, 121, 14156–14165. (62) González, D.; Heras-Domingo, J.; Pantaleone, S.; Rimola, A.; Rodríguez-Santiago, L.; Solans-Monfort, X.; Sodupe, M. Water Adsorption on MO2 (M = Ti, Ru, and Ir) Surfaces. Importance of Octahedral Distortion and Cooperative Effects. ACS Omega 2019, 4, 2989–2999. (63) Heras-Domingo, J.; Sodupe, M.; Solans-Monfort, X. Interaction between Ruthenium Oxide Surfaces and Water Molecules. Effect of Surface Morphology and Water Coverage. J. Phys. Chem. C 2019, 123, 7786–7798. (64) Wilson, H. F.; Tang, C.; Barnard, A. S. Morphology of Zinc Oxide Nanoparticles and Nanowires: Role of Surface and Edge Energies. J. Phys. Chem. C 2016, 120, 9498–9505. (65) Tang, C.; Spencer, M. J. S.; Barnard, A. S. Activity of ZnO Polar Surfaces: An Insight from Surface Energies. Phys. Chem. Chem. Phys. 2014, 16, 22139–22144. (66) Feigl, C. A.; Barnard, A. S.; Russo, S. P. Size- and Shape-Dependent Phase Transformations in Wurtzite ZnS Nanostructures. Phys. Chem. Chem. Phys. 2012, 14, 9871. (67) Manna, L.; Wang; Cingolani, R.; Alivisatos, A. P. First-Principles Modeling of Unpassivated and Surfactant-Passivated Bulk Facets of Wurtzite CdSe: A Model 37 System for Studying the Anisotropic Growth of CdSe Nanocrystals. J. Phys. Chem. B 2005, 109, 6183–6192. (68) Mishra, A. K.; Roldan, A.; de Leeuw, N. H. CuO Surfaces and CO2 Activation: A Dispersion-Corrected DFT+ U Study. J. Phys. Chem. C 2016, 120, 2198–2214. (69) Hu, M.; Zhang, J.; Wang, W.-D.; Qin, Y.-X. Ab-Initio Density Functional Theory Study of a WO3 NH3 -Sensing Mechanism. Chinese Phys. B 2011, 20, 082101. (70) Nolan, M. Enhanced Oxygen Vacancy Formation in Ceria (111) and (110) Surfaces Doped with Divalent Cations. J. Mater. Chem. 2011, 21, 9160. (71) Skorodumova, N. V.; Baudin, M.; Hermansson, K. Surface Properties of CeO2 from First Principles. Phys. Rev. B 2004, 69, 075401. (72) Paier, J.; Penschke, C.; Sauer, J. Oxygen Defects and Surface Chemistry of Ceria: Quantum Chemical Studies Compared to Experiment. Chem. Rev. 2013, 113, 3949–3985. (73) Geysermans, P.; Finocchi, F.; Goniakowski, J.; Hacquart, R.; Jupille, J. Combination of (100), (110) and (111) Facets in MgO Crystals Shapes from Dry to Wet Environment. Phys. Chem. Chem. Phys. 2009, 11, 2228. (74) Suchitra; Pan, J.; Waghmare, U. V. High Tunability of the Work Function of (001) Surface of ReO 3 with O-Vacancies: First Principles Analysis. J. Appl. Phys. 2014, 116, 034304. (75) Bendavid, L. I.; Carter, E. A. First-Principles Predictions of the Structure, Stability, and Photocatalytic Potential of Cu2O Surfaces. J. Phys. Chem. B 2013, 117, 15750–15760. (76) Walsh, A.; Catlow, C. R. A. Structure, Stability and Work Functions of the Low Index Surfaces of Pure Indium Oxide and Sn-Doped Indium Oxide (ITO) from 38 Density Functional Theory. J. Mater. Chem. 2010, 20, 10438. (77) Pantaleone, S.; Ugliengo, P.; Sodupe, M.; Rimola, A. When the Surface Matters: Prebiotic Peptide-Bond Formation on the TiO 2 (101) Anatase Surface through Periodic DFT-D2 Simulations. Chem. - A Eur. J. 2018, 24, 16292–16301. (78) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (79) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. (80) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561. (81) Grimme, S. Accurate Description of van Der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463–1473. (82) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (83) van der Walt, S.; Colbert, S. C.; Varoquaux, G. The NumPy Array: A Structure for Efficient Numerical Computation. Comput. Sci. Eng. 2011, 13, 22–30. (84) Oliphant, T. E. Python for Scientific Computing. Comput. Sci. Eng. 2007, 9, 10– 20. (85) Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (Pymatgen): A Robust, Open-Source Python Library for Materials Analysis. Comput. Mater. Sci. 2013, 68, 314–319. 39 (86) PyYAML https://pyyaml.org/. (87) Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (Pymatgen): A Robust, Open-Source Python Library for Materials Analysis. Comput. Mater. Sci. 2013, 68, 314–319. (88) Togo, A.; Tanaka, I. Spglib: a software library for crystal symmetry search http://arxiv.org/abs/1808.01590. (89) The Gibbs Surface Is Defined as the Ratio between Center-Surface Distance and the Interlayer Distance plus Half of the Interlayer Distance (s = Floor(d/L) + 0.5L Where d Is the Center-Surface Distance and L Is the Interlayer Separation. (90) Rodríguez, A.; Rodríguez-Fernández, R.; A. Vázquez, S.; L. Barnes, G.; J. P. Stewart, J.; Martínez-Núñez, E. Tsscds2018: A Code for Automated Discovery of Chemical Reaction Mechanisms and Solving the Kinetics. J. Comput. Chem. 2018, 39, 1922–1930. (91) Pietrucci, F.; Andreoni, W. Graph Theory Meets Ab Initio Molecular Dynamics: Atomic Structures and Transformations at the Nanoscale. Phys. Rev. Lett. 2011, 107, 085504. (92) Honkala, K.; Hellman, A.; Remediakis, I. N.; Logadottir, A.; Carlsson, A.; Dahl, S.; Christensen, C. H.; Nørskov, J. K. Ammonia Synthesis from First-Principles Calculations. Science 2005, 307, 555–558. (93) Wöll, C. The Chemistry and Physics of Zinc Oxide Surfaces. Prog. Surf. Sci. 2007, 82, 55–120. (94) Diebold, U.; Koplitz, L.; Dulub, O. Atomic-Scale Properties of Low-Index ZnO Surfaces. Appl. Surf. Sci. 2004, 237, 336–342. 40 TOC Graphic BCN Models ✓ q Crystal Structure Stoichiometric q Surface Energies NP models q Size Range 41