Structure, electronic properties and vibrational spectra of (MgF2)n clusters through a combination of genetic algorithm and DFT-based approach

In this article, we look at the option of using a stochastic optimisation technique, namely genetic algorithm (GA) in association with density functional theory (DFT) to find out the global minimum structures of (MgF2)n clusters with the range of n being between 2 and 10. To confirm whether the structures are indeed the acceptable ones, we go on to evaluate several properties like IR spectroscopic modes, vertical excitation energy, cluster formation energy, vertical ionisation potential and the HOMO–LUMO gap. We stress on the fact that an initial estimation of structure using GA, on two empirical potentials (with and without inclusion of polarisation), leads to a very quick convergence to structures which are quite close to the structures obtained from quantum chemical calculations done from the outset, such as using a DFT calculation. The general structural trend of these systems to form three-dimensional networks is also clear from our study. The lowest energy isomers of these clusters show preference for four-membered Mg2F2 and six-membered Mg3F3 rings. In the IR spectra of (MgF2)n clusters, a blueshift of the Mg–F symmetric stretch and a redshift of asymmetric Mg–F stretching as n increases are obtained.


Introduction
A correct evaluation of geometrical structures of clusters is important as it leads to the correct estimation of chemical and physical properties of such systems. The study of atomic and molecular clusters is an important area of research for both experimentalists and theoreticians in the last few decades. Nanoscale clusters is an important topic in current scientific research because of the strong dependence of their electronic properties on their size and structure. MgF 2 being an important material with interesting optical properties, (MgF 2 ) n clusters have been studied by theoretical [1][2][3][4][5][6][7][8][9][10] methods. Bulk MgF 2 have been studied by experimental methods [11][12][13][14][15] also. Detailed experimental study [16,17] has been performed on MgF 2 trimer. Such clusters are too small to characterise and their geometry can be analyzed by photoelectron spectroscopy in combination with ab initio calculation [18][19][20][21][22][23]. This technique requires a global minimum search which is based on the assumption that the most energetically stable isomers are likely to form in an experiment and will be the one to be probed in the experimental scenario.
Finding the global minimum of the potential energy surface (PES) [24,25] of an N atomic cluster is difficult because the number of local minima can be very large even for moderate values of N. It must be emphasised here that the lowest energy structure or the global minimum is of paramount importance, but the other low-lying structures * Corresponding author. Email: pinakc@rediffmail.com in addition to the global minimum have their importance too. In a real experimental situation, if the energy barrier between the global and the low-lying structures is not high, the structure which can be sensed experimentally at a finite temperature and in a given finite time is generally an average over all these near-energy conformations. Deciphering the correct structure would essentially mean finding the correct set of interatomic distances for which the potential energy is a minimum. One effective strategy to search out the global minimum in a multiple minima surface is to use non-deterministic stochastic techniques that use concepts of random walk and Markov chains to surmount energy barriers separating different energy conformations and ultimately finding out the global minimum.
For a method to be truly global, it should have the ability to search out the global minimum irrespective of the initial starting point, and in case of being trapped in a local minimum, it should have the ability to jump out of the local attractive basins and proceed on to the global geometry. Nature provides us an insight into powerful strategies and techniques which, if adopted into a mathematical procedure, can really produce efficient and competitive techniques to search really tough surfaces. The way natural evolution works, with selection of fitter individuals followed by occasional crossover of genetic material between relatively fitter individuals and a provision for mutation have over time led to the creation of better individuals. These concepts if incorporated into an algorithm which goes by the name of genetic algorithm (GA) [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] can be an efficient problem solver. There are also a group of time-tested techniques which perform a global search on a PES like Simulated Annealing [41][42][43][44][45][47][48][49] and Basin Hopping [50][51][52][53][54][55][56]. We, in the present communication, would like to use GA as a method of our choice. A mathematical equivalent of selection, crossover and mutation have found good use as a truly global and stochastic search algorithm, which is able to find out low-energy structures of complex large molecules and clusters. GA has been used to solve a variety of problems ranging from structure elucidation of atomic and molecular clusters, solution of Schrodinger equation and predicting geometries of polymeric materials. Cluster systems ranging from Lennard-Jones cluster, Morse cluster and metallic clusters modelled by Gupta potential have been studied, and work of several leading research groups needs a mention. The work of Deaven et al. [57,58], Pullan et al. [59][60][61], Niesse and Mayne [62][63][64], Johnston et al. [65], Michaelian et al. [66][67][68], Chaudhury et al. [69,70] and Hartke et al. [71][72][73][74] deserve special attention.
Francisco et al. [2] used the Basin Hopping technique to search the global minimum of (MgF 2 ) n . Jansen et al. [4] investigated (MgF 2 ) 3 and (MgF 2 ) 4 using Simulated Annealing. The trimer of MgF 2 has been studied in some detail using theoretical [75][76][77][78][79][80] and experimental [81,82] methods. Most of these studies have only considered dimers and trimers of MgF 2 molecule; however, Franscisco et al. [2] reported larger (MgF 2 ) n using the Basin Hopping methodology. The current work is on global optimisation employing GA in conjunction with Density Functional Theory to study (MgF 2 ) n cluster up to n = 10. We would like to see how the IR spectroscopic modes, HOMO (highest occupied molecular orbital)-LUMO (lowest unoccupied molecular orbital) gap, vertical ionisation potential, vertical electron affinity and vertical excitation energy change with increase in cluster size. It must be mentioned that Johnston et al. have implemented a GA + DFT (density functional theory) scheme, where GA is explicitly coupled with electronic structure calculation to find out global minimum structures like an impressive application towards study of Au-Ag nanoalloys [83] and bismuth-doped tin clusters [84].

The methodology 2.1. The potential for the system
In order to describe the interactions that exist between all the atoms constituting a cluster of the type (MgF 2 ) n , where n is the number of Mg atoms, we have used the potential energy function [85] consisting of a two-body coulomb-plus-Buckingham potential. The potential includes the attractive interactions of both the dispersive and coulomb type and also keeps a provision for an exponential repulsive term, which prevents unusually close approach of the constituents making up the cluster. The total energy of the cluster is given by where q i and q j are the charges on ions i and j, respectively, and r ij is the distance between the ions 'i' and 'j'. The parameters are given in Table 1. The charges are Mg = 1.66|e| and F = −0.83|e| and ρ = 0.215Å. We have also obtained structures with an empirical potential which includes the polarisation contribution on F − ion and this has been discussed in detail in Section 3.

Genetic algorithm (GA)
Search surfaces which are extremely rugged demand parallel approaches to be used if the critical points with very low energy, including the global minimum, needs to be found out. Nature can give us guidelines to model such procedures. Natural processes like extremely efficient genetic evolution or the cooperative behaviour of insects can help mathematicians to model optimisation procedures which are highly potent as problem solvers. The tendency to switch over to these so-called natural algorithms have seen an increasing usage in recent years. The way how insects like ants or bees locate their food or shelter has led to the growth of two algorithms, popularly Ant Colony Optimisation (ACO) and Particle Swarm Optimisation (PSO). The techniques are based on the methods used by a colony of insects, involving co-operative behaviour among them, to locate the shortest path towards its target. These algorithms, generally being applicable in a large variety of problems, have found usage in solving problems of structural chemistry where essentially the problem reduces to an in-depth search on a multiple minima surface separated by high energy barriers. These algorithms have an in-built power of coming out of locally attractive basins and continue its search towards a deeper minimum. This is possible because of randomness being an integral part of these search strategies as opposed to the more conventional gradient-based methods which use only the reduction of gradient as the main goal of search. We, in this work, have used a method whose working principle is inspired by the principles of conventional genetics.
Having defined the potential energy of interaction between the constituents of the clusters, the goal now would be to use an appropriate search strategy, which will unubiquitously be able to find out the minimal points supported by the potential energy function and so the structures of the (MgF 2 ) n clusters. As we have already emphasised, searching out minimal points on rugged multidimensional surfaces can be a non-trivial task for any search algorithm. We, in this present problem, have used GA [86] as the search algorithm. This is a stochastic scheme, based on the biological processes of crossover and mutation and the acceptance of the Darwinian concept of Survival of the Fittest. GA is different from most optimisation techniques in its working principles. While most techniques start from a single trial guess solution and progressively improves it, GA follows a different methodology. It generates randomly a large pool of starting trial solutions and subjects it to certain mathematical operations, namely selection, crossover and mutation, which closely mimic the equivalent operations in natural genetics. The starting solution pool as a whole is improved iteratively in successive steps, commonly called generations in GA terminology. As the entire pool of solution is subjected to the changes, as opposed to a single solution, the process to begin with has a large amount of potentially good information to explore and improve upon. This gives GA a high degree of efficiency, which other techniques lack. The essential steps and the working principles in a GA are described in our earlier publications [69,70]. The details of the parameters used in GA simulation are given in Table 2.

Computational details
The optimised coordinates obtained from GA simulation are used as input for geometry optimisation using DFT. In our strategy, all the global and local structures for each cluster size are fully optimised at the B3LYP/6-311++G(d,p) level of theory. All of the quantum chemical calculations reported in this paper were obtained using GAUSSIAN 09 [87] suite of programs. DFT calculations are performed to find out the lowest lying structures of neutral (MgF 2 ) n clusters using the hybrid B3LYP functional [88][89][90] and the 6-311++G(d,p) basis set as incorporated in Gaussian 09 program [87]. All calculations are performed with default DFT grid sizes and convergence criteria. Care was taken so that there are no imaginary frequencies in the vibrational spectrum, which clearly shows that the structures are indeed minimum on the potential energy hypersurface. In order to understand the relative stabilities of these clusters, we have calculated the formation energy E f = (E n /n) − E 1 of the lowest lying structures. E 1 and E n are the total energies of (MgF 2 ) 1 and (MgF 2 ) n clusters after the DFT calculation using the GA-optimised structures as inputs. TDDFT/B3LYP calculations are carried out to evaluate vertical excitation energies to the lowest singlet states. The vertical ionisation potentials (IP v ) and vertical electron affinity (EA v ) are also calculated by single-point B3LYP/6-311++G(d,p) calculations at optimised neutral cluster geometries for singly positively charged and negatively charged clusters. where are the total energies of positively charged, negatively charged and neutral clusters, respectively, without employing zero-point correction. The second difference in energy ( 2 E) is given by , where E n + 1 and E n−1 are the total energies of (MgF 2 ) n + 1 and (MgF 2 ) n−1 clusters. HOMO-LUMO gaps (E g ) is calculated by the following relation:

Structures of (MgF 2 ) n clusters
The monomer of MgF 2 is found to have a linear structure. The computed Mg-F bond length and F-Mg-F angle are 1.78Å and 180.1 • , with the experimental bond length being 1.77Å [91]. For MgF 2 dimer, the global minimum possesses D 2h symmetry. The geometry is characterised by a planar four-membered ring with two bridging fluorine and two terminal fluorine atoms and is in agreement with earlier reported calculations of Francisco et al. [2]. The global geometry of MgF 2 trimer is found to have a double-bridged D 2d structure containing two terminal fluorine atoms, and is also in good agreement with quantum chemical calculations of Francisco et al. The lowest energy structure of (MgF 2 ) 4 has a C 2h symmetry possessing six bridging two-fold coordinated fluorine atoms and two terminal fluorine atoms disposed in a trans fashion. The geometry is highlighted by the presence of two three-fold coordinated Mg atoms and two four-fold coordinated Mg atoms. The global minimum structure of the pentamer is of C s nature. The structure possesses eight bridging fluorine atoms and one three-fold coordinated fluorine atom. Five four-fold coordinated magnesium atoms with one terminal fluorine atom is the noticeable element in the structure. The lowest energy geometry of (MgF 2 ) 6 is of C 2 point group. The geometry is marked by the presence of six four-fold coordinated Mg atoms and one four-fold coordinated fluorine atom in addition to two terminal fluorines. A C s symmetry structure with two terminal fluorines is found to be the most stable structure for (MgF 2 ) 7 . The structure is made up of six four-fold coordinated Mg and one three-fold coordinated Mg. Eleven bridging fluorines with one three-fold coordinated fluorine is also a key element in the structure. The global minimum of (MgF 2 ) 8 is of C s symmetry with the presence of one terminal fluorine. The structure consists of seven four-fold coordinated magnesium and one three-fold coordinated magnesium in addition to 15 bridging fluorines. Francisco et al. reported a C 2v structure as the global minimum which is 0.7 kcal/mole less stable than our structure. For (MgF 2 ) 9 , the global minimum has C 1 point group with 18 bridging fluoride atoms and the presence of nine four-fold coordinated magnesium atoms. The global conformer of the decamer (MgF 2 ) 10 has a very symmetric cage-like structure with a C 2v symmetry. The two terminal fluorines are oriented in a cis fashion with 17 bridging fluorine atoms associated with one tetra-coordinated fluorine atom, and this is the key element in the structure of this size. The structures of global and local minimum (MgF 2 ) n clusters are given in Figures 1 and 2. The relative energies of local structures are given in kcal/mole in Figures 1  and 2.

Spectroscopy of (MgF 2 ) n clusters
The estimation of vibrational frequencies are carried out to confirm that all of the evaluated geometries are indeed minima on the model PES. The vibrational frequencies can be used for comparing experimental identification of these clusters wherever experimental data are available. The IR spectra of (MgF 2 ) n clusters are shown in Figure 3. The F-Mg-F bending mode for the monomer is located at 153 cm −1 , which is in agreement with the experimentally determined value of 160 ± 3 cm −1 [92]. For MgF 2 dimer, an intense band at 762 cm −1 is found to be the asymmetric stretching of the Mg-F terminal bonds. The dimer shows an intense peak also at 264 cm −1 , which again matches with the experimental observation of Neelamraju et al. [5].
The trimer is seen to have a peak at 769 cm −1 contributed by the asymmetric stretching of two terminal Mg-F bonds. The symmetric stretching mode is located at 486 cm −1 and matches well with the experimentally assigned peaks by Lesiecki et al. [93] and Hauge et al. [94]. Bands around 700 cm −1 and higher are mainly due to contributions from asymmetric stretching of the Mg-F bonds. Bands between 400 and 600 cm −1 are a result of Mg-F symmetric stretches. Bonds with lower coordinations are stronger than those with higher coordination, so they give higher vibrational frequencies. The infrared absorption peaks around 580-630 cm −1 are a result of Mg-F-Mg bridging bonds. In the IR spectra of (MgF 2 ) n clusters, a blueshift of the symmetric stretch and a redshift of asymmetric stretching as n increases are found. It is also observed that bonds formed as a result of greater coordination occur at lower frequencies.
The frequency of asymmetric stretching band is found to decrease from 770 cm −1 [(MgF 2 ) 3 ] to 615 cm −1 [(MgF 2 ) 10 ], which is explainable as there is an increase in Mg-F bond length of (MgF 2 ) n clusters with increasing n. Experimental results for these systems are available for n = 2 and 3 only, and for n > 3, our results are in proximity with the trends obtained from theoretical calculations of Neelamraju et al. [5].

Electronic properties of (MgF 2 ) n clusters
To have an idea about how the stabilities of the clusters of various sizes vary, we have proceeded to calculate the formation energy E f (= E n /n − E 1 ) and the nucleation en- Here, E 1 , E N and E N + 1 are the lowest determined energies for the MgF 2 monomer, (MgF 2 ) n and (MgF 2 ) n + 1 clusters, respectively.   stable in comparison to a single MgF 2 . E N depicts the pattern of the growth of (MgF 2 ) n clusters from (MgF 2 ) 1 and (MgF 2 ) n−1 units and is found to have an even-odd oscillation. From Figure 4(b), it is observed that the nucleation energy is less negative for even-sized clusters. The growth of the odd n clusters have been found to be easier in comparison to even ones, and the quickest nucleation process has been found to be an important factor in the production of (MgF 2 ) 9 from (MgF 2 ) 8 and (MgF 2 ) 1 . The second difference ( 2 E) in total energies of the clusters with cluster size (n) is depicted in Figure 4(c) and it has been evaluated using the standard expression [ 2 E = 2E n − (E n + 1 + E n−1 )]. 2 E is a measure of relative stability of the specific clusters (n) with respect to its immediate neighbours (n − 1 and n + 1). Positive (negative) 2 E signifies an unstable (stable) cluster with respect to its immediate neighbours [95]. 2 E shows a clear even-odd oscillation with the increase in n. 2 E is negative for even clusters and this proves the greater stability of even-sized clusters compared to the odd ones. The 2 E exhibits a clearly noticeable minimum at (MgF 2 ) 8 and indicates its highly enhanced stability compared to the immediate neighbours. The energy gap (E g ) between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), the vertical ionisation potential (IP v ) and vertical electron affinity (EA v ) are useful and necessary characteristics, which throw light into the electronic properties of these clusters. The HOMO-LUMO gap is considered to be an important indicator for deciphering electronic stability of clusters. It represents the tendency of a molecule to take part in a chemical reaction to a large extent. A large HOMO-LUMO gap is interpreted as the sign of enhanced chemical stability. The variation of HOMO-LUMO gap with an increased cluster size shows an irregular pattern. With the increase in cluster size, the HOMO-LUMO gap gradually approaches the band gap (9.5 eV) [96] of bulk (MgF 2 ) material. The variation of vertical ionisation potential (VIP) with cluster size is shown in Figure 4(e). The plot shows a downward surge at n = 6 and then a gradual increase. This may be attributed to the fact that the larger clusters having a higher fluorine content offer a greater reluctance to the knocking off an electron. The vertical electron affinity (VEA) shows an upward trend initially with the gradual increase in cluster size. The reason for this behaviour is that the clusters which have an excess of fluorine atoms can take up excess electrons easily and thus show high electron affinity values. The vertical excitation energy (E ve ), i.e., the energy needed to go from the ground singlet to first excited singlet state shows a decreasing trend for smaller size clusters with a clearly noticeable minimum at n = 5 and then an increasing trend though not in a very regular way. The calculated electronic properties are given in Table 3.

Discussion on the quality of potential used for (MgF 2 ) n clusters
The potential, which we have used to calculate the preoptimised structures before using them as inputs for subsequent quantum chemical calculation contain coulomb interactions (between point charges) and a repulsive part which grows exponentially at shorter distances. Logically, the question that might arise is whether this choice is sufficient or can be modified. Obviously, the next term which can be included is due to polarisation forces, especially the contribution from the polarisation of F − ion. The Mg ion has low polarisability and is expected not to contribute much. Therefore, one can modify the empirical potential energy expression by adding the term where − − → µ ind k is the induced dipole moment of the k-th ions (the k sites are F − only) and −−→ E coul k is the electric field present at ion 'k' because of the charges only. and where − → Q kj and ↔ T kj are the corresponding propagators arising out of coulomb and dipole interactions. With the inclusion of this V pol term in the empirical potential [Equation (1)], we have proceeded on to check if the global minimum structures do change or not compared to where V pol was not present. We have presented the 10 reevaluated global minimum structures (MgF 2 ) n (n = 1-10) as supplementary material. A comparison of any two structures of same n shows that there is a shortening of Mg-F bonds by a slight amount, which is expected as an additional attractive term that is now present in the modified empirical potential. However, the geometry and shape of the cluster remains unaltered which is also quite evident. At this point it must be mentioned that, since the structures obtained with the inclusion of polarisability term does not alter the structures obtained from the empirical potential compared to the situation where polarisability was not included, the final converged structure with the quantum chemical calculation remains the same. In Figure 3, the infrared spectrum of the studied systems has been depicted. Since the final converged structures with the two sets of inputs from search with empirical potentials (without polarisation and with polarisation) remain the same, we do not observe any change in the final infrared spectrum. The same is true for Figure 4(a)-(f). However, one can get an estimate of the pattern of change in the formation energy (E f ) even after calculation with the empirical potential. This is shown as inset figure in Figure 4(a). The variation pattern is in a similar line with the quantum chemical result [not shown as inset in Figure 4(a)]. However, the absolute values of the formation energies are obviously different as they are evaluated at two different levels of theory.

Conclusion
A hybrid GA-DFT method has been used to generate the lowest energy structures of (MgF 2 ) n clusters with n = 1-10. We have studied the structures, spectroscopy, growth behaviour and electronic properties of neutral (MgF 2 ) n clusters in the size range n = 1-10. The interconnection between electronic and structural properties has been studied by considering the variations of clusterformation energies (E cf ), HOMO-LUMO gap (E g ), vertical ionisation potential (IP v ) and vertical excitation energy (E ve ) with cluster size (n). According to the calculated E cf , E N and 2 E values, even n clusters are found to have higher stability than odd ones. From the second-order difference of total energy ( 2 E), we find that the n = 5, 7 and 9 cluster is thermodynamically unstable with respect to disproportionation to smaller and larger clusters. In general, the lowest lying structures of n ≤ 9 have at least one or two terminal fluorine. However, the structures n = 9 and 10 do not have any monovalent fluorine atom in the calculated lowest energy structure. The coordination of Mg and F atom is an important criterion to explore the stability of the clusters. With the increase in size, even n clusters tend to form compact structures that are more stable with larger HOMO-LUMO gaps. There is a very close correspondence of the computed infrared modes to the experimental frequencies.
For lowest lying structures, characteristic infrared peaks above 700 cm −1 are due to terminal Mg-F bonds and below 680 cm −1 due Mg-F-Mg bridging bonds. The bands within 400-600 cm −1 are due to four-fold coordinated Mg atoms. Characteristic peak close to 630 cm −1 can be useful for probing surface bridging Mg-F-Mg structures. Our vibrational spectroscopic calculations closely follow the experimentally available IR peaks for n = 2 and 3 and for n > 3, where only theoretical results are available, we also notice a close correspondence.