Abstract

We present a set of Gupta potentials fitted against highest-level ab initio data for interactions of the alkaline earth metals: beryllium, magnesium, and calcium. Reference potential energy curves have been computed for both pure and mixed dimers with the coupled cluster method including corrections for basic set incompleteness and relativistic effects. To demonstrate their usability for the efficient description of high-dimensional complex energy landscapes, the obtained potentials have been used for the global optimization of 38- and 42-atom clusters. Both pure and mixed compositions (binary and ternary clusters) were investigated. Distinctive trends in the structure of the latter are discussed.

1. Introduction

Metallic clusters have become over the years a subject of intense study, both theoretical as well as experimental [1]. Interest stems from the distinct properties they reveal when compared to the bulk phase and how these may change as a function of the cluster size. Different compositions (in binary, ternary, and higher mixtures) can also lead to new chemical and physical phenomena. Nanoalloys are a prime example of how both factors can be combined for material design and application in catalysis [2, 3]. The computational study of their structures is a challenging task for two interlacing reasons. On the one hand, the number of local minima is considered to scale exponentially with the cluster size, making the search for the global minimum NP-hard [4]. This property reflects back on all algorithms designed to explore the energy landscape of such systems. On the other hand, a suitable theoretical description of the interactions in play is required. It needs to be accurate enough to properly describe the energy landscape for a wide range of bonding patterns. It should also be amenable to computation, meaning that the computation of several hundred many-body interactions can be carried out in a sensible time frame. This is even more important since multiple thousands of these computations are required for a proper sampling of the energy landscape.

One of the most successful approaches to the study of metallic clusters has been the combination of fitted potentials with global optimization algorithms [5–8]. The former are usually obtained by fitting experimental data or electronic structure results to an analytic expression. The brute force use of quantum mechanical methods is impractical due to the computational cost, particularly linked to its scaling relative to the system size. Even semiempirical methods may be too costly as the prefactors are high enough to hinder a proper sampling of conformational space.

In this work, we have made use of correlated wave function methods to calculate the two-body interaction potential of alkaline earth metals (Be, Mg, and Ca). Emphasis has been placed on obtaining converged energy profiles relative to basic set, relativistic, and electronic correlation effects. The high-level reference data thus obtained was mapped to a two-body Gupta-type potential, [9] which in turn could be used to explore the structure of pure, binary, and ternary clusters. A few comments should be made about this choice of approach. First of all, it follows a bottom-to-top rationale that no information about nano- or macroscopic materials is used in the fit. It is purely based on first principles results that no empirical information (aside from the form of the chosen potential) has been included. This can certainly be seen as an advantage, since it allows us to improve the description in a systematic way. However, since the reference data has been computed with computationally demanding methods, it is not possible to benchmark the fit by repeating calculations for a selected test set of clusters. In fact, some of the terms included in the energy expression would be hard to obtain even for a 3-atom system. The advantages and disadvantages of our choice are later discussed in the text.

2. Methods and Techniques

Both for the cluster structure optimization as well as the potential fit, the OGOLEM framework for global optimization was used. Its features have been introduced in a series of publications [10–12]. Therefore, we will restrict ourselves to a brief discussion of the relevant features. The OGOLEM framework is loosely based on genetic algorithms as described in [13] replacing the generation-based scheme with the more efficient genetic pool scheme. While standard generation-based schemes feature serial bottlenecks at the end of every generation, a pool-based scheme removes this constraint through constant updates of a genetic population, allowing for a more efficient parallelization of the algorithm. As a side effect, elitism is a built-in feature of any genetic pool scheme, therefore removing the need to define additional criteria for it. Since the genetic pool contains all current solution candidates, parent individuals are chosen from it (father based on ranked fitness, mother randomly) and subject to the usual genetic operations: crossover and mutation. The crossover operator used for the global potential fit is a one-point genotype operator accompanied by a genotype mutation (probability 5%). For the cluster structure optimization, our implementation of a phenotype operator [11] is used, again accompanied by a genotype mutation (probability 5%). It should be noted here that no explicit exchange mutation (as, e.g., proposed in [6] and applied in [11]) was used. The phenotype implementation already includes some internal exchange which proved effective enough for lightly mixed clusters as targeted in this study.

In the case of cluster structure optimization, the solution candidates are then subject to a graph-based collision and dissociation detection. Should a candidate structure show either, it will be rejected and does not enter the subsequent local optimization. In the case of the potential fits, no such restriction is applied. Finally, it is attempted to add the fitter of the two locally optimized individuals to the genetic pool. This operation is only successful if it does not violate the fitness-based diversity criterion. After a defined number of these iterations, a converged solution pool is obtained, containing the global minimum candidate. In the case of cluster structure optimizations, such candidate is only accepted if four independent runs yield the same individual.

3. Global Fit of Potentials

All two-body interactions of beryllium, magnesium, and calcium have been fitted against highest-level ab initio data. The numerical data will be published elsewhere [27]. To obtain the highest possible accuracy at a still affordable computational footprint, different levels of theory based on wave function methods are combined as follows:𝐸inter=𝐸∞HF+Ξ”πΈβˆžCCSD(T)+Δ𝐸rel+Δ𝐸Q,(1) where 𝐸∞HF is the CBS[3,4,5]-extrapolated HF/aug-cc-pCVXZ [14–16] energy as proposed by Feller [17], Ξ”πΈβˆžCCSD(T) is the CBS[4,5]-extrapolated correlation energy using the CCSD(T)/aug-cc-pCVXZ (𝑋=𝑄,5) level of theory with the π‘‹βˆ’3 formula, Δ𝐸rel is a relativistic correction using a Douglas-Kroll Hamiltonian at the CCSD(T)/aug-cc-pCVTZ-DK [15] level of theory, and Δ𝐸Q is the quadruples contribution to the correlation energy obtained with CCSDT(Q)/aug-cc-pVTZ with the frozen core approximation in place. All calculations were performed with the Molpro2010.1 program package [18]. The CCSDT(Q) runs were carried out by the MRCC program [19, 20] interfaced to the latter.

The quality of this data set is high enough to reproduce the experimental dissociation energy of 11.1 kJ molβˆ’1 and equilibrium distance of 2.45 Å for the beryllium dimer [21] and can be expected to be of similar quality for the other interactions. Additionally, it provides a consistent data set for all pairs. For the latter property, the inclusion of all electrons in the Ξ”πΈβˆžCCSD(T) term calculation and the inclusion of relativistic effects are of particular importance.

The Gupta potential [9] used is of the regular form ξ‚Έξ‚΅π‘ŸπΈ(π‘Ž,𝑏)=𝐴(π‘Ž,𝑏)β‹…expβˆ’π‘(π‘Ž,𝑏)π‘Žπ‘π‘Ÿ0βˆ’ξƒŽ(π‘Ž,𝑏)βˆ’1ξ‚Άξ‚Ήπœ’(π‘Ž,𝑏)2ξ‚Έξ‚΅π‘Ÿβ‹…expβˆ’2β‹…π‘ž(π‘Ž,𝑏)π‘Žπ‘π‘Ÿ0,(π‘Ž,𝑏)βˆ’1ξ‚Άξ‚Ή(2) where π‘Ÿπ‘Žπ‘ is the distance between atoms π‘Ž and 𝑏, and 𝐴(π‘Ž,𝑏), 𝑝(π‘Ž,𝑏), π‘Ÿ0(π‘Ž,𝑏), πœ’(π‘Ž,𝑏), and π‘ž(π‘Ž,𝑏) are the parameters to be fitted against the reference.

Due to the rigid nature of the Gupta potential, a weighting of data points was necessary to guarantee a good fit. This weighting followed the rationale that an exact reproduction of the depth and position of the minimum is most important. A good reproduction of the attractive part of the potential was the second target, and less focus was placed on reproducing the repulsive part. We consider these to be reasonable design principles, reflecting the standard demands on potentials. Used weighting factors are tabulated in the supplementary information (see Supplementary Table in the Supplementary Material available online at doi: 10.1155/2012/648386).

The derived potentials are depicted in Figures 1 and 2 with the numerical values of the parameters to four digits precision available in Table 2. Perhaps one of the most striking features, upon inspection of the figures, is the difficulty in describing the weak-bonding regime. Some of the potential curves show a close to linear profile on approaching the minimum. This is the case for the Be-Be interaction and, less drastically, for the Be-Mg interactions. In the former case, a clear platteau is visible. Under the constraints of the potential form chosen, it is not possible to correctly reproduce this behavior without significantly affecting the description of the minimum. Nevertheless, all fitted potentials accurately describe the position and depths of the minimum correctly and are in overall good agreement with the reference. The Mg-Ca and Ca-Ca fits reproduce extremely well the reference data. Numerical information on the fitting quality can be obtained from Table 1. It should be noted though that the depth of the potential needs to be taken into account. The average deviation of 0.14 kJΒ·molβˆ’1 for the Ca-Ca interaction (minimum depth approximately 11 kJΒ·molβˆ’1) is less severe than the average deviation of 0.11 kJΒ·molβˆ’1 for the Mg-Mg interaction (minimum depth approximately 5 kJΒ·molβˆ’1). Further enhancements in the description would ultimately require another potential type, either another rigid potential more suitable for these interactions or a more flexible potential form. Both Morse potentials and LJ-type potentials were found to be unsuitable to overcome this principle problem. In a recent study by Li et al. [22], the Tang-Toennies potential model was used to fit experimental data of homogeneous alkaline earth dimers. The attractive part of the Be-Be interaction could not be perfectly described in this case either.

Further enhancement to the potential would also be possible through parametrization of three-body terms. These would have to be computed at a lower level of theory due to the large number of points needed and the size increase in the system. The computation of quadruple excitations is particularly costly and would be hard to perform in systems other than dimers. A possible approach would be to add an effective 3-body term in agreement with experimental structural data or by using simulation results at a lower level. Caution should be taken in computing such a term from three-atom systems for two reasons. First of all, it is expected that basic set superposition effects (BSSE/BSIE) can contaminate the potential. Most importantly, we note that many-body stabilization is overestimated when considering only 3-body interactions [25]. To illustrate the BSSE/BSIE problem, we compare in Figure 3 the energy profile for the Mg dimer computed at the CCSD(T)/aug-cc-pCVTZ level (CCSD(T)/AVTZ) and the energy obtained from the first two terms in (1) (CCSD(T)/CBS). The difference between the two sets of data is exclusively due to differences in the basic set. The use of a triple-zeta quality basic set leads to a clear overestimation of the well depth. The CCSD(T)/AVTZ level of theory predicts the equilibrium distance at 7.4 π‘Ž0 with a dissociation energy of 5.1 kJΒ·molβˆ’1 in contrast to the CCSD(T)/CBS prediction of 7.6 π‘Ž0 and 4.0 kJΒ·molβˆ’1, respectively. This amounts to an error of approximately 20% in the dissociation energy. If one were to estimate three-body terms with the triple-zeta basis, an overestimation will be expected. The basic functions of a third atom can contribute to the basic space of the neighboring dimer, resulting in a biased potential. Only close-to-CBS values could be used for correctly estimating 3-body contributions.

In general, we expect that the inclusion of 3-body terms should amount to an overall compression of the structure which would, in turn, induce local structural changes [25]. This could, however, be balanced by even higher-order terms in the many-body expansion. Work in this direction is underway.

4. Cluster Structure Optimization

To demonstrate the real-world applicability of the derived potentials, they have been used in the global optimization of medium-sized alkaline earth clusters. We focused on clusters of 38 alkaline earth atoms since this size typically exhibits the most interesting structural behaviour in the medium size regime [8]. To check whether the observed structural trends are specific to this cluster size, similar compositions in 42 atom clusters have been optimized. The structural data will be available from the Cambridge Cluster Database [26] after publication. All global minimum candidate structures are depicted in Figures 4 and 5.

The homogeneous clusters show icosahedral structural motifs. Depending on the atom in play, the structure varies slightly. While Be38, Ca38, and Ca42 possess mirror plane symmetry and seem to be magic numbers, Be42, Mg38, and Mg42 lack a number of atoms in defined positions, which is clear through visual inspection. It should be noted that no stable fcc structure could be located for any of the alkaline earth metals.

The same principle motifs hold true for the binary compositions. Common features are icosahedral substructures and real or pseudo mirror plane symmetry. Additional structural motifs can be observed for all binary clusters. First of all, a segregation of atom types can be observed in the form of the well-known core-shell structures [8] for all clusters containing beryllium. Beryllium forms an icosahedral core which can be easily explained with the potential profiles. The Be-Be interaction exhibits a deep and narrow minimum at a short distance. In contrary, the Mg-Mg and Ca-Ca interactions are both either not as deep (magnesium) or not as narrow (both magnesium and calcium). The formation of core-shell structures is also supported by the shape of the Be-Mg and Be-Ca potentials. In both cases, the minimum is located at longer distances than the Be-Be equilibrium distance and is not as deep as the Be-Be one. Obviously, the system must maximize the number of Be-Be contacts for an energetically low structure, which is only the case for a small icosahedral beryllium core.

A segregation of atom types can also be observed for the Mg/Ca binary compositions albeit not in the form of core-shell structures. Again, the potentials provide evidence for this behaviour. The Ca-Ca interactions possess a deeper minimum than the Mg-Ca interaction which in turn is slightly deeper than the Mg-Mg interaction. The system must therefore maximize the number of Ca-Ca contacts, followed by the number of Mg-Ca contacts. Since the equilibrium distance of the Ca-Ca is longer than the Mg-Ca and Mg-Mg one, a core-shell structure would require a very high Mg/Ca ratio. As can be seen from Figure 4(q), even a 29 : 9 ratio is not sufficiently high for such behaviour. In any other ratio, calcium forms the icosahedral backbone of the structure, with the magnesium atoms literally melting on that backbone, as can be seen, for example, in Figures 4(l) and 4(r). The resulting structures may probably be best described as Janus particles [8], possessing both magnesium and calcium character on the surface. Closely related is the ball-and-cup structure found, for example, in Figure 4(l).

The same design principles hold true when moving to ternary compositions. In the most simple case, when substituting single atoms, the binary cluster structure is slightly distorted but remains overall unchanged. This can be, for example, clearly seen in the transformation from the binary Be21 Mg21 (Figure 4(h)) to the ternary Be20 Mg20 Ca2 (Figure 5(j)) cluster. Once the composition contains more atoms of the third species, the cluster structure is again subject to the principle rules that have been formulated earlier. Beryllium forms a small icosahedral core, with magnesium and calcium segregating around it. This behaviour is most pronounced in the Be13 Mg15  Ca10 (Figure 5(e)) and Be6  Mg16  Ca16 (Figure 5(i)) cluster structures. In the earlier cluster, the beryllium core is large enough in comparison to the number of magnesium and calcium atoms to allow forming two half-shells around the core. In the latter, the core is small enough so that the calcium atoms form the shell and magnesium atoms remain at the surface. This ordering is due to the dissociation energy of the Be-Ca interaction being higher than the one of the Be-Mg interaction.

It is possible to conclude that alkaline earth clusters in the studied size regime seem to obey well-defined and rational building rules when using the Gupta model. A possible fault, and one which will be addressed in later work [27], is the problematic description of the beryllium atom. It is unclear how the deviations in the fit can influence the cluster structures. This, however, requires a more flexible functional form than the Gupta potential.

5. Conclusions

Gupta potentials for all bimetallic interactions involving beryllium, magnesium, and calcium are derived from highest-level ab initio data using global optimization techniques. All potentials reproduce the position and depths of the minimum correctly. The potentials have been subsequently used for the global optimization of medium-sized cluster structures, namely, up to ternary 42 atom clusters.

The structures obtained reveal several systematic trends. Clusters containing beryllium will form beryllium cores surrounded by a shell of the other atoms in play. Magnesium and calcium segregate, forming a calcium backbone with magnesium on the surface.

Acknowledgment

The authors gratefully acknowledge the financial support from the German Excellence Initiative through the Free Floater Research Group program of the University of GΓΆttingen.

Supplementary Materials

The electronic supplementary material contains the tabulated weighting factors used for the global fit of the Gupta-potentials.

  1. Supplementary Table