Economical computational tools suitable for estimations of electronic structure and molecular geometry of transition metal complexes (TMC) are highly in demand. The molecular mechanics (MM)  both itself and in molecular dynamics framework is intensely used in calculations of proteins and other polyatomic organic molecules. During the last decade several attempts were made to apply the conventional MM scheme to the metal ion complexes with organic ligands Hay,Rappe,Landis95,Landis98,Comba,CombaCCR,Lehmann,Sabolovic,Hambley,Zimmer. The main problem here is that in TMC several electronic states may occur in a narrow energy range close to its ground state. Sometimes, the potential energy surfaces (PES) corresponding to different electronic terms of the metal ion d-shell intersect which results in spin transitions SpinTran. In organic molecules this problem normally does not appear and MM description is valid since the electronic excited states are well separated on the energy scale from the respective ground state. In these cases a single quantum state of electronic system suffice for description of a molecule. Clearly, this is not mandatory true for TMCs.
Also, within MM it is hard to get an adequate modeling of coordination sphere, in particular, to account for the flexibility of coordination polyhedron. The most straightforward way is to describe deformations of valence angles in coordination sphere (angles with metal atom as their center) with potential functions more sophisticated than harmonic potentials. Also a so-called 'points-on-a-sphere' (POS) approach was proposed [,]. It suggests the shape of coordination polyhedron to be ultimately dictated by the inter-ligand van der Waals-like interactions. Recently, it has been shown [,] that it may be further improved by considering not the inter-ligand interaction (described through common nonbonding 6-12 or 6-exp potentials) but repulsion of effective interaction centers placed somewhere on the coordination bonds, as it is suggested by well-known and quite successfull qualitative theories by Gillespie  and Kepert . This approach recently brought insight of coordination geometries diversity. It allows a proper description of many cases of significant distortion in coordination geometry (for discussion and examples see [,]). However, being an MM method it is unable (and obviously not designed) to describe spin states of coordination compound, which is necessary to discuss magnetic properties, as well as to provide correct estimates for energetics of a large number of important chemical and biochemical processes where coordination compounds take part.
In order to incorporate electronic effects of the partially filled d-shell in TMC's it was proposed in [,,,,] to include the energy of the d-shell as a separate contribution. It is done in variance with works [,,] where the accent is put on estimating the spectral characteristics of the d-shell at the geometry assessed with use of an MM treatment. Including the energy of the d -shell explicitly allows to account for electronic structure influence on the geometry of TMC. These are quantum effects specific for the open d-shell which appear due to possible degeneracies of different electronic terms of the latter at certain complex geometries. Experimentally this would correspond to the Jahn-Teller complexes and to spin active complexes. However, the ligand field energy in Refs. Deeth1,Deeth2,CombaCCR,CombaMMAOM1,CombaMMAOM2 depends only on the distance between metal ion and ligand donor atoms, which seems to be an oversimplified picture since the effects of lone pair orientation on the ligand field must be taking into account.
Promising methods for quantitative estimates on TMCs are proposed in the framework of the hybrid quantum mechanics/molecular mechanics (QM/MM) methodology [,,]. Most of their applications belong to organometallic realm or to that of complexes of heavy (second and third transition row) metals, without addressing different spin states within a calculation. In the frame of these attempts rather large part of the TMC is treated by QM method leaving to MM only the periphery of the molecule.
In the works [,] we proposed and tested a hybrid QM/MM description of TMCs targeted at first transition row metal complexes. It is based on the general approach of  to description of molecular electronic structure and potential energy, which makes it possible to apply the QM description to that part of molecule in which electronic terms are close on the energy scale, and to use the MM description for the part where the electronic states remain distant in energy. In Refs. TchIJQC,DarhPletTchJCC a combination of the local version  of the effective Hamiltonian of crystal field (EHCF) method , EHCF(L), with the MMGK procedure  has been proposed and implemented for a series of Fe(II) complexes with nitrogen containing ligands. Within the EHCF(L) approach the effective crystal field is presented as a sum of lone pair contributions (see below) which allows for detailed description of the ligand field dependence not only on metal-ligand distance but also on lone pair orientation with respect to metal ion. In the present paper, we report further improvement of the EHCF(L) approach which allows the application of the enhanced hybrid method to molecular and electronic structure estimates for the series of Fe(II) and Co(II) complexes of low- and high-spin.
The paper is organized as follows: in the next Section we briefly review the basic features of the EHCF(L) method [,]. Next we describe an improved EHCF(L) approach taking into account the ligand polarization in TMC. The last Section provides the parameterization and application of EHCF(L)/MMGK approach to calculations of several complexes.
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 composition and structure. In this section we review working approximations based on the EHCF(L) theory  and propose the improved EHCF(L) method taking into account the ligand polarization.
The concept of separating electron variables is to be employed when a hybrid
QM/MM method is developed : electrons have to be divided into
groups, some of the groups whose excited electronic states are accessible in
the experiment are treated by a QM method whereas the behavior of other
groups whose excited electronic states lay high in energy (and are not
accessible in experiment) 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
SouTchMis it is shown that the d-shell can be described by the effective
QM Hamiltonian Hdeff:
In our paper  we have derived and tested a local version of
the EHCF method EHCF(L). It was shown that the splitting parameter 10Dq
can be estimated with the error not exceeding 0.1 eV (this accuracy compares
to that of the EHCF method itself) by the formula:
The resonance integrals bmL in eq. (6) can be expressed through the tL vector of 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 L :
Then, introducing the quantities:
The expression eq. (9) defines the ell¢L parameters in terms of the quantities which can be calculated within the EHCF(L) method. Their relation with the standard angular overlap model (AOM)  is described in details in Ref. . There eqs. (9), (10) have been used to calculate the values of the es and ep parameters for a series of octahedral complexes with nitrogen containing ligands. That calculation was in a good agreement with experimental 10Dq values (within 10% accuracy).
In the previous Section we reviewed the EHCF(L) theory which allows to estimate the crystal field in terms of local electronic structure parameters (ESP) of the ligands. By this method it can be done for arbitrary geometry of the complex, which is prerequisite for developing a hybrid QM/MM method.
The natural way to go further with this technique is to apply the
perturbation theory to obtain estimates of the l-system Green's function
entering eqs.(6) and/or (9). We assume for the purposes
of the present work that the bare Green's function for the l-system has a
The Coulomb interaction between the ligands themselves and between each of them and the metal ion when turned on does not break the block diagonal structure of the bare Green's function G00l. Then the approximate Green's function for the l-system conserves the form eq. ( 11) but with the poles of Green's function matrices G0L corresponding to the orbital energies of the ligand molecules in the Coulomb field induced by the central ion and by other ligands (L¢ ¹ L ) rather than to those of the free ligands. This estimate for the ligand orbital energies is sufficiently accurate to enter eq. (7).
In the following Subsection we consider an implementation of this approach taking into account the Coulomb field effects and thus allowing to express the Green's function of the l-system in terms of the Green's functions of separate ligands.
The simplest picture the influence of the central ion on the surrounding
ligands reduces to that of the Coulomb field affecting the positions of the
poles of the Green's function (orbital energies) of the free ligand. The
form of the CMO's of each ligand is left unchanged which is a picture of the
rigid ligands' MO's (RLMO) Ref. . According to
GreenMonster, the effect of the Coulomb field upon the orbital energies is
These formulae comprise the RLMO model of the electronic structure of the l-system of the TMC. The RLMO procedure has been implemented in the program suite ECFMM 1.0 . Its application to analysis of molecular geometries are described in Refs. [,]. Despite satisfactory results in geometry and spin states description of some iron(II) complexes, the main conclusion is that the crystal field is reproduced with too large error due to overestimated repulsion of d-electrons from the ligands lone pairs. Due to that in Refs. TchIJQC,DarhPletTchJCC we had to change the Racah parameters differently for complexes with pyridine-like and amino nitrogen donor atoms, where in the latter corresponding error was small. In the following Subsection we propose an improved model taking into account polarization effects in the ligand sphere which is able to describe both types of ligands within unified parameterization.
The model of electronic structure of TMC which considers the metal ion as a point charge equal to its oxidation degree or formal charge, is habitually called the sparkle model . Within models of that type semiempirical SCF calculation is performed for the ligands of the complex placed in the electrostatic field induced by the central ion with its formal charge ('sparkle'). Thus, the charge redistribution occurring in the ligands is obtained by performing a standard SCF procedure for them.
Within models of the sparkle family the effect of the external Coulomb field does not reduce to the renormalization of the orbital energies as it is within the RLMO model (see above). The electron distribution also changes when the ligand molecules are put into the field. That means that the density matrix of the system varies and, accordingly, the effective point charges QN residing on the atoms of the molecule in eq. (15) change. We will describe this situation by means of polarizability concept . According to it the difference between polarized and non-polarized effective charge on atom A is:
Though procedures of that sort are admitted in modern MM schemes directed to the systems with significant charge redistribution  we consider such a procedure to be too resource consuming and restrict ourselves by several lower orders with respect to P. Then the term Pdh0 corresponds to the first order perturbation by the Coulomb field induced by the metal ion and bare (non-polarized) ligand charges. The second order term corresponds to the perturbation due to the Coulomb field induced by the mutually (first order) polarized charges:
The atom-atom mutual polarizability matrix P has a block-diagonal
In order to evaluate PL we consider the mutual atomic
orbital polarizabilities PabL (0):
Finally, according to the general formulae given, say, in , the matrix element of the bare orbital mutual polarizability entering eq. (21) is given by:
This is the method for construction the renormalized polarizability matrix PL for the ligand L . Form of the total matrix P for the whole TMC is given by eq. (20). With use of this matrix, we can obtain renormalized atomic charges by the eq. ( 18).
In this Subsection we formulated the perturbative version of the Sparkle approximation for the Green's function G0l of the l-system. It satisfies the requirement imposed above that the Green's function of the l-system must be expressed in terms of those of the free ligands. Then, eqs. (16)-(19) comprise the perturbative form of the Sparkle model of the l-system's electronic structure. As we show in the results Section below, it yields effective atomic charges of sufficient precision using the point charges of the free ligand as a zero approximation. The charges thus obtained are used for calculation of the Sii(f) term according to eq.(14) and for renormalizing the orbital energies by eq.(13). The proposed procedure improves the junction between the EHCF(L) method playing role of the QM procedure and the MM part, as shown below, where details of the calculations performed within this approximation are given (Section 4).
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:
In order to obtain the effective Hamiltonian for the d-shell used in eq.(26) the electronic structure parameters (ESP's) of the l-system must be used in eqs. (1) - (4), (9), (10). These ESP's are condensed in the l-system Green's function. In the previous Section we presented general formulae which comprise the perturbation approach to evaluation of the Green's function of the l-system using those of the separate free ligands as a zero approximation. Estimates of the l-system Green's function following the prescriptions of the above Section can be performed for arbitrary molecular geometry. Inserting this approximate form of the l-system Green's function into the EHCF(L) formulae eqs. (9), (10) 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.
The RLMO and Sparkle models represent the Green's function G0l eqs. (11) - (12) including only the ligand MO energy shifts by Sii(f) eq.(14) calculated with use of the effective atomic charges. The latter are the charges for either free ligands or those polarized by metal ion and other ligands for the RLMO or Sparkle models, respectively. That all comprises the two versions of the hybrid EHCF(L)/MM approach to evaluation of the PES of TMC's. The RLMO approach was thoroughly investigated in [,] so in the present paper we focused on the perturbative Sparkle version of the model.
The EHCF(L)/MMGK method described above in general terms is a specific case of a general hybrid scheme involving QM and MM components which both require extensive 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.
The d-shell parameters are taken from the original EHCF method Ref. SouTchMis without changes. These are the specific exponents of atomic d-orbitals, d-electron core attraction parameter Udd for each metal atom. The Coulomb repulsion of d-electrons is characterized by three parameters: gdd, and the Racah parameters B and C. In the general theoretical setting of the EHCF method the Racah parameters must be taken standard for the free ions as tabulated, say, in Ref. . Pragmatically, however, the values specific for the complex are used in order to reach better agreement between theoretical and experimental spectra [,,,]. In the context of the present study directed towards uniform description of a wide range of complexes with many different ligands only the single values of the Racah parameters common for all complexes of a given metal ion make sense. For the complexes of Fe(II) the Racah parameters B and C for the free Fe2+ cation are used like in [,,,]. In Souexp1,Souexp2,Souexp3,Souexp4 where the EHCF method has been employed for electronic spectra calculations of some Co2+ complexes, various values of the Racah parameters have been used. For example, in case of the complex CoCl64- these values were B=780 cm-1, C=3432 cm-1 , whereas for the complexes Co(H2O)62+ and Co(NH3)62+ values of B=850 cm-1, C=3935 cm-1 and B=885 cm-1, C=4099 cm-1 were used respectively . The free Co2+ ion Racah parameters are B0=971 cm-1, C0=4366 cm-1 . For the considered set of Co(II) complexes, the specific single set of these parameters for Co2+ ion is used both for the low- and high-spin complexes with different coordination numbers and geometries which are somewhat reduced as compared to the free ion values.
Organic part of a molecule and metal ion coordination sphere (leaving out effects of the d-shell) is described in the present hybrid procedure in terms of the MMGK method [,].
Within it the total conformation energy of a molecule is:
- the energy of bond stretching (except metal - donor atom bonds);
Bonding interaction of metal valence 4s- and 4p-subshells with ligands
is currently modeled within the MM part of the combined EHCF(L)/MM scheme
through the Morse potential:
The arrangement of the donor atoms around the metal is dictated by mutual
repulsion between the effective centers lying on the M-L bonds on the
distance reff from the metal ion. This term implicitly partially
accounts for the electronic effects in the coordination sphere which could
not be described within the standalone EHCF formalism (which gives only the d-shell energy). The energy of the 'bond repulsion' in the coordination
sphere is then:
The MM parameters for organic part of molecule were primarily taken from the CHARMM force field [,], while specific metal-dependent parameters must be fitted within different versions of the EHCF(L)/MM method separately (compare Ref.  and the present work). We note once again that the single parameters set is used for any spin state of the TMI under consideration.
Since the EHCF(L)/MMGK approach is a specific case of a general QM/MM scheme, where the entire system 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 interacton is habitually termed as a 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) DarhTchAOM,TchIJQC theories. 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 certain type of 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. These two kinds of parameters corresponding, respectively, to the d-l interaction itself and to the l-system electronic structure (ESP's) 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 gsd 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 bML characteristic for a metal - donor atom pair, scaling the resonance interaction. These parameters are transferred from the original EHCF  to the EHCF(L)/MM without change. The orbital exponents necessary for calculating the overlap integrals employed throughout the parameterizing the resonance integrals are also taken from  as they are there.
ESP's of the l-system. RLMO model
The ESPs of the l-system required for the calculation of the effective Hamiltonian eq.(1) are the one-electron densities (effective charges), orbital energies, MO LCAO expansion coefficients. The original EHCF  method employs the CNDO approximation  in order to estimate these quantities. They are calculated for arbitrary molecular geometry by the approximate SCF procedure extended to the entire l -system. The local version of the EHCF employed in the present work additionally requires a set of expansion coefficients for each local state (LMO's) related to the LP involved in the complex formation, as parameters. The expansion coefficients of the LP [`(c)]La, where a runs over the donor atom AO's having the dominating contribution to the LP, and calculated within the ligand fixed coordinate frame (LFCF) are treated as ESPs of the l-system as well. These parameters cLi, [`(c)]Ls, and [`(c)]Lp proposed in  are calculated separately for the free ligand molecules and are fed to the EHCF(L)/MM procedure as parameters.
In the RLMO approximation Ref.  the orbital energies are
estimated perturbatively which is more economical from the computational
point of view but requires a larger number of parameters. Within the RLMO
model the electronic structure of the free ligand prototype is supposed to
be unchanged during the complex formation. Thus, for the EHCF(L)
calculations we use the charge distribution calculated for the free ligand
itself i.e. the effective point charges which are found from the CNDO
calculation on the free ligand (QA0, A Î L ) and consider them
as ESP's for the l-system. Also the orbital energies of the ligand MO's
having non-zero contribution to the LP of the donor atom calculated for the
free ligand are to be fed to the EHCF(L)/MM procedure. They are used to
estimate the positions of the poles of the Green's function in the Coulomb
field of the charges within the complex according to the formula eq. (
13) with use of the partial densities of the i-th MO's on the atom A
of the ligand L :
ESP's of the l-system. Perturbative Sparkle model
The Sparkle model of the electronic structure of the l-system includes ESPs which are the same as in the RLMO model. The difference with the latter is that the effective charges used in the Sparkle model in eqs. (17), (18), and (19) are renormalized due to polarization by the Coulomb interaction between the ligands themselves and with the metal ion. By eq. (20) the matrix of mutual atomic polarizabilities for the whole l-system is constructed from the polarizability matrices of the free ligands eq. (24) calculated preliminary. All the mutual polarizabilities between atoms of the different ligands are neglected since they have too small values. It was checked by direct test calculation of the whole polarizability matrix of separate ligands.
Then it is possible to calculate the polarized charges in all the required orders by eq. (19) starting from the bare ones. In fact the second order in eq. (19) and full summation of the perturbation series by eq. (18) give very close results. In most cases the second or even the first order polarization suffice for our purposes (see the next Section).
In our present study the basic procedure for treating PES of TMC within a general QM/MM-like framework is constructed. To summarize, we reformulated in the local form (i.e. in terms of the effective field increments induced by the lone pairs of the ligands) the semi-empirical EHCF theory which previously allowed to calculate with quite good accuracy the crystal field induced by the ligand's on the TMC's d-shells. This gave us explicit formulae for the crystal field matrix expressed through the ESP's of the free ligands and the procedure to calculate them. In the framework of our approach the crystal field matrix is calculated for arbitrary arrangement and orientation of the ligands around the central TMI.
In the present work we constructed a procedure combining the EHCF(L) approach and the specific form of the MM (MMGK Ref. ) by eq. (26). It is implemented within the ECFMM 1.1 package MMECFPackage which allows gradient minimization for the energy eq. ( 26). The package allows also to consider ligands or their fragments also as rigid bodies. As a consequence, the number of geometry variables considerably decreases which allows to speed up the minimization. Technically, the ligand geometries employed within the rigid body scheme are first pre-optimized with use of the MM potentials only, and in the further calculations their internal geometry is fixed. First and second orders for charges eq. (19) as well as the exact formula eq. ( 18) of the Perturbative Sparkle model were implemented in the ECFMM package. Parameters fitted within EHCF(L)-PS/MMGK model (where PS stands for Perturbative Sparkle) for pairs metal atom-donor atom, where metal is Fe(II) or Co(II) and donor atom types (MMGK) are NA and N3, that is sp2- and sp3-hybridized nitrogen, are present in Table 1. Names and Cambridge Crystal Structure Data Bank (CCSDB) codes of the calculated complexes are listed in Table 2.
As it is stated in Section 2.1 the effective charges and orbital energies of the l-system are needed to estimate the EHCF. In the EHCF method Ref. SouTchMis these quantities are calculated with use of the semiempirical SCF procedure (CNDO). We compared the free ligand charges and the charges calculated within the SCF sparkle estimates and perturbative sparkle estimates according to procedure described in Section 2.3.2. The results are presented in Tables and . Numeration of corresponding ligand atoms can be found on Fig. 1.
The data of Table show that the free ligand effective atomic charges (column 1) are quite strongly affected by the interaction with metal ion described by the SCF sparkle procedure (column 2). The difference between the results of the SCF sparkle model (2) and the charges obtained by the first order PS model (column 3) described in Section 2 (eqs. ( 17), (18)) is however small especially for charges of the peripheral atoms. Charges on the donor atoms in (3) are slightly larger than these in (2), but the difference is also small. On the other hand, the difference between SCF sparkle (2) and free ligand (1) charges are close to the difference between (3) and (1). One can see that whatever method of taking into account the charge renormalization gives rather close results when it goes about the charges on all the ligand atoms including the donor ones.
The values of the orbital energies are obtained from eq. (14) that is the first-order correction to the orbital energy for the MO in the field induced by all the charges of the complex except those in the ligand under consideration itself. The charges in eq. (14) are taken from the first order PS model. Orbital energies are close to those obtained within the sparkle SCF scheme, which is illustrated by the data presented in Table . Thus, such an approximation for the orbital energies of the ligands is shown to simulate the results of the SCF sparkle model.
In conclusion, the first order polarized point charges as well as the orbital energies eq. (14) using these charges of eq. ( 19) are obtained to be in fair agreement with those from a semiempirical SCF procedure corresponding to the sparkle model. Moreover, the difference between charges on the donor atoms and on peripheral atoms indeed, obtained within the sparkle SCF and PS models, are in all cases smaller than 0.05 which results in less than 2 % possible error in the EHCF estimates for the crystal field.
The methodology described above was applied to 30 complexes of Fe2+ listed in Table 2 together with the ligand names, relevant Cambridge Crystal Structure Data Bank (CCSDB) reference codes, and the experimental spins of the ground states. The ligands are shown in Fig. 1. The series contains compounds with monodentate and polydentate ligands of both low- and high-spin ground states.
Experimental geometries of the above complexes were taken from the CCSDB. Hydrogen atoms were added where necessary. The complexes 23-30 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 (at least in terms of molecular geometry).
As a test, we calculated the 10Dq parameter for octahedral complex 11 as a function of the metal-nitrogen distance with use of the EHCF(L)-PS, and by the original EHCF procedure. It was found that for the ''interesting'' range of the interatomic separations (about 2 Å ) either first or second perturbation orders of the EHCF(L)-PS model employed in the present work fairly coincide with the standard EHCF curve which for the purposes of the present paper is considered as the exact one. Thus, it is not necessary to use the renormalized Racah parameters, which we used in Refs. TchIJQC,DarhPletTchJCC. So, unlike the RLMO model , we are able to use a single parameter set on the metal d-shell for ligands containing both the NA (sp2) and N3 (sp3) nitrogen atom types.
Parameterizations for the NA and N3 atom types were performed separately. The proposed EHCF(L)-PS method was initially applied to a test set of the Fe2+ complexes. The test set comprised complexes 1, 2, 12, 13 (for NA) and 7, 19, 20 (for N3). Structures of the test set with different spin states (singlet and quintet) were optimized by analytical gradient procedure starting from the experimental structures with rigid (preliminary optimized with pure MMGK) ligands. Optimization is performed until root-mean-squared (RMS) energy gradient is smaller than 0.1 kcal· mol-1·Å . Criteria for good parameterization are the correct ground spin state and small difference in molecular geometry of coordination sphere (bond lengths and valence angles on metal atom). Obtained parameters are presented in Table 1. One can see that the rigidity of the metal-ligand bond, measured by D0a2 for the N3 type atom parameter is only slightly smaller than that for the NA type atom. Geometry optimization for the rest of series of the investigated complexes was performed through the same scheme, but with use of parameters already defined.
Below we consider results of our calculations of Fe(II) complexes of Table 2. Metal-donor atom bond lengths and averaged valence angles centered on metal atom as well as their RMS from experimental crystal structure are present in Tables 5, 6 for different spin states of the metal and compared with experimental values. The possible terms or spin states considered in calculations are quintet (5T2g) and singlet (1A1g) prototypes for octahedral coordination of the Fe2+ ion.
Calculated geometries of this series of the complexes, in general, agree rather well with the experimental data. We especially notice that complexes with ligands containing different types of donor atoms (NA and N3) in the single molecule are calculated correctly. However, making parameterization we should keep in mind that crystal structure can be a result of not only metal interaction with the ligands but also of intermolecular interactions among the crystal neighbors. Thus, major deviations from experimental geometry may be a result of crystal surrounding influence, especially, that of counter ions (see discussion in ). Analysis of obtained results shows that geometry structure details are in good agreement with crystal data (rms bond lengths is within 0.06 Å ), despite the fact that in the complexes 4-6 with low-spin experimental structure our method gives wrong ground-spin state but correct geometry for the experimental low-spin state. For the high-spin complexes 11-22 as well as for the low-spin complexes 1-3 and 7-10 we obtain both correct ground spin state and acceptable geometry. It can be concluded that current parameterization of EHCF(L)-PS method is somewhat biased towards the high-spin states. This manifests itself in the fact that the calculated high-spin state equilibrium geometries correspond to noticeably shorter Fe-N bond lengths than the experimental ones. Technically, the reason may be the stiffness of the Morse potential employed to model the Fe-N MM energy increment. For the complexes with spin isomers we obtained similar result, since in all the cases the high-spin form of the isolated molecular metal containing complex cations have lower energy. Analysis of effects of counteranions upon the spin forms in the spin-crossover compounds is given in Ref. . The geometry of both low- and high-spin isomers are calculated correctly with the rms smaller than 0.07 Å . If only the low-spin molecules are considered the rms does not exceed 0.04 Å .
To conclude, we notice that it is the first time when calculation is performed for such a wide range of Fe(II) complexes (27 individual molecules) of different ground state spins within single parameterization, and reproduces ground state spin as well as geometry of the crystal structure with reasonable accuracy.
Within the proposed improvement of EHCF(L) method, we also calculated the Co2+ complexes with different geometry of coordination polyhedra. The results of calculations are given in Table 7. We tried to have maximal diverse test calculation, thus selected octahedral 31-33, tetrahedral 34 and pyramidal 35 high-spin (quartet with the prototype state 4T1g) complexes and the low-spin square pyramidal 38 and square planar complexes 36-37 and 39-40 (doublet with the prototype state 2T1g). The entire set of the Co2+ complexes was used for parameterization of the NA and N3 atom types. With use of the Racah parameters given in Table 7 the correct ground spin states are reproduced for the whole set of the compexes calculated at their experimental geometries. This parameterization allows to obtain even more precise results than for the Fe2+ complexes with the rms for the bond lenghts less than 0.05 Å .
In the high-spin complexes 31-35 structure and spin states are correctly predicted with low deviation with rms smaller than 0.05 Å in bond lengths. Interestingly, the complex 32 is known to be near the spin crossover point in solution . According to our calculations, energy difference between quartet and doublet minima for this complex is the smallest among all the complexes 31-35 and is equal to 3.7 kcal/mol, while in other cases it is more than 10 kcal/mol (see Table 7).
Remarkably enough that the same parameterization allows to reproduce the ground state spin for the low-spin compounds 36-40 provided the calculation is performed with due caution. Particularly, in complexes 37,39 we had to include a water molecule in the calculation according to the experimental crystal structure of the complex. Its position was fixed during the optimization procedure since it is needed for correct calculation of the crystal field. On the contrary, the unit cell of the complex 36 contains two ClO4- counter ions which were considered in  as extremely weakly bonded (the distance Co-O(ClO4-) is 2.409 Å ). By our calculation, we optimized both the complex ion 36 itself and that including two ClO4- anions fixed at their positions in the unit cell and found no difference in geometry and spin of the ground state, so with and without counteranions the result of optimization is the same. For these low-spin complexes 36,37,39 average rms for bond lenghts is 0.04 Å . For the rest of the Co2+ low-spin complexes series, the square pyramidal 38 and the square planar 40, geometry was correctly predicted with the rms smaller than 0.03 Å in the bond lengths.
The results show that the proposed method can be used for precise calculation on geometry and spin states of Co2+ complexes with different coordination numbers and coordination patterns. The proposed methodology thus covers in a uniform way different ground state spins and even coordination numbers of the cobalt(II) complexes.
On the basis of the above analysis it can be stated that the concert usage of the EHCF(L)-PS 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 QM/MM-like description for the PES of different spin states of the iron(II) and cobalt(II) complexes with nitrogen-containing ligands is achieved with use of the single spin-independent parameterization specific for each metal atom and MM type of the donor atom.
The used EHCF(L) procedure allows for a detailed description of the d-shell energy as a function of composition and geometry of the ligand sphere, taking into account the correlation of electrons in the d-shell by using the full CI wave function for them. This allows to handle correctly reaction of the d-shell to subtle changes of the crystal field induced by the surrounding and by this to be sure that the spin intersection point is correctly located.
On the other hand a due attention is paid to reproducing the dependence of the crystal field itself on the tiny ligand geometry and ESP's variations. Explicit form for the crystal field matrix elements reproduce their dependence not only on interatomic separations but also on all kinds of valence angles.
The authors gratefully acknowledge Dr. I.V. Pletnev for valuable discussions and providing the wrapping MMGK suite MMPC1030 . The authors are grateful to Profs. T.W. Hambley and M. Zimmer for sending reprints of their recent works. The usage of the CCSDB is supported by the RFBR grant No 99-07-90133. This work has been supported by the RAS through the grant No 6-120 dispatched by the Young Researchers Commission of RAS.
|Metal||Atom||D0, kcal/mol||a, Å-1||r0, Å||A, kcal· Å 6 /mol||deff|
|no.||ligand name||ground state||CCSD refcode||Ref.|
|no.||ligand name||ground state||CCSD refcode||Ref.|
|Bipyridine in 1|
|Terpyridine in 6|
|Pyridine in 11|
|Methyl-bipyridine in 12|
|Ethylenediamine in 20|
|Bipyridine in 1|
|Pyridine in 11|
1Corresponding author. Address for correspondence: Karpov Institute of Physical Chemistry, Vorontsovo Pole 10, 105064 Moscow, Russia