M.B. Darkhovskii ^1,3, I.V. Pletnev ^2,3 ,
and
A.L.Tchougréeff^1,3
1 L.Y. Karpov Institute of Physical Chemistry
² Chemistry Department, M.V. Lomonosov Moscow State University
³ Center for Computational Chemistry at the
M.V. Keldysh Institute for Applied Mathematics

Low- and high-spin iron (II) complexes studied
by effective crystal field method combined with molecular mechanics.

Abstract

1  Introduction

Exploring potential energy surfaces (PES) of organic molecules is one of the main application fields of molecular mechanics (MM) []. Conventional MM scheme is, however, hardly applicable to the metal ion complexes where the main problem is posed by an adequate modeling of coordination sphere, in particular, accounting for the flexibility of coordination polyhedron and, more generally, for the interplay of steric and electronic factors. Existing MM schemes consider any atom in a molecule as one having a definite ideal ('strain-free') stereochemistry, which is allowed to be distorted only slightly. While being reasonable in organic chemistry, this construction is not valid when applied to most coordination compounds; it precludes any analysis of strongly distorted structures and/or structural rearrangements.

Physical pre-condition of validity of the MM description for organic molecules is that their electronic excited states are well separated from the respective ground state on the energy scale []. In these cases a single (ground) quantum state of electronic system suffice to describe a molecule and MM becomes valid. By contrast, behavior of the metal d-shell is mostly quantum, since several electronic states may be observed in a narrow energy range close to its ground state. Sometimes the PES's corresponding to different electronic terms of the d-shell intersect which results in spin transitions [,]. The natural way to theoretically handle the situation is to integrate a relevant QM description for the central atom and its closest surrounding (first coordination sphere) with the MM.

The problem of MM description of transition metal complexes (TMC) is thus a special case of a more general problem of constructing hybrid quantum mechanics/molecular mechanics (QM/MM) methods. In Ref. [] a general approach to description of molecular electronic structure and potential energy of such combined systems was suggested with proper formulation of the QM/MM junction part of the total energy, which consists of the field, polarization and resonance (covalent) contributions to the intersystems' interaction.

An approach fitting the above scheme was proposed and tested in Ref. TchDarh for the PES's of the spin-crossover Fe(II) complexes. It allows to calculate the energies of lowest states of the d-shell in the crystal field induced by the ligands. These energies are combined with MM energy of the ligand sphere in order to get the estimate of the total energy of the TMC. The physical preconditions of validity of the effective Hamiltonian for crystal field (EHCF) method [] used to evaluate the crystal field at the d-shells are opposite to those of the standard MM schemes. This approach is physically relevant for the situations when the d-shell of the TMC is a relatively isolated electron group which fundamentally retains the system of multiplets characteristic for the free transition metal ion in the respective oxidation state. This situation may be described by contrasting it to, say, methane molecule whose electronic structure manifests no trace of the multiplet system of the free carbon atom: in the hybrid states of carbon many states even of different electron count are sufficiently mixed []. Another validity precondition is that the lower excitation states of the entire TMC are those of the d-shell itself. These two limiting conditions single out the complexes we are going (and are able) to treat. From the point of view of an experimentalist they can be characterized as those fitting into the crystal field theory. We call them rather arbitrarily the ''Werner-type TMCs''. In fact this approach concerns largely the complexes of the first transition row metals in their common oxidation states II and III. The features of their electronic structure and of the spectrum of their electronic states sketched above sufficiently distinguish them from other types of transition metal compounds, like first row carbonyls or complexes of second and third transition row. This explains the reasons why the concepts of Ref. [] extracted from analysis of numerical experiments performed largely on complexes of heavy transition metals do not apply to TMC considered in the present paper in full. However, some important elements of the bonding picture in TMCs described in Frenking work for the Werner-type complexes as well which will be discussed below.

In the EHCF approximation the energy of the d-shell depends both on the molecular geometry (position and orientation of the ligands with respect to the central metal ion) and on the electronic structure of the ligands. The method [] is computationally intensive since it requires calculation of electronic structure of the whole ligand sphere at each step of geometry optimization. In Ref. [] a combination of the local version (Ref. []) of the EHCF method (EHCF(L)) with the special version of the MM scheme - MMGK procedure (Ref. PletMel,RazPlet) has been proposed and implemented with parameters estimated for a series of Fe(II) complexes with nitrogen containing ligands. In the present paper, we report on further improvements of the approach DarhRazPletTch and on application of the hybrid method to studying the geometry of a series of Fe(II) complexes extended and optimized with use of a procedure different from that of Ref. []. A sufficient point is that the complexes considered here have different values of the ground state total spin.

The paper is organized as follows: in the next Section we briefly review the basic features of the EHCF(L) method DarhRazPletTch,DarhTchRJPC,SouTchMis allowing to determine the crystal field in terms of the Green's function of the lone pairs of the ligands bound to the metal atom. It is shown that the Green's function can be divided into separate contributions from the lone pairs of the free ligands with perturbative corrections coming from the Coulomb interaction of each ligand with the rest of the complex. The last Section provides the applications of EHCF(L)/MMGK to some problems close to ''real-world'' ones.

2  Hybrid EHCF(L) - MM model

The key point for incorporation of transition metal ions (TMI) into MM is to estimate the energy of the d-shell as a function of the ligand sphere's composition and geometry. In this Section we review the working approximation based on the EHCF(L) theory Ref.[] performing this task and the EHCF(L)/MM junction procedure for the MMGK version of the MM method.

2.1  Basics of EHCF(L)

The EHCF method is based on a general concept of separating electron variables. The same concept applies when a hybrid QM/MM method is to be developed []: electrons have to be divided into groups, some of the groups whose excited electronic states are accessible in the experiment are treated quantum mechanically whereas the behavior of other groups whose excited electronic states lay high in energy are modeled with use of MM. In a TMC comprising one TMI and ligands around it the basis of valence atomic orbitals (AO's) containing the 4s-, 4p-, and 3d-AO's of the metal atom (for a first transition row element) and those of the ligand atoms, is according to Ref. [] divided into the d-system which contains only 3d-orbitals of the TMI and the l-system which contains 4s-, 4p-AO's of the TMI and the valence AO's of the ligand atoms. In the EHCF method [] it is shown that the effective QM Hamiltonian H_d^eff for the d-shell has the form:

Hdeff = å mns  Umneffdms+dns+ 1 2 å mnrh  å st  (mn | rh)dms+drt+dhtdns

(1)

where the d-electron Coulomb interaction term is inherited from the free metal ion and the effective core attraction parameters U_mn^eff contain contributions from the Coulomb and the resonance interaction of the d- and l-systems:

Umneff = dmnUdd+Wmnatom+Wmnfield+Wmncov,

(2)

where

Wmnatom = dmn( å a Î s,p  gmaPaa)

(3)

is the repulsion of electrons in the d-shell from those in the 4s- and 4p-AO's of the metal;

Wmnfield = å L  QLVmnL

(4)

is the Coulomb interaction of d-electrons with the net charges Q_L on the ligand atoms, having the standard crystal field theory form []. The covalence part:

Wmncov = å L  å L Î L   bmLbnLGLLadv(Ad)

(5)

ultimately comes from the resonance interaction between the d- and l-systems. In eq.(5) the index L enumerates the ligands, the subscripts L enumerate the one-electron local states referring to the lone pairs (LP's), residing on the donor atom. According to Ref. TchPCCP in the general QM/MM context the formulae eqs.(3)-( 5) represent the exact form of the intersubsystem junction, which depend on the electron structure parameters (ESPs) of the ligands. The ESP's are the one-electron densities P_aa; (the effective charges are derived from them), the orbital energies e_i and the expansion coefficients of the one-electron local states (see below) involved in the resonance interaction with the d-shell. The latter ESP's are condensed in [,] the advanced Green's function G^adv for the ligands given by

GLLadv(e) = -\dsum\limitsi\dfracniciL2e-(gdi-ei)

(6)

where c_iL is the coefficient of the local state (LP) expansion over the canonical MO's (CMO's), which are obtained by the max Y⁴ localization procedure []; g_di is the interaction energy between d-electron and electron on the i-th CMO, and e_i is the i-th CMO energy.

The resonance integrals b_mL in eq.(5) are expressed through the vector t^L formed by the resonance integrals between the metal d-AO's and the L-th LMO taken in the diatomic coordinate frame (DCF) related to the ligand (donor atom) L :

bmL = å \Sb lL Î L \endSb RlmL tl L.

(7)

The coefficients R_lm^L form a unitary matrix R^L transforming d-orbitals from the global (laboratory) coordinate frame (GCF) to the DCF.( The latter is chosen in such a way that its z-axis is the straight line connecting the metal atom with the ligand donor atom). With use of these quantities the covalence part eq.(5) can be rewritten in a rotationally invariant form:

Wmncov = å L  å ll¢  RlmL ell¢L Rl¢mL ell¢L = å L Î L   tlLGLLadv(Ad)tl¢L+,

(8)

where the quantities ell¢L form a eL matrix. The expression eq. (8) allows to calculate the ell¢L parameters with use of the EHCF(L) method. Their relation with the standard angular overlap model (AOM) parameters is described in details in Refs. [,]. The EHCF(L) method has been explored in Refs. DarhTchRJPC,DarhRazPletTch. It has been shown to be able to reproduce the splitting parameter 10Dq of d-shell levels with a 10% precision DarhTchRJPC.

2.2  Total energy in hybrid EHCF(L) / MM model

The total energy of a TMC in its n-th electronic state in the EHCF(L)/MM approximation is taken as in Ref.[] where it is shown to be:

En = EL+Edeff(n) Edeff(n) = áYnd|Hdeff|Yndñ

(9)

The term E_d^eff(n) is the d-shell energy calculated as the n-th eigenvalue of the effective d-shell Hamiltonian eq.(1). The ligand energy E_L is replaced by E_MM, the MM energy of the ligands. In the present work we assume that the effective d-shell Hamiltonian is estimated by the EHCF(L) method described in the previous Section. The contribution E_d^eff(n) apparently is not a MM-like 'force field' and has a different nature.

2.3  Approximate treatment of the ligand electronic structure parameters. Rigid ligands' MO's (RLMO) model

For evaluation of the effective Hamiltonian for the d-shell, the ESP's of the l-system must be inserted in eqs. (1)-(4), (5), (8). These ESP's are contained in the l-system Green's function. In the present Section we discuss general formulae which comprise perturbative evaluation of the Green's function of the l-system using those of the separate free ligands as a zero approximation. Inserting this approximate form of the l-system Green's function into the EHCF(L) formulae eqs. (5), (8) yields the required estimate for the crystal field acting on the d-shell of a central TMI in terms of the separate increments of the lone pairs for each molecular configuration of the TMC. Within this model the influence of the central ion on the surrounding ligands reduces to that of the Coulomb field DarhRazPletTch. In terms of the corresponding Green's function, the Coulomb field affects the positions of the poles of the Green's function (orbital energies) of the free ligand. According to [], the effect of the Coulomb field upon the orbital energies can be represented by:

(GL )-1 = (G0L )-1-S(f)

(10)

where G₀^L is the Green's function for the free ligand with poles corresponding to the free ligand's MO's energies:

(G0L (e))LLadv = å i Î L   lim d® 0+  (ciLL )2ni e-eL i(0)+id

(11)

and the self-energy term S^(f) is due to the external (ion) Coulomb field. The perturbed Green's function G^L within the first order has the same form as G₀^L eq. (11) but with the poles e_{L i}:

eL i = eL i(0)+Sii(f)

(12)

The self-energies S_ii^(f) are taken as those of pure electrostatic interaction between the partial electron densities and point atomic charges:

Sii(f) » \dsum\limitsN Î L riNL dhN

(13)

where r_iN^L is the partial electron density of the i-th CMO of the ligand L on the N-th atom of the ligand:

riNL = å a Î N  ( ciaL ) 2

(14)

(here c_ia^L are the MO LCAO coefficients of the free ligand) and the core Hamiltonian perturbation dh_N is:

dhN = -e2( \dfrac(ZM-nd)RN+\dsum\limits\Sb L ¢ ¹ L N¢ Î L ¢\endSb \dfracQN¢RNN¢)

(15)

The perturbations dh_aa of the core Hamiltonian matrix element in the AO basis are taken to be equal to the corresponding atomic quantity dh_N. The quantities Q_N remain fixed through calculation and do not depend on the complex structure.

These formulae and approximations comprise the RLMO model of the electronic structure of the l-system of the TMC. Despite its empirical implementation it reflects some important features of bonding characteristic for TMCs Frenking, namely the Coulomb interaction between metal ion and ligands. The work on including further elements of electronic structure into this scheme is under way now. The RLMO procedure of estimating the parameters of the l-system electronic structure has been implemented in the program suite EHCF(L)/MMGK []. The results of its application to analysis of molecular geometries of Fe(II) complexes are given below.

2.4  Parameters used in the EHCF(L)/MMGK approaches

The above EHCF(L)/MMGK method with the RLMO procedure (EHCF(L) (RLMO) / MMGK) for the ligands' ESP's in general terms is (as it has been mentioned in the Introduction) a specific case of a general hybrid scheme including QM and MM components which both require special parameterization. The entire set of parameters consists of three subsets. In our case these are the subsets related to the QM description of the d-shell, the parameters of the MM part and those relevant to the junction between the MM and QM subsystems.

2.4.1  d-shell parameters

The d-shell parameters are taken from the EHCF method Ref. [] without changes. These are the atom specific exponents of d-orbitals and d -electron core attraction parameter U_dd. The Coulomb repulsion of d-electrons is characterized by three parameters: g_dd, and the Racah parameters B and C which are either specific for the complex (if known) or are standard for the free ions tabulated, say, in Ref. []. All these parameters are described in details in Ref. []. In the present paper we use the standard ionic values of B and C, though may be uniformly scaled for all complexes considered (see below).

2.4.2  MM parameters

Organic part of a molecule and metal ion coordination sphere (leaving out effects of the d-shell) in the present hybrid procedure is described in the frame of the MMGK method [,]. The latter is the most straightforward attempt to include metal ions in the MM model and retain the MM formulation itself may be significant by modifying the MM 'force fields' formalism for metal ions. The angular deformations in coordination sphere are described with potential functions more sophisticated than harmonic ones. In this framework, the so-called 'points-on-a-sphere' (POS) scheme was proposed [,,,]. It suggests the shape of coordination polyhedron to be ultimately dictated by the inter-ligand van der Waals - like interactions. Recently, this approach has been criticized and improved [,] by considering not the inter-ligand interaction but repulsion of effective interaction centers placed at the coordination bonds, as it is suggested by well-known and quite successful qualitative theories by Gillespie [] and Kepert []. The Gillespie-Kepert MM (MMGK) as described in [,] allows a proper description of many cases of significant distortion in coordination geometry.

Within the MMGK method the total conformation energy of a molecule is:

EMM = å Eb+ å Eang+ å Etors+ å Enb+ å Eimp+ å Erep,

(16)

where the energy terms are

Eb = 0.5Kr(r-r0)2

- the energy of bond stretching (except metal - donor atom bonds);

Eang = 0.5Kq (q-q0)2

- the energy of valence angle bending; note that the valence angles involving the metal ion as vertex are not considered as they are described through the Gillespie-Kepert term (see below eq.());

Etors = 0.5V0(1+cos[n(f+y)])

- the energy of torsion interaction;

Enb = eij æ ç è r0 rij ö ÷ ø 12   -2eij æ ç è r0 rij ö ÷ ø 6

-the energy of non-bonded interaction;

Eimp = 0.5Kimpd2

-the energy of improper torsion (out-of-plane) interaction.

Metal-ligand bonds stretching is modeled with Morse potential:

Eb = D0[e-a(r-r0)-1]2.

(17)

in order to adequately reproduce large variations of bond lengths occurring in different TMC spin states.

According to the MMGK model, the arrangement of the donor atoms around the metal is dictated by repulsion of the effective centers lying on the M-L bonds at a distance of r_eff from the metal ion. This term implicitly (partially) accounts for the electronic effects in the coordination sphere which are not covered by the standalone EHCF formalism (which gives only the d-shell energy) since the former are pertinent to the s-p interactions with the lone pairs. The energy of the 'bond repulsion' in the coordination sphere is:

Eijrep = Aij/Rij6,

(18)

where

Rij2 = ri,eff2+rj,eff2-2ri,effrj,effcos( ^ XiMXj   ); ri,eff = R(M-Xi)di,eff;

and R(M-X_i) is the actual M-X_i bond length; A_ij, d_eff,i, d_eff,j are the GK force field parameters, characterizing the repulsion energy (A) and positions of the repulsion centers (d_eff).

The MM parameters for organic part of molecule were primarily taken from the CHARMM force field [,], while metal-specific parameters are fitted within different versions of the EHCF(L)/MM method (see Ref. DarhRazPletTch and the present work).

Parameters of the Morse potential fitted in the present paper (see below for details of the fitting procedure) to be used with the EHCF(L)/MM approximate scheme are given in Table l. At the time being we limit ourselves to fit only the parameter set of the Morse potential for the sp²-hybridized nitrogen (atom type 8 within MMPC force field).

Full set of the MMGK force field parameters is available upon request.

2.4.3  Junction parameters

Since the EHCF(L)/MMGK procedure is a specific case of a general QM/MM approach and the entire complex is divided into two parts, namely the d-shell and the l-system, their interaction requires separate attention (within the standard EHCF model this interaction ultimately results in the d-shell splitting). In the QM/MM context the intersystem interaction is habitually termed as junction. Not like in other hybrid QM/MM schemes the form of the junction in the present EHCF(L)/MM scheme is not taken ad hoc but is given by the EHCF [] and EHCF(L) TchDarh,DarhTchRJPC,DarhRazPletTch theories represented by eqs.( 3)-(5). The precise numerical values of the junction-related quantities are calculated on the basis of the theory reviewed above. An important component of this theory is that some electronic structure underlying the MM-part of the system is assumed. Parameters characterizing this implied electronic structure of the l-system are used in order to estimate the intersystem junction (interaction). These two kinds of parameters corresponding, respectively, to the d-l interaction itself and to the l-system electronic structure are characterized below.

d-l interaction parameters  

In the original EHCF theory the specific parameters describing the interaction between the d- and l-systems were fit to reproduce the d-level splitting for octahedral complexes with a specific donor atom. The set of the intersystem interaction parameters includes the g_sd and [`(g)]_pd parameters of the Coulomb interaction between the d-shell and transition metal valence s- and p-electrons. These parameters are taken from the Oleari's work [], the valence state ionization potentials for the d-shell and the donor atoms are taken from BoehmGleiter, and the dimensionless factors b_ML characteristic for a metal - donor atom pair, scaling the resonance interaction, are transferred from the original EHCF [] to the EHCF(L)/MM without change. The orbital exponents necessary for calculating the overlap integrals employed in parameterizing the resonance integrals are also taken from [] as is.

Electronic structure parameters of the l-system.  

The electronic structure parameters of the l-system required for the calculation of the effective Hamiltonian eq.(1) by eqs.( 3)-(5) are the one-electron densities (effective charges), orbital energies, MO-LCAO expansion coefficients. The original EHCF SouTchMis method employs the CNDO approximation [] in order to estimate these quantities. They are calculated for each molecular geometry by the approximate SCF procedure extended to the entire l-system. The local version EHCF(L) reduces the computation effort by the price of requiring additional parameters: (i) a set of expansion coefficients for each local state (LMO's) related to the LP involved in the complex formation; (ii) the orbital energies e_{L i}⁽⁰⁾ of the ligand MO's having non-zero contribution to the lone pair of the donor atom; (iii) a set of partial densities r_iN^L of the i-th MO's on the atom N for each ligand L eq.(14); and (iv) the effective point charges (Q_A⁰, A Î L ). The expansion coefficients of the LP over MO's c_Li^L , and the coefficients [`(c)]_La, where a runs over the AO's of the donor atom having the dominating contribution to the LP, are calculated within the ligand fixed coordinate frame (LFCF). These calculations are done separately for the free ligand molecules with use of the semiempirical CNDO procedure which is known to be parameterized to decently reproduce the dipole moments of organic molecules in terms of effective atomic charges. Since the charge redistribution between the ligands and the metal s- and p-orbitals is neglected in the RLMO model, the occupancies of the metal AOs P_aa in eq.(3) remain zero.

3  Results and Discussion

In the present study a procedure for treating PES of TMC within a general QM/MM-like framework is applied to a series of complexes. In Ref. DarhRazPletTch it was shown how to reformulate the semi-empirical EHCF theory which previously allowed to calculate with acceptable accuracy the crystal field induced by the ligand's on the TMC's d-shells in the local form (i.e. in terms of the crystal field increments induced by the lone pairs of the ligands) . This gave us explicit formulae for the crystal field matrix expressed through the ESP of the free ligands and a procedure to calculate them as functions of molecular geometry.

3.1  Implementation

In our previous work [] a procedure combining the EHCF(L) QM approach and the MMGK MM procedure by eq.(9) was implemented in the MMECF 1.0 package [] which allows both gradient minimization and minimization without derivatives [,] for the molecular energy. The package also allows to consider either whole ligands or their fragments as rigid bodies. Technically the ligand geometries employed while parameter fitting (see below) within the rigid body scheme were first pre-optimized with use of 'pure' MM, and remained fixed in the further calculations.

3.2  Spin states and geometry of iron (II) complexes

The methodology described above was applied to twenty two Fe²⁺ complexes listed in Table 2 together with relevant Cambridge Crystal Structure Data Bank (CCSDB) reference data and the spins of the ground states. The ligands are shown in Fig.1. The series contains compounds with monodentate and polydentate ligands, both low- and high-spin ground states. As we notice above, there is a parameter set for nitrogen N(sp²) atom only, so we select from the CCSDB database complexes with ligands containing this type of donor atom strictly.

Experimental geometries of the above complexes were taken from the CCSDB. Hydrogen atoms were added where necessary. The complexes XVII-XXII exhibit spin-crossover and crystal structures for both low- and high-spin states are known which allows detailed comparison of results of our calculations with experiment.

Since the original EHCF method provides an accurate description for the crystal field in TMC's itself [], as well as its dependence on fine geometry variations [,], we first fitted the MM parameters describing the Fe-N interactions (it is that part of the MM model which is affected by introducing the EHCF(L) d-shell energy). The calibration has been performed for the set of molecules I-II, VIII-X by trial-and-error procedure with a twofold restriction: (i) to ensure that experimental total spin of ground states is properly predicted, and (ii) to reproduce experimental geometry of this spin isomer.

Within such a setting, the complexes I - II, VIII-IX have been already studied in Ref.[] with use of the standard Racah parameters for free Fe²⁺ (B₀ = 917 cm^-1, C₀ = 4040 cm^-1) [] and within the rigid ligand approximation. The general conclusion was that the shape of the potential energy profiles is in agreement with experimental data, but some important details are wrong. In order to clarify this situation, the 10Dq parameter for octahedral [Fe(Py)₆]²⁺ complex as a function of the metal-nitrogen distance with use of the EHCF(L)(RLMO) and by the original EHCF procedure has been calculated in []. It was observed that for the ''interesting'' range of the Fe-N distances (near 2 Å , where most of Fe(II)-N distances lie) the 10Dq value in EHCF(L)(RLMO) is 1.5-2.2 times larger than in the EHCF which is shown to fairly coincide with experimental values at experimental geometry of the complex []. In other words, the d-shell splitting in the RLMO model is too strong as compared to the experiment (the latter being well reproduced by the EHCF method).

Analysis of this situation with use of the Tanabe-Sugano diagrams performed in Ref.[] has shown that since the EHCF(L)(RLMO) model systematically overestimates the splitting parameter, the correct position of the crossection point of the low- and high-spin terms can be recovered by scaling the free ion B₀ and C₀ values by a factor h. As it was shown in Ref.[], the best agreement between the calculated and experimental geometries in the EHCF(L)(RLMO) model is reached at h = 1.5 (B and C Racah parameters of 1430 and 6040 cm^-1, respectively).

In the present work parameterization of the MM part of the method was performed for the specified value of h for complexes I-II, VIII-X with rigid ligands using gradient minimization. The parameters obtained for the EHCF(L)(RLMO) model (Table 1) slightly differ from those of Ref. DarhRazPletTch due to different minimization procedure (minimization without derivatives) used in that work.

Below we consider results of our calculations with use of the parameters given in Table 3. All the complexes have been studied by the rigid ligands optimization using analytic gradients. The geometry was optimized starting from the experimental structures until the root-mean-squared (RMS) energy gradient was less than 0.1 kcal· mol^-1·\.Å ^-1. Calculated geometries of this series of the complexes, in general, agree rather well with the experimental data.

In the low-spin (S=0) complexes I-VII maximal RMS difference between the calculated and experimental metal-nitrogen bond lengths is 0.04 Å and maximal difference for the bond angles in coordination sphere is 8^o. The low-spin ground states of the complexes are reproduced in all cases besides the last complex VII where the energy difference between low- and high-spin forms is very small.

We performed full optimization of ligands geometry for the singlet complexes I and II and found certain improvement of results as compared to rigid ligands optimization. Notably, for the complex I the difference between the Fe-N bond lengths of central and terminal pyridine rings (which is inherent feature of this structure attainable to mis-direction of terminal nitrogen lone pairs) obtained in the calculation (0.11 Å ) only slightly differs from the experimental value (0.12 Å ). The Fe-N bond lengths themselves are in good agreement with experiment with the difference about 0.03-0.04 Å .

A numerical experiment was performed on the low-spin complex V. Corresponding ligand, bpmy, may show two binding modes which differ by the nitrogen atoms used for coordination (marked with * or ** in Fig.1). Here, one can test whether the proposed theory correctly predicts an actual binding mode. The results presented in Table 4 show that the ground state total energy for the (*) isomer is lower than that for the (**) isomer (notice the same ground state spins in two complexation modes). Thus, the (**) complexation mode is thermodynamically far less preferable; this fairly agrees with experiment [] where only the low-spin (*) mode is detected.

Experimentally, strongly distorted octahedral [Fe(Py)₆]²⁺ complex VIII has the high-spin ground state. Our EHCF(L)(RLMO)/MMGK calculation correctly reproduces this fact. However, the Fe-N bond lengths are significantly underestimated. This also happens for the majority of high-spin complexes of the series; such a behavior perhaps comes from too rigid relation between potential depth, elasticity and anharmonism of the used Morse form of the potential (which does not allow to fit these parameters independently). At the same time, calculated valence angles at the metal ion in VIII-XV are in good agreement with experiment which indicates correct description of ligand-ligand repulsion by the MMGK model. An important feature is that the torsional angles between the pyridine rings in the calculated structure VIII differ significantly from the experimental values. With this in mind one may speculate with good confidence on strong nonbonded interactions in the crystal environment which are responsible for observed discrepancies. Indeed, the close contacts in experimental structure of VIII lead to highly asymmetric environment for the cation VIII. The negatively charged Fe(CO)₆^4- ion in this structure is close to three pyridines which are more distant from the Fe²⁺ ion than other three, while nothing should break the symmetry for the isolated structure of VIII.

Calculation on tris-(2-methyl-bipyridine) complex IX also yields the high-spin ground state, in agreement with experiment. The calculated lengths of the Fe-N bonds for the methylated ring are shorter than the experimental ones which is a common feature for the high-spin complexes.

The complex X with six 1-methyl-imidazole ligands is another example of successful reproducing high-spin ground state. However, the Fe-N bond lengths are also underestimated here. Manifestation of the Jahn-Teller effect in the experimental geometry where two axial bonds are by 0.03 Å  longer than four equatorial ones is not reproduced in the calculation. Very probable reason for the latter drawback is the absence of the p-overlap between the ligand p-system and the t_2g orbitals of the metal, which is the physical source of the Jahn-Teller distortion in the Fe(II) complexes in the present model. Here only the lone pairs of the ligands are taken into account in eq.(5). It was pointed out in paper DarhTchRJPC that taking into account some p-orbitals in eq.(5) gives considerable values of the e_p parameters. It is still to be tested whether taking into account localized p-orbitals on the donor atoms contributes to improve the above result.

In the complex XVI the spin-crossover is observed though the structure of the low-spin form is only known. Indeed, the high-spin ground state in this case has a lower energy; however, the calculated difference between the singlet and quintet states is only about 2.5 kcal/mol which can indicate a potential existence of spin-crossover. In the work [] the complex XVI was calculated by quasi-relativistic INDO method which also get optimal geometry featuring the Fe-N distances about 2.02 and 2.13 Å . These values may be characterized as indicating to the high-spin state. That small difference indicates also that the interaction with counterions and other components of the crystal structure may be important to estimate the true ground state of a TMC in the crystal (see below).

Both low- and high-spin calculated geometries of spin-crossover complex XVII-XVIII [] are is in general agreement with the experiment, but there are large difference between the calculated (2.23 Å ) and experimental (2.09 Å ) Fe-N distances for the central pyridine nitrogen. Important structural details nevertheless reproduced in our calculation are (i) the slight (7^° ) twist between the terminal (with the phenyl substituents) and central pyridine rings and (ii) stacking of the phenyl substituent from one tpy with the central pyridine ring belonging to another tpy ligand.

The calculation results require a special caution while compared with experimental data if the complexes manifest both spin isomers. In general such a situation is only possible if the energy difference between two spin forms is small. However, the last condition refers not to the relative energies of the isolated isomer molecules but to those in the appropriate environment (crystal or liquid one). The renormalization of the relative isomer energy incurred by interaction of the spin active moiety with its surrounding must be considered separately.

The complexes XVII-XVIII provide an interesting example of the influence of crystal package upon the structure and the spin state of the complex. The counterion in the high-spin isomer XVII is PF₆^- while for the low-spin isomer XVIII it is ClO₄^-. It is concluded from the studies on organic Bechgaard salts that internal (''chemical'') pressure in the crystals of (TMTSF)₂X is higher if X= ClO₄^- while X=PF₆^- [] corresponds to a lower internal pressure. The larger pressure according to the LeChatelier principle favours the smaller unit cell volumes, smaller metal-nitrogen distances, and thus makes the low-spin form more stable SpinTran1.

The influence of counterion on relative stability of the low- (XVII) and high-spin (XVIII) forms of the ionic complex can be inerpreted in terms of the following simple thermodynamic model.

Let us consider a crystal containing both low- and high-spin TMC ions and their counterion (anion) both modeled as charged spheres of different radii interacting via Coulomb attraction and exchange repulsion forces. The free energy dG of a spin-crossover crystal (per unit cell) can be written analogously to the expression taken from []:

dG = (dh-Tds)x+kT(xlnx+(1-x)ln(1-x))- ae2 r + c (r-R)m

(19)

where dh and ds are intramolecular spin-crossover enthalpy and entropy change, respectively, T is the temperature, x is the fraction of the high-spin form (0 £ x £ 1), a is the Madelung constant which is of order of unity, c is a parameter of the intermolecular interaction potential, and r is the distance between anion and cation charge centers. The moieties forming the crystal are in fact bulk organic molecules and thus the exchange repulsion between the component molecules in fact depends on the distance y between the surfaces of the spheres representing them in the model, rather than on the distance between their centers. The Coulomb attraction on the other hand depends namely on the distance r between the centers. These two quantities are related by the condition y = r-R where R is the sum of the radii of the spheres representing the cation and the anion moieties forming the crystal. The transition from a low to high-spin form in a cation from the point of view of such a crystal is a spontaneous increase of the radius of its modeling sphere by a small quantity dr_HL. For the crystal with the x fraction of the high-spin cations we have on the average:

R = R0+xdrHL = rA+rL+xdrHL

where r_A and r_L are the effective radii of the anion and the low spin form of the cation, respectively. Inserting we get for the last two terms of eq.(19) as contribution from potential energy of anion-cation interaction:

dE = - ae2 y+R0+xdrHL + c ym

Then, taking into account that y << R₀ as well as dr_HL << R₀, we can simplify it so that:

dE = - ae2 R0 (1- 1 R0 (y+xdrHL))+ c ym

and the enthalpy dh is effectively renormalized due to interaction by the positive quantitity:

dh¢ = dh+ ae2 R02 drHL

Thus, the effective relative energy of the high-spin isomer is lower for larger r_A (and thus R₀). In other words, the larger anions (like PF₆^-) make high-spin form be easier accessible in the crystal. It is in accordance with experimental facts cited above.

In any case one should expect that going from the vacuum to the condensed medium would shift the preference towards low-spin form due to internal pressure. In all the examples where the proposed method predicts the ground spin state erroneously as compared to the experiment it is always the high-spin state which is predicted so that it possibly can be corrected by taking into account the effects of the crystalline environment for which actual measurements are performed.

In the complex XXI (high-spin form) a good agreement was obtained in metal-nitrogen distances despite the large (comparing to other low-spin molecules) difference of these distances in the corresponding low-spin form XXII. It is understandable as in XXI distances are not so long (as for other high-spin molecules) and excessive rigidity of the Morse potential does not influence the minimum position too much.

Despite the difference in the metal-donor atom bond distances one can see acceptable accuracy of the hybrid EHCF(L)(RLMO)/MMGK scheme. Notably, the differences of calculated energies of the spin isomers, under suggested parameterization (and with the above stated precautions), are large enough to explain the fact that only one ground spin state is experimentally observed for those considered compounds where it actually happens. The hybrid EHCF(L)(RLMO)/MMGK method correctly predicts the relative stability of different isomers of TMCs as well as the influence of the ligand substituents on the spin of the ground state as, for example, in the case of complexes II (low-spin with bpy) and X (high-spin with methylated bpy).

General comparison of the RMS deviations between the experimental and calculated distances for a series of complexes shows that for the low- and high-spin complexes considered separately they are respectively 0.03 and 0.05 Å . In Refs.[,] authors give better values of the RMS deviations for the Ni(II) complexes (0.01-0.02 Å ). Nevertheless, our calculations show reasonable correspondence with experiment not only in geometries but also in the energy differences between the low- and high-spin forms. Values of the RMS Fe-N distances obtained here can be interpreted in the light of work [] where distribution of these distances in the six-coordinated Fe(II) complexes was analyzed using the facilities of the CSD database. Two major maxima were found in that distribution: one about 1.9 Å  and another 2.2 Å  which can be identified with the distances characteristic for the low- and high-spin complexes. Our calculations generally reproduce this picture rather well.

In the work [] a new force field for TMC is proposed and good accuracy in correspondence between experimental and calculated structures is claimed. The data set used for parameterization in Ref.[] is, however, limited to low-spin organometallic compounds, at least for Fe(II) complexes. Other examples of calculations on different Fe(II) spin states in complexes with bulk organic ligands present in the literature (see Refs. [,]) are performed by pure MM methods with different parameters for the low- and high-spin iron ions.

4  Conclusion

In the present paper we develop the EHCF(L) methodology [] and the general approach to development of the hybrid QM/MM methods TchPCCP in application to the problem of transition metal incorporation into MM. On this route we reconsidered the EHCF(L)(RLMO)/MMGK approach which resulted in a unified parameterization for different spin states of Fe(II) complexes with nitrogen ligands. On the basis of the performed analysis it can be stated that the concert usage of the EHCF(L)(RLMO) procedure as a QM component for describing the geometry dependence of the d-shell energy together with the MMGK procedure as the MM component for describing the ligand energy a unified MM-like description for the PES of different spin states of iron(II) complexes is achieved. Remarkably it uses the single spin-independent parameterization of the MM-part of the system.

5  Acknowledgments

The work is performed with partial financial support through the grant No 6-120 dispatched to ALT by the Young Researchers Commission of RAS. The usage of the CCSDB is supported by the RFBR grant No 99-07-90133. 

References

[]

Burkert, U., Allinger, T. Molecular Mechanics. ACS, Washington, D.C. 1982.

[]

Darhovskii; M.B., Tchougréeff, A.L. Chem. Phys. Reps. 1999, 18, 73.

[]

König, E.; Ritter, G.; Kulshreshtha, S.K. Chem. Rev. 1985, 85, 219.

[]

Gütlich, P. Struct. Bond. 1981, 44, 83.

[]

Tchougréeff, A.L. Phys. Chem. Chem. Phys. 1999, 1, 1051.

[]

Soudackov, A.V.; Tchougréeff, A.L.; Misurkin, I.A. Theor. Chim. Acta. 1992, 83, 389.

[]

C.A. Coulson. Valence. Oxford Univ. Press. 1961.

[]

Frenking, G.; Frölich, N. Chem. Rev. 2000, 100, 717.

[]

Darhovskii, M.B.; Razumov, M.G.; Pletnev, I.V.; Tchougréeff, A.L. Int. J. Quant. Chem. 2002, 88, 588.

[]

Darhovskii, M.B.; Tchougréeff, A.L. Russ. J. Phys. Chem. 2000, 75, 112.

[]

Lever, A.B.P. Inorganic electronic spectroscopy. Amsterdam, Elsevier Science, 1984.

[]

Abrikosov, A.A.; Gor'kov, L.P.; Dzyaloshinskii, I.E. Methods of quantum field theory in statistical physics. Dobrosvet, Moscow, 1998 [in Russian], 2nd edition.

[]

http://qcc.ru/~ netlab: Entry ECFMM.

[]

Pletnev, I. V.; Melnikov, V. L. Koord. Khim. 1997, 23, 205 [in Russian].

[]

Razumov, M. G.; Melnikov, V. L.; Pletnev, I. V. J. Comput. Chem. 2001, 22, 38.

[]

Perkins, P.G.; Stewart, J.J.P. J. Chem. Soc. Faraday Trans.(II). 1982, 78, 285.

[]

Hambley, T. W.; Hawkins, C. J.; Palmer, J. A.; Snow, M. R. Aust. J. Chem. 1981, 34, 45.

[]

Reichert, D. E.; Hancock, R. D.; Welch, M. J. Inorg. Chem. 1996, 35 (24), 7013-7020.

[]

Cundari, T. R.; Fu, W.; Moody, E. W.; Slavin, L. L.; Snyder, L. A.; Sommerer, S. O. J. Phys. Chem. 1996, 100 (46), 18057-18064.

[]

Johnson, D. W.; Brubaker, G. R. Inorg. Chim. Acta. 1986, 119, 131.

[]

Gillespie, R. J.; Hargittai, I. The VSEPR Model of Molecular Geometry; Allyn and Bacon: 1991.

[]

Kepert, D. L. Inorganic Stereochemistry; Springer: Berlin, 1982.

[]

Tchougréeff, A.L. Chem. Phys. Reps. 1997, 16, 62.

[]

Tchougréeff, A.L. The effective crystal field methodology as used to incorporate transition metals into molecular mechanics, In Molecular modeling and dynamics of bioinorganic systems (eds. Banci L., Comba P.): NATO ASI workshop, Kluwer, 1997. P.217.

[]

Burton, V.; Deeth, R.; Kemp, C.; Gilbert, P. J. Am. Chem. Soc. 1995, 8407.

[]

Deeth, R.; Paget, V. J. Chem. Soc. Dalton Trans. 1997, 537.

[]

Gerloch, M.; Wooley, R.G. Struct. Bond. 1983, 371.

[]

Tchougréeff, A.L. Chem. Phys. Reps. 1998, 17 (6), 62.

[]

Tchougréeff, A.L.; Soudackov, A.V.; Misurkin, I.A.; Bolvin, H.; Kahn, O. Chem. Phys. 1995, 193, 19.

[]

Pople, J.A.; Beveridge, D.L. Approximate molecular orbital theory, McGraw-Hill, 1970.

[]

Shaeffer, C.E. Struct. Bond. 1968, 5, 68.

[]

Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; Dates, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187.

[]

Smith, J. C.; Karplus, M. J. Am. Chem. Soc. 1992, 114, 801.

[]

Di Sipio, L.; Tondello, E.; De Michelis, G.; Oleari, L. Chem. Phys. Lett. 1971, 11, 287.

[]

Böhm, M.C.; Gleiter, R. Theor. Chim. Acta 1981, 59, 127.

[]

Nelder, J. A.; Mead, R. Computer J. 1965, 7, 308.

[]

Brent, R. P. Algorithms for minimization without derivatives. Prentice Hall, Englewood Cliffs, 1973.

[]

Laine, P.; Gourdon, A.; Launay, J.-P. Inorg. Chem., 1995, 34, 5156.

[]

Dick, S. Z. Kristallogr. - New Crystal Structures, 1998, 213, 356.

[]

Kucharski, E.S.; McWhinnie, W.R.; White, A.H. Aust. J. Chem., 1978, 31, 53.

[]

De Munno, G.; Poerio, T.; Viau, G.; Julve, M.; Lloret, F. Angew. Chem., Int. Ed. Engl., 1997, 36, 1459.

[]

Breu, J.; Range, K.-J.; Herdtweck, E. Monatsh. Chem., 1994, 125, 119.

[]

Boca, R.; Baran, P.; Dlhan, L.; Fuess, H.; Haase, W.; Renz, F.; Linert, W.; Svoboda, I.; Werner, R. Inorg. Chim. Acta, 1997, 260, 129.

[]

Onggo, D.; Scudder, M.L.; Craig, D.C.; Goodwin, H.A. Aust. J. Chem., 2000, 53, 153.

[]

Childs, B.J.; Cadogan, J.M.; Craig, D.C.; Scudder, M.L.; Goodwin, H.A. Aust. J. Chem., 1997, 50, 129.

[]

Doedens, R.J.; Dahl, L.F. J. Am. Chem. Soc., 1966, 88, 4847.

[]

Onggo, D.; Hook, J.M.; Rae, A.D.; Goodwin, H.A. Inorg. Chim. Acta, 1990, 173, 19.

[]

Seel, F.; Lehnert, R.; Bill, E.; Trautwein A. Z. Naturforsch., Teil B, 1980, 35, 631.

[]

Hibbs, W.; Arif, A.M.; van Koningsbruggen, P.J.; Miller, J.S. CrystEngComm., 4,1999.

[]

Boinnard, D.; Cassoux, P.; Petrouleas, V.; Savariault, J.-M.; Tuchagues, J.-P. Inorg. Chem., 1990, 29, 4114.

[]

Lorente, M.A.M.; Dahan, F.; Sanakis, Y.; Petrouleas, V.; Bousseksou, A.; Tuchagues, J.-P. Inorg. Chem., 1995, 34, 5346.

[]

Sugiyarto, K.H.; Craig, D.C.; Rae, A.D.; Goodwin, H.A. Aust. J. Chem., 1995, 48, 35.

[]

Reger, D.L.; Little, C.A.; Rheingold, A.L.; Lam, M.; Concolino, T.; Mohan, A.; Long, G.J. Inorg.Chem., 2000, 39, 4674.

[]

Constable, E.C.; Baum, G.; Bill, E.; Dyson, R.; van Eldik, R.; Fenske, D.; Kaderli, S.; Morris, D.; Neubrand, A.; Neuburger, M.; Smith, D.R.; Wieghardt, K.; Zehnder, M.; Zuberbuhler, A.D. Chemistry-A European Journal, 1999, 5, 498.

[]

Holland, J.M.; McAllister, J.A.; Zhibao Lu; Kilner, C.A.; Thornton-Pett, M.; Halcrow, M.A. Chem.Comm., 2001, 577.

[]

Stassen, A.F.; de Vos, M.; van Koningsbruggen, P.J.; Renz, F.; Ensling, J.; Kooijman, H.; Spek, A.L.; Haasnoot, J.G.; Gutlich, P.; Reedijk, J. Eur. J. Inorg. Chem., 2000, 2231.

[]

Linert, W.; Konecny, M.; Renz, F. J. Chem. Soc. Dalton Trans., 1994, 1523.

[]

Jerome, D.; Schulz, H.J. Adv. Phys., 1982, 31, 443.

[]

Tchougréeff, A.L.; Darhovskii, M.B. Int. J. Quant. Chem., 1996, 57, 903.

[]

Simaan, J.; Poussereau, A.; Blondin, G.; Girerd, J.-J.; Defaye, D.; Philouze, C.; Guilhem, G.; Tchertanov, L. Inorg. Chim. Acta, 2000, 299, 221.

[]

Gajewski, J.J.; Gilbert, K.E.; Kreek, T.W. J. Comput. Chem., 1998, 19, 1167.

[]

Hancock, R.D.; Weaving, J.S.; Marques, H.M. J. Chem. Soc. Chem. Commun., 1989, 1176.

[]

Marques, H.M.; Munro, O.Q.; Grimmer, N.E.; Levendis, D.C.; Marsicano, F.; Pattrick, G.; Markoulides, T. J. Chem. Soc. Faraday Trans., 1995, 91, 1741.