# Excited states computation of models of phenylalanine protein chains: TD-DFT and composite CC2/TD-DFT protocols

Marine Lebel, Thibaut Very, Eric Gloaguen, Benjamin Tardivel, Michel Mons, Valérie Brenner

► To cite this version:

Marine Lebel, Thibaut Very, Eric Gloaguen, Benjamin Tardivel, Michel Mons, et al.. Excited states computation of models of phenylalanine protein chains: TD-DFT and composite CC2/TD-DFT protocols. International Journal of Molecular Sciences, 2022, 23 (2), pp.621. 10.3390/ijms23020621. cea-03524523

HAL Id: cea-03524523

https://cea.hal.science/cea-03524523

Submitted on 13 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Excited States Computation of Models of Phenylalanine Protein Chains: TD-DFT and Composite CC2/TD-DFT Protocols

Marine Lebel, Thibaut Very †, Eric Gloaguen, Benjamin Tardivel, Michel Mons D and Valérie Brenner *ID

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Lebel, M.; Very, T.; Gloaguen, E.; Tardivel, B.; Mons, M.; Brenner, V. Excited States Computation of Models of Phenylalanine Protein Chains: TD-DFT and Composite CC2/TD-DFT Protocols. *Int. J. Mol. Sci.* **2022**, *23*, 621. https://doi.org/10.3390/ijms23020621

Academic Editor: Isabelle Callebaut

Received: 16 November 2021

Accepted: 3 January 2022

Published: 6 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

LIDYL, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette, France; laiguillon32@hotmail.fr (M.L.); thibaut.very@idris.fr (T.V.); eric.gloaguen@cea.fr (E.G.); benjamin.tardivel@cea.fr (B.T.); michel.mons@cea.fr (M.M.)

* Correspondence: valerie.brenner@cea.fr

† Present address: IDRIS-CNRS, Orsay CEDEX F-91403, France.

**Abstract:** The present benchmark calculations testify to the validity of time-dependent density functional theory (TD-DFT) when exploring the low-lying excited states potential energy surfaces of models of phenylalanine protein chains. Among three functionals suitable for systems exhibiting charge-transfer excited states, LC-ωPBE, CAM-B3LYP, and ωB97X-D, which were tested on a reference peptide system, we selected the ωB97X-D functional, which gave the best results compared to the approximate coupled-cluster singles and doubles (CC2) method. A quantitative agreement for both the geometrical parameters and the vibrational frequencies was obtained for the lowest singlet excited state (a ππ* state) of the series of capped peptides. In contrast, only a qualitative agreement was met for the corresponding adiabatic zero-point vibrational energy (ZPVE)-corrected excitation energies. Two composite protocols combining CC2 and DFT/TD-DFT methods were then developed to improve these calculations. Both protocols substantially reduced the error compared to CC2 and experiment, and the best of both even led to results of CC2 quality at a lower cost, thus providing a reliable alternative to this method for very large systems.

**Keywords:** electronic excitations; protein modeling; molecular simulation; time-dependent density functional theory; photophysics

## 1. Introduction

Understanding the electronic dynamics of bio-relevant systems usually requires the exploration of potential energy surfaces (PES) of their low-lying excited states, in particular in order to localize non-equilibrium geometries such as conical intersections (CI) [1–5]. This can be accomplished by performing non-adiabatic dynamic simulations considering all the low-lying significant excited states for the dynamics of these systems. A good compromise between accuracy and computational times in simulations on medium-sized systems is to use a single reference method for electronic structure calculations, providing that the excited states are obviously dominated by single excitations. Among the single reference methods, the time-dependent density functional theory (TD-DFT) [6], which describes excited states within response theory, exhibits a favorable cost-performance ratio, but even to only perform a qualitative exploration of the PES, the functional has to be judiciously chosen according to the nature of the excited states considered [7–11]. In this context, we developed an original, innovative computational strategy in order to document the conformer selective dynamics of capped peptide building blocks (including the phenylalanine (Phe) residue), serving as models of proteins [12–19]. Gas-phase isolation was investigated first, enabling cross-checking between the experimental data and relevant quantum chemistry methods in order to validate the theoretical approach. The challenge in such calculations is that these systems exhibit specificities that orient the choice of the approach: (i) their size (medium-size systems where the smallest one, a capped peptide with one residue, already contains at least thirty atoms); (ii) their lack of symmetry, (iii) their great
flexibility due to the non-covalent interactions that govern their structure; and (iv) their
multiple singlet close-low-lying excited states featuring very different natures, locally
excited states on one or several peptide bonds (nπ\*CO) or on one or several phenyl rings
(ππ\*) and even charge transfer (CT) states from the backbone to a phenyl ring (nπ\*) [17–20].
This justifies the “multiscale” character of our strategy, which consists of, first, non-adiabatic
dynamics simulations within the TD-DFT framework to provide hints about the critical
motions that drive the deactivation and then refine them at two better levels of theory. The
two better levels of theory are (i) the standard approximate coupled-cluster singles and
doubles method (CC2) [21–25] and (ii) a multireference configuration interaction (MRCI)
method [26,27].

The first key point in the development of this strategy was to demonstrate, by compar-
ison with multireference configuration interaction (MRCI) calculations on such systems, the
validity of the CC2 method. The investigated geometries corresponded not only to the equi-
librium geometry of the initially excited state accessible from the Franck–Condon region
but also to the equilibrium geometries of all close-low-lying excited singlet states of these
systems. In addition, relevant geometries along with the energy profile of the deactivation
mechanisms, such as the CI between excited states [17,19], were also considered.

In the second key step, an extensive benchmark of the CC2 method (structure, energetic and vibrational frequencies of the first ππ\* excited state) was performed on a series of capped peptides of increasing size and containing several residues, by comparison with experimental data in order to assess the accuracy of this method [18].

Finally, the third key point that remained to be solved consisted in assessing the ability
of the TD-DFT method to model the PES of these systems qualitatively. This is one of
the objectives of the present work. We first focused on the choice of the functional to be
used for performing the calculations (including energetics as well as the first and second
derivatives of the energy) on a series of capped peptides of increasing size and containing
several different residues. These were chosen for their ability to adopt in the ground state
prototypical secondary structures of proteins [15]. Due to the presence of CT states in
these systems, benchmarking calculations were performed using three long-range corrected
hybrid functionals, LC-ωPBE [28], CAM-B3LYP [29], and ωB97X-D [30] that are suited
for the calculations of such states [8,9,31]. These calculated data were directly compared
with both experimental spectroscopic data and results of the more refined, correlated,
single reference CC2 method [18]. The second objective of this work, moreover, was to
evaluate composite CC2/TD-DFT protocols for the 0–0 excitation energies calculations
and to determine if such protocols can become, for large and very large-sized systems,
an alternative to more refined single reference methods such as CC2, which become not
feasible in this case. Indeed, it has recently been demonstrated on several series of organic
and inorganic molecules that, while the vertical excitation energies are sensitive to the
choice of the method, CC2 or TD-DFT, the difference in zero-point vibrational energy
(ZPVE) corrections between the excited and ground states as well as the excited state
relaxation energy is much less sensitive to it [32–35]. Therefore, it would be desirable to
develop and assess such composite protocols capable of yielding 0–0 excitation energies of
CC2 quality at a lower cost. In this regard, we investigated two types of protocols: the first
one that does not use CC2 for both the geometry optimizations and frequency calculations
and the second one that uses it only for the geometry optimizations.

**2. Results and Discussion**

2.1. *Selection of the Functional in TD-DFT Calculations*

The selection of the functional was carried out in three steps on the reference system, the N-acetyl-phenylalaninyl-amide (Ac-Phe-NH₂, Fa in short), and more specifically on its four conformers (Fa A-D) identified in the experiment (Figure 1) and characterized by CC2 calculations [18,36]. First, at the DFT-D optimized geometry (B97-D2/TZVPP) of the ground state, the TD-DFT/cc-pVDZ five lowest excited states of the four conformers and their vertical excitation energies were analyzed by comparison with those obtained
at the CC2/cc-pVDZ level. Second, the geometry optimization of the lowest excited state accessible from the Franck–Condon region, a singlet $\pi\pi^*$ excited state localized on the phenyl ring, was performed at the TD-DFT/cc-pVDZ level for each of the four conformers and the quality of the optimized geometries obtained was assessed by comparison with the CC2 ones. Finally, for each of the four conformers, the adiabatic ZPVE-corrected excitation energies and the IR signature, the *amide A* vibrational frequencies (NH stretches), of this lowest $\pi\pi^*$ excited state were evaluated at the TD-DFT/cc-pVDZ level and compared to CC2/cc-pVDZ results as well as to experimental data.

### 2.1.1. Nature and Vertical Excitation Energies for the Five Lowest Excited States of Fa Conformers

In the ground state, four conformers (A–D) with two different types of backbone folding were observed and assigned to conformations lying in an energy range of 0–6.5 kJ/mol (DFT-D level, Figure S1a,b, and Table S1). At the CC2 level, whatever the conformers, the first and the fourth excited states were locally $\pi\pi^*$ excited states centered on the phenyl ring (see Table 1, Table S2a–c, and Figure 2). The second and the third were locally $n\pi^*_{CO}$ excited states, each one centered on a peptide bond; the first peptide bond, that of the N-terminal side, and the second peptide bond, that of the C-terminal side. For Fa A, the lowest $n\pi^*_{CO}$ excited state was localized on the first peptide bond, whereas for Fa B-D, which exhibits a folded backbone, the lowest $n\pi^*_{CO}$ excited state was localized on the second peptide bond. The fifth state for Fa A-C corresponded mainly to an $n\pi^*$ charge transfer (CT) excited state involving an electronic charge transfer from the lone pairs of the backbone (one or two peptide bonds) to the phenyl ring. In the case of Fa D, the fifth state was mainly a locally $n\pi^*_{CO}$ excited state (74%), and the first state with a dominant CT character was the ninth state.

<table><caption>**Table 1.** Nature and CC2/cc-pVDZ vertical excitation energies of the five lowest excited states of Fa B at the B97-D2/TZVPP optimized geometry of the ground state.</caption><thead><tr><th>Fa B</th><th>E<sub>vert</sub><sup>[a]</sup> (eV)</th><th>Nature of the State</th><th>NTOs (%, Occupied→Virtual)<sup>[b]</sup></th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>5.259</td><td>ππ*</td><td>1–53%<br>2–47%</td></tr><tr><td>S<sub>2</sub></td><td>5.781</td><td>nπ<sup>*</sup><sub>CO(2)</sub></td><td>1–98%</td></tr><tr><td>S<sub>3</sub></td><td>5.905</td><td>nπ<sup>*</sup><sub>CO(1)</sub></td><td>1–99%</td></tr><tr><td>S<sub>4</sub></td><td>6.534</td><td>ππ*</td><td>1–73%<br>2–25%</td></tr><tr><td>S<sub>5</sub></td><td>6.787</td><td>n<sub>(1,2)</sub> [π<sup>*</sup>π<sup>*</sup><sub>CO(1)</sub>] + ππ*</td><td>1–80%<br>2–12%</td></tr></tbody></table>

<footnote>[a] The theoretical values are given with the number of significant digits obtained in the experiment. [b] In the event of several couples of natural transition orbitals (NTOs), only those whose contribution to the wave function is greater than 10% are reported.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Ac-Phe-NH₂ (Fa) and Ac-X-Phe-NH₂ with X = Gly, Phe, Gln (GFa, FFa, and QFa, respectively) conformers and their experimental 0–0 transition energies (eV). Arrows indicate non-covalent intramolecular bonds (hydrogen bonds (C₅, C₇, and C₁₀), NH-π bonds, and phenyl ring-ring arrangements). Letters A, B, ... refer to the notation adopted in the experimental reports to designate the conformers observed, which correspond to different types of backbone folding and non-covalent interactions (see the Methods section for a detailed description). In the case of several S₀ conformers, which present both similar backbone folding and non-covalent networks, a single quotation mark is put to distinguish them, such as B and B'. In the case of several S₁ conformers whose conformations have been obtained from the geometry optimization of the same S₀ conformer but correspond to different excitation such as excitation via each chromophore, a number is put in subscript to distinguish them, such as A₁ and A₂.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Contours ($\pi\pi^* \pm 0.0015$ au, $n\pi^*_{CO(2)} \pm 0.03$ au, and CT $\pm 0.007$ au) of the CC2 electron density difference between excited and ground states, calculated at the B97-D2/TZPP optimized geometry of the Fa B ground state. A density increase (decrease) is indicated in blue (red).</figcaption></figure>

**LC-ωPBE:** The nature of the four low-lying excited states was well reproduced for all the four conformers. However, among the ten lowest low-lying excited states, no CT excited state was identified for any of the four conformers, and only states with minor CT contributions were identified. In particular, the fifth excited state corresponded to a pure locally $\pi\pi^*$ excited state centered on the phenyl ring for Fa C and mainly to a $\pi\pi^*$ excited state (58%, 48%, and 55%) for Fa A, B, and D (Table S3a–d). We decided then to study the influence of the $\omega$ parameter of the functional on the existence of CT states and, in particular, on the CT contribution in these low-lying excited states. In order to achieve this, we performed a series of calculations at different values of $\omega$ for A and B conformers: the standard value $0.4 \text{ bohr}^{-1}$ and four surrounding values, i.e., 0.45, 0.35, 0.30, and $0.25 \text{ bohr}^{-1}$. The results highlighted that a $\omega$ value of $0.30 \text{ bohr}^{-1}$ was the best compromise allowing to obtain for the fifth excited state a CT contribution in accordance with the CC2 data, without damaging the nature of the other valence excited states, in particular the fourth state, which corresponded to a $\pi\pi^*$ excited state (Table S3e). However, the mean absolute error (MAE) between the CC2 and TD-DFT levels on the vertical excitation energies of these valence states significantly increased for this $\omega$ value (Table S3f).

**CAM-B3LYP and ωB97X-D:** The nature of all excited states (Tables S4a–d and S5a–d) were well reproduced by the two functionals except the fifth state of Fa D, which corresponded mainly to a CT state, in contrast with the locally $n\pi^*_{CO}$ excitation feature found at the CC2 level. The mean absolute (MAE) and mean signed errors (ME) on the vertical excitation energies were similar: 0.21 and −0.13 eV, respectively, for CAM-B3LYP and 0.19 and −0.10 eV, respectively, for ωB97X-D. In addition, if the vertical excitation energies of the lowest $\pi\pi^*$ excited states were overestimated (an MAE of 0.19 eV with an ME of +0.19 eV for both functionals), all the others excited states including the higher $\pi\pi^*$ excited states were underestimated (an MAE and ME of 0.24 and −0.24 eV, respectively, for CAM-B3LYP; 0.19 eV and −0.19 eV, respectively, for ωB97X-D).

## 2.1.2. Geometry Optimization of the Lowest Singlet $\pi\pi^*$, IR Signature, and 0–0 Excitation Energy

In view of the poor results obtained on the nature and vertical excitation energies for the Fa A-D five lowest excited states with the LC-ωPBE functional, the geometry optimization of the lowest excited state (S₁), a $\pi\pi^*$ excited state, was only performed with the CAM-B3LYP and ωB97X-D functionals.

Selected characteristic geometrical parameters of the lowest $\pi\pi^*$ excited state optimized geometry of Fa A-D are reported in Table S6a–c for the CC2, CAM-B3LYP, and ωB97X-D levels, and a comparison of the whole set of geometries is depicted in Figure 3 and Figure S6. Whereas only relatively small deviations from CC2 were obtained with the ωB97X-D functional for the four Fa A–D conformers, the CAM-B3LYP functional gave a significant deviation for Fa A. Indeed, the interactions that govern the arrangement of the backbone above the phenyl ring and, more particularly, the NH...π interaction differed strongly in this conformer as if the dispersion forces were not quite correctly taken into account by the CAM-B3LYP functional. In fact, the NH₂ group was shifted out of the phenyl ring, and the NH...π interaction was strongly reduced. With the ωB97X-D functional, all
the discrepancies for the Fa A-D optimized geometries fell within a range from −8 to 6° for the dihedral angles with a mean absolute deviation (MAD) of 3° and a range from −0.04 to 0.14 Å for the intramolecular distances with a MAD of 0.05 Å. Similar discrepancies were obtained with CAM-B3LYP for the Fa B–D conformers, i.e., within a range from −8 to 6° for the dihedral angles with a MAD of 4° and a range from −0.06 to 0.23 Å for the intramolecular distances with a MAD of 0.10 Å. On the contrary, the geometry of the Fa A exhibited large discrepancies, i.e., into a range from 4 to 26° for the dihedral angles with a MAD of 17° and into a range from −0.43 to 0.75 Å for the intramolecular distances with a MAD of 0.33 Å.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Optimized geometries of the S₁ state for Fa A and B. Comparison of the CC2/cc-pVDZ geometry (blue structures) with that obtained with ωB97XD/cc-pVDZ (red structures) and with CAM-B3LYP/cc-pVDZ (green structures). For the sake of comparison, phenyl rings were overlapped.</figcaption></figure>

On the contrary, the optimized geometry of the ground state of Fa A (S₀, see Table S6a–c) obtained with CAM-B3LYP did not present such large deviations. The discrepancies from 0 to 5° for the dihedral angles with a MAD of 5° and from −0.05 to 0.12 Å for the intramolecular distances with a MAD of 0.07 Å, were similar to those obtained for Fa B–D ground states, i.e., from −9 to 5° for the dihedral angles with a MAD of 4° and from −0.06 to 0.22 Å for the intramolecular distances with a MAD of 0.11 Å. Furthermore, they were comparable to those obtained with ωB97X-D for the Fa A–D ground states, i.e., from −8 to 2° for the dihedral angles with a MAD of 2° and from −0.04 to 0.11 Å for the intramolecular distances with a MAD of 0.05 Å. These results led to a root mean square deviation (RMSD) of the ωB97X-D intramolecular distances relative to CC2 for the Fa A–D ππ* excited state (S₁) conformers almost equal to that obtained for the Fa A–D ground state (S₀) conformers, i.e., 0.06 Å and 0.07 Å, respectively. These very close values were small but larger by one order of magnitude than those obtained for the covalent bonds. Moreover, no difference was observed between RMSD in the ground (S₀) and in the first ππ* excited state (S₁). This was in contradiction with what was reported for the covalent bond lengths by us in this work as well as by others in a benchmark set of 20 medium-sized aromatic organic compounds [37] for which the RMSD(S₁) for these two functionals, despite being small, 0.009 Å, was larger than their RMSD(S₀), i.e., 0.006 Å. These results demonstrate that non-covalent bonds governed by relatively weak interactions, such as dispersion forces, were more sensitive to the method than the covalent bonds but, on the other hand, do not depend significantly on the nature of the state, ground or excited.

In order to go further in the comparison of geometries, the harmonic vibrational frequencies, in particular, the three NH stretch frequencies (NHₚhe, NH₂sym, and NH₂anti.) of the *amide A* region (Table S7) and their shifts between the S₁ and the S₀ states (Table 2) were calculated. Once again, except for Fa A with the CAM-B3LYP functional, the CC2 values were relatively well reproduced by the two functionals, especially if we first consider the order of magnitude as well as the sign of the shift. Indeed, except for CAM-B3LYP Fa A, when a large experimental shift was obtained, the two functionals gave large theoretical shifts close to that obtained by the CC2 method, and a similar behavior was obtained for small shifts. Regarding the comparison with the experiment, DFT/TD-DFT and CC2 predict the shift sign properly but tend to overestimate the frequency shifts, in line with a general overestimate of the harmonic frequencies in the *amide A* region. We then carried out for the series of capped peptides (see the following section) the protocol developed for the CC2 calculations in order to obtain the mode-dependent corrected frequencies from
the optimal harmonic frequency mode-dependent (NH<sub>Phe</sub>, NH<sub>2sym</sub>, and NH<sub>2anti</sub>) linear scaling functions [18].

<table><caption>**Table 2.** Theoretical shifts (harmonic x and mode-dependent corrected y, i.e., x or x/y) of the *amide* A region frequencies of the lowest $\pi\pi^*$ excited state optimized geometry (S<sub>1</sub>) relative to the ground state optimized geometry (S<sub>0</sub>) for the Fa A–D conformers obtained at the CC2/, CAM-B3LYP and ωB97X-D/ /cc-pVDZ levels, together with the corresponding available experimental data (cm<sup>-1</sup>).</caption><thead><tr><th></th><th>Δν<sub>S1/S0</sub> (cm<sup>-1</sup>)</th><th>CC2 [17]</th><th>CAM-B3LYP</th><th>ωB97X-D</th><th>Experiment</th></tr></thead><tbody><tr><td rowspan="3">Fa A</td><td>NH<sub>Phe</sub></td><td>+11/+9</td><td>-6</td><td>+10</td><td>-1</td></tr><tr><td>NH<sub>2 sym.</sub></td><td>-15/-10</td><td>-29</td><td>-9</td><td>-9</td></tr><tr><td>NH<sub>2 anti.</sub></td><td>-18/-9</td><td>-22</td><td>-8</td><td>-6</td></tr><tr><td rowspan="3">Fa B</td><td>NH<sub>Phe</sub></td><td>-44/-37</td><td>-44</td><td>-25</td><td></td></tr><tr><td>NH<sub>2 sym.</sub></td><td>-5/-4</td><td>-5</td><td>-6</td><td></td></tr><tr><td>NH<sub>2 anti.</sub></td><td>-1/-1</td><td>-2</td><td>-3</td><td></td></tr><tr><td rowspan="3">Fa C</td><td>NH<sub>Phe</sub></td><td>-44/-37</td><td>-38</td><td>-44</td><td>-24</td></tr><tr><td>NH<sub>2 sym.</sub></td><td>-3/-1</td><td>-8</td><td>-2</td><td>-1</td></tr><tr><td>NH<sub>2 anti.</sub></td><td>-1/-1</td><td>-1</td><td>+2</td><td>-1</td></tr><tr><td rowspan="3">Fa D</td><td>NH<sub>Phe.</sub></td><td>0/0</td><td>-1</td><td>-1</td><td></td></tr><tr><td>NH<sub>2 sym.</sub></td><td>-2/-2</td><td>0</td><td>-7</td><td></td></tr><tr><td>NH<sub>2 anti</sub></td><td>0</td><td>9</td><td>-1</td><td></td></tr></tbody></table>

Finally, if the conformers stability order compared to both CC2 and experiment (Tables 3 and 4) was well reproduced by the two functionals, we observed a systematic similar strong overestimation of the adiabatic ZPVE-corrected excitations energies. A MAE of 0.40 eV (0.41) with a ME of +0.40 (0.41) eV for ωB97X-D (CAM-B3LYP) compared to CC2 and a MAE of 0.50 eV (0.51) with a ME of +0.50 (+0.51) eV for ωB97X-D (CAM-B3LYP) compared to the experiment were obtained whereas the MAE (ME) for CC2 [18] compared to experiment was of 0.11 (+0.11) eV only. Such a similar behavior was previously reported for a set of 66 organic medium- and large-sized molecules containing aromatic heterocyclic compounds and aromatic (aliphatic) hydrocarbons as well as substituted aromatic hydrocarbons with (MAE, ME) of (0.32, +0.30) eV for ωB97X-D/def2-TZVP ((0.33, +0.30) for CAM-B3LYP)/def2-TZVP) compared to experiment, whereas the CC2/def2-TZVP (MAE, ME) was (0.11, +0.09) eV [33].

<table><caption>**Table 3.** Adiabatic ZPVE-corrected excitation energies of the lowest $\pi\pi^*$ excited state (S<sub>1</sub>) of the Fa A-D conformers obtained at the CC2/aug-(N,O,π)-cc-pVDZ//CC2/cc-pVDZ, the CAM-B3LYP/aug-(N,O,π)-cc-pVDZ//CAM-B3LYP/cc-pVDZ and ωB97X-D/aug-(N,O,π)-cc-pVDZ//ωB97X-D/cc-pVDZ levels, together with the experimental 0–0 transition energies.</caption><thead><tr><th>ΔE<sub>adia</sub> (eV)</th><th>Fa A</th><th>Fa B</th><th>Fa C</th><th>Fa D</th></tr></thead><tbody><tr><td>CC2/aug(N,O,π)-cc-pVDZ//CC2/cc-pVDZ [17]</td><td>4.754</td><td>4.770</td><td>4.767</td><td>4.791</td></tr><tr><td>CAM-B3LYP/aug(N,O,π)-cc-pVDZ//CAM-B3LYP/cc-pVDZ</td><td>5.168</td><td>5.167</td><td>5.142</td><td>5.192</td></tr><tr><td>ωB97XD/aug(N,O,π)-cc-pVTZ//ωB97XD/cc-pVDZ</td><td>5.140</td><td>5.161</td><td>5.160</td><td>5.174</td></tr><tr><td>Experiment [a]</td><td>4.650</td><td>4.663</td><td>4.653</td><td>4.666</td></tr></tbody></table>

<footnote>[a] The experimental values take into account the number of significant digits obtained in the experiment, and the theoretical values are given with the same number of digits.</footnote>
<table><caption>**Table 4.** MAE and ME (eV) on the adiabatic ZPVE-corrected excitation energies of the lowest $\pi\pi^*$ excited state ($S_1$) between DFT/TD-DFT, P1 or P2 and CC2 or Experiment. For DFT/TD-DFT, only the results obtained with ($\omega$B97X-D) are reported. For the definition of P1 and P2, see the following section and Methods section. Bs designed all the systems contained in the benchmark set (see Methods section).</caption><thead><tr><th>Calculation Systems Level</th><th></th><th>MAE CC2</th><th>ME CC2</th><th>MAE Exp.</th><th>ME Exp.</th></tr></thead><tbody><tr><td>DFT/TD-DFT</td><td>Fa A-D</td><td>0.40</td><td>+0.40</td><td>0.50</td><td>+0.50</td></tr><tr><td>DFT/TD-DFT</td><td>Bs</td><td>0.39</td><td>+0.35</td><td>0.50</td><td>+0.43</td></tr><tr><td>P1</td><td>Bs</td><td>0.15</td><td>+0.14</td><td>0.33</td><td>+0.22</td></tr><tr><td>P2</td><td>Bs</td><td>0.03</td><td>+0.02</td><td>0.14</td><td>+0.11</td></tr></tbody></table>

## 2.2. 0–0 Excitation Energies and IR Signature of the $\pi\pi^*$ for the Series of Capped Peptides

The results obtained on the reference system, the Fa A-D conformers, in particular, those highlighting that the dispersion forces are only partially taken into account with CAM-B3LYP, led us to select only the $\omega$B97X-D functional in order to perform the calculations on the series of capped peptides.

First of all, as for the reference system, the nature and the optimized geometry of the lowest $\pi\pi^*$ excited state were well reproduced compared to those obtained at the CC2 level (Tables S8a–c and S9). In particular, FFa A led to two low-lying $\pi\pi^*$ excited states FFa $A_1$ and FFa $A_2$, for which the excitations were localized on either of the phenyl rings whereas FFa C led to one low-lying $\pi\pi^*$ excited state for which the excitation was delocalized on the two phenyl rings (Figure 4 and Figure S9 and Table S9). Moreover, these qualitative differences were in perfect adequacy with the experimental data [18]. Finally, the MAD obtained for the characteristic geometrical parameters of $S_0/S_1$ states were very similar to that obtained for the reference system, 3–4° for the dihedral angles and 0.04–0.05 Å for the intermolecular distances to compare to 3–2° and 0.05–0.05 Å.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Couples of NTOs at the $\omega$B97XD/cc-pVDZ level for FFa $A_1$, $A_2$, and C. Only couples whose contribution to the wave function is greater than 10% are drawn (see Table S9 for the values of contribution).</figcaption></figure>

Regarding the energetics (Tables 4 and 5) and compared to CC2, a systematic strong overestimation of the adiabatic ZPVE-corrected excitation energies (MAE of 0.41 eV with an ME of +0.41 eV) similar to that obtained for the reference system was observed for all the conformers except that of FFa C, which was underestimated (−0.26 eV). The total MAE (reference system plus the series of capped peptides) was equal to 0.39 eV with a total ME of +0.35 eV. A similar strong overestimation was observed when we compared to the experiment: a total MAE of 0.50 eV with a total ME of +0.43 eV. Again, the FFa C adopted a behavior opposite to that of all the other conformers, and its adiabatic ZPVE-corrected excitation energy was strongly underestimated at 0.41 eV. In this context, it is important to remind that this exception concerns the FFa C conformer for which the 0–0 transition was not directly measured but extrapolated by comparison with the vibronic progression
of an analogous system containing only one phenylalanine [18]. Finally, the errors at the ωB97X-D level on the ZPVE-corrected adiabatic excitations energies compared to both CC2 and experiment seems systematic along with the series without any significant size effects.

<table><caption>**Table 5.** Adiabatic ZPVE-corrected excitation energies of the lowest ππ* excited state (S₁) of the GFa A and B', FFa A₁, A₂ and C and QFa A and C conformers obtained at the CC2/aug-(N,O,π)-cc-pVDZ//CC2/cc-pVDZ and ωB97X-D/aug-(N,O,π)-cc-pVDZ//ωB97X-D/cc-pVDZ levels, together with the experimental 0–0 transition energies.</caption><thead><tr><th>ΔE<sub>adia</sub> (eV)</th><th>GFa A</th><th>GFa B'</th><th>FFa A<sub>1</sub></th><th>FFa A<sub>2</sub></th><th>FFa C</th><th>QFa A</th><th>QFa C</th></tr></thead><tbody><tr><td>CC2<br>aug(N,O,π)-cc-pVDZ//<br>cc-pVDZ [17]</td><td>4.754</td><td>4.771</td><td>4.714</td><td>4.729</td><td>4.479</td><td>4.803</td><td>4.793</td></tr><tr><td>ωB97XD<br>aug(N,O,π)-cc-pVTZ//<br>cc-pVDZ</td><td>5.178</td><td>5.178</td><td>5.154</td><td>5.149</td><td>4.219</td><td>5.207</td><td>5.162</td></tr><tr><td>Experiment [a]</td><td>4.648</td><td>4.644</td><td>4.648</td><td>4.658</td><td>4.630</td><td>4.662</td><td>4.658</td></tr></tbody></table>

The harmonic and mode-dependent corrected frequency shifts of the NH stretch vibrations between the ground and the ππ* excited states were calculated for all the conformers and reported in Table 6 together with those at the CC2 level and the experimental ones when available. The mode-dependent corrected shifts have been calculated from the mode-dependent corrected frequencies, which have been determined from the optimal harmonic frequency mode-dependent (NH, NH₂sym, and NH₂anti) linear scaling functions (ν<sub>exp.</sub> = aν<sub>theo.</sub> + b). The linear scaling functions were determined from both the 42 experimental *amide A* region frequencies available for S₀ and the 22 experimental *amide A* region frequencies available for S₁ (Tables S7a and S10b,d,f) and allowed us to take into account in these theoretical S₀ and S₁ frequencies both the method and basis set errors and anharmonicity effects (Figure S11). As followed for the CC2 level, this strategy of systematic correction of calculated harmonic frequencies provides reliable predictions of the theoretical S₁-S₀ frequency shifts in the *amide A* region for such systems, with an RMSD of 6 cm⁻¹ similar to that obtained at the CC2 level, i.e., 5 cm⁻¹ [18].

<footnote>[a] Same details as Table 3.</footnote>
<table><caption>**Table 6.** Theoretical shifts (harmonic x and mode-dependent corrected y, i.e., x or x/y) of the *amide* A region harmonic frequencies of the lowest $\pi\pi^*$ excited state optimized geometry (S₁) relative to the ground state optimized geometry (S₀) for the GFa A and B', FFa A₁, A₂, and C and QFa A, C conformers obtained at the CC2/ and ωB97X-D(x)/cc-pVDZ levels, together with the corresponding available experimental data (cm⁻¹).</caption><thead><tr><th></th><th>Δν<sub>S1/S0</sub> (cm<sup>-1</sup>)</th><th>CC2 [17]</th><th>ωB97X-D</th><th>Experiment</th></tr></thead><tbody><tr><td rowspan="4">GFa A</td><td>NH<sub>Gly</sub></td><td>-7/-6</td><td>-2/-2</td><td>-2</td></tr><tr><td>NH<sub>Phe</sub></td><td>-18/-15</td><td>-10*/-9</td><td>-18</td></tr><tr><td>NH<sub>2</sub> sym.</td><td>-5/-3</td><td>+3*/+2</td><td>+3</td></tr><tr><td>NH<sub>2</sub> anti.</td><td>0/0</td><td>+1/+1</td><td>-9</td></tr><tr><td rowspan="4">GFa B'</td><td>NH<sub>Gly</sub></td><td>0/0</td><td>-1/-1</td><td>+1</td></tr><tr><td>NH<sub>Phe</sub></td><td>-21/-18</td><td>-4/-4</td><td>-18</td></tr><tr><td>NH<sub>2</sub> sym.</td><td>+1/+1</td><td>+1/+1</td><td>+2</td></tr><tr><td>NH<sub>2</sub> anti.</td><td>0/0</td><td>+1/+1</td><td>+1</td></tr><tr><td rowspan="4">FFa A<sub>1</sub></td><td>NH<sub>Phe1</sub></td><td>-41/-35</td><td>-40/-36</td><td>-33</td></tr><tr><td>NH<sub>Phe2</sub></td><td>-9/-8</td><td>-7/-6</td><td>0</td></tr><tr><td>NH<sub>2</sub> sym.</td><td>-4/-3</td><td>+3/+2</td><td>-1</td></tr><tr><td>NH<sub>2</sub> anti.</td><td>-2/-1</td><td>0/0</td><td>0</td></tr><tr><td rowspan="4">FFa A<sub>2</sub></td><td>NH<sub>Phe1</sub></td><td>-5/-4</td><td>-6*/-5</td><td>-1</td></tr><tr><td>NH<sub>Phe2</sub></td><td>-34/-29</td><td>-21*/-19</td><td>-24</td></tr><tr><td>NH<sub>2</sub> sym.</td><td>-2/-1</td><td>+3/+2</td><td>-1</td></tr><tr><td>NH<sub>2</sub> anti.</td><td>-1/-1</td><td>0/0</td><td>0</td></tr><tr><td rowspan="4">FFa C</td><td>NH<sub>Phe1</sub></td><td>-12/-10</td><td>-9/-8</td><td>-1</td></tr><tr><td>NH<sub>Phe2</sub></td><td>-74/-63</td><td>-60*/-55</td><td>-24</td></tr><tr><td>NH<sub>2</sub> sym.</td><td>-30/-21</td><td>-21*/-15</td><td>-1</td></tr><tr><td>NH<sub>2</sub> anti.</td><td>-11/-6</td><td>-7/-6</td><td>0</td></tr><tr><td rowspan="5">QFa A</td><td>NH<sub>Gln</sub></td><td>-2/-2</td><td>-6/-5</td><td>-1</td></tr><tr><td>NH<sub>Phe</sub></td><td>-12/-10</td><td>-22/-20</td><td>-1</td></tr><tr><td>NH<sub>2</sub> sym./C-term</td><td>+5/+3</td><td>-2/-1</td><td>0/0</td></tr><tr><td>NH<sub>2</sub> anti./C-term</td><td>+8/+4</td><td>0/0</td><td>0/0</td></tr><tr><td>NH<sub>2</sub> sym./cChain</td><td>-2/-1</td><td>0/0</td><td>0/0</td></tr><tr><td rowspan="5">QFa C</td><td>NH<sub>Gln</sub></td><td>-9/-8</td><td>-9/-8</td><td>-1</td></tr><tr><td>NH<sub>Phe</sub></td><td>-16/-14</td><td>-32/-29</td><td>-1</td></tr><tr><td>NH<sub>2</sub> sym./C-term</td><td>-2/-1</td><td>-2/-1</td><td>0/0</td></tr><tr><td>NH<sub>2</sub> anti./C-term</td><td>0/0</td><td>0/0</td><td>0/0</td></tr><tr><td>NH<sub>2</sub> sym./cChain</td><td>-1/-1</td><td>-2/-1</td><td>0/0</td></tr><tr><td></td><td>NH<sub>2</sub> anti./Chain</td><td>+1/+1</td><td>+3/+2</td><td>0/0</td></tr></tbody></table>

## 2.3. Composite Protocols CC2-DFT/TDDFT and 0-0 Excitation Energies Calculations

The adiabatic ZPVE-corrected excitation energies were calculated according to two protocols for the full series of capped peptides (Tables 4 and 7). These two protocols combined DFT/TD-DFT and CC2 calculations. In the first protocol denoted hereafter P1, both the geometry optimizations and harmonic frequencies calculations were performed at the DFT/TD-DFT level. In the second protocol, P2, the geometry optimizations were performed at the CC2 level, whereas the harmonic frequencies calculations were performed at the DFT/TD-DFT level (see Methods section for details). A systematic overestimation was observed except for FFa C, like what was obtained at the CC2 level. The two protocols led to a decrease in the error compared to both CC2 and the experiment. However, this decrease was moderate for protocol P1, whereas protocol P2 exhibited a strong decrease, eventually providing an error compared to the experiment similar to that obtained at the CC2 level. Indeed, P1 gives an MAE of 0.15 eV with an ME of +0.14 eV compared to CC2 and an MAE of 0.33 eV with an ME of +0.22 eV compared to the experiment, whereas P2 gives an MAE of 0.03 eV with an ME of +0.02 eV compared to CC2 and an MAE of 0.14 eV

<footnote>* Coupled modes</footnote>
with an ME of +0.11 eV compared to the experiment. In addition, as for the CC2 method, no size effects were observed for the two protocols. Indeed, in the P1 protocol, the MAE for the capped peptides containing one residue was equal to 0.20 eV (ME of +0.20 eV) compared to CC2 and 0.31 eV (ME of +0.31 eV) compared to the experiment, whereas the MAE for the capped peptides containing two residues was equal to 0.23 eV (ME of +0.10 eV) compared to CC2 and 0.34 eV (ME of +0.17 eV) compared to experiment. The protocol P2, compared to CC2, provides an MAE of 0.02 eV (ME of +0.02 eV) for the capped peptides containing one residue and an MAE of 0.03 eV (ME of +0.02 eV) for the capped peptides containing two residues.

The difference between MAE and ME for the systems containing two residues is essentially due to the underestimation obtained for FFa C. In addition, a similar behavior was observed for the comparison with the experiment: the MAE was equal to 0.14 eV (ME of +0.14 eV) for the capped peptides containing one residue and equal to 0.14 eV (ME of +0.09) for the capped peptides containing two residues. The protocol P2 allows then to reproduce very well the CC2 adiabatic ZPVE-corrected excitation energies and then to obtain a quantitative agreement on these data compared to the experiment as for the CC2 method. This very good agreement can be explained by the low sensitivity of the ΔΔE₀₀ term, the difference of ZPVE corrections between the ground (S₀) and excited (S₁) states relative to the method, CC2 or DFT/TD-DFT.

<table><caption>**Table 7.** Adiabatic ZPVE-corrected excitation energies of the lowest $\pi\pi^*$ excited state (S₁) of the Fa A-D, GFa A and B', FFa A₁, A₂, and C and QFa A and C conformers obtained with the protocols P1 and P2, together with the CC2/aug-(N,O,π)-cc-pVDZ//CC2/cc-pVDZ ones and the experimental 0–0 transition energies.</caption><thead><tr><th>ΔE<sub>adia</sub> (eV)</th><th>CC2 [17]</th><th>DFT/TD-DFT<br>ωB97XD</th><th>P1</th><th>P2</th><th>Experiment <sup>[a]</sup></th></tr></thead><tbody><tr><td>Fa A</td><td>4.754</td><td>5.140</td><td>4.953</td><td>4.778</td><td>4.650</td></tr><tr><td>Fa B</td><td>4.770</td><td>5.161</td><td>4.972</td><td>4.798</td><td>4.663</td></tr><tr><td>Fa C</td><td>4.767</td><td>5.160</td><td>4.971</td><td>4.793</td><td>4.663</td></tr><tr><td>Fa D</td><td>4.791</td><td>5.174</td><td>4.983</td><td>4.809</td><td>4.666</td></tr><tr><td>GFa A</td><td>4.754</td><td>5.178</td><td>4.867</td><td>4.814</td><td>4.648</td></tr><tr><td>GFa B'</td><td>4.771</td><td>5.178</td><td>4.979</td><td>4.802</td><td>4.644</td></tr><tr><td>FFa A<sub>1</sub></td><td>4.714</td><td>5.154</td><td>4.962</td><td>4.771</td><td>4.648</td></tr><tr><td>FFa A<sub>2</sub></td><td>4.729</td><td>5.149</td><td>4.954</td><td>4.773</td><td>4.658</td></tr><tr><td>FFa C</td><td>4.479</td><td>4.219</td><td>4.023</td><td>4.451</td><td>4.630</td></tr><tr><td>QFa A</td><td>4.803</td><td>5.207</td><td>5.014</td><td>4.825</td><td>4.662</td></tr><tr><td>QFa C</td><td>4.793</td><td>5.162</td><td>4.965</td><td>4.779</td><td>4.658</td></tr></tbody></table>

Such behavior was already reported by Durbeej et al. for the set of 66 organic medium-sized and large molecules containing aromatic heterocyclic compounds and aromatic (aliphatic) hydrocarbons as well as substituted aromatic hydrocarbons [34]. Finally, these two protocols give results that, compared to the experiment, behave identically to the CC2 ones, i.e., systematic overestimation and no size effect. They can, therefore, be a reliable alternative to the CC2 method for large and very large-sized systems with two levels of accuracy, a high level with P2 and a moderate one with P1.

### 3. Methods

#### 3.1. Benchmark Set Composition

The benchmark set (Figure 1) consists of the four conformers (A–D) of a reference system, N-acetyl-phenylalaninylamide (Ac-Phe-NH₂, Phe: phenylalanine), and two conformers of three capped dipeptides, Ac-Gly-Phe-NH₂ (A and B', Gly: glycine), Ac-Phe-Phe-NH₂ (A and C) and Ac-Gln-Phe-NH₂ (A and C, Gln: glutamine). These systems adopt in their ground state prototypical secondary structural features of proteins [15] and present specific low-lying excited states, such as localized or delocalized states on the different significant

<footnote>[a] Same details as Table 3.</footnote>
groups of these systems (aromatic rings and peptide bonds) as well as charge transfer states from the backbone to the phenyl ring. Detailed information regarding these systems is summarized below. The L label refers to the preferential backbone orientation adopted by the natural amino acids of L configuration, as opposed to those of D configuration. The labels a, g+, or g- refer to the anti, gauche+, or gauche- orientation of the phenyl side chain relative to the backbone, the $\chi_i^1$ dihedral angle $N-C_\alpha-C_\beta-C_\gamma$ being close to +180° for a, +60° for g+, and -60° for g- (Figure S1b). Finally, the index n of the H-bond label $C_n$ indicates the number of atoms in the ring formed by the H-bond.

**The reference system, Fa:**[16,18,36,38] In the ground state, four conformers (A–D) with two different types of folding backbone were observed and assigned to conformations lying in an energy range of 0–6.5 kJ/mol: one β-strand extended conformation (A ($\beta_L$(a))) stemming from a $C_5$ H-bond and three γ-turn folded conformations stabilized by a $C_7$ H-bond and differing by the orientation of the phenyl side chain (B ($\gamma_L$(g+))), $C(\gamma_L$($g$-)) and $D(\gamma_L$($a$))). Three structures present an NH...π bond, the β-strand conformation A, and two γ-turns, B and C. The corresponding 0–0 excitation energy to the lowest $\pi\pi^*$ excited state was both theoretically and experimentally determined for each conformer, whereas only excited state IR (ESIR) spectra of A and C were measured.

**GFa:**[15,16,18,38] Among the five conformers observed and assigned to conformations lying in an energy range of 0–6 kJ/mol, two conformers were chosen, A and B'. They correspond to two different types of folding backbone: a 7-7 (double γ-turn) extended conformation with an L chirality of the turns (A ($7_L-7_L$($g$-))) and a β-turn folded conformation of types II' ($B'$(π-10II'($g$+))). These two families of backbone folding correspond to secondary structural features of proteins: the $\gamma_2$ ribbon and β-turn secondary structures, respectively stabilized by successive $C_7$ H-bonds (γ-turns) and $C_{10}$ H-bonds (β-turn). In addition, one structure, the β-turn, presents an NH...π bond. The corresponding 0-0 excitation energy to the lowest $\pi\pi^*$ excited state was both theoretically and experimentally determined for each conformer as well as their ESIR spectra.

**FFa:**[16,18,39] Among the three conformers observed and assigned to conformations lying in an energy range of 0–15 kJ/mol, two conformers, A and C, were chosen. They correspond to both different types of backbone folding and different relative orientations of the two phenyl rings: one β-turn type I conformation (an $\alpha_L-\gamma_L$ structure, A (π-π-10I($g$+,g+))) where the phenyl rings interact according to a T-shape arrangement and one $\beta_L-\gamma_L$ conformation (C ($5-\pi-7_L$($a$,$g$+))) where the phenyl rings interact according to a face-to-face arrangement. FFa A is characterized by a $C_{10}$ H-bond and C by a usual combination of $C_5$ and $C_7^{eq}$ H-bonds. In addition, these conformations present at least one NH-π bond. For the A conformer, the two intense transitions observed in the near UV spectrum were assigned, in accordance with the CC2 calculations and ESIR spectra, to the origin transition of each excited chromophore, $A_1$ and $A_2$. On the contrary, the experiment suggested that the C conformer exhibits only one intense transition in the near UV spectrum, which can be assigned to the lowest $\pi\pi^*$ excited state. CC2 calculations confirmed this hypothesis showing that this excited state was delocalized on the two chromophores. Only ESIR spectra of $A_1$ and $A_2$ were measured.

**QFa:**[18,40] Among the three conformers observed and assigned to conformations lying in an energy range of 0–7.5 kJ/mol, two conformers were chosen, A and C. They corresponded to the same type of folding backbone, i.e., a type I β-turn backbone, stabilized by a $C_{10}$ H-bond. This bond was combined to a side chain/main chain $C_7$ H-bond bridging the NH site of the first peptide bond to the oxygen atom of the Gln residue side chain CO-NH₂ group labeled $7^\epsilon$ and presenting a similar NH...π bond; A ($7^\epsilon$-π-10I($g$+))) and C ($7^\epsilon$-π-10I($g$+))). The two structures differed by the orientation of the CO-NH₂ group, i.e., the conformation of the Gln residue side chain. The corresponding 0–0 excitation energy to the lowest $\pi\pi^*$ excited state was determined both theoretically and experimentally for both conformers.
## 3.2. DFT/TD-DFT Calculations

DFT/TD-DFT calculations were carried out with the Gaussian Package [41]. Ground-state DFT and excited-state TD-DFT calculations were performed by using three long-range corrected hybrid functionals, LC-ωPBE [28], CAM-B3LYP [29], and ωB97X-D [30]. The difference between these long-range hybrid GGA functionals resides mainly in the percentage of short and long-range exact exchange and the explicit taken into account of a dispersion correction. The CAM-B3LYP functional comprises 19% exact exchange plus 81% Becke 1988 (B88) exchange interaction at short-range, and 65% exact exchange plus 35% B88 at long-range, the intermediate region being smoothly described through the standard error function with parameter 0.33. In the LC-ωPBE functional, the short-range and long-range exact exchange coefficient is equal to 1 with an ω equal to 0.40 bohr⁻¹. The ωB97X-D functional includes 22.2% short-range exact exchange and 100% long-range exact exchange controlled by an ω of 0.20 bohr⁻¹. In addition, this functional includes explicitly a dispersion correction via empirical potentials of the form −C₆/R⁶. In order to compare to our previous CC2 calculations,[18] we used the same combination of basis sets: (i) the Dunning's correlation-consistent polarized valence double-ξ (cc-pVDZ) basis set [42] for single-point energy calculations at the ground state of the reference system, the geometry optimizations of both ground state and the first ππ* excited states of the benchmark set as well as the harmonic frequencies calculations; (ii) this basis set was augmented with diffuse functions taken from the aug-cc-pVDZ [43] on each oxygen and nitrogen atom and only for one in two carbon atoms of the phenyl ring for the single-point energy calculations. Our previous work demonstrated on a series of capped peptides of increasing size and containing different residues that the CC2 method in combination with the aug(N,O,π)-cc-pVDZ//cc-pVDZ basis sets, proves to be very reliable for calculations of the 0-0 excitation energies of such systems compared to experiment, leading to a mean absolute error (MAE) of 0.10 eV and a mean signed error (ME) of +0.08 eV. In addition, this level of theory was further justified by the faster convergence of DFT/TD-DFT calculations, with respect to the size of the one-electron basis set, compared to conventional correlated methods such as CC2. All DFT and TD-DFT calculations were performed using the *ultrafine* grid. In the ground-state DFT calculations, the convergence criterion on the RMS density was 10⁻⁸. In the TD-FDT calculations, fifteen states were solved, and the convergence threshold was 10⁻³ au on the wave function and 10⁻⁶ au on the energy for single-point energy calculations and 10⁻⁶ au on the wave function and 10⁻⁸ au on the energy for geometry optimizations. The geometry of both ground and excited states was optimized until the residual mean force was smaller than 4.5 10⁻⁴ au, and the harmonic frequencies were analytically (DFT) or numerically (TD-DFT) determined. Harmonic frequencies calculations allow us to both verify that the optimized geometries correspond to true minima and calculate for each conformer and state the zero-point vibrational energy (ZPVE). In addition, taking advantage of a large amount of experimental data available (42 *amide A* region frequencies available for S₀ and 22 for S₁), we determined optimal harmonic frequency mode-dependent (NH, NH₂sym, and NH₂anti) linear scaling functions (νₚₐₗ. = aνₚₐₗ. + b) in order to correct the calculated S₀ and S₁ harmonic frequencies for method and basis set errors as well as anharmonicity effects.

## 3.3. *Composite CC2/DFT-TDDFT Protocols for 0-0 Excitation Energies*

We evaluated two composite protocols to calculate the adiabatic ZPVE-corrected excitation energies (ΔE₀₀) according to the approach developed in ref. [32]. In the first one, denoted hereafter P1, the CC2 [21–25] method was only used to calculate the vertical excitation energy calculation (ΔE_v), and both the geometry optimizations and harmonic frequencies calculations were performed with the DFT(S₀) and TD-DFT(S₁) methods. In the second protocol (P2), the CC2 method was used for both the vertical excitation energy
calculation and the geometry optimizations, but the harmonic frequencies calculations were performed with the DFT(S₀) and TD-DFT(S₁) methods:

$$ \Delta E_{00} = \Delta E_{v} + \Delta\Delta E_{ad} + \Delta\Delta E_{00} $$

$$ \Delta\Delta E_{ad} = \Delta E_{ad} - \Delta E_{v} $$

$$ \Delta\Delta E_{00} = \Delta E_{00} - \Delta E_{ad} $$

$$ P1: \Delta E_{00} = \Delta E_{v}(CC2) + \Delta\Delta E_{ad}(DFT/TD - DFT) + \Delta\Delta E_{00}(DFT/TD - DFT) $$

$$ P2: \Delta E_{00} = \Delta E_{v}(CC2) + \Delta\Delta E_{ad}(CC2) + \Delta\Delta E_{00}(DFT/TD - DFT) $$

where ΔE_v is the vertical excitation energy at the equilibrium geometry of the ground-state, ΔE_ad is the energy difference between the ground and excited states at their respective optimized geometry, and ΔE_00 is this energy difference, including the zero-point vibrational energy (ZPVE) corrections. A combination of two basis set was used for all the calculations. First, the geometry optimizations and frequency calculations were performed within the cc-pVDZ basis set, whatever the method used, and second, the single-point energy calculations for the determination of the vertical excitation energies were performed within the aug(N,O,π)-cc-pVDZ.

The CC2 computational details summarized hereafter were those of our previous benchmark [18]. CC2 calculations were performed with the TURBOMOLE package [44] using the resolution-of-identity (RI) approximation and associated auxiliary basis sets [45]. Frozen cores for the 1s electrons were employed, and all calculations were carried out in the C₁ point-group symmetry. Ten singlet states were considered, and D₁, D₂ diagnostics and % < $\bar{E}_1$ | E₁ > borthogonal norm were calculated in order to evaluate the capability of the CC2 method to properly describe the ground and excited states of such systems [25,46,47]. The convergence criterion used in single-point energy calculation was 10⁻⁸ au on the density for the HF calculation, 10⁻⁹ au for the RI-CC2 ground state energy for the iterative coupled-cluster methods, and 10⁻⁶ au for the convergence threshold for the norm of residual vectors in eigenvalue problems for the RI-CC2 excited states calculations. In the geometry optimization of both ground and lowest ππ* excited states, the convergence criterion used corresponded to a norm of the Cartesian gradient lower than 10⁻⁴ au. Finally, the harmonic frequencies were calculated by numerical differentiation of the analytic gradients using central differences and a step length of 0.02 au.

### 3.4. Characterization of the Excited States

Orbital-relaxed first-order properties were determined, in particular the density, at DFT or CC2 level for the ground states and at TD-DFT or CC2 level for the excited states. DFT/TD-DFT and CC2 differences between the density of excited states and that of the ground state were then calculated. In addition, a post-processing tool interfaced to TURBOMOLE and GAUSSIAN, Nancy_EX-2.0 [48], was used in order to analyze the density and character of the excited states and obtain, at both TD-DFT and CC2 level, the so-called natural transition orbitals (NTOs) [49,50] of each state. Instead of describing one excitation with multiple canonical spin orbitals couples, all the physical information on the nature of the electronic transition was gathered in one (sometimes two) couple(s) of NTOs.

### 4. Conclusions

Benchmarking calculations were performed to determine the ability of the TD-DFT method to model qualitatively the PES of bio-relevant systems such as capped peptides. In the first step, the most accurate functional, ωB97X-D, among three long-range corrected hybrid functionals have been selected on criteria based on energetics as well as the first and second derivatives of the energy of the low-lying excited states of the four conformers of the smallest capped peptide, Fa A-D. This highlights that the intramolecular dispersion interactions can play an important role in the geometry of prototypical secondary structures of proteins and must be taken into account explicitly. In the second step, this functional was validated on a series of capped peptides of increasing size and containing residues
of different nature. We obtained for the geometry optimization of both ground and $\pi\pi^*$ excited states as well as for the harmonic frequencies calculations a quantitative agreement compared to CC2 and experiment. However, the adiabatic ZPVE-corrected excitation energies were systematically strongly overestimated, and we then developed protocols combining CC2 and DFT/TD-DFT methods, and one of these protocols allowed us to obtain excitation energies of CC2 quality at the lowest computational cost (Table S12).

In conclusion, the TD-DFT method with a system-appropriate functional, the ωB97X-D for capped peptides, is well suited for a qualitative exploration of the PES of the low-lying excited states, a key step in non-adiabatic dynamics simulations performed to determine deactivation mechanisms in bio-relevant systems. In addition, computation of adiabatic ZPVE-corrected excitation energies of CC2 quality can be obtained by employing CC2 for both the single-point energy calculations (vertical excitation energies) and the geometry optimizations and DFT/TD-DFT for the calculation of the ZPVE corrections. We even obtained a moderate agreement using the CC2 method only for the single-point energy calculations and the DFT/TD-DFT for both the geometry optimizations and frequencies calculations. These protocols, therefore, constitute an alternative to the CC2 method for very large systems, i.e., when this computationally demanding method (see Table S12) is no longer applicable. Finally, the TD-DFT method exhibits a favorable cost-performance ratio, and its involvement in the calculations can be adapted according to the level of expected accuracy.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/ijms23020621/s1.

**Author Contributions:** Theoretical calculations, M.L., T.V., and V.B.; Experimental data, B.T., E.G., and M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work received financial support from the Agence Nationale de la Recherche (ANR) as Grant ANR-14-CE06-0019-01- ESBODYR. This work was granted access to the HPC facility of [TGCC/CINES/IDRIS] under the Grant 2015-t2015087412 and t2015087372, 2016-t2016087540, 2017-A0010807540 and 2020-A0070807540 awarded by GENCI (Grand Equipement National de Calcul Intensif) and to the CCRT High-Performance Computing (HPC) facility at CEA under the Grant CCRT2015/CCRT2016/CCRT2017/CCRT2020-p606bren.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study and not reported in the Supplementary Materials are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## References

1. El-Khoury, P.Z.; Scharpio, I.; Huntress, M.; Melaccio, F.; Gozen, S.; Frutos, L.M.; Olivucci, M. Computational Photochemistry and Photobiology. In *CRC Handbook of Organic Photochemistry and Photobiology*, 3rd ed.; CRC Press: Boca Raton, FL, USA, 2012; Volume 2, pp. 1029–1056.

2. Worth, G.A.; Cederbaum, L.S. Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection. *Annu. Rev. Phys. Chem.* **2004**, *55*, 127–158. [CrossRef] [PubMed]

3. Levine, B.G.; Martinez, T.J. Isomerization through conical intersections. *Annu. Rev. Phys. Chem.* **2007**, *58*, 613–634. [CrossRef] [PubMed]

4. Domcke, W.; Yarkony, D.R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. *Annu. Rev. Phys. Chem.* **2012**, *63*, 325–352. [CrossRef] [PubMed]

5. Yarkony, D.R. Nonadiabatic Quantum Chemistry-Past, Present, and Future. *Chem. Rev.* **2012**, *112*, 481–498. [CrossRef]

6. Runge, E.; Gross, E.K.U. Density-functional theory for time-dependent systems. *Phys. Rev. Lett.* **1984**, *52*, 997–1000. [CrossRef]

7. Marques, M.A.L.; Gross, E.K.U. Time-dependent density functional theory. *Annu. Rev. Phys. Chem.* **2004**, *55*, 427–455. [CrossRef]

8. Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. *Chem. Rev.* **2005**, *105*, 4009–4037. [CrossRef]
9. Casida, M.E.; Huix-Rotllant, M. Progress in Time-Dependent Density-Functional Theory. *Annu. Rev. Phys. Chem.* **2012**, *63*, 287–323. [CrossRef]

10. Maitra, N.T. Perspective: Fundamental aspects of time-dependent density functional theory. *J. Chem. Phys.* **2016**, *144*, 220901. [CrossRef]

11. Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. *J. Chem. Phys.* **2004**, *120*, 8425–8433. [CrossRef]

12. Malis, M.; Loquais, Y.; Gloaguen, E.; Biswal, H.S.; Piuzzi, F.; Tardivel, B.; Brenner, V.; Broquier, M.; Jouvet, C.; Mons, M.; et al. Unraveling the Mechanisms of Nonradiative Deactivation in Model Peptides Following Photoexcitation of a Phenylalanine Residue. *J. Am. Chem. Soc.* **2012**, *134*, 20340–20351. [CrossRef]

13. Malis, M.; Loquais, Y.; Gloaguen, E.; Jouvet, C.; Brenner, V.; Mons, M.; Ljubic, I.; Doslic, N. Non-radiative relaxation of UV photoexcited phenylalanine residues: Probing the role of conical intersections by chemical substitution. *Phys. Chem. Chem. Phys.* **2014**, *16*, 2285–2288. [CrossRef]

14. Loquais, Y.; Gloaguen, E.; Alauddin, M.; Brenner, V.; Tardivel, B.; Mons, M. On the near UV photophysics of a phenylalanine residue: Conformation-dependent pi pi\* state deactivation revealed by laser spectroscopy of isolated neutral dipeptides. *Phys. Chem. Chem. Phys.* **2014**, *16*, 22192–22200. [CrossRef]

15. Loquais, Y.; Gloaguen, E.; Habka, S.; Vaquero-Vara, V.; Brenner, V.; Tardivel, B.; Mons, M. Secondary Structures in Phe-Containing Isolated Dipeptide Chains: Laser Spectroscopy vs Quantum Chemistry. *J. Phys. Chem. A* **2015**, *119*, 5932–5941. [CrossRef]

16. Sohn, W.Y.; Brenner, V.; Gloaguen, E.; Mons, M. Local NH-pi interactions involving aromatic residues of proteins: Influence of backbone conformation and pi pi\* excitation on the pi H-bond strength, as revealed from studies of isolated model peptides. *Phys. Chem. Chem. Phys.* **2016**, *18*, 29969–29978. [CrossRef]

17. Ben Amor, N.; Hoyau, S.; Maynau, D.; Brenner, V. Low-Lying excited states of model proteins: Performances of the CC2 method versus multireference methods. *J. Chem. Phys.* **2018**, *148*, 184105. [CrossRef]

18. Dupuy, M.-S.; Gloaguen, E.; Tardivel, B.; Mons, M.; Brenner, V. CC2 Benchmark for Models of Phenylalanine Protein Chains: 0-0 Transition Energies and IR Signatures of the ππ\* Excited State. *J. Chem. Theory Comput.* **2020**, *16*, 601–611. [CrossRef]

19. Brenner, V.; Very, T.; Schmidt, M.; Gordon, M.; Hoyau, S.; Ben Amor, N. Model protein excited states: MRCI calculations with large active spaces vs CC2 method. *J. Chem. Phys.* **2021**, *154*, 214105. [CrossRef]

20. Serrano-Andres, L.; Fulscher, M.P. Theoretical study of the electronic spectroscopy of peptides. III. Charge-transfer transitions in polypeptides. *J. Am. Chem. Soc.* **1998**, *120*, 10912–10920. [CrossRef]

21. Christiansen, O.; Koch, H.; Jorgensen, P. The 2nd-order approximate coupled-cluster singles and doubles model CC2. *Chem. Phys. Lett.* **1995**, *243*, 409–418. [CrossRef]

22. Hattig, C.; Weigend, F. CC2 excitation energy calculations on large molecules using the resolution of the identity approximation. *J. Chem. Phys.* **2000**, *113*, 5154–5161. [CrossRef]

23. Hattig, C.; Kohn, A. Transition moments and excited-state first-order properties in the coupled-cluster model CC2 using the resolution-of-the-identity approximation. *J. Chem. Phys.* **2002**, *117*, 6939–6951. [CrossRef]

24. Hattig, C. Geometry optimizations with the coupled-cluster model CC2 using the resolution-of-the-identity approximation. *J. Chem. Phys.* **2003**, *118*, 7751–7761. [CrossRef]

25. Kohn, A.; Hattig, C. Analytic gradients for excited states in the coupled-cluster model CC2 employing the resolution-of-the-identity approximation. *J. Chem. Phys.* **2003**, *119*, 5021–5036. [CrossRef]

26. Bories, B.; Maynau, D.; Bonnet, M.L. Selected excitation for CAS-SDCI calculations. *J. Comput. Chem.* **2007**, *28*, 632–643. [CrossRef]

27. Ben Amor, N.; Bessac, F.; Hoyau, S.; Maynau, D. Direct selected multireference configuration interaction calculations for large systems using localized orbitals. *J. Chem. Phys.* **2011**, *135*, 014101. [CrossRef]

28. Vydrov, O.A.; Scuseria, G.E.; Perdew, J.P. Tests of functionals for systems with fractional electron number. *J. Chem. Phys.* **2007**, *126*, 154109. [CrossRef]

29. Yanai, T.; Tew, D.P.; Handy, N.C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). *Chem. Phys. Lett.* **2004**, *393*, 51–57. [CrossRef]

30. Chai, J.D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. *Phys. Chem. Chem. Phys.* **2008**, *10*, 6615–6620. [CrossRef] [PubMed]

31. Maitra, N.T. Charge transfer in time-dependent density functional theory. *J. Phys.-Condens. Matter* **2017**, *29*, 423001. [CrossRef] [PubMed]

32. Winter, N.O.C.; Graf, N.K.; Leutwyler, S.; Hattig, C. Benchmarks for 0-0 transitions of aromatic organic molecules: DFT/B3LYP, ADC(2), CC2, SOS-CC2 and SCS-CC2 compared to high-resolution gas-phase data. *Phys. Chem. Chem. Phys.* **2013**, *15*, 6623–6630. [CrossRef] [PubMed]

33. Fang, C.F.; Oruganti, B.; Durbeej, B. How Method-Dependent Are Calculated Differences between Vertical, Adiabatic, and 0-0 Excitation Energies? *J. Phys. Chem. A* **2014**, *118*, 4157–4171. [CrossRef] [PubMed]

34. Oruganti, B.; Fang, C.F.; Durbeej, B. Assessment of a composite CC2/DFT procedure for calculating 0-0 excitation energies of organic molecules. *Mol. Phys.* **2016**, *114*, 3448–3463. [CrossRef]

35. Loos, P.F.; Jacquemin, D. Evaluating 0-0 Energies with Theoretical Tools: A Short Review. *Chemphotochem* **2019**, *3*, 684–696. [CrossRef]
36. Chin, W.; Mons, M.; Dognon, J.P.; Piuzzi, F.; Tardivel, B.; Dimicoli, I. Competition between local conformational preferences and secondary structures in gas-phase model tripeptides as revealed by laser spectroscopy and theoretical chemistry. *Phys. Chem. Chem. Phys.* **2004**, *6*, 2700–2709. [CrossRef]

37. Wang, J.; Durbeej, B. How accurate are TD-DFT excited-state geometries compared to DFT ground-state geometries? *J. Comput. Chem.* **2020**, *41*, 1718–1729. [CrossRef]

38. Alaudin, M.; Vaquero-Vara, V.; Habka, S.; Tardivel, B.; Gloaguen, E.; Mons, M. Private communication, 2020.

39. Gloaguen, E.; Loquais, Y.; Thomas, J.A.; Pratt, D.W.; Mons, M. Spontaneous Formation of Hydrophobic Domains in Isolated Peptides. *J. Phys. Chem. B* **2013**, *117*, 4945–4955. [CrossRef]

40. Vaquero-Vara, V.; Sohn, W.Y.; Tardivel, B.; Brenner, V.; Gloaguen, E.; Mons, M. Private communication, 2020.

41. Gaussian, D.R.; Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; et al. Gaussian 16, Revision C.01, Gaussian, Inc.: Wallingford, CT, USA, 2016.

42. Dunning, T.H. Gaussian-basis sets for use in correlated molecular calculations.1. The atoms boron through neon and hydrogen. *J. Chem. Phys.* **1989**, *90*, 1007–1023. [CrossRef]

43. Kendall, R.A.; Dunning, T.H.; Harrison, R.J. Electron-affinities of the 1st-row atoms revisited-systematic basis-sets and wave-functions. *J. Chem. Phys.* **1992**, *96*, 6796–6806. [CrossRef]

44. TURBOMOLE v6.4-2012, v7.0-2015 & v7.2-2017, A Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007. Available online: http://www.turbomole.com (accessed on 12 November 2021).

45. Weigend, F.; Kohn, A.; Hattig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. *J. Chem. Phys.* **2002**, *116*, 3175–3183. [CrossRef]

46. Janssen, C.L.; Nielsen, I.M.B. New diagnostics for coupled-cluster and Moller-Plesset perturbation theory. *Chem. Phys. Lett.* **1998**, *290*, 423–430. [CrossRef]

47. Nielsen, I.M.B.; Janssen, C.L. Double-substitution-based diagnostics for coupled-cluster and Moller-Plesset perturbation theory. *Chem. Phys. Lett.* **1999**, *310*, 568–576. [CrossRef]

48. NAncy-EX. Available online: http://nancyex.sourceforge.net (accessed on 13 June 2014).

49. Etienne, T.; Assfeld, X.; Monari, A. Toward a Quantitative Assessment of Electronic Transitions' Charge-Transfer Character. *J. Chem. Theory Comput.* **2014**, *10*, 3896–3905. [CrossRef]

50. Etienne, T.; Assfeld, X.; Monari, A. New Insight into the Topology of Excited States through Detachment/Attachment Density Matrices-Based Centroids of Charge. *J. Chem. Theory Comput.* **2014**, *10*, 3906–3914. [CrossRef]