Global Optimization of Copper Clusters at the ZnO(10-10) Surface Using a DFT-based Neural Network Potential and Genetic Algorithms

I Introduction

NOTE: The following article has been submitted to The Journal of Chemical Physics. After it is published, it will be found at https://aip.scitation.org/journal/jcp.

Supported clusters, with the cluster and the support being either metals or oxides, are an active area of research Gates (1995); Santra and Goodman (2002); Wang (2004); Heiz and Bullock (2004); Yu, Porosoff, and Chen (2012) in surface science and catalysis. Accordingly, they have been studied extensively both with theoretical Cabria, López, and Alonso (2010); Robles and Khanna (2010); Tang, Yang, and Dai (2012) and experimental Heiz et al. (1997); Haas et al. (2000); Yamaguchi and Iglesia (2010); Peters et al. (2013) methods. A prominent example for such catalysts that is used in the industrial synthesis of methanol Bowker et al. (1988); Saito et al. (1996); Sá et al. (2010), consists of ZnO and Cu nanoparticles. Given the commercial importance of methanol as a solvent and as a reagent, this system has received a lot of attention in recent years, and many theoretical French et al. (2008); Cheng et al. (2014); Artrith, Hiller, and Behler (2013); Hellström and Behler (2016); Mora-Fonz et al. (2017) as well as experimental Dulub, Boatner, and Diebold (2002); Koplitz, Dulub, and Diebold (2003); Kroll and Köhler (2007); Behrens et al. (2012); Kuld et al. (2016); Kattel et al. (2017) studies are available with the aim of unraveling the structural details of this system.

Computer simulations can provide valuable information, complementary to experimental data, but theoretical studies of supported clusters are hampered by several challenges. The first one concerns the system size, as supported clusters can be very large containing hundreds or even thousands of atoms, and large supercells are needed for realistic structural models of the support to prevent spurious interactions between adsorbed clusters and other finite-size effects. This rules out the use of demanding electronic structure methods like density-functional theory (DFT) for all but the most simple systems, as most simulation techniques like molecular dynamics (MD) or Monte Carlo (MC) require the computation of thousands to millions of configurations.

Simple and thus computationally cheap empirical potentials would allow to study a large number of structures, but reliable atomistic potentials for supported clusters, which provide a direct functional relation between the atomic positions and the potential energy, are difficult to construct due to the complex and diverse bonding and interactions present in these systems, and consequently they are available only for a few systems Zhou et al. (2004); van Duin et al. (2010). Thus the second challenge is the lack of reliable potentials, which are able to cover the wide variety of atomic environments and bonding patterns present in the nanoparticle and the support.

A possible solution to this dilemma is offered by machine learning (ML) potentials Handley and Popelier (2010); Behler (2016); Botu et al. (2017), which can provide very accurate potential energy surfaces (PES) based on reference data sets for representative configurations obtained in electronic structure calculations. Once constructed, ML potentials allow to compute the energies and forces of large systems in a very efficient way enabling the study of realistic structural models, while the first-principles accuracy of the underlying reference method is essentially maintained. Starting with the seminal work of Doren and coworkers in 1995 Blank et al. (1995), the technology of ML potentials has made rapid progress in recent years, and many algorithms are available, like neural network-based approaches Lorenz, Groß, and Scheffler (2004); Behler, Lorenz, and Reuter (2007); Schütt et al. (2018); Unke and Meuwly (2019); Jiang and Guo (2013); Smith, Isayev, and Roitberg (2017); Ghasemi et al. (2015), Gaussian approximation potentials Bartók et al. (2010), moment tensor potentials Shapeev (2016), spectral neighbor analysis potentials Thompson et al. (2015), and many others. Here, we utilize high-dimensional neural network potentials (NNPs) proposed by Behler and Parrinello in 2007 Behler and Parrinello (2007); Behler (2015, 2017) to study small copper clusters supported on the (10$\bar{1}$0) surface of ZnO. High-dimensional NNPs of this type have been used successfully in the past for a series of related systems like copper Artrith and Behler (2012), brass nanoparticles Weinreich et al. (2020), zinc oxide Artrith, Morawietz, and Behler (2011), water at copper Kondati Natarajan and Behler (2017a) and zinc oxide surfaces Quaranta, Hellstrom, and Behler (2017), and also for the copper-zinc oxide system Artrith, Hiller, and Behler (2013).

Specifically, here we use NNPs to perform genetic algorithm (GA) Michalewicz and Janikow (1991); Deaven and Ho (1995) based global optimization (GO) searches of pure copper clusters containing between four and ten atoms supported on the low index zinc oxide (10$\bar{1}$0) surface. Other recent papers Kolsbjerg, Peterson, and Hammer (2018); Jennings et al. (2019) have shown the usefulness of such a combined NNP+GA approach. Combining both algorithms allows for generating a new potential for a target system “on-the-fly”, while at the same time requiring a minimal amount of electronic structure calculations. The generated NNP is then used to run GA searches, and good candidates from this GA search are used to further refine the potential. This procedure could, in the future, easily be adapted to run automated global optimization searches on a wide range of materials.

Previous GO searches of this or similar systems utilized either parametrized classical force fields Mora-Fonz et al. (2017), or relied on ab-initio calculations, which requires them to pre-optimize the clusters in an isolated configuration in vacuum Wan et al. (2018) before depositing them on the surface. The NNP allows us to easily include the influence of the support on the optimized clusters in the GO search, without needing to resort to isolated clusters. Additionally, it allows us to run much longer GA searches which provides for better sampling.

The paper is organized as follows: Section II summarizes the key concepts of high-dimensional NNPs and GA searches. In Section III we describe our computational setup for the reference DFT calculations, the procedure for generating the NNP data set, the parameters for the NNP construction and the details of the GA search. In Section IV we present our results addressing the structural properties of the clean ZnO support (IV.1), the global minimum and low-energy local minima from Cu4 to Cu10 (IV.2), the interface structures (IV.3), the geometric properties of the optimized clusters (IV.4) and the role of a possibly frozen support surface (IV.5). Finally, we draw our conclusions in Section V.

II Methods

Neural networks (NN) are one of the main methods in use in the realm of machine learning. They are capable of reproducing any well-behaved function Hornik (1991), such as the potential energy surface (PES) of a system of atoms, if a set of points where the function has been evaluated is available. In the present work we utilize fully connected, high-dimensional, feed-forward neural networks introduced by Behler and Parrinello Behler and Parrinello (2007); Behler (2017), which are constructed to reproduce reference data sets containing the energies and forces obtained from ab-initio calculations for a diverse set of structures.

In the Behler-Parrinello approach, the system for which the energy and forces are to be predicted is broken down into atomic environments, each of which consists of a central atom and all other atoms within a certain cutoff radius. Each of these subsystems is processed by a separate feed-forward NN, unique to each element but shared across atoms of the same element. Each of these elemental NNs possesses the same set of weight parameters, and predicts an atomic energy contribution as a function of the chemical environment. These atomic energies are then added to provide the total potential energy of the system. Analytic gradients to determine the atomic forces are readily available.

In the input layers of the atomic NNs the information about the atomic environments is provided using vectors of symmetry functions Behler (2011). These serve as structural fingerprints with the required translational, rotational and permutational invariance. In a training process the weight parameters and biases connecting the different nodes of the network are adjusted until the output matches the available ab-initio reference data. For this, it is necessary to sample the PES of the system, by performing single point calculations at different positions of this PES. Once the NNP has been trained, it is able to predict the energies and forces at a fraction of the computational costs of the underlying electronic structure method, and can be applied to much larger systems than are accessible by DFT. For all details about the method, the training process and the validation strategies for high-dimensional NNPs the interested reader is referred to a series of recent reviews of the method Behler (2014, 2015, 2017)

II.2 Genetic Algorithm

Genetic algorithms Michalewicz and Janikow (1991); Deaven and Ho (1995) are a class of global optimization algorithms. They rely on codifying the properties of the system in form of “genes” of a given configuration, assigning a fitness criterion, which in the case of global structural optimization often is the potential energy, but other target criteria are possible Jensen et al. (2014), and then selecting and crossing over these genes with the purpose of generating new, fitter candidates. Along the way, variability is introduced into the candidate pool in the form of mutations that directly affect this genome. The final goal of this GO algorithm is to quickly generate and evaluate candidates, exploring the configuration space until the fittest candidate corresponding to the global minimum is found.

GA can be utilized for optimizing different kinds of systems Bozkurt et al. (2018), and the procedure is rather standardized for structural optimization of atomic systems Vilhelmsen and Hammer (2012, 2014a); Huang et al. (2019); Buendía et al. (2017); Heydariyan et al. (2018). The genes in this case consist of the coordinates of the atoms of interest, and elemental composition for multi-elemental systems. Crossovers between candidates are achieved by cutting the system along random planes Deaven and Ho (1995) and recombining the generated parts while maintaining the system’s size and stoichiometry. After crossover events, the new configurations are minimized, in an approach similar to basin hopping Monte Carlo Wales and Doye (1997) (BHMC), and it is acknowledged that both algorithms in fact sample the same transformed PES configuration space Wales and Doye (1997).

Strictly speaking, a genetic algorithm requires that the properties of the system are encoded into strings that can then be mutated and combined, such as for example listing the elements at given positions of a defined crystal structure Jensen et al. (2014). Modifying directly the coordinates of the system as we do in this work is more akin to modifying the phenotype rather than the genotype of the configurations. This kind of algorithms are more accurately described as “evolutionary” algorithms. We nonetheless refer to it as a genetic algorithm throughout the text as this is the more commonly in-use designation, and evolutionary algorithms are a subset of genetic algorithms.

The fitness is usually evaluated through the potential energy since the goal is to find the lowest energy configuration. Mutations are applied as direct, large-scale modifications to the coordinates of the candidates, such as mirrorings through a plane, i.e., discarding one half of the cluster and mirroring the other half, which leads to new configurations if the original cluster was not symmetrical through that plane; or rotations, such as partial rotations of a section of a cluster, or rotations of a supported cluster over a support. This results in a rapid exploration of the PES by generating and evaluating thousands of structurally different candidates in quick succession.

Some GAs make use of a generational ansatz, where the algorithm proceeds in steps by fixing the fittest population, generating new candidates, evaluating them, and proceeding to recalculate and update the fittest population. This ansatz is ideal for ab-initio based searches, since there is a long wait time for evaluating the new population that can be partially reduced by parallelization within a generation. This has been shown not to be necessary for a successful GA search Vilhelmsen and Hammer (2014b), and our NNP based calculations are extremely fast compared to similarly sized ab-initio relaxations. For this reason we utilize a continuous approach, where the fittest population is updated every time a relaxation is finished, and new candidates are continually pulled from this population.

Mutations are the key for efficiently exploring the PES. Similar to Monte Carlo simulations, poor mutations will advance the system through an inefficient route, and only rarely new good structures will be generated. Mutations should thus result in big changes to the candidate, but at the same time generate “reasonable” structures, which is aided by the final geometry relaxation every candidate goes through. In our case we have made use of the following mutations, which are applied only to the atoms belonging to the cluster, while the atoms of the support only change positions in the subsequent local geometry optimizations:

Rattle mutation: Inspired by basin hopping Monte Carlo Wales and Doye (1997), all atoms in the cluster are displaced at random up to a maximum amount. The displacement is slightly biased in the positive $z$ direction to avoid atoms penetrating into the support. This move helps the cluster to expand across the interface with the support and also to grow in height.

Twist mutation: The whole cluster is twisted by a random angle around a vector that passes through its center of geometry (COG) and is perpendicular to the slab surface. We note that the COM is used here as an element-independent generalization instead of the also frequently used center of mass (COM) of the cluster, since the COM does not coincide with the COG for clusters with more than one element. The cluster is also lifted a small random amount in the $z$ direction to avoid repulsive overlaps with the topmost surface atoms. This move helps the cluster to align itself with the atomic patterns of the support.

Angular mutation: Several atoms of the cluster are picked based on their distance to the COG with the furthest away atoms being picked first, and displaced to positions on a hemisphere centered on the COG and whose radius is the average distance of all atoms in the cluster to the COG. This move has been inspired by Rondina et al. Rondina and Da Silva (2013), which describe other possible interesting moves that can be adapted from BHMC. This move helps the cluster grow out in a 3D mode, avoiding flat clusters.

Mirror and shift mutations: A random plane passing through the COG of the cluster is selected, and the cluster is mirrored through this plane. To avoid frequent overlaps at the mirror plane, the mirror image is shifted a small amount away from the plane at random, which increases the acceptance ratio of the move. This move helps to generate more symmetric configurations, which tend to have a lower energy.

Molecular dynamics mutation: A short MD simulation in the $NVT$ ensemble (constant number of particles, volume and temperature) is performed including cluster and support atoms at high enough temperatures ensuring that the cluster atoms are mobile, which can even be above the melting temperature of the cluster, and then the system is cooled down slowly to avoid completely amorphous configurations. This move is inspired by simulated annealing Kirkpatrick, Gelatt, and Vecchi (1983) and minima hopping Goedecker (2004) optimization procedures, and it is well-known that MD is good at finding transition paths over low energy barriers Goedecker (2004). The cluster is thus allowed to naturally jump over potential energy barriers as provided by the available kinetic energy. This takes advantage of the shape of the PES, instead of being completely random such as the rattle mutation. This move is more expensive to run, since it requires a couple of thousand of MD steps, and it would be extremely expensive to implement in an ab-initio based search, but it is straightforward to perform with the NNP approach. Even then, this move requires more computing time than the other simple mutation plus local relaxation moves. For this reason, this move is given a lower occurrence probability.

Interatomic distances, within the cluster as well as between cluster and support, are checked after every type of move to avoid too repulsive atomic overlaps, and if this check fails the mutation is repeated a number of times until successful, or discarded if too many failures accumulate. This is most relevant for ab-initio based GA searches, where long relaxations due to closely overlapping atoms would cost a lot of computation time. But, this is still important in our case, since unphysical atomic positions can give rise to configurations that are outside the training set of the NNP, leading to less reliable energy and force predictions due to extrapolation.

We find it necessary to remark that although the moves might bias the final cluster shape in certain directions, they only provide an avenue of change, since in the end the energy of the cluster as evaluated through a particular force field dictates whether the structures are accepted or not. This means that it is hard to over-bias the cluster, but on the other hand it is possible to not offer enough mutations, resulting in the cluster exploring only a part of the configuration space.

A key component to a successful GA run is maintaining structural diversity in the breeding population. Otherwise, this population can become overrun by a small number of very fit candidates, which leads to always breeding the same children configurations and the GA search stalling. This structural diversity can be achieved in a number of ways, such as penalizing candidates that have already been chosen for breeding Vilhelmsen and Hammer (2014b). In general, a way of quantifying structural similarity is needed. In this case we use an interatomic distance comparator Vilhelmsen and Hammer (2012), which is insensitive to rotations and translations of the compared atoms. A first filter is applied comparing the energies of the structures in question. Afterwards, the comparator value is calculated employing the interatomic distances. This comparator is quick to compute and easy to calculate, which makes it ideal for this application. In our experience it is effective at keeping the breeding population distinct.

The size of the actively breeding population is kept fixed at a predefined maximum number of candidates, which correspond to the fittest candidates that are also structurally different. The fitness criterion is calculated as Vilhelmsen, Walton, and Sholl (2012):

where $f_{i}$ is the fitness of a particular candidate, $E_{i}$, $E_{\textrm{min}}$ and $E_{\textrm{max}}$ correspond to the energies of the candidate and the minimum as well as maximum energies in the whole population. This creates a fitness distribution that falls sharply as we get away form the optimal candidate and has a constant low value for the worst candidates.

Candidates are then extracted at random from this population utilizing a roulette wheel algorithm Johnston (2003), where the following criterion is evaluated,

to decide if a candidate is kept, where $f_{i}$ is the fitness of the randomly chosen candidate, $f_{\textrm{max}}$ is the fitness of the fittest candidate, and ran generates a random number between 0.0 and 1.0. With this equation, fitter candidates will tend to get chosen more often since their ratio of $f_{i}/f_{\textrm{max}}$ will be closer to 1.0, but all candidates in the active population have a non-zero chance of being picked.

Other GO methods have been considered for this system. Grid-based Rahm and Erhart (2017) optimization methods are not ideal for our purpose, since the resulting atomic grid of the clusters is not known beforehand, and as seen in the Results sections, is can be strongly modified by the support. Instead, methods with flexible grids could be tested Paleico and Behler (2020). BHMC Wales and Doye (1997) is in our experience too slow, since structures are still connected through a Monte Carlo-like chain. Instead, GA offers the advantage of evaluating many structures in parallel. Of course, a number of BHMC simulations could be run in a trivial parallelization approach, but GA benefits from biasing future candidates with information from previous candidates, thanks to the crossover operations.

III Computational Details

The reference method we have used to generate the training set for parameterizing the NNP is density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965), as implemented in the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996); Kresse and Joubert (1999) version 5.4.4. This software performs periodic DFT calculations with plane wave basis functions treating the core electrons with PAW pseudopotentials Blöchl (1994); Kresse and Hafner (1994). We have utilized the pseudopotentials recommended by the VASP manual for each element. We have chosen the generalized gradient approximation (GGA) functional suggested by Perdew, Burke, and Ernzerhof (PBE) Perdew, Burke, and Ernzerhof (1996) to describe electronic exchange and correlation.

We have run convergence tests for all relevant settings of the DFT calculations, until energy and force values converged below the aimed NNP root mean squared errors (RMSE), i.e., 1 meV/atom for total energies and 100 meV/Bohr for the forces. For this purpose, we utilized a plane wave cutoff of 500 eV, and a k-point grid of $12\times 12\times 12$ for the 3.61 Å cubic unit cell of fcc Cu containing 4 atoms. For larger supercells the number of k-point has been scaled accordingly to ensure a similar level of convergence. For slab systems, at least 13 Å vacuum have been used in the $z$ direction perpendicular to the surface to minimize interactions with the periodic images of the slab. We applied a Gaussian smearing at the Fermi level, with a 0.1 eV $\sigma$ factor. An example input file for VASP with a complete list of the settings is given in the supporting information.

Although the NNP can be trained to DFT reference data including dispersion corrections Kondati Natarajan and Behler (2017b); Quaranta, Behler, and Hellström (2019); Morawietz et al. (2016), which are often problematic in GGA functionals, we have not done this in the present work. The reason is that in particular the dispersion-corrected PBE functional often yields overbinding, which is present to some extent also for the plain PBE functional, and too small lattice constants are obtained with commonly used correction schemes Grimme et al. (2010); Tkatchenko et al. (2012); Bučko et al. (2014), which is in agreement with previous reports of overbinding Antony and Grimme (2006) for certain exchange-correlation functionals. As we consider a good prediction of the bulk properties of copper and zinc oxide as important, we have chosen the uncorrected PBE functional providing the most accurate description. It should be noted that, as discussed below, the energy differences between local minima can be very small and the subtle energetic ordering can change depending on if and which dispersion correction scheme is used.

III.2 Generation and Composition of the Reference Data Set

The key to generating a reliable NNP is to produce a reference data set covering all atomic environments that are relevant for the chosen type of simulation. Since carrying out these electronic structure calculations is the most time-consuming step of the NNP construction, this data set should be kept as small as possible. This avoids wasting computational time on calculating ab-initio data on configurations that a given type of simulation will not visit, or the reverse, missing training data for configurations that emerge in a simulation. For example, Monte Carlo simulations are more likely to visit high energy configurations corresponding to less-favorable relative atomic positions when compared to molecular dynamics simulations, that generate trajectories guided by the interatomic forces.

It is thus important to tailor the structure generation method to the simulations being performed, and if possible to integrate both processes Artrith and Behler (2012). Several aspects need to be considered. Structures need to be generated covering a wide range of the PES, but not at random since this samples irrelevant high-energy regions. Structures also need to fit within the available computational resources, i.e., they should contain only up to about 200-300 atoms per structure in case of DFT, while also avoiding to compute identical or very similar environments repeatedly.

The protocol we have followed for the generation of the NNP in the present work can be summarized and categorized into the following phases, also presented in fig.1:

In phase I, structures are generated from known configurations, i.e., known crystal structures and surface models, for which either experimental data or theoretical information is available. In our case, this corresponds to Cu and ZnO bulk structures in their fcc and wurtzite crystal structures, respectively; low Miller index surfaces of both materials, and Wulff-shape nanoparticles. Additionally, the NNP includes some additional data for brass structures Weinreich et al. (2020). These base structures are then modified in a well-defined way. For example, the experimental lattice parameters and interatomic distances of the bulk structures are varied in a small range, to accommodate for possible DFT deviations from the experimental values, and to provide the NNP with more information around the equilibrium conditions. Additionally, thermal distortions are emulated by randomly displacing the atoms around their equilibrium positions, which is also important to break the symmetry of the atomic configurations. For slab structures, different numbers of layers are included. This results in a NNP that is applicable to a narrow range of conditions centered around the base structures, and is useful for finding simple structural and energetical parameters Artrith, Hiller, and Behler (2013), such as lattice constants, bulk moduli, and surface energies.

In phase II, $NVT$ and $NPT$ ensemble simulations are utilized. With these simulations the system is naturally allowed to explore the available configurations at given temperatures and pressures. Included in this is melting and subsequently cooling the systems to introduce realistic liquid and amorphous environments into the data set, which would be difficult to achieve with the procedural methods of phase I. Structures are sampled at regular intervals, e.g. every couple of picoseconds, or by utilizing an ensemble of NNPs Artrith and Behler (2012); Behler (2015) to scan specifically for new configurations. After phase II, the NNP is usually able to perform long $NVT$ and $NPT$ MD simulations without visiting unsampled regions of the PES. In the present case we require both high energy configurations, generated e.g. by GA mutations and crossings, which are not well sampled in a canonical ensemble, and close to equilibrium low energy structures. Other simulation methods that allow the NNP to explore reasonable sections of configuration space are also useful, such as Monte Carlo atom swaps Weinreich et al. (2020) that exchange elements, simulated annealing Kirkpatrick, Gelatt, and Vecchi (1983), and semi-grand canonical ensemble simulations Kofke and Glandt (1988).

Phase III aims to include high-energy structures and consists of two parts. First, random copper clusters are generated on different ZnO surfaces. Structures are checked for short interatomic distances and rejected if atoms are too close together to avoid unphysically large, spurious repulsive forces, which can lead to problems both in the DFT calculations and in the fitting procedure. We set this threshold at 20% below the nearest neighbor distance for the respective material, which is sufficient to cover the usual deviations from this distance observed in clusters and provides the required information about the relevant repulsive parts of the PES. This enables the NNP to perform the second part, where iterative GA searches are carried out. In this iterative approach, GA searches are run based on a given NNP for different cluster sizes, which is very fast. Both minimized and unrelaxed structures resulting from a mutation or crossover operation are then extracted and evaluated, again utilizing an ensemble of NNPs Artrith and Behler (2012); Behler (2015). With this procedure we obtained about 200-500 new structures per GA run. These structures then have been recalculated by DFT, a set of new NNPs has been generated for the extended data set, and this process has been repeated iteratively several times until convergence has been reached, which is defined here by two NNPs trained on the same data set agreeing on the predictions for the obtained clusters within the accepted energy and force tolerance.

Two problems need to be solved in this hierarchical approach. The first is to avoid recalculating structures that are already well represented in the data set. The use of ensembles of NNPs can assist the selection, since structures that are already included will exhibit small errors in all NNPs trained on this data set, while unknown structures will have very different errors for different NNPs. A complementary method based on the quick classification of SF sets and atomic environment descriptions has been employed here as well, and also other algorithms like CUR are in principle available to select the most distinct training structures Imbalzano et al. (2018). Of course, repeated environments cannot be avoided completely, since for example the presence of a slab support will often lead to redundant atomic environments, in particular in the bottom-most layers of the support, which hardly change while the cluster and the top layers are optimized.

The second problem to solve is extracting smaller structures from the very large systems used in the simulations, once atomic environments that need to be added to the DFT data set have been identified. This is a problem particular to the subset of structures derived from supported clusters, since other structures that we have used in training the NNP (bulk solids, slabs, etc.) are already intrinsically periodic and designed to fit within the available ab-initio calculations.

To avoid periodic images affecting the GA search, which can lead to self-interacting periodic configurations such as wires and sheets of copper atoms on the surface, the GA simulations are performed on large slab supports with up to a thousand atoms. These structures are too large for routine DFT calculations, and even if they would be feasible, they would offer very little new information with respect to the required computational cost since most support atoms are too far away from the cluster to notably interact. To reduce the system size for the DFT calculations, we need to extract the relevant atoms from the simulations Eckhoff and Behler (2019); Weinreich et al. (2020) while preserving the local atomic environment of these atoms. In the particular case of supported clusters we also need to maintain the periodicity of the extracted subsystems, which is required by the VASP code, which is not able to perform non-periodic calculations, and, more importantly, for obtaining realistic structural models without artificially merging ZnO surfaces at the edges of the periodic boundary condition.

To achieve this, the key realization is that rectangular cuts that preserve the periodicity of the underlying surface lattice will also tile periodically correctly. This assumption holds, of course, if the atoms in the slab have not been substantially distorted from their equilibrium positions. To account for possible small distortions due to cluster-support interactions, our cutting algorithm checks that i) the number of copper atoms in the cell and in particular the zinc-to-oxygen ratio has been preserved and ii) in the cut-out subsystem, atomic distances across the periodic boundary condition are reasonable. This sometimes requires centering the cutting surfaces on different atoms of the system until a cut fulfilling requirements i) and ii) has been found. Additionally, slightly adapting the size of the cutting box from the the size dictated by the ZnO lattice constants by 1-3% might be required. Taking into account these considerations, it is possible to obtain periodic subsystems with minimal mismatch at the cell boundaries that are suitable for routine DFT calculations. This procedure is shown in fig. 2.

III.3 Construction of the Neural Network Potential

The final data set obtained from the procedures described in the previous section consists of 73,136 structures covering about 5 million atomic environments. Therefore, in total 73,136 energies and about 15 million force components are available for training the NNP. The data set contains a wide range of configurations, including bulk, slab, cluster, and amorphous structures for the subsystems Cu, ZnO, and CuZn as well as supported cluster structures for the full ternary (Cu,Zn,O) system.

We have utilized our in-house program RuNNer Behler (2015, 2017) for the construction of the NNP, which is freely available under the GPL3 license. To choose the best architecture, we have optimized different NN architectures using the hyperbolic tangent activation function in the hidden layers and a linear function for the output layer. For simplicity we have chosen to use the same architecture for each element in the system. The parameters of the symmetry functions used to describe the atomic environments, which have been derived from previous work on ZnO and Cu Artrith, Hiller, and Behler (2013) with some functions removed without significant loss of accuracy, are given in the supporting information. The resulting errors are shown in Table 1. Based on these errors we have selected the architecture with 2 hidden layers and 15 nodes in each layer for the subsequent simulations. This architecture yields a 2.5 meV/atom root mean square error for the training energies and 59.7 meV/Bohr for the training forces, with very similar values for the testing data not included in the fitting process indicating the absence of significant overfitting. As further evidence of the quality of the potential, we present a predicted vs. DFT energy plot in our supporting information for the structures in our dataset.

III.4 Genetic Algorithm

Utilizing a NNP for the GA search allows us to perform the simulations in large cells, with 448 support atoms, with even larger cells also being possible. Each cell was generated to have a lateral diameter of at least 20 Å in each direction of the ZnO(10$\bar{1}$0) surface to provide enough separation of the clusters from their periodic images, resulting in a final size of 23.03 $\times$ 21.19 Å2, or a 6 $\times$ 4 supercell. Each slab was generated with enough layers to provide a thickness of at least 7 Å, which avoids atoms on top of the layer from “seeing” past the atoms on the bottom-most layer due to the NNP symmetry function cutoff. The vacuum added in the $z$ direction (at least 13 Å) ensures that the top-most atoms of the cluster do not interact with the bottom-most atoms of the slab through the PBC.

For each cluster size, the algorithm is seeded by generating and optimizing 15 clusters starting from random geometries. These clusters are initialized in an area in the center of the supercell to avoid the randomly placed atoms being spread too far apart, which leads to poor initial structures. Random initial structures too far from reasonable configurations might lead to extrapolations of the NNP since the underlying data set samples mostly the energetically favored compact clusters, while very high energy configurations, e.g. isolated atoms may not be covered in all situations with a very high precision, which is acceptable as they are physically less relevant. Indeed, after this initial random sampling of the PES the GA search quickly abandons these initial random structures due to their high energies.

The supporting clean slab is pre-optimized at the beginning of the GA search, and then every new candidate cluster is relaxed together with the supporting interface. We note that this is sometimes avoided in DFT-based GA simulations because of the added costs, but as we will show in Sec. IV.5, this can lead to different results. For the minimization we use the LBGFS algorithm (limited memory BFGS Liu and Nocedal (1989)), with a maximum force criterion of 0.005 eV/Å and a limit of 1000 evaluation steps, as implemented in the ASE library Larsen et al. (2017). This optimizer was chosen mainly for being the fastest for our particular setup and rather robust. Care should be taken that the structures are properly optimized, as we found that otherwise false LM can appear in the final candidate pool that relax into other LM when stricter criteria are used, particularly for small clusters.

The emerging active breeding population is kept at a maximum of 50 structures during the simulations. It is ensured to be structurally distinct by utilizing the interatomic distance comparator Vilhelmsen and Hammer (2012). This comparator has three criteria that need to be set. We use a minimum total energy difference of 0.1 eV. If the energy difference between two structures is larger, they are assumed to be different. A maximum correlation distance of 0.7 Å, and a maximum cumulative difference of 0.025 Å are also utilized for the comparison of interatomic distances. For the last two values, the larger the number, the less “strict” the comparator is. If a too-strict comparator setting would be used, small differences in structures within the minimization convergence threshold might be detected as significant, resulting in many clusters corresponding to the same structure. In our experience, these values are suitable to keep the breeding population structurally distinct (that is, varied).

The mutation moves we employed have a number of parameters that can be adjusted. We have utilized the following set of values:

60% of the atoms in the cluster are rattled, chosen at random.

Maximum displacements of up to 3 Å in the $x$ and $y$ directions are used. This value is chosen because the distance between sites on the ZnO surface is also in the order of 3 Å, and we want the cluster atoms to be able to jump to a neighboring site. Displacements in the negative $z$ direction towards the surface are limited to 1 Å to avoid cluster atoms burrowing into the support, while positive displacements remain the same as in the other directions. The probability is uniform in this range.

A vertical offset of 0.5 Å is added to the cluster after twisting it, to avoid cluster atoms overlapping directly with the support. This offset is removed by the geometry relaxation that takes place after the mutation.

between 25 and 50% of the atoms in the cluster are displaced by the move, starting with the atoms farthest away from the COG.

Atoms are positioned at random on a hemisphere centered on the COG of the cluster, with a radius that is the average distance of all the atoms in the cluster to the COG. This radius is then modified by a random amount between -0.5 and 3.0 Å.

The mutation can be biased to position atoms towards the top of the hemisphere and thus promote cluster growth in the $z$ direction by restricting the allowed spherical angles, but in our case we have not added this bias.

The angle between the mirroring plane and the surface plane can be biased to favor mirrorings parallel or perpendicular to the surface, but we have left this unbiased in our simulations.

Mirrored atoms are shifted to avoid overlaps close to the mirroring plane. This shift is a random value between 0.0 (that is, no shifting) and 1.5 Å. Since there is usually already a small gap at the mirror plane, this amount is sufficient to avoid serious overlaps given the usual Cu-Cu nearest neighbor distance ($d_{o}$) of approximately 2.55 Å.

The MD is run for 3000 steps with a 5 fs timestep, for a total of 15 ps.

The temperature is set to 1000 K, which enables to melt small Cu clusters (the bulk copper melting temperature is around 1360 K, but due to the high surface to volume ratio clusters have lower melting temperatures Cui et al. (2017); Weinreich et al. (2020)).

Each search for each cluster size is repeated 4 times, until 3000 structures have been evaluated in each run, for a total of 12.000 structures. Clusters at each cluster size are then merged into a single dataset before continuing with the analysis performed in Sec. IV.2. The same structures are usually found in each run (see fig. 3 a), described below), but the independent runs are kept as a consistency check, and we can expect this situation to be different for larger clusters where finding a given putative GM becomes more difficult Paleico and Behler (2020).

The GA search described here has been implemented around the already available functions in the ASE (Atomic Simulation Environment) Python library Larsen et al. (2017), version 3.17.0, with some extra functions as required. In particular all of the mutations needed to be either modified (rattle and mirror) or implemented from scratch (the rest) following the library’s templates. The GA library is based around a database module, which makes it easy to store and retrieve structures in an organized way and query the results of the GA search. Due to the large amount of structures that can be calculated with the NNP, the database needs to be purged regularly from about 500 stored structures onward keeping the fittest candidates, otherwise the updating process for the current population becomes excessively slow. Energy and force evaluations were provided by the n2p2 neural network library Singraber, Behler, and Dellago (2019) in combination with LAMMPS Plimpton (1995).

To illustrate the course of a typical simulation Figure 3 shows the results of running four GA searches for the GM of Cu10. In a) we can see how the GA search quickly finds the GM in fewer than 300 structures for all runs. The energy before this point decreases rapidly, and no lower energy structures are found in the following 2700 configurations generated in each run (for a total of 3000 structures per independent run).

Figure 3 b) shows how the different mutations (or no mutation, just pairing) affect the energy of the generated relaxed candidates, for two GA searches for Cu10. The energy is expressed relative to the parents’ average energy, as given by:

Just pairing with no mutation seems to be the most consistent operation, with energy changes clustering around the 0% mark. Angular mutation moves generate the highest energy configurations, which is usually a consequence of fewer atoms being attached to the surface. Molecular dynamic runs are very effective at generating lower energy configurations, with their data set the only one with a median energy change below zero. They still can generate very high energy configurations due to the occasional formation of separated isolated atoms. Rattle, mirror and twist mutations exhibit a similar behavior between them, with the twist mutation obtaining the largest energy gain.

IV Results

The structure of the support is very important for the obtained optimized clusters. For this reason, we have carefully checked the description of the bulk material as well as of the relaxed geometry of the ZnO surface provided by the NNP. The optimized lattice parameters of the bulk ZnO wurtzite structure we obtained are presented in Table 2. They are basically indistinguishable from the reference DFT results, and are only 2% above the experimental valuesLide (2009), which can be ascribed to the employed exchange correlation functional.

For the ZnO(10$\bar{1}$0) surface that we use in this work, Figure 4 shows all the relevant angles and distances as obtained from a geometry relaxation of a 4 layer, $(4\times 7)$ ZnO(10$\bar{1}$0) slab. We observe the known geometry changes due to relaxation from the introduction of vacuum interfaces, with layer distances shortening in particular close to the surface, and tilting of the ZnO dimers to an angle of 9.85 degrees. These values are very similar to the ones obtained from DFT calculations, and to previous results reported in the literature Meyer and Marx (2003).

IV.2 Global Optimization Results

In the next step, we have employed GA simulations to identify the GM and energetically favorable LM structures of a series of cluster sizes at the ZnO(10$\bar{1}$0) surface. Figure 5 shows the results of running the GA algorithm for Cu clusters between 4 and 10 atoms. Four separate runs have been carried out per system yielding 3000 structures each, and then the results of the runs for each cluster size have been combined and analyzed. Before doing this, as shown in Fig. 3 it was confirmed that different repeated runs tend to find the same GM although; since the process is stochastic, this is not necessarily always the case Vilhelmsen and Hammer (2014b). Still, most of the minima have been found multiple times, within a single run and also across runs.

We find a rich variety of structural patterns. First, we note that it is possible to detect a number of “geometric families” that are consistently present in the best ranked clusters. One important general trend is that Cu atoms tend to interact strongly with the oxygen atoms of the surface, preferring a direct adsorption, or the adsorption between the Zn and O of a ZnO dimer. Further, they can also be located between two surface oxygen atoms, like the Cu atom at the bottom of the image of LM 2 for Cu5 in Fig. 5.

For Cu4, we identify the simplest structural patterns, and here we can derive our second important principle: Even for small sizes the clusters already adopt a 3D shape, and do not lie flat on the surface. The same is observed for all other cluster sizes, with only a few rather flat configurations being found among the lowest energy clusters although they do appear at higher energies. Additionally, if somewhat trivial, the Cu atoms do indeed form clusters and do not spread across the support into individual adsorption sites. This clustering is well known from experimental results even at low coverages Dulub, Boatner, and Diebold (2002), and can be predicted from the relative binding and surface energies of the materials involved Henry (2005). Some separated clusters are observed in simulations of the larger clusters as a consequence of certain mutation moves, but these always correspond to high energy configurations. Finally, we notice that neither here nor for any other size do the clusters form extended periodic structures in one or two dimensions, as has been proposed in other studies Mora-Fonz et al. (2017). These kinds of structures were only ever observed in our simulations for very small periodic cells, which prevent the formation of non-interacting clusters due to rather high copper coverages.

Due to the small number of atoms in the Cu4 cluster, it appears to form only two main structural patterns at low energies. The first pattern is a truncated square pyramid with only 3 atoms present on the support, which is found in the GM in Fig. 5 and in LM 5. Notice that the GM and LM 5 are in fact very different with respect to the support, as the adsorption footprint on the ZnO(10$\bar{1}$0) surface is not symmetric because of the orientation of the ZnO dimers corresponding to the polar (0001) and (000$\bar{1}$) surfaces. In this direction, we have a short Zn-O distance within the dimers, and a long distance between different dimers. For both structures, the GM and LM 5, the three copper atoms attached to the surface are in the vicinity of oxygen atoms. However, in the GM two of them are in between two surface ZnO dimers and only one is found on top of a dimer, while in LM 5 the opposite is observed. This indicates a preferential direction, since rotating the cluster by 180∘ does not result in the same support environment.

The second pattern for the Cu4 clusters found in LM 1, LM 2 and LM 3 corresponds to small, more or less distorted tetrahedra. Perfectly stacked structures like LM 2 and 3 cannot easily be realized for bigger clusters, but are sometimes observed for specific cluster sizes at higher energies. Notice that also these structures exhibit the preferential directionality mentioned in the previous paragraph, and are elongated to some extent in the ZnO dimer direction. Due to their structural similarity, we cannot fully exclude that these structures are connected by only small energy barriers or even correspond to the same shallow minimum, which might also be dependent on the employed exchange correlation functional or subtleties of the NNP, and further investigations using e.g. DFT-based nudged elastic band calculations Henkelman, Uberuaga, and Jonsson (2000) might be used to finally clarify this question. However, such calculations would be very expensive to carry out with DFT in our simulation setup.

Finally, we notice that all the Cu4 clusters are seemingly contained within a $(1\times 1)$ surface cell. The clusters do not form flat 1D or 2D structures, but rather gain height in the $z$ direction while remaining confined within this box.

For Cu5, many of the structures are natural extensions of the Cu4 structures. In particular, the GM is the completed square pyramid with the missing atom added in one corner, while LM 1 is a distortion of this structure. LM 2 looks like a extended tetrahedron emerging from LM 1 of Cu4, and shows a preferred orientation. In LM 3 the distortion of LM 1 extends to both sides of the pyramid base making it very likely that this distortion represents true local minima in LM 1 and LM 3. The local minima LM 4 and LM 5 resemble distorted structures emerging from the GM of Cu4 by the attachment of an additional copper atom.

For the Cu6 clusters we notice a number of interesting developments. The GM now corresponds to the structure of LM 3 of Cu5 with an atom attached. We will call this pattern a saddle. For LM 3 we observe a first cluster with distorted five-fold symmetry, and a number of rather flat clusters occur that follow the corrugation of the surface in the cases of LM 2 and LM 5. More interestingly, LM 2 is the first cluster to extend beyond the $(1\times 1)$ surface cell, which will become the norm for larger clusters. This new extension is larger in the “short” direction of the surface, perpendicular to the direction of the Zn-O dimers. This is reasonable, since this distance is shorter and thus easier to “bridge over” as the cluster grows.

For Cu7, many of the clusters appear to be directly related to those at Cu6. The new GM is a direct evolution of the previous saddle-like GM of Cu6, with another related structure at LM 5. LM 1 and 4 are related to the five-fold symmetry cluster at the LM 3 of Cu6, and LM 3 is similar to the previously flat triangular cluster. Of these clusters, only 2 remain within the $(1\times 1)$ cell, which are interestingly the lowest and the highest energy structure, with the rest increasing in footprint.

For Cu8, the GM is a moderately distorted version of the cluster at the LM 1 of Cu7, with an extra atom added in the lower left corner of the image, but it can also be thought of being derived from LM 4 thus unifying two rather stable sub-patterns. For LM 1 and 5 we notice new modifications of the saddle, a new five-fold derived cluster emerges in LM 3, and a first cluster with a buckled hexagonal base pattern plus an attached atom at LM 4.

For Cu9, the saddle-based pattern once again becomes the GM. Additionally, at this size all the shown clusters are for the first time outside the $(1\times 1)$ cell.

Finally, for Cu10, by attaching an additional atom the previous LM 1 of Cu9 becomes the GM, with the saddle configuration falling to LM 2. However, it must be noted that in particular for Cu10 the energy differences are extremely small, and cannot be resolved with full confidence, neither by the NNP nor by the underlying DFT calculations.

These analyses of the structural and energetic trends leads us to a general conclusion: the shapes of these clusters are a compromise between the interaction between copper and the support, which would tend to favor spread out, flat clusters, and the cohesive energy of the copper clusters, which tends to favor three-dimensional clusters. Additionally, the granular nature, corrugation and symmetry of the support, with preferred adsorption sites, and the strain that it induces on the clusters, complicates this interaction. This is a known effect that can be observed even with simple potentials Eckhoff, Schebarchov, and Wales (2017). In Fig. 6 we plot the binding energy of the GM structures as a function of the cluster size. We define the binding energy as:

IV.3 Interface Structure

An interesting question is what local structure the clusters adopt at the contact with the interface. We can expect that for very large clusters, far away from the influence of the support, the central atoms will organize in the geometry of fcc bulk copper. However, right at the interface, a number of opposing trends influence the structure. First, it is known that finite solids such as nanoparticles tend to adopt shapes dominated by the need to minimize the surface energy, as described by Wulff Wulff (1901). This, in principle, applies to large structures where surface and bulk energetic terms dominate, that are in equilibrium, with only one element, and more importantly, considering solid-vacuum interfaces, although attempts have been made to extend Wulff’s theory to interfacing materials Henry (2005) and alloys Ringe, Van Duyne, and Marks (2011). Since here we are addressing very small clusters, which in addition are strongly influenced by the ZnO interface, we have to take into consideration its effect. The degree of interaction between the cluster atoms and the support can lead to new preferred geometries. This can be due to low-energy interaction sites dictating a given geometry, and specific surface structures favoring these interactions. This can be interpreted in various forms, such as strong metal-support interactionTauster (1987); Pan et al. (2017); Figueiredo et al. (2019), or the theory of lattice mismatch Eckhoff, Schebarchov, and Wales (2017); Koda et al. (2016).

The Cu-ZnO system presents such an inherent lattice mismatch between the two materials. If we take one low index surface from each material, it is not possible to produce a coincident supercell without straining, i.e., expanding, contracting, or even reconstructing, at least one of the two subsystems. Such lattice mismatches can induce defects and strain in structures, which are known to enhance catalytic properties Mavrikakis, Hammer, and Nørskov (1998); Amakawa et al. (2013); Bissett et al. (2013), and it is still present to some extent for finite structures such as clusters. It is partially reduced by losing the requirement for a periodic structure, which allows for lateral changes of the interatomic distances, but the atoms in the interface area are still guided by the underlying geometry of the support to adopt regular positions.

We note that the clusters studied in this work are very small and thus regular surface facets corresponding to the low-Miller index surfaces in a Wulff-like picture cannot be expected, but still small structural features geometrically resembling these surfaces may emerge, which we will now investigate. The interfacial atoms of the minimized clusters of Fig. 5 are shown in Fig. 7. For this purpose we define interfacial atoms as those within a 1.6 Å vertical distance of the lowermost Cu atom of the respective cluster. This interval was chosen as it is thick enough to potentially contain a monolayer of all the low index Cu surfaces. The ideal (100) and (111) monolayers are flat, but the (110) monolayer has a thickness of about 1.275 Å, i.e., one-half of the nearest neighbor distance, $d_{0}$. The distance criterion has been chosen about 0.3 Å larger than this value to accommodate also distorted positions of Cu atoms on top of the support, with vertical distances falling in the range of 1.3 to 1.6 Å. Additionally, those atoms further away from the substrate are displayed in a perspective view with larger radii, providing a sense of the 3D structure of these interface atoms. Identifying surface features at the interface in such a heterogeneous system is not trivial, but we attempt an analysis by a comparison with the well-known low Miller index surfaces of copper, (111), (100) and (110), which are known to be the most stable ones Tyson and Miller (1977) for the material in vacuum, in the given order Tran et al. (2016).

For Cu4 and Cu5 we observe a variety of patterns, including some triangular configurations. The GM of Cu5 is reminiscent of the (100) surface, organized in approximately square units of Cu atoms. This surface, which is the second most stable out of the three low index cuts in vacuum, only appears for the smallest clusters. In LM 3 of Cu5 we observe the first example of a structural pattern that repeats for many of the larger clusters: five interfacial atoms, structured as a rectangle with an atom in the center that is slightly above the plane formed by the other atoms. Both, the buckled (110) and the planar (111) surfaces of fcc crystals are somewhat related. The (111) surface is usually described as a planar layer of hexagonally arranged atoms, but can also be reinterpreted as a rectangle plus a central atom-in-plane pattern. Its hexagons are completely symmetrical with equivalent bonds. The rectangle has two characteristic distances: a short one equal to $d_{0}$, the nearest neighbor distance of 2.55 Å, and a long one equal to $\sqrt{3}a_{0}$, the lattice constant of bulk copper, yielding 4.42 Å. The (110) surface is usually interpreted as the rectangle plus central atom pattern, with the central atom belonging to a higher or deeper layer resulting in a buckled and rather open surface, but it can actually also be reinterpreted as asymmetrical tilted hexagons, i.e., they are longer in one dimension than the other one. The distances in the rectangle consist of: a short one, which is again equal to $d_{0}$, 2.55 Å, and a long one equal to $a_{0}$, 3.64 Å. This means that actually the (111) and (110) surfaces are in fact quite similar, with the key difference being whether the hexagonal pattern is flat or buckled and the related consequences for the interatomic distances. We can smoothly convert from one surface to the other by flattening the (110) configuration, which causes the other Cu atoms to move away and increases the length of the long side of the rectangle.

As can be seen in Fig. 7, at the interface we obtain a mixture of both behaviors, with clusters showing a structural intermediate of both surfaces. Some interfaces are flatter with larger distances in the long direction (closer to (111)), while in others we observe the opposite configuration (closer to (110)). In fact, this can vary within one cluster size or even one structure, such as the GM at Cu10, which exhibits two parts with the two different extremes. In effect we observe structures that are in a continuum between the two extremes of Cu (110) and (111), which is a consequence of the cluster shape at the interface adapting to the ZnO surface geometry, while still subtle distortions resulting from the Cu-Cu interactions are present. Due to the small size of the clusters they are rather flexible in shape and none of the atoms has the full copper environment of the bare metal surfaces in vacuum. This structural behavior is expected to change with an increasing number of copper atoms, finally resulting in copper interface structures more similar to one of the ideal copper surfaces.

IV.4 Properties

In the next step, we have analyzed a number of geometric properties of the obtained clusters to identify possible structural trends. Specifically, we have looked at:

The fractional coordination number Iannuzzi, Laio, and Parrinello (2003) of the atoms in the cluster. The fractional coordination number converts the usual integer-valued coordination number into a continuous value. For an atom within a certain distance of the central atom, like $r_{0}=2.55$ Å, the nearest neighbor distance in bulk copper, the coordination number contribution is 1. For atoms beyond, smaller values which rapidly decrease with increasing distance from the coordinated atom are used. In effect atoms beyond the nearest neighbor distance are considered as partially coordinating. Here, this avoids rapidly changing coordination numbers due to cluster atoms being slightly beyond the nearest neighbor distance, which frequently occurs for small metal clusters.

The height of the cluster, as given by the difference in $z$ coordinates of the top and bottom-most atoms in the cluster.

The maximum distance of any atom to the center of mass (COM) of the cluster. This shows the maximum extent of the cluster.

The average nearest neighbor distances, for each atom and for the whole cluster, with neighbors in the cluster up to 3.0 Å away.

The volume of the cluster as obtained from the QHull convex hull algorithm library Barber, Dobkin, and Huhdanpaa (1996). The convex hull algorithm returns the smallest convex polyhedron that contains a collection of points, i.e., the position of the cluster atoms. From this one can derive the properties of a polyhedron such as volume and surface area, but also other properties, for instance if an arbitrary point is inside or outside the constructed polyhedron.

The area of the atoms at the interface, also using QHull. Interface atoms are chosen following the definition in Sec. IV.3, i.e., within 1.6 Å vertical distance of the lowest cluster atom.

Some of these properties may have catalytic consequences, such as the coordination number of the cluster atoms, while others characterize the overall geometry of the clusters. We have plotted the results in Fig. 8, showing the values for the GM, and the average value for the GM plus the first 5 LM, with the bars indicating the minimum and maximum values of the property within this group.

In Fig. 8 a) we plot the fractional coordination number, and we can observe that this quantity is not yet saturated for clusters of size 10, since we are still in a rather small size regime that does not allow to approach the fcc bulk value of 12. The cluster height in b) does not show any clear trend for the different GMs, and there is quite a spread within each class, but the average shows a steadily increasing value in spite of the diversity of the cluster geometries. As expected, the average nearest neighbor distance in c) slowly approaches the equilibrium value (2.55 Å in fcc Cu) from below, i.e., most cluster atoms are closer than in a bulk environment. This is related to the pattern imposed by the ZnO support as explained in the previous section but also a general property of small metal clusters, as the remaining bonds in undercoordinated atoms tend to be stronger and shorter.

Fig. 8 d) shows that the clusters slowly increase in radius, and is the only property where the GMs show a common behavior in that all of them are at the average or below indicating a more compact shape. The plot of interface area with the support in e) shows a rather slow growth past 6-7 atoms, which implies that the extra atoms coming into the clusters are not used to substantially increase the footprint of the cluster on the support. This might be related to the fact that the cluster needs enough atoms to bridge the sometimes significant gaps between oxygen adsorption sites on the surface. Finally, f) shows the interface area, i.e., footprint on the support, to volume ratio of the cluster, which supports the analysis of e). Note that most of the GM lie at the lower end of this range, that is, they tend to optimize their surface to volume ratio, but this trend may not necessarily extrapolate to larger clusters.

Unfortunately, it does not seem possible to easily differentiate the GMs by utilizing just the calculated properties, although some general trends can be discerned, such as GM usually being more compact (Fig. 8 d), and with a minimized interface surface to volume ratio (see Fig. 8 f). We ascribe this to the overall very small cluster size in this study.

Finally, we have also colored the clusters by those properties that can be expressed on a per-atom basis, as seen in Figs. 9 for the fractional coordination number and 10 for the shift from the fcc copper nearest neighbor distance. As expected from the property plots, we observe a gradual increase in the coordination numbers with cluster size. We observe the first fully surrounded/coordinated Cu atom with 8 other nearest neighbor Cu atoms at the center of the GM of Cu9, (see Fig. 9 Cu9 a). This is the first size that can present a full coordination sphere, and this happens to be the putative GM at this size. The only other fully coordinated cluster appears as a LM for Cu10. Although it would be expected for fully coordinated structures to be energetically favored, this does not seem to be the case at these small cluster sizes.

In Fig. 10 we similarly observe a slow approach to the expected nearest neighbor distance of 2.55 Å. Notice that only two clusters present positive deviations from this nearest neighbor distance (i.e., a value above 2.55 Å), LM 2 of Cu4 and LM 3 of Cu9.

IV.5 Frozen Surface Optimization

An interesting question addresses the influence of the surface beyond offering merely a substrate for the cluster. In particular in ab-initio based searches Davis, Horswell, and Johnston (2016); A. Hussein, A. Davis, and L. Johnston (2016) the surface is sometimes kept fixed to save on computational time during the multiple minimizations required. In other simulations the surface is kept frozen in first, coarse minimizations to save on resources and then the best clusters from these pre-optimizations are further optimized with more expensive methods and flexible surfaces. This could lead to biased results if the initial filtering employs a frozen support. But what happens, if we do not allow the surface to relax to the presence of the cluster?

To analyze this, we first investigate the effect of the cluster on the surface. As an example, Fig. 11 shows the distortions in the ZnO slab for the GM of Cu10. In this case, some of the atoms below the cluster have been slightly displaced from their layer, being pulled towards the cluster. Fig. 11 a) shows in a top view the lateral distortion on the surface of the slab. We can see that very small lateral distortions reach all the way to the periodic boundary even for our very large surface cell, which shows another reason why utilizing large enough optimization cells is important. Smaller cells might inhibit or alter these distortions, and thus change the energetic ranking of the clusters. Further we can see that these distortions do not just follow the symmetry of the cluster but can exhibit pronounced directionalities close to the cluster due to the asymmetric orientation of the ZnO surface dimers. b) plots this effect in a cross section for the interior of the slab. In this case, we can see that the atoms are modified down to the bottom-most layer. This effect is also found for other clusters to varying degrees. This shows the importance of utilizing multiple support layers, as some ab-initio based GA searches reduce the size of the support down to even two layers only Davis, Horswell, and Johnston (2016); A. Hussein, A. Davis, and L. Johnston (2016); Kolsbjerg, Peterson, and Hammer (2018); Bazhenov and Honkala (2019); Engel, Francis, and Roldan (2019) to save on computational cost, or freeze the bottom-most layers for the same purpose. These extra layers are not very expensive for the NNP, since its computation cost scales well with system size.

Overall, this effect might seem to be still rather small, but it clearly demonstrates that there is some notable interaction between the support and the cluster that requires the support to slightly deform. As this might have consequences for the identification of the minimum energy structures, we have repeated the GA search for the Cu10 cluster to show how freezing the atoms of the surface affects the results of the search. The obtained GM and the first 5 LM for the frozen surface are shown in Fig. 12, in the same format as Fig. 5 for the flexible surface. We note that in this case the energetic ordering of the clusters is different, and new cluster shapes appear that are not present in the flexible case. In fact, the GM of the flexible case was not observed among the energetically lowest 50 minima of the frozen surface search. In contrast, a structure similar to LM 2 of the original search (Fig. 5 10 c)) and to the GM at cluster size 9 (Fig. 5 9 a)) now become the GM and LM 1. We can thus see that the flexibility of the support can have a large impact on the relative stability of the studied clusters and on the identified LM, as overall the energy difference are quite subtle and not very large. Even though some of the correct GM and LM structures might reappear in a subsequent second optimization in which the frozen substrate constraint is lifted, at least if there are no sizeable barriers, usually only the best intermediate structures are further refined, and thus the best overall structures may be missed in such a second step.

V Conclusions

We have performed a genetic algorithm global optimization search on the most stable structures of copper clusters with 4 to 10 atoms on the (10$\bar{1}$0) surface of zinc oxide employing a high-dimensional neural network potential constructed for this purpose using DFT reference data. This search, which essentially provides results of first-principles quality at a fraction of the computational costs, can be performed in very large surface cells demonstrating the need to use a flexible support to obtain meaningful low-energy structures. For the obtained global and local minima we have extracted important structural trends that aid in the understanding of the interaction between both materials. We have shown that a number of structural families are consistently present among the low energy configurations, often exchanging their ranking as the cluster size changes. A detailed analysis of various geometric properties of these clusters has been performed. In general, the GM aim to minimize their size and their interface area to volume ratio. The energy differences within a given cluster size are often very small, and thus clusters are expected to interconvert at moderately elevated temperatures. In general, we find that the clusters prefer to interact with the substrate through the surface oxygen atoms, with the footprint of the cluster limited to a rather small surface area evidencing the competition between support-cluster and cluster-cluster interaction. With respect to the interfacial structure of the clusters and the support, we observe that at the interface the clusters arrange themselves in a continuum between Cu (111) and (110) structures, already starting at cluster size 5, and becoming more prevalent at larger cluster sizes.

In summary, these results establish a complex picture of the interaction between the cluster and its support. The cluster needs to 1. interact with the support at specific sites, and within this constraint 2. minimize its energy. This kind of strong dependence between cluster morphology and support properties is already well known in the literature. Studies are available both for theoretical coarse models Henry (2005); Eckhoff, Schebarchov, and Wales (2017), and more specific atomistic simulations and experimental results Tauster (1987); Pan et al. (2017); Figueiredo et al. (2019), and is known as strong metal-support interaction.

VI Supplementary Material

As described in the text, supplementary material is provided:

Two tables with the parameters for the symmetry functions of the neural network potential.

A plot of predicted vs. DFT energy for the structures in our dataset for the chosen potential.

VII Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References