Challenges in Multireference Perturbation Theory for the Calculations

Subscriber access provided by Kaohsiung Medical University

Quantum Electronic Structure

Challenges in Multireference Perturbation Theory for the Calculations of the g-Tensor of First Row Transition Metal Complexes Saurabh Kumar Singh, Mihail Atanasov, and Frank Neese J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00513 • Publication Date (Web): 01 Aug 2018 Downloaded from http://pubs.acs.org on August 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Challenges in Multireference Perturbation Theory for the Calculations of the g-Tensor of First Row Transition Metal Complexes Saurabh Kumar Singh†, Mihail Atanasov†§, Frank Neese†*

AUTHOR ADDRESS †

Department

of

Molecular

Theory

and

Spectroscopy,

Max-Planck

Institute

for

Kohlenforschung, Kaiser Wilhelm-Paltz-1, Mülheim an der Ruhr, Germany. §

Institute of General and Inorganic Chemistry, Bulgarian Academy of Sciences, Akad. Georgi

Bontchev Street 11, 111i3 Sofia, Bulgaria

KEYWORDS: ab initio calculations, density functional theory, g-tensor, CASSCF, NEVPT2, Zeeman Interaction, spin-Hamiltonian parameters, electron paramagnetic resonance, transition metal complexes

1 Environment ACS Paragon Plus

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract: In this article, we have studied thirty-four S = 1/2 complexes of first row transition metal complexes in d1, d5, d7 and d9 configurations in an attempt to determine the intrinsic accuracy of the scalar-relativistic complete active space self-consistent field (CASSCF) and N-electron valence perturbation theory (NEVPT2) methods with respect to predicting molecular gvalues. CASSCF calculations based on active spaces which contain only metal-based orbitals largely overestimate the g-values compared to experiment and often fail to provide chemically meaningful results. Incorporation of dynamic correlation by means of the NEVPT2 method significantly improves the transition energies, with a typical error relative to experiment of 2000-3000 cm-1. As a result, a lowering in the g-shift by almost an order of magnitude is obtained relative to the CASSCF results. However, the g-shifts are still overestimated compared to experiment since CASSCF leads to an overly ionic description of the metalligand bond, and hence to spin-orbit coupling matrix elements that are too large. Inclusion of the double d-shell along with appropriate bonding counterparts to the anti-bonding d-orbitals in the active space led to the correct trends in the g-values for all studied complexes, with the linear regression coefficient (R) equal to 0.93 over the whole data set. Various technical aspects of the calculations such as the influence of relativity, importance of picture change effects, solvation effects, and comparison between second-order perturbation and effective Hamiltonian based theories have also been systematically studied. Additionally, g-tensor calculations were performed with five popular density functional theory (DFT) methods (B3LYP, M06L, M06, TPSSh, and PBE0) to compare with wavefunction (WF) methods. Our results suggest that WF based methods are remarkably better than DFT methods. However, despite the fact that WF theory has come a long way in computing the properties of large, open-shell transition metal complexes, methodological work is still necessary for truly high accuracies to be reached.


Page 2 of 40

Page 3 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction: Open-shell transition metal ions are ubiquitous in the fields of catalysis and biology.1 Among their most important features are their peculiar magnetic properties, such as magnetic bistability, which give them numerous potential applications in materials chemistry.2 A range of sophisticated experimental techniques such as electron paramagnetic resonance (EPR), magnetic circular dichroism (MCD), Mössbauer spectroscopy, and inelastic neutron scattering (INS) are among the useful tools which allow insight into the geometric and electronic structure of these systems. Among these methods, EPR is perhaps the most widely used technique.3 In EPR, one observes the interaction of the intrinsic electric and magnetic moments with an external magnetic field, which can be described by:

r r Hˆ spin = β BgS

(1)

r

r

where β is the Bohr magneton, B is the external magnetic field, S is the electronic spin operator, and g is the 3 x 3 matrix described as the Zeeman coupling matrix (g-matrix). The g-matrix parameterizes the Zeeman effect arising from the interaction of the total electronic magnetic dipole moment with an external magnetic field. In reality, the observed gvalues for paramagnetic systems, for either main group or metal complexes, deviate from the free electron values, and one may refer to g-shifts similar to the chemical shifts in NMR spectroscopy. Deviations from the free electron value arise from the presence of spin-orbit coupling, and hence EPR contains a wealth of information regarding the electronic structure of the system. Theoretical tools play an indispensable role in the interpretation of EPR experiments. Quantum chemical methods are not only used to understand trends in experimental data but also to extract the maximum chemical information, including gaining insight into the structure of unknown species, developing magneto-structural correlations, and 3 Environment ACS Paragon Plus


correlating electronic structure to the reactivity of complexes. Thus, it is not surprising that the prediction of EPR spectra from first-principles calculations has a long history and significant general interest.4 Early attempts to compute g-values were based on either the Hückel Approximation5 or semi-empirical methods such as INDO6 and CNDO.7 These methods have been employed on large polynuclear complexes with reasonable accuracy,6d but several crude approximations and the limits of empirical parameterization preclude the general applicability of these methods. In earlier attempts to extract g-values from first-principles calculations, Moores and McWeeny employed restricted Hartree-Fock (HF) calculations along with a perturbation approach to treat spin-orbit coupling.8 In the context of the ab initio framework, Vahtras et al. employed multiconfigurational self-consistent field (MC-SCF) methods for the calculation of g-values and this method has been applied to several real systems.9 However, this method lacks the treatment of dynamic correlation. In order to capture dynamic correlation to a reasonable extent, the MC-SCF method with second order perturbative corrections was employed to extract g-values using the equations provided by Gerloch and McMeeking.10 Meanwhile, the calculation of g-tensors using multireference configuration interaction (MRCI) approach, pioneered by Grein and coworkers, incorporates both static and dynamic correlation and was extensively applied for the calculation of g-tensors in small molecules.11 Both perturbational and variational approaches were employed to connect electronic structure and spin-Hamiltonian parameters.12 10a 10b 9b Density functional theory (DFT) based methods have been regarded as attractive alternatives to post-HF methods as they approximate electron correlation at lower computational cost. In the DFT framework, the g-values are calculated as a first-order property with a two-component approach or as a second-order property with a one-component approach.13 Numerous Kohn-Sham methods are available for calculating g-tensors using different density functionals such as local spin-density approximations (LSDA), generalized 4 Environment ACS Paragon Plus

Page 4 of 40

Page 5 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


gradient approximation (GGA), hybrid functionals, as well as double hybrid functionals with an admixture of exact exchange.14 15

16

The results are quite impressive for organic radicals.

However, these methods are known to significantly underestimate g-values in the case of transition metal ion complexes, for which the quality depends critically on the choice of functional and SOC treatment. 14-17 There has been a growing interest in computing the g-tensors of transition metal complexes using wavefunction based ab initio methods. Complete active space self-consistent field (CASSCF) based methods followed by CASPT2 or NEVPT2 perturbational treatment of dynamic correlation have been found to be the most reliable in computing the energy gaps between the ground and excited states of transition metal complexes and typically perform much better than common DFT methods. Hence, this approach has been one of the preferred choices for the computation of g-values.18 Regardless of the accurate description of excitation energies, there are several cases where CASSCF based methods fail to provide chemically meaningful g-values and this is due to the poor description of metal-ligand covalency at the CASSCF level.19 It was found that in many cases metal-ligand covalency could be recovered by incorporating some of the ligand-to-metal charge-transfer (LMCT) or metal-to-ligand charge-transfer states (MLCT) – thus, conscientious and thoughtful input on the end of the user is needed to obtain chemically meaningful g-values.18b, 19b, 20 Analyzing the performance of CASSCF/NEVPT2 methods in computing the g-value of transition metal complexes demands a careful testing of several important factors such as the choice of active space, impact of basis sets, relativity, solvation effects, and the influence of low-lying excited states. In the present contribution, we have assembled a test set of thirty-four spectroscopically well-characterized complexes of first-row transition metal complexes as a benchmark set for the CASSCF/NEVPT2 method. In the first section, the theoretical background of the spin Hamiltonian (SH) parameters of S=1/2 complexes is briefly described, 5 Environment ACS Paragon Plus


followed by an exposition of the computational methodology and an outline of the ligand field origin and interpretation of the g-matrix. In the second part, the importance of the active space and basis sets, the impact of relativity, and the influence of the low-lying excited states on computed g-values is analyzed. Finally, DFT results using five different functionals, i.e. B3LYP, M06, M06L, TPSSh and PBE0 are reported and compared to the wavefunction results.


Page 6 of 40

Page 7 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Theoretical Background and Computation: Theory of the g-Matrix As discussed in detail elsewhere,21 to second order, the effective Hamiltonian that describes the magnetic sublevels of the electronic ground state can be written as, eff SM H MM Hˆ 1 Ψ 0SM ′ − ′ = δ MM ′ E0 + Ψ 0

∑∆

−1 k

Ψ 0SM Hˆ 1 Ψ KS ′M ′′ Ψ SK′M ′′ Hˆ 1 Ψ 0SM ′

(2)

KS ′M ′′

Here, S represents the spin of the system, while M represents the magnetic quantum number and the label K identifies the different eigenstates of the unperturbed Hamiltonian Hˆ 0 (assumed to be the Born-Oppenheimer (BO) Hamiltonian). The Hˆ 1 term contains the spin and orbital Zeeman effect as well as the spin-orbit coupling ( Hˆ SOC + Hˆ ZE ). The denominators ∆K represent the non-relativistic or scalar relativistic excitation energies (Ek – E0). A necessary requisite for the validity of eq (2) is the condition that the ground state spin-multiplet of dimension 2S+1 must be energetically sufficiently isolated from the remaining excited states. This treatment yields the g-tensor in the following form: ( RMC) (GC) (OZ / SOC) gKL = δKL ge(SB) +∆gKL +∆gKL +∆gKL

(3)

Here, ∆gRMC, ∆gGC and ∆gOZ/SOC represent the shift in the g-value due to the relativistic mass correction (RMC), gauge correction (GC) and Zeeman/spin-orbit contribution (OZ/SOC) ( SB )

respectively. The first term g e

represents the first order free-electron g-value contribution

to the g-matrix. The isotropic relativistic mass correction21d is given by

( RMC ) g KL = δ KL

α 2 1 ge 2 S 2

Ψ 0SS

∑∇ s

2 i zi

Ψ 0SS

i


(4)


Page 8 of 40

Here, ge is the free electron value (2.00239), α is the fine structure constant, while S is the total spin of the ground state. The gauge correction21a to the g-value is given by

( GC ) g KL =

1 Ψ 0SS S

∑ ξ (r ){r r − r iA

iA i

r }s zi Ψ 0SS

(5)

iA; K i ; L

i, A

where ξ(riA) represents the radial part of the effective charge reduced SOC operator, riA = ri – RA is the position of ith electron relative to the nucleus A and szi is the z-component of the ith electron. This term is a crude and simple approximation to the more general ‘diamagnetic’ gtensor contribution that contains one- and two-electron terms.11a The approximation is warranted because this term is small. The largest contribution to the g-value is the Zeeman interaction/spin-orbit contribution (OZ/SOC) and is represented by

( OZ / SOC ) g KL

  Ψ 0SS  1 =− ∆ b−1  S b( S =S )  b SS + Ψ 0 

∑

∑

li ; K Ψ bSS

Ψ bSS

i

∑z

∑

z L ;i s z ;i Ψ 0SS

i

s

K ;i z ;i

i

Ψ

SS b

Ψ

SS b

∑l

i;L

Ψ

SS 0

i

      

(6)

Here, Ψ0 represents the ground state wavefunction, while Ψb represents an excited state of SS

SS

the same multiplicity. li;K is the Kth component of angular momentum operator for electron i (referred to the global origin) while li;L is the Lth component of the angular momentum operator for the electron i. It is well known that the g-tensor in this formulation is not gauge invariant. Here we follow Luzanov et al.,22 who proposed to choose the center of electronic charge as a center of common origin.21a In order to evaluate the second-order contribution to the g-tensor exactly, one would have to compute all eigenstates of the BO Hamiltonian explicitly, a task which is clearly impossible. Hence, one must either truncate the sum-over-states (SOS) in an approximate treatment of the BO problems or rely on linear response theory in most practical 8 Environment ACS Paragon Plus

Page 9 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


applications.9b In the linear response language, the second order contribution to the g-matrix can be evaluated as

( OZ / SOC ) g KL =−

1 S

∑ µυ

α −β ∂Dµυ

∂BK

µ zL υ

(7)

α −β where Dµυ is the spin-density matrix in terms of the orbital basis sets of the calculations. It is

important to stress that both the sum-over-states and linear response formalisms are only valid if the approximation of a well-isolated ground state is satisfied. For transition metal ions with orbitally degenerate ground states or heavy elements with very large SOC interactions, this condition is not fulfilled. Instead, one relies on an infinite order treatment of the SOC in terms of the quasi-degenerate perturbation theory (QDPT),23 where the excited states are computed explicitly and perturbations are treated variationally by diagonalizing the matrix of

Hˆ BO + Hˆ SOC in the basis of states of Ψ SM I (here, I represents all the roots and M stands for the magnetic quantum number). The matrix elements of basis set of Ψ I

SM

in QDPT treatment can

be written as ′ ′ Ψ SM Hˆ BO + Hˆ SOC Ψ SM = δ IJ δ SS ′δ MM ′ E I( S ) + Ψ SM Hˆ SOC Ψ SM I J I J

(8)

It is important to point out that for well-isolated ground states the truncated sum-over-states and QDPT treatments lead to nearly identical results since the SOS formalism can be derived from QDPT. In QDPT, the g-tensor is defined by the equations first put forward by the Gerloch and McMeeking10b and in an ab initio context by Bolvin.10a



Computational Methodology All calculations are performed using the ORCA 4.0 program.24 Geometry optimizations employed the BP86 functional25 and the all-electron def2-TZVP basis sets26 for all atoms (metal and ligand). Dispersion corrections are treated by Grimme’s atom-pairwise dispersion correction approach as implemented in ORCA.27 Tight SCF (1x10-8 Eh) and tight geometry optimization criteria were used during the geometry optimization. The computation of the numerical Hessian with no negative eigenvalue indicates that the optimized geometries are minima on the potential energy surface. Since most of the complexes are either tri- or dianionic in nature, we have employed the Conductor-like Polarizable Continuum Model (CPCM) to completely screen the charge density of the complexes.28 The optimized geometries of all complexes are provided in the ESI. The resolution of identity (RI) approximation 29 was used along with the corresponding auxiliary basis sets to speed up the calculations. For the calculation of properties we have employed the CASSCF method.30 Minimal active spaces comprising only d-electrons consisted of CAS(9,5) for Cu(II) and Ni(I) complexes, CAS(7,5) for low-spin Co(II), Fe(I) and Ni(III) complexes, CAS(1,5) for V(IV), Cr(V) and Mn(VI) complexes, and CAS(5,5) for low-spin Fe(III) and Mn(II) complexes. To obtain the correct electronic picture and precise g-values, we have systematically varied the active space; details of the active space selection are provided in the results and discussion section. Scalar relativistic effects are incorporated by means of the Douglas-Kroll-Hess (DKH) approximation.31 Here, we have employed the DKH adapted version of the def2TZVP basis sets for all atoms.32 To capture the effect of the dynamic correlation, NEVPT233 was employed on top of the converged CASSCF wavefunction. All CASSCF calculations were converged to the numerical precision of the computer (“EXTREMESCF” in ORCA). As an initial guess to the CASSCF calculations, we have utilized quasi-restricted orbitals (QROs) from DFT calculations.34 10 Environment ACS Paragon Plus

Page 10 of 40

Page 11 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To investigate the influence of the basis set choice, we have performed calculations with the DKH adapted versions of the def2 basis sets that include DKH-def2-TZVP, DKHdef2-TZVPP and DKH-def2-QZVPP.32 Relativistic contracted atomic natural orbitals (ANORCC) basis sets of double-ζ (ANO-RCC-DZP), triple-ζ (ANO-RCC-TZP), and quadruple-ζ (ANO-RCC-QZP) are also used in the calculations.35 The SOC operator is considered as an effective reduced one-electron operator treated by the spin-orbit mean-field (SOMF) approximation.36 The importance of picture change effects on the SOC operator was investigated according to the formalism of ref 14a (see ESI for details). In order to compute the g-values using DFT methods, calculations were performed using the five different functionals PBE0,37 B3LYP,25a M06,38 M06L,39 and TPSSh.40 The computed g-values and description of input files are provided in the ESI. In all calculations, the DKH Hamiltonian retaining one-center terms was employed to account for scalar relativistic effects. To fairly compare with ab initio results, we have employed the DKH adapted version of the def2-TZVP basis sets for all atoms. Both the integration grids and integration accuracy were increased to 7 for transition metal centers, and 6 for all other elements (in ORCA nomenclature).



Ligand Field Theory of g-matrix: Ligand field theory has served as an important tool to understand and interpret the optical and EPR properties of transition metal complexes. Before we begin our discussion on the ab initio computed g-values in transition metal complexes of dN configuration, we briefly summarize the ligand field interpretation of g-values. It is well known from the theoretical foundations of the SH parameters that g-shifts in transition metal complexes solely emerge from d-d excitations. The ligand field expression for the shift in the g-value is given as,41

g KL ≈ g eδ KL +

ζ ion 

∑

∆ i−,1bci2cs2 d i lKM d s d s lLM d i −  S0  i→s

 ∆ −s ,1bcs2ca2 d s lKM d a d a lLM d s  (9) s→a 

∑

where ζion is the one-electron spin-orbit coupling constant of the metal ion, S0 the spin of the ground state, di, ds and da the corresponding doubly occupied, singly occupied, and virtual orbitals, respectively, and c is the MO coefficient of the corresponding orbital. The ∆i,b and ∆s,b are the excitation energies for the promotion of an electron from di to ds and ds to da orbitals respectively. The contribution to g-values emerges only from excited states with the same multiplicity as that of the ground state.41b It is evident from eq. 9 that excitations from doubly occupied molecular orbitals (DOMOs) to singly occupied molecular orbital (SOMOs) always yield a positive contribution to the g-shift, while SOMO to virtual orbital (VMO) excitations lead to a negative contribution to the g-shift. The applicability and accuracy of the above equation critically depends on: (a) correct excitation energies (affecting the denominator) and (b) a proper description of the metal-ligand covalency (both central field and symmetry restricted covalency effects will contribute to the g-shifts41b). For the d9 configuration (i.e. Cu(II) and Ni(I)) complexes in square planar ligand field environment with ligands placed along –x and –y axes, the unpaired electron resides in the metal dx2-y2 based molecular orbital. It is readily derived in a ligand field framework that the contributions to the g-matrix are: 12 Environment ACS Paragon Plus

Page 12 of 40

Page 13 of 40 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


g zz = g e + 8∆ xy → x 2 − y 2 ζ M α12 β12 g xx = g yy = g e + 2∆ yz / xz→ x 2 − y 2 ζ M α12γ 12

(10) (11)

where α1, β1 and γ1 are the MO coefficients of the dx2-y2, dxy and dxz/dyz orbitals, while the various ∆’s are the transition energies in equation (10-11). For predominantly ionic bonding α2i, β2i, γ2i ~ 1, while significant covalent delocalization leads to α2i (and to a much lesser extent β2i, γ2i) < 1. A classic example that can be explained by these ligand field expressions are the g-values of the D4h [CuCl4]2- complex (g|| = 2.232. and g⊥ = 2.049). Closely related expressions for the g-tensor hold for d1 complexes, where the unpaired electron resides in the metal dxy based MO. However, the signs of the all g-shifts are reversed in accord with the discussion above. Low-spin d7 complexes, such as those of Co(II) and Ni(III), are treated in a tetragonally distorted environment with a weak ligand field on the z-axis. This leads to a vacant dx2-y2 orbital followed by the dz2 orbital as the SOMO. The resulting expressions for the g-tensor may be expressed as41a

g zz = g e

(12)

6δ12γ 12ζ M g xx = g e + ∆ yz → z 2

(13)

6δ12γ 22ζ M g yy = g e + ∆ xz→ z 2

(14)

There is no contribution to gzz as the dz2 → dx2-y2 transition does not contribute to the g-value. In the case of low-spin d5 complexes, the ground state is 2T2g in an octahedral ligand field environment and in case of degeneracy of three t2g orbitals, the simple perturbation techniques are not applicable. In case of low-symmetry, i.e. distorted octahedron, if we 13 Environment ACS Paragon Plus


assume that the unpaired electronic resides in the dxy orbital and ζM

Challenges in Multireference Perturbation Theory for the Calculations

Recommend Documents