In silico characterization of nanoparticles

Abstract

Nanoparticles (NPs) make for intriguing heterogeneous catalysts due to their large active surface area and excellent and often size-dependent catalytic properties that emerge from a multitude of chemically different surface reaction sites. NP catalysts are, in principle, also highly tunable: even small changes to the NP size or surface facet composition, doping with heteroatoms, or changes of the supporting material can significantly alter their physicochemical properties. Because synthesis of size- and shape-controlled NP catalysts is challenging, the ability to computationally predict the most favorable NP structures for a catalytic reaction of interest is an in-demand skill that can help accelerate and streamline the material optimization process. Fundamentally, simulations of NP model systems present unique challenges to computational scientists. Not only must considerable methodological hurdles be overcome in performing calculations with hundreds to thousands of atoms while retaining appropriate accuracy to be able to probe the desired properties. Also, the data generated by simulations of NPs are typically more complex than data from simulations of, for example, single crystal surface models, and therefore often require different data analysis strategies. To this end, the present work aims to review analytical methods and data analysis strategies that have proven useful in extracting thermodynamic trends from NP simulations.

---

1 Introduction

Deployment of heterogeneous catalysts in nanoparticulate form has various benefits. Not only do nanoparticles (NPs) promise high mass activity due to the more favorable surface-to-volume ratio compared to catalysts with grain sizes in the micro- or millimeter range, NPs can also express heightened or outright different catalytic properties compared to the bulk material as a result of finite- and quantum-size effects and due to the availability of a plethora of chemically distinct surface reaction sites. However, synthesis of shape- and size-controlled NPs is challenging, and Ostwald ripening as well as other degradation effects can impact the longevity of NP catalysts.¹ In order to streamline material optimization cycles, interest has therefore been growing in computational approaches that can predict favorable catalyst candidate structures for a reaction of interest. More and more, computational science is asked to establish structure-activity relationships for model NP catalysts in the 1 to 5 nm range and to investigate degradation mechanisms.

Historically, computational investigation of NP properties was usually carried out by performing density functional theory (DFT) calculations of single crystal model surfaces that correspond to the surface facets of a NP of interest.^2–6 However, results from surface models are often not transferable to NPs for several reasons. Firstly, the number of chemically different surface sites on a NP is larger than on single crystal surface models. In particular, the highly reactive, undercoordinated edge and vertex sites are hard to represent using a surface model.⁷ Furthermore, when using stepped model surfaces to mimic the lower coordination of NP edge sites, one finds that such models contain both convex and concave surface structures while NPs typically only contain one type.⁸ Another disparity between surface models and actual NPs is that quantum-size effects can affect the properties of clusters and small NPs in non-systematic ways.^9,10 In fact, not all properties observed for NPs will converge to the bulk limit as the system size increases.¹¹ Surface models are therefore unsuitable to study certain NP properties, even if the surface model is used as a stand-in for very large NPs. Finally, investigation of catalyst degradation processes using surface models and DFT can be arduous given the computational limitations with regards to system size and time scale that such processes typically occur at.

Nevertheless, the surface model approach remained the default option for a long time to model the behavior of extended facets of NPs because accurate simulations of NPs with hundreds or thousands of atoms are computationally demanding. Recent advances in computational methodology have changed this circumstance. Force field-based methods such as ReaxFF^12,13 or semi-empirical methods such as xTB¹⁴ and Sutton-Chen type potentials^15–18 have brought down costs for calculations of NPs significantly while retaining reasonable accuracy. Similarly, highly scalable codes enable electronic structure calculations of model sizes previously thought inaccessible using DFT. Advances in the field of hybrid simulations such as the QMMM ansatz also push the limits of system size.^19–23

With computational costs for NP simulations brought down, one important question remains: how best to extract useful information from the model systems and simulations? Analyzing properties and the reactivity of NPs is not as straight forward as analyzing a single crystal surface for example. The abundance of different binding sites will naturally lead to noisier results. Oftentimes, molecular dynamics (MD) or Monte Carlo (MC) based simulations are used to study NP model systems. Such simulations can create large data sets that require specialized data analysis strategies to extract the desired information.

This work aims to review methods of data analysis that we found useful in our research in this field^24–26 to extract information about thermodynamic properties and structure-activity relationships from simulations of NP catalysts. Exemplary studies applying the presented methods are highlighted when appropriate, without the intention to be exhaustive in this regard.

Many of the methods discussed in this review will be illustrated using atomic configurations of 3 nm octahedral Pt NP structures in various states of oxidation. These structures are taken from an openly accessible data set available under the https://doi.org/10.5281/zenodo.6322004 which our group recently published alongside a study comparing electrochemical oxidation trends of Pt NPs of various shapes.²⁶ Additionally, we have made available a git repository of a Jupyter book that serves as a tutorial on how to implement many of the methods presented below using Python3: https://bjk24.gitlab.io/in-silico-review/intro.html. The code is written to directly interact with the data set for oxidized octahedral NPs. More information on the required setup to follow the tutorial and execute the code for yourself can be found on the website linked above.

Of course, almost all of the presented methods are already very efficiently implemented in existing software packages for daily use. However, such programs are typically closed-source or at least hard to decipher due to radical optimization for speed. We therefore encourage readers, who we anticipate are scientists first and programmers second, to run and modify our more straight forward test script suite to familiarize themselves with the implementation of these methods and with code-driven data analysis in general.

The remainder of this review is structured as follows: Section 2 shortly summarizes open-source software solutions that enable users to generate atomic configurations of NPs. Section 3 is concerned with analytical methods aimed at extracting information from individual structures while Section 4 presents methods aimed at analyzing entire data sets. Some concluding remarks are given in Section 5. All figures original to this work used the Visual Molecular Dynamics (VMD) software²⁷ and the Tachyon renderer²⁸ to visualize atomic configurations. If not indicated otherwise, Pt atoms are grey while oxygen atoms are red.

2 Generating nanoparticle model systems

Many commercial computational chemistry software suites allow users to generate NP configurations. Here, three open-source alternatives are mentioned to reduce the activation barrier for newcomers.

• Atomic simulation environment (ASE): ASE is a Python3-based framework for molecular and materials simulations that interfaces with many commercial and open-source computational codes.²⁹ The ase.cluster class provides various means for generating NPs based on user-provided information about the lattice type and constants as well as Miller indices of the surfaces that should be cut. Pre-configured methods exist to generate typical NP shapes such as octahedra. A method for Wulff constructions exists as well which requires users to supply the formation energy of different types of surface facets. The method will then create a NP shape that minimizes the overall surface energy.

• WulffPack: WulffPack is a powerful python3 package used to generate NPs using Wulff construction.³⁰ WulffPack features built in methods to visualize the generated NPs and interfaces with ASE to export atomic coordinates.

• Nanocut: nanocut is a python3 package that can be used to generate NPs by providing cell and Miller index information. The program can also generate periodic surfaces and crystals, nanotubes, and nanowires.³¹ The ability to chain sequences of cutting procedures allows users to realize unusual shapes. At the time of writing, development appears to have halted. Nanocut may still be worth a look for users who find other alternatives too restrictive as compatibilty issues should be quick to fix due to the open-source nature of the project.

3 Analyzing single structures

3.1 Quantifying the surface coverage of nanoparticles

In computational surface science, coverage is typically given in terms of monolayers, i.e. with respect to the number of atoms within the top layer of the underlying periodic surface model. This is often a convenient measure because it directly relates to single-crystal experiments for example. In case of NPs, however, a monolayer is somewhat ill-defined due to the inhomogeneous nature of the available surface sites. Experimental data on NPs, for example from X-ray spectroscopy methods, can give insight into the relative occurrence of different chemical species or oxidation states in a material. Hence, giving coverage as a ratio of x = number of adatoms/number of substrate atoms is found to be more useful.^24,32–34 Senftle and co-workers for example have characterized the degree of oxidation of a Pd NP via its O : Pd ratio.³² This method of designating coverage will be used throughout this review.

3.2 Coordination number-based approaches

3.2.1 Classical coordination numbers. The classical coordination number (CN) of a given atom is calculated by counting its nearest-neighbor atoms. CNs have been used to classify adsorption positions on model surfaces for decades.³⁵ The usage of CNs to address adsorption sites on NPs is advantageous compared to the commonly used way of referring to adsorption sites on extended surfaces by means of the underlying surface composition (FCC, HCP, threefold, etc.) simply due to the sheer number of different available binding sites.

Moreover, CNs have been shown to be useful as a descriptor for catalytic activity. In 2009, Jiang and co-workers showed that the CNs of surface atoms on small NPs linearly correlate with metal d-band centers,³⁶ priming CNs as descriptors for reactivity.³⁷ Mpourmpakis and co-workers exploited this relationship further: they used CNs to characterize the interactions between small gold clusters (Au₁₆ to Au₄₅) and CO molecules.³⁸ From these results, the group extrapolated a transferable concept that uses the CN and the bonding angle of an adspecies on the particle surface as descriptors for the binding energy. This approach therefore allows for the prediction of NP reactivity without requiring expensive electronic structure calculations. Certainly, this concept cannot be generalized to any surface and adspecies as it neglects various aspects such as surface rearrangement and crowding effects. To address some of these limitations, the group introduced a refined descriptor approach based on the relationship between CO binding energy and surface CN that takes into account interdependence of NP structure and adsorbate coverage in a follow-up study.³⁹

Two interesting applications are highlighted in the following to further illustrate the usefulness of CNs for catalysis research. Using O atoms and CO molecules as probes, Kleis et al. employed CNs to identify the point of transition at which the adsorption properties on cuboctahedral Au NPs will converge to the extended surface behavior.¹⁰ By relating adsorption energy values to CNs, the group was able to show that for particles with more than 561 atoms (diameter of ca. 2.7 nm), the adsorption energy on a (111)-indexed Au NP facet becomes invariant against the extent of the facet and is close to the number obtained for an extended Au(111) model surface. In another example, Ouyang et al. chose a Au(997) surface to simulate a NP system and used CNs to study O and NO adsorption.⁶ The group identified the relative activity of different adsorbed oxygen species towards NO reduction which they found to be “island O” > “terrace O” > “O near step-edge sites”.

A noteworthy limitation of CNs is their low sensitivity for microscopic features. For example, the ABC-layered structure of a Pt(111) surface gives rise to two different threefold adsorption sites, HCP and FCC. It is well known that these two sites interact differently with certain adspecies such as oxygen; however, since CNs only take into account the first coordination sphere, adatoms adsorbed at HCP and FCC sites are assigned the same CN. Furthermore, when CNs are used with NPs, they do not capture edge or kink site effects. As shown in Fig. 1a, all facet atoms barring those located exactly on the edges and vertices are assigned values of the corresponding single crystal surface.

Fig. 1 Comparison of (a) classical coordination number analysis and (b) generalized coordination number analysis of a 3 nm octahedral Pt NP. Top: Outside view. Bottom: Cross section. Different colors indicate different coordination numbers.

The study by Kleis et al. discussed above established that a certain particle size is required before facet atoms express properties analogous to the corresponding single crystal surface.¹⁰ CNs as a categorization tool are too coarse-grained to capture such effects. This limitation is further exemplified in the cross section view in Fig. 1a. The cross section illustrates that CNs immediately converge to the bulk value of an FCC crystal at the first layer below the surface. This step-like change does not reflect the more gradual – albeit still fast – decay of surface properties with increasing depth.

3.2.2 Generalized coordination numbers. Generalized coordination numbers were introduced by Calle-Vallejo and co-workers in 2014.⁴⁰ The of a specific atom constitutes the average of the CNs of its nearest-neighbor atoms, see eqn (1):Here, CN(j) are the classical coordination numbers of neighboring atoms j and CN_max is a weighing factor. CN_max defines the largest possible coordination number (12 for an FCC crystal, 9 for a BCC crystal, etc.). Contrary to CNs, are distinct for different types of three-fold adsorption sites or bridge positions on ABC-layered surfaces. Similarly, give a more nuanced picture of the different adsorption sites on NPs, see Fig. 1b.

Calle-Vallejo and co-workers have shown the usefulness of in several studies. The group used to study the interplay between adsorption energetics and adsorbate-induced core–shell deformation.⁴¹ They show that the deformation cost is site-independent and can be assumed constant for a specific adspecies. In 2015, the group used to devise a relationship for adsorption properties of oxygen species on transition-metal surfaces that does not require expensive electronic structure calculations.⁴² The relationship they propose relies on knowledge of the valence of the adatom, which can be obtained from electron-counting rules, and the chemical environment of the active site, which they represent via. The approach can be used to predict surface features that satisfy the energetic requirements of each intermediate of a given reaction mechanism, which Calle-Vallejo et al. confirmed against benchmark DFT calculations. The group successfully applied this approach to inform the experimental modification of a Pt(111) surface which produced a material with threefold enhanced ORR activity.⁴³ Tymoczko et al. further used this approach to explain the heightened experimental HER activity of a Cu–Pt(111) near-surface alloy.⁴⁴

Recently, Calle-Vallejo et al. used to clarify why simulation results from single crystal model surfaces are oftentimes not transferable to NPs.⁸ In particular, they note that Pt surfaces contain both convex and concave sites while typical NP systems like cuboctahedra contain only convex sites that are less active towards the ORR. They find that active sites with promote ORR activity, where corresponds to Pt atoms in an extended Pt(111) surface.

are ill-defined in case of alloy catalysts for which not only the number but also the elemental nature of neighboring atoms are important to describe the properties of a specific site. Similarly, cannot be used to study systems for which strain is the dominant factor governing chemical reactivity.

3.2.3 Orbitalwise coordination numbers. In order to overcome the limitations of outlined above, Ma and Xin extend the coordination number approach to orbitals. Orbitalwise coordination numbers (CN^α, where α refers to the orbital α = s, d)⁴⁵ are based on the electronic structure by way of moment characteristics of the projected density of states. Ma and Xin define CN^α in eqn (2) aswhere M^α_2,1 is the second moment from the moments theorem⁴⁶ and t^α,∞_nn is the sum over all two-center hopping integrals of an α = s or d electron to valence orbitals of nearest-neighbor atoms. M^α_2,1 is sensitive to perturbations introduced by the local chemical environment of a reactive site, thus accounting for effects such as material strain.

Ma and Xin tested the CN^α descriptor against for CO adsorption on small Au clusters of 38 to 181 atoms, see Fig. 2.

Fig. 2 Comparison of CO adsorption energy trends on Au clusters and extended surface obtained with DFT and a (a) - or (b) CN^α-based descriptor method. Reprinted with permission from Ma and Xin, Phys. Rev. Lett., 2017, 118, 036101. Copyright 2017 by the American Physical Society.

The CN^α-based model, which in this case relied on s electron hopping integrals (CN^s) exclusively, is found to correlate better with DFT-PBE results (R² = 0.94) than the model based on (R² = 0.75). Wang et al.⁴⁷ used this approach to study larger Au NPs optimized using the embedded-atom method (EAM).⁴⁸ Adsorption energy results on individual sites were evaluated using s-orbital dependent coordination numbers. Wang et al. further included a microkinetics model to predict CO oxidation rates as a function of particle size. They found improved agreement between experiments and the CN^α based model over using classical or .

3.2.4 Adjusted coordination numbers. Even though offer greater resolution than classic CNs, they can still be too coarse-grained for highly disordered systems such as transition metal oxides for which different oxidation states of surface metal atoms, surface defects, and other non-trivial electronic properties complicate the picture. Adjusted coordination numbers (ACNs) were proposed by Fung and co-workers in 2017 as a descriptor for transition metal oxides in particular.⁴⁹ Here, the ACN of an oxygen atom located at the surface of the system is calculated based on its own coordination number, CN_O, and the sum over the coordination numbers of nearest-neighbor metal centers, , as presented in eqn (3):The constant scaling factor λ is related to the partial charge of the surface oxide. Therefore, conducting a partial charge analysis (see Section 3.3) of all relevant solid-state oxide crystal structures is necessary to calibrate the method for each new material. ACNs offer higher resolution than classic CNs or and are therefore advantageous for describing the local chemical environment of low-coordination, surface-bound species such as oxides in more detail. In their publication, Fung et al. calibrated ACNs for bulk Co₃O₄, V₂O₃, NiO, and C–H bond activation barriers.⁴⁹

3.3 Partial charge analysis

3.3.1 Mulliken. Mulliken proposed the first scheme for partial charge analysis of molecules in 1955.^50–53 The Mulliken scheme is based on the theory of representing molecular wavefunctions as linear combinations of atomic orbitals (LCAO) and is used to quantify the bonding character of the molecular orbital of a pair of two atoms. From LCAO theory, the molecular orbital of a pair of atoms can be described as⁵⁴

ψ_i = c_ijϕ_j + c_ikϕ_k, (4)

where ψ_i is the molecular LCAO wavefunction, ϕ_j,k are the wavefunctions of atoms j and k, and c_ij and c_ik are coefficients of the atomic orbitals. Squaring the wavefunction in eqn (4) gives the probability density. Normalization of the probability density gives the expression

1 = c_ij² + c_ik² + c_ijc_ikS_jk, (5)

where S_jk is the overlap integral over the orbitals of atoms j and k. Mulliken defined c_ij² and c_ik² as the atomic-orbital populations and c_ijc_ikS_jk as the overlap population. The latter is <0 for antibonding, 0 for nonbonding, and >0 for bonding molecular orbitals. These quantities can be arranged in matrix form, the population matrix, where atomic-orbital populations constitute the diagonal and overlap populations the off-diagonal elements. Since applying this type of analysis to each pair of atoms in a system can lead to a large amount of data, different ways of post-treatment have been devised. The net population matrix, for example, is the sum of all population matrices for occupied orbitals; here, the diagonal elements give the total charge in each atomic orbital.

Mulliken charges suffer from a major limitation. Computational codes using the LCAO approach represent atoms through basis sets of different makeup and size. Basis sets may include varying numbers of basis functions to represent an orbital. The total atomic charge extracted via Mulliken's approach depends on the basis set size and has no basis set limit.⁵⁵ Mulliken occupancies typically do not sum up to the total number of electrons in the system and can even be negative in rare cases.⁵⁵ Finally, Mulliken analysis cannot be used with plane wave calculations directly.

3.3.2 Hirshfeld. While Mulliken charge analysis applies a basis-set based partitioning scheme, the Hirshfeld approach (1977) is based on separation of the electron density.⁵⁶ Here, the charge density is divided at each point between atom pairs, proportional to their respective contributions based on the free atoms. This strategy led to the byname of “stockholder” partition analysis. The partitioning procedure is based on earlier works by Politzer and Harris, refined by Hirshfeld to ensure that atomic fragments are always well defined.^56,57 Comparison of the molecular electron density at a site r̄, ρ^mol(r̄), to the charge density of the “promolecule”, ρ^pro(r̄), which is constructed from the individual non-interacting atomic densities, results in the deformation density ρ_d(r̄), as shown in eqn (6):⁵⁸ is the sum over the spherically averaged ground state electron densities of all free atoms α, centered at the location of nucleus R̄_α. For an uncharged molecule, the effective charge of an atom q_α can be obtained using eqn (7),Here, the weighing function ω_α(r̄) denotes the relative share that atom α contributes to the density at site r̄.

De Proft and co-workers extensively tested the Hirshfeld scheme on a series of functionalized organic molecules.⁵⁹ They note that while overall trends for charges are in good agreement with measurements, the magnitude of charges and dipoles is frequently underestimated. This conclusion is corroborated by other benchmark studies.^60–62 Furthermore, Hirshfeld charges can be affected unfavourably by atomic fragments exerting electron withdrawal or excess electrons in an atom's vicinity, as often observed by contradictory charges on hydrogen atoms.^58,63 Saha et al. proposed refinements to overcome these limitations.⁵⁸

3.3.3 Natural population analysis. Natural population analysis (NPA) was introduced in the late 1980s by the Weinhold group^55,64,65 as a reaction to the strong basis set dependence of the Mulliken method. NPA uses natural atomic orbitals (NAOs)⁶⁶ which are based on natural orbitals (NOs) as introduced by Löwdin in the 1950s.^67,68 NAOs do not depend on the basis set due to an occupation-dependent weighing factor which makes sure that unoccupied Rydberg states do not significantly impact the results of the method. Thus, unlike NOs, NAOs always remain close in shape to the original single-atom orbitals and electrons are localized close to the atomic centers for the purpose of counting occupations.

Unlike in case of the Mulliken approach, NPA produces only physically meaningful, non-negative occupations and the overall charge will sum up to the total number of electrons in the system. Weinhold and co-workers showed that NPA can significantly improve results over Mulliken and Löwdin analysis by varying the exponent of a 2p basis function in the SCF calculation of H₂ and relating the resulting energy changes to charge results obtained from the different partial charge analysis methods.⁶⁵

The NPA and Mulliken schemes share the limitation of not being directly applicable to plane wave calculations. Furthermore, dipole moments are overestimated significantly using NPA, which is discussed to be a result of the far-reaching tails of (diffuse) basis functions.^62,69

3.3.4 Bader (atoms-in-molecules). Bader charge analysis, also referred to as the atoms-in-molecules (AIM) method, was proposed in 1990^70,71 and has been implemented in various codes over the years.^72–81 Bader analysis differs from the previously described methods in that it uses the electron density, thereby allowing both LCAO and plane waves methods to be used. Some implementations are, moreover, completely independent of the nuclear coordinates.^79–81

Fundamentally, Bader analysis separates the system by creating curved surfaces that run through minima in the charge density. The volume segments enclosed by these surfaces are referred to as Bader regions. The total charge of each Bader region is obtained by integration. Bader charge analysis has been shown to systematically overestimate charges and dipole moments.^60,62,69

3.3.5 Voronoi deformation density. The Voronoi Deformation Density (VDD) was introduced by Guerra and co-workers in 2004.^62,82 Similar to the Bader method, VDD integrates the electron density over atomic domains to calculate the charge. In case of VDD, these atomic domains are defined by Voronoi polyhedra.⁸³ Furthermore, VDD utilizes the deformation density which is defined as the density change going from Hirshfeld-like promolecules to the final electron density. The charge of atom A is then obtained via integration:Results for VDD and Hirshfeld are often similar due to their shared methodological roots. Guerra et al. benchmarked VDD against Hirshfeld, Bader, and NPA, noting that Bader or NPA overestimate charges.⁶² A recent benchmark study by Choudhuri et al. suggest that Hirshfeld analysis tends to underestimate charges which likely affects VDD as well.⁶⁰ Similar to AIM methods, VDD results are consistent against the basis set size.

3.3.6 Charge model 5. Charge Model 5 (CM5) was developed by Marenich and co-workers in 2012⁸⁴ and constitutes the latest iteration in the CMx (x = 1–5 and 4 M) series of methods.^85–90 CM5 is a semiempirical method that involves fitting charges pre-calculated using the Hirshfeld method with parameters from a test set of dipole moments and charges to overcome the limitations of the base method as discussed in Section 3.3.2. CM5 was benchmarked with promising results.^60,84,91 The CM5 method suffers, like any other (semi-) empirical method, from limited transferability to systems outside the training set.

3.3.7 Density derived electrostatic and chemical analysis. Density derived electrostatic and chemical analysis (DDEC) has been developed over the last decade by Manz and co-workers, with the latest iteration of the approach being DDEC6.^61,92–96 The DDEC approach relies on separation of the electronic density in the vein of the VDD and AIM methods. DDEC6 was tested on pure metals, ice crystals, a large protein system, and water and ozone molecules, as well as a set of heteroatom-encapsulated endohedral C₆₀ complexes. In a recent benchmark, Choudhuri and Truhlar compared DDEC5 with their CM5 method as well as with Bader and Hirshfeld charges for various metal dioxides, sulfides, selenides, and metal complex crystal structures.⁶⁰ They point out that for the same metal center, DDEC6 charges tend to be systematically bigger or smaller than CM5 charges depending on the type of ligand. For example, DDEC6 charges were found to be 20–25% bigger than CM5 charges for MO₂, MCO₃, M(CN)₆, and M(NCNH)₂ structures but smaller by about the same amount for MS₂ and MSe₂ structures. They hypothesize that this behavior may be due to sensitivity of DDEC6 towards the diffuseness of the charge density on the metal centers and advocate for the use of partial charge methods that are validated against dipole moments.

3.3.8 Charge equalization methods. All partial charge analysis methods introduced so far require electronic structure calculations to obtain an optimized charge density of the system under investigation. Although electronic structure calculations can be performed for surprisingly large systems using highly optimized and scalable DFT codes, it is still more common to treat NP model systems using force field methods. In the following, two methods are presented that are constructed to work on the basis of classical potential terms.

The electronegativity equalization method (EEM) was first introduced by Mortier et al. (1985)^97–99 based on Sanderson's principle of electronegativity equalization states.¹⁰⁰ With its easily adaptable equations and transferable parameters, the model provides an efficient analysis of charge distributions at a low computational cost. EEM has been validated and successfully applied for applications including organic molecules,^101–109 inorganic solids,^110–116 biomolecular systems,^117–122 metal-organic frameworks,^123–125 and electrochemical systems.^126–128 Further improvements and formalisms evolved from the original EEM principle. Derivative methods include the atom-bond electronegativity equalization (ABEEM),¹²⁹ charge equilibration (QEq and CHEQ),^130,131 fluctuating charges (FlucQ),¹³² chemical potential equalization (CPE),¹³³ and split charge equilibration (SQE) approaches.¹³⁴

Our group recently used EEM partial charge analysis to characterize oxidized cuboctahedral platinum NPs with to identify the prevalent type of oxide in this system which at first glance appears amorphous.²⁴ An exemplary EEM partial charge distribution is shown in Fig. 3. To provide context for the obtained partial charge distribution, 3 nm spherical particles cut from bulk PtO, α-PtO₂, β-PtO₂, and Pt₃O₄ are analyzed and the maxima of the oxygen partial charge distributions are shown as dotted lines in Fig. 3. Comparison to the reference particles indicates that the oxygen atoms in the oxidized octahedral NP most likely correspond to Pt₃O₄.

Fig. 3 Distribution of EEM partial charges on oxygen atoms for a fully oxidized (p_O = 1 mbar, 500 K) 3 nm octahedral Pt NP obtained using the ReaxFF-GCMC simulation approach. Dotted reference lines indicate the oxygen partial charge distribution maxima obtained from spherical 3 nm reference particles cut from bulk oxide materials. The analysis indicates that the oxide composition of the particle is most closely related to a Pt₃O₄ reference particle. Adapted from Kirchhoff et al.²⁴

Notably, EEM-based approaches suffer from cubic scaling of the dipole polarizability with the system size^135,136 and often predict fractional molecular charges for systems with well-separated molecules.¹³⁷ Fractional charges result in a conductor-like charge distribution throughout the system and an ill-defined description of oxidized and reduced species.

Verstraelen et al. addressed these limitations by proposing the atom-condensed Kohn–Sham density functional theory approximated to second-order (ACKS2) scheme.^138,139 Here, the EEM formalism is expanded by additional variables and quadratic energy terms to overcome the unrealistic long-range charge smearing. ACKS2 provides a more realistic description of charge delocalization. All parameters of the ACKS2 model have direct physical meaning and can be obtained from DFT calculations.^138,140

Recent applications of the ACKS2 approach, among others, are the pseudoclassical treatment of explicit electrons in reactive force field simulations,^141–143 proton-transfer reactions in subcritical and supercritical water systems,¹⁴⁴ and charge-transfer processes in electrolytes of Li-ion batteries.^143,145,146

3.4 Radial atomic density distributions

Radial atomic density distributions visualize the distribution of a particular chemical species in a complex heteroatomic material. Senftle and co-workers used it to illustrate different stages of oxidation of a Pd NP,³² hydrogen uptake of a Pd NP,³³ as well as Pd carbide formation.³⁴ To obtain such a distribution, atoms of the respective chemical species are counted in volume segments radially outward from the center of the system. In Fig. 4a, radial atomic density distributions are used to show the distribution of oxygen atoms in barely oxidized (blue), surface-oxidized (orange), and fully oxidized (green) 3 nm cuboctahedral Pt NPs. While oxygen atoms are mostly located in the outer shell of the particle at surface oxidation conditions (blue and orange), oxidation all the way to the core can be observed beyond this point (green). Broadening of the green distribution indicates expansion of the NP as a result of oxidation.

Fig. 4 (a) Radial oxygen distributions of a barely oxidized (blue), surface oxidized (orange), and fully oxidized (green) 3 nm octahedral Pt NPs obtained from ReaxFF GCMC calculations. (b) Integrated, normalized radial oxygen density distributions of the same structures.

Integration over the volume segments and normalization of the results gives the distributions shown in Fig. 4b. Compared to Fig. 4a, this way of presentation is less noisy and allows for large amounts of data to be presented at once. For example, Fig. 4b presents integrated distributions for the three aforementioned barely oxidized (blue), surface-oxidized (orange), and fully oxidized (green) structures. While Fig. 4b still retains information about the oxidation-dependent expansion of the particles, the regular radial atomic density distribution in Fig. 4a can be easier to read depending on what information is most important to visualize; for example, it is immediately obvious from the absolute Frequency value in Fig. 4a that the highest density of O atoms in the configuration obtained at 400 K is found between 12 and 14 Å, whereas this information is contained in the slope of the corresponding data in Fig. 4b.

3.5 Oxidation state analysis

When atomic (partial) charges are not available, for example in force field calculations, the oxidation state of an atom in a metallic NP can be roughly estimated by the number of its metal-metal bonds. A density of oxidation states can thus be obtained by showing the number of atoms within a specific oxidation state (i.e. number of metal-metal bonds) as a function of the radial distance r from the particle center. Hong and van Duin¹⁴⁷ used this technique to monitor the formation of different types of oxides during a ReaxFF-driven reactive molecular dynamics (MD) simulation of an aluminium NP in the presence of oxygen. Fig. 5 shows an application of the method to a strongly oxidized 3 nm octahedral Pt NP.

Fig. 5 A simplified oxidation state analysis of a strongly oxidized 3 nm octahedral Pt NP. The analysis shows that the particle has some core–shell character, with a strongly oxidized exterior (0–74% Pt neighbors, red and green) but some remaining metallic states at the center (100% Pt neighbors, blue). Different oxidation states of Pt atoms are estimated by the ratio between the total number of nearest-neighbor atoms and the number of nearest-neighbor oxygen atoms. Note that in this evaluation, undercoordinated Pt atoms at edges and kink sites with no oxygen neighbors still count as fully metallic (100% Pt neighbors). Depending on the desired information, the ratio may also be calculated with respect to the maximum number of bonds in the bulk crystal to account for the unsaturated nature of edge and kink sites.

3.6 Radial distribution function

A radial distribution function (RDF, g(r)) is obtained by calculating distances between all pairs of atoms in a system and presenting them as a volume-normalized density (ρ = N/V). Usually, the density is obtained by counting all distances that fall into r + dr sphere segments originating from the center of the system. The RDF is therefore entirely dependent on the atomic configuration of the system.

The g(r) can be measured in X-ray^148,149 or neutron scattering experiments.^149–152 Theory can play an integral role in the analysis of scattering data. For example, hybrid quantum mechanical/molecular mechanical (QM/MM) simulations have been successfully used to predict RDFs and assist in the interpretation of the experimental result.^149,153,154

Going back to our guinea pig, the 3 nm octahedral Pt NPs, the g(r) can be used to quantify the degree of amorphicity in a series of increasingly oxidized structures. Fig. 6 illustrates RDFs of an almost clean octahedral particle (blue) as well as of a surface-oxidized (orange) and of a fully oxidized particle (green).

Fig. 6 RDFs of clean (blue), surface-oxidized (orange), and fully oxidized (green) 3 nm octahedral Pt NPs obtained from ReaxFF-GCMC simulations. Significant broadening of the peaks is observed as long-range ordering of the particles decreases as a result of progressing oxidation.

The RDFs become more and more liquid-like as oxidation proceeds and long-range ordering is lost.

A recent example of RDFs applied to NP systems was published by Thomä et al. who experimentally measured the RDF response of the water solvation shell around faceted iron oxide NPs.¹⁵⁵ The group used calculated RDFs of different bonding situations of water to assign the correct structures to the measured bands. In another recent computational study, Zeng et al. used RDFs to monitor the structure and response to heating of a core–shell aluminium/aluminium oxide NP and its surrounding water shell over the course of a ReaxFF reactive MD simulation.¹⁵⁶

3.7 Common neighbor analysis

Common neighbor analysis (CNA) was introduced by Honeycutt and Andersen in 1987.¹⁵⁷ CNA seeks to identify reoccurring structural patterns in atomic configurations. CNA is based on earlier works by Blaisten-Barojas¹⁵⁸ and Haymet¹⁵⁹ and was further developed and first implemented into code by Clarke, Faken, and Jónsson.^160,161 In CNA, each pair of atoms in a structure is assigned a triple of indices, jkl.‡ The three indices represent the number of nearest neighbors common to both atoms (j), the number of bonds between common neighbors (k), and the number of bonds in the longest chain formed by the k bonds between common neighbors (l). These triples are unique for characteristic structural motifs such as FCC, HCP, BCC, or icosahedral arrangements of atoms. Large structures can therefore be characterized based on the occurrence of smaller sub-structures with fingerprint triples. While relatively insensitive to small displacements of atoms, the method will quickly highlight larger changes in a crystal structure caused, for example, by local melting, and can therefore be used to analyze phase transitions. In their original publication, Honeycutt and Andersen annealed small Lennard-Jones clusters of 13 to 309 atoms in MD simulations and used CNA to identify the dominant structural feature (FCC, HCP, or icosahedral) of the relaxed structures.

CNA can be understood as a partitioning of the RDF of the system. The RDF, g(r), of a structure can be represented as a sum of RDFs of each jkl triple, g_jkl(r), so that

This renders g_jkl(r) an interesting tool for the interpretation of the whole g(r). For example, by integrating the first peak in g_jkl(r), the number of pairs N_jkl with jkl can be determined asN is the total number of atoms in the system and r_c is the integration cutoff. By choosing integration limits accordingly, other peaks can be analyzed as well. The method is sensitive to changes to r_c since the cutoff radius determines which atoms are considered nearest-neighbors and, therefore, which triples are found.

As exemplary applications, CNA has been used to describe solid-liquid phase transitions such as melting of small clusters,^162,163 strain and deformation in metallic systems,¹⁶⁴ and properties of glasses.¹⁶⁵ CNA routines can be found in many computational tools, such as the open-source Atomic Simulation Environment (ASE).²⁹

3.8 X-Ray diffraction spectrum

(Powder) X-ray diffraction (XRD) spectra of NPs can be calculated from a given structure using the Debye scattering equationwhere I is the resulting scattering intensity, q is the scattering vector with the incidence angle θ and the incidence wavelength λ, i and j are atomic indices, N is the total number of atoms in the system, f_i and f_j are scattering factors specific to the chemical species of atoms i and j, and r_ij is the distance between atoms i and j. The scattering factors f can be calculated from tabulated values^166,167 using the relationshipwhere a_i, b_i, and c are the tabulated parameters (i = 1–4), θ is the incidence angle and λ is the incidence wave length.

Fig. 7 compares powder XRD spectra of barely oxidized (blue), surface oxidized (orange), and fully oxidized (green) 3 nm octahedral Pt NPs. Band broadening is observed already for the barely oxidized system as a result of its nanometer size.¹⁶⁸ As oxidation proceeds, bands broaden even further and the overall intensity I reduces, indicating that the NP becomes increasingly disordered. Similar to RDFs, with which this approach has in common the dependence on atomic pair distances, the calculated powder XRD spectra can therefore be used to quantify the degree of amorphicity in a series of structures aside from determining the characteristic diffraction peaks.

Fig. 7 Simulated powder XRD spectra of clean (blue), surface-oxidized (orange), and fully oxidized (green) 3 nm octahedral Pt NPs. Assuming Cu-K_α incidence radiation.

Examples for successful use of simulated XRD spectra in literature are abundant. To name a few, Chiche and co-workers showed that calculated XRD spectra can be used in combination with measurements to reveal detailed structural information.¹⁶⁹ By comparing simulated and experimental XRD spectra of MgO and bohemite, the group showed that the simulated spectra contain detailed information about size and shape of the studied 3–7 nm NPs. They note that XRD in such cases outperforms TEM measurements which fail to provide morphological details for particles in this size range. Naicker et al. used simulated XRD spectra to monitor the characteristic rutile, anatase, and brookite bands during MD simulations of TiO₂ NPs to study the surface rearrangement.¹⁷⁰ Vogel and co-workers cross-referenced simulated and measured powder XRD spectra to characterize the reversible uptake of H in Pd NPs.¹⁷¹ Similarly, Gilbert and co-workers used simulated powder XRD spectra in combination with extended X-ray absorption fine structure (EXAFS) and wide-angle X-ray scattering (WAXS) spectroscopy techniques to study the increase in crystallinity in ZnS NPs upon hydration.¹⁷²

3.9 Solvent-accessible surface area and interaction heatmaps

The Solvent-Accessible Surface Area (SASA) method can be used to create heatmaps of the interaction strength of probe atoms or molecules with a nanoparticle structure. The SASA method is rooted in biomolecular applications and computational protein structure prediction^173,174 and was originally proposed by Richmond in 1984.¹⁷⁵ To calculate the SASA, the van der Waals (vdW) radius of each surface atom in a structure is extended by a certain distance, for example by the vdW radius of the probe.¹⁷⁴ By systematically placing probe atoms (e.g. oxygen, hydrogen, or a water molecule) on the extended vdW surface, an interaction heatmap visualizing the strength of interaction of the structure with the probe is obtained.

Fig. 8 shows the interaction heatmap for an oxygen probe atom with an oxidized 3 nm octahedral Pt NP.

Fig. 8 Interaction heatmap for an oxygen probe atom with a fully oxidized 3 nm octahedral Pt NP. Equidistant interaction points on the expanded vdW surface of the NP are used. Red: strong interaction; white: weak interaction.

By assigning interaction points to the respective nearest surface atoms and averaging the interaction energy results assigned to each atom, the interaction strength of each surface atom can be estimated. This approach is exemplified in Fig. 9 for a pristine 3 nm octahedral Pt NP, using an oxygen probe atom. In this specific case, the NP edges are revealed to interact most strongly with the oxygen probe, suggesting that adsorption on these sites is preferred during early stages of the oxidation.^24,25 The SASA heatmap approach therefore constitutes an intuitive way of visualizing surface reactivity even for complex systems.

Fig. 9 Interaction heatmap of an oxygen probe atom with a pristine 3 nm octahedral NP. Interaction strengths for the surface atoms are calculated by averaging over interaction points on the SASA nearest to each surface atom. Interaction strength: green > blue > red.

3.10 Two-phase thermodynamics

The Two-Phase Thermodynamics (2PT) method enables the calculation of the density of states (DoS) of liquids from molecular dynamics trajectories, thus allowing the identification of translational, rotational and intramolecular vibrational modes.^176–178 The 2PT method can aid in the interpretation of complex data obtained using surface infrared or Raman spectroscopy methods, in particular for studying the solid-electrolyte interface.^179,180

In the 2PT scheme, the velocity autocorrelation function (VACF) is calculated via the time-dependent atomic velocities, giving access to the DoS function. The total DoS is then separated into a diffusive gas-like component and a solid-like component. The components are treated separately via hard sphere thermodynamics and quantum statistics, enabling calculation of the entropy and the heat capacity.

The 2PT scheme has been used successfully in literature. In a reactive MD based study, Cheng et al. identified intermediates of the CO₂ reduction reaction by their characteristic bands in the vibrational DoS.¹⁸⁰ The group included explicit solvation and the electrode potential in their simulations. Pascal and co-workers evaluated the translational, rotational and vibrational entropic contributions of confined water using the 2PT method from MD simulations. They studied the wetting behavior of hydrophilic carbon nanotubes and report that water in the nanotube is found to be more stable despite the confinement.¹⁸¹ Persson and co-workers introduced a spatially resolved 3D variant of the 2PT method.¹⁷⁹ Using MD simulations, the group studied local contributions to the solvation enthalpy, entropy, and free energy of hydrophilic and hydrophobic small molecular solutes under the influence of ions and the hydrogen bonding network.

4 Analyzing sets of calculations

4.1 Normalized formation energy

To study the stability of different NP structures as a function of particle size, the total energy per atom, E_atom, can be set in relation with the the NP size to obtain stability trends. Such analysis was performed, for example, by Huang et al. on Pt NPs of various shapes from 20 to 100 nm.¹⁸ Another useful variant of this analysis is to show E_atom as a function of the cubic root of the number of atoms in the system, N^−1/3. Kirchhoff et al. used this type of analysis to study the stability of cuboctahedral, octahedral, spherical, dodecahedral, and cubic Pt NPs in a size range of ca. 1–10 nm.²⁴ The results are reproduced in Fig. 10. Displaying the energy per atom as a function of N^−1/3 linearizes the data. This allows for intuitive comparison of the stability of, in this example, different NP shapes. Differences between the particles are most distinct at small sizes and values converge to the bulk E_atom for the ReaxFF force field that they were obtained with, which is 5.77 eV atom⁻¹.¹⁸² This linearized presentation also reveals slight alterations in the facet composition. For example, small differences in the relative size of (111) and (100) facets of cuboctahedral particles show as deviations from the otherwise linear behavior since the surface-to-volume ratio is not constant anymore. The same observation can be made for slightly trunctated octahedral or dodecahedral particles or when the edges of the cubic particle start to round off, for example as a result of surface reactions that take place.

Fig. 10 Stability of cuboctahedral, octahedral, spherical, dodecahedral, and cubic Pt NPs as a function of particle size obtained from ReaxFF calculations. N: number of atoms in a particle. The key diameter of the spherical particle, d_sphere, is given on top. Bottom: Illustration of NP structures of ca. 3 nm size. Adapted from Kirchhoff et al.²⁴

4.2 Adsorption isotherms and isobars

Adsorption isobars and isotherms can reveal information about the structure of a system in the presence of an adspecies at different thermodynamic conditions. Adsorption isobars or isotherms are constructed by showing the surface coverage as a function of either the partial pressure of the adspecies or the system temperature while keeping constant the other respective other value. For example, Senftle et al. constructed hydrogen adsorption isotherms to explore the uptake of hydrogen in Pd as a function of Pd cluster size.³³

A data set for this type of analysis can be achieved in different ways. One route is to manually generate structures of increasing coverage and diverse arrangements of an adspecies x and to calculate the thermodynamic stability of each structure at various different chemical potentials μ_x(p,T) in order to find the most stable coverage for each condition. Alternatively, schemes like the grand-canonical Monte Carlo (GCMC) algorithm^32–34 can be used to generate diverse data sets of structures at certain μ(p,T) conditions using a stochastic sampling approach.

Returning to the data set of oxidized 3 nm octahedral Pt particles, which was generated using a reactive force field and the GCMC scheme at oxygen chemical potentials μ_O(p,T) corresponding to 200–1200 K at UHV (p_O = 10⁻¹⁰ mbar) and NAP (p_O = 1 mbar),²⁶ adsorption isobars can be constructed as shown in Fig. 11. The coverage convention introduced in Section 3.1 is used. The individual data points shown in Fig. 11 correspond to the structures with the most favorable energy of formation at the corresponding μ_O(p,T).

Fig. 11 Adsorption isobars for a 3 nm octahedral Pt particles oxidized at UHV (p_O = 10⁻¹⁰ mbar) and NAP (p_O = 1 mbar) conditions and 200–1200 K temperature. The degree of oxidation is quantified as the quotient between oxygen and platinum atoms in the system (x_O = N_O/N_Pt). Adapted from Kirchhoff et al.²⁶

Fig. 11 reveals that the oxidation process in case of this particular model system can be broadly divided into three steps: phase 1, surface adsorption at high temperatures marked by a slow increase of the N_O/N_Pt ratio; phase 2, quick increase of N_O/N_Pt marking full oxidation of the particle; and phase 3, decoration of the oxidized particles by additional dioxygen species on the surface.²⁴ This trend is most obvious in case of the more aggressive NAP conditions.

4.3 Adsorbate phase diagrams

The thermodynamic adsorption phase diagram can be regarded as a generalized case of an adsorption isotherm and isobar (Section 4.2). To construct a phase diagram, the energy of formation is displayed as a function of the chemical potential μ(p,T) of the adspecies. The formerly separated quantities of temperature and pressure are thus combined in this case. An exemplary phase diagram for adsorption of oxygen on Pt(111) published by Fantauzzi and co-workers¹⁸³ is given in Fig. 12. This type of analysis can be used to predict the approximate coverage under given thermodynamic conditions and can thus be used to supplement experimental assessments¹⁸⁴ and predict new structures.¹⁸³

Fig. 12 Phase diagram for the adsorption of O on a Pt(111) model system. The data set was obtained using ReaxFF-GCMC simulations at a MC temperature of 650 and 600 K. The colored boxes (a–i) correspond to the lines in the phase diagram. Pt atoms are shown in shades of grey depending on the layer or light blue when buckled (only for Θ_O ≤ 1.13 ML). The phase diagram shows five thermodynamically stable phases: A clean Pt(111), B p(2 × 2) adsorbed atomic oxygen, C a low- and a high-coverage surface oxide phases (1 and 2) and D α-PtO₂ bulk oxide. Reprinted with permission from Fantauzzi et al., Angew. Chem., Int. Ed., 2017, 56, 2594–2598. Copyright 2017 by John Wiley & Sons.

Adsorbate phase diagrams can quickly become complex as the number of variables increases, for example in case of co-adsorption of multiple chemical species. Ferguson and co-workers have therefore developed a visualization tool entitled Surface Phase Explorer (SPE) to facilitate visualization of multidimensional phase diagrams.¹⁸⁵

Finally, phase diagrams can be brought into the realm of electrochemistry using the Extended Ab Initio Thermodynamics (EAITD) approach.^3,186,187 This approach makes use of the computational standard hydrogen electrode and therefore assumes that the reaction H₂ ↔ 2H⁺ + 2e⁻ is in equilibrium (assuming pH 0, 298 K, 1 atm pressure).¹⁸⁸ The SHE can be applied to all electrochemical reactions that involve proton transfer reactions. For our exemplary data set of oxidized 3 nm octahedral Pt NPs, we can construct a thermodynamic cycle based on the assumption that all oxygen atoms must come from the water splitting reaction, H₂O ↔ H₂ + ½O₂, to introduce H₂ into the reaction. The potential-dependent energy of formation for each oxidized particle, ΔE^system_F(Δϕ), can then be calculated as

ΔE^system_F(Δϕ) = E^system_tot − E^ref_tot − N_O(μ_O + 2eΔϕ) with (14)

μ_O = μ_H2O − μ_H2, (15)

where E^system_tot is the computed total energy of the oxidized system, E^ref_tot is the total energy of a reference system (here: the pristine NP without adatoms), N_O is the number of oxygen adatoms in the oxidized system, μ_O, μ_H2O, and μ_H2 are the chemical potentials of oxygen, water, and hydrogen, respectively, e is the elementary charge, and Δϕ is the electrode potential. The resulting electrochemical phase diagram, which shows the thermodynamically most stable oxide structures at SHE conditions, is given in Fig. 13. The coverage convention introduced in Section 3.1 is used to classify the oxidized structures.

Fig. 13 Electrochemical phase diagram depicting the most stable adsorbate structures on a 3 nm octahedral Pt NP. Adsorbate structures are characterized by their O : Pt ratio, x_O, as outlined in Section 3.1. Adapted from Kirchhoff et al.²⁶

5 Conclusions

The present work reviewed analytical methods useful for extracting thermodynamic information from computational simulations of NP systems. As computational methods continue to improve and costs for high-performance computing are brought down by the year, simulations of large NP model systems instead of stand-ins such as periodic surface models will become the norm. Many established analytical methods, such as partial charge analysis or radial distribution functions, can easily be scaled for application to NP systems. On the other hand, the more complex nature of NPs requires new methods and approaches to analyze structure-activity relationships in particular. The authors hope that the present work will serve as a starting point for computational researchers working with NP structures.

Author contributions

BK: conceptualization, methodology, data curation, visualization, writing – original draft. CJ: writing – review & editing. DG: writing – review & editing. DF: conceptualization, supervision. HJ: supervision. TJ: writing – review & editing, supervision, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

BK acknowledges funding from the Professional Graduate and Training Center Ulm (ProTrainU). Further, the authors gratefully acknowledge support by the DFG (german science foundation) through the collaborative research centers SPP-2248 (project ID 441209207) as well as SFB-1279 (project ID 316249678), SFB-1316 (project ID 327886311), TRR-234 (364549901). The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the DFG through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster).

References

1. J. C. Meier, C. Galeano, I. Katsounaros, J. Witte, H. J. Bongard, A. A. Topalov, C. Baldizzone, S. Mezzavilla, F. Schüth and K. J. J. Mayrhofer, Beilstein J. Nanotechnol., 2014, 44–67.
2. N. Tian, Z.-Y. Zhou and S.-G. Sun, J. Phys. Chem. C, 2008, 112, 19801–19817.
3. P. Kaghazchi, F. C. Simeone, K. A. Soliman, L. A. Kibler and T. Jacob, Faraday Discuss., 2009, 140, 69–80.
4. G. A. Tritsaris, J. Greeley, J. Rossmeisl and J. K. Nørskov, Catal. Lett., 2011, 141, 909–913.
5. T. Zhu, E. J. M. Hensen, R. A. van Santen, N. Tian, S.-G. Sun, P. Kaghazchi and T. Jacob, Phys. Chem. Chem. Phys., 2013, 15, 2268.
6. G. Ouyang, K.-J. Zhu, L. Zhang, P.-F. Cui, B.-T. Teng and X.-D. Wen, Appl. Surf. Sci., 2016, 387, 875–881.
7. I. V. Yudanov, R. Sahnoun, K. M. Neyman and N. Rösch, J. Chem. Phys., 2002, 117, 9887–9896.
8. F. Calle-Vallejo, M. D. Pohl, D. Reinisch, D. Loffreda, P. Sautet and A. S. Bandarenka, Chem. Sci., 2017, 8, 2283–2289.
9. Y. Volokitin, J. Sinzig, L. J. de Jongh, G. Schmid, M. N. Vargaftik and I. I. Moiseevi, Nature, 1996, 384, 621–623.
10. J. Kleis, J. Greeley, N. A. Romero, V. A. Morozov, H. Falsig, A. H. Larsen, J. Lu, J. J. Mortensen, M. Dułak, K. S. Thygesen, J. K. Nørskov and K. W. Jacobsen, Catal. Lett., 2011, 141, 1067–1071.
11. F. Viñes, J. R. B. Gomes and F. Illas, Chem. Soc. Rev., 2014, 43, 4922–4939.
12. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
13. Y. Han, D. Jiang, J. Zhang, W. Li, Z. Gan and J. Gu, Front. Chem. Sci. Eng., 2016, 10, 16–38.
14. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009.
15. A. P. Sutton and J. Chen, Philos. Mag. Lett., 1990, 61, 139–146.
16. Y. Qi, T. Çağin, W. L. Johnson and W. A. Goddard, J. Chem. Phys., 2001, 115, 385–394.
17. N. Tian, Z.-Y. Zhou and S.-G. Sun, J. Phys. Chem. C, 2008, 112, 19801–19817.
18. R. Huang, Y.-H. Wen, Z.-Z. Zhu and S.-G. Sun, J. Mater. Chem., 2011, 21, 11578–11584.
19. A. Warshel and M. Levitt, J. Mol. Biol., 1976, 103, 227–249.
20. B. T. Thole and P. T. van Duijnen, Theor. Chim. Acta, 1980, 55, 307–318.
21. M. J. Field, P. A. Bash and M. Karplus, J. Comput. Chem., 1990, 11, 700–733.
22. X. Lu, D. Fang, S. Ito, Y. Okamoto, V. Ovchinnikov and Q. Cui, Mol. Simul., 2016, 42, 1056–1078.
23. A. O. Dohn, Int. J. Quantum Chem., 2020, 120, e26343.
24. B. Kirchhoff, L. Braunwarth, C. Jung, H. Jónsson, D. Fantauzzi and T. Jacob, Small, 2019, 16, 1905159.
25. L. Braunwarth, C. Jung and T. Jacob, Top. Catal., 2020, 1647–1657.
26. B. Kirchhoff, C. Jung, H. Jónsson, D. Fantauzzi and T. Jacob, J. Phys. Chem. C, 2022, 126, 6773–6781.
27. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
28. J. Stone, MSc thesis, Computer Science Department, University of Missouri-Rolla, 1998.
29. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys. Cond. Matter, 2017, 29, 273002.
30. J. Rahm and P. Erhart, J. Open Source Softw., 2020, 5, 1944.
31. F. Uekermann, S. Fiedler and B. Aradi, nanocut, version 12.12; 2019, https://bitbucket.org/aradi/nanocut, accessed March 3, 2022.
32. T. P. Senftle, R. J. Meyer, M. J. Janik and A. C. T. V. Duin, J. Chem. Phys., 2013, 139, 044109.
33. T. P. Senftle, M. J. Janik and A. C. T. van Duin, J. Phys. Chem. C, 2014, 118, 4967–4981.
34. T. P. Senftle, A. C. T. van Duin and M. J. Janik, Catal. Commun., 2014, 52, 72–77.
35. I. J. Robertson, M. C. Payne and V. Heine, Europhys. Lett., 1991, 15, 301–306.
36. B. Hammer and J. K. Nørskov, Surf. Sci., 1995, 343, 211–220.
37. T. Jiang, D. J. Mowbray, S. Dobrin, H. Falsig, B. Hvolbæk, T. Bligaard and J. K. Nørskov, J. Phys. Chem. C, 2009, 113, 10548–10553.
38. G. Mpourmpakis, A. N. Andriotis and D. G. Vlachos, Nano Lett., 2010, 10, 1041–1045.
39. M. G. Taylor, N. Austin, C. E. Gounaris and G. Mpourmpakis, ACS Catal., 2015, 5, 6296–6301.
40. F. Calle-Vallejo, J. I. Martínez, J. M. García-Lastra, P. Sautet and D. Loffreda, Angew. Chem., Int. Ed., 2014, 53, 8316–8319.
41. F. Calle-Vallejo, P. Sautet and D. Loffreda, J. Phys. Chem. Lett., 2014, 5, 3120–3124.
42. F. Calle-Vallejo, D. Loffreda, M. T. M. Koper and P. Sautet, Nat. Chem., 2015, 7, 403–410.
43. F. Calle-Vallejo, J. Tymoczko, V. Colic, Q. H. Vu, M. D. Pohl, K. Morgenstern, D. Loffreda, P. Sautet, W. Schuhmann and A. S. Bandarenka, Science, 2015, 350, 185–189.
44. J. Tymoczko, F. Calle-Vallejo, W. Schuhmann and A. S. Bandarenka, Nat. Commun., 2016, 7, 10990.
45. X. Ma and H. Xin, Phys. Rev. Lett., 2017, 118, 036101.
46. J. P. Gaspard and F. Cyrot-Lackmann, J. Phys. C: Solid State Phys., 1973, 6, 3077–3096.
47. S. Wang, N. Omidvar, E. Marx and H. Xin, Phys. Chem. Chem. Phys., 2018, 20, 6055–6059.
48. M. S. Daw and M. I. Baskes, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 6443–6453.
49. V. Fung, F. F. Tao and D.-e Jiang, J. Phys. Chem. Lett., 2017, 8, 2206–2211.
50. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
51. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1841–1846.
52. R. S. Mulliken, J. Chem. Phys., 1955, 23, 2338–2342.
53. R. S. Mulliken, J. Chem. Phys., 1955, 23, 2343–2346.
54. T. J. Zielinski, E. Harvey, R. Sweeney and D. M. Hanson, J. Chem. Ed., 2005, 82, 1880.
55. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
56. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
57. P. Politzer and R. R. Harris, J. Am. Chem. Soc., 1970, 92, 6451–6454.
58. S. Saha, R. K. Roy and P. W. Ayers, Int. J. Quantum Chem., 2009, 109, 1790–1806.
59. F. D. Proft, C. V. Alsenoy, A. Peeters, W. Langenaeker and P. Geerlings, J. Comput. Chem., 2002, 23, 1198–1209.
60. I. Choudhuri and D. G. Truhlar, J. Chem. Theory Comput., 2020, 16, 5884–5892.
61. T. A. Manz and N. G. Limas, RSC Adv., 2016, 6, 47771–47801.
62. C. F. Guerra, J.-W. Handgraaf, E. J. Baerends and F. M. Bickelhaupt, J. Comput. Chem., 2004, 25, 189–210.
63. E. R. Davidson and S. Chakravorty, Theor. Chim. Acta, 1992, 83, 319–330.
64. A. E. Reed and F. Weinhold, J. Chem. Phys., 1985, 83, 1736–1740.
65. A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899–926.
66. A. E. Reed and F. Weinhold, J. Chem. Phys., 1983, 78, 4066–4073.
67. P.-O. Löwdin, Phys. Rev., 1955, 97, 1474–1489.
68. P.-O. Löwdin, Phys. Rev., 1955, 97, 1490–1508.
69. F. Martin and H. Zipse, J. Comput. Chem., 2005, 26, 97–105.
70. R. F. W. Bader, Atoms in molecules: a quantum theory, Clarendon Press, Oxford; New York, 1990.
71. R. F. W. Bader and C. F. Matta, J. Phys. Chem. A, 2004, 108, 8385–8394.
72. P. L. A. Popelier, Theor. Chim. Acta, 1994, 87, 465–476.
73. B. B. Stefanov and J. Cioslowski, J. Comput. Chem., 1995, 16, 1394–1404.
74. P. L. A. Popelier, Comput. Phys. Commun., 1998, 108, 180–190.
75. B. Silvi and C. Gatti, J. Phys. Chem. A, 2000, 104, 947–953.
76. P. L. A. Popelier, Theor. Chem. Acc., 2001, 105, 393–399.
77. N. O. J. Malcolm and P. L. A. Popelier, J. Comput. Chem., 2003, 24, 437–442.
78. N. O. J. Malcolm and P. L. A. Popelier, J. Comput. Chem., 2003, 24, 1276–1282.
79. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
80. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, J. Comput. Chem., 2007, 28, 899–908.
81. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204.
82. F. M. Bickelhaupt, N. J. R. van Eikema Hommes, C. Fonseca Guerra and E. J. Baerends, Organometallics, 1996, 15, 2923–2931.
83. G. Voronoi, J. Reine Angew. Math., 1908, 1908, 198–287.
84. A. V. Marenich, S. V. Jerome, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2012, 8, 527–541.
85. J. W. Storer, D. J. Giesen, C. J. Cramer and D. G. Truhlar, J. Comput.-Aided Mol. Des., 1995, 9, 87–110.
86. J. Li, T. Zhu, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. A, 1998, 102, 1820–1831.
87. P. Winget, J. D. Thompson, J. D. Xidos, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. A, 2002, 106, 10707–10717.
88. J. M. Brom, B. J. Schmitz, J. D. Thompson, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. A, 2003, 107, 6483–6488.
89. C. P. Kelly, C. J. Cramer and D. G. Truhlar, Theor. Chem. Acc., 2005, 113, 133–151.
90. R. M. Olson, A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2007, 3, 2046–2054.
91. J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, J. Chem. Theory Comput., 2014, 10, 2802–2812.
92. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2010, 6, 2455–2468.
93. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2012, 8, 2844–2867.
94. N. G. Limas and T. A. Manz, RSC Adv., 2016, 6, 45727–45747.
95. T. A. Manz, RSC Adv., 2017, 7, 45552–45581.
96. N. G. Limas and T. A. Manz, RSC Adv., 2018, 8, 2678–2707.
97. W. J. Mortier, K. Van Genechten and J. Gasteiger, J. Am. Chem. Soc., 1985, 107, 829–835.
98. W. J. Mortier, S. K. Ghosh and S. Shankar, J. Am. Chem. Soc., 1986, 108, 4315–4320.
99. K. A. Van Genechten, W. J. Mortier and P. Geerlings, J. Chem. Phys., 1987, 86, 5063–5071.
100. R. T. Sanderson, Science, 1951, 114, 670–672.
101. T. Verstraelen, V. Van Speybroeck and M. Waroquier, J. Chem. Phys., 2009, 131, 044127.
102. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
103. W. J. Mortier, S. K. Ghosh and S. Shankar, J. Am. Chem. Soc., 1986, 108, 4315–4320.
104. Z.-Z. Yang and C.-S. Wang, J. Phys. Chem. A, 1997, 101, 6315–6321.
105. P. Bultinck, R. Vanholme, P. L. A. Popelier, F. De Proft and P. Geerlings, J. Phys. Chem. A, 2004, 108, 10359–10366.
106. P. Bultinck, W. Langenaeker, P. Lahorte, F. De Proft, P. Geerlings, M. Waroquier and J. P. Tollenaere, J. Phys. Chem. A, 2002, 106, 7887–7894.
107. Y. Zhong, G. L. Warren and S. Patel, J. Comput. Chem., 2008, 29, 1142–1152.
108. J. E. Davis, G. L. Warren and S. Patel, J. Phys. Chem. B, 2008, 112, 8298–8310.
109. R. S. Varekova, Z. Jirouskova, J. Vanek, S. Suchomel and J. Koca, Int. J. Mol. Sci., 2007, 8, 572–582.
110. K. A. Van Genechten, W. J. Mortier and P. Geerlings, J. Chem. Phys., 1987, 86, 5063–5071.
111. G. O. A. Janssens, B. G. Baekelandt, H. Toufar, W. J. Mortier and R. A. Schoonheydt, J. Phys. Chem., 1995, 99, 3251–3258.
112. A. Ostadhossein, A. Rahnamoun, Y. Wang, P. Zhao, S. Zhang, V. H. Crespi and A. C. T. van Duin, J. Phys. Chem. Lett., 2017, 8, 631–640.
113. D. Akbarian, D. E. Yilmaz, Y. Cao, P. Ganesh, I. Dabo, J. Munro, R. Van Ginhoven and A. C. T. van Duin, Phys. Chem. Chem. Phys., 2019, 21, 18240–18249.
114. R. Heidler, G. O. Janssens, W. J. Mortier and R. A. Schoonheydt, Microporous Mater., 1997, 12, 1–11.
115. K. Smirnov and D. Bougeard, Chem. Phys., 2003, 292, 53–70.
116. J. Wen, T. Ma, W. Zhang, G. Psofogiannakis, A. C. van Duin, L. Chen, L. Qian, Y. Hu and X. Lu, Appl. Surf. Sci., 2016, 390, 216–223.
117. C.-M. Ionescu, S. Geidl, R. Svobodová Vařeková and J. Koča, J. Chem. Inf. Model., 2013, 53, 2548–2558.
118. Y. Ouyang, F. Ye and Y. Liang, Phys. Chem. Chem. Phys., 2009, 11, 6082–6089.
119. A. Baeten and P. Geerlings, THEOCHEM, 1999, 465, 203–207.
120. J. L. Banks, G. A. Kaminski, R. Zhou, D. T. Mainz, B. J. Berne and R. A. Friesner, J. Chem. Phys., 1999, 110, 741–754.
121. T. Verstraelen, E. Pauwels, F. De Proft, V. Van Speybroeck, P. Geerlings and M. Waroquier, J. Chem. Theory Comput., 2012, 8, 661–676.
122. S. Patel, A. D. Mackerell Jr. and C. L. Brooks III, J. Comput. Chem., 2004, 25, 1504–1514.
123. C. E. Wilmer and R. Q. Snurr, Chem. Eng. J., 2011, 171, 775–781.
124. L. Huang, K. L. Joshi, A. C. T. V. Duin, T. J. Bandosz and K. E. Gubbins, Phys. Chem. Chem. Phys., 2012, 14, 11327–11332.
125. L. Huang, T. Bandosz, K. L. Joshi, A. C. T. van Duin and K. E. Gubbins, J. Chem. Phys., 2013, 138, 034102.
126. A. Ostadhossein, S.-Y. Kim, E. D. Cubuk, Y. Qi and A. C. T. van Duin, J. Phys. Chem. A, 2016, 120, 2114–2127.
127. Y. K. Shin, M. Y. Sengul, A. S. M. Jonayat, W. Lee, E. D. Gomez, C. A. Randall and A. C. T. V. Duin, Phys. Chem. Chem. Phys., 2018, 20, 22134–22147.
128. M. Raju, P. Ganesh, P. R. C. Kent and A. C. T. van Duin, J. Chem. Theory Comput., 2015, 11, 2156–2166.
129. Z.-Z. Yang and C.-S. Wang, J. Phys. Chem. A, 1997, 101, 6315–6321.
130. A. K. Rappe and W. A. Goddard, J. Phys. Chem., 1991, 95, 3358–3363.
131. Y. Zhong and S. Patel, J. Phys. Chem. B, 2010, 114, 11076–11092.
132. S. W. Rick, S. J. Stuart and B. J. Berne, J. Chem. Phys., 1994, 101, 6141–6156.
133. D. M. York and W. Yang, J. Chem. Phys., 1996, 104, 159–172.
134. R. A. Nistor, J. G. Polihronov, M. H. Müser and N. J. Mosey, J. Chem. Phys., 2006, 125, 094108.
135. G. Lee Warren, J. E. Davis and S. Patel, J. Chem. Phys., 2008, 128, 144110.
136. R. A. Nistor and M. H. Müser, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 104303.
137. J. Cioslowski and B. B. Stefanov, J. Chem. Phys., 1993, 99, 5151–5162.
138. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, J. Chem. Phys., 2013, 138, 074108.
139. T. Verstraelen, S. Vandenbrande and P. W. Ayers, J. Chem. Phys., 2014, 141, 194114.
140. P. Gütlein, L. Lang, K. Reuter, J. Blumberger and H. Oberhofer, J. Chem. Theory Comput., 2019, 15, 4516–4525.
141. M. M. Islam, G. Kolesov, T. Verstraelen, E. Kaxiras and A. C. T. van Duin, J. Chem. Theory Comput., 2016, 12, 3463–3472.
142. B. Evangelisti, K. A. Fichthorn and A. C. T. van Duin, J. Chem. Phys., 2020, 153, 104106.
143. M. M. Islam and A. C. T. van Duin, J. Phys. Chem. C, 2016, 120, 27128–27134.
144. H. Manzano, W. Zhang, M. Raju, J. S. Dolado, I. López-Arbeloa and A. C. T. van Duin, J. Chem. Phys., 2018, 148, 234503.
145. M. J. Hossain, G. Pawar, B. Liaw, K. L. Gering, E. J. Dufek and A. C. T. van Duin, J. Chem. Phys., 2020, 152, 184301.
146. K. A. O’Hearn, M. W. Swift, J. Liu, I. Magoulas, P. Piecuch, A. C. T. van Duin, H. M. Aktulga and Y. Qi, J. Chem. Phys., 2020, 153, 084107.
147. S. Hong and A. C. van Duin, J. Phys. Chem. C, 2015, 119, 17876–17886.
148. E. B. Sirota, H. D. Ou-Yang, S. K. Sinha, P. M. Chaikin, J. D. Axe and Y. Fujii, Phys. Rev. Lett., 1989, 62, 1524–1527.
149. A. K. Soper, J. Phys.: Condens. Matter, 2007, 19, 335206.
150. N. S. Gingrich and L. Heaton, J. Chem. Phys., 1961, 34, 873–878.
151. J. L. Yarnell, M. J. Katz, R. G. Wenzel and S. H. Koenig, Phys. Rev. A: At., Mol., Opt. Phys., 1973, 7, 2130–2144.
152. F. Li and J. S. Lannin, Phys. Rev. Lett., 1990, 65, 1905–1908.
153. A. O. Dohn, E. Biasin, K. Haldrup, M. M. Nielsen, N. E. Henriksen and K. B. Møller, J. Phys. B: At., Mol. Opt. Phys., 2015, 48, 244010.
154. A. O. Dohn, E. Biasin, K. Haldrup, M. M. Nielsen, N. E. Henriksen and K. B. Møller, J. Phys. B: At., Mol. Opt. Phys., 2016, 49, 059501.
155. S. L. J. Thomä, S. W. Krauss, M. Eckardt, P. Chater and M. Zobel, Nat. Commun., 2019, 10, 995.
156. H. Zeng, X. Cheng, C. Zhang and Z. Lu, J. Phys. Chem. C, 2018, 122, 9191–9197.
157. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem., 1987, 91, 4950–4963.
158. E. Blaisten-Barojas, Kinam, 1984, 6A, 71–84.
159. A. D. J. Haymet, Chem. Phys. Lett., 1984, 107, 77–82.
160. A. S. Clarke and H. Jónsson, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 3975–3984.
161. D. Faken and H. Jónsson, Comput. Mater. Sci., 1994, 2, 279–286.
162. C. Cleveland, W. Luedtke and U. Landman, Phys. Rev. Lett., 1998, 81, 2036–2039.
163. C. L. Cleveland, W. D. Luedtke and U. Landman, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 5065–5077.
164. M. R. Sørensen, M. Brandbyge and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 3283–3294.
165. Y. Q. Cheng and E. Ma, Prog. Mater. Sci., 2011, 56, 379–473.
166. D. T. Cromer and J. B. Mann, Acta Crystallogr., 1968, 24, 321–324.
167. P. J. Brown, A. G. Fox, E. N. Maslen, M. A. O’Keefe and B. T. M. Willis, International Tables for Crystallography, American Cancer Society, 2006, pp. 554–595.
168. C. F. Holder and R. E. Schaak, ACS Nano, 2019, 13, 7359–7365.
169. D. Chiche, M. Digne, R. Revel, C. Chanéac and J.-P. Jolivet, J. Phys. Chem. C, 2008, 112, 8524–8533.
170. P. K. Naicker, P. T. Cummings, H. Zhang and J. F. Banfield, J. Phys. Chem. B, 2005, 109, 15243–15249.
171. W. Vogel, W. He, Q.-H. Huang, Z. Zou, X.-G. Zhang and H. Yang, Int. J. Hydrogen Energy, 2010, 35, 8609–8620.
172. B. Gilbert, H. Zhang, F. Huang, J. F. Banfield, Y. Ren, D. Haskel, J. C. Lang, G. Srajer, A. Jürgensen and G. A. Waychunas, J. Chem. Phys., 2004, 120, 11785–11795.
173. E. Durham, B. Dorr, N. Woetzel, R. Staritzbichler and J. Meiler, J. Mol. Model., 2009, 15, 1093–1108.
174. B. Lee and F. Richards, J. Mol. Biol., 1971, 55, 379-IN4.
175. T. J. Richmond, J. Mol. Biol., 1984, 178, 63–89.
176. S.-T. Lin, M. Blanco and W. A. Goddard, J. Chem. Phys., 2003, 119, 11792–11805.
177. S.-T. Lin, P. K. Maiti and W. A. Goddard, J. Phys. Chem. B, 2010, 114, 8191–8198.
178. T. A. Pascal, S.-T. Lin and W. A. Goddard III, Phys. Chem. Chem. Phys., 2011, 13, 169–181.
179. R. A. X. Persson, V. Pattni, A. Singh, S. M. Kast and M. Heyden, J. Chem. Theory Comput., 2017, 13, 4467–4481.
180. T. Cheng, A. Fortunelli and W. A. Goddard, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 7718–7722.
181. T. A. Pascal, W. A. Goddard and Y. Jung, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 11794–11798.
182. D. Fantauzzi, J. Bandlow, L. Sabo, J. E. Mueller, A. C. T. V. Duin and T. Jacob, Phys. Chem. Chem. Phys., 2014, 16, 23118–23133.
183. D. Fantauzzi, S. Krick-Calderón, J. E. Mueller, M. Grabau, C. Papp, H.-P. Steinrück, T. P. Senftle, A. C. T. van Duin and T. Jacob, Angew. Chem., Int. Ed., 2017, 56, 2594–2598.
184. D. Fantauzzi, J. E. Mueller, L. Sabo, A. C. T. V. Duin and T. Jacob, ChemPhysChem, 2015, 16, 2797–2802.
185. G. A. Ferguson, V. Vorotnikov, N. Wunder, J. Clark, K. Gruchalla, T. Bartholomew, D. J. Robichaud and G. T. Beckham, J. Phys. Chem. C, 2016, 120, 26249–26258.
186. T. Jacob, J. Electroanal. Chem., 2007, 607, 158–166.
187. S. Venkatachalam and T. Jacob, Density functional theory applied to electrocatalysis, John Wiley & Sons, Ltd., 2010.
188. S. Trasatti, Pure Appl. Chem., 1986, 58, 955–966.

---

Footnotes

† An online, Jupyter book based tutorial is available which outlines how to implement many of the methods presented here: https://bjk24.gitlab.io/in-silico-review/intro.html. The code is built to interact with a recently published data set of oxidized NP structures (DOI: https://doi.org/10.5281/zenodo.6322004).²⁶ All figures in this manuscript involving the oxidized 3 nm octahedral NPs were obtained using these scripts and data set. We recommend interested readers to execute and modify the provided code examples to deepen their understanding of these various evaluation methods.

‡ Note that use of quadruples is proposed in the original publication by Honeycutt and Andersen.¹⁵⁷ We focus on the more commonly used triples representation introduced by Clarke and Jónsson in 1993.¹⁶⁰

---