We present a study on the structural, electronic, and magnetic properties of Fen(n  =  2  −  20) clusters by performing density functional tight binding (DFTB) calculations within a basin hopping (BH) global optimization search followed by density functional theory (DFT) investigations. The structures, total energies and total spin magnetic moments are calculated and compared with previously reported theoretical and experimental results. Two basis sets SDD with ECP and 6-31G** are employed in the DFT calculations together with BLYP GGA exchange-correlation functional. The results indicate that the offered BH-DFTB/DFT strategy collects all the global minima of which different minima have been reported in the previous studies by different groups. Small Fe clusters have three kinds of packing; icosahedral (Fe9−13), centered hexagonal antiprism (Fe14−17, Fe20), and truncated decahedral (Fe17(2), Fe18−19). It is obtained in a qualitative agreement with the time of flight mass spectra that the magic numbers for the small Fe clusters are 7, 13, 15, and 19 and with the collision induced dissociation experiments that the sizes 6, 7, 13, 15, and 19 are thermodynamically more stable than their neighboring sizes. The spin magnetic moment per atom of Fen(n = 2 − 20) clusters is between 2.4 and 3.6 μB for the most of the sizes. The antiferromagnetic coupling between the central and the surface atoms of the Fe13 icosahedron, which have already been reported by experimental and theoretical studies, is verified by our calculations as well. The quantitative disagreements between the calculations and measurements of the magnetic moments of the individual sizes are still to be resolved.

Clusters, as the aggregates of a few to thousands of atoms or molecules, constitute an important intermediate state between microscopic atoms and molecules and macroscopic condensed matter since the equilibrium arrangement of atoms and electronic structure of the clusters differ fundamentally from the bulk.1,2 Transition metal nanoclusters are being studied extensively by theoretical and experimental methods since their structural, electronic, and magnetic properties which may revolutionize many fields of nanotechnology such as heterogeneous catalysis and magnetic storage, are unique. In particular, Fe clusters are one of the most interesting examples due to the magnetic properties and importance in magnetic storage devices.

There are a number of studies on free small Fe clusters in the literature.3–44 Time-of-flight mass spectra of free Fen(n = 2 − 40) clusters produced by a laser vaporization source reveal the sequence of magic numbers n = 7, 13, 15, 19, 23.13 Fragmentation energies of mass selected Fe n 0 / + clusters were obtained for n = 2 − 19 using collision-induced dissociation where the largest bond energies are found for the Fe13 and Fe 13 + clusters.8 The Fe-Fe bond length was measured for Fe2 dimers trapped in argon and neon matrixes as (1.87 ± 0.13) Å3 and (2.02 ± 0.02) Å,5 respectively. Vibrational frequencies were obtained for gas-phase Fe2,4 and Fe 2 ,7 and for matrix-isolated Fe3.14 The geometry of a neutral Fe17 cluster was obtained in a salt.25 Stern-Gerlach experiments serve to determine the total magnetic moments of isolated transition metal clusters in the gas phase. Significantly enhanced magnetic moments in the range 2.5-5.4 μB per atom, much larger than the bulk value of 2.07 μB,45 have been observed for small Fen clusters in these experiments.6,9,10,12,18,34 Separations of the total magnetic moment into its spin and orbital constituents were recently done for Fe n + clusters using X-ray magnetic circular dichroism spectroscopy.36,42,44 It was found that the orbital magnetic moments are smaller by a factor of 2.5-20 than the spin magnetic moments.

The close relationship between structural stability and magnetic properties of Fen clusters motivated many theoretical studies as well.11,17,20,22–24,27–33,37–40,43 Density functional theory (DFT) with the local spin density approximations (LSDA), generalized gradient approximations (GGA) or hybrid methods and higher level ab initio calculations have been performed for small iron clusters. A recent theoretical discussion and the relevant literature on the ground state Fe2 dimer can be found in Ref. 32. An almost complete list of references on the DFT calculations of small iron clusters can be found in Ref. 38 including the ones considering the noncollinearity of local spin magnetic moments. After all of these studies, there are still important controversies about the ground state atomic arrangements, magnetic moments, noncollinear magnetism, and the stability spectrum of small Fe aggregates, which are still waiting to be solved. One of the main reasons for these controversies is the difficulty of obtaining the true ground state atomic arrangement of the clusters with relatively higher sizes.

As a nanocluster size is reduced to a few tens of atoms, the properties of the clusters are strongly dependent on cluster size and thus, a major focus in this context is to search their ground state structures. Although the experimental techniques such as chemisorption on cluster surfaces, Raman spectra, and mobility measurements that can lead to information about cluster structure have been progressed significantly, these are generally not sufficient to conclusively identify the geometries of size-specific clusters. In principle, quantum theoretical first principles methods such as DFT can be used to directly obtain different isomeric forms of the clusters. However, identifying the lowest energy structure of a cluster of even modest size using DFT is a computationally demanding task. It is a well known fact that the number of possible atomic configurations for a cluster increases exponentially with cluster size which makes the cluster global optimization problem a so called NP-hard problem. Since the computationally demanding first principle techniques like DFT are not practical for extensive surveys of the energy surfaces of the clusters other than having a few atoms, an alternate method that retains much of accuracy of DFT but is considerably more efficient is required. The DFT-based tight binding method (DFTB)46 is such a technique which combines the efficiency of the tight binding approach with the transferability of DFT. This method has been applied successfully to investigate a wide variety of clusters, molecules, and condensed matter systems.19,21,26,35,41

The aim of this work is to investigate the global minima of Fen(n = 2 − 20) clusters more extensively and efficiently in order to be able to get rid of the inconsistencies related to the atomic arrangement of these clusters. For this aim, we have developed a python code which combines the self-consistent charge DFTB and a basin-hopping (BH) global optimization algorithm.47 The BH Monte-Carlo (MC) minimization is one of the most reliable algorithms to search for the lowest energy structure of atomic clusters. It transforms the complex energy landscape into a collection of basins, and explores them by hopping, which is achieved by random MC moves and acceptance/rejection using the Metropolis criterion. After getting the lowest energy isomers of the clusters, we have calculated the total energies, magnetic moments, and vibrational frequencies of the first 10 isomers of each size through the spin-polarized DFT within the generalized gradient approximation (GGA). The transition states are separated from the local minima through the frequency calculations. Thus, the whole procedure is abbreviated as BH-DFTB/DFT calculations. Noncollinear effects are not considered in the present work since it is not the primary focus of the study and it is reported that geometric distortions restore the collinearity of the local spin magnetic moments for Fen(n = 2 − 15) clusters.23 All the iron cluster structures reported previously in the literature were considered in this work to resolve some of the contradictory findings regarding possible ground state structures.

The computational procedure illustrated in Fig. 1 schematically consists of two phases: the BH-DFTB phase and the DFT phase. The aim of the BH-DFTB phase is to search the potential energy surface of the iron clusters as efficient and complete as possible, which is necessary to get the true global minimum. In order to achieve this aim, we first compose a pool of initial configurations of the iron clusters which are used as the seeds of several BH global optimization runs for the search of the potential energy surface described by the DFTB method. The initial configurations included in the pool are of three classes: The first class of coordinates are obtained from the previous classical calculations such as the ones reported by Doye et al.48 and by Elliott et al.49 where the interatomic forces are described by the Morse and Finnes-Sinclair potentials, respectively. As a second class of initial configurations, the previously reported DFT results28,38 have been included in the pool. Finally as the third class of coordinates, we have prepared a number of initial configurations having different molecular point group symmetries and bulk like crystallographic symmetries. In the BH-DFTB phase, we have performed several BH MC minimization runs starting from each of these initial coordinates, where a thousand MC steps have been taken in each minimization runs. Since different total spin states may form different basins on the potential energy landscape, different spin states have been checked before starting any of the MC runs by performing local minimizations for each initial configuration. The MC runs are performed for the optimized total spin magnetic moments.

FIG. 1.

Schema of the procedure employed in the present study.

FIG. 1.

Schema of the procedure employed in the present study.

Close modal

In the traditional tight binding methods, the matrix elements of the Hamiltonian are treated empirically, whereas within DFTB these are calculated with DFT derived basis functions and atomic potentials. In addition, DFTB method can take into account possible charge transfers between different atoms and can include spin polarization. A single point energy calculation in DFTB is at least two orders of magnitude faster than in DFT, depending on cluster size,19 which allows a more thorough search of cluster configuration space than would be practical in DFT. We have used the Fe-Fe interaction parameters which have been calibrated for the spin-polarized self-consistent-charge DFTB method by Zheng et al.50 These parameters work quite well for the geometries of organic molecules including O, N, C, H with the first row transition metal elements Sc, Ti, Fe, Co, and Ni. In this first phase of the calculations, the default convergence criteria of the DFTB+ code have been employed.46 

After getting all the low-lying isomers of the iron clusters at the end of the BH-DFTB calculations, we started the second (DFT) phase of our procedure. In this phase, the first 10 isomers of each size of the iron clusters having the lowest total DFTB energies are chosen and local DFT geometry optimizations imposing no symmetry constraints are performed. Their DFT energies, spin multiplicities, and vibrational frequencies at the optimized geometries are calculated with a GGA xc-functional. Although DFTB is a good and fast approximation to DFT, the potential energy surfaces of iron clusters described by DFTB are not exactly the same with the ones described by DFT. In Table I, we present the orders of iron clusters isomers obtained by DFTB and DFT calculations, where the number of isomers increases with the decrease of the binding energies. The table indicates that the orders of isomers change significantly, even though all of the chosen DFTB isomers are still the low-lying structures of DFT calculations. We, then, compared our DFT calculations with the previous DFT results reported in the literature by reoptimizing them with the basis set and xc-functional chosen in the present study. During this comparison, we never got a case where a previously reported structure has a lower energy than the one we obtained with the described procedure. By the end of this check, we finally identified the lowest energy geometric, electronic, and magnetic structures of iron clusters, which are reported in Section III.

TABLE I.

Orders of isomers for DFTB and DFT calculations.

n DFTB DFT n DFTB DFT n DFTB DFT n DFTB DFT
13  17 
13  17 
10  13  17 
13  17 
13  17  10 
10  13  17 
13  17 
13  17 
13  10  17 
10  10  13  10  17  10 
10  14  18 
10  14  18  10 
10  14  18 
10  14  18 
10  14  10  18 
10  14  18 
10  14  18 
10  14  18 
10  10  14  18 
    10  10  14  10  18  10 
11  10  15  19 
11  15  19 
11  15  19 
10  11  15  19 
11  15  19 
11  15  19 
11  15  10  19 
11  15  19 
11  15  19 
10  11  10  15  10  19  10  10 
10  12  16  20 
12  16  20 
12  16  20 
12  10  16  20 
12  16  20  10 
12  16  20 
12  16  20 
12  16  20 
12  16  10  20 
10  12  10  16  10  20  10 
n DFTB DFT n DFTB DFT n DFTB DFT n DFTB DFT
13  17 
13  17 
10  13  17 
13  17 
13  17  10 
10  13  17 
13  17 
13  17 
13  10  17 
10  10  13  10  17  10 
10  14  18 
10  14  18  10 
10  14  18 
10  14  18 
10  14  10  18 
10  14  18 
10  14  18 
10  14  18 
10  10  14  18 
    10  10  14  10  18  10 
11  10  15  19 
11  15  19 
11  15  19 
10  11  15  19 
11  15  19 
11  15  19 
11  15  10  19 
11  15  19 
11  15  19 
10  11  10  15  10  19  10  10 
10  12  16  20 
12  16  20 
12  16  20 
12  10  16  20 
12  16  20  10 
12  16  20 
12  16  20 
12  16  20 
12  16  10  20 
10  12  10  16  10  20  10 

NWChem 6.1 program package51 has been used to perform geometry optimizations and total energy calculations by DFT. Two basis sets have been employed: the first one is the Stuttgart relativistic, small core 199752 basis sets (SDD) and effective core potentials (ECP), where the outer most 16 electrons of free Fe atom (3s23p63d64s2) are treated explicitly, and the corresponding contracted Gaussian basis functions are (6s5p3d1f). The second one is 6-31G** all electron Gaussian basis set53 with contracted basis functions of 5s4p2d1f. Becke’s GGA exchange functional54 and Lee-Yang-Parr correlation functional55 (BLYP) has been employed in the calculations. Justification of the choice of the xc-functional and the basis set is discussed in the next section by giving a comparison between the results obtained for the Fe2 dimer. In the calculations with SDD basis set, the default convergence criteria, which are 1 x 10−6 Hartree for energy, and 5 x 10−4 Hartree/a0 for energy gradient, are used, while the energy criteria for the calculations with 6-31G** basis set is relaxed to 1 x 10−5 Hartree.

For the optimum choice of the xc-functional and the basis set to be used in the DFT phase, we have performed a number of calculations with different functionals and basis sets for the magnetic moment, binding energy, bong length, and vibrational frequency of Fe2 dimer and compared the obtained results with the available experimental data in Table II. Five basis sets with ECPs (LANL2DZ,56 LANL2TZ,56 CRENBL,57 SDD,52 QZVP58) plus the 6-31G** all electron basis set53 and four xc functionals (two GGA: BLYP,54,55 MPWPW91,59–61 a meta-GGA: TPSS,62 and a hybrid: PBE063) have been tried.

TABLE II.

The comparison of the experimental data with the calculated results by employing the chosen xc-functionals and basis sets for the properties of Fe2 dimer.

Magnetic Moment (μB) Binding Energy (eV)
Experimental6 = 6.5 ± 1 Experimental8 = 1.14 ± 0.10
Basis Sets BLYP MPWPW91 TPSS PBE0 BLYP MPWPW91 TPSS PBE0
LANL2DZ ECP (3s3p2d)  1.87  1.91  1.55  0.31 
LANL2TZ ECP (5s5p3d)  1.51  1.71  1.63  -0.09 
CRENBL ECP (7s6p6d)  1.44  1.60  1.47  1.10 
SDD ECP (6s5p3d1f)  1.51  1.69  1.58  0.91 
QZVP ECP (11s6p5d3f1g )  2.27  2.33  1.92  0.98 
6-31G** All Elec. (5s4p2d1f)  1.69  2.10  1.75  -0.31 
  Bond Length (Å)   
  Experimental3 = 1.87 ± 0.13  Vibrational freq. (cm−1
  Experimental5 = 2.02 ± 0.02  Experimental4 = 299.6 
Basis Sets  BLYP  MPWPW91  TPSS  PBE0  BLYP  MPWPW91  TPSS  PBE0 
LANL2DZ ECP (3s3p2d)  2.04  2.03  2.20  2.14  406.0  435.6  229.1  391.9 
LANL2TZ ECP (5s5p3d)  2.04  2.02  2.29  2.13  395.4  422.5  293.8  386.6 
CRENBL ECP (7s6p6d)  2.02  2.01  2.00  2.09  404.5  436.4  418.0  397.9 
SDD ECP (6s5p3d1f)  2.02  2.01  2.28  2.11  393.0  427.2  284.7  389.4 
QZVP ECP (11s6p5d3f1g )  2.02  2.00  2.00  2.09  447.2  576.6  410.6  549.3 
6-31G** All Elec. (5s4p2d1f)  2.03  2.02  2.01  1.99  416.1  533.1  435,0  455.4 
Magnetic Moment (μB) Binding Energy (eV)
Experimental6 = 6.5 ± 1 Experimental8 = 1.14 ± 0.10
Basis Sets BLYP MPWPW91 TPSS PBE0 BLYP MPWPW91 TPSS PBE0
LANL2DZ ECP (3s3p2d)  1.87  1.91  1.55  0.31 
LANL2TZ ECP (5s5p3d)  1.51  1.71  1.63  -0.09 
CRENBL ECP (7s6p6d)  1.44  1.60  1.47  1.10 
SDD ECP (6s5p3d1f)  1.51  1.69  1.58  0.91 
QZVP ECP (11s6p5d3f1g )  2.27  2.33  1.92  0.98 
6-31G** All Elec. (5s4p2d1f)  1.69  2.10  1.75  -0.31 
  Bond Length (Å)   
  Experimental3 = 1.87 ± 0.13  Vibrational freq. (cm−1
  Experimental5 = 2.02 ± 0.02  Experimental4 = 299.6 
Basis Sets  BLYP  MPWPW91  TPSS  PBE0  BLYP  MPWPW91  TPSS  PBE0 
LANL2DZ ECP (3s3p2d)  2.04  2.03  2.20  2.14  406.0  435.6  229.1  391.9 
LANL2TZ ECP (5s5p3d)  2.04  2.02  2.29  2.13  395.4  422.5  293.8  386.6 
CRENBL ECP (7s6p6d)  2.02  2.01  2.00  2.09  404.5  436.4  418.0  397.9 
SDD ECP (6s5p3d1f)  2.02  2.01  2.28  2.11  393.0  427.2  284.7  389.4 
QZVP ECP (11s6p5d3f1g )  2.02  2.00  2.00  2.09  447.2  576.6  410.6  549.3 
6-31G** All Elec. (5s4p2d1f)  2.03  2.02  2.01  1.99  416.1  533.1  435,0  455.4 

The accurate description of the electronic structure of transition metal dimers is not a straight forward task as one can see in the study of Angeli and Cimiraglia,32 who have reported that the 9 Σ g and 7Δu states of the iron dimer are the most probable candidates for the assignment of the ground state according to their multi reference perturbation theory calculations. This difficulty manifests itself in the present study in the way that while GGA functionals, when they are employed with ECP basis sets, predict the total spin magnetic moment of the iron dimer as 6 μB, the hybrid functional PBE0 estimates it as 8 μB for all the basis sets used here. The all electron 6-31G** basis set calculations result in 8 μB without depending on the xc-functional choice. On the other hand, the total spin magnetic moment obtained by the meta-GGA functional TPSS depends on the choice of the basis set. The result shuttles between 6 and 8 μB as the number of gaussian functions included in the calculations increases, where the reason of this shuttling is out of the scope of this study.

In accordance with the previous theoretical calculations64 the hybrid functional PBE0 tends to overestimate bond distance with the exception of all electron case and thus underestimate binding energy. Both of the GGA functionals BLYP and MPWPW91 produce better bond distances than the others two, where the choice of basis set plays a minor role since the bond distance decreases slightly with the increasing number of gaussian functions used in the basis sets. A similar oscillatory behavior of TPSS observed in the spin magnetic moment can be found in the calculated bond distances as well. On the other hand, the binding energies obtained by the GGA functionals are not in better accord with experiment for they significantly overestimate the iron dimer binding energy. Interestingly, the largest basis set employed in the present study (QZVP) and the all electron basis set even enhance the deviation from the reference value, which is in contrast with the expectations. We do not find any advantage of using TPSS against GGA functionals in terms of the dimer binding energy. Although, PBE0 produces closer results to the experimentally measured value, it is clear that this is due to the overestimated bond distances. Vibrational frequency calculations do not present a solid reason for the choice of a specific combination of basis sets and xc-functionals either. Nearly all of the frequency calculations overestimate the frequency substantially with a few exceptions. Thus, we have decided to employ the BLYP functional together with SDD and 6-31G** basis sets in the DFT phase of our calculations. The reason for the employment of the all electron basis set in addition to the SDD basis set will be more clear, when we discuss the binding energies and magnetic moments in the next sections.

We present the optimized geometries of the lowest energy structures and a few low-lying isomers of Fen(n = 3 − 20) clusters along with some typical bound distances obtained by the BLYP/SDD choice of the functional and the basis set in Fig. 2 and Fig. 3, respectively. The point group symmetries, binding energies (Eb), total spin magnetic moments (μs), average bond distances (da), HOMO-LUMO gaps for the spin down electrons (Egap), and the lowest (wl) and the highest (wh) vibrational frequencies of these structures are given in Table III for both of the choices of BLYP/SDD and BLYP/6-31G**. The morphologies of the ground state arrangements of the clusters obtained by BLYP/6-31G** is the same with the ones given in Fig. 2, but the bond distances differ. The average bond distances of BLYP/6-31G** calculations can be found in Table III. The point group symmetry of each cluster was determined using a tolerance of 0.1 Å between the atomic positions in the computed structure and the symmetrize reference structure. Many of our ground state geometries experience stabilizing Jahn-Teller distortions that reduce their symmetry. The lowest energy structure of Fe3 is obtained as a scalene triangle, although in many of the previous calculations isosceles triangles were identified as the ground state geometry.20,27,28,37,43 Kim et al.43 reported that the isosceles triangle is only 0.0026 eV/atom lower in energy than the scalene triangle. In our calculations the isosceles triangle is 0.06 eV/atom (with BLYP/SDD) higher in energy than the scalene one.

FIG. 2.

Geometrical configurations corresponding to the lowest total energy states of Fen(n = 3 − 20) clusters. Bound distances are given in Å.

FIG. 2.

Geometrical configurations corresponding to the lowest total energy states of Fen(n = 3 − 20) clusters. Bound distances are given in Å.

Close modal
FIG. 3.

Geometrical configurations corresponding to the low-lying isomers of Fen(n = 3 − 20) clusters. Bound distances are given in Å.

FIG. 3.

Geometrical configurations corresponding to the low-lying isomers of Fen(n = 3 − 20) clusters. Bound distances are given in Å.

Close modal
TABLE III.

Point group symmetries, binding energies (Eb in eV/atom), total spin magnetic moments (μs in μB/atom), average bond distances (da in Å), HOMO-LUMO gaps for spin down electrons (Egap in eV), the lowest (wl) and the highest (wh) vibrational frequencies (in cm−1) for the obtained lowest energy structures of Fen(n = 2 − 20) clusters. Whenever they differ, BLYP/SDD and BLYP/6-31G** calculations are given as the first and the second values, respectively.

n Symmetry Eb μs da Egap wl- wh
Dh  0.76 / 0.43  3.00 / 4.00  2.02 / 2.03  0.62 / 1.02  393.0 / 416.1 
Cs / C2v  1.12 / 1.23  3.33  2.28 / 2.26  0.78 / 0.26  9.2 / 155.5 - 335.5 / 381.3 
C3v / C3  1.43 / 1.75  3.50  2.38 / 2.29  0.71 / 0.40  46.5 / 21.5 - 344.4 / 406.8 
C2v / Cs  1.70 / 2.15  3.20  2.39 / 2.30  0.59 / 0.26  52.6 / 46.7 - 341.9 / 406.3 
D4h  1.93 / 2.44  3.33  2.46 / 2.37  0.48 / 0.51  104.3 / 72.3 - 316.0 / 394.1 
C1 / Cs  2.05 / 2.62  3.14  2.44 / 2.35  0.73 / 0.35  73.2 / 102.2 - 315.6 / 386.5 
C2v  2.10 / 2.70  3.00  2.45 / 2.35  0.46 / 0.29  104.0 / 123.4 - 318.4 / 368.9 
C2v  2.13 / 2.79  2.89  2.44 / 2.35  0.55 / 0.37  86.7.1 / 76.7 - 324.9 / 381.5 
10  C1 / C3  2.17 / 2.87  3.00  2.46 / 2.37  0.45 / 0.40  -70.4 / -7.4 - 319.5 / 252.7 
11  C1 / C2v  2.21 / 2.92  3.09  2.48 / 2.39  0.53 / 0.26  45.0a / 30.9 - 345.5a / 399.0 
12  Cs  2.25 / 2.97  3.17  2.50 / 2.41  0.38 / 0.28  40.0 / 20.1 - 326.0 / 395.2 
13  C1  2.29 / 3.06  3.38 / 2.62  2.52 / 2.40  0.34 / 0.20  43.8 / -75.2 - 339.2 / 375.8 
14  C1 / C6  2.34 / 3.09  3.29  2.52 / 2.44  0.40 / 0.33  59.9 / 62.8 - 279.4 / 328.1 
15  D3  2.41 / 3.19  3.33 / 3.02  2.53 / 2.44  0.49 / 0.39  64.6 / 95.1 - 316.4 / 363.5 
16  D3h  2.43 / 3.22  3.25  2.54 / 2.45  0.37 / 0.34  24.1 / 65.3 - 286.0 / 328.2 
17  S4 / Cs  2.46 / 3.23  3.18  2.54 / 2.46  0.28 / 0.23  26.7a / -48.7 - 276.8a / 328.7 
18  Cs / C2v  2.47 / 3.27  3.11  2.53 / 2.47  0.23 / 0.25  68.5 / 82.4 - 276.4 / 319.7 
19  C1  2.48 / 3.30  3.05  2.54 / 2.45  0.44 / 0.50  54.6 / 80.1 - 284.5 / 327.0 
20  C1  2.48 / 3.22  3.00  2.53 / 2.44  0.21 / 0.19  69.1a / 42.6 - 367.8a / 386.7 
n Symmetry Eb μs da Egap wl- wh
Dh  0.76 / 0.43  3.00 / 4.00  2.02 / 2.03  0.62 / 1.02  393.0 / 416.1 
Cs / C2v  1.12 / 1.23  3.33  2.28 / 2.26  0.78 / 0.26  9.2 / 155.5 - 335.5 / 381.3 
C3v / C3  1.43 / 1.75  3.50  2.38 / 2.29  0.71 / 0.40  46.5 / 21.5 - 344.4 / 406.8 
C2v / Cs  1.70 / 2.15  3.20  2.39 / 2.30  0.59 / 0.26  52.6 / 46.7 - 341.9 / 406.3 
D4h  1.93 / 2.44  3.33  2.46 / 2.37  0.48 / 0.51  104.3 / 72.3 - 316.0 / 394.1 
C1 / Cs  2.05 / 2.62  3.14  2.44 / 2.35  0.73 / 0.35  73.2 / 102.2 - 315.6 / 386.5 
C2v  2.10 / 2.70  3.00  2.45 / 2.35  0.46 / 0.29  104.0 / 123.4 - 318.4 / 368.9 
C2v  2.13 / 2.79  2.89  2.44 / 2.35  0.55 / 0.37  86.7.1 / 76.7 - 324.9 / 381.5 
10  C1 / C3  2.17 / 2.87  3.00  2.46 / 2.37  0.45 / 0.40  -70.4 / -7.4 - 319.5 / 252.7 
11  C1 / C2v  2.21 / 2.92  3.09  2.48 / 2.39  0.53 / 0.26  45.0a / 30.9 - 345.5a / 399.0 
12  Cs  2.25 / 2.97  3.17  2.50 / 2.41  0.38 / 0.28  40.0 / 20.1 - 326.0 / 395.2 
13  C1  2.29 / 3.06  3.38 / 2.62  2.52 / 2.40  0.34 / 0.20  43.8 / -75.2 - 339.2 / 375.8 
14  C1 / C6  2.34 / 3.09  3.29  2.52 / 2.44  0.40 / 0.33  59.9 / 62.8 - 279.4 / 328.1 
15  D3  2.41 / 3.19  3.33 / 3.02  2.53 / 2.44  0.49 / 0.39  64.6 / 95.1 - 316.4 / 363.5 
16  D3h  2.43 / 3.22  3.25  2.54 / 2.45  0.37 / 0.34  24.1 / 65.3 - 286.0 / 328.2 
17  S4 / Cs  2.46 / 3.23  3.18  2.54 / 2.46  0.28 / 0.23  26.7a / -48.7 - 276.8a / 328.7 
18  Cs / C2v  2.47 / 3.27  3.11  2.53 / 2.47  0.23 / 0.25  68.5 / 82.4 - 276.4 / 319.7 
19  C1  2.48 / 3.30  3.05  2.54 / 2.45  0.44 / 0.50  54.6 / 80.1 - 284.5 / 327.0 
20  C1  2.48 / 3.22  3.00  2.53 / 2.44  0.21 / 0.19  69.1a / 42.6 - 367.8a / 386.7 
a

PBE0/SDD.

There is an agreement in the literature17,20,21,23,27,28,37,38,40,43 that for the iron clusters with n = 4 − 8, tetrahedron, trigonal bipyramid, octahedron, (distorted) pentagonal bipyramid (PBP), and bidisphenoid (bicapped octahedron) are their lowest energy structures, respectively, although the point group symmetries of these structures change from study to study. Our findings are consistent with the literature morphologically in this size range. All the other properties presented in Table III can be compared with the available literature. We found that the most stable isomer of Fe9 is a bicapped PBP in C2v symmetry, which is in agreement with the findings of Ma et al.,28 Cervantes-Salguero and Seminario,37 Yuan et al.,40 and Kim et al.43 but differs from the other DFT calculations,17,21,23,38 where a tri-capped trigonal prism (D3h) is reported as the lowest energy structure. The tri-capped trigonal prism has 12.5 meV/atom high energy than the bicapped PBP in our calculations.

We identified the most stable Fe10 structure as a tri-capped PBP structure, which is in agreement with Ma et al.,28 Cervantes-Salguero and Seminario,37 Yuan et al.,40 and Kim et al.,43 but differs from the other previous DFT calculations by Dieguez et al.,17 Köhler et al.,21 Rollmann et al.,23 and Gutsev et al.,38 who suggested a bicapped square antiprism. The energy difference between the tri-capped PBP and bicapped square antiprism is 20.8 meV/atom (with BLYP/SDD). Our finding that the lowest energy Fe11 structure is a tetra-capped PBP with C1 symmetry is consistent with all of the previously reported DFT calculations17,21,28,38,40,43 cited in this work. Following an icosahedral growth pattern, the most stable Fe12 structure adopts a penta-capped PBP with Cs symmetry, which is a distortion of C5v, in accord with all of the previous DFT studies.17,21,28,38,40,43 Our finding that the distorted icosahedral structure (C1) is the most stable atomic arrangement of Fe13 is in a morphological agreement with all of the previous DFT studies.17,21,23,28,38,40,43 However the reported symmetries are contradictory since Ma et al.28 reported D2h symmetry for their icosahedral structure, while Yuan et al.40 and Kim et al.43 reported D5d, D3d symmetries, respectively. A more symmetric icosahedron (Ci) is found to be 3.14 meV/atom (with BLYP/SDD) higher in energy.

Similar to many of the previous works, the most stable Fe14−17 clusters have been obtained as structures based on a centered hexagonal antiprism (HAP). These structures can be described as single capped centered HAP (Fe14), bicapped centered HAP or twisted, double hexagonal bipyramid (Fe15), tri-capped (Fe16), and tetra-capped centered HAP (Fe17), where the capping atoms are on the same hexagonal surface. Kim et al.43 reported C6v, D6d, C2v, and D2 symmetries for these structures, whereas we obtained C1, D3, D3h, and S4 symmetries (with BLYP/SDD), respectively, due to our tight tolerance criteria. Although Gutsev et al.38 and Yuan et al.40 suggested icosahedral structures for Fe14, they obtained the same HAP morphologies for Fe15−17 clusters. The second isomer that we obtained for Fe14 (see Fig. 3) is morphologically similar to the lowest energy structure of Gutsev et al.38 but it has 6.64 meV/atom (with BLYP/SDD) higher energy than the single capped centered HAP. Yuan et al.40 reported D3, D3h, and D2d symmetries for Fe15−17, respectively. Dieguez et al.17 reported similar HAP structures for Fe14−15 clusters with C2v (Fe14) and D6h (Fe15) symmetries. Rollmann et al.23 and Ma et al.28 obtained the same bicapped centered HAP for Fe15, while Ma et al.28 reported D6d symmetry for this cluster.

Our calculations indicate that the most stable structures of Fe18−19 clusters are based on a 13-atom truncated decahedron (TDE). The lowest energy structures of these clusters can be described as a penta-, and hexa-capped TDE with Cs, and C1 symmetries (with BLYP/SDD), which can be obtained with small distortions of D5h, and D4d symmetries, respectively. These structures were reported before by Köhler et al.,21 Gutsev et al..38 and Yuan et al.40 The symmetries given by Yuan et al.40 are C2v for Fe18 and D4d for Fe19. Ma et al.28 suggested a bi-icosahedral structure with Ih symmetry for Fe19. One should realize that the second isomer of Fe17 is a tetra-capped TDE (see Fig. 3).

Finally, our calculations suggest that the most stable Fe20 structure (C1) belongs to the HAP packing instead of icosahedral or TDE packings, which is a bicapped, double, centered HAP missing two atoms or twisted triple hexagonal bipyramid missing two atoms. Köhler et al.,21 Gutsev et al.,38 and Yuan et al.40 reported a capped bi-icosahedron for Fe20, where Yuan et al.40 identified the C2v symmetry for that structure. We identified the capped bi-icosahedron as the second isomer which has 8.8 meV/atom (with BLYP/SDD) higher energy than the lowest energy structure.

We define the total binding energy of an Fen cluster as E(n) = nEtot(Fe atom) − Etot(Fen) and the binding energy per atom as Eb(n) = E(n)/n. The second difference in energy is defined as Δ2E(n) = E(n + 1) + E(n − 1) − 2E(n). The dependence of our Eb calculations and previously reported results21,28,37,38,40,43 of Eb(n) and Δ2E(n) on the size (n) of the iron clusters is displayed in Fig. 4(a) and Fig. 4(b), respectively. Fig. 4(a) shows that while our Eb results with 6-31G** are very similar with the ones obtained by Gutsev et al.38 and Ma et al.28 especially after certain sizes, our Eb results with SDD are in good accord with the ones by Cervantes-Salguero and Seminario37 who studied Fen(n = 2 − 10) clusters using the OPBE/TZV level of DFT. The other previous Eb calculations40,43 performed by employing different GGA functionals and basis sets predicted sightly higher energies than these calculations. The DFTB results reported by Köhler et al.21 are significantly higher than all of the binding energies obtained by DFT calculations.

FIG. 4.

Dependence of the binding energy per atom and second difference in energy on the size of the iron clusters.

FIG. 4.

Dependence of the binding energy per atom and second difference in energy on the size of the iron clusters.

Close modal

Abundances in time of flight (TOF) mass spectra, called magic atomic numbers, are associated with particularly stable structures that are less likely than average to lose or acquire additional atoms in collisions. The second difference in the binding energies is the corresponding quantity to look for the magic numbers since this difference is a measure of the relative binding energy of an n-atom cluster with respect to clusters with n + 1 and n − 1 atoms. Thus, the peaks of this function indicate relatively more stable structures that can be associate with magic numbers of atoms. A TOF spectra13 revealed the sequence of magic numbers n = 7, 13, 15, 19, and 23 for iron clusters, which are indicated by dotted vertical lines in Fig. 4. One can see how sensitive are the peaks of Δ2E by comparing the plots of the data by Kim et al.43 and Yuan et al.40 Their Eb values are almost the same for the whole size range of n = 2 − 17. However, the small differences in the energies at the sizes n = 6 − 8 produce a peak in the Δ2E plot of the data by Kim et al.43 but a dip in the plot of the data by Yuan et al.40 for n = 7. Thus, as Kim et al.43 catches the experimentally observed magic number of 7, Yuan et al.40 misses it. A similar case can be found between our results with SDD and the ones by Cervantes-Salguero and Seminario.37 Although these two sets of Eb data are consistent with each others, as we observe a small peak in our Δ2E function at n = 7 (the peak is more pronounced in our calculations with 6-31G**), there is a dip in their function at the same size. The TOF peaks at the sizes n = 13, 15, and 19 are observed in many of the previous calculations21,38,40,43 along with our own results.

It is possible to further analyze the relative stabilities of iron clusters by comparing the fragmentation energies with those given by the collision induced dissociation experiments,8 which is defined as ΔE(n) = E(n) − E(n − 1) and displaced in Fig. 5 as a function of cluster size. Dissociative experimental observation of Lian et al.8 coincides well with the TOF mass spectra of Sakurai and coworkers13 since the local maxima in the dissociation experiment occur at n = 6, 7, 13, 15, and 19 similar to that of the mass spectra. In general, while Gutsev et al.’s,38 Yuan and coworkers’s,40 and our own calculations with 6-31G** overestimate the fragmentation energies, our findings with SDD underestimate them especially for the size range of n ≥ 7. All of the theoretical calculations presented in Fig. 5 agree with the experiment that the sizes of n = 6, 13, 15, and 19 are thermodynamically more stable than their neighboring sizes due to their relatively higher dissociation energies with the following exceptions: our BLYP/SDD calculations do not predict relative thermodynamic stability of Fe13 and Fe19 clusters correctly, as Yuan and coworkers40 did not predict it for Fe19. There are also some common inconsistencies between most of the theoretical works and the experiment. The only DFT study predicting the sharp dip at the size of 8 observed in the experiment is our BLYP/6-31G** results. In addition, as all of the theoretical works including our own BLYP/SDD but excluding BLYP/6-31G** methods estimate a dip at n = 9, the experiment finds a peak contrarily. Finally, our BLYP/6-31G** findings suggest an extra peak at n = 10.

FIG. 5.

Dependence of the fragmentation energy on the size of the iron clusters.

FIG. 5.

Dependence of the fragmentation energy on the size of the iron clusters.

Close modal

In Fig. 6, we compared experimental18,36,44 and theoretical21,28,38,40,43 spin magnetic moments reported in the literature along with our own calculations. We have included in the figure both the results of Stern-Gerlach experiments by Knickelbein,18 who performed magnetic deflection measurements on neutral Fen clusters without discriminating between spin and orbital contributions, and the x-ray spectroscopies by Niemeyer et al.36 and Meyer et al.,44 who presented size dependent spin and orbital magnetic moments of Fe n + cations. Except the Knickelbein’s measurements on Fe10−12 clusters, which are substantially higher than the others, the most of the magnetic moments displayed in the graph are between 2.4 and 3.6 μB where one can hardly find a consensus on the magnetic moment of an individual size. Although Knickelbein18 argued that the high moments of Fe10−12 might be due to the unquenched orbital moments, more recent findings of Niemeyer et al.36 and Meyer et al.44 are not supporting this argument since the measured orbital moments are not as high as what it was argued.

FIG. 6.

Spin magnetic moments per atom for the iron clusters.

FIG. 6.

Spin magnetic moments per atom for the iron clusters.

Close modal

A common feature of all the experimental studies is the clear dip observed in the magnetic moment of Fe13 cluster, which is related to the antiferromagnetic coupling between the central and the surface atoms of the Fe13 icosahedron. Our BLYP/SDD calculations together with the ones by Ma et al.,28 Gutsev et al.,38 and Kim et al.43 do not predict the antiferromagnetic coupling. That is why we do not include our magnetic moment results with SDD in Fig. 6. On the other hand our 6-31G** results and the ones by Köhler et al.21 and Yuan et al.40 can identify the antiferromagnetic coupling at size 13. The theoretical calculations of these three groups are close to the experimental findings of Knickelbein18 and Niemeyer et al.36 at this size, while the spin magnetic moment measurements of Meyer et al.44 are somehow lower than the others. We have found in our BYLP/6-31G** calculations of Fe13 that the central and the surface Fe atoms has -2.22 and 3.02 μB per atom spin magnetic moment, respectively, which can be compared with the values of -1.71 and 3.06 μB per atom reported by Yuan et al.40 obtained in their VASP calculations. Thus, the total spin magnetic moment of the icosahedral Fe13 cluster is calculated as 2.62 μB per atom, which is very close to the experimental prediction of 2.4 μB per atom by Knickelbein18 and Niemeyer et al.36 Another common characteristics of the experimental works is the small peak at size 14, which is predicted in our 6-31G** calculations correctly.

In the present study, we offer a systematic method called BH-DFTB/DFT for the search of the geometric, electronic, and magnetic properties of any atomic or molecular clusters and applied it for the investigation of the properties of small, free Fen(n = 2 − 20) clusters. The aim of the method is first to search the potential energy surface of the clusters as extensively and efficiently as possible, to identify their low-lying isomers and than to calculate the properties with a more accurate method. The general idea is not new. Here we are just presenting an example of this general idea and show how it works on the example of the small Fe clusters. Small Fe clusters are studied in the literature relatively largely than some other elements. Thus, the most of the structures obtained here have already been reported in the literature. However, we believe that such an extensive search of the potential energy surface of clusters should be a standard in the investigation of the properties of these clusters whenever possible in order to decrease the number of inconsistencies in the literature. The obtained results indicate that the strategy employed in the present work collects all the global minima of which different minima have been reported in the previous studies by different groups.

Small Fe clusters have three kinds of packing; icosahedral (Fe9−13), HAP (Fe14−17, Fe20), and TDE (Fe17(2), Fe18−19). The energy differences between the low-lying isomers of these structures are sometimes so small that different xc-functionals and basis sets employed in the DFT calculations may produce different order of isomers, which can be considered as another reason of the present inconsistencies. Many of the theoretical works including the present study are in a qualitative agreement with the TOF mass spectra that the magic numbers for the small Fe clusters are 7, 13, 15, and 19 and with the collision induced dissociation experiments that the sizes 6, 7, 13, 15, and 19 are thermodynamically more stable than their neighboring sizes. There is no quantitative agreement between neither the experimental nor the theoretical calculations of the magnetic moment calculations. The spin magnetic moment per atom is between 2.4 and 3.6 μB for the most of the sizes considered here. The antiferromagnetic coupling between the central and the surface atoms of the Fe13 icosahedron, which is previously reported by some experimental and theoretical studies, is verified by our calculations in this work as well. The quantitative disagreements between the calculated and measured values of the magnetic moments of the individual sizes are still to be resolved.

This research is supported by the Scientific and Technical Research Council of Turkey (TUBITAK) grant number 114F297. The numerical calculations reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).

1.
P.
Jena
and
J.A.W.
Cattleman
,
Nanoclusters: A Bridge Across Disciplines
(
Elsevier
,
Oxford
,
2010
).
2.
J.
Zhao
,
X.
Huang
,
P.
Jin
, and
Z.
Chen
,
Coord. Chem. Rev.
289-290
,
315
(
2015
).
3.
P.A.
Montano
and
G.K.
Shenoy
,
Solid. State Commun.
35
,
53
(
1980
).
4.
M.
Moskovits
and
D.P.
DiLella
,
J. Chem. Phys.
73
,
4917
(
1980
).
5.
H.
Purdum
,
P.A.
Montano
,
G.K.
Shenoy
, and
T.
Morrison
,
Phys. Rev. B
25
,
4412
(
1982
).
6.
D.M.
Cox
,
D.J.
Trevor
,
R.L.
Whetten
,
E.A.
Rohlfing
, and
A.
Kaldor
,
Phys. Rev. B
32
,
7290
(
1985
).
7.
D.G.
Leopold
and
W.C.
Lineberger
,
J. Chem. Phys.
85
,
51
(
1986
).
8.
L.
Lian
,
C.X.
Su
, and
P.B.
Armentrout
,
J. Chem. Phys.
97
,
4072
(
1992
).
9.
I.M.L.
Billas
,
J.A.
Becker
,
A.
Chatelain
, and
W.A.
de Heer
,
Phys. Rev. Lett.
71
,
4067
(
1993
).
10.
I.M.L.
Billas
,
A.
Chatelain
, and
W.A.
de Heer
,
Science
265
,
1682
(
1994
).
11.
X.G.
Gong
and
Q.Q.
Zheng
,
J. Phys.: Condens. Matter
7
,
2421
(
1995
).
12.
I.M.L.
Billas
,
A.
Chatelain
, and
W.A.
de Heer
,
J. Magn. Magn. Mater.
168
,
64
(
1997
).
13.
M.
Sakurai
,
K.
Watanabe
,
K.
Sumiyama
, and
K.
Suzuki
,
J. Chem. Phys.
111
,
235
(
1999
).
14.
T.L.
Haslett
,
K.A.
Bosnick
,
S.
Fedrigo
, and
M.
Moskovits
,
J. Chem. Phys.
111
,
6456
(
1999
).
15.
A.
Hirt
,
D.
Gerion
,
I.M.L.
Billas
,
A.
Chatelain
, and
W.A.
de Heer
,
Phys. Rev. B
62
,
7491
(
2000
).
16.
P.B.
Armentrout
,
Annu. Rev. Phys. Chem.
52
,
423
(
2001
).
17.
O.
Dieguez
,
M.M.G.
Alemany
,
C.
Rey
,
P.
Ordejon
, and
L.J.
Gallego
,
Phys. Rev. B
63
,
205407
(
2001
).
18.
M.B.
Knickelbein
,
Chem. Phys. Lett.
353
,
221
(
2002
).
19.
P.
Bobadova-Parvanova
,
K.A.
Jackson
,
S.
Srinivas
,
M.
Horoi
,
C.
Köhler
, and
G.
Seifert
,
J. Chem. Phys.
116
,
3576
(
2002
).
20.
G.L.
Gutsev
and
C.W.
Bauschlicher
, Jr.
,
J. Phys. Chem. A
107
,
7013
(
2003
).
21.
C.
Köhler
,
G.
Seifert
, and
T.
Frauenheim
,
Chem. Phys.
309
,
23
(
2005
).
22.
S.
Li
,
M.M.G.
Alemany
, and
J.R.
Chelikowsky
,
Phys. Rev. B
73
,
233404
(
2006
).
23.
G.
Rollmann
,
P.
Entel
, and
S.
Sahoo
,
Comput. Mater. Sci.
35
,
275
(
2006
).
24.
M.L.
Tiago
,
Y.
Zhou
,
M.M.G.
Alemany
,
Y.
Saad
, and
J.R.
Chelikowsky
,
Phys. Rev. Lett.
97
,
147201
(
2006
).
25.
I.A.
Gass
,
C.J.
Milios
,
M.
Evangelisti
,
S.L.
Heath
,
D.
Collision
,
S.
Parsons
, and
E.K.
Brechin
,
Polyhedron
26
,
1835
(
2007
).
26.
C.
Köhler
and
T.
Frauenheim
,
J. Comput. Theor. Nanos.
4
,
264
(
2007
).
27.
S.
Yu
,
S.
Chen
,
W.W.
Zhang
,
L.
Yu
, and
Y.
Yin
,
Chem. Phys. Lett.
446
,
217
(
2007
).
28.
Q.M.
Ma
,
Z.
Xie
,
J.
Wang
,
Y.
Liu
, and
Y.C.
Li
,
Solid State Commun.
142
,
114
(
2007
).
29.
R.
Singh
and
P.
Kroll
,
Phys. Rev. B
78
,
245404
(
2008
).
30.
K.
Boufala
,
L.
Fernandez-Seivane
,
J.
Ferrer
, and
M.
Samah
,
J. Magn. Magn. Mater.
322
,
3428
(
2010
).
31.
D.R.
Roy
,
R.
Robles
, and
S.N.
Khan
,
J. Chem. Phys.
132
,
194305
(
2010
).
32.
C.
Angeli
and
R.
Cimiraglia
,
Mol. Phys.
109
,
1503
(
2011
).
33.
S.
Datta
,
M.
Kabir
, and
T.
Saha-Dasgupta
,
Phys. Rev. B
84
,
075429
(
2011
).
34.
X.
Xu
,
S.
Yin
,
R.
Moro
,
A.
Liang
,
J.
Bowlan
, and
W.A.
de Heer
,
Phys. Rev. Lett.
107
,
057203
(
2011
).
35.
Y.
Okamoto
,
F.
Kawamura
,
Y.
Ohta
,
A.J.
Page
,
S.
Irle
, and
K.
Morokuma
,
J. Comput. Theor. Nanos.
8
,
1755
(
2011
).
36.
M.
Niemeyer
,
K.
Hirsch
,
V.
Zamudio-Bayer
,
A.
Langenberg
,
M.
Vogel
,
M.
Kossick
,
C.
Ebrecht
,
K.
Egashira
,
A.
Terasaki
,
T.
Möller
,
B. v.
Issendorff
, and
J.T.
Lau
,
Phys. Rev. Lett.
108
,
057201
(
2012
).
37.
K.
Cervantes-Salguero
and
J.M.
Seminario
,
J. Mol. Model.
18
,
4043
(
2012
).
38.
G.L.
Gutsev
,
C.A.
Weatherford
,
P.
Jena
,
E.
Johnson
, and
B.R.
Ramachandran
,
J. Phys. Chem. A
116
,
10218
(
2012
).
39.
G.L.
Gutsev
,
C.A.
Weatherford
,
K.G.
Belay
,
B.R.
Ramachandran
, and
P.
Jena
,
J. Chem. Phys.
138
,
164303
(
2013
).
40.
H.K.
Yuan
,
H.
Chen
,
A.L.
Kuang
,
C.L.
Tian
, and
J.Z.
Wang
,
J. Chem. Phys.
139
,
034314
(
2013
).
41.
W.
Song
,
M.
Jiao
,
K.
Li
,
Y.
Wang
, and
Z.
Wu
,
Chem. Phys. Lett.
588
,
203
(
2013
).
42.
A.
Langenberg
,
K.
Hirsch
,
A.
Lawicki
,
V.
Zamudio-Bayer
,
M.
Niemeyer
,
P.
Chmiela
,
B.
Langbehn
,
A.
Terasaki
,
B. v.
Issendorff
, and
J.T.
Lau
,
Phys. Rev. B
90
,
184420
(
2014
).
43.
E.
Kim
,
A.
Mohrland
,
P.F.
Weck
,
T.
Pang
,
K.R.
Czerwinski
, and
D.
Tomanek
,
Chem. Phys. Lett.
613
,
59
(
2014
).
44.
J.
Meyer
,
M.
Tombers
,
C.
van Wüllen
,
G.
Niedner-Schatteburg
,
S.
Peredkov
,
W.
Eberhardt
,
M.
Neeb
,
S.
Palutke
,
M.
Martins
, and
W.
Wurth
,
J. Chem. Phys.
143
,
104302
(
2015
).
45.
C.T.
Chen
,
Y.U.
Idzerda
,
H.-J.
Lin
,
N.V.
Smith
,
G.
Meigs
,
E.
Chaban
,
G.H.
Ho
,
E.
Pellegrin
, and
F.
Sette
,
Phys. Rev. Lett.
75
,
152
(
1995
).
46.
B.
Aradi
,
B.
Hourahine
, and
Th.
Frauenheim
,
J. Phys. Chem. A
111
,
5678
(
2007
).
47.
D.J.
Wales
and
J.P.K.
Doye
,
J. Phys. Chem. A
101
,
5111
(
1997
).
48.
J.P.K.
Doye
,
D.J.
Wales
, and
R.S.
Berry
,
J. Chem. Phys.
103
,
4234
(
1995
).
49.
J.A.
Elliot
,
Y.
Shibuta
, and
D.J.
Wales
,
Philos. Mag.
89
,
3311
(
2009
).
50.
G.
Zheng
,
H.A.
Witek
,
P.
Bobadova-Parvanova
,
S.
Irle
,
D.G.
Musaev
,
R.
Prabhakar
,
K.
Morokuma
,
M.
Lundberg
,
M.
Elstner
,
C.
Köhler
, and
Th.
Frauenheim
,
J. Chem. Theory Comput.
3
,
1349
(
2007
).
51.
M.
Valiev
,
E.J.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T.P.
Straatsma
,
H.J.J.
van Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T.L.
Windus
, and
W.A.
de Jong
,
Comput. Phys. Commun.
181
,
1477
(
2010
).
52.
M.
Dolg
,
U.
Wedig
,
H.
Stoll
, and
H.
Preuss
,
J. Chem. Phys.
86
,
866
(
1987
).
53.
V.
Rassolov
,
J.A.
Pople
,
M.
Ratner
, and
T.L.
Windus
,
J. Chem. Phys.
109
,
1223
(
1998
).
54.
A.D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
);
A.D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
55.
C.
Lee
,
W.
Yang
, and
R.G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
56.
P.J.
Hay
and
W.R.
Wadt
,
J. Chem. Phys.
82
,
270
(
1985
);
P.J.
Hay
and
W.R.
Wadt
,
J. Chem. Phys.
82
,
284
(
1985
);
P.J.
Hay
and
W.R.
Wadt
,
J. Chem. Phys.
82
,
299
(
1985
).
57.
M.M.
Hurley
,
L.F.
Pacios
,
P.A.
Christiansen
,
R.B.
Ross
, and
W.C.
Ermler
,
J. Chem. Phys.
84
,
6840
(
1986
).
58.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
59.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
108
,
664
(
1998
).
60.
Y.
Zhao
and
D.G.
Truhlar
,
J. Phys. Chem. A
109
,
5656
(
2005
).
61.
J.P.
Perdew
,
J.A.
Chevary
,
S.H.
Vosko
,
K.A.
Jackson
,
M.R.
Pederson
,
D.J.
Singh
, and
C.
Fiolhais
,
Phys. Rev. B
46
,
6671
(
1992
).
62.
J.
Tao
,
J.
Perdew
,
V.
Staroverov
, and
G.
Scuseria
,
Phys. Rev. Let.
91
,
146401
(
2003
).
63.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
64.
F.
Furche
and
J.P.
Perdew
,
J. Chem. Phys.
124
,
044103
(
2006
).