Elena
Molteni
*ab,
Giuseppe
Mattioli
a,
Paola
Alippi
a,
Lorenzo
Avaldi
a,
Paola
Bolognesi
a,
Laura
Carlini
a,
Federico
Vismarra
cd,
Yingxuan
Wu
cd,
Rocio Borrego
Varillas
d,
Mauro
Nisoli
cd,
Manjot
Singh
e,
Mohammadhassan
Valadan
ef,
Carlo
Altucci
ef,
Robert
Richter
g and
Davide
Sangalli
ab
aIstituto di Struttura della Materia-CNR (ISM-CNR), Area della Ricerca di Roma 1, Via Salaria km 29.300, CP 10, Monterotondo Scalo, Roma, Italy. E-mail: elena.molteni@mlib.ism.cnr.it
bDipartimento di Fisica, Universita' degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy
cDipartimento di Fisica, Politecnico di Milano, Piazza Leonardo da Vinci, 32, Milano, Italy
dCNR-Istituto di Fotonica e Nanotecnologie, Piazza Leonardo da Vinci, 32, Milano, Italy
eDipartimento di Scienze Biomediche Avanzate, Universita' degli Studi di Napoli Federico II, via Pansini 5, I-80131, Napoli, Italy
fIstituto Nazionale Fisica Nucleare (INFN), Sezione di Napoli, Napoli, Italy
gSincrotrone Trieste, Area Science Park, Basovizza, Trieste, Italy
First published on 19th November 2021
The electronic energy levels of cyclo(glycine–phenylalanine), cyclo(tryptophan–tyrosine) and cyclo(tryptophan–tryptophan) dipeptides are investigated with a joint experimental and theoretical approach. Experimentally, valence photoelectron spectra in the gas phase are measured using VUV radiation. Theoretically, we first obtain low-energy conformers through an automated conformer–rotamer ensemble sampling scheme based on tight-binding simulations. Then, different first principles computational schemes are considered to simulate the spectra: Hartree–Fock (HF), density functional theory (DFT) within the B3LYP approximation, the quasi-particle GW correction, and the quantum-chemistry CCSD method. Theory allows assignment of the main features of the spectra. A discussion on the role of electronic correlation is provided, by comparing computationally cheaper DFT scheme (and GW) results with the accurate CCSD method.
DKPs can be considered as precursors of longer oligopeptides, they have been detected, e.g., in meteorites,5 and, as chiral molecules, they can catalyze enantioselective reactions.6 They may therefore have had a role in the synthesis of the first biomolecules on Earth, and in the emerging of the so-called “homochirality of life”, i.e. the fact that naturally occurring proteins are composed of L-amino acids.5 The formation of cyclo-dipeptides in simulated prebiotic conditions has been studied in ref. 6. The correlation between chirality and self-organization of cyclo-dipeptides has also been investigated in the literature.4
Biomolecules are usually studied in solution, often in water, in order to mimic physiological conditions. The presence of the solvent can significantly affect their chemical and physical properties and indeed the analysis of the effects of type of solvent, pH and temperature on e.g. protein stability and conformation is an interesting field of research. On the other hand, in spite of the challenge to evaporate these molecules intact, gas phase characterization opens the way to determining the intrinsic properties of the molecule without environmental interferences. It also allows to investigate fundamental processes like ultra-fast charge migration via pump and probe experiments.7 Photoinduced charge transfer through molecules is involved in a variety of mechanisms and phenomena, including charge separation in photosynthetic centers, charge generation in organic solar cells, photoprotection and photodamage.8–12
Experimental works on cyclo-dipeptides or linear dipeptides have already been reported in the literature, e.g., the study of the correlation between photoelectron spectra and structures in c-PhePhe, c-TyrPro and c-HisGly13 or in c-GlyGly, c-PhePro and c-LeuPro,14 the laser UV photoionization of the linear TrpTyr and TrpTrp dipeptides15 and the investigation on aromatic cyclo-dipeptides, such as c-TrpTrp, in methanol solution as building blocks for self-assembled biocompatible supramolecular structures with photoluminescence properties.16 The electron mobility and dissociation of several different short linear peptides in the gas phase via UV and IR multiphoton absorption,17 as well as the fragmentation of c-AlaAla, a process relevant to the understanding of the early stages of evolution of life,18 have attracted some interest.
Here a joint experimental and theoretical study of the electronic properties of c-GlyPhe, c-TrpTyr, and c-TrpTrp has been performed. The valence photoemission spectra (PES) in gas phase have been measured using synchrotron radiation and interpreted with the help of ab initio simulations. In our computational protocol, stable geometries of DKP in gas phase are extracted through a large-scale conformational search carried out using tight-binding (TB) simulations. Density functional theory (DFT) calculations are then used to compute the gas phase electronic properties of the most stable conformers. The electronic energy levels and the resulting densities of states (DOS) are compared with the measured PE spectra.
For each of the three dipeptides, a few lowest energy conformers found with this method were then subjected to a geometry optimization run, and, subsequently, to self-consistent electronic structure calculations: both calculations were performed within DFT22,23 using either the Quantum ESPRESSO (QE)24,25 or the ORCA26 suites of programs. There is a generally close agreement on the energy ordering of molecular configurations between tight-binding and ab initio calculations. A very detailed investigation on this point has been performed in the case of c-TrpTrp, and is discussed in the ESI.†
In the case of QE, simulations of the same molecules have been performed in the framework of a plane-waves/pseudopotential method. Initial structures selected by the CREST algorithm have been accommodated in large cubic supercells (40 a.u., i.e. ≈21.167 Å, side length) and fully optimized. Then, in view of the more computationally demanding GW calculations, also a face-centered cubic (FCC) cell has been used, yielding the same distance between molecules in replicated cells with a smaller volume and, therefore, a lower number of G-vectors. Total energies have been calculated using norm-conserving Troullier–Martins atomic pseudopotentials,27 a plane-wave basis set, and the B3LYP hybrid exchange–correlation functional,28,29 (and also the M062X functional30 in some cases) including the D3 pairwise dispersion correction for van der Waals (vdW) interactions.31 Satisfactorily converged results have been achieved by using cutoffs of 90 Ry on the plane waves and of 360 Ry on the electronic density, respectively. H 1s2 (2 electrons), C 2s2-2p2 (4 electrons), N 2s2-2p3 (5 electrons), and O 2s2-2p4 (6 electrons) electrons have been treated as valence electrons. All of the inner shell electrons are embedded in the pseudopotentials. The Makov–Payne correction to the total energy was also computed and an estimate of the vacuum level was calculated, which allowed us to properly align electronic energy levels.32 For geometry relaxation runs the Broyden–Fletcher–Goldfarb–Shanno (BFGS) scheme was used.33
In the case of ORCA, DFT simulations have been performed in an all-electron localized-basis-set framework.26,34 Kohn–Sham orbitals have been expanded on large, Gaussian-type def2-QZVPP basis sets.35,36 The corresponding def2/J basis has been also used as an auxiliary basis set for Coulomb fitting in a resolution-of-identity/chain-of-spheres (RIJCOSX) framework. Molecular geometries have been fully re-optimized and their properties investigated by using the same B3LYP functional used in the case of QE, including the D3 pairwise dispersion correction.
Reference values for molecular ionization energies have been obtained using ORCA, in the same Gaussian-type orbitals (GTO) framework discussed above but through a different wavefunction approach. As introduced above, the simulations are based on the equation of motion – coupled cluster theory (EOM-CC), which extends a primarily ground-state method like the single-reference coupled-cluster one to the calculation of excited/ionized states as well as to the corresponding excitation/ionization energies.41 Within the ORCA implementation, core and valence ionization energies, as well as electron attachment energies for the mapping of unoccupied electronic states, are calculated by grouping electron pairs in domain-local pairs of natural orbitals (DLPNO), and the chain-of-sphere approximation is used to speed up the calculation of exchange-like integrals with four virtual labels (COSX).42,43 The present implementation makes the treatment of large systems (2196 basis functions in the case of the largest c-TrpTrp) possible, with a loss of accuracy lower than 0.01 eV with respect to canonical algorithms.
Fig. 1 Experimental photoemission spectra of the c-GlyPhe, c-TrpTyr, and c-TrpTrp dipeptides. The feature at −12.6 eV in the c-TrpTrp spectrum is due to residual water in the sample. |
All amino acids except glycine (Gly) have a chiral center: each chiral molecule can exist in two chemically indistinguishable forms, enantiomers, labeled “S” and “R”, which can only be identified through their different interaction with circularly polarized electromagnetic radiation, yielding opposite circular dichroism spectra. Out of the possible enantiomers or diastereoisomers resulting from the combinations of the “S” or “R” forms of the two constituent amino acids, we have performed our calculations – for each of the three investigated dipeptides (c-GlyPhe, c-TrpTyr, c-TrpTrp) – on the form on which experiments were performed, i.e. the “S” enantiomer for c-GlyPhe, and the “S,S” diastereoisomer for both c-TrpTyr and c-TrpTrp.
TB ΔE (mHa) | TB pop. (%) | ORCA ΔE (mHa) | |
---|---|---|---|
B3LYP + vdW | |||
conf1 | 0.00 | 84.3 | 0.00 |
conf2 | 2.40 | 6.6 | 3.47 |
conf3 | 2.87 | 4.0 | 3.92 |
We use cyclo(GlyPhe), which is the smallest of the three peptides under study, to compare the electronic properties with different levels of approximation for the exchange–correlation (xc) energy: HF, B3LYP and M062X, with the EOM-CCSD method, and with the GW approximation. This was done in order to investigate the sensitivity of these properties to the chosen method and xc functional, and to find the method and/or functional which yields the best trade off between agreement with experimental data and computational cost. In particular, we focused on the comparison between measured and theoretically estimated – with different approaches - PES.
We inspect the electronic densities of states (DOS) for the GlyPhe conformer 1 (see Fig. 3, panel (a)) resulting from DFT B3LYP (magenta curve in the figure) and HF (yellow curve) occupied electronic levels, and we compare the results with those obtained by calculating ionization energies using the EOM-CCSD method (blue) and the GW scheme (green curve). The latter are used as theoretical reference, and compared with the gas-phase experimental photoemission spectrum (black) of GlyPhe. In the EOM-CCSD DOS reported here each energy root has been weighed with its percentage contribution from single transitions, while in the GW simulation the quasi-particle poles εn are weighed by the renormalization factor Zn.
We can observe that the B3LYP calculation provides a better description of the energy distribution of occupied electronic states with respect to the HF one, provided that a rigid shift of −2.5 eV is applied‡ to the B3LYP levels in order to align the B3LYP HOMO level with the EOM-CCSD at −9.25 eV and thus to the experimental IE. We also report that the B3LYP IE obtained via total energy difference is EN−10 − EN0 = 8.57 eV, while the corresponding value obtained using the meta-GGA global hybrid M062X functional is 9.25 eV, thus matching the CCSD value. HF tends to overbind all main features of the PE spectrum by ≈−2 eV. This is why a rigid shift in the opposite direction is applied to the plot of the occupied HF levels. However, apart from this rigid shift, the agreement between the distribution of HF levels and the PE spectrum is not good, and HF tends to overestimate the energy splitting of the levels close to the HOMO of ≈+2 eV. As a consequence, the unshifted εHFHOMO = 9.29 eV is very close to the CCSD result, but this is the result of two errors with opposite sign. For better clarity the unshifted HOMO/IE are reported in Table 2 and later discussed.
The EOM-CCSD method, taking into account dynamical electronic correlation, is able to yield a PE spectral function (blue curve) in good agreement with the experimental PE spectrum, despite starting from the Hartree–Fock calculation which by itself badly describes the energy distribution of electronic levels. The EOM-CCSD PE energies and the B3LYP DOS are in a surprising and satisfactory agreement with each other in the analyzed energy range, although the CCSD scheme is numerically much more demanding. The GW corrections (green curve) on top of the B3LYP energies lead to a final result which is even closer to the EOM-CCSD one, in particular in the width of the first DOS feature (the one at ≈−10 eV) and in the detailed shape of the various peaks/features in the −11 eV to −19 eV range. GW also partially cancels the underestimation of the HOMO binding energy. A residual underestimation of −1.47 eV remains. The GW levels in Fig. 3 are then shifted by such value. Finally GW and EOM-CCSD also agree for the renormalization of the PE peaks intensity, determined by the percentage contribution from single transitions in EOM-CCSD and by the renormalization factor Zn in GW. In the EOM-CCSD approach, contributions from double transitions are negligible up to ≈−19 eV (an energy close to the double ionization threshold49), and the mixing between single and double transitions is very small in the whole analyzed energy range, (i.e., up to ≈−25 eV). With the exception of two energy roots, with percentage contributions from single transitions of 80%, ≈10% respectively, the contribution from single transitions to a given energy root is either larger than 90%, or below 0.5%. Similarly the renormalization factors Zn of the GW calculation are all above 0.75 up to ≈−23 eV.
In panel (b) of Fig. 3 we compare DFT B3LYP occupied energy levels and DOS of the two lowest energy conformers of the GlyPhe peptide. Although the geometries of these two conformers are quite different (Fig. 2), the energy positions of their electronic levels (inset of Fig. 3(b)) display minor mutual differences only. Once a Lorentzian broadening is applied to these electronic levels in order to obtain the densities of states of the two conformers, these latter appear hardly distinguishable from each other, being both in good agreement with the features in the PE spectrum.
We now focus on empty electronic levels. The inset of Fig. 3, panel (a), shows the position of unoccupied electronic energy levels of GlyPhe conformer 1 obtained within Hartree–Fock (HF) and B3LYP, compared with the inverse PE energies obtained within the GW method on top of B3LYP and within the EOM-CCSD method. For the GW method we only consider the correction to the QE B3LYP bound empty states. No energy shift is applied to the empty states besides the vacuum level correction computed within HF (0.18 eV on the cubic cell) and B3LYP (0.18 eV on the cubic cell; 0.25 eV on the FCC cell). This is somehow equivalent to assume that the shifts applied to the occupied states correspond to an opening of the band gap, i.e. the energy difference between the IE and minus the EA. EOM-CCSD gives a negative EA, i.e. the first root has positive energy of 1.15 eV (fifth column in Table 2). In contrast B3LYP gives εKSn < 0 for n = LUMO and 5 bound empty states in total. The B3LYP EA obtained via total energy difference is instead negative EN0 − EN+10 = −0.81 eV. At the HF level the first unoccupied state localized on the molecule, obtained from the simulation with localized basis set, is at εHFLUMO ≈ 2.3 eV. We refer to such state as a resonance; in the ESI† we show the spatial localization of this state. In the figure inset, HF resonances from ORCA (with localized basis set) are highlighted in red in the dense list of states obtained with QE using the plane-waves basis set, which also captures vacuum states. In a similar way, B3LYP empty levels obtained with ORCA are highlighted in brown in the list of B3LYP empty states obtained with QE. In the B3LYP case we obtain respectively five (two) bound empty states with QE (ORCA), but without a significant variation in the energy of the first unoccupied state (−0.75 eV with QE B3LYP, −0.70 eV with ORCA B3LYP). The GW correction to the bound B3LYP empty levels also gives a negative EA with εGWLUMO = 1.21 eV, quite close to the EOM-CCSD result. GW corrections to the delocalized vacuum states are not considered. Similarly to the occupied levels case, in the EOM-CCSD approach, the contribution from single transitions is larger than 90% for all the inverse PE energy roots, and in GW the renormalization factors Zn are larger than 0.9 for all the computed empty levels.
HF eigen. | B3LYP Etot | M062X Etot | GW eigen. | CCSD roots | Exp. | |
---|---|---|---|---|---|---|
IE | 9.29 | 8.57 | 9.25 | 7.94 | 9.25 | 9.54 ± 0.02 |
-EA | 2.30 | 0.81 | 0.91 | 1.21 | 1.15 | |
gap | 11.59 | 9.38 | 10.16 | 9.15 | 10.4 |
We underline that in the positive energy range it is not possible to compare the energy levels distribution between the different approaches due to the presence of delocalized states from the continuum. Continuum states are captured when using the plane-waves basis set. On the other hand the approach based on localized basis set only captures resonances, although missing their delocalized contribution (resonances have positive energy, as a consequence they must be in part delocalized) and thus possibly introducing an error in the resulting energy. A good description of resonances using a plane-waves basis set would require very large boxes or possibly some lifetime to accelerate convergence. Still the comparison between results with the two basis sets is interesting to be investigated also with the finite boxes we use (see ESI†).
The overall picture is that GlyPhe has a negative EA of 1.15 eV, a IE of 9.25 eV and a gap, or difference between IE and EA, of 10.4 eV. In Table 2 we report the IE and the EA from the different approaches, assuming that HF and B3LYP eigenvalues for the HOMO and the LUMO are a zero order approximation. B3LYP eigenvalues, not shown in the table, underestimate both the IE (6.91 eV) and EA (0.75 eV), and as a consequence the gap, while HF does the opposite.50 B3LYP and M062X total energies and GW eigenvalues still underestimate the gap and the IE and EA, but the error is much smaller. In particular, the meta-GGA global hybrid M062X functional yields vertical IE and EA values in better agreement than B3LYP with those obtained with the CCSD method, and therefore a IE value which is closer to the experimental one (also vertical). The different results are mostly affected by how the electron exchange or Fock term enter in the description. HF has the full exchange term which opens too much the electronic gap, while B3LYP has only a fraction of the exchange which is not enough. Both GW and CCSD add extra terms which lead to a screening of the exchange (in GW via the RPA screening and in CCSD via the inclusion of terms which mimic the IP screening) leading to the correct inclusion of this effect. The dynamical corrections included both within GW and CCSD do not play a significant role since all poles have dominant single character (/are well described within the QP approximation). However CCSD also includes extra terms, neglected in GW, that give a sizable effect. These extra terms are known to be important in small molecules, but their relevance usually decreases while increasing the system size (up to the limit of extended systems where GW gives a very good physical description). Finally, while the use of screened exchange and extra terms is important for a correct electronic gap, the energy levels distribution is much less sensitive as already discussed.
In order to shed light on similarities between the EOM-CCSD on top of HF methods and the GW on top of B3LYP scheme, we also analyzed the spatial localization of some B3LYP occupied and empty electronic states in the energy region of the frontier orbitals, encompassing the gap, as shown in Fig. 4 for the GlyPhe conformer 1. Each of the EOM-CCSD ionization energies can be interpreted as a linear combination of ionizations of the subsiding Hartree–Fock levels. We thus compare the spatial localization of the HOMO B3LYP bound electronic states with that of the weighed superposition of the HF states mostly contributing to the EOM-CCSD ionization energies.
As shown in Fig. 5, the agreement is very good with the first EOM-CCSD ionization energy mixing HOMO and HOMO−2 HF orbitals to lead to a final wave-function which is very similar to the B3LYP HOMO. We did the same analysis considering also the following 2 EOM-CCSD PE energies which we label here IE-1 and IE-2 (ESI†). IE-1 mixes HOMO, HOMO−1 and HOMO−2 HF orbitals leading to a wave-function very close to the B3LYP HOMO−2 state. IE-2 mixes HOMO−5, HOMO−2, HOMO−1 and HOMO HF orbitals, leading to a wave-function similar (but with some differences) to the B3LYP HOMO−5 state. In conclusion, from the analysis of both electronic densities of states and spatial localization of molecular orbitals, we can thus conclude that DFT B3LYP calculations provide a reasonably good description of the electronic properties, if compared to the EOM-CCSD approach, but at a much lower computational cost. This is why we take the liberty of using the wording “molecular orbitals” also for the B3LYP ones.
Fig. 5 GlyPhe conformer 1: Hartree–Fock orbitals (left) contributing to the CCSD ionization energy (center), and the B3LYP HOMO (right). |
In Fig. 6 we show the five lowest energy conformers of the “S,S” diastereoisomer of the cyclo(TrpTyr) peptide as obtained by a tight-binding conformational search using the CREST code. Conformers are labeled (1, 2,…) according to their tight-binding energy ordering (column 1 of Table 3). For TrpTyr the energy ordering of the conformers obtained through the tight-binding run changes upon performing a DFT B3LYP geometry relaxation, with either the ORCA or the QE code, although their geometries remain substantially unaltered. In particular, the energy ordering obtained after performing ORCA B3LYP relaxations (column 5 of Table 3) is the following (using conformer labels based on their tight-binding energy ordering): 1, 5, 4, 2, 3. Thus conformer 1 remains the lowest energy geometry, while conformers 5 and 4, in which (as in conformer 1) the side-chain rings of the tyrosine and tryptophan amino acids are closer to each other, become energetically favored with respect to conformers 2 and 3. Based on these results, we decided to focus on conformers 1, 5 and 4 for further calculations on electronic properties.
Fig. 6 Geometries of the five lowest energy conformers of the cyclo(TrpTyr) peptide. Color codes as in Fig. 2. |
TB ΔE (mHa) | TB pop. (%) | ORCA ΔE (mHa) | |||
---|---|---|---|---|---|
B3LYP | vdW | Sum | |||
conf1 | 0.00 | 30.4 | 0.00 | 0.00 | 0.00 |
conf2 | 0.38 | 20.3 | 0.32 | 2.52 | 2.84 |
conf3 | 0.57 | 16.6 | 0.59 | 2.57 | 3.16 |
conf4 | 1.14 | 0.09 | 0.26 | 1.79 | 2.06 |
conf5 | 1.49 | 0.06 | 1.49 | −0.49 | 1.00 |
Intramolecular weak interactions such as the van der Waals forces can have an important role in molecule stability, and the chosen method for their computational treatment can alter the energy ordering of molecule conformers. π–π and CH–π interactions can contribute to stabilizing specific conformers, in molecules containing aromatic rings. Indeed, out of the five low-energy geometries obtained through the tight-binding conformational search, conformers 1, 5 and 4 appear to be stabilized to some extent by π–π interactions between the phenol ring of tyrosine and the indole ring of tryptophan, while the geometries of conformers 2 and 3 suggest a possible CH–π interaction between the Cβ group of tyrosine and the indole ring of tryptophan. Interestingly, the dispersive vdW contribution (column 4) to the computed total energy of c-TrpTyr has an important role in determining the DFT B3LYP energy ordering of its conformers (see column 5 of Table 3). We recall that a similar vdW correction is present in the TB simulation. More in detail: the ORCA B3LYP calculation (third column in the table) tends to energetically favor conformer 4 with respect to the energy ordering obtained from the tight-binding conformational search (first column) as long as we do not consider the vdW contribution to the DFT total energy. The vdW contribution to the DFT energy (fourth column) is negative for conformer 5 only – and zero or positive for the other conformers – thus strongly favoring conformer 5. As a final result, the energy ordering of conformers according to the total B3LYP energy including the vdW contribution is, as already mentioned, 1, 5, 4, 2, 3. The same final energy ordering is also obtained using the D4 vdW dispersion correction51,52 in the ORCA B3LYP calculation, instead of the D3 discussed so far: the only remarkable difference is that the ORCA calculation with the D4 dispersion correction yields a positive vdW contribution for conformer 5; this is not sufficient to change the energy ordering of the conformers, but it brings the energy of conformer 5 very close to that of conformer 4 (differing by 1 meV only).
In Fig. 7, panel (a), we show the electronic densities of occupied states for TrpTyr conformer 1 (i.e. the most populated one according to both tight-binding conformational search and B3LYP DFT) obtained either from DFT B3LYP (magenta curve), or from a valence photoemission spectrum (blue) obtained through a EOM-CCSD calculation. Also for this cyclo-dipeptide, as already observed for GlyPhe, the DFT B3LYP and EOM-CCSD densities of states are in good agreement with each other and with the experimental photoemission spectrum (black curve), once a −2.5 eV energy shift is applied to the B3LYP DOS. Also here (as in the GlyPhe case) the EOM-CCSD energy roots have been weighed with their percentage contribution from single transitions. For TrpTyr this percentage is higher than 90 in all the analyzed energy range, which is indeed shorter than the one considered for GlyPhe (since considering the same number of energy roots for a larger molecule).
As for the sensitivity of the DFT B3LYP densities of states to molecule conformation, panel (b) of the same figure shows the DFT B3LYP DOS for the three lowest energy conformers of TrpTyr according to DFT B3LYP, i.e. conformers 1, 5 and 4 according to tight-binding labeling. Also for TrpTyr we can only observe minor differences among the densities of states of different conformers, such as a shoulder at ≈8 eV ionization energy, which is most pronounced in conformer 5, less in conformer 1, and almost absent in conformer 4. As such, also for TrpTyr the densities of states of all the three analyzed conformers are in good agreement with the experimentally measured photoemission spectrum of the molecule (black curve in the figure).
The spatial localization of electronic states of TrpTyr lowest energy conformer (conformer 1) is shown in Fig. 8. Frontier orbitals only are considered, calculated within DFT B3LYP with either the plane-wave pseudopotential code QE or the localized basis all-electron code ORCA. As expected, the two codes give almost identical results, apart from a switch between HOMO−2 and HOMO−1 which are very close in energy (Table 4).
QE (eV) | ORCA (eV) | |
---|---|---|
LUMO+3(*) | −0.3325 | +0.2208 |
LUMO+2(*) | −0.4689 | +0.0528 |
LUMO+1 | −0.6307 | −0.3518 |
LUMO | −0.7573 | −0.5886 |
HOMO | −8.2152 | −8.0903 |
HOMO−1 | −8.7344 | −8.5768 |
HOMO−2 | −8.7602 | −8.6369 |
From Table 4 we can also observe that the ORCA B3LYP calculation yields two bound empty states only (column 2); on the other hand, in the QE B3LYP calculation we find seven bound empty states, i.e. three further states in addition to the ones appearing in column 1 of the table. As already observed for the c-GlyPhe dipeptide, in the ORCA B3LYP calculation with localized bases empty levels are generally less “dense” in energy with respect to those obtained from the plane-waves QE B3LYP calculation, while the energy distribution of occupied states is very similar with the two approaches. This trend can be rationalized by recalling that DFT calculations with plane-waves basis sets are able to also capture continuum(-like) delocalized empty electronic states, as opposed to localized basis codes, which only capture resonances, possibly introducing an error in their energies.
Also for the “S,S” form of TrpTrp, as already found for TrpTyr, the energy ordering of conformers as obtained through the initial TB conformational search is not maintained in QE DFT B3LYP (see columns 1 and 3 of Table 5). In this case not even the lowest energy conformer is conserved: labeling conformers according to their TB energy order, within DFT B3LYP with the QE code the lowest energy geometry is conformer 3, followed by conformers 2, 1, 4 and 5 (their geometries are shown in Fig. 9). Also for this dipeptide, some of the five TB lowest energy conformers, i.e. conformers 2, 1 and 5, appear likely to be stabilized to some extent by π–π interactions between the aromatic rings of the two constituent amino acids – here two tryptophan indoles – and others (conformers 3 and 4) suggest CH–π interactions between the Cβ group of one tryptophan and the indole ring of the other.
TB ΔE (mHa) | TB pop. (%) | QE ΔE (mHa) | |
---|---|---|---|
B3LYP + vdW | |||
conf1 | 0.00 | 35.1 | 0.00 |
conf2 | 0.52 | 20.3 | −0.15 |
conf3 | 1.21 | 19.5 | −0.41 |
conf4 | 1.51 | 14.3 | 1.74 |
conf5 | 2.53 | 4.8 | 5.00 |
Fig. 9 Geometries of the five lowest energy conformers of the cyclo(TrpTrp) peptide. Color codes as in Fig. 1. |
A more detailed analysis (see ESI†), exploring the dependence of the DFT total energy ordering of the five lowest energy conformers from the tight binding conformational search on the computational details of the DFT calculation, such as cell size, localized bases (ORCA) vs. plane waves (QE), type of vdW treatment, exchange–correlation functional, pseudopotentials, has shown that the energy ordering among conformers 1, 2 and 3 can change in some cases, but these three conformers are always lower in energy than conformers 4 and 5. Moreover, the energy differences among conformers 1, 2 and 3 are in most cases within the value of kBT at room temperature, and lower than those between these three conformers and the other two (conformers 4 and 5). This supports our choice to analyze conformers 1, 2 and 3 as lowest energy/most populated conformers of this peptide.
The agreement between the electronic densities of states obtained either from DFT B3LYP energy levels (green curve in panel (a) of Fig. 10) or from a valence photoemission spectrum (blue curve) obtained through a EOM-CCSD calculation is quite good, and both curves (obtained for the lowest energy conformer according to DFT B3LYP) reproduce reasonably well the experimental photoemission spectrum (black) of TrpTrp in the analyzed energy range. The EOM-CCSD energy roots have been weighed with the percentage contribution from single transitions, which is also in this case (as in TrpTyr) higher than 90 for all the obtained ionization energies.
A comparison between DFT B3LYP densities of states of the three lowest energy conformers of the TrpTrp dipeptide (panel (b) of Fig. 10) shows a negligible conformational sensitivity of the DOS also for this molecule: the three DOS curves are almost indistinguishable, despite the rather different geometries of the analyzed conformers (Fig. 9). Also in this case, the calculated density of states alone would not be sufficient in order to draw any conclusions on the presence of the different possible conformers in the experimental sample.
In Fig. 11 we report the spatial localization of QE DFT B3LYP highest occupied and lowest unoccupied electronic states of the TrpTrp lowest energy conformer (conformer 3). Interestingly, molecular orbitals ranging from HOMO−3 to LUMO+1 are alternately localized on the indole ring of either of the two tryptophan amino acids. This is likely due to structural differences responsible for a different chemical environment of the two otherwise identical indole groups. Molecular orbitals localized on the diketopiperazine ring resulting from peptide cyclization appear at lower energies (HOMO−4 to HOMO−7 levels).
Indeed, the IE is given by the difference between EN0 (the ground state total energy with N electrons) and EN−10 (the ionized state with N − 1 electrons), so that, when comparing the IE values of different molecules, the observed trends will depend non-trivially on the combined effect of variations in EN0 and in EN−10. In particular, the trend observed for the IE values of the three investigated cyclo-dipeptides c-GlyPhe, c-TrpTyr, c-TrpTrp, i.e. lower IEs for larger aromatic system, can be rationalized by considering that a larger aromatic group will make the ionized molecule more stable (lowering of EN−10) due to charge delocalization. When comparing the IEs of these cyclo-dipeptides with those of the corresponding single amino acids, instead, the above-mentioned stabilizing effect of the DKP ring on the neutral molecule (lowering of EN0) can explain the higher values of IE found for cyclo-dipeptides with respect to single amino acids.
Similar trends are predicted by B3LYP and CCSD calculated densities of electronic states (colored curves in Fig. 3, 7 and 10). From the numerical simulations we can analyze the orbitals which determine the ionization potential and contribute to the first peak in the DOS, in the energy range between 8–10 eV. To this end we plot in Fig. 12, for the case of the TrpTyr, the electronic density obtained from the orbitals in the selected energy range. Similar informations can be obtained for GlyPhe and TrpTrp, by looking at the orbitals in Fig. 4 and 11. From the plots we clearly see that the ionization potential is determined by “pz-like” orbitals in the aromatic rings and that these orbitals provide a large contribution to the whole first feature (consisting of two substructures, which merge for GlyPhe only) in the DOS. The second substructure in this first DOS feature shows contributions also from the orbitals on the DKP ring (highlighted in blue in Fig. 12). In the case of phenylalanine, the energy of the ionization of the phenyl ring approaches the one of the DKP ring,14 and the DOS shows a set of six states very close to each other. As a consequence the highest occupied orbitals of the phenyl side chain are mixed with the DKP ones, and produce the single broader feature observed in the experiment as well as in the calculated DOS for GlyPhe. After these frontier orbitals, theory predicts a region void of states consistently with the gap observed in the experimental spectra. Then a region with a high density of states (29 in the case of TrpTyr as shown in Fig. 12) follows. These electronic states lie in the planes of the aromatic rings, primarily due to the network of sp2 orbitals of the carbon atoms in these rings.
Fig. 12 Bottom left: [−20, −5] eV part of the plot (see Fig. 7 panel (a)) showing the B3LYP calculated DOS of TrpTyr conformer 1. Other panels: plots showing the density of states, for TrpTyr conformer 1, integrated on selected energy ranges, corresponding to specific regions of the DOS plot. |
All the above considerations are based on the analysis of the B3LYP orbitals, which is simple to perform. The good quality of the B3LYP description is confirmed by the agreement with the experimental data, but also by the more refined calculations using the CCSD approach. The latter confirm that the signatures of PES are mostly related to single particle features, with only a very weak satellite appearing in GlyPhe, i.e. the smallest of the molecules. This is in contrast to what often happens in small molecules where mixing with double and higher order excitations can have an important role. The theoretical description is reminiscent of what happens in extended systems which involve mostly s and p electrons. In such cases DFT provides a reasonable description of the electronic density, and only the eigenvalues need to be corrected. Indeed, for GlyPhe only, we also performed calculations within the state of the art DFT + GW scheme, usually employed to describe extended systems, and found good agreement with the CCSD scheme. The key ingredient of the GW scheme is the screening, whose role grows in importance with the system size, while making other correlation effects included in schemes like CCSD less important. GW is not on top of the CCSD results, but it is reasonably close, and it also confirms that the PES can be well described in terms of single-particle features.
This can be further understood in terms of the available phase space. When the system size is increased, so is the number of single particle orbitals, which then provide enough phase space to reach a good description of the electronic properties of the molecules. Observations related to the system size can be done also for the HOMO–LUMO gap and its trend with the molecule size. The HOMO–LUMO gap provides an indication of the stability of the activated molecule towards further chemical reactions. In our case we found the HOMO–LUMO gap tends to decrease as the size of the side chain group becomes larger. This is similar to what happens for example in solid state physics when moving from smaller to bigger nano-structures, to the bulk limit. The bigger the system, the more available phase space to relax, the smaller the electronic gap.
Moreover, we also used the molecules to analyze different computational schemes, in particular comparing an accurate quantum chemistry scheme, based on CCSD, with the state of the art approach for extended systems, GW on top of DFT. The goal was to validate the computationally cheaper DFT results; but also to draw considerations related to the importance of dynamical corrections (i.e. multiple excitations or satellites), electronic screening, and the system size. We consider CCSD as the reference scheme for these molecules, but we also observe that their size is approaching the point where the DFT + GW scheme can also be considered a good alternative.
Footnotes |
† Electronic supplementary information (ESI) available: Text and figures about unbound electronic states in the continuum: Fig. S1: calculated DOS for both occupied and empty electronic states of c-GlyPhe, obtained with different computational methods; Fig. S2: spatial localization of “resonance” states and of completely delocalized states. Text and tables (Tables S1–S4) on the energy ordering of c-TrpTrp conformers and its dependence on several computational details. Fig. S3 ans S4: spatial localization of EOM-CCSD and of B3LYP orbitals of c-GlyPhe. Tables S5 and S6: calculated vertical and adiabatic ionization energies and electron affinities for the cyclo-dipeptides under study. See DOI: 10.1039/d1cp04050b |
‡ This is an extra shift on top of the vacuum level correction which we compute within the Makov–Payne and Martyna–Tuckerman (MT) schemes. For the box size used both schemes give ≈0.25 eV. |
This journal is © the Owner Societies 2021 |