Theoretical Foundations of the Models Implemented in LUMPAC

Process of Geometry Optimization

The potential energy surface (PES) consists of a surface of the calculated energy () as a function of the geometric parameters of a molecule (q), , where n is the number of geometric parameters. A steady point in the PES is defined by a flat point on the surface, being mathematically represented by . If , the steady point is a transition state. In contrast, if , the steady point corresponds to a minimum of energy, usually a local minimum. The lowest of the local minima, the global minima, usually defines the most stable ground state geometry of the molecule.

In general, the geometry optimization procedure consists in supplying an input molecular structure with geometrical parameters hoped to be as close as possible to the steady point desired. This reasonable geometry is then submitted to an algorithm which systematically alters the atomic positions until a minimum steady point is reached, defined by its geometric parameters .

The computational chemistry methods that are applied to perform geometry optimizations may be divided roughly into two groups: molecular mechanics methods and quantum methods. In molecular mechanics methods, roughly speaking, the molecule is treated as set points, each representing an atom, interconnected by springs, representing the chemical bonds. Such methods are very fast because the potentials are classic and no electronic wave functions are present. In contrast, quantum methods attempt to solve, for the entire molecular system, the famous SchrΓΆdinger eigenvalue equation , where is the Hamiltonian operator, is the wave function (the eigenvector), and E is the total energy of the system (the eigenvalue).

The quantum computational methods are usually divided into four groups: i) the Hartree-Fock methods which solve the self-consistent field SchrΓΆdinger equation under the independent particle approximation, calculating all possible two-electron integrals; ii) the post Hartree-Fock methods, which also take into consideration the contributions due to electron correlations; iii) the semiempirical quantum chemical methods, which are based on the Hartree-Fock method, though some integrals are replaced by parameters, adjusted to reproduce experimental data during its process of development; and, finally, iv) the methods based on density functional theory (DFT), which consider the electron density as the fundamental entity, and not the wave function as all other quantum methods do.

The method of choice for performing geometry optimizations depends mainly on the number of atoms present in the desired system.

Normally, the standard procedure for geometry optimization consists in building the chemical structure by adding atoms in arbitrary positions and by subsequently connecting them according to their chemical bonds. The next step is to previously optimize the geometry by using less costly computational methods. Molecular mechanics methods are quite fast methods and some of them have parameters available for almost all elements. As a result, although molecular mechanics methods do not provide accurate geometries, such methods are usually applied to transform the drawn geometries into reasonably good starting structures for refinement by computationally more accurate and usually more expensive calculations [1].

The ground state geometries of lanthanide complexes can be calculated by two different quantum chemical based approaches: i) DFT or ab initio methods with effective core potentials (ECP) for treating the lanthanide ions, or ii) semiempirical methods. In 2006 [2] and 2011 [3] two papers were published in order to compare geometries predicted by semiempirical methods with those predicted by ab initio and DFT methods, with crystallographic data as reference. In contrast to what would be expected, the results showed that by enlarging the size of the basis set, or by including electron correlation, or both, deviations of the predicted coordination polyhedrons with respect to the crystallographic ones consistently increased, reducing the quality of the results. And among all abinitio methods evaluated, the method RHF/STO-3G using the MWB core effective potential was the most efficient for predicting the coordination polyhedron of lanthanide complexes [2]. This result confirms that the Sparkle models, which demand a much lower computational effort, have higher accuracy in calculations and modeling, when compared to ab initio/ECP ones [3].

In LUMPAC, the Sparkle models will be used in the geometry optimization step because such methods have an excellent capability of geometry prediction and also a considerably low computational cost. As a result, we will now provide a detailed description of these models. The procedure of development of the Sparkle model consists in parameterizing a semiempirical Hamiltonian, such as AM1 or PM3, for example, in which the lanthanide ion is replaced by a +3e point charge. This point charge is subjected to a repulsive potential , where the parameter quantifies the size of the ion. This mathematic entity is called βsparkleβ. As the bond between Ln^{3+} and atom ligands has high ionic character, the Sparkle model has been consistently proven to be adequate.

The first Sparkle model, named SMLC (Sparkle Model for the Calculation of Lanthanide Complexes), was developed by Andrade and coworkers in 1994 [4]. This version was parameterized for the AM1 semiempirical model with only one experimental structure in the parameterization set: the tris(acetylacetonate)-(1,10-phenanthroline) of europium (III). When this Sparkle model version was evaluated for a representative test containing 96 europium complexes, the SMLC model lead to errors of approximately 0.68 for LnL, lanthanideligand atom, distances. However, the second parameterized version of the Sparkle model [5], SMLC II, published in 2004, included Gaussian functions in the core-core repulsion energy. The errors for LnβL distances decreased from 0.68 to 0.28 when tested with the same europium structures set.

A new and much more sophisticated parameterization scheme was then carried out within AM1 for the Sparkle model in 2005 and was initially developed for Eu^{3+}, Gd^{3+}, and Tb^{3+}[6]. This new version of the model was called Sparkle/AM1. The main changes consisted in the application of more sophisticated statistical techniques, both in the selection of the most representative training sets as well as in the validation of the parameters obtained. In the Sparkle/AM1 model development of the three parameterized ions, more than 200 different crystallographic structures were used together with a new response function for minimization in the parameterization procedure. These changes made it possible to decrease the errors for LnβL distances from 0.28 to 0.09 in europium complexes. Test sets of gadolinium complexes (70 structures) and of terbium complexes indicated errors of approximately 0.07 . Then, the Sparkle/AM1 model was generalized for all types of ligands and parameterized for all 15 trivalent lanthanide ions [7-14]. Currently, the Sparkle models are also parameterized for the following semiempirical models: PM3 [15-21], PM6 [22], PM7 [23] and RM1 [24].

The geometry optimization of lanthanide complexes has a great importance for studying the luminescent properties of the system. All published Sparkle models (Sparkle/AM1, Sparkle/PM3, Sparkle/PM6, Sparkle/PM7 and Sparkle/RM1) are available in MOPAC2012 [25]. The choice of which of the available Sparkle models is to be used must be based mainly on the capability of the underlying semiempirical method, either AM1, PM3, PM6, RM1 or PM7 to correctly describe the specific ligands involved. Nevertheless, many tests performed by our group suggest Sparkle/RM1 to be the version that presents the best overall results.

Excited States Calculation

The singlet and triplet excited states of the organic part can be calculated by using methods based on time-dependent density functional theory (TD-DFT) [26] or by the semiempirical INDO/S method [27, 28]. In 2001, Gorelsky and Lever [36] compared these two methodologies for the ground and excited states calculations of Ru(II) complexes. The electronic spectra obtained by these two different methods showed excellent agreement with each other. However, even today, the TD-DFT method is still inappropriate to treat complexes with more than 100 atoms, due to its high demand of computational resources.

Santos and coworkers evaluated the accuracy of the semiempirical INDO/S method in comparison with TD-DFT ab initio results in studies of lanthanide complexes [29]. The results showed that triplet state energies calculated by the semiempirical method presented errors similar to those obtained by TD-DFT methodology, with the advantage of being hundreds of times faster.

In this context, the geometries optimized by the Sparkle models are used to calculate the singlet and triplet excited states by using the configuration interaction simple (CIS) of INDO/S, which has an accuracy of about 1000 cm^{-1} [27, 28]. This method is implemented in ZINDO [30] and ORCA [31] programs. In this procedure, a point of charge +3e represents the lanthanide ion [32].

Intensity Parameters Calculation

The intensity parameters, ( = 2, 4, and 6), are calculated by Judd-Ofelt theory [33, 34]. According to this theory, the central ion is affected by the nearest neighbor atoms, through a static electric field also referred as crystal or ligand field. Judd and Ofelt described, in independent works, the importance of the electric dipole mechanism for the 4f4f transitions from the mixing of a ground state 4f^{N} configuration with excited state configurations of opposite parity through the odd terms of the ligand field Hamiltonian. All 4f orbitals have the same parity, that is , where = 3 for lanthanide ions. Then, the mixing are between the 4f orbitals plus higher-n orbitals, such as the 5d orbital, which presents = 2, and has an opposite parity to that of the f orbital.

The intensity parameters describe the interaction between the lanthanide and ligand atoms, and are calculated by Eq. ().

One aspect which is very important for the possible application of this theory is to know the values that each of the rank variables , , and may assume in relation to each other. As can been seen in Eq. (), for example, when is equal to 2, will be equal to 1 and 3, whereas the values of will be equal to 0, 1, ..., .

The parameters are calculated by:

The first term, , refers only to the forced electric dipole (ED) contribution, and are given by Eq. ().

The term corresponds to the difference of energy between the ground state barycenters and the first excited state configuration of opposite parity. The radial integrals, , were taken from reference [35], with an extrapolation for the quantity . The values of radial integrals for Eu^{3+} ion are = 0.9175 a.u., = 2.0200 a.u., = 9.0390 a.u., and = 110.0323 a.u. The terms are numeric factors associated with each lanthanide ion and are estimated from radial integrals of Hartree-Fock calculations [36]. The values of are ; ; ; ; , and , for Eu^{3+} ion [36].

The second term of Eq. (), , refers only to the dynamics coupling (DC) contribution and is given by Eq. (). This contribution is complementary to the one from the Judd-Ofelt static electric field model, and was firstly considered by Mason and coworkers [37]. The dynamics coupling mechanism, which is more important than the electric dipole mechanism for some transitions, is due to the high gradient of the electromagnetic field generated by the ligands when they interact with an incident external field. The DC mechanism depends on the nature of both ligands and on the coordination geometry, and explains the hipersensitivity in 4f4f transitions [36].

The quantity is a shielding field due to 5s and 5p filled orbitals of lanthanide ions, which have a radial extension larger than those of the 4f orbitals [36]. The values of are 0.600, 0.139 and 0.100. is a tensor operator of rank ( = 2, 4, and 6) with values = -1.366, = 1.128, and = -1.270 for lanthanide ions. is the Kronecker delta function. As such, is equal to 0 when is different from the .

The parameters ( = 1, 3, 5, and 7), given by Eq. (), are the so-called odd-rank ligand field parameters and contains a sum over the surrounding atoms.

are the conjugated spherical harmonics. As can be observed in Eq. (), the spherical harmonics depend on the spherical coordinates of the j ligand atoms.

The term present in Eq. () according to the Simple Overlap Model (SOM) [38, 39] developed by prof. Oscar Malta (UFPE, Brazil), formalizes that crystal field Hamiltonian and is adequately calculated as a function of the charge density between the lanthanide ion and the ligand atoms.

The SOM model assumes two postulates [38]: i) the 4f energy potential is generated by charges, uniformly distributed in a small region located around the mid-points of the lanthanideβligand chemical bonds; and ii) the total charge in each region is equals to , where the parameter is proportional to the magnitude of the total overlap between the lanthanide ion and the ligand atoms. Figure 1 shows a sketch of the effective charges for a hypothetical complex (LnL_{3}). The vector represents the position of the ligand atoms, and the vector represents the position of the electron of the central metal ion.

Figure 1. Graphical representation of the Simple Overlap Model.

In other words, the term introduces a correction to the crystal field parameters of the point charge electrostatic model (PCEM), , such that . This way, this correction confers a degree of covalency to the point charge model from the inclusion of the parameter , since PCEM treats the metal-ligand atom bonds as a purely electrostatic phenomenon.

The effective charges are assumed to be at positions defined at the distances given by . The factor , given by Eq. (), indicates that the effective charges may not be exactly at . The plus sign in Eq. () is used when the barycenter of the overlap region is displaced towards the ligand, which happens in the case of oxygen and fluorine coordinating atoms. The minus sign is used when this barycenter is displaced towards the central ion, as is the case of nitrogen and chlorine coordinating atoms.

The overlap between 4f orbitals and the valence orbitals of the ligands, , is calculated by Eq. ().

where is a constant equal to 0.05 and is equal to 3.5 for the lanthanides. is the smallest among all lanthanideβligand atom distances.

The parameters ( = 1, 3, 5, and 7), like the parameter , also depends on the coordination geometry and on the chemical environment around the lanthanide ion, and is given by Eq. ().

The limitations in the intensity parameters calculation consist in determining the quantities, and . As a result, it is necessary to use the experimental intensity parameters.

The charge factors and polarizabilities, used in and calculations, respectively, are adjusted to reproduce the experimental intensity parameters. During the adjustment procedure, the intensity parameters calculated () from the optimized geometry, obtained from Sparkle model, are compared with the experimental intensity parameters (). The response function () is defined by Eq. ().

Emission Radiative Rate Calculation

The emission radiative rate (), taking into account the magnetic dipole and forced electric dipole mechanisms, is given by Eq. ():

where is the difference of energy between the ^{5}D_{0} and ^{7}F_{J} states (in cm^{-1}), is the Planck constant, is the degeneracy of the initial state, and is the refractive index of the medium, usually assumed to be equal to 1.5. (Eq. ()) and in Eq. () are the magnetic dipole and forced electric dipole mechanisms, respectively.

The squared matrix elements , , and are equal to 0.0032, 0.0023, and 0.0002, respectively, for Eu^{3+} [40].

where is the electron mass. The matrix elements that appear in Eq. () above are determined according to the intermediate coupling mechanism. The ^{5}D_{0}^{7}F_{1} transition is the only one that does not have contributions from the electric dipole mechanism and are quantified theoretically as esu^{2} cm^{2} [41].

The ^{5}D_{0}^{7}F_{J} transitions ( = 0, 3, and 5) are forbidden by magnetic dipole and forced electric dipole mechanisms, that is, their contributions are equal to 0.

The contributions of each transition to the emission radiative rate are calculated by Eq. (), and are named branching ratios ().

Energy Transfer Rates Calculation

The theoretical model used to calculate the energy transfer rate between the organic ligands and the lanthanide ion was developed by Malta and coworkers [42, 43]. According to this model, the energy transfer rates, , are given by the sum of two terms:

The term , given by Eq. (), corresponds to the energy transfer rate obtained from the multipolar mechanism.

The quantities are the electric dipole contributions to the intensity parameters, taking into account only the contributions of the parameters . are reduced matrix elements of the tensor operators . The parameters Ξ³_{Ξ»} are calculated by Eq. ().

In Eq. (), is the total angular momentum quantum number of the lanthanide ion. is the degeneracy of the initial state of the ligand, and a specifies the spectroscopic term of the 4f orbitals. is the dipole strength associated with the transitions in the ligands. The quantity , calculated by Eq. (), corresponds to the temperature-dependent factor and contains a sum of Frank Condon factors.

The factor (Eq. ()) is the ligand state bandwidth-at-half-maximum (in cm^{-1}), and is the difference of energy between the donor and acceptor states involved in the energy transfer process. For lanthanide complexes, the energy donor states correspond to the singlet and triplet excited states, whereas the acceptor states correspond to the lanthanide ion excited states. Typical values of are in the range of 10^{12}10^{13} erg^{-1}.

The second term of Eq. (), refers to the energy transfer rates obtained from the exchange mechanism and are calculated by Eq. ().

In Eq. (), ( = -1, 0, 1) is the spherical component of the spin operator of electron in the ligand. The is the component of its dipole operator, and is the total spin operator of the lanthanide ion. Typical values of the squared matrix element of the coupled dipole and spin operators lie in the range 10^{-34}10^{-36} esu^{2} cm^{2} [44].

Malta proposed some corrections for the energy transfer rates equations in 2008 [44]. The first one corresponds to the addition of a shielding factor, to the first term of Eq. () . This contribution had been initially neglected for dipole-dipole mechanism. The second one corresponds to the replacement of the quantity in Eq. () by the overlap integral. This last correction causes a change of three orders in the magnitude of the energy transfer rate calculated by the exchange mechanism. Nevertheless, as these rates are still much higher than the radiative and non-radiative rates, the general conclusions for the theoretical quantum yield obtained from previous work remain valid [44].

The energy transfer rates depend on the distance difference between donor and acceptor states involved in the process of energy transfer. This distance is known as , and for its determination, it is necessary to estimate the molecular orbital coefficients of the atom () that contributes to the ligand states (triplet or singlet). It is important to know the distance from atom () to the lanthanide ion. The quantities and are calculated by data obtained from excited states calculations using semiempirical INDO/S method. This way, the is given by Eq. ().

The energy back-transfer rates () are obtained by multiplying the transfer rate () by the Boltzmann factor, , considering the room temperature. refers to the energy difference between the donor and acceptor levels, and is Boltzmann constant.

The most important transfer channels for systems containing europium ion, according to Malta [45], are shown in Figure 2.

Figure 2. Transfer channels involved in the energy transfer rate processes of systems containing europium ion.

The total angular momentum selection rules, , of the lanthanide 4f states, are complementary. The europium excited states that are more likely to accept energy from ligands through the direct Coulomb interaction mechanism are ^{5}D_{2}, ^{5}L_{6}, ^{5}G_{6}, and ^{5}D_{4}. The energy transfer from ligand excited states to the ^{5}D_{1} level is allowed by the exchange interaction mechanism (Eq. ()). Although the energy transfer to the ^{5}D_{0} level, in principle, is forbidden by direct interaction or exchange mechanism, the selection rule can be relaxed by a mix of the total angular moments (βs) [43, 45].

Emission Quantum Yield Calculation

The emission quantum yield, given by Eq. (), is defined as the ratio between the emitted and absorbed light intensities.

where is the ^{5}D_{0} level population. and correspond to the S_{0} singlet level population and absorption rate, respectively. The normalized population levels, , are obtained from the appropriate rate equations given by Eq. ().

From Eq. (), or represent the transfer rates between and states, or and states.

The ^{5}D_{0} level population depends on the non-radiative emission rate , which still cannot be theoretically calculated. However, can be quantified via Eq. () from the and the experimental lifetime (). Because of this, the theoretical emission quantum yield depends on the experimental lifetime.

The normalized populations of the states involved in the process of energy transfer are obtained from diagonalizing the matrix that contains the rate equations showed in Figure 3. As can been in Figure 3, the matrix is assembled from energy transfer and back-transfer rates, and . The matrix diagonal contains the transfer channels responsible for the energy depopulation of the states in the matrix columns. The channels in red (Figure 2) are not normally included in the energy transfer rates diagrams (Figure 2) due to the non-resonance condition presented between some ligand and the europium excited states.

The emission quantum yield is then calculated from the population given by Eq. ().