- Email: [email protected]

QM/MM calculations of EPR hyperfine coupling constants in blue copper proteins Seongho Moona,*, Serguei Patchkovskiib, Dennis R. Salahuba,c,* a

De´partment de Chimie, Faculte´ des Arts et des Sciences, Universite´ de Montre´al, Montre´al, Que., Canada H3C 3J7 b Steacie Institute for Molecular Science, National Research Council of Canada, Ottawa, Ont., Canada K1A 0R6 c University of Calgary, 2500 University Drive NW, Calgary, Alta., Canada T2N 1N4 Received 4 November 2002; accepted 16 December 2002

Abstract We examine the hyperfine structure of the paramagnetic active sites in two oxidized blue copper proteins (azurin and stellacyanin), using a combined quantum mechanics/molecular mechanics (QM/MM) approach. In our implementation, the QM region (treated with density functional theory (DFT)) is truncated by one-electron carbon capping potentials. The effects from surroundings are modeled by including the electric field, due to the MM domain, in the DFT Hamiltonian. The hyperfine structure of the active sites is well described by a small QM region, embedded in MM partial charges. Because of the simplicity of the implementation, the present approach enables us to investigate the hyperfine parameters of large, realistic models of biological active sites. q 2003 Elsevier B.V. All rights reserved. Keywords: Azurin; Stellacyanin; Electron paramagnetic resonance; Qm/mm

1. Introduction Modern high-field magnetic spectroscopy techniques, such as electron paramagnetic resonance (EPR), electron nuclear double resonance (ENDOR), electron spin echo envelope modulation (ESEEM), and nuclear magnetic resonance (NMR), allow the study of the hyperfine structure of oxidized blue copper proteins [1 – 9]. These experimental methods provide a wealth of data on individual nuclei within the active center of the protein. At the same time, notwithstanding these data, complete hyperfine * Corresponding authors. E-mail addresses: [email protected] (S. Moon), [email protected] (S. Patchkovskii).

structures of the active sites are not yet available for any blue copper protein [6,7]. Recently, several theoretical studies for the hyperfine parameters of the copper active sites in azurin and its mutants were carried out using density functional theory (DFT) calculations on model systems with about 100 atoms [10,11]. The calculation results indicate that the electron spin density is largely localized on the copper center and the directly bound aminoacids. DFT calculations for truncated model systems provide useful insights into the electronic structure and magnetic properties of the copper active site. At the same time, it should be noted that the copper active site is surrounded by both the protein backbone and solvent molecules. This

0166-1280/03/$ - see front matter q 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0166-1280(03)00306-3

288

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

environment changes from one system to another (e.g. azurin and stellacyanin, Figs. 1 and 2). Despite tremendous advances in computer technology and computational techniques, models containing around 100 atoms, using basis sets of moderate quality, still represent substantial undertakings. Extension of the model size, to treat additional backbone residues and environment effects at the DFT level, would be computationally expensive. The problem is compounded by the pronounced sensitivity of the hyperfine parameters to the basis set size [12]. The main goal of the present work is to find a way of reducing the size of the DFT model required for an adequate description of the blue copper active site, and at the same time include the effects from surroundings at a reasonable computational cost.

Fig. 2. The QM part in stellacyanin from Cucumis sativus; terminal carbon atoms were replaced by capping carbon atoms in the QM/MM calculations.

Fig. 1. The QM part in wild type azurin from Pseudomonas aeruginosa; terminal carbon atoms were replaced by capping carbon atoms in the QM/MM calculations.

Because saturated carbon atoms of the protein backbone attenuate delocalization of the spin density [9], it should be possible to treat the immediate environment of the copper center with quantum mechanical (QM) techniques, and its further surroundings with molecular mechanics (MM). In combined QM/MM approaches, two issues are critical for the success of the simulations: separation of the model into the QM and MM parts, and handling of the long-range interactions between the subsystems. Here, we adopt the capping potential approach [13] for capping the dangling bonds in the QM region, created at the subsystem boundary. In this technique, bonds are terminated by one-electron ‘carbon’ atoms, carrying empirically parameterized effective core potentials (ECP). This approach is implemented using conventional ECP expansions, and can be incorporated in most DFT programs with minimal code modification. Additionally, treatment of the interactions between the subsystems can be simplified greatly [13]. In general, the interactions (Coulomb, induction, repulsive) between the subsystems can be described by introducing additional one-electron terms in the Hamiltonian of the QM region [14]. Recently, Cui et al. have shown that the electrostatic effect of the environment can be represented by MM partial charges.

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

This description appears to be adequate in NMR chemical shift calculations, as long as the MM part is not directly connected to the NMR nucleus, and polarization effects dominate interactions between the QM and MM subsystems [15]. Here, we adopt a similar approach to the hyperfine parameter calculations in two copper blue proteins: azurin and stellacyanin. The rest of this contribution is organized as follows: Section 2 outlines the theory of the hyperfine parameters in a QM/MM framework. Section 3 provides computational details. Section 4 presents the results of our calculations. Finally, Section 5 provides the conclusions, and outlines the directions for future work.

(see Ref. [18] for the explicit expressions and discussion). Indices u and v denote Cartesian components of vectors and tensors (u; v ¼ x;y; or z). Provided that standard exchange-correlation functionals, independent of the current density, are used, spin densities ra and rb are the unperturbed ground state densities [19,20]. The hyperfine tensor for a magnetic nucleus N is given by the second-order response of the electronic energy EKS with respect to magnetic moments mN and me : In the uncoupled DFT approximation this leads to the expression (up to first order) [17,18,28]: " # occ X ›2 EKS 0 ANuv ¼ ¼ kf0is lh^ 11 Neiuv lfis l ›mNu ›meiv m ¼m ¼0 is N

e

ð3Þ

2. Theory 2.1. Hyperfine interactions in DFT framework In the presence of nuclear and electronic magnetic moments mN and me ; the non-relativistic, spin-unrestricted Kohn – Sham (KS) electronic energy [16] is given by [17] EKS ðm~N ; m~e Þ ¼

occ X

nis hfis ðm~N ; m~e Þj 2 12 72

is

þ h^ 0 ðm~N ; m~e Þjfis ðm~N ; m~e Þi þ J½ra þ rb

ð þ EXC ½ra ; rb þ vð~rÞðra þ rb Þd~r ð1Þ and h^ 0 ðm~N ; m~e Þ ¼

289

X

nuc X

mNu h^ 10 Nu þ

u,x;y;z N

þ

meiu h^ 01 eiu

u,x;y;z i

nuc elec X X X u;v,x;y;z N

X

eX lec

mNu h^ 11 Neiuv meiv

ð2Þ

i

where nis is the occupation number of the KS molecular orbital (MO) fis of spin s (a or b), 72 is the Laplacian operator, J and EXC are the Coulomb and exchange-correlation energies, vðrÞ is the nuclear potential, ra and rb are the spin-up and spin-down electron densities, operator h^ 0 includes contributions to one-electron Hamiltonian, which are linear and bilinear in the electronic and nuclear magnetic moments

where f0is are unperturbed KS MO ði ¼ 1…; Nocc ; s; s ¼ a; bÞ; and h^ 11 Neiuv is the perturbation operator of the hyperfine tensor. This operator arises from a contact mechanism due to the finite probability of finding electrons directly at the nucleus, and the long-range spin dipole –dipole interactions. Additional terms will appear in Eq. (3) if higher order terms [17] and relativistic effects [21] are taken into account. These effects can be neglected for the light nuclei, examined here. The explicit form of the operator h^ 11 Neiuv is given by [17,18] 8p r r d 2 3r r h^ 11 dðriN Þ^siv þ iN iN uv 5 iNu iNv s^ iv Neiuv ¼ 3 riN

ð4Þ

where the first and second terms represent the Fermi contact (isotropic) and the dipole – dipole (anisotropic) interactions, respectively. After the necessary substitutions and units conversion, the Fermi contact and dipole – dipole interaction terms in non-relativistic, spin-unrestricted DFT are given by [17,18] Aiso ðNÞ ¼ gbgN bN kSz l21 occ D E 8p X £ f0is dðriN Þ^siz f0is 3 is ¼

basis X a2b 4p gbgN bN kSz l21 Pmv kxm ldðriN Þlxv l 3 mv

ð5Þ

290

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

and Tuv ðNÞ ¼ gbgN bN kSz l21 occ X 0 riN riN duv 2 3riNu riNv £ fis s^ iz f0is r5

3. Computational details

iN

is

basis X

1 b ¼ gbgN bN kSz l21 Pam2 v 2 mv r r d 2 3r r £ xm iN iN uv 5 iNu iNv xv riN

ð6Þ

where g is the free electron g-value (2.00231931), gN is the g-value of nucleus N; b is the Bohr magneton, bN is the nuclear magneton, kSZ l is the expectation value of the z-component of the total electronic spin, d b is the Dirac delta function, Pam2 is the spin density v matrix, s^Z is the one-electron spin operator for the zdirection, and x denotes atomic basis functions. 2.2. QM/MM: electrostatic contributions to hyperfine tensor The final forms of the isotropic and dipolar hyperfine interactions indicate that the spin density matrix is the only part which can be affected by the MM part. In the simple electronic embedding model of the QM region, the QM electron density is polarized by the MM partial charges. Using the capping carbon approach [13], the magnetic field-free KS energy of the QM region embedded in MM electrostatic fields can be written as follows EKS ðtotalÞ¼EKS ðQMÞþEKS ðelectrostatics;QM=MMÞ * + occ ð X 1 2 ¼ nis fis 227i fis þ vð~rÞðraþrb Þd~r is þJ½raþrb þEXC ½ra ;rb

* + occ X X qJ X eff VK ðriK Þ2 þ nis fis f r is is J[MM iJ K[CAP ð7Þ where VKeff is the ECP of capping carbons and qJ are partial charges of the MM atoms. The last two terms in Eq. (7) are independent of magnetic fields, and vanish upon differentiation of the total KS energy (Eq. (3)). As a result, VKeff and qJ =riJ affect the hyperfine coupling tensors only indirectly, through the changes to the spin density of the QM subsystem.

In this work, hyperfine tensors were calculated in the QM/MM framework. The QM part was treated by DFT, using the deMon program [22]. The MM part was treated with an empirical AMBER 1994 force field [23]. The capping carbon approach was adopted for the boundary atoms. Spin-unresticted Kohn – Sham calculations were performed, with Becke_exchange [24] and Perdew_correlation functionals [25]. A double zeta plus polarization (DZVP2) basis set, augmented by ECPs, p-functions, (16s10p5d)/ (9s6p4d), was used for copper [12]. The IGLO-III basis set was used for the main group atoms [26]. Molecular structures of azurin and stellacyanin were taken from X-ray diffraction results (4AZU and 1JER files from the Protein Data Bank). The remaining hydrogen atoms were generated using the AMBER Xleap program [27], without reoptimizing the geometry. For azurin, four aminoacids (His46, His117, Cys112, and Met121) were included in the QM calculation. Weakly bound Gly45 was not considered in the QM calculation, but was treated as part of the MM subsystem. For stellacyanin, four aminoacids (His46, His94, Cys89, and Gln99) were involved in the QM calculation. In both cases, the aminoacids of the QM region were truncated at the a-carbon. The only side-chains were treated at the QM level (Figs. 1 and 2). Two approaches were used for terminating the polypeptide chain. In the hydrogen link atom approach, hydrogen atoms were added to the MM side of the broken C – C covalent bond of the chain to satisfy the valency of the QM ˚ ). region (here, the C – H bond length is set to 1.09 A In the capping carbon approach, broken bonds were terminated by one-electron carbon atoms, carrying empirically parameterized ECPs. All the remaining atoms of both proteins were treated as electronically innocent bystanders, using MM partial charges. Solvent effects were treated by including the consensus crystallographic water molecules in the MM subsystem.

4. Results and discussion We begin by examining the influence of MM partial charges on the hyperfine structure and

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

Mulliken atomic spin population of the QM region. To this end, we employ two model systems, labeled as ‘C-Cap QM’ and ‘QM/MM’. The first model (C-Cap QM) consists of the pure QM subsystem, capped with carbon ECPs. It contains 42 atoms. The QM/MM model additionally incorporates MM partial charges (2140 and 1929 MM atoms for azurin and stellacyanin, respectively). Results, obtained with these two models, are compared to the experimental data [6 – 9], and to published QM results for an extended model, capped by hydrogen link atoms (H-Link L-QM) [10,11]. In the latter case, four aminoacids, directly coordinated to the copper center, were supplemented by the atoms of the neighboring

291

aminoacids, for a total of 117 and 93 QM atoms in azurin and stellacyanin, respectively. DFT results for the extended model of stellacyanin, calculated in the present work, were obtained for a slightly smaller model, where the peptide chain of the axial aminoacid (Gln99) was not elongated. Tables 1 and 2 compare calculated hyperfine parameters of azurin and stellacyanin for these model systems. Mulliken spin populations of the two models, studied presently, are shown in Table 3. From Tables 1 and 2, it can be seen that the capping pseudopotential technique gives similar results to the hydrogen link atom approach, but at a significantly lower computational cost (in terms of numbers of

Table 1 Calculated hyperfine parameters of azurin in MHz

His46

a

c

H-link L-QM1a

H-link L-QM2b

Exp.c

2109 147 117 2263

2106 144 114 2258

2105 152 138 2290

2207 193 154 2347

2158

, 2 100

Nd1

Aiso T1 T2 T3 Aiso T1 T2 T3 Aiso Aiso Aiso

15.6 21.7 21.4 3.1 0.60 20.14 20.04 0.18 0.68 0.89 0.41

15.6 21.7 21.4 3.1 0.56 20.14 20.04 0.17 0.63 0.96 0.33

18.7 22.2 21.9 4.1 0.62 20.18 20.05 0.23 0.87 1.30 0.34

20.3 22.0 21.6 3.6 0.7 20.2 0.0 0.2 0.9 0.4 0.3

16.1

19.1 22.1 21.7 3.8 1.01 20.19 20.05 0.24 0.68 1.48 0.70

18.5 22.0 21.6 3.7 0.93 20.18 20.05 0.23 0.62 1.43 0.66

22.0 22.4 21.9 4.3 1.21 20.19 20.07 20.25 0.83 1.77 0.77

22.6 22.2 21.7 3.9 1.0 20.2 20.1 0.3 1.5 0.4 0.4

22.4

Hd2 H11 H12

Aiso T1 T2 T3 Aiso T1 T2 T3 Aiso Aiso Aiso

Hb1 Hb2

Aiso Aiso

Nd1

N12

b

QM/MM

Aiso T1 T2 T3

Hd2 H11 H12

Cys112

C-cap QM

Cu

N12

His117

H-link QM

34 60

34 64

UB1LYP results of the large model taken from Ref. [10]. BP86 results of the large model taken from Ref. [11]. ENDOR, ESEEM and NMR results taken from Refs. [5–7,9].

30 54

28 46

0.9

0.79 1.32 0.60

1.5

2.00 0.79 1.06 28 32

18.1 21.3 20.8 1.2 0.87 20.17 20.07 0.23 1.49 1.06/1.48 0.56 25.1 21.5 21.1 2.7 1.30 20.32 20.04 0.35 1.61 1.45/1.02 27 28

292

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

Table 2 Calculated hyperfine parameters of stellacyanin in MHz

His46

a

c

M-QM/MMa

H-link L-QMb

273 176 113 2289

270 172 112 2284

253 190 100 2289

253 198 106 2304

259 175 93 2263

Nd1

Aiso T1 T2 T3 Aiso T1 T2 T3 Aiso Aiso Aiso

16.93 22.15 21.75 3.90 1.04 20.18 20.03 0.22 1.56 1.26 0.95

16.63 22.12 21.73 3.85 1.02 20.18 20.03 0.21 1.51 1.28 0.89

16.49 22.06 21.66 3.72 0.96 20.17 20.04 0.21 1.48 1.25 0.76

17.75 22.24 21.82 4.05 1.04 20.18 20.04 0.22 1.65 1.32 0.81

18.30 22.45 21.97 4.42 1.04 20.21 20.03 0.24 1.46 1.12 0.77

Hd2 H1 1 H1 2

Aiso T1 T2 T3 Aiso T1 T2 T3 Aiso Aiso Aiso

23.92 23.10 22.49 5.59 1.71 20.23 20.06 0.29 1.82 1.53 1.69

23.54 23.04 22.45 5.49 1.64 20.22 20.06 0.29 1.74 1.50 1.67

25.33 23.28 22.63 5.91 1.68 20.25 20.08 0.32 1.79 1.61 1.75

26.43 23.50 22.82 6.32 1.75 20.26 20.08 0.34 1.98 1.72 1.81

26.79 23.39 22.71 6.10 1.61 20.26 20.06 0.32 1.56 1.36 1.38

Hb1 Hb2

Aiso Aiso

Nd1

N1 2

b

QM/MM

Aiso T1 T2 T3

Hd2 H1 1 H1 2

Cys89

C-cap QM

Cu

N1 2

His94

H-link QM

32 38

33 40

36 38

28 28

22 27

Exp.c

1.53/1.78 0.70/1.09 0.62

1.77/1.51 1.40/1.01 13 16

The peptide chain of Cys89 is elongated by the atoms of neighboring aminoacids. This work. NMR results taken from Ref. [9].

electrons that have to be explicitly treated in the QM subsystem). Potentially, the capping carbon approach also provides a better description of spin delocalization by retaining carbon atom character at the domain boundary. The qualitative difference between the protein environment of the copper center in azurin and stellacyanin can be seen clearly from the simplified ribbon structures, shown in Figs. 1 and 2. In azurin, the aminoacids directly coordinated to the copper center are well screened by the protein backbone. The only exception is the d2 and 12 positions of His117, which are exposed to the solvent. On the other hand, in stellacyanin, parts of His46 and His94 opposite

from the copper center are easily accessible for the solvent molecules. These characteristic structural features lead to quite different trends in the QM/ MM results for the hyperfine structure and Mulliken spin population. For the copper nucleus, the anisotropic hyperfine tensors of azurin are greatly affected by the MM partial charges (changes of 8 – 32 MHz), but the isotropic hyperfine constant is not (change of 1 MHz). On the other hand, both the dipolar and isotropic constants of stellacyanin copper are greatly affected by MM partial charges (changes of 5 –18 and 17 MHz for the dipolar and isotropic constants, respectively). As can be seen from Table 3,

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

293

Table 3 Comparison of mulliken spin populations for azurin and stellacyanin models, with different approaches for the dangling bond capping Azurin, stellacyanin

Azurin

H-link QM C-cap QM QM/MM

Stellacyanin

H-link QM C-cap QM QM/MM

a b

Cys112, Cys89 Cu S

Hb1

Hb2

His46, His46 Nd1

His117, His94 Nd1

0.236a 20.016b 0.230 20.014 0.267 20.013

0.661 0.020 0.666 0.021 0.608 0.019

0.014 0.013 0.018 0.017 0.016 0.016

0.026 0.024 0.029 0.027 0.024 0.023

0.031 0.008 0.033 0.008 0.047 0.009

0.037 0.009 0.037 0.009 0.044 0.011

0.270 20.008 0.264 20.009 0.278 20.006

0.602 0.015 0.611 0.016 0.583 0.013

0.012 0.011 0.012 0.012 0.012 0.012

0.017 0.016 0.019 0.018 0.017 0.017

0.042 0.011 0.044 0.011 0.044 0.009

0.057 0.012 0.057 0.013 0.062 0.015

Total spin populations (s-, p- and d-type orbitals). Spin populations of s-type orbitals.

the unpaired electron spin populations (p- and dtypes) on the copper atom in azurin and stellacyanin are increased by 0.038 and 0.017, respectively, under the influence of MM charges. This illustrates the stronger influence of the MM charges on the anisotropic tensors of the copper in azurin. (We note that dipolar tensors are mainly determined by the interaction between the nuclear spin and anisotropic spin-density distribution, provided by the valence p and d-orbitals on copper). At the same time, the isotropic results indicate that, in azurin, s-type populations on copper are not affected by MM charges, in contrast to stellacyanin. In order to test the convergence of the QM/MM calculations with the QM subsystem size for the copper center in stellacyanin, we also examined an alternative QM region for the QM/MM calculation in this system (M-QM/MM). In this model (53 QM atoms), some atoms of neighboring residues attached to Cys89 were transferred from the MM part into the QM region. Because Cys89 sulfur interacts strongly with the copper (Table 3), accurate handling of Cys89 may be expected to be more important, compared to His46 and His95. From the results (see M-QM/MM in Table 2), it is seen that the hyperfine structures of the copper are not affected much by the elongation of Cys89 chain. By inference, extension of the QM subsystem to the direct neighbours of the histidine ligands, is likely to be even less important.

For the histidine nitrogen, in azurin, our QM/MM results for the hyperfine structure are in a good agreement with the previously published results for larger, pure QM models, and the experimental data. Nitrogen atoms, coordinated directly to the copper center (in d1 position), carry larger spin populations compared to the remote 12 nitrogen atoms. As a result, calculated isotropic hyperfine parameters of dnitrogen atoms (His46: 15.6 MHz; His117: 18.5 MHz in C-Cap QM; 18.7 and 22.0 MHz in QM/MM; vs. 18.1 and 25.1 MHz in experiment) are somewhat more sensitive to MM charges, than those of 1nitrogen atoms (His46: 0.56 MHz; His117: 0.93 MHz in C-Cap QM; 0.62 and 1.21 MHz in QM/MM; vs. 0.87 and 1.30 MHz in experiment). The absolute difference (3.1 and 3.5 MHz in d1 position; 0.06 and 0.28 MHz in 12 position) of Aiso values between the pure QM (C-Cap QM) and QM/MM models indicates relatively strong effects from the MM subsystem. In stellacyanin, the QM/MM calculations give nitrogen hyperfine parameters similar to the QM results (L-QM, this work) for larger models. As can be seen from Table 2 and Fig. 2, the effects from the protein backbone and solvent molecules are small, compared to azurin. At the same time, the results from the extended QM/MM model (M-QM/MM) are closer to the results from the large QM region (L-QM), and to the experiment. This indicates that additional shortrange interactions from neighbors of Cys89

294

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

(e.g. induction, exchange repulsion, charge transfer), completely neglected in standard QM/MM models, should be treated explicitly. For the cysteine b-hydrogen, in azurin, the QM/MM Aiso values (b1: 30 MHz; b2: 54 MHz) represent a significant improvement relative to the pure QM (C-Cap QM) results (34 and 64 MHz). The values are close to the previously published QM results [10,11] for larger models (28 and 46 MHz in LQM1; 28 and 32 MHz in L-QM2), and to the experiment (27 and 28 MHz). On contrary, in stellacyanin, there is no significant difference between the pure QM (C-Cap) (b1: 33 MHz; b2: 40 MHz) and QM/MM (36 and 38 MHz) results and the QM/MM results are far from the QM values for a larger model (b1: 22 MHz; b2: 27 MHz), and from experiment (13 and 16 MHz). In both cases, there are still significant discrepancies between the QM/MM and experimental values for b-hydrogens. Again, these discrepancies arise from the omission of short-range interactions of the cystein residue with the protein backbone. Treating these interactions explicitly within the DFT region improves the results considerably (M-QM/ MM—b1: 28 MHz; b2: 28 MHz in stellacyanin). For the histidine hydrogen, the QM/MM results are in better agreement with experiment than the pure QM (C-Cap) treatment of the radical center. In azurin, QM/MM Aiso values represent a significant change (0.11 – 0.34 MHz) relative to the pure QM (C-Cap) results, except 12 hydrogen of His46. Following the trend seen for other hyperfine parameters, in stellacyanin, the QM/MM results deviate much less from the pure QM results (by 0.02 –0.13 MHz). In azurin, the hyperfine parameters of histidine hydrogens should be affected mainly by the protein backbone (d2: 0.24 MHz; 11: 0.34 MHz in His46, 11: 0.34 MHz in His117), while those of the d2 and 12 hydrogens in His117 should be affected mainly by solvent molecules (d2: 0.21 MHz; 12: 0.11 MHz). In stellacyanin, an influence of the solvent may be expected for most hydrogens, but the effect (0.02 – 0.13 MHz) is less pronounced than in azurin.

models of paramagnetic sites in biological systems. The technique employs empirically parameterized carbon pseudopotentials [13] for capping bonds at the domain boundary. Long-range interactions between the QM and MM subsystems are treated by adding an electrostatic term, due to the MM charges, to the oneelectron part of the QM Hamiltonian. The approach is applied to the calculation of hyperfine parameters of two blue copper proteins (azurin and stellacyanin). The new QM/MM approach shows two significant advantages: Firstly, it can be implemented easily without code modification. Capping potentials can be treated as conventional ECPs, and point charges can be handled for energy calculations by most DFT programs. Secondly, hyperfine parameters of paramagnetic proteins are well described by a small QM model surrounded by MM partial charges. This is in part due to the use of the capping potentials, which reduce disturbances of the electronic structure of the QM subsystem, arising from broken bonds at the domain boundary. At the same time, care should be taken not to place the domain boundary too close to the nuclei of interest. In the present study, b-hydrogens of the cysteine residue (directly coordinated to the copper center) were placed too close to the domain boundary used in our initial QM/MM model. Adequate treatment of such centers requires that additional short-range interactions at the boundary (induction, exchange repulsion, and charge transfer) should be included. In this work, the desired effects were achieved by an elongation of the polypeptide chain at the cysteine site. In summary, our hybrid QM/MM method provides computationally efficient access to the hyperfine parameters of the interesting local (QM) part in a protein, while at the same time including the essential effects from its surroundings. This approach can be extended to the calculation of other magnetic properties.

Acknowledgements 5. Conclusion We developed a new QM/MM approach for the calculations of hyperfine parameters using realistic

The authors would like to thank Gino A. DiLabio for providing carbon capping potentials and for useful discussions of QM/MM theories. We are grateful to NSERC and to CIPI for financial support.

S. Moon et al. / Journal of Molecular Structure (Theochem) 632 (2003) 287–295

References [1] J.W.A. Coremans, O.G. Poluektov, E.J.J. Groenen, G.W. Canters, H. Nar, A. Messerschmidt, J. Am. Chem. Soc. 116 (1994) 3097. [2] E.I. Solomon, M.J. Baldwin, M.D. Lowery, Chem. Rev. 92 (1992) 521. [3] J.E. Roberts, T.G. Brown, B.M. Hoffman, J. Peisach, J. Am. Chem. Soc. 102 (1980) 825. [4] J.E. Roberts, J.F. Cline, V. Lum, H. Freeman, H.B. Gray, J. Peisach, B. Reinhammar, B.M. Hoffman, J. Am. Chem. Soc. 106 (1984) 5324. [5] M.M. Werst, C.E. Davoust, B.M. Hoffman, J. Am. Chem. Soc. 113 (1991) 1533. [6] J.W.A. Coreman, O.G. Poluektov, E.J.J. Groenen, G.W. Canters, H. Nar, A. Messerschmidt, J. Am. Chem. Soc. 118 (1996) 12141. [7] J.W.A. Coremans, O.G. Poluektov, E.J.J. Groenen, G.W. Canters, H. Nar, A. Messerschmidt, J. Am. Chem. Soc. 119 (1997) 4726. [8] I. Bertini, S. Ciurli, A. Dikiy, R. Gasanov, C. Luchinat, G. Martini, N. Safarov, J. Am. Chem. Soc. 121 (1999) 2037. [9] I. Bertini, C.O. Fernandez, B.G. Karlsson, J. Leckner, C. Luchinat, B.G. Malmstrom, A.M. Nersissian, R. Pierattelli, E. Shipp, J.S. Valentine, A.J. Vila, J. Am. Chem. Soc. 122 (2000) 3701. [10] A.R. Jaszewski, J. Jezierska, Chem. Phys. Lett. 343 (2001) 571. [11] M. Swart. Ph.D Thesis, Materials Science Center, Rijksuniversiteit Groningen (RuG), Nijenborgh 4, 2002. [12] M. Munzarova, M. Kaupp, J. Phys. Chem. A 103 (1999) 9966.

295

[13] G.A. DiLabio, M.M. Hurley, P.A. Christiansen, J. Chem. Phys. 116 (2002) 9578. [14] M.S. Gordon, M.A. Freitag, P.B. Banyopadhyay, J.H. Jensen, V. Kairys, W.J. Stevens, J. Phys. Chem. A 105 (2001) 293. [15] Q. Cui, M. Karplus, J. Phys. Chem. B 104 (2000) 3721. [16] R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, 16, Oxford University Press, Oxford, 1989. [17] D.L. Beveridge, in: G.A. Segal (Ed.), Semiempirical Methods of Electronic Structure Calculation, Plenum Press, New York, 1977, Chapter 5. [18] V.G. Malkin, O.L. Malkina, L.A. Eriksson, D.R. Salahub, in: J.M. Seminario, P. Polizer (Eds.), Modern Density Functional Theory: A Tool for Chemistry, Theoretical and Computation Chemistry, vol. 2, Elsevier, Amsterdam, 1995, pp. 9, Chapter 9. [19] A.M. Lee, N.C. Handy, S.M. Colwell, J. Chem. Phys. 103 (1995) 10095. [20] H. Fukui, Magn. Res. Rev. 11 (1985) 205. [21] P. Pyykko, Adv. Quantum Chem. 11 (1978) 354. [22] A. St-Amant, D.R. Salahub, Chem. Phys. Lett. 169 (1990) 387. [23] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, P. Kollman, J. Am. Chem. Soc. 117 (1995) 5179. [24] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [25] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [26] W. Kutzelnigg, U. Fleischer, M. Schindler, NMR-Basic Principles and Progress, vol. 23, Springer-Verlag, Heidelberg, 1990. [27] The X-leap program is distributed with the AMBER (version 6.0) program. (http://www.amber.ucsf.edu/amber/amber. html) [28] R. McWeeny, Methods of Molecular Quantum Mechanics, 2nd ed., Academic Press, London, 1989, Chapter 11.