Study of structures and thermodynamics of CuNi nanoalloys using a new DFT-fitted atomistic potential

Shape, stability and chemical ordering patterns of CuNi nanoalloys are studied as a function of size, composition and temperature. A new parametrization of an atomistic potential for CuNi is developed on the basis of ab initio calculations. The potential is validated against experimental bulk properties, and ab initio results for nanoalloys of sizes up to 147 atoms and for surface alloys. The potential is used to determine the chemical ordering patterns of nanoparticles with diameters of up to 3 nm and different structural motifs (decahedra, truncated octahedra and icosahedra), both in the ground state and in a wide range of temperatures. The results show that the two elements do not intermix in the ground state, but there is a disordering towards solid–solution patterns in the core starting from room temperature. This order–disorder transition presents different characteristics in the icosahedral, decahedral and fcc nanoalloys.


Introduction
In recent years, nanoscale metallic alloys have become a widely studied topic at the frontier between condensed matter physics, chemistry and materials science. Such heterogeneous nanoparticles (NPs), often referred to as nanoalloys, have been shown to exhibit magnetic, catalytic and optical properties that can be different from those of bulk systems and of the elemental NPs. The size-, shape-and composition-dependence of these properties is not at all trivial but it is extremely important, since its comprehension would allow tailored design of nano-sized systems for specific applications in many technological fields, from industrial catalysis, to data storage and to medical research. 1 CuNi in particular is a system that has gained increasing interest for what concerns bulk surface alloys 2,3 and nanoalloys, 4 because of its catalytic properties. It has been shown that bulk surface CuNi alloys have a higher reactivity for carbon oxide reduction than elemental nickel or copper. They have also been used to catalyze the growth of mono-and bi-layer graphene with a fine control of thickness and uniformity. 3,5 The interaction between graphene layers and CuNi nanoparticles in nanocomposites has been shown to develop electrical and photoresponse properties making them promising candidates as future photodetectors. 6 Regarding biomedical applications, CuNi magnetic nanoparticles are employed in studies involving selective tissue hyperthermia. 4,7 CuNi nanoalloys exhibit remarkable magnetic properties. 8 Interest in small CuNi nanoalloys derived also from the possibility of producing chiral isomers. 9 From a thermodynamic point of view such a system is weakly miscible, presenting a bulk miscibility gap which extends up to 630 K. 10 The two elements present a small lattice mismatch and a small positive enthalpy of mixing. Such features suggest that, in clusters of small sizes, strain and pressure related driving forces, recently found to dominate segregation effects in NPs, 11 should be much weaker. Moreover CuNi NPs should present smaller miscibility gaps, shifted down to and below room temperature, drastically influencing the synthesis of such clusters. 12 Recent work on thermodynamic models applied to NPs in the 4-10 nm size range showed results confirming such finite-size effects on the phase diagram. 13 The study of the structure and thermodynamics of nanoalloys can in principle be performed at the ab initio level. However this approach is limited either to the study of a small number of isomers with sizes of several hundred atoms, 14 or to the full optimization of clusters of quite small sizes (well below 50 atoms). 15 Moreover, ab initio methods are out of question when simulating growth processes, whose time scale is of several microseconds at least. 16,17 On the other hand, thermodynamic modeling, which uses macroscopic concepts to study nanoscale objects, becomes less and less reliable when nanoparticle diameters decrease below B5 nm. For these reasons, the development of atomistic models is often necessary when studying nanoparticles and nanoalloys. Atomistic models can still provide reliable microscopic descriptions of these systems, being at the same time much less computationally expensive than ab initio calculations. Atomistic models thus allow a satisfactory exploration of the energy landscape for nanoparticles counting up to a few thousand atoms. They can be used to simulate phase transitions in such nanoparticles, and the growth on realistic time scales. The main problem with atomistic models is the validation of their accuracy, which must be performed case-by-case against experimental data, when available, and against ab initio calculations.
In this paper, a new atomistic model for the CuNi system is developed and tested, in the framework of a potential based on the second-moment approximation of the tight binding potential (SMTB), 18 also known as the Gupta potential. 19 SMTB potentials are nowadays used in the modeling of several nanoalloys, specifically to study problems of great complexity, such as global minimum searches, 20,21 phase transitions [22][23][24] and kinetic or growth processes, 17,25 since the cost of ab initio simulations is still too high for the extensive investigation required for these studies. However computational modeling of catalytic processes is not possible at the atomistic level, hence there is a need to develop a potential as compatible as possible with ab initio methods, to select structures and compositions for further ab initio studies.
Parameter sets for elemental copper and nickel have already been published in the work of Cleri and Rosato. 26 Such parameter sets were fitted to experimental bulk quantities, and developed with the aim of reproducing bulk properties as their main application. We have tried to develop a parametrization starting from Cu-Cu and Ni-Ni parameters of ref. 26, and fitting the heterogeneous Cu-Ni interactions on experimental solubility energies of impurities. The strategy of fitting on experimental values of solubility energies was successfully employed in the case of Ag-Ni, Ag-Cu and Ag-Pd heterogeneous interactions. 16,17,27 However, in the case of Cu-Ni, this strategy does not work very well for two reasons: (a) the experimental values of the solubility energies of impurities present quite large uncertainties, with a large difference between the values that can be found in the literature; 28,29 (b) the Ni potential of ref. 26 is not really satisfactory, 30 because it produces very high surface energies and the energy difference between hcp and fcc bulk phases is too large (0.05 eV per atom against the experimental value of 0.01 eV 31 ). This large difference between fcc and hcp phases would artificially enhance the energies of structural defects such as twin planes and stacking faults. These defects, while not important in bulk studies, can greatly influence the shape of clusters, where the crystalline symmetries are usually broken. Moreover, such a large energy difference destabilizes the hcp phase experimentally observed for Ni nanoparticles adsorbed on the MgO(100). 32,33 To avoid these problems, we have adopted the strategy of fitting all parameters (for both homogeneous and heterogeneous interactions) with the aid of DFT calculations. Therefore our potential will be a DFT-based one. After fitting the potential, we checked it against further DFT data (not used in the fitting), to the available experimental results and to recent thermodynamic calculations. 13 The potential is then used to study the chemical ordering in different structural motifs, both in the limit T -0 (by global optimization searches) and at finite temperatures (by Monte Carlo simulations) with the aim of singling out phase transitions. The results shown in the following section are thus a contribution to determine the nanoscale phase diagram in the size range between 10 2 and 10 3 atoms, in which DFT calculations are too cumbersome and thermodynamic modeling is not fully reliable.

DFT calculations
Calculations were performed within the Density Functional Theory approach implemented in the Quantum-Espresso package 34 with available ultrasoft pseudopotentials. 35 For exchange and correlation, the spin unrestricted Generalized Gradient Approximation in the Perdew-Burke-Ernzerhof 36 implementation was used. For the selection of the plane waves, energy cutoffs of 30 and 240 Ry were chosen to describe the wave function and the electronic density, respectively, and the sampling of the first Brillouin zone was performed at the G point.

Atomistic calculations
Atomistic calculations were performed by means of the SMTB potential, whose fitting procedure is discussed in Section 3.
Here we describe the calculations performed by using that potential.
2.2.1 Global optimization method. The optimal chemical ordering pattern within a given structural motif has been obtained via basin-hopping (BH) global optimization searches. This algorithm consists of Metropolis Monte Carlo simulations in which a local minimization is performed after each move. 37 Depending on the type of moves involved, this method can be used for both full structural optimization and for chemical ordering optimization within a given structural motif. Since in this work we are most interested in the latter task, in our simulations we use only exchange moves, where the identity of two atoms of different species is swapped. Both random and tailored moves were used. 14,38 In our simulations the structure of the cluster is fixed (apart from local relaxations), and the moves allow us to find the best homotops 39 (i.e. the best chemical ordering patterns within a given geometry) at the specified size and composition. For each size and composition we have performed at least 5 independent runs of 5 Â 10 5 moves.
2.2.2 Monte Carlo simulations. Since we are interested also in the stability of such chemical ordering at finite temperature, we have performed Monte Carlo (MC) simulations in a range of increasing temperatures. During these simulations only random exchanges have been used, and local minimization have been performed after each move to allow stress relaxation. It has been checked that this procedure reproduces the results of standard Monte Carlo simulations (without local minimization, but with random local displacements) in known cases. 40 Depending on the size of the cluster the simulations consisted in 1.0-2.5 Â 10 6 single Monte Carlo steps. The Monte Carlo simulations of bulk systems, which were used in the potential fitting procedure, have been performed under periodic boundary conditions in a box containing 1296 atoms.
2.2.3 Molecular dynamics simulations. Molecular dynamics (MD) simulations have been performed in the NVT canonical ensemble using a homemade code for a preliminary study of nanoparticle melting. The program uses the Velocity-Verlet integrator for the solution of the equations of motion coupled with an Andersen thermostat. The annealing rate has been set to 2 K ns À1 , with a time-step of 5 fs.

Potential fitting and validation
The SMTB potential used in the atomistic calculations has the following form; the potential energy of the system depends on the relative distances between atoms r ij and it is written as the sum of single-atom contributions E i : The parameters p ij , q ij , A ij , x ij , and r 0ij depend on the two atomic species of the pair ij only. As for r 0ij , we take the equilibrium distances in the bulk crystals at zero temperature for homogeneous pairs, and the average of the two equilibrium distances for heterogeneous pairs. r c is an appropriate cutoff distance. In this work we choose to put r c equal to the third-neighbor distance in the respective bulk solids for homogeneous pairs, whereas for heterogeneous pairs the r c is again the arithmetic average. The potential is then linked to zero at the fourthneighbor distance by a fifth-order polynomial function in such a way that the resulting potential is continuous with continuous derivatives for all interatomic distances. As stated in the Introduction, the parameters of the atomistic model are fully fitted by using DFT data. For each type of interaction (Cu-Cu, Ni-Ni, and Cu-Ni), there are four parameters (p, q, A, x) to be fitted. 26 For the homogeneous Cu-Cu and Ni-Ni interactions, the parameters are fitted on DFT data of lattice constant a, bulk cohesive energy per atom E c , bulk modulus B, and energy difference per atom DE hcp-fcc between hcp and fcc bulk phases. This fitting strategy has proved to be effective in developing atomistic potentials that are accurate for describing cluster structures, as discussed in ref. 30 and 33. The parameters for the mixed interactions were often chosen to be the mean values of those of the elemental interactions. While this has been proven to be rather successful in systems with a tendency to form solid solutions in the bulk, such as PdPt, 41 we judged that this procedure is a poor approximation for our system, which presents a miscibility gap. For this reason, only p and q have been taken as the arithmetic mean values of the pure parameters, while a and x have been fitted to reproduce the single impurity solubility energy E s of a Ni atom in a Cu bulk (and vice versa) as obtained in DFT simulations. The reference values are shown in Table 1. This procedure has been adopted for weakly mixing systems, such as Ag-Ni and Ag-Cu, with quite good results. 14,17 To reduce the computational cost of the simulations, a cutoff r c on interatomic distances has been adopted.

Ground-state validation
The potential has been validated through comparison to DFT results about clusters and surfaces. In the case of bulk crystal surfaces we checked pure-metal surface energies and the vertical relaxation of Ni layers (from 1 to 4) on Cu crystals and vice versa.
In the case of clusters we compared the energy ordering of different homotops for some different structures.
Atomistic and DFT results about the energies of (111) and (100) bulk crystal surfaces in Cu and Ni are shown in Table 2. The agreement between the two types of calculations is quite good.
Let us now consider the structural relaxation in the z-direction of Ni-layers on Cu crystals (and vice versa). In our DFT calculations, we fixed the position of the two lowest Cu(Ni) layers in a 5-layer slab to simulate an infinite bulk, allowing the relaxation of the other three Cu(Ni) layers. The further Ni(Cu) layers on top are all allowed to relax. For the atomistic simulations we just simulate a large slab of Cu(Ni) fcc bulk with the Ni(Cu) layers on top (2300-2600 atoms in total) and removed the periodic boundary conditions along the z axis. DFT and SMTB results are in good agreement, as can be seen in Fig. 1.
In the case of clusters, we considered three structures: the 55-atom truncated octahedron and icosahedron, and the 147-atom octahedron. For each structure, several different compositions were examined. Some of them are shown in Fig. 2. Since a thorough search of the Potential Energy Surface (PES) minima at the DFT level for these sizes is still computationally too expensive, the structures have been taken from a basin-hopping search at the atomistic level and afterwards relaxed by DFT. The comparison of the energy ordering shows that the preferential site for a single impurity is correctly identified in all cases, for both Ni impurity in Cu clusters and vice versa.
However, in general the agreement is not perfect, so two comments should be made. First, the energy differences between homotops in SMTB calculations are always smaller than those at the DFT level. Second, the agreement between the two methodologies, while quite good from a qualitative point of view, is subject to failure when considering the fine details: such properties as the copper segregation in the surface and the aggregation of nickel atoms in the core are very well reproduced, but the ordering of low-lying homotops is often not the same at DFT and SMTB levels.
The most evident failure is to reproduce the character of the incomplete copper capping in the case of very high nickel concentration. At the SMTB level the copper atoms occupy the less coordinated surface positions, that is all vertices first, then the edges, then the (100) facets and finally the (111) facets. This characteristic is common to all cluster shapes and sizes. There is a clear indication that this is not the case at the DFT level. For example, two representative homotops of the Cu 12 Ni 43 icosahedron have been considered. In the first homotop, twelve Cu atoms occupy all vertices. In the second homotop, they form a compact patch on one side of the surface. While at the SMTB the potential prefers the first homotop, the homotop with the compact patch is preferred at the DFT level.

Finite temperature validation
To check whether this parametrization can be trusted for predicting the correct chemical ordering patterns at finite temperature, we validated it against a well-known experimental property, that is, the transition in the bulk phase-diagram between the segregated phase and the solid solution. For the Cu-Ni systems several phase diagrams have been proposed, but there is a good agreement that at 50% Ni concentration the miscibility gap is up to a temperature of the order of 630 K, and at 25% Ni concentration this temperature is around 500 K. 10 The results of the MC simulations are shown in Fig. 3. The significant percentage of heterogeneous bonds at 0 K is due to the interface between the Cu and Ni segregated phases, and is just a spurious effect of the finite size of the simulation box. The transition temperature is given by the position of the inflexion point in the order parameter curve. For the 50% Ni concentration, the transition between a well-segregated system and a quasi-random mixing can be found around 550 K. At 25% Ni concentration, the transition temperature is B400 K. These results are in reasonable agreement with the experimental data on bulk alloys, even though the transitions are at somewhat lower temperatures. As we will see below, our transition temperatures for nanoparticles are however higher than those given by the atomistic models used in the literature. 12 The mixing enthalpy of the bulk alloy has been calculated by MC simulations at different compositions at a temperature of 973 K, which is well above the miscibility gap, but still in the solid range. This temperature has been considered in the literature, both in experiments and in simulations. 12,28,29

Optimal chemical ordering patterns
Here we determine the optimal chemical ordering according to our atomistic model for three magic-number nanoalloys: the 434-atom decahedron, the 561-atom icosahedron and the 586-atom fcc truncated octahedron. For each structure, several compositions are considered. A selection of results is shown in Fig. 5. The optimal chemical ordering patterns of these nanoalloys bear some resemblance to those of other non-mixing systems in which one of the elements has a lower surface energy and has a larger atomic radius (such as AgNi, AgCu, AgCo and AuCo 11 ), even if the difference in the size of nickel and copper atoms is quite small. For all sizes and structures the clusters show a complete surface segregation of copper atoms, given that the number of copper atoms is sufficient to cover all nickel atoms. Moreover, for all intermediate compositions, the structures show a very clear quasi-Janus chemical pattern, 42 with an offcentre nickel core with a monolayer coating of copper on one side. The main difference with AgNi, AgCu, AgCo and AuCo is that the lattice mismatch in CuNi is so small that the quasi-Janus clusters have almost no surface distortions. However, a small mismatch is sufficient to trigger the typical instabilities that lead to low-symmetry chemical ordering patterns 43 such as the quasi-Janus patterns, because of strain release. Quasi-Janus patterns have been found also in simulations of another system with a small mismatch, PtIr, but the driving force in that case was not related to strain release. 44 In the icosahedron, the best impurity location for a nickel atom is the central position, also for all other Ih magic-sizes, that is, 55, 147 and 309. This is expected since the Ni atom is (slightly) smaller and the centre of an icosahedron is its most compressed site. For the 561-atom cluster specifically, upon increasing the Ni concentration we first found a filling of the 55-atoms core, then a morphological instability appears (see Fig. 5, second row), similar to those already found in the literature, 11 in which the core takes a low-symmetry off-centre shape.
Nickel-rich clusters show incomplete surfaces where the copper atoms occupy positions on the surface in a precise pattern from the lowest coordinated sites, i.e. vertices, to the highest coordinated sites, that is the (111) facets. This kind of decoration, though, is not in agreement with the DFT results at very small sizes, where the tendency is found to be in favour of an incomplete but segregated surface capping, so that this atomistic result should therefore be taken with some care.
The subvertex site of the decahedral structure is the preferential location for a single nickel impurity. As much as 7 atoms are found to form a small crown around the five-fold axis just under the vertex. Larger inclusions then follow a very clear and unique pattern: if we consider the decahedron as divided in five slices (the original five tetrahedra composing the structure) by its 5-fold axis, the Ni atoms occupy one slice at a time, showing magical composition numbers of 63, 104, 145, 186 and 212 nickel atoms, as shown in Fig. 6. It should be noted however that at each compositions homotops that do not completely fill these slices are very close in energy to the optimal homotop, with DE B 0.03 eV.
In the truncated octahedra, nickel atoms always form subsurface (hence off-centre) aggregates, for all compositions. This pattern has been already found in several systems. 11,45 However in CuNi the difference with respect to AgCu, AgNi or AgCo is that the surface monolayer of copper atoms has neither holes nor deformations. 11

Thermodynamic stability of segregation
A complete segregation is not experimentally observed in several studies, for example in ref. 13, where some degree of intermixing is usually found. This intermixing may be caused by temperature effects. For this reason, in the following we study the temperature dependence of segregation by MC simulations. We consider again the 434-atoms Dh, the 561-atoms Ih and 586-atoms fcc TO structures. The temperature range considered is from 0 K (the ground states) up to 800 K, well above the mixing temperature of the bulk alloy, but well below the melting temperatures of the clusters (the latter will be briefly studied in Section 4.3). We decided to limit our investigation to clusters with 25% Ni composition. Since all surfaces are entirely made of copper, at this composition the core has roughly 50% nickel content and the effects of segregation and disordering are therefore more interesting.
All three structures show a monotonic increase of the percentage of heterogeneous bonds with temperature, indicating a smooth transition from complete segregation to quasi-random mixing. The transition temperature is identified as the inflexion point of the curve, and it is approximately 350 K for all the motifs considered. This transition temperature is higher than those found in simulations using other atomistic models. 12 It is interesting to comment that the percentage of nickel atoms in the surface is negligible even at the highest simulated temperature, meaning that the average atomic occupation for a surface atom is c surf Ni o 0.005. The value of the order parameter for completely random mixing has been therefore calculated as that of a completely randomly mixed core with a pure copper surface.
Although the three structures present similar features, some differences are present. The degree of segregation in the icosahedral structure is higher than that of the other two structures even at temperatures above the transition, as can be well seen in Fig. 7, where the average occupations for clusters at 600 K are shown. In contrast, the TO core present a perfect solid solution, and for the decahedral cluster only a weak nickel preferential occupation can be seen along the 5-fold axis. These effects are clearly related to the different internal pressure profiles.
While inside a TO the pressure is everywhere homogeneous (except for a slight compression of the subsurface atoms), the Ih motif has a highly compressed core and the Dh clusters present a somewhat compressed 5-fold symmetry axis. The strain release associated with the insertion of a smaller nickel atom in such location evidently has an effect on the average occupation even at temperature well above the transition.
The transition temperatures between ordered and intermixed patterns can be compared to the results of the work by Guisbiers et al., 13 in which phase diagrams for 4 nm-sized nanoparticles of several shapes were calculated by thermodynamic modeling, finding transitions at 200, 400 and 400 K for Dh, Ih and fcc shapes respectively. The results for Ih and fcc nanoalloys are in good agreement with our findings, considering that our clusters have a smaller size, circa 2.5 nm, and thus are expected to present somewhat lower transition temperatures. In the case of Dh, we do not find any significant temperature difference with the other motifs.

Melting
In order to evaluate the melting range of the structures, we have performed a preliminary MD study. For each structure, several molecular dynamics simulations have been made, starting from both intermixed and segregated seeds. At these sizes pure MD simulations are often not enough to obtain precise evaluations of the melting temperature from caloric curves, as shown in ref. 46, and more sophisticated algorithms need to be implemented. Nevertheless our results are sufficient to locate the melting range with some confidence, so that they allow us to make two comments: (a) all structures melt in the range around 1000-1100 K, which is in very good agreement with the results obtained by Guisbiers et al. by means of thermodynamic calculations 13 (b) all melting ranges are well above the highest temperature at which we performed MC runs at fixed structures for the study of the intermixing transition. Fig. 6 The shape of the nickel core for clusters of the 434-atom decahedral structure for five magic compositions corresponding to the filling of the five ''slices''. Ni atoms are shown as large gold spheres, copper as small red spheres. Fig. 7 Top row: average local nickel occupations for Dh, Ih and TO clusters at 25% Ni composition at 600 K. Bottom panel: percentage of heterogeneous bonds for the same clusters as a function of temperature: Dh (+), TO (Â) and Ih (*) clusters.

Conclusions
In this paper, a new DFT-based parametrization of the SMTB atomistic potential for CuNi has been developed. The new potential is found to describe well all the properties analysed, those regarding both infinite (bulk) and semi-infinite (surface) systems and those regarding finite-size clusters. Specifically, the potential is shown to correctly reproduce structural properties, such as surface vertical relaxation in thin films and impurity locations in clusters, and thermodynamic properties, such as miscibility gap and mixing enthalpy in bulk alloys. These results make this parametrization particularly suitable to study structural properties of nanoalloys at finite temperature, but also to perform thorough explorations of the nanoalloy energy landscape which can be used as the starting point of subsequent ab initio studies. This new parametrization has then been used to study the chemical ordering pattern for several cluster shapes, sizes and compositions, both at the level of ground-state search and in finite-temperature simulations. Concerning the optimal chemical ordering patterns, the system shows qualitative trends which are similar to those of other weakly miscible systems, but with some important differences. First of all, the CuNi system present very little lattice mismatch, thus showing quasi-Janus patterns with no surface reconstruction and almost no distortion, at variance in the cases of, e.g., AgCu or AgNi. As the temperature is increased, the chemical ordering pattern changes: while the surface segregation of copper atoms is almost complete even at quite high temperatures, the internal segregation is rapidly disordered, with a transition around 350 K for an overall composition of 25% in Ni. This finite-temperature effect, which is relevant even at room temperature, has thus to be taken into account in the synthesis of such nanoalloys.
The intermixing transition has different characteristics for the different structural motifs. This is related to the intrinsic pressure profile which derives from the non-crystalline nature of icosahedra and decahedra. The more compressed sites, such as those along the five-fold axis for the Dh and especially the core sites for the Ih, retain a higher nickel occupation even at temperatures well above the order-disorder transition.
In conclusion, we have shown that atomistic modelling can be useful in reproducing the nanoscale phase diagram of CuNi nanoalloys, thus being helpful in designing appropriate synthesis methods for these nanoalloys, which present a variety of core-shell or solid-solution patterns.