We present a study on the structural, electronic, and magnetic properties of Fe_{n}(*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 (Fe_{9−13}), centered hexagonal antiprism (Fe_{14−17}, Fe_{20}), and truncated decahedral (Fe_{17(2)}, Fe_{18−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 Fe_{n}(*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 Fe

_{13}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.

## I. INTRODUCTION

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 Fe_{n}(*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 Fe_{13} and Fe$ 13 + $ clusters.^{8} The Fe-Fe bond length was measured for Fe_{2} 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 Fe_{2},^{4} and Fe$ 2 \u2212 $,^{7} and for matrix-isolated Fe_{3}.^{14} The geometry of a neutral Fe_{17} 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 Fe

_{n}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 Fe_{n} 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 Fe_{2} 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 Fe_{n}(*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 Fe_{n}(*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.

## II. COMPUTATIONAL METHOD

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 results^{28,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.

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.

n
. | DFTB . | DFT . | n
. | DFTB . | DFT . | n
. | DFTB . | DFT . | n
. | DFTB . | DFT . |
---|---|---|---|---|---|---|---|---|---|---|---|

5 | 1 | 2 | 9 | 1 | 2 | 13 | 1 | 2 | 17 | 1 | 4 |

5 | 2 | 8 | 9 | 2 | 1 | 13 | 2 | 9 | 17 | 2 | 6 |

5 | 3 | 5 | 9 | 3 | 10 | 13 | 3 | 6 | 17 | 3 | 9 |

5 | 4 | 9 | 9 | 4 | 3 | 13 | 4 | 7 | 17 | 4 | 5 |

5 | 5 | 1 | 9 | 5 | 9 | 13 | 5 | 1 | 17 | 5 | 10 |

5 | 6 | 10 | 9 | 6 | 8 | 13 | 6 | 5 | 17 | 6 | 7 |

5 | 7 | 7 | 9 | 7 | 4 | 13 | 7 | 3 | 17 | 7 | 8 |

5 | 8 | 6 | 9 | 8 | 6 | 13 | 8 | 8 | 17 | 8 | 2 |

5 | 9 | 4 | 9 | 9 | 5 | 13 | 9 | 10 | 17 | 9 | 3 |

5 | 10 | 3 | 9 | 10 | 7 | 13 | 10 | 4 | 17 | 10 | 1 |

6 | 1 | 1 | 10 | 1 | 1 | 14 | 1 | 6 | 18 | 1 | 1 |

6 | 2 | 2 | 10 | 2 | 5 | 14 | 2 | 2 | 18 | 2 | 10 |

6 | 3 | 7 | 10 | 3 | 6 | 14 | 3 | 1 | 18 | 3 | 7 |

6 | 4 | 6 | 10 | 4 | 2 | 14 | 4 | 9 | 18 | 4 | 3 |

6 | 5 | 5 | 10 | 5 | 3 | 14 | 5 | 10 | 18 | 5 | 4 |

6 | 6 | 3 | 10 | 6 | 4 | 14 | 6 | 8 | 18 | 6 | 2 |

6 | 7 | 8 | 10 | 7 | 7 | 14 | 7 | 5 | 18 | 7 | 6 |

6 | 8 | 9 | 10 | 8 | 8 | 14 | 8 | 4 | 18 | 8 | 9 |

6 | 9 | 4 | 10 | 9 | 10 | 14 | 9 | 7 | 18 | 9 | 5 |

6 | 10 | 10 | 9 | 14 | 10 | 3 | 18 | 10 | 8 | ||

7 | 1 | 3 | 11 | 1 | 10 | 15 | 1 | 1 | 19 | 1 | 1 |

7 | 2 | 4 | 11 | 2 | 1 | 15 | 2 | 2 | 19 | 2 | 4 |

7 | 3 | 9 | 11 | 3 | 2 | 15 | 3 | 3 | 19 | 3 | 8 |

7 | 4 | 10 | 11 | 4 | 3 | 15 | 4 | 4 | 19 | 4 | 9 |

7 | 5 | 1 | 11 | 5 | 5 | 15 | 5 | 5 | 19 | 5 | 5 |

7 | 6 | 2 | 11 | 6 | 8 | 15 | 6 | 6 | 19 | 6 | 3 |

7 | 7 | 8 | 11 | 7 | 4 | 15 | 7 | 10 | 19 | 7 | 2 |

7 | 8 | 5 | 11 | 8 | 7 | 15 | 8 | 9 | 19 | 8 | 6 |

7 | 9 | 6 | 11 | 9 | 9 | 15 | 9 | 8 | 19 | 9 | 7 |

7 | 10 | 7 | 11 | 10 | 6 | 15 | 10 | 7 | 19 | 10 | 10 |

8 | 1 | 10 | 12 | 1 | 4 | 16 | 1 | 2 | 20 | 1 | 5 |

8 | 2 | 5 | 12 | 2 | 1 | 16 | 2 | 4 | 20 | 2 | 2 |

8 | 3 | 9 | 12 | 3 | 3 | 16 | 3 | 3 | 20 | 3 | 4 |

8 | 4 | 1 | 12 | 4 | 10 | 16 | 4 | 6 | 20 | 4 | 1 |

8 | 5 | 2 | 12 | 5 | 9 | 16 | 5 | 5 | 20 | 5 | 10 |

8 | 6 | 3 | 12 | 6 | 8 | 16 | 6 | 7 | 20 | 6 | 9 |

8 | 7 | 6 | 12 | 7 | 7 | 16 | 7 | 8 | 20 | 7 | 6 |

8 | 8 | 8 | 12 | 8 | 5 | 16 | 8 | 9 | 20 | 8 | 7 |

8 | 9 | 7 | 12 | 9 | 2 | 16 | 9 | 10 | 20 | 9 | 8 |

8 | 10 | 4 | 12 | 10 | 6 | 16 | 10 | 1 | 20 | 10 | 3 |

n
. | DFTB . | DFT . | n
. | DFTB . | DFT . | n
. | DFTB . | DFT . | n
. | DFTB . | DFT . |
---|---|---|---|---|---|---|---|---|---|---|---|

5 | 1 | 2 | 9 | 1 | 2 | 13 | 1 | 2 | 17 | 1 | 4 |

5 | 2 | 8 | 9 | 2 | 1 | 13 | 2 | 9 | 17 | 2 | 6 |

5 | 3 | 5 | 9 | 3 | 10 | 13 | 3 | 6 | 17 | 3 | 9 |

5 | 4 | 9 | 9 | 4 | 3 | 13 | 4 | 7 | 17 | 4 | 5 |

5 | 5 | 1 | 9 | 5 | 9 | 13 | 5 | 1 | 17 | 5 | 10 |

5 | 6 | 10 | 9 | 6 | 8 | 13 | 6 | 5 | 17 | 6 | 7 |

5 | 7 | 7 | 9 | 7 | 4 | 13 | 7 | 3 | 17 | 7 | 8 |

5 | 8 | 6 | 9 | 8 | 6 | 13 | 8 | 8 | 17 | 8 | 2 |

5 | 9 | 4 | 9 | 9 | 5 | 13 | 9 | 10 | 17 | 9 | 3 |

5 | 10 | 3 | 9 | 10 | 7 | 13 | 10 | 4 | 17 | 10 | 1 |

6 | 1 | 1 | 10 | 1 | 1 | 14 | 1 | 6 | 18 | 1 | 1 |

6 | 2 | 2 | 10 | 2 | 5 | 14 | 2 | 2 | 18 | 2 | 10 |

6 | 3 | 7 | 10 | 3 | 6 | 14 | 3 | 1 | 18 | 3 | 7 |

6 | 4 | 6 | 10 | 4 | 2 | 14 | 4 | 9 | 18 | 4 | 3 |

6 | 5 | 5 | 10 | 5 | 3 | 14 | 5 | 10 | 18 | 5 | 4 |

6 | 6 | 3 | 10 | 6 | 4 | 14 | 6 | 8 | 18 | 6 | 2 |

6 | 7 | 8 | 10 | 7 | 7 | 14 | 7 | 5 | 18 | 7 | 6 |

6 | 8 | 9 | 10 | 8 | 8 | 14 | 8 | 4 | 18 | 8 | 9 |

6 | 9 | 4 | 10 | 9 | 10 | 14 | 9 | 7 | 18 | 9 | 5 |

6 | 10 | 10 | 9 | 14 | 10 | 3 | 18 | 10 | 8 | ||

7 | 1 | 3 | 11 | 1 | 10 | 15 | 1 | 1 | 19 | 1 | 1 |

7 | 2 | 4 | 11 | 2 | 1 | 15 | 2 | 2 | 19 | 2 | 4 |

7 | 3 | 9 | 11 | 3 | 2 | 15 | 3 | 3 | 19 | 3 | 8 |

7 | 4 | 10 | 11 | 4 | 3 | 15 | 4 | 4 | 19 | 4 | 9 |

7 | 5 | 1 | 11 | 5 | 5 | 15 | 5 | 5 | 19 | 5 | 5 |

7 | 6 | 2 | 11 | 6 | 8 | 15 | 6 | 6 | 19 | 6 | 3 |

7 | 7 | 8 | 11 | 7 | 4 | 15 | 7 | 10 | 19 | 7 | 2 |

7 | 8 | 5 | 11 | 8 | 7 | 15 | 8 | 9 | 19 | 8 | 6 |

7 | 9 | 6 | 11 | 9 | 9 | 15 | 9 | 8 | 19 | 9 | 7 |

7 | 10 | 7 | 11 | 10 | 6 | 15 | 10 | 7 | 19 | 10 | 10 |

8 | 1 | 10 | 12 | 1 | 4 | 16 | 1 | 2 | 20 | 1 | 5 |

8 | 2 | 5 | 12 | 2 | 1 | 16 | 2 | 4 | 20 | 2 | 2 |

8 | 3 | 9 | 12 | 3 | 3 | 16 | 3 | 3 | 20 | 3 | 4 |

8 | 4 | 1 | 12 | 4 | 10 | 16 | 4 | 6 | 20 | 4 | 1 |

8 | 5 | 2 | 12 | 5 | 9 | 16 | 5 | 5 | 20 | 5 | 10 |

8 | 6 | 3 | 12 | 6 | 8 | 16 | 6 | 7 | 20 | 6 | 9 |

8 | 7 | 6 | 12 | 7 | 7 | 16 | 7 | 8 | 20 | 7 | 6 |

8 | 8 | 8 | 12 | 8 | 5 | 16 | 8 | 9 | 20 | 8 | 7 |

8 | 9 | 7 | 12 | 9 | 2 | 16 | 9 | 10 | 20 | 9 | 8 |

8 | 10 | 4 | 12 | 10 | 6 | 16 | 10 | 1 | 20 | 10 | 3 |

NWChem 6.1 program package^{51} 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 1997^{52} basis sets (SDD) and effective core potentials (ECP), where the outer most 16 electrons of free Fe atom (3*s*^{2}3*p*^{6}3*d*^{6}4*s*^{2}) are treated explicitly, and the corresponding contracted Gaussian basis functions are (6*s*5*p*3*d*1*f*). The second one is 6-31G** all electron Gaussian basis set^{53} with contracted basis functions of 5*s*4*p*2*d*1*f*. Becke’s GGA exchange functional^{54} and Lee-Yang-Parr correlation functional^{55} (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 Fe_{2} 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/*a*_{0} 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.

## III. RESULTS AND DISCUSSION

### A. The choice of xc-functional and basis set

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 Fe_{2} 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} QZVP^{58}) plus the 6-31G** all electron basis set^{53} and four xc functionals (two GGA: BLYP,^{54,55} MPWPW91,^{59–61} a meta-GGA: TPSS,^{62} and a hybrid: PBE0^{63}) have been tried.

. | Magnetic Moment (μ)
. _{B} | Binding Energy (eV) . | ||||||
---|---|---|---|---|---|---|---|---|

. | Experimental^{6} = 6.5 ± 1
. | Experimental^{8} = 1.14 ± 0.10
. | ||||||

Basis Sets . | BLYP . | MPWPW91 . | TPSS . | PBE0 . | BLYP . | MPWPW91 . | TPSS . | PBE0 . |

LANL2DZ ECP (3s3p2d) | 6 | 6 | 6 | 8 | 1.87 | 1.91 | 1.55 | 0.31 |

LANL2TZ ECP (5s5p3d) | 6 | 6 | 8 | 8 | 1.51 | 1.71 | 1.63 | -0.09 |

CRENBL ECP (7s6p6d) | 6 | 6 | 6 | 8 | 1.44 | 1.60 | 1.47 | 1.10 |

SDD ECP (6s5p3d1f) | 6 | 6 | 8 | 8 | 1.51 | 1.69 | 1.58 | 0.91 |

QZVP ECP (11s6p5d3f1g ) | 6 | 6 | 6 | 8 | 2.27 | 2.33 | 1.92 | 0.98 |

6-31G** All Elec. (5s4p2d1f) | 8 | 8 | 8 | 8 | 1.69 | 2.10 | 1.75 | -0.31 |

Bond Length (Å) | ||||||||

Experimental^{3} = 1.87 ± 0.13 | Vibrational freq. (cm^{−1}) | |||||||

Experimental^{5} = 2.02 ± 0.02 | Experimental^{4} = 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) . | ||||||
---|---|---|---|---|---|---|---|---|

. | Experimental^{6} = 6.5 ± 1
. | Experimental^{8} = 1.14 ± 0.10
. | ||||||

Basis Sets . | BLYP . | MPWPW91 . | TPSS . | PBE0 . | BLYP . | MPWPW91 . | TPSS . | PBE0 . |

LANL2DZ ECP (3s3p2d) | 6 | 6 | 6 | 8 | 1.87 | 1.91 | 1.55 | 0.31 |

LANL2TZ ECP (5s5p3d) | 6 | 6 | 8 | 8 | 1.51 | 1.71 | 1.63 | -0.09 |

CRENBL ECP (7s6p6d) | 6 | 6 | 6 | 8 | 1.44 | 1.60 | 1.47 | 1.10 |

SDD ECP (6s5p3d1f) | 6 | 6 | 8 | 8 | 1.51 | 1.69 | 1.58 | 0.91 |

QZVP ECP (11s6p5d3f1g ) | 6 | 6 | 6 | 8 | 2.27 | 2.33 | 1.92 | 0.98 |

6-31G** All Elec. (5s4p2d1f) | 8 | 8 | 8 | 8 | 1.69 | 2.10 | 1.75 | -0.31 |

Bond Length (Å) | ||||||||

Experimental^{3} = 1.87 ± 0.13 | Vibrational freq. (cm^{−1}) | |||||||

Experimental^{5} = 2.02 ± 0.02 | Experimental^{4} = 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 \Sigma g \u2212 $ 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

*μ*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.

_{B}In accordance with the previous theoretical calculations^{64} 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.

### B. Geometrical Configurations

We present the optimized geometries of the lowest energy structures and a few low-lying isomers of Fe_{n}(*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 (*E _{b}*), total spin magnetic moments (

*μ*), average bond distances (

_{s}*d*), HOMO-LUMO gaps for the spin down electrons (

_{a}*E*), and the lowest (

_{gap}*w*) and the highest (

_{l}*w*) 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 Fe

_{h}_{3}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.

n
. | Symmetry . | E
. _{b} | μ
. _{s} | d
. _{a} | E
. _{gap} | w- _{l}w
. _{h} |
---|---|---|---|---|---|---|

2 | D_{∞h} | 0.76 / 0.43 | 3.00 / 4.00 | 2.02 / 2.03 | 0.62 / 1.02 | 393.0 / 416.1 |

3 | C / _{s}C_{2v} | 1.12 / 1.23 | 3.33 | 2.28 / 2.26 | 0.78 / 0.26 | 9.2 / 155.5 - 335.5 / 381.3 |

4 | C_{3v} / C_{3} | 1.43 / 1.75 | 3.50 | 2.38 / 2.29 | 0.71 / 0.40 | 46.5 / 21.5 - 344.4 / 406.8 |

5 | C_{2v} / C _{s} | 1.70 / 2.15 | 3.20 | 2.39 / 2.30 | 0.59 / 0.26 | 52.6 / 46.7 - 341.9 / 406.3 |

6 | D_{4h} | 1.93 / 2.44 | 3.33 | 2.46 / 2.37 | 0.48 / 0.51 | 104.3 / 72.3 - 316.0 / 394.1 |

7 | C_{1} / C _{s} | 2.05 / 2.62 | 3.14 | 2.44 / 2.35 | 0.73 / 0.35 | 73.2 / 102.2 - 315.6 / 386.5 |

8 | C_{2v} | 2.10 / 2.70 | 3.00 | 2.45 / 2.35 | 0.46 / 0.29 | 104.0 / 123.4 - 318.4 / 368.9 |

9 | C_{2v} | 2.13 / 2.79 | 2.89 | 2.44 / 2.35 | 0.55 / 0.37 | 86.7.1 / 76.7 - 324.9 / 381.5 |

10 | C_{1} / C_{3} | 2.17 / 2.87 | 3.00 | 2.46 / 2.37 | 0.45 / 0.40 | -70.4 / -7.4 - 319.5 / 252.7 |

11 | C_{1} / C_{2v} | 2.21 / 2.92 | 3.09 | 2.48 / 2.39 | 0.53 / 0.26 | 45.0^{a} / 30.9 - 345.5^{a} / 399.0 |

12 | C _{s} | 2.25 / 2.97 | 3.17 | 2.50 / 2.41 | 0.38 / 0.28 | 40.0 / 20.1 - 326.0 / 395.2 |

13 | C_{1} | 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 | C_{1} / C_{6} | 2.34 / 3.09 | 3.29 | 2.52 / 2.44 | 0.40 / 0.33 | 59.9 / 62.8 - 279.4 / 328.1 |

15 | D_{3} | 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 | D_{3h} | 2.43 / 3.22 | 3.25 | 2.54 / 2.45 | 0.37 / 0.34 | 24.1 / 65.3 - 286.0 / 328.2 |

17 | S_{4} / C _{s} | 2.46 / 3.23 | 3.18 | 2.54 / 2.46 | 0.28 / 0.23 | 26.7^{a} / -48.7 - 276.8^{a} / 328.7 |

18 | C / _{s}C_{2v} | 2.47 / 3.27 | 3.11 | 2.53 / 2.47 | 0.23 / 0.25 | 68.5 / 82.4 - 276.4 / 319.7 |

19 | C_{1} | 2.48 / 3.30 | 3.05 | 2.54 / 2.45 | 0.44 / 0.50 | 54.6 / 80.1 - 284.5 / 327.0 |

20 | C_{1} | 2.48 / 3.22 | 3.00 | 2.53 / 2.44 | 0.21 / 0.19 | 69.1^{a} / 42.6 - 367.8^{a} / 386.7 |

n
. | Symmetry . | E
. _{b} | μ
. _{s} | d
. _{a} | E
. _{gap} | w- _{l}w
. _{h} |
---|---|---|---|---|---|---|

2 | D_{∞h} | 0.76 / 0.43 | 3.00 / 4.00 | 2.02 / 2.03 | 0.62 / 1.02 | 393.0 / 416.1 |

3 | C / _{s}C_{2v} | 1.12 / 1.23 | 3.33 | 2.28 / 2.26 | 0.78 / 0.26 | 9.2 / 155.5 - 335.5 / 381.3 |

4 | C_{3v} / C_{3} | 1.43 / 1.75 | 3.50 | 2.38 / 2.29 | 0.71 / 0.40 | 46.5 / 21.5 - 344.4 / 406.8 |

5 | C_{2v} / C _{s} | 1.70 / 2.15 | 3.20 | 2.39 / 2.30 | 0.59 / 0.26 | 52.6 / 46.7 - 341.9 / 406.3 |

6 | D_{4h} | 1.93 / 2.44 | 3.33 | 2.46 / 2.37 | 0.48 / 0.51 | 104.3 / 72.3 - 316.0 / 394.1 |

7 | C_{1} / C _{s} | 2.05 / 2.62 | 3.14 | 2.44 / 2.35 | 0.73 / 0.35 | 73.2 / 102.2 - 315.6 / 386.5 |

8 | C_{2v} | 2.10 / 2.70 | 3.00 | 2.45 / 2.35 | 0.46 / 0.29 | 104.0 / 123.4 - 318.4 / 368.9 |

9 | C_{2v} | 2.13 / 2.79 | 2.89 | 2.44 / 2.35 | 0.55 / 0.37 | 86.7.1 / 76.7 - 324.9 / 381.5 |

10 | C_{1} / C_{3} | 2.17 / 2.87 | 3.00 | 2.46 / 2.37 | 0.45 / 0.40 | -70.4 / -7.4 - 319.5 / 252.7 |

11 | C_{1} / C_{2v} | 2.21 / 2.92 | 3.09 | 2.48 / 2.39 | 0.53 / 0.26 | 45.0^{a} / 30.9 - 345.5^{a} / 399.0 |

12 | C _{s} | 2.25 / 2.97 | 3.17 | 2.50 / 2.41 | 0.38 / 0.28 | 40.0 / 20.1 - 326.0 / 395.2 |

13 | C_{1} | 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 | C_{1} / C_{6} | 2.34 / 3.09 | 3.29 | 2.52 / 2.44 | 0.40 / 0.33 | 59.9 / 62.8 - 279.4 / 328.1 |

15 | D_{3} | 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 | D_{3h} | 2.43 / 3.22 | 3.25 | 2.54 / 2.45 | 0.37 / 0.34 | 24.1 / 65.3 - 286.0 / 328.2 |

17 | S_{4} / C _{s} | 2.46 / 3.23 | 3.18 | 2.54 / 2.46 | 0.28 / 0.23 | 26.7^{a} / -48.7 - 276.8^{a} / 328.7 |

18 | C / _{s}C_{2v} | 2.47 / 3.27 | 3.11 | 2.53 / 2.47 | 0.23 / 0.25 | 68.5 / 82.4 - 276.4 / 319.7 |

19 | C_{1} | 2.48 / 3.30 | 3.05 | 2.54 / 2.45 | 0.44 / 0.50 | 54.6 / 80.1 - 284.5 / 327.0 |

20 | C_{1} | 2.48 / 3.22 | 3.00 | 2.53 / 2.44 | 0.21 / 0.19 | 69.1^{a} / 42.6 - 367.8^{a} / 386.7 |

^{a}

PBE0/SDD.

There is an agreement in the literature^{17,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 Fe_{9} is a bicapped PBP in *C*_{2v} 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 (*D*_{3h}) 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 Fe_{10} 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 Fe_{11} structure is a tetra-capped PBP with *C*_{1} symmetry is consistent with all of the previously reported DFT calculations^{17,21,28,38,40,43} cited in this work. Following an icosahedral growth pattern, the most stable Fe_{12} structure adopts a penta-capped PBP with *C _{s}* symmetry, which is a distortion of

*C*

_{5v}, in accord with all of the previous DFT studies.

^{17,21,28,38,40,43}Our finding that the distorted icosahedral structure (

*C*

_{1}) is the most stable atomic arrangement of Fe

_{13}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

*D*

_{2h}symmetry for their icosahedral structure, while Yuan

*et al.*

^{40}and Kim

*et al.*

^{43}reported

*D*

_{5d},

*D*

_{3d}symmetries, respectively. A more symmetric icosahedron (

*C*) is found to be 3.14 meV/atom (with BLYP/SDD) higher in energy.

_{i}Similar to many of the previous works, the most stable Fe_{14−17} clusters have been obtained as structures based on a centered hexagonal antiprism (HAP). These structures can be described as single capped centered HAP (Fe_{14}), bicapped centered HAP or twisted, double hexagonal bipyramid (Fe_{15}), tri-capped (Fe_{16}), and tetra-capped centered HAP (Fe_{17}), where the capping atoms are on the same hexagonal surface. Kim *et al.*^{43} reported *C*_{6v}, *D*_{6d}, *C*_{2v}, and *D*_{2} symmetries for these structures, whereas we obtained *C*_{1}, *D*_{3}, *D*_{3h}, and *S*_{4} 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 Fe_{14}, they obtained the same HAP morphologies for Fe_{15−17} clusters. The second isomer that we obtained for Fe_{14} (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 *D*_{3}, *D*_{3h}, and *D*_{2d} symmetries for Fe_{15−17}, respectively. Dieguez *et al.*^{17} reported similar HAP structures for Fe_{14−15} clusters with *C*_{2v} (Fe_{14}) and *D*_{6h} (Fe_{15}) symmetries. Rollmann *et al.*^{23} and Ma *et al.*^{28} obtained the same bicapped centered HAP for Fe_{15}, while Ma *et al.*^{28} reported *D*_{6d} symmetry for this cluster.

Our calculations indicate that the most stable structures of Fe_{18−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 *C _{s}*, and

*C*

_{1}symmetries (with BLYP/SDD), which can be obtained with small distortions of

*D*

_{5h}, and

*D*

_{4d}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

*C*

_{2v}for Fe

_{18}and

*D*

_{4d}for Fe

_{19}. Ma

*et al.*

^{28}suggested a bi-icosahedral structure with

*I*symmetry for Fe

_{h}_{19}. One should realize that the second isomer of Fe

_{17}is a tetra-capped TDE (see Fig. 3).

Finally, our calculations suggest that the most stable Fe_{20} structure (*C*_{1}) 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 Fe_{20}, where Yuan *et al.*^{40} identified the *C*_{2v} 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.

### C. Binding Energies

We define the total binding energy of an Fe_{n} cluster as *E*(*n*) = *nE _{tot}*(Fe atom) −

*E*(Fe

_{tot}_{n}) and the binding energy per atom as

*E*(

_{b}*n*) =

*E*(

*n*)/

*n*. The second difference in energy is defined as Δ

_{2}

*E*(

*n*) =

*E*(

*n*+ 1) +

*E*(

*n*− 1) − 2

*E*(

*n*). The dependence of our

*E*calculations and previously reported results

_{b}^{21,28,37,38,40,43}of

*E*(

_{b}*n*) and Δ

_{2}

*E*(

*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

*E*results with 6-31G** are very similar with the ones obtained by Gutsev

_{b}*et al.*

^{38}and Ma

*et al.*

^{28}especially after certain sizes, our

*E*results with SDD are in good accord with the ones by Cervantes-Salguero and Seminario

_{b}^{37}who studied Fe

_{n}(

*n*= 2 − 10) clusters using the OPBE/TZV level of DFT. The other previous

*E*calculations

_{b}^{40,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.

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 spectra^{13} 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 Δ_{2}*E* by comparing the plots of the data by Kim *et al.*^{43} and Yuan *et al.*^{40} Their *E _{b}* 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 Δ

_{2}

*E*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

*E*data are consistent with each others, as we observe a small peak in our Δ

_{b}_{2}

*E*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 calculations

^{21,38,40,43}along with our own results.

### D. Fragmentation Energies

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 coworkers^{13} 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 Fe_{13} and Fe_{19} clusters correctly, as Yuan and coworkers^{40} did not predict it for Fe_{19}. 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.

### E. Magnetic Moments

In Fig. 6, we compared experimental^{18,36,44} and theoretical^{21,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 Fe_{n} 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 Fe_{10−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 Knickelbein

^{18}argued that the high moments of Fe

_{10−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.

A common feature of all the experimental studies is the clear dip observed in the magnetic moment of Fe_{13} cluster, which is related to the antiferromagnetic coupling between the central and the surface atoms of the Fe_{13} 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 Knickelbein^{18} 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 Fe_{13} 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

*μ*per atom reported by Yuan

_{B}*et al.*

^{40}obtained in their VASP calculations. Thus, the total spin magnetic moment of the icosahedral Fe

_{13}cluster is calculated as 2.62

*μ*per atom, which is very close to the experimental prediction of 2.4

_{B}*μ*per atom by Knickelbein

_{B}^{18}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.

## IV. CONCLUSIONS

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 Fe_{n}(*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 (Fe_{9−13}), HAP (Fe_{14−17}, Fe_{20}), and TDE (Fe_{17(2)}, Fe_{18−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 Fe

_{13}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.

## ACKNOWLEDGMENTS

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).