Hybrid Local Molecular Orbital: Molecular Orbital Calculations for

Subscriber access provided by University of South Dakota

Quantum Electronic Structure

Hybrid Local Molecular Orbital: Molecular Orbital Calculations for Open Shell Systems Milica Feldt, and Ricardo A. Mata J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00727 • Publication Date (Web): 28 Aug 2018 Downloaded from http://pubs.acs.org on August 30, 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 35 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

Hybrid Local Molecular Orbital: Molecular Orbital Calculations for Open Shell Systems Milica Feldt∗,†,‡ and Ricardo A. Mata∗,† †Institut für Physikalische Chemie, Universität Göttingen, Tammannstrasse 6, D-37077 Göttingen, Germany ‡Current address: Division of Quantum Chemistry and Physical Chemistry, Department of Chemistry, Katholieke Universiteit Leuven, Celestijnenlaan 200F, 3001 Leuven, Belgium E-mail: [email protected]; [email protected] Phone: +32-(0)16-325011; +49-(0)551-3923149

Abstract In this work we discuss first applications of the local molecular orbital: molecular orbital (LMOMO) scheme on open-shell systems, focusing on the advantages of isolating the orbital space of (or near) metal centers. We have used as benchmark ligand exchange reactions, discussing the multi-reference character observed in local and canonical calculations, the impact of the local domain approximations and the convergence of the hybrid scheme. After drawing some conclusions on how to build a selection for high level regions, we applied the method to the rate determining steps in a nitrite reductase catalysed reaction step. We have been able to obtain an overall description of the proton, electron transfers energetics occurring in the system which provide a picture for the first time in consistence with the experimental findings. The use of the hybrid scheme is particularly useful in the calculation of the electronic affinities which had previously only been calculated at lower levels of theory. These results show

1

ACS Paragon Plus Environment

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

the potential of such local orbital approaches to deal with specific correlation effects confined to a metal center or their direct vicinity.

1

Introduction

Understanding the reactivity of complexes involving transition metals is a long enduring subject. They are essential in modern synthetic catalysis, new materials, energy applications, but also in nature where over one third of all proteins contain metal ions as essential prosthetic groups. 1 Different quantum chemical methods can be used for assessing such systems and assessing their role in catalysis, but one is usually forced to a hard compromise between accuracy and computational cost. Wave function approaches are often necessary, but these methods can be expensive, particularly since both static and dynamical electron correlation has to be taken into account. This is not a completely lost cause for single reference methods, as some works have recently shown, 2,3 but one is still required to include a large number of configurations to at least qualitatively capture the properties of such systems. However, some of these effects are relatively local in nature, even those arising from static correlation. The biggest obstacles are linked to the metal centers and their direct vicinity, foremost the metal d-electrons and the strongest interaction spaces. This observation raises the possibility of dealing with the problem through local approaches, which restrict the configuration space and can be used to single out sections of a molecular complex. In the Local Molecular Orbital:Molecular Orbital (LMOMO) scheme, 4,5 localized orbitals are used to split the system into different regions. This allows for high accuracy in regions where bond breaking/formation takes place, while the remaining environment is described at a low level. Coupled cluster and Møller-Plesset perturbation (MP2) approaches can be combined in a single calculation, without resource to model systems. 6 It was previously shown that this method provides excellent results when applied to closed-shell systems. 4,5,7 On the example of Mo-dependent enzymes, we were able to show that by treating the metal orbitals at the

2


Page 2 of 35

Page 3 of 35 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


CCSD(T) level, while calculating the rest of the system with MP2, one could reproduce the higher level results for a reaction pathway at a fraction of the cost. 5 Since open shell metal centers commonly occur in enzymes and catalytical complexes, with the closed-shell cases being more of an exception, it is of utmost importance to extend this approach to open-shell systems, albeit the added obstacles. We report here for the first time calculations of the LMOMO approach on open-shell complexes. The major challenges are the onset of static correlation and the potential loss of electron localisation when looking to unpaired electron systems. In this paper, we will first evaluate the LMOMO accuracy on the example of small benchmark models. We make use of binding energies as a reference set, comparing full coupled-cluster to LMOMO results. By considering such ligand exchange reactions, we are able to analyse strong correlation effects in these open-shell systems. After assessing the accuracy for selected models, we apply the method to the calculation of the reaction mechanism of the copper nitrite reductase. 8 The latter involves the calculation of more demanding properties such as electron affinities.

2

Method

The LMOMO approach has been presented in another publication, 4 so we refrain from a lengthy introduction. The total correlation energy can be expressed as a sum of orbital P pair energies ij ij , irrespective of the single reference method in question. Hidden in the aforementioned ij ’s is the dependence on the virtual orbitals which are used to build the many configurations included in the calculation. These could be projected atomic orbitals (PAOs), 9 orbital specific virtuals (OSVs) 10 or pair natural orbitals (PNOs), 11 just to name a few examples. The only important characteristic is that both sets of orbitals (occupied and virtual) are in some way spatially local. In this work, we made use of projected-atomic orbitals (PAOs), as applied in the original work of Pulay, but the following description is independent of this choice.

3



Local correlation methods commonly make distinctions on how the pair energies are computed depending on the distance between the orbitals i and j. For example, in order to reach linear scaling, pairs which are very far apart need to be at some point neglected, otherwise only a quadratic asymptotic scaling is possible. In the case of LMOMO, one distinguishes the pair energies not on their distance but depending on the absolute position of the orbitals. A region of interest is defined, an active site, where the pair energies located there are calculated at the highest level. This could be where a bond is being broken or formed, possibly including the next neighboring atoms. Only pairs with orbitals i and j with both i and j located in the active region will be computed at the full coupled cluster level. Other regions can be defined where lower level methods are used for the respective pair energies, such as CCSD, MP2 or simply neglecting the pairs altogether (Hartree-Fock). A graphical representation of the scheme is shown in Figure 1.

k

"kl

i l

"ij

A

"jl

j

E

Figure 1: Schematic representation of the system partitioning in LMOMO between active site (A) and environment (E). The orbitals i and j are located in the active site and their pair energy εij will be calculated at the highest level of theory. Between k and l, located in the environment, only a lower level approach will be used. The pair energy εjl will usually be computed at the lowest level, but can be easily included into the high level calculation. All calculations involving the active site will include the effect of the amplitudes computed for the environment.

The advantage of this method in comparison to QM/QM methods based on a subtractive approach like ONIOM 12,13 is that there is no need for model systems to be built. The par4


Page 4 of 35

Page 5 of 35 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


titioning (into active site and environment) is carried out at the orbital level. This greatly reduces the artifacts in the calculation due to caps and other corrections, and allows for the choice of very small regions, as small as a single atom. It also allows for a straightforward highly accurate embedding of the wave function calculation within the active site. To illustrate this, let us consider what is perhaps the simplest doubles residual equation, the orbital-invariant MP2 expression:

ij ij Rab = Kab +

X

X kj ik ij fkj . fcb − fik Tab + Tab fac Tcbij + Tac

c

k

(1)

ij Equation (1) shows that the amplitude Tab will couple with all other amplitudes featuring

three orbitals in common with this double excitation. Through the other residuals, ultimately all amplitudes couple, as expected. In ONIOM or other subtractive approaches, the residuals for the active region will only feature orbital pairs located there. In LMOMO, the amplitudes of pairs in the environment are included, even if they are not recomputed. For example, in a LCCSD(T):LMP2 calculation, the LCCSD(T) residuals will feature LMP2 amplitudes from orbital pairs in the lower level region. The LMP2 amplitudes of the environment are also calculated including LMP2 amplitudes from the active region. The corresponding pair energies are later replaced by the higher level result. The same will happen with mixed pairs, where i is in the active region and j in the environment. This is the main reason why LMOMO can be applied with such small active site definitions. It should be noted that some embedding is recovered in an ONIOM like expression

E = EHL (A) − ELL (A) + ELL (A+E),

(2)

whereby EX (Y) is the energy of system Y (A=active region, F=full system) computed at the X level of theory (HL=high level, LL=low level). However, the embedding effect is only included by the low level, as far as it is described in the ELL (A+E) calculation. Similar approaches have been used by other authors in the last years. The method was

5



always close in spirit to DFT-in-DFT, 14,15 whereby the calculation of a DFT subregion is carried out embedded in the density of an environment. The extensions of Manby and Miller et al. 16,17 to WFT-in-DFT and Hegely et al. 18 to WFT-in-WFT result in very similar procedures. The embedding, as mentioned above, is somewhat different. The WFT-inDFT methods are fundamentally different, given that the embedding is provided by a DFT potential. The Kalláy WFT-in-WFT approach which is based on the cluster-in-molecule approach 19 makes use of domains where the embedding is exact, but these corrections are obtained in a subtracting fashion. These will converge asymptotically, but it is unclear how they will compare. Most likely this will also depend on the system under investigation, with aromatic systems being somewhat problematic. 18 An alternative multilevel cluster-inmolecule approach was proposed by Li and Piecuch. 20 The Manby and Miller approach uses projectors to separate the two regions, but the reference wave function for regions A and E are calculated separately (albeit exactly coupled). 16,17 In LMOMO, the partitioning is made on top of a single reference (HF) calculation. The multilevel approach of Neese and coworkers 21 makes use of a similar scheme as in LMOMO, but extends to the use of varying accuracy thresholds among the regions. This is a natural extension when dealing with PNOs. There is no information on how the coupling between the different regions is carried out in their published manuscript. In the Molpro implementation of open-shell LMOMO the single occupied orbitals always need to be included in the high-level region. Although this is not a strict requirement, it stays in line with the algorithmic structure of the LUCCSD(T0) code of Liu and Werner, whereby the latter are indexed separately from the doubly occupied orbitals. 22 Beyond that, we have removed the non-linear contributions to the coupled cluster residuals, as mentioned in the original paper, 4 Equations (2-3), by truncating the respective terms in the doubles residual to the union of orbital domains in the active region. If only the orbitals in the active region are correlated at the coupled cluster level, the computational effort of the respective residuals will become asymptotically independent of the environment system used.

6


Page 6 of 35

Page 7 of 35 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


The LMOMO code for open-shell systems is available starting from Molpro 2015.

2.1

Computational Details

All calculations on the model systems were carried out using mainly the aug-cc-pVTZ basis set 23–25 and corresponding fitting basis. This basis set was not used for H in the calculations of the ligand exchange reactions where cc-pVTZ basis set 23 was used. Furthermore, in calculations containing Cu aug-cc-pVTZ-PP basis set 26 in combination with the effective core potential ECP10MDF 27 was applied. Density fitting was applied in all calculations in combination with corresponding basis sets, except in the case of metals. In the case of metals the def2-QZVPP/JKFIT basis set 28 was used for JKFIT. In the case of nitrite reductase the Dunning cc-pVTZ basis 23,24 set was used for all atoms except for Cu. In the case of Cu the aug-cc-pVTZ-PP basis set 25 in combination with the ECP10MDF 27 effective core potential was used. The density fitting approximation was also used with corresponding basis sets for all atoms except for the metal. Auxiliary basis sets used for Cu are def2-QZVPP/JKFIT 28 and aug-cc-pVTZ-PP/MP2FIT 29 for the Coulomb and exchange part, respectively. In all local calculations the Pipek-Mezey localization scheme 30 was used for orbitals, while orbital domains were determined according to the natural population analysis (NPA) criteria 31 with TN P A = 0.03. The distance criteria was used for the classification of the orbital pairs. Strong pairs were defined within the distance of 3 bohr and close pairs within 5 bohr. Strong orbital pairs were fully included in the coupled cluster part and close pairs were treated exclusively at the MP2 level. All calculations were carried out with a development version of Molpro 2012.2. 32 All geometry optimizations were done using the ORCA program package. 33 In the case of the small open-shell test system the geometry optimizations were performed using B3LYP 34,35 with the def2-TZVP basis set 36,37 starting from the structures from Reference 38. The desired coordination geometry was obtained by constraining the L-M-L angles during the geometry 7



Page 8 of 35

optimizations, except in the case of the [MII (H2 O)2 (SH)(CH3 COO)(Im)]·H2 O (SP) complex. In this case, the starting structure was in an octahedral coordination geometry and the optimization did not converge during the constrained optimization. In the case of the relaxed optimization, one water significantly moved away from the copper and at the end of the optimization was located in the second coordination shell. The rest of the complex slipped into a square pyramidal geometry. In the case of nitrite reductase, the geometry optimizations were performed using B3LYP-D3 34,35,39,40 with the def2-SVP basis set. 36,37 The RIJCOSX approximation 41 was used in combination with the def2-SVP/JK auxiliary basis set 28 to speed up the calculations. In all optimizations the α-carbons were constrained to emulate the strain of the backbone, yet allowing a certain amount of structural flexibility. In the case of CuNiR all shown energies represent the free energies including a cluster model of the active site and were calculated as:

∆Gtot = ∆E(elec) + ∆G(solv;) + ∆G(ZPE) + ∆G(therm);

(3)

where ∆E(elec) is the electronic energy computed with different methods, ∆G(solv;) is the solvation energy which is calculated at the LRMP2 level and represents the difference between COSMO calculations with dielectric constants of 4 and 1. ∆G(ZPE) is the zero point energy (ZPE) and ∆G(therm) is the thermal correction to the Gibbs free energy. The last two terms were obtained from frequency calculations using the B3LYP-D3 method with def2-SVP basis set. It was suggested that entropy effects can be sensitive to the freezing of certain atoms during the geometry optimizations, 8 therefore enthalpies are provided in the Supporting Information. The entropy contributions do not change the energies by more than 2.5 kcal/mol as can be observed by comparing the two sets of values. The observations and conclusions later made in the text are independent of the entropy contributions.

8


Page 9 of 35 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


3

Results and discussion

3.1

Open-Shell Test Systems

To assess the accuracy of the LMOMO method small benchmark systems were used to calculate binding energies. These systems were previously used in a study of Gutten et al., 38 where they investigated the influence of the optimization procedure, the solvation effects and the accuracy of the calculated electronic energy on the results. They observed that RI-MP2 results were generally in good agreement with UCCSD(T). Contrarily, DFT was in stark disagreement with the coupled cluster reference. They made use of seven different functionals: PBE-D, B3LYP-D, TPSS-D, M06, BHLYP, M052X and MPWB1K. Reasonable results were obtained in the case of uncharged ligands, while in the case of charged ligands results deviated by more than 10 kcal/mol. Furthermore, they also concluded that open shell systems have larger errors compared to closed shell systems. Based on these observations, we concluded that this would be a suitable set for benchmark calculations. Different levels of theory were used to calculate binding energies, e.g. UCCSD(T), LUCCSD(T0) and LMOMO. As test systems, four different complexes (Figure 2.) were used with three different metal ions (Fe2+ , Cu2+ and Mn2+ ). These three metals were taken in their high-spin state which is assumed to be their electronic ground state. The ligands were chosen to be of varied size, but also to represent all three common metal-binding atoms in proteins (O, S, N). The ligand exchange energies were calculated corresponding to the process: {L1 + . . . + Ln }c−2 + [M(H2 O)n ]2+ −→ [MLn ]c + nH2 O

(4)

where the {Li }c−2 notation is used to denote the overall charge of the (non-interacting) ligands. In almost all cases the number of ligated water molecules before the reaction was equal to the total number of ligands bound to the metal center after the reaction. It should be pointed out that in the case of the SP complex the starting complex was [MII (H2 O)6 ]2+ in octahedral coordination geometry and the SP complex was in a square pyramidal confor9



Page 10 of 35

mation with one water molecule in the second shell. In this study we calculated the ligand exchange energy ∆Eexc as: ∆Eexc (Ln ) = E([MLn ]c ) + nE(H2 O)− X − E([M(H2 O)n ]2+ ) − ( E(Li ))c−2

(5)

i

Figure 2: Model complexes: (A) [MII (CH3 S)(H2 O)]+ in linear coordination geometry (LI complex); (B) [MII (H2 O)2 (H2 S)(NH3 )]2+ in tetrahedral coordination geometry (TH complex); (C) [MII (CH3 S)(NH3 )(H2 O)(CH3 COO)] in square planar coordination geometry (SQ complex); (D) [MII (H2 O)2 (SH)(CH3 COO)(Im)]·H2 O in square pyramidal coordination geometry with one water in the second coordination shell (SP complex).

Single point calculations were carried out on all reactant and product states using different levels of theory. In the case of the square pyramidal systems (SP) we could not obtain the full UCCSD(T) values. In the LMOMO calculations only orbitals containing the metal center 10


Page 11 of 35 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


Table 1: Ligand exchange energies (in kcal/mol) of the studied metal ions with model sites, ∆Eint (Ln ), (for the Reaction (L1 + ... + Ln )c−2 + 2+ c [M(H2 O)n ] → [MLn ] + nH2 O), Calculated Using Various ab Initio Methods and 2+ [M(H2 O)n ] as References. (a) [MX2 ]+ stands for the [MII (CH3 S)(H2 O)]+ complex. (b) [MX4 ] stands for the [MII (CH3 S)(NH3 )(H2 O)(CH3 COO)] complex. (c) 2+ II 2+ + [MY4 ] stands for the [M (H2 O)2 (H2 S)(NH3 )] complex.(d) [MX6 ] stands for the II [M (H2 O)2 (SH)(CH3 COO)(Im)]·H2 O complex. coord complex LI [CuX2 ]+(a) [FeX2 ]+ [MnX2 ]+ TH [CuX4 ]2+(b) [FeX4 ]2+ [MnX4 ]2+ SQ [CuY4 ](c) [FeY4 ] [MnY4 ] SP [CuX6 ](d) [FeX6 ] [MnX6 ]

UCCSD(T) LUCCSD(T0) LMOMO LRPM2 RMP2 –295.4 –296.3 –295.9 –318.3 –318.3 –265.8 –263.3 –263.8 –259.7 –259.6 –256.8 –255.2 –255.6 –252.9 –252.9 –17.4 –17.5 –17.8 –15.4 –15.4 –6.6 –6.9 –7.5 –6.4 –6.4 –5.2 –5.5 –5.9 –5.0 –5.0 –377.9 –374.9 –375.5 –370.3 –371.7 –367.9 –364.9 –365.2 –363.6 –364.8 –354.0 –351.3 –351.6 –350.2 –351.5 –348.7 –348.6 –344.1 –346.5 –325.0 –325.9 –324.6 –327.6 –322.2 –322.7 –321.7 –324.4

were treated at the coupled cluster level and the rest of the system at the LRMP2 level. In Table 1 the results obtained using different wavefunction methods are shown. First, we compare canonical and local methods to estimate how much the PAO-based local approximations effect the results. The influence of the domain approximation can be estimated by comparing canonical and local MP2 results. In general, there is a reasonable agreement between the two, with errors below 1.5 kcal/mol. The difference is only slightly larger in the case of the square planar coordination geometry (SQ). It has been observed in other systems that this difference can be strongly related to the inherent removal of correlation basis set superposition effects (BSSE). 42 Further proof would require an estimate of counterpoise corrections for the reactions featured. This is, however, hard to define for the complexes. To circumvent this definition problem, we used the gCP correction featured in the HF-3c method, 43 as a quick tool to estimate the relative weight of BSSE in each of the complexes. The latter is a fitted formula depending on the atoms and their distances to provide interand intramolecular corrections. In Figure 3 we plot the differences between local and canon11



ical methods as a function of the latter correction. The plots reveal a reasonable correlation taking into account the approximative nature of the gCP formula. For example, one observes that the square planar complexes will be the most affected, in line with the domain errors.

3

F e S Q

2

[k c a l/m o l]

C u L I

C u S Q

M n L I

F e T H

M n S Q

1

c c -lc c

F e L I

0 M n T H

E

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

Page 12 of 35

-1 C u T H

-2 -1

0

E

1 g C P

2

3

[k c a l/m o l]

Figure 3: Comparison between gCP and the difference between LCCSD(T0) and CCSD(T).

On the other hand, the difference between canonical and local coupled cluster results is larger, albeit no more than 3 kcal/mol. One should note that such errors (measured as relative or even absolute) are similar to the ones reported in the first LUCCSD(T0) benchmarks. 22 In the latter study, however, the investigated systems did not contain any metal centers. Some of the differences can be attributed to the pair approximation. In previous work it was shown that error cancellation between domain and pair approximation is important and can lead to the good agreement between local and canonical coupled cluster. 42 Since the domain errors are very small in these calculations, this cancellation did not occur.

12


Page 13 of 35 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


It should be noted that the larger errors occur in the case of ionic ligands, since in that case the charge is not as well localized as in uncharged ligands. Nevertheless, one finds a fair agreement between local and canonical results, particularly as relative errors. Furthermore, we wish to compare the LRMP2 and LUCCSD(T0) results. The difference between these two methods is less than a few kcal/mol in almost all cases. The only exception is for the linear complex of [CuII (CH3 S)(H2 O)]+ where the difference is formidable, reaching 20 kcal/mol. The same difference was observed in the canonical calculations and in the previous study. 38 On the other side, the difference between LMOMO and LUCCSD(T0) methods is more or less constant and does not amount to more than 1 kcal/mol independent of the investigated system. As mentioned above, DFT gave results which were not in good agreement with the UCCSD(T) values depending on the charge of the ligand. 38 Therefore, we can conclude that the LMOMO method performs better than DFT methods independent of the charge of the ligand. Furthermore, we compare the overall performance of the LRMP2 and LMOMO methods. Normalized Gaussians are plotted in Figure 4, representing the error distributions. The center of the Gaussian represents the average difference between LRMP2 or LUCCSD(T0):LRMP2 and LUCCSD(T0). The width represents the root mean square deviation (RMSD). As it can be seen, an average deviation for both methods is below 1 kcal/mol. However, the main disadvantage of LRMP2 is that the results are much more erratic than those obtained with the LMOMO method as it is seen from the width of the Gaussians.

13



Page 14 of 35

-0 .4 0

L U C C S D (T 0 ):L R M P 2

-0 .0 6

-6

-4

-2

0

L R M P 2

2

4

6

E [k c a l/m o l] Figure 4: Normalized Gaussians centered at the average difference, the width represents the RMSD. LUCCSD(T0) was taken as a reference.

One could in general question the applicability of coupled cluster for this set of reactions, considering that such transition metal complexes could carry some multireference character. The chosen suite of single reference methods would be then ill-advised. One of the ways to determine the magnitude of the multireference character is to analyse the T1 diagnostic as provided by a CCSD calculation. 44 The T1 diagnostic 45 is the Forbenius norm of the single-excitation amplitude vector divided by the square root of the number of electrons correlated: ||t1 || T1 = √ nelec

(6)

and its value should be smaller than the recommended threshold of 0.05 in order to allow 14


Page 15 of 35

the treatment with single reference methods. 44 T1 values were obtained from CCSD(T) calculations as well as from local and LMOMO calculations and they are shown in Figure 5. By looking at the CCSD results one can see that the T1 values do not exceed 0.05, with the largest value in the case of Fe square-planar complex of 0.042. It should be pointed out that all water complexes have small T1 diagnostics. Among the metals, complexes containing Mn have the smallest values, while those containing Fe show the largest values. All in all, we can conclude that these systems either do not have multireference character or the character is weak (in case when T1 is larger than 0.02). The largest deviations between MP2 and coupled cluster also seem to bear no correlation to the T1 values. This would lead us to conclude that the chosen set of reactions would be indeed adequate for the proposed benchmark. C C L U L M L M

0 .0 5

S D C C S D O M O O M O - c a lc

0 .0 3

0 .0 2

1

d ia g n o s tic

0 .0 4

T 0 .0 1

]

2+

]

2+

2+

]

2+

]

2+

]

]

2+

2+

2+

2+

2+

2+

2+

6 6

) 6 ]

[M

n( H

2

e( H [F

u( H [C

n( H

2

2

[M

O

O

nX

) 6 ]

eX

) ] 6

O [C 2

[M

[F

uX

) 4 ]

O

O

nY

) 4 ]

eY

2

[M

u( H

e( H [F

[M

[C

n( H

2

[M

e( H [F

2

2

[C

O

O

[F

) 4 ]

4

4

nX

O

eX [F

2

u( H [C

[C

n( H [M

uY

) ]

] 4

) ]

] 4

) ] O

uX

4

2

4

) ]

2

2

[M

O 2

e( H [F

u( H

2

[C

[F

O

O

2

nX

) ]

2

eX

2

uX

2

) ]

]

6

4

4

4

]

]

]

2+

2+

+ 2+

+

+

0 .0 0

[C

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


Figure 5: T1 diagnostic for systems in this study calculated with different methods. The horizontal line at 0.05 represent commonly used thresholds for T1 diagnostic. 44 The horizontal line at 0.03 represent our recommended threshold for local correlation methods. CCSD, LUCCSD and LMOMO values are obtained from Molpro calculations. LMOMO-calc is recalculated LMOMO value to account for the correct number of electrons as described in text. [MX2 ]+ stands for the II + II [M (CH3 S)(H2 O)] complex. [MX4 ] stands for the [M (CH3 S)(NH3 )(H2 O)(CH3 COO)] complex. [MY4 ]2+ stands for the [MII (H2 O)2 (H2 S)(NH3 )]2+ complex. [MX6 ]+ stands for the [MII (H2 O)2 (SH)(CH3 COO)(Im)]·H2 O complex. T1 values can be also obtained from local calculations as well as from the hybrid scheme.

15



If we look at the T1 diagnostic from the LUCCSD(T0) calculations, we see that the results are mainly underestimated compared to the CCSD(T) results. In the case of small T1 values the difference is insignificant. However, in the case of larger T1 values this difference is found to be around 0.02. Therefore, we would suggest to lower the acceptance threshold for LUCCSD(T0). Instead of the recommended value of 0.05 for the canonical approach we would make use of a more cautious value of 0.03. Furthermore, LMOMO results are in good agreement with LUCCSD(T0) results and they always have smaller value. The largest difference is 0.009. However, in the case of LMOMO the number of all electrons was used to calculate the T1 diagnostic, but instead only the number of electrons in the high level region should be used. When this is corrected for, nelec was replaces by nregion in Equation (6), the difference becomes smaller than 0.005, with only elec one outlier ([FeII (H2 O)2 (H2 S)(NH3 )]2+ complex) where the difference is 0.008. From all this we can conclude that the T1 diagnostic obtained either directly from LMOMO calculations, or corrected for the number of electrons in high level region, can be used to estimate the multireference character of the system if the threshold of 0.03 is used.

3.2

Nitrite Reductase

After benchmarking the hybrid LMOMO scheme on the small model systems, we will demonstrate its usefulness in a computational application featuring the copper nitrite reductase (NiR). The latter reduces the substrate NO− 2 (nitrite) to NO through a multi-step electron and proton transfer process. The question, however, is how the different events are sequenced. To clarify the sequence of events leading to substrate conversion, we used a cluster representation of the active site. This included the T2 copper, which is coordinated to three histidines (His94, His129 and His300-II) as well as Asp92 in the second coordination shell. Two water molecules were also included in the model, since in some of the crystal structures water molecules are present around Asp92. To model this active site the crystal structure of the blue Alcaligenes xylosixidans NiR 46 is used (Figure 6). 16


Page 16 of 35

Page 17 of 35 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


Figure 6: X-ray structure of the catalytically active T2-copper site in NiR from Alcaligenes xylosoxidans. 46

The beginning of the reaction mechanism of copper nitrite reductase was investigated. Three different processes can occur at this stage: protonation, reduction or substrate (NO2 ) binding. In order to know the exact sequence of steps, the different alternatives have been investigated. All possible combinations of these processes were calculated leading to six different mechanisms (Figure 7 and Table 2). The initial and the final states are the same in all of them. In the initial state, the T2 copper is in its oxidized state and Asp92 is deprotonated (OS). In the final state, the T2 copper is in its reduced state and Asp92 is protonated (RS). In order to determine which mechanism is most favorable proton and electron affinities, as 17



well as the binding energy of NO− 2 were calculated. The optimized structures are shown in Figures 8 and 9 for the complexes with oxidized and reduced copper, respectively, using the same nomenclature as in the previous study by Siegbahn and coworkers. 8 The latter depends on the protonation state and whether the substrate is already coordinated. Given the presence of various titratable residues near the T2 site, and not only Asp92, there are several different proposal for the proton transfer. It is generally agreed upon that Asp92 will serve as a proton donor, but this can easily couple to other residues, such as His249. Several different pathways have been proposed and studied by Lintuluoto and Lintuluoto 47,48 looking at the pKa values of neighboring residues and the electron affinities of the T2 site. The latter values are very sensitive to the computational protocol used. A simple choice of effective dielectric can tip the scales between the T1 and T2 sites. 48 All previous works have been based on DFT methods so one has to question how much the results might be influenced by the method itself. Figure 7: Schematic representation of all possible mechanisms for the first three steps of Copper Nitrite Reductase.

In our study we used the same model as Siegbahn and coworkers. 8 In this model NO− 2 coordinates to copper through one of the oxygens (Figure 8 and 9). It should be noted that in a recent theoretical study it was shown that the nitrite could also bind in a bidentate 18


Page 18 of 35

Page 19 of 35 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


Table 2: All possible mechanisms for the first three steps of Copper Nitrite Reductase. NO−

H+

e−

(1)

2 OS −−−→ SS −−→ P1S −→ RS

(2)

2 OS −−−→ SS −→ P1ES −−→ RS

(3)

2 P1S −→ RS OS −−→ P1WS −−−→

(4)

2 OS −−→ P1WS −→ P1EWS −−−→ RS

(5)

2 OS −→ OES −−→ P1EWS −−−→ RS

(6)

2 OS −→ OES −−−→ P1ES −−→ RS

NO−

e−

NO−

H+

e−

H+

e−

H+

NO−

e−

H+

e−

NO−

NO− H+

fashion to the copper. 47,49 Even though the NO binding mode is in agreement with a previous theoretical study by McDouall et. al 50 and our own calculations, the NO− 2 binding mode is in disagreement with our study and that of Siegbahn. 8 Li and coworkers 49 found that this bidentate way of NO− 2 binding occurs only in the case when both residues Asp92 and His249 are protonated. Lintuluoto and coworkers 47 observed that the coordination could depend on the details of the active site model. We have opted for the Siegbahn description and keep with the oxygen coordination through all the stationary points.

Figure 8: Optimized structures of the T2 site with copper in the oxidized state.

19



Page 20 of 35

Figure 9: Optimized structures of the T2 site with copper in the reduced state.

To calculate the proton affinities, pKa values should be calculated and compared with the proton donor. From the relation:

pKa =

∆Go,pKa solv RT ln10

(7)

the pKa is obtained. In the above equation ∆Go,pKa is the free energy of deprotonation and solv RT ln10 = 5.743 kJ/mol. The free energy of deprotonation is defined as the free energy change of the reaction: + HAsolv → A+ solv + Hsolv .

(8)

The value of −1119.4 kJ/mol 51 has been used for the free solvation energy of H + . Furthermore, to determine the electron affinities, the redox potential of different intermediates was calculated using the equation: ∆Go,red solv E =− nF 0

(9)

where ∆Go,red solv is the free energy of reduction, F is the Faraday constant and n is the number of electrons being transferred. The free energy of reduction is defined as the free energy change of the reaction: Oxsolv + e− → Redsolv 20


(10)

Page 21 of 35 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


with Oxsolv and Redsolv being the oxidized and reduced species in solution. To obtain the redox potential which are experimentally relevant a reference electrode must be considered. For this to be done, previously published experimental or theoretical values for the absolute reduction potential of the reference electrode should be used. In this study, the value of 4.281 V obtained by Isse and coworkers 51 was applied.

3.2.1

Calculation of proton affinities

We will start by calculating the proton affinities using LRMP2, LUCCSD(T0) and LMOMO methods. The multireference character was also checked for the CuNiR system. It was seen that this system is of single reference character with the T1 diagnostic from LUCCSD(T0) calculations not exceeding 0.016, well below our suggested threshold value of 0.03 for local coupled cluster. Similar values were obtained with corrected LMOMO diagnostics (see ESI) with the difference to LUCCSD(T0) values not exceeding 0.005. In the LMOMO calculations three different region selections were used (Figure 10) for the calculation of all properties. In the first selection only the metal center was correlated at the high level leading to only 14 orbitals treated at the coupled cluster level. This selection will be denoted as R1. In the second selection (R2 selection) the first neighbors of the metal as well as the substrate were treated at the high level, while in the third one the position where the proton binds was also included (R3). The total number of valence orbitals treated at the high level was 34 and 36 for R2 and R3, respectively. This is a third of the number of valence orbitals in the full LUCCSD(T0) calculation (98 valence orbitals).

21



Figure 10: Selection of high level regions (represented in red) on the example of the SS structure.

The results for the proton affinities are shown in Table 3. Differences between the methods (LRMP2, LMOMO-LUCCSD(T0):LRMP2 and DFT) and the LUCCSD(T0) reference are shown in Figure 11. It can be seen that the difference between R1 and the full calculation is only a few kcal/mol and by increasing the high level region the error is slightly decreasing. In the end, the largest error is less than 3 kcal/mol and it is comparable with the error obtained from the LRMP2 results. We are left to conclude that in the calculations of proton affinities no significant improvement is obtained by treating a part of the system at the coupled cluster level since the LRMP2 results were already close to the target values. Furthermore, DFT results from Reference 8 were also compared with LUCCSD(T0) results. B3LYP in combination with the LACV3P**+ basis set was used in this case. 8 The latter values did not include thermal corrections, but we observe in our own calculations that this correction is no larger than 2 kcal/mol (see ESI). This will not change the results significantly, therefore we compared the values directly. As it can be seen from Figure 11, DFT gives results which are more than 10 kcal/mol away from the coupled cluster results. Even considering an error bar of 2 kcal/mol, the conclusion is the same.

22


Page 22 of 35

Page 23 of 35 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


Table 3: Proton affinities calculated with LUCCSD(T0):LRMP2 in comparison with the LUCCSD(T0) and LRMP2 methods. DFT results are taken from Reference 8. All values are in kcal/mol.

OS SS OES P1ES

3.2.2

LRMP2 LUCCSD(T0):LRMP2 LUCCSD(T0) R1 R2 R3 246.5 245.7 246.5 247.1 248.5 272.6 271.9 272.4 272.9 274.6 277.7 277.0 278.0 278.7 279.3 301.4 300.7 301.0 301.1 304.0

DFT 266.5 284.2

Calculation of electron affinities

Moving to the results on the electron affinities, the picture changes significantly. One can observe that the LRMP2 results deviate strongly from the coupled cluster values, up to 27 kcal/mol (Table 4 and Figure 11). However, in Table 4 it can be seen that the treatment of only the metal orbitals at the high level is enough to almost fully correct the LRMP2 results. By increasing the size of the high level region an accuracy of about 2 kcal/mol can be reached. For the calculation of electron affinities DFT methods were also found unsuitable. Even though the electron affinity of the P1S state was found in perfect agreement with the coupled cluster results the correct order of the electron affinities for these states was not obtained and the errors can vary between 0.5 kcal/mol and 35 kcal/mol (Table 4 and Figure 11). This clearly highlights the difficulty in obtaining these sets of values within a reasonable accuracy range. Table 4: Electron affinities calculated with LUCCSD(T0):LRMP2 in comparison with LUCCSD(T0) and LRMP2 methods. DFT results are taken from Reference 8. All values are in kcal/mol.

OS SS P1WS P1S

LRMP2 LUCCSD(T0):LRMP2 LUCCSD(T0) R1 R2 R3 125.5 102.2 101.4 101.0 100.8 103.1 79.0 78.4 78.5 76.3 156.7 133.5 132.9 132.6 131.5 131.9 107.8 107.1 106.6 105.6

23


DFT 111.1 40.1 106.2


2 0

L R M L M O L M O L M O B 3 L

P 2 M O -R 1 M O -R 2 M O -R 3 Y P

1 0

5

D G

e rr

[k c a l/m o l]

1 5

0

-5 O S

S S

O E S

P 1 E S

3 0

2 0

0

-1 0

e rr

[k c a l/m o l]

1 0

D G

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

Page 24 of 35

L R M L M O L M O L M O B 3 L

-2 0

-3 0

-4 0 O S

S S

P 1 W S

P 2 M O -R 1 M O -R 2 M O -R 3 Y P P 1 S

Figure 11: TOP: Difference between LRMP2, LMOMO and DFT results with the LUCCSD(T0) results for the proton affinities. DFT results are taken from Reference 8 and lack thermal correction to free energy. BOTTOM: Difference between LRMP2, LMOMO and DFT results with the LUCCSD(T0) results for the electron affinities. DFT results are taken from Reference 8 and lack thermal correction to free energy.

24


Page 25 of 35 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


3.2.3

Mechanism of CuNiR

Given that the LMOMO scheme accurately reproduces the full results, we will try to answer the question which mechanism nitrite reductase follows, by using the LMOMO R3 results. For that purpose the pKa values were calculated. A negative value (–14.88) was obtained in the OS state meaning that the protonation is extremely unlikely to occur in this state. The pKa values of the SS and OES state are similar and amount to 3.91 and 8.13. The largest pKa value was obtained for the P1ES state (24.43). All these values are very similar or above the pKa value of 3.9 for the σ-carboxy group of aspartic acid, which was taken as a reference, since the protonation occurs at this position. Therefore, four of the proposed mechanisms are possible, mechanisms (1), (2), (5) and (6). Thus, by only looking into the proton affinities the final conclusion about the mechanism through which the reaction occurs cannot be made. Next, we will take a look at the redox potentials. The redox potential of one state (SS state) was found to have a negative value (–0.88 V) meaning that the reduction should not occur in this step. The values for the other three states are 0.10 V, 0.34 V and 1.47 V for OS, P1S and P1WS, respectively. Obtained redox potentials should be compared with the redox potential of the T1 copper site (247 mV) since we know that the electron is transferred from the latter. 52,53 Comparing these redox potentials we can conclude that also the mechanism in which the OS state is involved is not likely. This leaves us with three possible mechanisms: (1), (3) and (4). Combining the results for the pKa values and redox potentials it can be seen that the reaction could happen only through mechanism (1). The other mechanisms should be rejected due to the too low proton and/or electron affinities. These results are in agreement with experiment, since there is experimental observations 54 that nitrite does not bind to the reduced copper site. In the end, we calculated the energetic of the reaction pathway (1). The results for all three regions and also for LRMP2 and LUCCSD(T0) are shown in Table 5. The difference 25



Page 26 of 35

between the LRMP2 and LUCCSD(T0) is significantly smaller in the case of the SS and P1S states, than in the case of RS state. This is to be expected, since the reduction occurs in the last step, which leads to a significant change in the electronic configuration. It is known that MP2 fails to describe such changes. Therefore, the RS state is our best choice to test our hybrid scheme. As it was the case with electron affinities, the treatment of the copper center at the LUCCSD(T0) level is already enough to improve LRMP2 results and reach the energy obtained by LUCCSD(T0). Increasing the region improves the results even further and in the case of the largest region, the error is below 1 kcal/mol in comparison with 22.7 kcal/mol in the case of LRMP2. In the case of other structures the difference is slightly above 1 kcal/mol (max 1.5 kcal/mol), but always smaller than the difference between the LUCCSD(T0) and LRMP2 results. Table 5: Comparison of the LUCCSD(T0):LRMP2 method with LRMP2 and LUCCSD(T0) methods for the mechanism (1). All values are in kcal/mol. LRMP2 OS SS P1S RS

4

0.00 –35.8 –308.4 –440.2

LUCCSD(T0):LRMP2 LUCCSD(T0) R1 R2 R3 0.00 0.00 0.00 0.00 –35.4 –37.5 –37.5 –37.2 –307.3 –309.9 –310.4 –311.8 –415.1 –416.9 –417.0 –417.5

Conclusions

In this work the extension of the hybrid LMOMO scheme to open shell systems and its first applications were discussed. The results of the hybrid scheme were compared to full coupled cluster for ligand exchange reactions. The latter energies are extremely sensitive to the treatment of correlation and allow for an optimal benchmark of the approximations used in our scheme. It was shown that even for total reaction energies exceeding 300 kcal/mol, LMOMO results deviated less than 1 kcal/mol independent of the metal, ligand and coordination. After a satisfactory accuracy was established in the case of the model systems, LUCCSD(T0)

26


Page 27 of 35 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


and LMOMO methods were applied to the calculation of the first three steps in the reaction mechanism of CuNiR. Of particular interest was the calculation of proton and electron affinities in tentative intermediates, a key to understanding the sequence of proton/electron transfer steps leading to product formation. It was shown that proton affinities were largely independent of the level of theory used. On the other hand, electron affinities are much harder to calculate, with LRMP2 significantly overestimating the computed affinities relative to coupled cluster (about 25 kcal/mol in all calculated states). The LMOMO approach was able to reproduce the full results even when treating limited orbital spaces around the metal. For the systems under study the only viable option is the application of high level correlated wave function methods. The results from the literature show that DFT can lead to large errors, not providing even the correct qualitative picture. 8,38 The LMOMO scheme provides a framework to best apply such correlated methods in the investigation of challenging properties in bioinorganic chemistry, while keeping with the accuracy of full calculations.

Acknowledgement M. F. was supported through a PhD scholarship from the International Research Training Group 1422.

Supporting Information Available T1 diagnostic for systems studied in this work. Proton and electron affinities without entropy correction. Examples of the input file. Cartesian coordinates of all structures used in this work. This material is available free of charge via the Internet at http://pubs.acs.org/.

27



5

Keywords

Electronic structure calculations, Coupled Cluster, Copper Nitrite Reductase

References (1) Andreini, C.; Bertini, I.; Cavallaro, G.; Holliday, G. L.; Thornton, J. M. Metal ions in biological catalysis: from enzyme databases to general principles. J. Biol. Inorg. Chem. 2008, 13, 1205–1218. (2) Kats, D.; Manby, F. R. Communication: The distinguishable cluster approximation. J. Chem. Phys. 2013, 139, 021102. (3) Bulik, I. W.; Henderson, T. M.; Scuseria, G. E. Can Single-Reference Coupled Cluster Theory Describe Static Correlation? J. Chem. Theory Comput. 2015, 11, 3171–3179. (4) Mata, R. A.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106. (5) Andrejić, M.; Mata, R. A. Local Hybrid QM/QM Calculations of Reaction Pathways in Metallobiosites. J. Chem. Theory Comput. 2014, 10, 5397–5404. (6) Andrejić, M.; Mata, R. A. Study of ligand effects in aurophilic interactions using local correlation methods. Phys. Chem. Chem. Phys. 2013, 15, 18115–18122. (7) Li, J.; Andrejić, M.; Mata, R. A.; Ryde, U. A Computational Comparison of Oxygen Atom Transfer Catalyzed by Dimethyl Sulfoxide Reductase with Mo and W. Eur. J. Inorg. Chem. 2015, 2015, 3580–3589. (8) De Marothy, S. A.; Blomberg, M. R. A.; Siegbahn, P. E. M. Elucidating the mechanism for the reduction of nitrite by copper nitrite reductase–A contribution from quantum chemical studies. J. Comput. Chem. 2007, 28, 528–539. 28


Page 28 of 35

Page 29 of 35 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


(9) Pulay, P. Localizability of Dynamic Electron Correlation. Chem. Phys. Lett. 1983, 100, 151–154. (10) Yang, J.; Chan, G. K.-L.; Manby, F. R.; Schütz, M.; Werner, H.-J. The orbital-specificvirtual local coupled cluster singles and doubles method. J. Chem. Phys. 2012, 136, 144105. (11) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (12) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H.-B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678–5796. (13) Bates, D. M.; Smith, J. R.; Tschumper, G. S. Efficient and Accurate Methods for the Geometry Optimization of Water Clusters: Application of Analytic Gradients for the Two-Body:Many-Body QM:QM Fragmentation Method to (H2 O)n , n = 3–10. J. Chem. Theory Comput. 2011, 7, 2753–2760. (14) Jacob, C. R.; Neugebauer, J.; Visscher, L. A flexible implementation of frozen-density embedding for use in multilevel simulations. J. Comput. Chem. 2008, 29, 1011–1018. (15) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F. A Simple, Exact DensityFunctional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564–2568. (16) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F. Accurate and systematically improvable density functional theory embedding for correlated wavefunctions. J. Chem. Phys. 2014, 140, 18A507. (17) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F. Density functional theory embedding for correlated wavefunctions: Improved methods for open-shell systems and transition metal complexes. J. Chem. Phys. 2012, 137, 224113. 29



(18) Hégely, B.; Nagy, P. R.; Ferenczy, G. G.; Kállay, M. Exact density functional and wave function embedding schemes based on orbital localization. J. Chem. Phys. 2016, 145, 064107. (19) Rolik, Z.; Kállay, M. A general-order local coupled-cluster method based on the clusterin-molecule approach. J. Chem. Phys. 2011, 135, 104111. (20) Li, W.; Piecuch, P. Multilevel Extension of the Cluster-in-Molecule Local Correlation Methodology: Merging Coupled-Cluster and Møller-Plesset Perturbation Theories. J. Phys. Chem. A 2010, 114, 6721–6727. (21) Sparta, M.; Retegan, M.; Pinski, P.; Riplinger, C.; Becker, U.; Neese, F. Multilevel Approaches within the Local Pair Natural Orbital Framework. J. Chem. Theory Comput. 2017, 13, 3198–3207. (22) Liu, Y. Linear Scaling High-spin Open-shell Local Correlation Methods. Ph.D. thesis, Institut für Theoretische Chemie der Universität Stuttgart, 2011. (23) Dunning, J. T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (24) Kendall, R. A.; Dunning, J. T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (25) Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3d elements Sc-Zn. J. Chem. Phys. 2005, 123, 064107. (26) Peterson, K. A.; Puzzarini, C. Systematically convergent basis sets for transition metals. II. Pseudopotential-based correlation consistent basis sets for the group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements. Theor. Chem. Acc. 2005, 114, 283–296. 30


Page 30 of 35

Page 31 of 35 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


(27) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials for group 11 and 12 atoms: adjustment to multi-configuration Dirac-Hartree-Fock data. Chem. Phys. 2005, 311, 227–244. (28) Weigend, F. Hartree-Fock exchange fitting basis sets for H to Rn. J. Comput. Chem. 2008, 29, 167–175. (29) Hill, J. G.; Platts, J. A. Auxiliary basis sets for density fitting-MP2 calculations: Nonrelativistic triple-ζ all-electron correlation consistent basis sets for the 3d elements Sc-Zn. J. Chem. Phys. 2008, 128, 044104. (30) Pipek, J.; Mezey, P. G. A fast intrinsic localization procedure applicable for abinitio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys. 1989, 90, 4916–4926. (31) Mata, R. A.; Werner, H.-J. Local correlation methods with a natural localized molecular orbital basis. Mol. Phys. 2007, 105, 2753–2761. (32) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, version 2012.1, a package of ab initio programs. 2012; see www.molpro.net. (33) Neese, F. The ORCA program system. Wiley Interdisciplinary Reviews: Computational Molecular Science 2, 73–78.

31



(34) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (35) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (36) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (37) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials. Theor. Chem. Acc. 1997, 97, 119–124. (38) Gutten, O.; Beššeová, I.; Rulíšek, L. Interaction of Metal Ions with Biomolecular Ligands: How Accurate Are Calculated Free Energies Associated with Metal Ion Complexation? J. Phys. Chem. A 2011, 115, 11394–11402. (39) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (40) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (41) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ’chain-of-spheres’ algorithm for the Hartree-Fock exchange. Chem. Phys. 2009, 356, 98–109. (42) Mata, R. A.; Werner, H.-J. Calculation of smooth potential energy surfaces using local electron correlation methods. J. Chem. Phys. 2006, 125, 184110.

32


Page 32 of 35

Page 33 of 35 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


(43) Kruse, H.; Grimme, S. A geometrical correction for the inter- and intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems. J. Chem. Phys. 2012, 136, 154101. (44) Jiang, W.; DeYonker, N. J.; Wilson, A. K. Multireference Character for 3d TransitionMetal-Containing Molecules. J. Chem. Theory Comput. 2012, 8, 460–468. (45) Lee, T. J.; Taylor, P. R. A diagnostic for determining the quality of singleâĂŘreference electron correlation methods. Int. J. Quantum. Chem. 1989, 36, 199–207. (46) Ellis, M. J.; Dodd, F. E.; Sawers, G.; Eady, R. R.; Hasnain, S. S. Atomic Resolution Structures of Native Copper Nitrite Reductase from Alcaligenes xylosoxidans and the Active Site Mutant Asp92Glu. J. Mol. Biol. 2003, 328, 429–438. (47) Lintuluoto, M.; Lintuluoto, J. M. DFT Study on Nitrite Reduction Mechanism in Copper-Containing Nitrite Reductase. Biochemistry 2016, 55, 210–223. (48) Lintuluoto, M.; Lintuluoto, J. M. DFT Study on Enzyme Turnover Including Proton and Electron Transfers of Copper-Containing Nitrite Reductase. Biochemistry 2016, 55, 4697–4707. (49) Li, Y.; Hodak, M.; Bernholc, J. Enzymatic Mechanism of Copper-Containing Nitrite Reductase. Biochemistry 2015, 54, 1233–1242. (50) Periyasamy, G.; Sundararajan, M.; Hillier, I. H.; Burton, N. A.; McDouall, J. J. W. The binding of nitric oxide at the Cu(i) site of copper nitrite reductase and of inorganic models: DFT calculations of the energetics and EPR parameters of side-on and end-on structures. Phys. Chem. Chem. Phys. 2007, 9, 2498–2506. (51) Isse, A. A.; Gennaro, A. Absolute Potential of the Standard Hydrogen Electrode and the Problem of Interconversion of Potentials in Different Solvents. J. Phys. Chem. B 2010, 114, 7894–7899. 33



(52) Olesen, K.; Veselov, A.; Zhao, Y.; Wang, Y.; Danner, B.; Scholes, C. P.; Shapleigh, J. P. Spectroscopic, Kinetic, and Electrochemical Characterization of Heterologously Expressed Wild-Type and Mutant Forms of Copper-Containing Nitrite Reductase from Rhodobacter sphaeroides 2.4.3. Biochemistry 1998, 37, 6086–6094. (53) Ghosh, S.; Dey, A.; Sun, Y.; Scholes, C. P.; Solomon, E. I. Spectroscopic and Computational Studies of Nitrite Reductase: Proton Induced Electron Transfer and Backbonding Contributions to Reactivity. J. Am. Chem. Soc. 2009, 131, 277–288. (54) Strange, R. W.; Murphy, L. M.; Dodd, F. E.; Abraham, Z. H.; Eady, R. R.; Smith, B. E.; Hasnain, S. Structural and kinetic evidence for an ordered mechanism of copper nitrite reductase-Edited by A. R. Fersht. J. Mol. Biol. 1999, 287, 1001 – 1009.

34


Page 34 of 35

Page 35 of 35 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


Graphical TOC Entry

35


Hybrid Local Molecular Orbital: Molecular Orbital Calculations for

Recommend Documents