Burda
email: burda@karlov.mff.cuni.cz

Matěj Pavelka

Jaroslav V Burda

Computational Study of Redox Active Centers of Blue Copper Proteins: A Computational DFT Study

Keywords: DFT calculations, plastocyanin, blue copper proteins, copper complexes

Introduction

Every sixth entry in the protein data bank (PDB) contains a metal cofactor. Such statistics demonstrate the importance of the metal ions for catalytic processes in a living cell. Hence, investigation of metal sites -such as copper redox centers -is essential for determination of molecular mechanisms in biochemistry. The Cu ions often exhibit interesting spectral properties, which originate from the unusual geometric and electronic structures that are imposed through their interactions with protein environment. Such proteins provide many functions: electron transfer, oxidation-reduction processes, oxygen transport and insertion, and so forth.

Blue copper proteins forms the so-called Type 1 Cu proteins [START_REF] Cowan | Inorganic Biochemistry: An Introduction[END_REF][START_REF]ProMiS: Prosthetic groups and Metal Ions in Protein Active Sites[END_REF]. These Cu proteins have an intense absorption band near 600 nm in the oxidized Cu(II) state. This transition is assigned with the S(cysteine)-Cu Ligand-to-Metal-Charge-Transfer (LMCT). Structures of the active sites usually contain a 4-coordinated Cu ion, though a coordination number of five was found in azurins. Structural motif consists of the arrangement Cu(I)/Cu(II):(His)2CysX, where X is Met or Gln. Examples of this motif can be found in pseudoazurin, rusticyanin, plastocyanin, mavicyanin, auracyanin, stellacyanin, umecyanin, and amicyanin. Interestingly, reduction of the Cu(II) ion to Cu(I) in type 1 proteins causes minimal structural changes. It results in a low activation barrier for the Cu(I)→Cu(II) redox process and thus a rapid electron Solomon [START_REF] Solomon | Spectroscopic Methods in Bioinorganic Chemistry: Blue to Green to Red Copper Sites[END_REF].

There is a large number of works investigating biological activity of the copper proteins using both experimental and theoretical approaches. By means of UV-VIS and EPR spectroscopy, Cu-center of azurin was studied [4]. Spectroscopic tools in combination with DFT calculations were used to investigate the role of amino acid in axial position to the copper complex and its influence on a reduction potential [START_REF] Palmer | Spectroscopic Studies and Electronic Structure Description of the High Potential Type 1 Copper Site in Fungal Laccase: Insight into Effect of the Axial Ligand[END_REF].

Published experimental works include a vast range of techniques, e.g. fluorescence spectroscopy [START_REF] Santra | Novel Interaction between Glutamate and the Cu2+ /DMABN/beta-CD Complex[END_REF], X-ray absorption near-edge structure (XANES) spectra [START_REF] Shimizu | Ligand field effect on the chemical shift in XANES spectra of Cu(II) compounds[END_REF] [START_REF] Barrett | Atomic Resolution Crystal Structures, EXAFS, and Quantum Chemical Studies of Rusticyanin and Its Two Mutants Provide Insight into Its Unusual Properties[END_REF], combination of EPR and electron-nuclear double resonance (ENDOR) techniques [START_REF] Manikandan | Structure of Copper(II)-Histidine Based Complexes in Frozen Aqueous Solutions As Determined from High-Field Pulsed Electron Nuclear Double Resonance[END_REF], EPR and UV-VIS spectra [START_REF] Sigman | Design and Creation of Cu(II)-Binding Site in Cytochrome c Peroxidase that Mimics the Cub-heme Center in Terminal Oxidases[END_REF], and many other. The basic aspects of a copper coordination in the protein environment are summarized in reviews [START_REF] Gray | Copper coordination in blue proteins[END_REF][START_REF] Solomon | Electronic Structures of Active Sites in Copper Proteins: Contribution to Reactivity[END_REF][START_REF] Solomon | Electronic Structures of Metal Sites in Proteins and Models: Contribution to Function in Blue Copper Proteins[END_REF][START_REF] Mirica | Stack: Structure and Spectroscopy of Copper-Dioxygen Complexes[END_REF]. Charge transfer (CT) dynamics of the blue copper proteins was investigated using a pumpprobe [START_REF] Book | Ultrafast Pump-Probe Studies of Excited-State Charge Transfer Dynamics in Blue Copper Proteins[END_REF] and resonance Raman spectroscopy [START_REF] Fraga | Charge-Transfer Dynamics in Plastocyanin, a Blue Copper Protein, from Resonance Raman Intensities[END_REF]. Some other experimental studies dealing with these proteins should be also mentioned [START_REF] Holland | Heterocumulene Metathesis by Iridium Guanidinate and Ureylene Complexes: Catalysis Involving Reversible Insertion To Form Six-Membered Metallacycles[END_REF][START_REF] Holland | Three-Coordinate Cu(II) Complexes: Structural Models of Trigonal-Planar Type 1 Copper Protein Active Sites[END_REF][START_REF] Taylor | Modelling the impact of geometric parameters on the redox potential of blue copper proteins[END_REF]. Theoretical studies of the copper interactions with amino acids have been reported recently [START_REF] Raffa | Ab initio model studies of copper binding to peptides containing a His-His sequence: relevance to the beta-amyloid peptide of Alzheimer's disease[END_REF][START_REF] Sabolovic | Why are Cu(II) amino acid complexes not planar in their crystal structures? An ab initio and molecular mechanics study[END_REF][START_REF] Sabolovic | Modeling anhydrous and aqua copper(II) amino acid complexes: A new molecular mechanics force field parametrization based on quantum chemical studies and experimental crystal data[END_REF][START_REF] Katz | Copper-binding motifs: Structural and theoretical aspects[END_REF][START_REF] Rulíšek | Theoretical studies of the interactions of transition metals with biomolecules[END_REF][25][START_REF] Rulíšek | Theoretical Studies of Metal Ion Selectivity. 3. A Theoretical Design of the Most Specific Metal-Binding Sites of Metalloproteins for Selected Transition Metal Ions (Co2+, Ni2+, Cu2+, Zn2+, Cd2+, and Hg2+)[END_REF][START_REF] Bertran | The Different Nature of Bonding in Cu+ -Glycine and Cu2+ -Glycine[END_REF][START_REF] Shoeib | A Study of Silver (I) Ion-Organonitrile Complexes: Ion Structures, Binding Energies, and Substituent Effects[END_REF][START_REF] Caraiman | Investigations of the gas-phase reactivity of Cu+ and Ag+ glycine complexes towards CO, D2O and NH3[END_REF][START_REF] Breza | On the reduced form of (imidazole-N-3)(Nsalicylidene-alaninato-O,N,O ') copper(II)[END_REF]. A lot of computational effort was devoted to the examination of copper proteins [START_REF] Shuku | Molecular orbital analysis of active site of oxidized azurin: Dependency of electronic properties on molecular structure[END_REF][START_REF] Olsson | Frozen density functional free energy simulations of redox proteins: Computational studies of the reduction potential of plastocyanin and rusticyanin[END_REF][START_REF] Olsson | The influence of axial ligands on the reduction potential of blue copper proteins[END_REF][START_REF] Olsson | On the relative stability of tetragonal and trigonal Cu(II) complexes with relevance to the blue copper proteins[END_REF][START_REF] Corni | Role of the electronic properties of azurin active site in the electron-transfer process[END_REF][START_REF] Warshel | Computer Modeling of Chemical Reactions in Enzymes and Solutions[END_REF][START_REF] Prabhakar | Theoretical study of the mechanism for the oxidative half-reaction of copper amine oxidase (CAO)[END_REF][START_REF] Yoshizawa | Conversion of methane to methanol at the mononuclear and dinuclear copper sites of particulate methane monooxygenase (pMMO): A DFT and QM/MM study[END_REF][START_REF] Van Den Bosch | Calculation of the redox potential of the protein azurin and some mutants[END_REF][START_REF] Siegbahn | Advances in Chem. Physics[END_REF]. The first theoretical spectra of plastocyanin were computed by Solomon's group [START_REF] Penfield | [END_REF]42] and many studies were reported recently, e.g. [START_REF] Sinnecker | QM/MM calculations with DFT for taking into account protein effects on the EPR and optical spectra of metalloproteins. Plastocyanin as a case study[END_REF][START_REF] Lacroix | Electronic structure of the perturbed blue copper site in nitrite reductase: Spectroscopic properties, bonding, and implications for the entatic/rack state[END_REF][45][START_REF] Randall | Spectroscopic and Electronic Structural Studies of Blue Copper Model Complexes. 1. Perturbation of the Thiolate-Cu Bond[END_REF][START_REF] Lacroix | Spectroscopic and geometric variations in perturbed blue copper centers: Electronic structures of stellacyanin and cucumber basic protein[END_REF][START_REF] Randall | Electronic structure contributions to electron transfer in blue Cu and Cua[END_REF][START_REF] Guckert | ELECTRONIC-STRUCTURE OF THE REDUCED BLUE COPPER ACTIVE-SITE -CONTRIBUTIONS TO REDUCTION POTENTIALS AND GEOMETRY[END_REF]. Calculation of the EPR parameters were published for aqua-ligated Cu(II) complexes at the modified B3LYP level [START_REF] De Almeida | Modeling of EPR parameters of copper(II) aqua complexes[END_REF]. Similarly also azurin and some other blue copper proteins were explored [START_REF] Remenyi | Density Functional Study of EPR Parameters and Spin-Density Distribution of Azurin and Other Blue Copper Proteins[END_REF]. In the case of plastocyanin, QM/MM method was successfully applied [START_REF] Sinnecker | QM/MM calculations with DFT for taking into account protein effects on the EPR and optical spectra of metalloproteins. Plastocyanin as a case study[END_REF]. Force-field potentials were developed for Cu-ligand interactions in plastocyanin [START_REF] Ando | Excited State Potentials and Ligand Force Field of a Blue Copper Protein Plastocyanin[END_REF] and applied in the study of LMCT dynamics for this protein [START_REF] Ando | Ligand-to-Metal Charge-Transfer Dynamics in a Blue Copper Protein Plastocyanin: A Molecular Dynamics Study[END_REF]. An analogous attempt to construct an empirical force-field for the oxidized form of blue [START_REF] Comba | A New Molecular Mechanics Force Field for the Oxidized Form of Blue Copper Proteins[END_REF]. Same authors studied coordination bonds in amicyanin and rusticyanin using QM/MM (IMOMM) method [START_REF] Comba | Hybrid quantum mechanics/molecular mechanics studies of the active site of the blue copper proteins amicyanin and rusticyanin[END_REF].

A great attention was also paid to examination of the "small inorganic complexes" in order to determine coordination geometries and electronic properties of various copper compounds. Investigation of these small models can provide a deeper insight into and an easier interpretation of the Cu(I)/Cu(II) behavior. Interactions of the both Cu cations with molecules like water, ammonia, or hydrogen-sulfide were intensively studied by many researchers . A more detailed literature discussion on this topics can be found in our previous papers where these complexes were inspected [START_REF] Burda | Theoretical model of copper Cu(I)/Cu(II) hydration. DFT and ab initio quantum chemical study[END_REF][START_REF] Pavelka | Theoretical model of copper Cu(I)/Cu(II) mixed aqua-amine complexes. DFT and ab initio quantum chemical study[END_REF][START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF].

This study focuses on models of active mononuclear centers, which are present in blue copper proteins (e.g. plastocyanin in Figure 1). This study gives a comparison of electronic properties for general tetra-coordinated complexes with the active centers of blue copper proteins (within the fixed-geometry from pdb structures).

Computational Details

Active sites of the blue copper proteins (Scheme 1) were studied at the DFT level of theory. In the Type A family of the blue proteins (with Met residue), the following protein structures from the PDB database [START_REF] Berman | The Protein Data Bank[END_REF] were considered: amicyanin (1AAC), auracyanin (1QHQ), plastocyanin (1KDI and 1KDJ), and rusticyanin (1A3Z and 1RCY). In the case of the Type B centers (with Gln residue), the following proteins were analyzed: mavicyanin (1WS8), stellacyanin (1JER), and umecyanin (1X9R and 1X9U). These X-ray crystal structures were obtained with resolution in range from During the second optimization step (labeled "Constrain optimization phase II"), the copper cations could move freely (Scheme 2b) in the fixed framework of frozen ligands (resulting from the previous optimization phase). In the last step, the full optimization of the whole complex was allowed leading to the reduction of all the Type A centers to one final structure and analogous reduction also occurred in the case of the B centers. Since the second optimization step did not bring any important changes in coordination, it is not considered in further text.

Quantum chemical calculations were performed at the DFT level using the B3LYP functional. For the H, C, O, N atoms, the 6-31+G(d) basis set was applied.

The copper core electrons were described by the Christiansen Averaged Relativistic

Effective Pseudopotentials (AREP) [START_REF] Hurley | Ab Initio Relativistic Effective Potentials with Spin-Orbit Operators[END_REF]. A consistent basis set was adopted for the Reduced centers with the Cu + cation were represented by a closed shell singlet electronic ground state, while the oxidized Cu(II) complexes, where 3d 9 electron configuration is present, possess doublet ground states. Hence, some attention had to be devoted to the correct initial guess for the SCF procedure and the appropriate character of the wave function was checked in the end of every calculation.

All optimizations were preformed in vacuo as well as with inclusion of the solvation effects. The Conductor-like screening model (COSMO) [82,[START_REF] Klamt | Cosmo -a New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient[END_REF] energies. In these calculations, Basis Set Superposition Error (BSSE) corrections and corrections on the deformation energies [START_REF] Boys | The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors[END_REF] were considered according to the equation eq. ( 1)

) ( Stab BSSE complex fragment deform E E E E - ∆ = - - ∑ ∑
where E complex represents the total energy of the whole complex and leads to incorrect values of these energies. The correct dispersion and repulsion terms were taken from the single-point calculations of the corresponding deformation corrections (for further details see ref. [START_REF] Zimmermann | Cisplatin Interaction with Cysteine and Methionine in Aqueous Solution: Computational DFT/PCM Study[END_REF]).

In addition, total coordination energies (∆ ∆ ∆ ∆E TCE ) were also determined. In this case only two fragments were considered in eq. ( 1): a) all the interacting ligand molecules treated as one subsystem and b) the central Cu ion as another one. These values can be understood as the binding energies of the cation with a pre-formed ligand shell. In our previous work [START_REF] Burda | Theoretical model of copper Cu(I)/Cu(II) hydration. DFT and ab initio quantum chemical study[END_REF][START_REF] Pavelka | Theoretical model of copper Cu(I)/Cu(II) mixed aqua-amine complexes. DFT and ab initio quantum chemical study[END_REF][START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF], the ∆ ∆ ∆ ∆E TCE energies were labeled ∆ ∆ ∆ ∆E Stex .

The difference between ∆ ∆ ∆ ∆E Stab and ∆ ∆ ∆ ∆E TCE estimates the energy, which must be invested to form the ligand shell arrangement in the absence of the ion. However, the actual interligand repulsion is larger in the presence of the cation due to the additional polarization effects [START_REF] Šponer | Stabilization of the Purine.Purine.Pyrimidine DNA Base Triplets by Divalent Metal Cations[END_REF]. The ∆ ∆ ∆ ∆E BE energies of the Cu-L ligand bonds were obtained using eq. ( 1)

where partitioning according to Cu-L bond must be used. The ∆ ∆ ∆ ∆E BE values were considered without the contributions of the deformation corrections.

In order to describe redox activity of the studied complexes, vertical and adiabatic ionization potentials (IP) were calculated for the reduced Cu(I) centers according to formula:

(2)

( ) ( ) Cu II Cu I IP E E = - .
In the case of the vertical IP, the E Cu(II) term represents the energy of a (2+) charged system calculated in the Cu(I) optimized structure. For the adiabatic IP, the E Cu(II)

energy was determined from the Cu(II) optimized structure.

The charge and spin distributions were analyzed in terms of partial atomic charges and spin densities were obtained by the Natural Population Analyses (NPA) [START_REF] Reed | Natural population analysis[END_REF]. Spin isodensities (ρ s = 0.01e/Å 3 ) were plotted together with selected Molecular Orbitals (MO).

Main axes of a diagonalized g-tensor were estimated in order to investigate the EPR spectrum and consequently the behavior of the unpaired electron for all studied Cu(II) complexes.

All the quantum chemical calculations were performed using Gaussian 03

program package [START_REF] Frisch | Self-Consistent Molecular Orbital Methods 25. Supplementary Functions for Gaussian Basis Sets[END_REF]. NBO v. 5.0 program from Wisconsin University [START_REF] Weinhold | NBO 5.0 Program University of Wisconsin[END_REF] was used to evaluate the Natural Bond Orbital (NBO) characteristics and AIMPAC program was used Bader's AIM [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF] (Atoms in Molecules) analyses . For visualization of geometries, spin densities, and molecular orbitals, graphical programs Molden v.

4.4 [START_REF] Schaftenaar | Molden: a pre-and post-processing program for molecular and electronic structures[END_REF] and Molekel v. 4.3 [START_REF] Flükiger | Development of the molecular graphics package MOLEKEL and its application to selected problems in organic and organometallic chemistry[END_REF][START_REF] Portmann | MOLEKEL: An Interactive Molecular Graphics Tool[END_REF] were used. of blue copper proteins were optimized in consequent steps with selected coordinates frozen in order to enlighten the role of the protein constraints.

Results and Discussion

Protein structures -Constrained optimizations

The averaged coordination distances for the oxidized Cu(II) centers of the both Types are compiled in Table 1. Bond lengths obtained from the constrained optimizations correspond directly to the PDB distances from the X-ray structures.

In average, the Cu(II)-S(cysteine) bonds are about 2.2 Å long. For a comparison, the Cu-

X (S, O, N) distances (in Å) for [Cu(H 2 S) m (H 2 O) n (NH 3 ) k ] 2+/+
cations from our previous study [START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF] are presented in Table 2 together with values obtained from the CSD database by Katz et al. [START_REF] Katz | Copper-binding motifs: Structural and theoretical aspects[END_REF]. Significantly longer Cu(II)-S distances (about 2.4 Å) were found in the 4-coordinated complexes with electroneutral H 2 S ligands. In the models of blue copper proteins, the shorter Cu-S(Cys) bond can be explained by a stronger enhancement of the dative bond by electrostatic interaction between the Cu 2+ cation and the present Cys model (thiolate [S(CH 3 )] 1-).

In the case of Met model (thioether group S(CH 3 ) 2 ), only a weaker monopole-dipole (µ = 1.6 D) interaction is present. This fact correlates with the visibly longer Cu-S(Met) bond (ca 2.89 Å). Glutamine exhibits coordination distances of about 2.15 Å.

It is caused by smaller radii of the oxygen atom (in comparison with more diffuse sulfur) and by a larger dipole moment (µ = 4.0 D) of the acetyl-amine group.

Corresponding bond lengths of the Cu(II) aqua complexes from our previous study [START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF] were determined to amount in average to 2.01 Å. The difference can be partially explained by the forced geometry of the protein centers, where the axial position of the 4 th ligand is disfavored (cf. below). When oxidation state is altered, the redox centers of blue copper proteins exhibit small changes in the copper-ligand distances, meaning that qualitatively the same geometry is conserved. This can be considered to be a consequence of the protein "matrix". For the previously studied [START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF] 4-coordinated inorganic copper complexes, transition from the Cu(II) to Cu(I) states is connected with pronounced structural changes from square-planar to tetrahedral ligand arrangements. When such systems are allowed to be even more relaxed, they transform to 2-coordinated structures as illustrated in Figure 3. This fact is important especially for the electron transferring proteins. Small structural changes of the active centers -basically changes of the Cu-L(Met/Gln) distances -minimize the reorganization energy λ λ λ λ within the redox process, and thus enable a rapid electron transfer.

Fully optimized structures in vacuo

When the models of reduced centers are optimized without any restrictions in vacuo, interesting changes can be noticed (see Figures 4a and4b). The Cu atom, which lies above S, N, N-plane in non-optimized Cu(I) structures, moves so that it is situated within this plane. complexes [START_REF] Pavelka | Copper Cation Interactions with Biologically Essential Types of Ligands: A Computational DFT Study[END_REF].

The full optimization has marginal effect on the Cu-S(Cys) distances in the case of the A centers. However in the B centers, the Cu-S(Cys) bonds elongate. Such behavior can be considered to be a result of the ligand-donation competition between the carbonyl and thiolate groups. The donation effect from a more distant methionine was found negligible in the A centers. In the B centers, the full optimization leads to a substantially different coordination arrangement (see Figure 4d). The relaxed complex possesses a geometry halfway between square planar and tetrahedral arrangements, which has some impact on electronic properties as will be demonstrated later.

A lot of studies was devoted to investigation of protein constraints on the active centers of blue copper proteins. An earlier paper [94] suggested that there are no constraints. However, such a proposal was later revised [START_REF] Olsson | The influence of axial ligands on the reduction potential of blue copper proteins[END_REF][START_REF] Olsson | On the relative stability of tetragonal and trigonal Cu(II) complexes with relevance to the blue copper proteins[END_REF]45]. The role of axial ligand has also been widely studied [START_REF] Palmer | Spectroscopic Studies and Electronic Structure Description of the High Potential Type 1 Copper Site in Fungal Laccase: Insight into Effect of the Axial Ligand[END_REF][START_REF] Olsson | The influence of axial ligands on the reduction potential of blue copper proteins[END_REF][START_REF] Randall | Spectroscopic and Electronic Structural Studies of Blue Copper Model Complexes. 1. Perturbation of the Thiolate-Cu Bond[END_REF][START_REF] Lacroix | Spectroscopic and geometric variations in perturbed blue copper centers: Electronic structures of stellacyanin and cucumber basic protein[END_REF][START_REF] Guckert | ELECTRONIC-STRUCTURE OF THE REDUCED BLUE COPPER ACTIVE-SITE -CONTRIBUTIONS TO REDUCTION POTENTIALS AND GEOMETRY[END_REF][START_REF] Randall | Spectroscopic and Electronic Structural Studies of Blue Copper Model Complexes. 2. Comparison of Three-and Four-Coordinate Cu(II)-Thiolate Complexes and Fungal Laccase[END_REF][START_REF] Solomon | Electronic structure of the oxidized and reduced blue copper sites: Contributions to the electron transfer pathway, reduction potential, and geometry[END_REF]. The section above provides some answers to the questions about constraints imposed by the protein backbone

Relaxed models in solvents

The studied redox centers are influenced by their environment. Such effects were examined at the QM/MM level [START_REF] Olsson | Frozen density functional free energy simulations of redox proteins: Computational studies of the reduction potential of plastocyanin and rusticyanin[END_REF][START_REF] Sinnecker | QM/MM calculations with DFT for taking into account protein effects on the EPR and optical spectra of metalloproteins. Plastocyanin as a case study[END_REF][START_REF] Datta | Theoretical determination of the standard reduction potential of plastocyanin in vitro[END_REF]. Polarized continuum models were also used for this kind of complexes in ref. [START_REF] Sinnecker | QM/MM calculations with DFT for taking into account protein effects on the EPR and optical spectra of metalloproteins. Plastocyanin as a case study[END_REF]. Given the model used in our calculations, the subtler differences within the individual families of (both Type A and B) proteins cannot be reproduced. However, the influence of solvent effects can be still estimated in the PCM regime through the screening of electrostatic interactions. These effects can be important for more accurate energy analyses.

In case of the oxidized Cu(II) centers, solvation effects exhibit a lower influence on the fully optimized geometries than in the reduced Cu(I) complexes.

Screened electrostatic interaction of the thiolate Cu-S(Cys) bond results in a mild elongation with an increasing value of the dielectric constant. Since the donation effect is very distance-sensitive (it is linked to the overlap of Cu and S AOs), the cysteine donation weakens correspondingly. Simultaneously, shortening of the Cu-N(His) bonds occurs increasing the donation competition. Interestingly, the distance of the Cu-O(Gln) and Cu-S(Met) bonds also shorten. Such a behavior causes a dative competition between the Cys and all other (neutral) ligands. In the Gln case, the strengthening of the Cu-O(Gln) bond is in accord with the hard-soft-acid-base (HSAB) principle (oxygen is closer to the relatively hard Cu 2+ cation than the sulfur atom [98]).

Energetics

Energy analyses and determination of electron properties in next the section were performed at the B3LYP/6-311++G(2df,2pd) level. 3 and Figure 5).

Stabilization energies

Ligand bonding energies

Bonding energies ∆ ∆ ∆ ∆E BE of all the ligands in the studied redox centers are compiled in densities together with the EPR spectra estimation.

Charge and spin density distribution

The electron spin densities in terms of partial atomic spin densities (ρ ρ ρ ρ s obtained by the NPA method) are summarized in Table 7. Since there is a strong dative interaction between Cu 2+ and cysteine model [S(CH 3 )] -, most of the spin density (80 -90 %) is located on the Cu-S(Cys) bond as it can be observed in

Molecular orbitals analysis

In this part, only behavior of single occupied molecular orbitals (SOMO) in the oxidized centers will be discussed. The analysis revealed three types of SOMO, which are illustrated in Figure 7. The first type (7a) represents a typical SOMO occurring in the A centers, which was also presented in some other works [45,[START_REF] Lacroix | Spectroscopic and geometric variations in perturbed blue copper centers: Electronic structures of stellacyanin and cucumber basic protein[END_REF][START_REF] Guckert | ELECTRONIC-STRUCTURE OF THE REDUCED BLUE COPPER ACTIVE-SITE -CONTRIBUTIONS TO REDUCTION POTENTIALS AND GEOMETRY[END_REF]. In this type, the 3p π orbital of S(Cys) is combined with In the third type of molecular orbital (7c), SOMO is dominantly composed of copper 3d AO's so that the lobes aim to all four Cu-ligand bonds.

EPR spectrum of oxidized centers

EPR spectra represent an important tool for investigation of the Cu(II) complexes, since the method is based on magnetic properties of the unpaired electrons. The key tensor property called the g-factor describes how a local magnetic field is induced by the external magnetic field. It approaches the value of 2.0023 for a single free electron. Similar values should be obtained for systems with delocalized unpaired electron. Otherwise, it can markedly rise. Main axes of the diagonalized g-tensor were determined for all the studied Cu(II) complexes (Table 9). Experimental values are listed as well. In case of the protein imposed structures, g ║ and g ┴ values in vacuo were estimated to be about 2.12 and 2.05 a.u., respectively. This anisotropic splitting of the g-factor is close to the experimental values (2.23 and 2.05) for the A centers of the blue copper proteins [START_REF] Solomon | Electronic Structures of Metal Sites in Proteins and Models: Contribution to Function in Blue Copper Proteins[END_REF][START_REF] Solomon | Electronic structure of the oxidized and reduced blue copper sites: Contributions to the electron transfer pathway, reduction potential, and geometry[END_REF]. Considering the solvent effects in the COSMO regime slightly improves the calculated g ║ values. However, in order to obtain a better agreement with the experiment, the hybrid functional needs to be modified by fitting the ratio of the DFT exchange-correlation and the Hartree-Fock exchange terms as discussed e.g. in review [START_REF] Solomon | Electronic Structures of Metal Sites in Proteins and Models: Contribution to Function in Blue Copper Proteins[END_REF]. On the other hand, the B3LYP functional successfully reproduces higher g ║ values for the B centers, even in vacuo. The fully optimized A complex exhibits similar behavior as the protein structures from the database. In case of the B center, full optimization reveals a different trend, where g ║ is very similar to g ┴ = 2.08. Such behavior can be explained by pronounced geometry changes in the fully optimized B center with a more competitive Gln ligand as mentioned in the previous geometry discussion and demonstrated in Figure 4d. In this way, extensive redistribution of unpaired electron occurs as it was already illustrated by plots of spin densities in Figure 6.

For the solvated complexes, the distinction between g ║ and g ┴ increases with dielectric constant in the A centers as a reaction to the shortening of Cu-N bonds and concentrating electron density in N,N,S plane. On the contrary, the difference between both g components decreases in the B centers since a simultaneous increase of the Cu-S(Cys) and decrease of Cu-O(Gln) distances causes formation of a more symmetrical (isotropic) arrangement (cf. Table 1).

Conclusions

Models of mononuclear centers of blue copper proteins were investigated in the present study. Two families of these proteins were examined: the Type A centers with methionine side chain as the fourth residue and the Type B centers with glutamine residue.

  interesting review on these complexes was published recently by

  made earlier based on DFT calculations

2 .Scheme 1 : 2 :

 212 Scheme 1: The 4-coordinated Cu(I)/Cu(II) redox center of blue copper proteins. Residue R is either methionine or glutamine.

  valence electrons. Double-ζ pseudoorbitals of Cu were augmented by the diffuse and polarization functions (α s = 0.025, α p = 0.35, α d = 0.07, and α f = 3.75). Similarly, pseudoorbitals of the sulfur atom were extended by analogous functions with exponents: α s = 0.077, α p = 0.015, and α d = 0.50.

  Structures of both oxidized and reduced centers of the A and B Types (see Figure 2)

4 :

 4 Scheme 4a, ligands with oxygen and sulfur atoms usually do not have dipole moment

  that an asymmetrical arrangement of the two Cu-N(His) distances is apparent in the PDB structures.

Scheme 5 :

 5 Scheme 5a) lead to an assumption that reduced centers without protein constraints

  oxidized nor reduced centers have optimum geometries. The extent of these forces can be estimated from the stabilization energies (cf. Section 3.2.1). b) When no constraint is applied, the axial copper-ligand bond elongates/shortens depending on reduced/oxidized state of metal, regardless the type of 4 th residue (Gln or Met).

  in coordination in reduced centers correspond to the fact that the Cu(I) complexes prefer geometry with a lower-coordination. Analogously, models of the fully relaxed oxidized Cu(II) centers tend to form four "equivalent" bonds as can be seen from unification of the Cu-N(His) bonds and shortening of the Cu-L distances of the 4 th ligand.

16 For

 16 the solvated Cu(I) centers, the Cu-S(Cys) distances elongate and similarly Cu-O(Gln) reflecting the reduced electrostatic enhancement in (implicit) solvent.

  Cu-L bonding energy, which results in lower ∆ ∆ ∆ ∆E TCE energy for fully optimized system. Inclusion of the environment effects influences the stabilization energy as the screening decreases the electrostatic contribution to the bonding interaction between the Cu(I)/Cu(II) cation and ligands. Especially remarkable influence can be seen in the negatively charged model of cysteine. Therefore, the obtained ∆ ∆ ∆ ∆E Stab values are pronouncedly lowered when compared to in vacuo calculations. The dielectric constant, proportionally lowers the stabilization energies. An interesting behavior was found when relaxation of the protein constraints (estimated as a difference of the ∆ ∆ ∆ ∆E Stab energies of constrained and fully optimized structures) is considered. These values are drawn in Figure 5. This Figure illustrates several trends: a) In both types of the centers, the reduced Cu(I) complexes are less stable under protein constraints than the oxidized complexes. b) The B centers exhibit a higher energy relaxation after full optimization. c) The relaxation energies are most pronounced in the gas phase calculations, where stabilization increases after full optimization up to 21 kcal/mol (cf. Table

Figure 6 .

 6 The B centers exhibit a higher portion of spin density ρ ρ ρ ρ s on Cu and lower ρ ρ ρ ρ s on the sulfur atom of the cysteine model when compared to the A centers. It shows that in oxidized Type B centers, always lower donation occurs from cysteine sulfur than in Type A centers, leaving larger spin density on copper cation. This is linked with a stronger donation ability of the Gln residue (than Met). These analyses are in good accord with other computational works -similar picture was already obtained in study of De Kerpel et al.[100]. Partial atomic charges δ δ δ δ obtained by the NPA are compiled in Table 8 for the copper cation and all the donating atoms. For both Cu(I) and Cu(II) complexes, atomic charges closely reflect the donation from the individual ligands. The partial charges give an insight on how much of electron density from the given atom is approximately involved in the (dative) bonding (cf. Tables 8). It can be noticed that Type B centers exhibit in oxidized state lower donation from cysteine -in accord with spin densities (and AIM analysis) but (generally) slightly larger donation occurs in case of histidine ligands, which is also supported by larger electron densities in Cu-N critical points in AIM analysis. In the case of Cysteine (δ δ δ δ(S) = -0.74 e for isolated [S(CH 3 )] -in vacuo), the partial charge is reduced (in absolute value) to -0.57 e in (I) centers. It shows that polarization by the Cu cation causes shifting of the electron density towards the cation in all the ligands. Part of this electron density is then used for metal donation. In case of the Cu(II) complexes, the effect is even more pronounced resulting in the reduction of partial charge of sulfur atom as low as -0.13 e. This situation can be recognized as a transfer of almost the whole electron from Cys to Cu(II). The same description follows from the discussion of the electron spin densities (see above). In the Met residue, positive δ δ δ δ(S) values were found on the S atom in all cases. In the reduced centers, the sulfur partial charges do not change substantially compared to an isolated molecule, since the Cu(I)-S(Met) distances are very long. In the Cu(II) complexes, the partial charges in vacuo are lowered by about 0.04 e compared to the isolated S(CH 3 ) 2 molecule reflecting some shortening of the Cu-S(Met) bond in oxidized centers as well as their positive ∆ ∆ ∆ ∆E BE bond energies. Results for the thioether ligand can be also compared with our previous study[79], where the [Cu(H 2 S) n ] 2+/+ complexes (n = 1 -4) with neutral hydrogen sulfide displayed δ δ δ δ(S) = -0.1 and -0.2 e for the Cu(II) and Cu(I) systems, respectively. These more negative values (in contrary to Met) are product of a much shorter coordination distances in these inorganic complexes. Partial charge on the Cu atom was estimated to be about 0.96 and 0.69 e for the oxidized and reduced A centers in vacuo, respectively. For the Cu(II)/Cu(I) cations in the B centers, slightly higher values δ δ δ δ = 1.04/0.74 e were obtained. The difference can be explained by a weaker Cu-Cys donation in the B complexes, which represents the main contribution to the partial charge reduction of the copper cations. When solvation effects are considered, the Cu partial charge increases with a decreasing dative interaction of the ligands as a consequence of a bond elongation. It is a result of the more effective screening of electrostatic interactions, which enhances especially in the case of negatively charged model of the Cysteine side chain.

  of copper. The second type (7b) contains nonbonding admixture of S(Cys)3p z orbital perpendicular to the of Cu. The fully optimized A centers and the protein structures of B centers possess this type of SOMO.

Figure 5 ). 6 .

 56 Figure 5). Higher relaxation energy is released within the full optimization for the

Figure 1 :

 1 Figure 1: PDB structure of plastocyanin (1KDI).

Figure 2 :

 2 Figure 2: Models of a) Type A and b) Type B centers of blue copper proteins. Family of Type A centers includes following peptides: amicyanin, auracyanin, plastocyanin, and rusticyanin. In case of the Type B centers, following proteins were considered: mavicyanin, stellacyanin, and umecyanin.

Figure 3 :

 3 Figure 3: Example of 4-coordinated [Cu(H 2 S) 2 (NH 3 ) 2 ] 2+ complex, where transition from Cu(II) to Cu(I) state goes along with pronounced structural change from squareplanar to tetrahedral ligand arrangement. When the system is even more relaxed, it evolves to 2-coordinated structure.

Figure 4 :

 4 Figure 4: Structures of reduced a) Type A and b) Type B centers after the full optimization in vacuo. Optimized structure of oxidized centers of c) Type A and d) Type B are presented too.

Figure 5 :

 5 Figure 5: The relaxation energies (in kcal/mol), determined as relative differences of ∆ ∆ ∆ ∆E Stab energies of constrained and fully optimized structures for: _■_ oxidized Type A centers, --♦-oxidized Type B centers, -.-▲-.-reduced Type A centers, and …▼… reduced Type B centers.

Figure 6 :

 6 Figure 6: Plots of spin isodensities (ρ s = 0.01 e/Å 3 ) for selected Cu(II) complexes. Oxidized a) Type A (plastocyanin) and b) Type B (umecyanin) centers with imposed protein geometry (constrained opt.) and c) Type A and d) Type B centers after full optimization are considered.

Figure 7 :

 7 Figure 7: Three different schemes of SOMO with examples. The a) scheme occurs in Type A centers (opt. phase I), b) scheme is in Type A (full opt.) and Type B (opt. phase I) complexes, and c) can be found in fully optimized Type B center.

  Having obtained the optimized structures, further analyses were performed at the higher computational level. More accurate 6-311++G(2df,2pd) basis set was used for the H, C, N, and O atoms. At this level, the Stuttgart ECP pseudopotentials were used for both the Cu and S atoms. The basis sets on the copper/sulfur atoms were consistently enlarged by spd/sp diffuse functions and 2fg/2df polarization functions (α f = 1.00, 0.26, α g = 0.66/ α d = 0.92, 0.29, α f = 0.57).

	F r o Energy analysis consists of evaluation of stabilization ∆ ∆ ∆ ∆E Stab and bond ∆ ∆ ∆ ∆E BE
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	7

was used to simulate two environments: water (with relative permittivity ε = 78) and protein-like environment (ε = 4).

Table 1 .

 1 While the thiolate Cu-S(Cys) bond remains practically unchanged, the

	3.07

In the Type B center, two additional H-bonds are formed between the carbonyl group and histidines (d(O..H(His)) = 2.16 and 2.41 Å). Substantially weaker hydrogen bonding [d(S(Met)..H(His)) = 2.96 Å and d(S(Cys)..H(His)) = thioether Cu-S(Met) and carbonyl Cu-O(Gln) bonds elongate to 4.4 and 3.7 Å, respectively. The difference between the Cu-N(His1) and Cu-N(His2) distances

Table 4 .

 4 Partitioning of the system into the ligand and the remaining part of the complex provides a useful insight into the strength of individual dative bonds. The

	Cu-S(Cys) coordination clearly dominates in all investigated complexes with the
	∆ ∆ ∆ ∆E BE energies in vacuo around 240 and 130 kcal/mol for the oxidized and reduced
	centers, respectively. Such high bonding energies are products of a dative

coordination enhanced by a strong electrostatic interaction between the Cu 2+/+ cation and cysteine model [S(CH 3 )] 1-. For the Cu(II) complexes, ∆ ∆ ∆ ∆E BE of the Cu-N(His)

In the reduced centers, the Cu-L coordination bonds in protein imposed structures are noticeably weaker than in the Cu(II) case, mainly due to the reduction of electrostatic contributions. Actually, the ∆ ∆ ∆ ∆E BE energies of the Cu-Met/Gln bonds exhibit negative values in our model. This can be explained by the fact that the ∆ ∆ ∆ ∆E BE energies also contain an interligand repulsion, which can be stronger than weak Cu(I)-L attractive interaction. For the fully optimized complexes in vacuo (without protein constraints) no negative bonding energies were obtained since the repulsive ligandligand interactions were minimized. In the geometry discussion, it was proposed that the reduced centers tend to lower coordination numbers. A partial support for this

Table 5 .

 5 It can be noticed that the IP Electron distributions of the redox centers of blue copper proteins were investigated in terms of the NPA partial charges, AIM and molecular orbital (MO) analysis, and spin

	values decrease rapidly with an increasing relative permittivity. In the constrained
	optimizations, lower values of both vertical and adiabatic IP for the B centers can be
	seen regardless of the type of the environment. Optimized structures in all three
	environments display an inverted order of IP (IP(type A) < IP(type B)) for vertical
	transitions and for the adiabatic IP adiab energies in vacuo and protein-like

neighborhood. It can be noticed that structural relaxation leads to increase of IP values as one could expect according to larger stabilization of reduced complexes from Figure 5. The more pronounced stabilization of reduced Type B centers (especially in gas phase and protein-like environment) is responsible for the inverted order of IP in these optimized (relaxed) complexes.

Experimentally obtained redox potentials E 0 for the studied proteins are presented in Table

6

. There are also recent calculations on the T-1 blue-copperproteins reduction potential, which treats larger model and reproduce experimental results fairly accurately

[START_REF] Li | Determinants of the Relative Reduction Potentials of Type-1 Copper Sites in Proteins[END_REF]

. Although our calculations are not able to reveal E 0 values for individual proteins, it is at least partially possible to compare the relative differences between the A and B families. The difference of 0.14 eV obtained from experiment is in good agreement with the difference IP adiab (type A) -IP adiab (type B)

Table 1 .

 1 This feature can be explained by a markedly different, (near-tetrahedral, i.e. more symmetrical) geometry of this complex. Coordination distances for oxidized Cu(II) and reduced Cu(I) complexes of the A and B centers of blue copper proteins (in Å). Bond lengths for structures obtained in both constrained and full optimization are presented. The optimizations were performed in vacuo, protein-like, and water environment. Bold indicates the averaged PDB distances, which are not changed in the optimization phase I.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 2 .

 2 Average

			Cu-X (S, O, N) bond lengths (in Å) for small
	[Cu(H 2 S) m (H 2 O) n (NH 3 ) k ] 2+/+ models[79] and corresponding values obtained from
	CSD database by Katz et al.[23].
	Cu(II)		Small models			CSD
		4-coord. 5-coord. 6-coord. 4-coord. 5-coord.	6-coord.
	Cu-N	2.03	2.09	2.28	1.98	2.03	2.34
	Cu-O	2.01	2.06	2.13	1.93	2.07	2.36
	Cu-S	2.40	2.49	2.41	2.28	2.43	2.72
	Cu(I)		Small models			CSD
		2-coord. 3-coord. 4-coord. 2-coord. 3-coord.	4-coord.
	Cu-N	1.91	2.02	2.12	1.90	1.98	2.04
	Cu-O	1.88	2.04	2.21	1.84	2.14	2.05
	Cu-S	2.19	2.29	2.35	2.17	2.26	2.33

Table 3 .

 3 ∆ ∆ ∆ ∆E stab stabilization and ∆ ∆ ∆ ∆E TCE sterically corrected stabilization energies (in kcal/mol) for oxidized Cu(II) and reduced Cu(I) complexes of the A and B centers of blue copper proteins. Energies are presented for structures obtained by all optimization schemes with and without environment effects.

			Oxidized Type A center	Oxidized Type B center
			in vacuo protein-like water in vacuo protein-like water
	∆ ∆ ∆ ∆E stab	Constr. opt Full opt.	589.9 596.4	330.3 336.3	256.7 263.3	595.4 602.0	332.3 339.0	258.2 265.1
	∆ ∆ ∆ ∆E TCE	Constr. opt Full opt.	616.9 619.7	372.6 380.4	298.7 309.8	631.2 623.9	383.9 384.5	308.1 308.4
			Reduced Type A center	Reduced Type B center
			in vacuo protein-like water in vacuo protein-like water
	∆ ∆ ∆ ∆E stab	Constr. opt Full opt.	229.4 242.0	125.1 132.9	96.9 104.6	229.4 250.9	123.9 138.4	95.0 106.9
	∆ ∆ ∆ ∆E TCE	Constr. opt Full opt.	254.6 254.9	165.9 166.2	137.5 134.9	258.3 254.3	169.3 164.9	140.6 131.2
	Table 4. ∆ ∆ ∆ ∆E BE bond energies (in kcal/mol) for oxidized Cu(II) and reduced Cu(I)
	complexes of the A and B centers in various environments.
			Oxidized type A center	Oxidized type B center
			in vacuo protein like water in vacuo protein like water
		Constr. opt 242.6

Table 5 .

 5 Vertical and adiabatic ionization potentials (in eV) representing transition from reduced to oxidized state. Values are presented for both the A and B centers obtained in vacuo as well as in protein-like and water environment.

			Vertical ionization		Adiabatic ionization
		Center in vacuo protein-like water in vacuo protein-like water
	Constr. opt.	A B	5.03 4.85	4.29 4.13	4.15 4.00	4.97 4.73	4.23 4.04	4.08 3.93
	Full opt.	A B	5.92 6.03	5.08 5.27	4.92 5.06	5.22 5.36	4.30 4.39	4.33 4.13

Table 6 .

 6 Redox potentials (in mV) for selected blue copper proteins[START_REF] Solomon | Electronic Structures of Active Sites in Copper Proteins: Contribution to Reactivity[END_REF].

	Page 35 of 45		Molecular Physics
		protein	redox potential
		amicyanin	261
		auracyanin	240
	Type A	plastocyanin	370
		rusticyanin	680
		averaged	388
		mavicyanin	285
	Type B	stellacyanin umecyanin	184 283
		averaged	251

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 7 .

 7 Spin electron density (in e) localized on copper, cysteine model, and the rest of complex obtained by the NPA method for oxidized Cu(II) complexes of the A and B centers.

					Spin density
				Cu	S(Cys) N,N,S/O
		in vacuo	Constr. opt Full opt.	0.34 0.38	0.56 0.50	0.09 0.12
	Type A center	Constr. opt protein-like Full opt. Constr. opt water Full opt.	0.40 0.44 0.43 0.48	0.48 0.39 0.44 0.32	0.12 0.18 0.13 0.20
		in vacuo	Constr. opt Full opt.	0.39 0.53	0.51 0.30	0.10 0.17
	Type B center	Constr. opt protein-like Full opt.	0.47 0.56	0.41 0.25	0.13 0.19
		water	Constr. opt Full opt.	0.50 0.57	0.36 0.23	0.14 0.19

Table 8 .

 8 Partial

			Oxidized Type A center	Oxidized Type B center
			in vacuo protein-like water in vacuo protein-like water
	Cu	Constr. opt 0.961 Full opt. 0.976	1.016 1.032	1.043 1.059	1.042 1.166	1.110 1.197	1.140 1.208
	S(Cys)	Constr. opt -0.130 Full opt. -0.209	-0.225 -0.312	-0.272 -0.174 -0.387 -0.350	-0.292 -0.424	-0.346 -0.457
	N(His)	Constr. opt -0.634 Full opt. -0.639	-0.650 -0.661	-0.655 -0.625 -0.646 -0.645	-0.651 -0.639	-0.656 -0.658
	N(His)	Constr. opt -0.652 Full opt. -0.632	-0.657 -0.642	-0.661 -0.644 -0.665 -0.628	-0.643 -0.655	-0.649 -0.643
	S(Met)/	Constr. opt 0.190	0.184	0.183 -0.738	-0.755	-0.770
	O(Gln)	Full opt.	0.225	0.237	0.258 -0.762	-0.770	-0.775
			Reduced Type A center	Reduced Type B center
			in vacuo protein-like water in vacuo protein-like water
	Cu	Constr. opt 0.697 Full opt. 0.648	0.712 0.687	0.715 0.686	0.737 0.675	0.750 0.693	0.752 0.698
	S(Cys)	Constr. opt -0.573 Full opt. -0.561	-0.647 -0.633	-0.678 -0.570 -0.659 -0.546	-0.652 -0.624	-0.684 -0.654
	N(His)	Constr. opt -0.594	-0.613	-0.626 -0.590	-0.623	-0.635
		Full opt.	-0.598	-0.622	-0.628 -0.586	-0.623	-0.635
	N(His)	Constr. opt -0.604	-0.630	-0.643 -0.603	-0.616	-0.631
		Full opt.	-0.582	-0.622	-0.637 -0.603	-0.614	-0.651
	S(Met)/	Constr. opt 0.223	0.205	0.197 -0.676	-0.708	-0.735
	O(Gln)	Full opt.	0.194	0.185	0.184 -0.718	-0.732	-0.742

atomic charges (in e) on copper and coordinated atoms S(Cys), N(His), S(Met), and O(Gln) obtained by the NPA method.

Table 9 .

 9 EPR g-factors for all studied Cu(II) centers (in a.u.). Experimental values taken from ref.[START_REF] Solomon | Electronic Structures of Active Sites in Copper Proteins: Contribution to Reactivity[END_REF].

	Page 37 of 45			Molecular Physics
		Exp.	in vacuo	protein-like	water
			Constr. Full opt. Constr. Full opt. Constr. Full opt.
	Type A	g ║ 2.226 2.115 g ┴ 2.045 2.052	2.131 2.056	2.131 2.058	2.137 2.061	2.138 2.061	2.143 2.065
	Type B	g ║ 2.297 2.123 g ┴ 2.051 2.058	2.085 2.079	2.140 2.066	2.090 2.082	2.147 2.070	2.120 2.073

URL: http://mc.manuscriptcentral.com/tandf/tmph 38 Captions to Figures:

Acknowledgement

This study was supported by MSM 0021620835 grants. We are grateful for the access to the excellent computational resources of the Charles University Meta-Centrum in Prague. A special thank must be given to the KFCHO department computer cluster administrator Dr. Šimánek.


Tomáš Boublık
email: boublikt@ujep.cz

e w O n l y Distribution functions of a heteronuclear hard dumbbell at a hard wall. Pure fluid and mixtures with hard spheres

Keywords: Inhomogeneous system, Distribution function, Heteronuclear Dumbbell, Residual chemical potential, Geometric quantities of enlarged dumbbell Inhomogeneous system, Distribution function, Heteronuclear Dumbbell, Cavity function, Residual chemical potential, Geometric quantities 2

Distribution functions (angle-averaged site/wall correlation functions) in the inhomogeneous systems of a heteronuclear hard dumbbell (of the site diameter ratio σ 2 /σ 1 = 0.5 and L * = 0.625) were studied theoretically by applying a simple method to determine the background (cavity) correlation function from the residual chemical potentials of the hard dumbbell, hard wall and the corresponding combined body. Pure fluid and approximately equimolar mixtures of dumbbells with hard spheres were considered at several densities. Theoretical distribution functions for the most important site/wall distances were calculated and compared with simulation data; fair agreement was found. Advantage of the proposed method is its simplicity (the residual chemical potential depends, besides the packing fraction, only on the geometric characteristics of the hard body -its volume, surface area and mean radius); the method provides quick information on the density profiles in different inhomogeneous systems.

1. Introduction Knowledge of the density profiles in the inhomogeneous systems of hard bodies represents an important input for theories of adsorption, catalysis or behavior of nano-systems. The problem of determining the density profiles or distribution functions has been tackled by simulation technique [1]- [6] or theoretically by employing the solution of the modified Orstein-Zernike equation [7]- [8] or by applying the density functional theory [9]- [11]. Solution of both types of theoretical expressions have been obtained by different numerical methods; some input data (for the corresponding homogeneous systems) are used.

An alternative to these rather involved approaches represents a method, based on the determination of the background (cavity) correlation function (bcf) from the residual chemical potential of a pair of interacting hard bodies and of the corresponding "combined body"(CB), resulting from the overlap of the particles of the considered pair. The formulation of the method was given by Meeron and Siegert [12] and further developed by Ballance and Speedy [13], Speedy [14], Labík et al. [15]- [16] and the present author [17]- [18]. In the original approaches, dealing with the hard sphere (hs) system, the characteristic geometric quantities -volume, V , surface area, S, mean radius, R (mean curvature integral, divided by 4π), for the dumbbell with the site-site distances larger than a sum of radii of the both interacting hs's were obtained as a prolongation of formulas for V and S holding for the shorter site-site distances; in the case of R, the value calculated for the corresponding convex body, enveloping the given hard dumbbell (hd) was used. Such approximations yielded very precise values of the distribution functions for distances slightly larger than the contact one, but diverge at larger hs's distances.

This fact and difficulties with the definition of the combined body in the case of inhomogeneous systems led us to the introduction of a new CB -the enlarged hard dumbbell (ehd), i.e. a hard body originated when the probe hard sphere is rolled in contact with the sites of the original hd [19]. This new model of CB makes it possible to determine bcf in the most important range of hb distances and determine bcf of inhomogeneous systems.

The approach was recently applied to inhomogeneous systems of hs near a planar hard wall (hw), hs in a hard pore, homonuclear hb/wall systems with hb in the parallel and perpendicular orientations (with respect to the wall plane), full site-wall distribution function and to ternary hs systems with hs's differing considerably in their size [19]- [22]. It was found that within the considered approach, rdf's for different pairs of components might be easily calculated just by modifying the probe hard sphere radius and the non-sphericity parameters of the hs system under consideration (characterized by the composition of different hs's,).

In this contribution we study an application of the simple method to determine the site/wall distribution functions in the pure heteronuclear dumbbell system and in the hhd-hs mixtures. In the first part of the paper a brief outline of the method to determine the background correlation function is given and the way of determining the geometric characteristics of different ehd/wall pairs is summarized. In the next section the details of the determination of geometric quantities of the heteronuclear dumbbell and of the corresponding combined hard bodies are discussed. Next part brings a comparison of the calculated distribution functions with the simulation results (taken from figures of the original paper [23]). Conclusions summarize the main traits of the approach and possible ways of further improvement.

Theory

The present theory is grounded on the relation between the background (cavity) correlation function [12] (bcf) Y ,

Y (x) = exp[u(x)/k B T ]g(x)
(where g is the pair distribution function, u -pair potential, x -distance) and the residual chemical potentials of a pair of interacting hard molecules and of the "combined body" (CB) (originated by overlapping of the given molecular pair), infinitely diluted in a system containing both the molecules i, j. For Y it holds true

k B T ln Y (x) = ∆µ i + ∆µ j -∆µ CB (x)
(where T -temperature and k B -Boltzmann constant). In the case of hard bodies and distances x larger than the contact one, the pair potential u(x) = 0 and bcf equals to the correlation function, g(x).

The residual chemical potential of hard spheres [or non-spherical convex hard bodies (hb)] can be obtained as a derivative of the residual Helmholtz function, ∆A, which is a function of the packing fraction, y, and nonsphericity parameters. For the fused hard sphere model or ehd's, no general/exact expression stands at our disposal; thus, expressions derived for ∆A of hard convex body systems have to be used.

The choice of the expression for the residual Helmholtz function depends on the system under consideration. In the case of the inhomogeneous systems, which are related to the colloidal limit of the hard body mixture, the adequate choice of the expression for ∆A is that, which meets the self-consistent constrains (SC) (see Ref [24]- [25]). The following equation meets fairly well the SC conditions and represents thus the proper choice of the relation for ∆A (for a pure hb fluid)

∆A N k B T = -ln(1 -y) + 3αy (1 -y) + βy 2 (3/2 -5y/12) (1 -y) 2 (1)
Parameters of non-sphericity, α and β for non-spherical hard body systems are defined as

α i = R i S i /3V i β i = Q i S 2 i /9V
S i /4π ≤ Q ≤ R 2 i holds true.
For mixtures of components k = i, j..., the non-sphericity parameters are

α s = rs/3ρv β s = qs 2 /9ρv 2 (2)
where

r = ρΣx k R k , s = ρΣx k S k , q = ρΣx k Q k and v = y = ρΣx k V k
and the extended expression, for the residual chemical potential, corresponding to (1) reads as

β∆µ = -ln(1 -y) + y (1 -y) [3α(R * + S * ) + V * ]+ y 2 (1 -y) 2 [β(Q * + 2S * )(3/2 -y/3) + 3αV * ] + y 3 (1 -y) 3 (7/3 -y/6)βV * (3) Here R * , S * , V * , Q * are dimensionless quantities, R * = R/R ref , S * = S/S ref , V * = V /V ref , Q * = Q/Q ref .
The reference quantities in the case of pure fluids are those of the considered system, thus

R * i , S * i , V * i , Q * i , equal to one, whereas R CB * = R CB /R i , S CB * = S CB /S i , etc. For mixtures of components k = i, j.. R ref = Σx k R k , generally X ref = x i X i
, where x i denotes here the mole fraction of component i. Defining

∆r = R * i + R * j -R CB * , ∆s = S * i + S * j -S CB * , ∆v = V * i + V * j -V CB * ∆q = Q * i + Q * j -Q CB * one can write ln Y = -ln(1 -y) + y (1 -y) [3α(∆r + ∆s) + ∆v]+ y 2 (1 -y) 2 [β(∆q + 2∆s)(3/2 -y/3) + 3α∆v] + y 3 (1 -y) 3 (7/3 -y/6)β∆v (4)

Pure hard dumbbells

As mentioned before, eq (4) was derived for hard convex bodies. In order to describe the behavior of non-convex bodies -diatomics, (triatomics etc.) and combined bodies -we adopt a model of the enlarged hard dumbbell (ehd) (generally enlarged fused hard body), originates when the probe hard sphere is rolled over hard dumbbell (in general case over a fused hard sphere body). The enlarged hard body possesses a smooth surface area (without any discontinuities and cusps). The idea of ehd is especially helpful in the case of the description of inhomogeneous systems of hard bodies near a hard wall.

In order to determine the geometric quantities of ehd, we employed the method due to Connolly [26,27], in which the ehd surface area and volume were obtained as sums of the convex and saddle contributions (no concave part was present). For the heteronuclear hard dumbbell, studied in this work, with the site-site ratio p = σ 2 /σ 1 = 0.5, the reduced site-site distance L * = l/σ 1 = 0.625 (and radii r i = σ i /2), the ehd surface area is given by two convex and two saddle contributions, see Figure 1. The surface area of the studied ehd is

S 1 = 2πr 2 1 (1 + sin θ 1 ) + 2πr 2 2 (1 + sin θ 2 ) + 2πr p [X(θ 1 + θ 2 ) -r p (sin θ 1 + sin θ 2 )]
(5) Here r p is the probe radius, X = (r 1 + r p ) cos θ 1 and

θ i = arcsin[d i /(r i + r p )]. Then S 1 /πσ 2 1 = [(1 + sin θ 1 ) + p 2 (1 + sin θ 2 )+ r * p (1 + r * p )(θ 1 + θ 2 ) cos θ 1 -r * 2 p (sin θ 1 + sin θ 2 )]/2 (6)
(where r * p = r p /r 1 ). The ehd volume of the considered dumbbell is similarly

V 1 /(πσ 3 1 /6) = 1 2 [(1 + sin θ 1 + cos 2 θ 1 sin θ 1 /2) + p 3 (1 + sin θ 2 + cos 2 θ 2 sin θ 2 /2)] + 3 4 r * p (1 + r * p ) cos θ 1 [(1 + r * p ) cos θ 1 (sin θ 1 + sin θ 2 ) -r * p (cos θ 1 + cos θ 2 )] + r * 3 p 4 (cos 2 θ 1 sin θ 1 +cos 2 θ 2 sin θ 2 +2 sin θ 1 +2 sin θ 2 )- 3 4 r * 2 p (1+r * p ) cos θ 1 (θ 1 +θ 2 ) (7)
The mean radius is given by sums of the convex and saddle contributions, obtained from eq (5) by differentiation with respect to r 1 resp. r 2 and minus derivative of the saddle parts of surface areas with respect to r p ; thus

R 1 /r 1 = 1 2 [(1+sin θ 1 )+p(1+sin θ 2 )+r * p (sin θ 1 +sin θ 2 )-(1+r * p ) cos θ 1 (θ 1 +θ 2 )] (8)
7 The structure of a pure hard diatomic near a hard wall is usually determined as the limiting case of a pair of hard dumbbell/hard sphere (with an infinitely large radius and concentration ≈ 0). To determine the distribution function within our approach we consider instead of (quantities of) the infinitely large hard sphere a hard disk of the radius r d (see Figure 2a for the perpendicular orientation of ehd with respect to the wall), with the zero volume and zero mean radius and with the surface area equal to S d = πr 2 d . The geometric quantities of the combined body possess then following forms

S CB /πσ 2 1 = S 1 /πσ 2 1 - 1 2 [(1-sin θ)-r * p (1+r * p )(π/2+θ) cos θ+r * 2 p (1+sin θ)] (9) 2R CB /σ 1 = R 1 /r 1 - 1 2 [(1-sin θ)+(1+r * p )(π/2+θ) cos θ/2-r * p (1+sin θ)] (10) and V CB /(πσ 3 1 /6) = V 1 /(πσ 3 1 /6) - 1 2 (1 -sin θ -cos 2 θ sin θ/2) + 3 4 r * p (1 + r * p ) cos θ[(1 + r * p ) cos θ(1 + sin θ) -r * p cos θ sin θ -r * 2 p (1 + r * p )(π/2 + θ) cos θ] + r * 3 p 2 (cos 2 θ sin θ + 2 sin θ + 2) ( 11
)
Here θ is an angle defined as θ = arcsin(2z/(1 + r * p ) -1) and z is the reduced distance (related to σ 1 ) of the site from the disk.

For distances z > [r * p + (1 + 2r * p )]/2 the combined body decomposes into two parts ending in cusps, see Figure 2b, and the integrals which determine the contribution of the saddle parts to the CB surface area or volume have to be considered for bounds φ -θ, [where r * p cos φ = (1 + r * p ) cos θ]. Thus, one has to subtract contributions for angles (-φ, φ) -i.e. δR, δS, δV , from the above expressions. These contributions are

δR/r 1 = -φ(1 + r * p ) cos θ/2 + r * p sin φ ( 12
) δS/πσ 2 1 = φr * p (1 + r * p ) cos θ -r * 2 p sin θ ( 13
)
δV /(πσ

3 1 /6) = 1 2 [3r * p (1 + r * p ) 2 cos 2 θ sin φ -3r * 2 p (1 + r * p ) cos θ(cos φ sin φ + φ) + r * 3 p sin φ(cos 2 φ + 2)] (14)
The above expressions are valid when the structure of the larger site/wall is considered. In the case of smaller site/wall functions, the probe sphere for small values of distance z contact -in special orientations, see Figure 3 -the larger instead of the smaller site. Then in eqs ( 9)- (11) it is necessary to consider the reduced quantities of the larger site, equal in all the cases to one [instead of 2R 1 /σ 1 , S 1 /(πσ 2 1 ) and V 1 /(πσ 3 1 /6)]. Functional Q (CB) is approximated as

Q (CB) = S (CB) (4πR (CB)2 /S (CB) ) 1/4 .

Mixtures

In this paper we study a special case of hard body mixtures, namely mixtures of the heteronuclear hard dumbbells and hard spheres (whose diameter is equal to the diameter of the larger hhb site, σ hs = σ 1 ).

The hard body/wall correlation function of a non-spherical particle can be evaluated from eqs ( 9)-( 11)which are applicable with the following changes: (i) the reference quantities R ref , S ref , V ref , Q ref for a binary mixtures are defined as X ref = Σx i X i (where x i denotes the mole fraction and X a single geometric quantity). (ii) The non-sphericity parameters α, β in the above equations (for pure fluid systems) should be replaced by parameters α s , β s defined in eq (2). (iii) the reduced probe radius r * p is determined as a weighted average from those of single components (see the discussion in Ref [21]).

It is obvious that in the case of hard sphere/hard wall correlation function, one has to consider instead of R 1 /(σ 1 /2), S 1 /πσ 2 1 , V 1 /(πσ 3 /6) and

Q 1 /(σ 1 /2) 2 in our case the hs values R 2 /(σ 1 /2) = S 2 /πσ 2 1 = V 2 /(πσ 3 1 /6) = Q 2 /(σ 1 /2) 2 = 1.

Results and discussion

Firstly we consider a pure heteronuclear hard dumbbell (hhd) system. The studied dumbbell is characterized by the reduced site-site distance L * = l/σ 1 = 0.625 and the ratio of the site diameters p = σ 2 /σ 1 = 0.5. Because of the low p-value, the perpendicular orientation of hhd to the hard wall is decisive for determination of the correlation function of the larger site. For this orientation the background correlation function is evaluated first, and consequently multiplied by the probability P of the general orientation of the dumbbell axis with respect to the normal of the wall. It is evident that for the given distance of the larger site from the wall, z, orientations with angles ζ > π/2 + arcsin((z -r 2 )/L * ) are forbidden. Thus, in the simplest approach we take as the probability function the ratio of acceptable to all The calculated correlation functions for packing fractions y = 0.407 and 0.356 are compared with simulation data of Nezbeda et al. [23] in Figure 4. It is obvious that the geometric method yields a fair prediction of the correlation function for both the studied systems with exception of the distances close to the contact value (where, however, experimental values taken from figures in Ref [23] are rather uncertain).

Evaluation of the background correlation function of the smaller site of the considered dumbbell (with diameter σ 2 = 0.5σ 1 )/wall in the perpendicular orientation is more complex, because for the lowest z values, the probe sphere with r p /r 1 = 0.9 (i.e. only slightly smaller than the larger site of hhd) contacts the larger site instead of the smaller site situated closer to the wall; for parallel orientation and majority of general orientations the combined body is determined mainly by the site -probe contacts. To this end, we evaluate bcf's for the larger site/wall situation, determined, however, for wall-site distance z + d av with the average distance d av resulting from d max = L * and that of a minimum, where the larger site is in contact with the wall. To get distribution function g W HD2 , we multiply the obtained bcf-value by the probability ratio of allowed and total orientations. Dependence of the distribution function of the smaller hhd site/wall on the distance z, of the center of site 2 from the wall is compared with the experimental data on the most important interval of z in Figure 5. Larger differences, observed in the central part of the graph are due to neglecting the smaller site-probe contributions to the combined body.

Next we consider approximately equimolar mixtures of hhd(1) and hard sphere(2) (hs) near a hard wall. To determine the hs/wall correlation function we calculate bcf; all the reduced geometric quantities of HS are equal to one, however, the reference quantities are given by prescription X = x i X i and instead of α α s has to be used. The reduced probe radius r p /r 1 is here r * p = 0.95. For the hs/wall function, there are no forbidden orientations, therefore g W HS = Y . Dependence of g W HS on the reduced distance z for mixtures with the packing fraction y = 0.407 and composition x 1 = 0.514 and y = 0.306, x 1 = 0.523 are compared with values taken from a figure in Ref [23]. Figure 6 shows a fair agreement of theory with simulation results. Evaluation of the correlation function g W HD1 resembles much the approach described in the previous part for pure hhd, however, with values of r the hhd+hs solution. The correlation function g W HD1 was again obtained by multiplying bcf for the given distance z by the probability function (same as in the case of pure hhd). Theoretical correlation functions calculated for the packing fractions and compositions mentioned above, agree well with experimental results; at z = 0.875 cusps are found, corresponding to a transition from the restricted to free rotation, see Figure 7.

Determination of the distribution function of the smaller site /hard wall in the mixture of hhd + hs follows the way used in the previous case, however, the combined body of the hhd+hard wall is taken as a result of contact between larger site(1) and the probe hard sphere, as in the case of pure hhd fluid. Comparison of the distribution function of smaller site/hard wall in the hhd+hs mixture with the packing fraction y = 0.358, and the mole fraction x 1 = 0.525 with simulation results are shown in Figure 8 for the interval of reduced distances z ∈ (0.25, 1). It is apparent that the effect of hard sphere (in mixture) on the course of the theoretical g W HD2 is small and the distribution function resembles that for pure hhd fluid; the same is true for experimental dependence of g vs z for a mixture and pure fluid.

Conclusions

The present study of heteronuclear hard dumbbell with the considerably different site diameters and relatively short reduced site-site distance brings in comparison with previously considered homonuclear dumbbell (Ref [20]) some simplifications and some new problems. Simplifications consist in fact that the distribution function determined on the basis of background correlation function for the perpendicular orientation of hhd with respect to the wall, multiplied by the probability of different general orientations is more accurate than in the case of homonuclear hard dumbbells. This is a result of the smaller effect of the other site on the range of acceptable orientations. On the other hand, determination of the the smaller site/hard wall bcf faces difficulties, as in the case of the perpendicular orientation and shortest distances z, the probe sphere (with r p ≈ r 1 ) contacts firstly the more removed larger site and only for larger distances z it contacts the smaller site, situated near the wall. For general orientations, interactions of the probe with the larger sphere of hhd are more important. Therefore we employ for the determination of the geometric quantities of the combined body the bcf for the removed larger site+hard probe situation. This fact causes some uncertainties within the procedure of determining the distribution function. In this and previous contributions we tried to show the potential of the simple method to determine distribution functions in inhomogeneous systems (on the basis of residual chemical potential of the hard fluid molecules, of the hard wall and of the combined body). For hard dumbbells and other non-spherical models, an improved but more involved version can be worked out, where the geometric quantities of the combined body for the single orientations are determined for each distance z and some averaging procedure to get the average site-wall distribution functions is applied.

4 Page

 4

6 Page

 6

9 Page

 9 orientation. For z ≥ z c = p/2 + L * the probability is equal to one; at values z = z c the correlation function discloses a cusp.

10 Page

 10 * p α s , β s and R ref , S ref , V ref , Q ref valid for the mixture and used here for in the case of

11 Page

 11

Figure 1 .

 1 Figure1. Geometry of the heteronuclear hard dumbbell of the reduced sitesite length L and hard sphere radii r 1 , r 2 and the probe radius r p (L * = 0.625, σ 2 /σ 1 = 0.5).

Figure 2 .

 2 Figure 2. Geometry of the hard dumbbell/hard wall in the perpendicular orientation. z-dumbbell-site reduced distance, θ i -characteristic angle, r pprobe sphere radius, a) for shorter distances, b) for z > [r p + (1 + 2r p ) 1/2 ]/2.

Figure 3 .

 3 Figure3. The same as in Figure2for the smaller site.

Figure 4 .

 4 Figure 4. Larger site/wall distribution function for packing fractions y = 0.407 (o, full line) and y = 0.356 (∆, dash line).

Figure 5 .

 5 Figure 5. Smaller site/wall distribution function for the system with the packing fraction y = 0.356.

Figure 6 .

 6 Figure 6. Approximately equimolar mixture of hard dumbbell+hard sphere. hard spher/wall distribution function, y = 0.407, x 1 = 0.514, (o, full line), y = 0.306, x 1 = 0.523 (∆, dash line).

Figure 7 .

 7 Figure 7. The same as in Figure6but for the site/wall distribution function.

Figure 8 . 14 Page

 814 Figure 8. Smaller site/wall distribution function for the hhd-hs mixture (y = 0.358, x 1 = 0.525) .

2 i 5 Page 6 of 30 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	where V i , S i , R i and Q i are the hcb geometric characteristics -volume, surface
	area, mean radius (exactly the mean curvature integral, divided by 4π); Q
	is a quantity of dimension [l 2 ], for which inequality
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgement

The work has been financially supported by grant No. IAA 400 720 710 of the Grant Agency of the Academy of Sciences of the Czech Republic.