Quantitative prediction of electronic absorption spectra of copper (II)–bioligand systems: Validation and applications

The visible region of the electronic absorption spectra of Cu(II) complexes was studied by time-dependent density functional theory (TD-DFT). The performance of twelve functionals in the prediction of absorption maxima (λ max ) was tested on eleven compounds with different geometry, donors and charge. The ranking of the functionals for λ max was determined in terms of mean absolute percent deviation (MAPD) and standard deviation (SD) and it is as follows: BHandHLYP > M06 ≫ CAM-B3LYP ≫ MPW1PW91~B1LYP~BLYP > HSE06~B3LYP > B3P86~ω-B97x- D ≫ TPSSh ≫ M06-2X (MAPD) and BHandHLYP > M06~HSE06 > ω-B97x-D~CAM-B3LYP~MPW1PW91 > B1LYP~B3LYP > B3P86 > BLYP ≫ TPSSh ≫ M06-2X (SD). With BHandHLYP functional the MAPD is 3.1% and SD is 2.3%, while with M06 the MAPD is 3.7% and SD is 3.7%. The protocol validated in the first step of the study was applied to: i) calculate the number of transitions in the spectra and relate them to the geometry of Cu(II) species; ii) determine the coordination ofaxial water(s); iii) predict theelectronic spectra of thesystems where Cu(II) is bound to human serum albumin (HSA) and to the regions 94–97 and 108–112 of prion protein (PrP). The results indicate that the proposed computational protocol allows a successful prediction of the electronic spectra of Cu(II) species and to relate an experimental spectrum to a specific structure.


Introduction
The molecular modeling of the bioinorganic structures, the electronic and spectroscopic properties of metalloenzymes and their model systems is a subject of growing interest over the last decade, due mainly to the effective development of computational methods. Among them, density functional theory (DFT) [1] and its time-dependent extension (TD-DFT) [2,3] have been frequently used because of their relative low computational cost and implementation in most of the available computational software. Moreover, the effect of the solvation of the chemical species can be treated routinely. Several studies and comprehensive reviews are available in this field focusing on the characterization of structural, functional and spectroscopic features of (transition metal)-bioligand systems [4][5][6][7][8]. With respect to the spectroscopic properties of these systems, electron paramagnetic resonance (EPR) [9][10][11][12][13] and 1 H and 13 C nuclear magnetic resonance (NMR) parameters [14][15][16][17] have been calculated with great accuracy by computational techniques and the effect of the pseudo-contact chemical shift on the NMR spectra originating from paramagnetic complexes was successfully simulated with the use of DFT methods [18][19][20][21]. The prediction of infrared (IR) and optical (UV-Vis, circular dichroism (CD) and magnetic circular dichroism (MCD)) spectra have been examined by several research groups. However, the interpretation of these spectra is not a trivial task since the spectroscopic behavior of a metal compound is influenced by various factors -such as the specific donors bound to the metal ion, the geometry and symmetry of the species formed -and the calculated properties are greatly dependent on the applied theoretical methods (i.e. functional and basis set for DFT) [22][23][24][25][26][27][28][29][30]. Therefore, finding robust and suitable methods to calculate accurately the spectroscopic parameters remains a great challenge for the chemists' community.
Copper, which is one of the most abundant transition metal ions in living organisms, has a crucial role in several biological processes. Copper containing metalloenzymes are involved in oxygenation reactions and oxygen transport [31][32][33][34], and in electron transfer reactions [35,36]. Moreover, the role of copper in neurodegenerative diseases is also extensively examined [37,38]. In all cases, copper interacts with the amino acid residues of a protein, although the coordination environment and geometry around the metal center may significantly differ. Consequently, the metal-peptide complexes are excellent choices for mimicking the active site of metalloenzymes and metalloproteins and, for this reason, these systems are thoroughly studied [39,40]. In addition to the biological species, inorganic copper complexes have been used as catalysts and several copper compounds showed anticancer activity [41,42].
Three oxidation states are characteristic for copper: +1 (configuration 3d 10 ), +2 (3d 9 ) and +3 (3d 8 ). Among those, +2 state is the most stable in biological systems -even if specific ligands are able to stabilize the +1 and +3 oxidation states [43,44] -and we will focus here on copper(II) species. In general, Cu(II) prefers tetragonally distorted octahedral geometries, even if it exhibits non-rigid structures such as compressed and elongated octahedral, tetra-and penta-coordinated geometries, which are accompanied by Jahn-Teller effect [45]. Often, this behavior is designed "plasticity".
The complexation of copper(II) can be easily studied by UV-Vis spectroscopy. In water, copper(II) exists as [Cu(H 2 O) 6 ] 2+ ion, which exhibits three electronic transitions that collapse into a single absorption band. The replacement of the coordinated water molecules by stronger donors shifts the absorption maxima (λ max ) to higher energies, and the value and molar absorptivity of λ max significantly depend on the nature of the coordinated atoms. Billo before and Sigel and Martin subsequently developed a general rule for estimating the value of λ max for simple copper(II) species and, particularly, for copper(II)-peptide complexes [46,47]. The empirical equation is based on the interpolation of a large number of Cu(II) electronic spectra (eq. 1) [47]: In Eq. (1), H 2 O, CO, COO -, imide-N, NH 2 and peptide-Ndenote the type of equatorial coordinating donors, while p, q, r, s, t indicate the number of these donors involved in the equatorial binding. In the presence of axial ligands, the estimated values of λ max through Eq. (1) can significantly differ from those experimentally observed [48].
The simulation of the electronic absorption spectra using TD-DFT was frequently discussed and a number of reviews are available in the literature [24,[49][50][51][52]. For the TD-DFT prediction of an electronic absorption spectrum, it is necessary to choose an appropriate combination of functional and basis set. The effect of solvation as specific correction for the excited states can also be considered and several comprehensive studies and reviews were published [24,26,53,54]. At this time, the literature data suggest that, for transition metal complexes, it is not possible to find a general computational condition to get a good agreement between an experimental and a calculated spectrum within the framework of TD-DFT approach. An investigation on four Ir(III) and seven Pt(II) complexes has shown that local functionals (GGA and meta-GGA) underestimate the excitation energies, whereas other functionals such as CAM-B3LYP and M06-2X overestimate these energies [55]; the best performance was obtained using hybrid functionals with the B3 exchange part independently of the correlation part (for example, B3LYP and B3PW91) [55].
The calculation of the electronic absorption spectra of Cu(II) (both for d-d transitions and ligand to metal charge transfers) and Zn(II) complexes (only for ligand to metal charge transfers) through 41 pure GGA, hybrid GGA, pure meta-GGA and hybrid meta-GGA functionals indicated that the hybrid approaches perform better than the others, with B1LYP giving the most accurate results [56]. In a recent publication, we investigated the electronic absorption spectra of square planar nickel(II) complexes using TD-DFT methods [57]. It was shown that the functionals with 25% HF exchange (HSE06, MPW1PW91, PBE0) performed better than the other ones; it was also demonstrated that the effect of the basis set on the calculated spectra is weaker than the effect of the functional and the use of triple-ζ type def2-TZVP set is enough to obtain results in agreement with the experiment [57].
In this work, a systematic study on the TD-DFT prediction of λ max in the visible region for several Cu(II) complexes is discussed. In the first section, a benchmark of eleven species will be considered and the performance of twelve functionals compared. In the second and third sections, the number of experimental and expected transitions and their relation with the metal complex geometry will be analyzed and the effect of the weak coordination of axial water(s) evaluated; in the two final sections, the best functional (BHandHLYP) will be used to calculate λ max for the binding of Cu(II) ion to the N-terminal region of human serum albumin (HSA) and to the regions 94-97 and 108-112 of prion protein (PrP).

Spectroscopic methods
UV-visible spectra of the Cu(II) complexes were recorded from 200 to 1000 nm on a PerkinElmer Lambda 25 or an Agilent Technologies Cary 8454 UV-Vis diode array spectrophotometer.

Density functional theory (DFT) and time-dependent density functional theory (TD-DFT) calculations
The geometries of all the complexes in the benchmark were optimized through Gaussian 09 (rev. D01) software [58] at DFT level of theory, using the hybrid Becke three-parameter B3P86 functional, which is often used because of its high degree of accuracy for structural prediction of transition metal species [59,60], combined with the tripleζ def2-TZVP basis set. The solvent effect was taken into account adopting the Polarizable Continuum Model (PCM) for water [61,62]. The adducts formed upon binding of Cu(II) ion to the region 1-5 of human serum albumin and to the regions 94-97 and 108-112 of prion protein were optimized with the functional B3LYP-D3, and the basis set 6-311g(d,p) for main group elements and SDD[ECP] plus f-polarization functions for Cu.

Statistical analysis
The statistical dispersion of the calculated values of λ max (λ max calcd ) from the experimental ones (λ max exptl ) for the j th complex of the benchmark is defined in terms of percent deviation (PD) and absolute percent deviation (APD). However, in the analysis of a large dataset it is necessary to refer to a mean absolute percent deviation (MAPD) [78], that was proposed previously as a criterion of quality [79,80]: where N is the number of compounds in the dataset. The mean percent deviation is denoted with MPD. Beside MAPD, the standard deviation (SD), related to the data dispersion, is also reported.

Validation: optimization of the structures, and analysis of the spectra and benchmark data
Eleven Cu(II) complexes with various ligands L were studied: [Cu  Table S1 of the Supplementary Material. As an example, Fig. 1 shows the overlap between the optimized geometries (in blue) and the XRD structures (in orange) of [Cu(bipy) 2 ] 2+ and [Cu (H -4 βAla-Gly-βAla-Gly)] 2-: the agreement is very satisfactory and this is obviously a good starting point to calculate accurately the electronic absorption spectra. The Cartesian coordinates of the optimized minima are available in Tables S2-S12 of the Supplementary material.
To ascertain the best computational protocol, the electronic absorption spectra of the eleven Cu(II) complexes were calculated varying the functional and using the basis set def2-TZVP. For each spectrum the visible range (400-800 nm) was taken into account because the d-d transitions between the copper orbitals fall in this region.
The number of d-d transitions is strongly influenced by the geometry and degree of distortion of the copper complexes. Cu(II) species, having a degenerate d 9 configuration in an octahedral crystal field, undergo Jahn-Teller effect and are characterized by elongated or compressed octahedral, square pyramidal or square planar structures, close to D 4h symmetry. In the analysis of the spectra, when only one transition is predicted, it can directly be compared to the experimental λ max value. When the number of excitations is higher than one, the spectrum generated by the superimposition of the corresponding bands is compared to the experimental one. In this work, the spectra were generated with Gabedit software [76] using a Gaussian shape.
The performance of the functionals was assessed using MAPD (Eq. (2)) and SD. The values of MAPD obtained for the eleven Cu(II) complexes are summarized in Fig 2. From an analysis of the data, it emerges that the best performance is reached with the functional BHandHLYP (50% HF exchange) and M06 global hybrid functional (27% HF exchange); their MAPDs are 3.1% and 3.7%, while the SD's are 2.3% and 3.7%, respectively. These functionals are followed by the long-rangecorrected CAM-B3LYP (short range HF exchange 19%, long range HF exchange 65%; MAPD = 9.7% and SD = 5.4%) and by MPW1PW91 (25% HF exchange; MAPD = 12.0% and SD = 5.5%) and B1LYP (25% HF exchange; MAPD = 12.4% and SD = 6.1%). All the other functionals tested give worse predictions and are not recommended to calculate an electronic absorption spectrum of a Cu(II) species.
The hybrid GGA BHandHLYP is a widely used functional, including 50% of HF and mixed LSDA/B88 exchange plus LYP correlation, particularly appropriate to describe hydrogen bonds, potential energy surfaces, transition states, excitation energies, and oscillator strengths with high quality results, comparable as other post-HF methods [81,82]. Its potential in predicting the electronic spectra of transition metal compounds were recently discussed for 'bare' V(IV) species [83,84]. The top performances reached by the meta hybrid-GGA M06 functional are particularly interesting considering its quite low capability of predicting the UV-vis spectra of Ni(II) [57] or EPR parameters of Cu(II) compounds [13], and the general low accuracy observed for functionals with a large number of empirical parameters in the calculation of the electronic density [85].
Considering MPD, the difference between the experimental and predicted λ max values is~0.0% and 0.2% when M06 and BHandHLYP are used, respectively. The use of the other functionals leads to an underestimation of λ max , from −7.5% (BLYP) to −9.3% (TPSSh). The only exception is obtained with M06-2X for which the difference is +15.9%. It is noticeable that the functional B1LYP, considered the best by Holland and Green [56], is outperformed by several functionals, and BHandHLYP, M06 and CAM-B3LYP work significantly better than it. The comparison between the experimental spectrum of [Cu (bipy) 2 ] 2+ and those calculated with various functionals is shown in Fig. 3. Looking at this figure, the functionals can be divided into three groups: one containing BHandHLYP and M06, another containing the functional with intermediate performance such CAM-B3LYP, B1LYP, B3LYP and B3P86, and the last group with the worse functional TPSSh.

Applications: correlation between geometry and type of spectrum
Due to the non-spherical symmetry and Jahn-Teller effect expected for a d 9 ion, the stereochemistry of Cu(II) complexes is characterized by non-rigid structures including several distorted geometries (phenomenon known as 'plasticity') [86,87]. Species with coordination numbers 4, 5 and 6, and geometries square planar, square pyramidal, bipyramidal trigonal, octahedral compressed and elongated, with all the possible distortions, have been reported [88]. The distortion degree determines the number of absorptions in the visible region. For hexacoordinated tetragonal structures, three transitions are expected from symmetry considerations, that collapses into a single experimental band if the distortion is small; the absorption becomes sharper and shifts to higher energy values with increasing the donor strength [89]. One experimental absorption is predicted for trigonal bipyramidal complexes, although the number of electronic transitions expected on the basis of ligand field theory is larger than those detected in the spectrum [88,89]. Theory also predicts one experimental band for tetrahedral species with small distortion [88,89]. Therefore, in contrast with other metal ions -for example oxidovanadium(IV) for which all the three or four d-d excitations are usually detected [90,91] -, for copper(II) the number of electronic transitions expected from symmetry considerations is larger than those experimentally detected in the spectrum. Therefore, it is not easy to determine the exact number of excitations below the broad absorptions and relate the experimental spectrum with the molecular structure of complexes. Computational methods may help to do this and interpret correctly a spectral signal. In the following paragraphs two examples will be discussed. In Fig. 4, the experimental spectra of [Cu(bipy) 2 ] 2+ and [Cu (H -2 GGG)]are given. The spectrum of [Cu(bipy) 2 ] 2+ is characterized by a broad band centered at 730 nm (Fig. 4a). This suggests a tetrahedral arrangement. The computational data indicate that three bands with different intensity are expected in the range 670-720 nm and this accounts for the large linewidth and asymmetry of the band, suggesting a distortion of the structure toward the square planar arrangement which causes the splitting of the e and t 2 orbitals of a tetrahedron. The predicted band at 716 nm is composed by three d-d transitions at -in increasing energy order -716.8, 698.3 and 670.7 nm, from MOs with prevalent Cu-d xy/xz , Cu-d z 2 and Cu-d xy/yz character to that based on Cud x 2 -y 2. This order follows the correlation diagram for the transformation of a regular tetrahedral to a square planar geometry, for which the metal energy order is: Cu-d xy < Cu-d z 2 < Cu-d xz~C u-d yz < Cu-d x 2 -y 2 [92]. This analysis also agrees with the X-ray structure of [Cu(bipy) 2 ] 2+ (ref. [93] and Fig. 1a), for which it has been proposed that the steric hindrance of hydrogens in 3,3′ position precludes the coplanar arrangement of the two bipy molecules [94]. In contrast, for a penta-coordinated structure two bands in the case of a regular trigonal bipyramidal arrangement, and three bands in the case of a distortion toward the square pyramid but at higher wavelengths, are expected [95]. The calculated transitions are listed in Table 1.
Another example is the spectrum of [Cu(H −2 GGG)] − (Fig. 4b). It is more symmetric than the one of [Cu(bipy) 2 ] 2+ , suggesting a less distorted geometry. The absorption maximum at 553 nm indicates a square planar or a very elongated octahedral geometry. The calculations allow us to determine the transitions from the MOs based on Cud xy and Cu-d z 2 orbitals to Cu-d x 2 -y 2, with comparable intensity, accounting for the shape of the experimental band. The predicted behavior is in agreement, in this case too, with the ligand field theory description of the d orbital energy order for a square planar geometry [92]. The absence of a third excitation in the visible region induces to exclude a rhombic distortion, which would remove the degeneration of the Cu-d xz and Cu-d yz orbitals, causing an increase of the transition number. Overall, for [Cu(H −2 GGG)] − , the results indicate that the geometry is almost regular, in line with the predicted structure with copper(II) ion in an almost perfect plane formed by the four donors, i.e. amino-NH 2 , amide-N − , amide-N − , and carboxylate-O − .

Applications: determination of the coordination of axial water
Due to the 'plasticity' of Cu(II) ion and its non-rigid stereochemistry [86,87], the coordination of axial water(s) is possible for a complex with four strong equatorial donors. For a system with two equatorial bidentate L or one tetradentate L' ligand and a monodentate X ligand (X can be a halide ion, a polyatomic anion such as ClO 4 -, NO 3 -, NCS -, etc., or a solvent molecule), the stoichiometry may vary from CuL 2 to CuL 2 X/CuL 2 X 2 and from CuL′ to CuL′X/CuL′X 2 , respectively. For example, with amino-N donors such as ethylenediamine (en) and trien,     2+ . In solution, it is not straightforward to determine if water is bound to copper and, eventually, how many molecules are coordinated to the metal ion since the effect of the axial ligands on the spectroscopic properties (λ max in an electronic spectrum or g and A in an EPR spectrum) is not easily predictable. Obviously, an electronic absorption spectrum changes in the absence or presence of axial water(s), with the change depending on the number of water (one or two) and on the strength of the bond (i.e., the distance Cu-OH 2 ). In these cases, TD-DFT methods can be useful to characterize the species in solution and prove if solvent coordination occurs.
In this section, the effect of the coordination of an axial water ligand to [Cu(trien)] 2+  The large change in the calculated absorption wavelength of the monohydrate complex is due to the short Cu-OH 2 distance expected for the penta-coordinated species (2.277 Å, Fig. 5). An analogous approach was applied for the bis-chelated complex of glycinate: the species [Cu(Gly) 2 ] and [Cu(Gly) 2     (Gly) 2 (H 2 O)] is due to the long Cu-OH 2 distance (2.376 Å experimental [98], and 2.434 Å calculated).

Applications: calculation of the spectrum of the Cu(II)-HSA adduct
Human serum albumin is involved in the transport of various metal ions such as Ni(II) [101], Zn(II) [102,103] and Cu(II) [104]. Two metal binding sites exist in its structure: the ATCUN (Amino Terminal Copper and Nickel) or site I in the N-terminal region that binds Cu(II) and Ni(II) through the donor set NH 2 , N -, N -, His3-N [104], and the multi-metal binding site (MBS) or site II to which Zn(II) and V(IV)O are bound through His67, Asn99, His247 and Asp249 residues [105,106]. The coordination of His3 in the fourth equatorial position plays a fundamental role and promotes the deprotonation of amides of Ala2 and His3 to yield a distorted square planar geometry. Bovine serum albumin (BSA) and rat serum albumin (RSA) behave in the same way. In contrast, albumin of other animals such as dog serum albumin (DSA) and chicken serum albumin (CSA), where His3 is replaced by Tyr3 (DSA) or Glu3 (CSA), does not bind Cu(II) through the ATCUN site [104]. The electronic spectrum of HSA is characterized by a broad band centred at 529 nm (Fig. 7); the λ max exptl is perfectly in agreement with that estimated by Eq. (1) for a (NH 2 , N − , N − , His3-N) coordination.
Exploiting the results of this work, we calculated the electronic absorption spectrum of the model complex containing the first five amino acids of human serum albumin (DAHKS-NHCH 3 ), a model which was already successfully used to simulate the EPR spectrum of Cu (II)-HSA adduct [13]. The Cartesian coordinates of this adduct are reported in Table S13 of the Supplementary Material. The optimized geometry is shown in Fig. 7a and is characterized by a square planar arrangement of the four equatorial donors with a slight distortion of the equatorial plane. The experimental and calculated spectra are compared in Fig. 7b. The agreement is rather good and the calculated value of λ max at 552 nm deviates only 4.3% from the experimental value of 529 nm (see Table 2). The predicted spectrum is broad due to the four transitions expected at 568.7 nm (oscillator strength f × 10 5 = 50), 555.3 nm (f × 10 5 = 10), 545.8 nm (f × 10 5 = 20) and 526.7 nm (f × 10 5 = 20). The large number of transitions compared to the systems with high symmetry is due to the distortion of the equatorial plane of the metal ion. The intense absorptions below 450 nm are assigned as metal to ligand or ligand to metal charge transfers (MLCT and LMCT).
The main d-d character of the absorptions in the range 500-600 nm has been ascertained from the analysis of the vertical transitions and the MOs involved in the transitions (Table 3). It must be highlighted that single-electron excitations originate from occupied MOs with metal character, lower in energy compared to the HOMO, to the LUMO orbital (which is the orbital Cu-d x 2 -y 2, bearing the unpaired electron). The MOs involved in the most relevant transitions are depicted in Fig. 8.

Applications: calculation of the spectrum of Cu(II)-PrP adducts
Prion proteins (PrP) are cellular proteins that undergo a misfolding in the C-terminal region (128-231) to scrapie form (PrP Sc ). In this latter form, the α-helical folded PrP gives β-sheet and generates aggregates [108]. The misfolding is related to several diseases known as transmissible spongiform encephalopathies, which include bovine spongiform encephalopathy in cattle, scrapie in sheep and goat and Creutzfeldt-Jakob disease in humans [109,110].
In this study, two sites were considered with the aim to calculate the electronic absorption spectrum: that containing the amino acids 94-97 modeled by the tetrapeptide Ac-GTHS-NHCH 3 , and that 108-112, modeled by the pentapeptide Ac-NMKHM-NHCH 3 .
The interaction of Cu(II) ion with Ac-GTHS-NHCH 3 can be described through the coordination of His96 that favors the deprotonation of the amide groups of His96, Thr95 and Gly94 to give (His96-N, N − , N − , N − ) donor set. The optimized structure is reported in Fig. 9 and the Cartesian coordinates in Table S14 of the Supplementary Material. The spectrum was simulated with BHandHLYP functional and basis set def2-TZVP. The experimental and the calculated maxima are at 555 nm and 585 nm, respectively, with a PD of 5.4%, in line with the result found with the model of HSA (Table 2).
In the system involving Cu(II) and fragment 108-112, two species are detected at pH 7.4: one of them is the adduct with (N − , N − , N − , His111-N) coordination [107,120,121], while for the second one two possible donor sets were proposed, (Met109-S, N − , N − , His111-N) [107] and (CO, N − , N − , His111-N) [120,121]; in both cases N − -s are the deprotonated backbone peptide nitrogens of Lys110 and His111. The two possibilities are given in Fig. 10. The Cartesian coordinates for these three adducts can be found in Tables S15-S17 of the Supplementary Material.
The value of λ max calcd for (N − , N − , N − , His111-N) coordination is very close to λ max exptl (555 nm vs. 550 nm with a PD of only 0.9%, Table 2). This allowed us to confirm the presence of this adduct in aqueous solution. Concerning the second species, two different calculations were carried out on the optimized structures with (Met109-S

Conclusions
The results presented in this study confirm and offer a new evidence of the potential of computational techniques in bioinorganic chemistry. In particular for copper, an element with many biological and pharmacological properties, both EPR spectra [13] and now electronic absorption spectra can be successfully predicted. Using B3LYP functional for EPR and BHandHLYP for electronic spectra, A z and λ max are predicted with MAPD of 8.6% and 3.1%, respectively. It is noteworthy   that, for Cu(II) complexes, the prediction of the optical spectra is better than that of electron paramagnetic resonance spectra, a finding which is in contrast with results for other metal ions such as V(IV) [11,12]. Thus, the assessment that NMR and EPR spectra can be calculated successfully, while the prediction of optical spectra is possible only qualitatively for metal complexes, must be reconsidered. In addition to the results presented here, we have recently demonstrated that the absorption spectra of square planar Ni(II) species can be predicted with a MAPD of 1.5% in the visible region using HSE06 and MPW1PW91 functionals [57]. Up today, there are no rules to determine the protocol (basis set and, mainly, functional) which gives the best agreement with the experimental data. In this situation, the computational procedure must be validated using a benchmark which includes as many as possible metal complexes with different geometries, coordination numbers and electric charges. The results of this manuscript suggest that functionals BHandHLYP and M06 perform much better than HSE06 and MPW1PW91 in the case of Cu(II); in contrast, with square planar Ni(II) complexes, HSE06 and MPW1PW91 are the best, while BHandHLYP and M06 give worse predictions [57]. This effect can be explained by considering the fact that the open shell systems, such as Cu(II) or V(IV) O complexes, are quite sensitive to the amount of Hartree-Fock exchange included in the functional, as already reported earlier [122]. Similar results were obtained for the prediction of hyperfine coupling constant in the EPR spectra. For example, for vanadium in the oxidation state +4, it was shown that the functional BHandHLYP is recommended for oxidovanadium(IV) complexes [11], while for 'bare' non-oxidovanadium(IV) species the double hybrid B2PLYP functional must be used [12]. As pointed out in a number of papers, the effect of the basis set is less important than the functional, and a triple-ζ polarized (for example, def2-TZVP) is enough to reach good predictions. It needs to be emphasized that, in this study, the functionals BHandHLYP and M06 outperformed B1LYP, which gave the best results for a limited number of Cu(II) species [56]. Accordingly, the latter one is not recommended if an accurate prediction is desired: specifically, the value of MAPD is 3.1% and 3.7% with BHandHLYP and M06, respectively, and as high as 12.4% with B1LYP. Finally, it must be highlighted that the procedure validated in this paper can be applied to predict the transitions and the electronic absorption spectra of Cu(II)-bioligand adducts. The results could be useful to confirm the geometry of a hypothesized adduct through the comparison between the experimental and calculated spectra and, in general, to relate the experimental spectral data to a specific structure of a metallospecies.

Declaration of competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.