Research Article  Open Access
Johannes M. Dieterich, Sebastian Gerke, Ricardo A. Mata, "A FirstPrinciplesBased Potential for the Description of Alkaline Earth Metals", Journal of Atomic and Molecular Physics, vol. 2012, Article ID 648386, 8 pages, 2012. https://doi.org/10.1155/2012/648386
A FirstPrinciplesBased Potential for the Description of Alkaline Earth Metals
Abstract
We present a set of Gupta potentials fitted against highestlevel 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 highdimensional complex energy landscapes, the obtained potentials have been used for the global optimization of 38 and 42atom 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 NPhard [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 manybody 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 twobody 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 highlevel reference data thus obtained was mapped to a twobody Guptatype 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 bottomtotop 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 3atom 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 generationbased scheme with the more efficient genetic pool scheme. While standard generationbased schemes feature serial bottlenecks at the end of every generation, a poolbased 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 builtin 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 onepoint 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 graphbased 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 fitnessbased 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 twobody interactions of beryllium, magnesium, and calcium have been fitted against highestlevel 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/augccpCVXZ [14–16] energy as proposed by Feller [17], is the CBS[]extrapolated correlation energy using the CCSD(T)/augccpCVXZ () level of theory with the formula, is a relativistic correction using a DouglasKroll Hamiltonian at the CCSD(T)/augccpCVTZDK [15] level of theory, and is the quadruples contribution to the correlation energy obtained with CCSDT(Q)/augccpVTZ 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 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 weakbonding regime. Some of the potential curves show a close to linear profile on approaching the minimum. This is the case for the BeBe interaction and, less drastically, for the BeMg 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 MgCa and CaCa 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 CaCa interaction (minimum depth approximately 11 kJ·mol^{−1}) is less severe than the average deviation of 0.11 kJ·mol^{−1} for the MgMg 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 LJtype potentials were found to be unsuitable to overcome this principle problem. In a recent study by Li et al. [22], the TangToennies potential model was used to fit experimental data of homogeneous alkaline earth dimers. The attractive part of the BeBe interaction could not be perfectly described in this case either.


(a)
(b)
(c)
(a)
(b)
(c)
Further enhancement to the potential would also be possible through parametrization of threebody 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 3body 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 threeatom 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 manybody stabilization is overestimated when considering only 3body 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)/augccpCVTZ 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 triplezeta 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^{−1} in contrast to the CCSD(T)/CBS prediction of 7.6 and 4.0 kJ·mol^{−1}, respectively. This amounts to an error of approximately 20% in the dissociation energy. If one were to estimate threebody terms with the triplezeta 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 closetoCBS values could be used for correctly estimating 3body contributions.
In general, we expect that the inclusion of 3body 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 higherorder terms in the manybody expansion. Work in this direction is underway.
4. Cluster Structure Optimization
To demonstrate the realworld applicability of the derived potentials, they have been used in the global optimization of mediumsized 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
(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
The homogeneous clusters show icosahedral structural motifs. Depending on the atom in play, the structure varies slightly. While Be_{38}, Ca_{38}, and Ca_{42} possess mirror plane symmetry and seem to be magic numbers, Be_{42}, Mg_{38}, and Mg_{42} 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 wellknown coreshell structures [8] for all clusters containing beryllium. Beryllium forms an icosahedral core which can be easily explained with the potential profiles. The BeBe interaction exhibits a deep and narrow minimum at a short distance. In contrary, the MgMg and CaCa interactions are both either not as deep (magnesium) or not as narrow (both magnesium and calcium). The formation of coreshell structures is also supported by the shape of the BeMg and BeCa potentials. In both cases, the minimum is located at longer distances than the BeBe equilibrium distance and is not as deep as the BeBe one. Obviously, the system must maximize the number of BeBe 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 coreshell structures. Again, the potentials provide evidence for this behaviour. The CaCa interactions possess a deeper minimum than the MgCa interaction which in turn is slightly deeper than the MgMg interaction. The system must therefore maximize the number of CaCa contacts, followed by the number of MgCa contacts. Since the equilibrium distance of the CaCa is longer than the MgCa and MgMg one, a coreshell 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 ballandcup 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_{21} Mg_{21} (Figure 4(h)) to the ternary Be_{20} Mg_{20} Ca_{2} (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_{13} Mg_{15} Ca_{10} (Figure 5(e)) and Be_{6} Mg_{16} Ca_{16} (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 halfshells 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 BeCa interaction being higher than the one of the BeMg interaction.
It is possible to conclude that alkaline earth clusters in the studied size regime seem to obey welldefined 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 highestlevel 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 mediumsized 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 Guptapotentials.
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. View at: Publisher Site  Google Scholar
 B. F. G. Johnson, “From clusters to nanoparticles and catalysis,” Coordination Chemistry Reviews, vol. 190–192, pp. 1269–1285, 1999. View at: Publisher Site  Google Scholar
 C. L. Bracey, P. R. Ellis, and G. J. Hutchings, “Application of coppergold alloys in catalysis: current status and future perspectives,” Chemical Society Reviews, vol. 38, no. 8, pp. 2231–2243, 2009. View at: Publisher Site  Google Scholar
 L. T. Wille and J. Vennik, “Computational complexity of the groundstate determination of atomic clusters,” Journal of Physics A, vol. 18, no. 8, pp. L419–L422, 1985. View at: Publisher Site  Google Scholar
 B. Hartke, “Global optimization,” WIREs Computational Molecular Science, vol. 1, no. 6, pp. 879–887, 2011. View at: Google Scholar
 R. L. Johnston, “Evolving better nanoparticles: genetic algorithms for optimising cluster geometries,” Dalton Transactions, no. 22, pp. 4193–4207, 2003. View at: Google Scholar
 L. O. PazBorbón, T. V. MortimerJones, R. L. Johnston, A. PosadaAmarillas, G. Barcaro, and A. Fortunelli, “Structures and energetics of 98 atom PdPt nanoalloys: potential stability of the Leary tetrahedron for bimetallic nanoparticles,” Physical Chemistry Chemical Physics, vol. 9, no. 38, pp. 5202–5208, 2007. View at: Publisher Site  Google Scholar
 L. O. PazBorbón, R. L. Johnston, G. Barcaro, and A. Fortunelli, “Structural motifs, mixing, and segregation effects in 38atom binary clusters,” Journal of Chemical Physics, vol. 128, no. 13, Article ID 134517, 2008. View at: Publisher Site  Google Scholar
 R. P. Gupta, “Lattice relaxation at a metal surface,” Physical Review B, vol. 23, no. 12, pp. 6265–6270, 1981. View at: Publisher Site  Google Scholar
 J. M. Dieterich and B. Hartke, “OGOLEM: global cluster structure optimisation for arbitrary mixtures of flexible molecules. A multiscaling, objectoriented approach,” Molecular Physics, vol. 108, no. 34, pp. 279–291, 2010. View at: Publisher Site  Google Scholar
 J. M. Dieterich and B. Hartke, “Compositioninduced structural transitions in mixed LennardJones clusters: global reparametrization and optimization,” Journal of Computational Chemistry, vol. 32, no. 7, pp. 1377–1385, 2011. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. Corevalence basis sets for boron through neon,” Journal of Chemical Physics, vol. 103, no. 11, pp. 4572–4585, 1995. View at: Google Scholar
 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, corevalence, and scalar relativistic basis sets for Li, Be, Na, and Mg,” Theoretical Chemistry Accounts, vol. 128, no. 1, pp. 69–82, 2011. View at: Publisher Site  Google Scholar
 J. Koput and K. A. Peterson, “Ab initio potential energy surface and vibrationalrotational energy levels of X,” Journal of Physical Chemistry A, vol. 106, no. 41, pp. 9595–9599, 2002. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 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 stringbased 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 coupledcluster theory,” Journal of Chemical Physics, vol. 123, no. 21, Article ID 214105, 2005. View at: Publisher Site  Google Scholar
 J. M. Merritt, V. E. Bondybey, and M. C. Heaven, “Beryllium dimercaught in the act of bonding,” Science, vol. 324, no. 5934, pp. 1548–1551, 2009. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 Jmol: an opensource java viewer for chemical structures in 3d, http://www.jmol.org/.
 Povray—the persistence of vision raytracer, http://www.povray.org.
 E. BlaistenBarojas and S. N. Khanna, “Development of a firstprinciples manybody potential for beryllium,” Physical Review Letters, vol. 61, no. 13, pp. 1477–1480, 1988. View at: Publisher Site  Google Scholar
 The cambridge cluster database, http://wwwwales.ch.cam.ac.uk/CCD.html.
 J. M. Dieterich, S. Fischmann, and R. A. Mata, In preparation.
Copyright
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.