Journal of Atomic, Molecular, and Optical PhysicsVolume 2012 (2012), Article ID 648386, 8 pagesdoi:10.1155/2012/648386

Research Article

A First-Principles-Based Potential for the Description of Alkaline Earth Metals

Johannes M. Dieterich, Sebastian Gerke, and Ricardo A. Mata

Institut für Physikalische Chemie, Universität Göttingen, Tammannstrasse 6, 37077 Göttingen, Germany

Received 14 October 2011; Accepted 5 December 2011

Academic Editor: David Wales

Copyright © 2012 Johannes M. Dieterich et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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: where is the CBS[]-extrapolated HF/aug-cc-pCVXZ [14–16] energy as proposed by Feller [17], is the CBS[]-extrapolated correlation energy using the CCSD(T)/aug-cc-pCVXZ () level of theory with the formula, is a relativistic correction using a Douglas-Kroll Hamiltonian at the CCSD(T)/aug-cc-pCVTZ-DK [15] level of theory, and 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⁻¹ 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 term calculation and the inclusion of relativistic effects are of particular importance.

The Gupta potential [9] used is of the regular form where is the distance between atoms and , and , , , , 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⁻¹ for the Ca-Ca interaction (minimum depth approximately 11 kJ·mol⁻¹) is less severe than the average deviation of 0.11 kJ·mol⁻¹ for the Mg-Mg interaction (minimum depth approximately 5 kJ·mol⁻¹). 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.

Table 1: Quality of the fit for all alkaline earth interactions. Distance regime in Å and absolute and average deviation in kJ/mol. Absolute deviation includes weights; average deviation is weight-free.

Pair

Abs. dev.
Avg. dev.

Be-Be
2.010.0
81
47.43
0.59

Be-Mg
2.515.0
126
25.71
0.20

Be-Ca
2.715.0
124
25.80
0.21

Mg-Mg
3.015.0
121
13.79
0.11

Mg-Ca
3.215.0
119
11.42
0.10

Ca-Ca
3.515.0
116
16.60
0.14

Table 2: Numerical values to four digits precision for the fitted Gupta potentials. All values in atomic units.

Parameter
Be-Be
Be-Mg
Be-Ca
Mg-Mg
Mg-Ca
Ca-Ca

1.7943
2.1964
1.3955
2.7232
2.5646
0.8815

4.2656
2.0473
3.5205
2.8024
2.0231
3.3835

2.0323
1.4388
2.3759
1.8582
1.7852
3.9567

−0.4088
0.2066
0.0913
−0.0221
−0.1453
−0.5749

2.7536
1.1548
1.6065
0.7203
0.9789
2.8282

(a)

(b)

(c)

Figure 1: Derived Gupta potentials for Be-Be, Be-Mg, and Be-Ca interactions.

(a)

(b)

(c)

Figure 2: Derived Gupta potentials for Mg-Mg, Mg-Ca and Ca-Ca interactions.

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 with a dissociation energy of 5.1 kJ·mol⁻¹ in contrast to the CCSD(T)/CBS prediction of 7.6 and 4.0 kJ·mol⁻¹, 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.

Figure 3: Mg-Mg interaction energy curves obtained at the CCSD(T)/CBS and CCSD(T)/AVTZ levels of theory.

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.

(a) 38/0/0

(b) 42/0/0

(c) 0/38/0

(d) 0/42/0

(e) 0/0/38

(f) 0/0/42

(g) 19/19/0

(h) 21/21/0

(i) 19/0/19

(j) 21/0/21

(k) 0/19/19

(l) 0/21/21

(m) 29/9/0

(n) 9/29/0

(o) 29/0/9

(p) 9/0/29

(q) 0/29/9

(r) 0/9/29

Figure 4: Global minimum candidate structures of homogenous and binary clusters of the alkaline earth metals beryllium (indigo), magnesium (yellow), and calcium (red). The caption X/Y/Z denotes the number of beryllium atoms X, of magnesium atoms Y, and calcium atoms Z. All graphics are obtained with Jmol [23] and POV-ray [24].

(a) 13/13/12

(b) 13/12/13

(c) 12/13/13

(d) 14/14/14

(e) 13/15/10

(f) 13/6/19

(g) 13/9/16

(h) 16/6/16

(i) 6/16/16

(j) 20/20/2

(k) 20/2/20

(l) 2/20/20

(m) 19/18/1

(n) 18/19/1

(o) 19/1/18

(p) 18/1/19

(q) 1/19/18

(r) 1/18/19

Figure 5: Global minimum candidate structures of ternary clusters of the alkaline earth metals: beryllium (indigo), magnesium (yellow), and calcium (red). The caption X/Y/Z denotes the number of beryllium atoms X, of magnesium atoms Y, and calcium atoms Z. All graphics are obtained with Jmol [23] and POV-ray [24].

The homogeneous clusters show icosahedral structural motifs. Depending on the atom in play, the structure varies slightly. While Be₃₈, Ca₃₈, and Ca₄₂ possess mirror plane symmetry and seem to be magic numbers, Be₄₂, Mg₃₈, and Mg₄₂ 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 Be₂₁ Mg₂₁ (Figure 4(h)) to the ternary Be₂₀ Mg₂₀ Ca₂ (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 Be₁₃ Mg₁₅  Ca₁₀ (Figure 5(e)) and Be₆  Mg₁₆  Ca₁₆ (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.

References

R. Ferrando, J. Jellinek, and R. L. Johnston, “Nanoalloys: from theory to applications of alloy clusters and nanoparticles,” Chemical Reviews, vol. 108, no. 3, pp. 845–910, 2008.

B. F. G. Johnson, “From clusters to nanoparticles and catalysis,” Coordination Chemistry Reviews, vol. 190–192, pp. 1269–1285, 1999.

C. L. Bracey, P. R. Ellis, and G. J. Hutchings, “Application of copper-gold alloys in catalysis: current status and future perspectives,” Chemical Society Reviews, vol. 38, no. 8, pp. 2231–2243, 2009.

L. T. Wille and J. Vennik, “Computational complexity of the ground-state determination of atomic clusters,” Journal of Physics A, vol. 18, no. 8, pp. L419–L422, 1985.

B. Hartke, “Global optimization,” WIREs Computational Molecular Science, vol. 1, no. 6, pp. 879–887, 2011.

R. L. Johnston, “Evolving better nanoparticles: genetic algorithms for optimising cluster geometries,” Dalton Transactions, no. 22, pp. 4193–4207, 2003.

L. O. Paz-Borbón, T. V. Mortimer-Jones, R. L. Johnston, A. Posada-Amarillas, G. Barcaro, and A. Fortunelli, “Structures and energetics of 98 atom Pd-Pt nanoalloys: potential stability of the Leary tetrahedron for bimetallic nanoparticles,” Physical Chemistry Chemical Physics, vol. 9, no. 38, pp. 5202–5208, 2007.

L. O. Paz-Borbón, R. L. Johnston, G. Barcaro, and A. Fortunelli, “Structural motifs, mixing, and segregation effects in 38-atom binary clusters,” Journal of Chemical Physics, vol. 128, no. 13, Article ID 134517, 2008.

R. P. Gupta, “Lattice relaxation at a metal surface,” Physical Review B, vol. 23, no. 12, pp. 6265–6270, 1981.

J. M. Dieterich and B. Hartke, “OGOLEM: global cluster structure optimisation for arbitrary mixtures of flexible molecules. A multiscaling, object-oriented approach,” Molecular Physics, vol. 108, no. 3-4, pp. 279–291, 2010.

J. M. Dieterich and B. Hartke, “Composition-induced structural transitions in mixed Lennard-Jones clusters: global reparametrization and optimization,” Journal of Computational Chemistry, vol. 32, no. 7, pp. 1377–1385, 2011.

N. Carstensen, J. M. Dieterich, and B. Hartke, “Design of optimally switchable molecules by genetic algorithms,” Physical Chemistry Chemical Physics, vol. 13, no. 7, pp. 2903–2910, 2011.

D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Kluwer Academic Publishers, 1989.

D. E. Woon and T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. V. Core-valence basis sets for boron through neon,” Journal of Chemical Physics, vol. 103, no. 11, pp. 4572–4585, 1995.

B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, and A. K. Wilson, “Gaussian basis sets for use in correlated molecular calculations. VII. Valence, core-valence, and scalar relativistic basis sets for Li, Be, Na, and Mg,” Theoretical Chemistry Accounts, vol. 128, no. 1, pp. 69–82, 2011.

J. Koput and K. A. Peterson, “Ab initio potential energy surface and vibrational-rotational energy levels of X,” Journal of Physical Chemistry A, vol. 106, no. 41, pp. 9595–9599, 2002.

D. Feller, “Application of systematic sequences of wave functions to the water dimer,” Journal of Chemical Physics, vol. 96, no. 8, pp. 6104–6114, 1992.

H.-J. Werner, P. J. Knowles, R. Lindh, et al., Molpro, version 2010.1, a package of ab initio programs, 2010, http://www.molpro.net

Mrcc, a string-based quantum chemical program suite written by M. Kállay, see also M. Kállay, P. R. Surján, Journal of Chemical Physics , vol.115, pp. 2945, 2001, http://www.mrcc.hu

M. Kállay and J. Gauss, “Approximate treatment of higher excitations in coupled-cluster theory,” Journal of Chemical Physics, vol. 123, no. 21, Article ID 214105, 2005.

J. M. Merritt, V. E. Bondybey, and M. C. Heaven, “Beryllium dimer-caught in the act of bonding,” Science, vol. 324, no. 5934, pp. 1548–1551, 2009.

P. Li, J. Ren, N. Niu, and K. T. Tang, “Corresponding states principle for the alkaline earth dimers and the van der waals potential of Ba2,” Journal of Physical Chemistry A, vol. 115, no. 25, pp. 6927–6935, 2011.

Jmol: an open-source java viewer for chemical structures in 3d, http://www.jmol.org/

Pov-ray—the persistence of vision raytracer, http://www.povray.org

E. Blaisten-Barojas and S. N. Khanna, “Development of a first-principles many-body potential for beryllium,” Physical Review Letters, vol. 61, no. 13, pp. 1477–1480, 1988.

The cambridge cluster database, http://www-wales.ch.cam.ac.uk/CCD.html

J. M. Dieterich, S. Fischmann, and R. A. Mata, In preparation