Ivana Antol

Zoran Glasovac

Davor Margetić

Rachel Crespo-Otero

Mario Barbatti

Insights on the Auxochromic Properties of Guanidinium Group

UV/Vis spectra of phenylguanidine (PHGU) in the gas phase and in acetonitrile have been simulated by TD-DFT calculations. Several DFT hybrid and long-range corrected functionals were tested with respect to CASPT2 gas phase calculations. Solvent effects were considered using polarizable continuum model (PCM) and compared with the measured data in acetonitrile.

Comparison with isoelectronic phenyl urea and related phenyltiourea was done as well. The PBE0 and long-range corrected CAM-B3LYP functionals were selected to investigate the effect of protonation on the excitation energies and absorption intensities of PHGU and several guanidine derivatives with different aromatic chromophoric groups (naphthyl, anthracenyl, quinolinyl, anthraquinonyl and coumarinyl). Also, the effect of complexation and specific interactions through hydrogen bonds with different anions was examined. It was shown that the protonation of guanidine subunit shifts the low energy absorption bands toward higher energies (hypsochromic shift). The shift is reduced upon complexation with anions. In phenylguanidine salts, max are correlated to the anion basicity and strength of H-bonding. The observed changes diminish upon increase of chromophoric size (naphthyl, anthracenyl). Theoretical predictions of UV/Vis spectra correlate well with experimentally measured spectra of selected guanidine derivatives and their salts.

INTRODUCTION

Hydrogen bonding is a ubiquitous phenomenon in nature, which is of great importance in functioning of living organisms. A vast area of supramolecular chemistry could be recognized as the "chemistry of hydrogen bonded systems". Recently, an overview of new discoveries in this research field has been published encompassing usage of hydrogen bonds from anion binding to construction of large polymeric structures. [START_REF]Hydrogen Bonded Supramolecular Structures[END_REF] Due to the mobility of the proton involved in hydrogen bond, this interaction can be used for regulation of electron transfer processes. For instance, guanidine-carboxylate salt bridge formed upon arginine-aspartate hydrogen bonding interaction regulates electron transfer in some biologically important redox processes. [START_REF] Puglisi | Role of RNA Structure in Arginine Recognition of TAR RNA[END_REF][START_REF] Pavletich | Zinc Finger-DNA Recognition: Crystal Structure of a Zif268-DNA Complex at 2.1 A[END_REF][START_REF] Berg | Zinc Finger Domains: From Predictions to Design[END_REF][START_REF] Howell | Functional Role of Aspartic Acid-27 in Dihydrofolate Reductase Revealed by Mutagenesis[END_REF][START_REF] Crane | Sulfite reductase structure at 1.6 A: evolution and catalysis for reduction of inorganic anions[END_REF] Nocera and co-workers also showed that electron transfer rates are significantly influenced by the directionality of the salt bridge dipole. [START_REF] Kirby | Significant Effect of Salt Bridges on Electron Transfer[END_REF] In such systems, electron transfer across the salt bridge is usually coupled with proton movement.

To achieve an efficient binding of anions, various proton donating functional groups have been used, among which ureas, [START_REF] Amendola | Anion Recognition by Hydrogen Bonding: Urea-Based Receptors[END_REF] thioureas, [START_REF] Li | Anion Complexation and Sensing Using Modified Urea and Thiourea-Based Receptors[END_REF] amides, [START_REF] Zhang | Selective Anion Sensing by a Tris-Amide CTV Derivative: 1H NMR Titration, Self-Assembled Monolayers, and Impedance Spectroscopy[END_REF][START_REF] Chmielewski | Anion Recognition by Neutral Macrocyclic Amides[END_REF][START_REF] Gupta | Anion Recognition Through Amide-Based Dendritic Molecule: A Poly(Vinyl Chloride) Based Sensor for Nitrate Ion[END_REF] amidines, [START_REF] Nieto | Pt(II)-Activated Coupling of Aminoethylferrocene with Benzonitrile. A Facile Access Route to a New Redox-Active Bis(ferrocenyl-amidine) Anion Sensor[END_REF] and guanidines [START_REF] Schmuck | Guanidinium Based Anion Receptors in[END_REF][START_REF] Schug | Noncovalent Binding between Guanidinium and Anionic Groups: Focus on Biological-and Synthetic-Based Arginine/Guanidinium Interactions with Phosph[on]ate and Sulf[on]ate Residues[END_REF][START_REF] Best | Abiotic Guanidinium Containing Receptors for Anionic Species[END_REF] are the most common. If connected to a suitable chromophore, the anion binding will result in observable change in color as in the case of fluoride anion/pyrrolidylamidothiourea based sensor. [START_REF] Evans | Anion Binding vs. Deprotonation in Colorimetric Pyrrolylamidothiourea Based Anion Sensors[END_REF] In this case, the effect is related to a deprotonation of the acidic thioureas by highly basic anions, like acetate or fluoride; while less basic anions, like nitrate, perchlorate or sulphate, give small, often negligible changes. [START_REF] Wu | Fluoride-Selective Colorimetric Sensor Based on Thiourea Binding Site and Anthraquinone Reporter[END_REF][START_REF] Jimenez | Selective Fluoride Sensing Using Colorimetric Reagents Containing Anthraquinone and Urea or Thiourea Binding Sites[END_REF][START_REF] Nishizawa | Anion Sensing by a Thiourea Based Chromoionophore via Hydrogen Bonding[END_REF] Guanidine functional group is especially interesting due to its acid/base properties. The specific Y-shaped structure of its protonated form allows efficient delocalization of the electrons, which results in high stability and, therefore, in its unusually high basicity. [START_REF] Margetić | Superbases for Organic Synthesis: Guanidines, Amidines, Phosphazenes and Related Organocatalysts[END_REF][START_REF] Hall | Relations between the Structure and Strength of Certain Organic Bases in Aqueous Solution[END_REF] Unlike urea and thiourea, guanidine can be easily protonated and, in this form, anion binding is assisted by coulombic charge-charge attractive forces. Therefore, we can expect that guanidines bind anions significantly stronger than ureas, thioureas or amidines. This latter feature of guanidinium cation has been widely used in design of strong and selective binders for anionic guests, [START_REF] Blondeau | Molecular Recognition of Oxoanions Based on Guanidinium Receptors[END_REF][START_REF] Houk | Abiotic Guanidinium Receptors for Anion Molecular Recognition and Sensing[END_REF] anion transporters, colorimetric probes, and sensors. [START_REF] Bose | Functionalized Guanidinium Chloride Based Colourimetric Sensors for Fluoride and Acetate: Single Crystal X-ray Structural Evidence of -NH Deprotonation and Complexation[END_REF] Moreover, the large scientific interest on guanidine species is motivated by its presence as a substructure in many natural compounds and commercial drugs. [START_REF] Sączewski | Biological Activities of Guanidine Compounds, 2008 -2012 Update[END_REF] A number of basic organocatalysts have been designed using guanidine as the pivotal building block. [START_REF] Ferris | Synthesis of Pyrimidines from Guanidine and Cyanoacetaldehyde[END_REF][START_REF] Robertson | Prebiotic Synthesis of Diaminopyrimidine and Thiocytosine[END_REF] Besides that, plate-like structure of the protonated guanidine shows two opposite tendencies: its hydrophobic upper and lower parts allow even cation-cation stacking with practically no repulsion, [START_REF] Vazdar | Like-Charge Ion Pairing in Water: An Ab Initio Molecular Dynamics Study of Aqueous Guanidinium Cations[END_REF] while its hydrophilic edge, defined by the directionality of N-H bonds, tends to form strong hydrogen bonds with polar molecules especially with those of similar geometry, like nitro or carboxylate groups. Although guanidine subunit is of the great potential, it also complicates the design of the receptor due to the increased number of the structural variables one should take care of.

Scheme 1. Selected aromatic guanidines

Development of modern quantum chemical calculation programs and methods allows us to

analyze structure of the anion receptors in detail and to gain insight on changes that occur upon anion binding. This information can be of great importance for the design of the efficient anion sensors. In this work, we explored anion binding properties of some simple aromatic guanidines (Scheme 1) and investigated structures and UV/Vis properties of the free receptor and its complexes with various anions. For this purpose, we carried out a comprehensive investigation of simulated and measured UV spectra for this class of molecules. Special emphasis is put on the orbital interactions between the guanidine subunit and chromophore identifying auxochromic properties of guanidine. These calculations have also served another purpose: to test how successfully DFT models can be employed for predicting UV/Vis properties of components for

QGU

the light-responsive systems, where photoenergy transfer is regulated through guanidine moiety and for the potential building blocks in construction of novel guanidine based anion sensors.

EXPERIMENTAL AND THEORETICAL METHODS

Computational details. Gas-phase ground state geometries were optimized at the B3LYP/6-311+G(d,p) level of theory. [START_REF] Becke | Density-Functional Thermochemistry. III. The Role of Exact Exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density[END_REF] Vibrational analyses were performed to verify the true minimum nature. Several orientations of the substituent group were tested and the lowest energy structure was further used for the calculations of the excited-state properties. The electronic absorption spectra of phenylguanidine (PHGU), were calculated by means of the time-dependent functional theory (TD-DFT) [START_REF] Runge | Density-Functional Theory for Time-Dependent Systems[END_REF] using different functionals: BMK, 33 PBE0, 34 M062X, [START_REF] Zhao | The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals[END_REF] B97XD, [START_REF] Chai | Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections[END_REF] LC-BLYP, [START_REF] Becke | Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior[END_REF][START_REF] Iikura | A Long-Range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals[END_REF] LC-PBE, [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Vydrov | Importance of Short-Range Versus Long-Range Hartree-Fock Exchange for the Performance of Hybrid Density Functionals[END_REF] and CAM-B3LYP. [START_REF] Yanai | A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP)[END_REF] For a recent review on the performance of different density functionals for calculation of the excitation energies in different molecules see ref. 42 . In addition to the default value of the  parameter in the long-range corrected functionals, calculations with =0.2 a0 -1 were performed. A tuning of  parameters for different DFT functionals were conducted for different systems and the optimized values were in the range 0.17-0.21 a0 -1 . [START_REF] Zhang | Theoretical Study of the Local and Charge-Transfer Excitations in Model Complexes of Pentacene-C60 Using Tuned Range-Separated Hybrid Functionals[END_REF] We thought it was worthwhile to test it herein on our molecules and, indeed, it was found that there is better agreement with Emax from measured spectra, if the decreased  = 0.2 a0 -1 was used not only for LC-BLYP and LC-PBE, but also for CAM-B3LYP. All calculations have been performed by using 6-311+G(d,p) basis set.

The gas-phase TD-DFT results were compared with corresponding complete active space perturbation theory to the second order (CASPT2) results based on the aug-cc-pVDZ basis set. [START_REF] Andersson | Second-Order Perturbation Theory with a CASSCF Reference Function[END_REF][START_REF] Andersson | Second-Order Perturbation Theory with a Complete Active Space Self-Consistent Field Reference Function[END_REF][START_REF] Dunning | Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen[END_REF] The active space consisted of 14 electrons and 11 orbitals. State averaging (SA) procedure at complete active space self-consistent field (CASSCF) and multi-state (MS)

CASPT2 were used for calculating vertical excitation energies with up to 10 roots. The CASPT2 vertical energies were calculated by using level shift of 0.3 a.u. [START_REF] Forsberg | Multiconfiguration Perturbation Theory with Imaginary Level Shift[END_REF] and default IPEA shift (0.25 a.u.). [START_REF] Ghigo | A Modified Definition of the Zeroth-Order Hamiltonian in Multiconfigurational Perturbation Theory (CASPT2)[END_REF] The same shifts have been used for the calculations of the excited states of guanidine. [START_REF] Antol | Guanidine and Guanidinium Cation in the Excited State -Theoretical Investigation[END_REF] Also, in all CASPT2 calculations, the core electrons (20 electrons, 10 orbitals) were kept frozen and the oscillator strengths were calculated with the CAS state interaction method (RASSI). [START_REF] Malmqvist | The CASSCF State Interaction Method[END_REF] The RASSI calculations were carried out following the MS-CASPT2 calculations.

Calculated UV spectra of the guanidines and their salts at the TD-DFT level were further refined by considering the molecules into acetonitrile through the polarizable continuum model (PCM).

The standard integral equation formalism (IEFPCM) and default atomic radii, [START_REF] Miertuš | Electrostatic Interaction of a Solute with a Continuum. A Direct Utilizaion of Ab Initio Molecular Potentials for the Prevision of Solvent Effects[END_REF] as implemented in Gaussian09, [START_REF] Frisch | Gaussian 09[END_REF] were used in conjunction with PBE0 and CAM-B3LYP ( = 0.2 a0 -1 ) functionals.

The binding energy between guanidines (protonated form) and different anions has been calculated as the difference between the total energy of the salt and the sum of the energies of corresponding isolated guanidinium cation and anion. The Boys-Bernardi counterpoise technique [START_REF] Boys | The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors[END_REF][START_REF] Simons | Walking on Potential Energy Surfaces[END_REF] has been used to correct the inherent basis set superposition error (BSSE).

The DFT and TDDFT electronic data were obtained using Gaussian 09 program. [START_REF] Frisch | Gaussian 09[END_REF] The CASPT2/CASSCF calculations were performed with Molcas 7.8 software, [START_REF] Aquilante | MOLCAS 7: The Next Generation[END_REF][START_REF] Veryazov | 2MOLCAS as a Development Platform for Quantum Chemistry Software[END_REF][START_REF] Karlström | MOLCAS: a Program Package for Computational Chemistry[END_REF] while Vega-ZZ [START_REF] Pedretti | a Versatile Program to Convert, Handle and Visualize Molecular Structure on Windows-based PCs[END_REF] and Molden 59 programs were used for visualization and geometry manipulations.

Experimental details. Phenylguanidine 60 and 1-naphthylguanidine were prepared by converting appropriate amines to benzoylthioureas, [START_REF] Rasmussen | Improved Procedures for the Preparation of Cycloalkyl-, Arylalkyl-, and Arylthioureas[END_REF] followed by guanylation using hexamethyldisilazane method [START_REF] Shinada | A New Entry for the Synthesis of N-Acyl-N′-substituted Guanidines[END_REF] and finally by debenzoylation and deprotonation in 30 % NaOH [START_REF] Rasmussen | Improved Procedures for the Preparation of Cycloalkyl-, Arylalkyl-, and Arylthioureas[END_REF][START_REF] Shinada | A New Entry for the Synthesis of N-Acyl-N′-substituted Guanidines[END_REF] (Scheme 2).

Synthetic details are given in the Supporting Information. Nitrate, acetate and formate salts were prepared by adding equimolar amount of the corresponding acid to the methanolic solution of the desired guanidine and evaporating to dryness. Formed crude salts were recrystallized from the acetonitrile.

UV/Vis spectra were recorded at 25°C using PG Instruments Ltd T80+ spectrophotometer with chemstation software. All spectroscopic studies described here were conducted in acetonitrile.

RESULTS AND DISCUSSION

Phenylguanidine. The study of auxochromic properties of guanidinium group was started by its attachment to the simplest possible aromatic ring, benzene. Firstly, the electronic states of phenyl-guanidine (PHGU) (Scheme 1) have been characterized in the gas phase and in solvent at the TD-DFT level of theory, varying the density functionals. The gas phase results were compared with CASPT2 calculations (the vertical excitation energies, oscillator strengths and leading configurations, as well as figures of active orbitals at CASPT2 are given in the Supporting Information), while the inclusion of surrounding solvent effects in the DFT calculations allows direct comparison of simulated and experimentally measured UV/Vis spectra. Acetonitrile solvent has been selected because specific solute-solvent interactions are not expected there and PCM [START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] provides a valid approximation for the solvent effects. To study UV/Vis spectra, the nonequilibrium PCM solutions for the TD-DFT calculations were applied. [START_REF] Cossi | Time-Dependent Density Functional Theory for Molecules in Liquid Solutions[END_REF] The calculations in solutions were performed using the PBE0 and CAM-B3LYP functionals. 1). The excitation energies calculated by TD-DFT methods are between 5.05 (LC-BLYP(=0.2)) and 5.54 eV (LC-BLYP); and the calculated CASPT2 values are in the same range.

Based on the case of PHGU as well as on additional data for akin PHTU and PHUR (details in Supporting Information) we highlight that the decrease of the  parameter in the long-range corrected functionals leads to a systematic decrease of excitation energies for all examined excited states. The decrease is the largest for the LC-BLYP and the smallest for CAM-B3LYP.

The comparison between measured and theoretically predicted excitation energies for the first bright state showed that the smallest deviations from the experimental value were obtained with the PBE0 functional. The Eexc for PHGU is in excellent agreement (the difference is less than 0.1 eV). The excitation energy for PHUR is overestimated by 0.22 eV, and the excitation energy of PHTU is underestimated by 0.12 eV. In the case of CAM-B3LYP all values are systematically overestimated by ca. 0.4 eV. The systematic overestimation is reduced to 0.2 eV by the replacement of the default  parameter with the value 0.2 a0 -1 , which is typically the expected shift between the vertical excitation and the band maximum. [START_REF] Barbatti | The UV Absorption of Nucleobases: Semi-Classical Ab Initio Spectra Simulations[END_REF] Therefore, it can be concluded that substitution effect on the absorption spectrum of PHGU is reasonably well reproduced by PBE0 and CAM-B3LYP(=0.2) methods.

Detailed inspection of the orbitals located at the two fragments in question (guanidine and chromophore) indicate the n- conjugative interaction between them. In a free neutral guanidine, the three highest occupied orbitals have electron density distributed within GU and nGU lobes, as

shown in Figure 1. The interaction with benzene orbital 2e1g is the most intense, leading to the destabilization of the HOMO. Since the energies of LUMO orbitals in benzene and in PHGU are almost equal, the HOMO-LUMO gap in PHGU is smaller than in benzene resulting in bathochromic shift. The guanidinium group is highly basic (pKa of guanidine in water: 13.6 kcal mol -1 ) [START_REF] Raczyńska | Consequences of Proton Transfer in Guanidine[END_REF] and it is easily protonated under mild conditions. The optimization of the structure of protonated PHGU predicts that the guanidinium group is more perpendicular with the plane defined by phenyl ring, with respect to parent PHGU, minimizing the interactions between two  systems. For instance, the C8-N7-C1-C2 dihedral angle is increased from 64.5° in PHGU to 78.8° in protonated PHGU.

The C1-N7 and N7-C8 bonds are elongated by 0.044 Å and 0.50 Å, respectively. The spectrum of protonated guanidine in acetonitrile solvent was calculated with the CAM-B3LYP(=0.2) and PBE0 functionals and it is convoluted with normalized Gaussian functions with 0.2 eV standard deviation. A comparison with neutral PHGU is given in Figure 2a and Figure 2b Protonation of guanidine subunit leads to removal of the nGU- interaction by converting nGU orbital to a -bond, thus diminishing interactions between fragments. Indeed, Mayer bond orders, calculated for the guanidine-chromophore junction bond in PHGU, decrease from ca 1.2 to ca 1.0 upon protonation. In line with that, second order perturbation analysis within the NBO basis shows electron density donation from nGU lone pair to the *(C1-C2) NBO located at the chromophore contributing with 16.9 kcal.mol -1 and this interaction is practically removed upon protonation. Additionally, protonation of the guanidine subunit reduces the extent of chr-GU interaction by stabilizing GU orbital due to an electrostatic effect and translation of guanidine orbitals from C1 to C3 symmetry group (if considered as a separate fragment). Besides that, torsional angle between chromophore and guanidine planes is practically 90°, indicating complete absence of conjugation. Consequently, the frontier orbitals as well as UV spectra of protonated forms look more like for the isolated chromophore.

To check the consistency of theoretical prediction, 5×10 -5 M solution of PHGU in acetonitrile was titrated by 5×10 -4 trifluormethanesulfonic acid. After addition of one equivalent of the acid, UV spectra did not change anymore, indicating complete protonation of the neutral guanidine and absence of any other protonation site. The spectra recorded as a titration starting and end point are shown in Figure 2c. The blue shift between the neutral and the protonated species is 21 nm (0.48 eV). Thus, the calculated spectra for the neutral and protonated phenylguanidine are in fair accordance with the measured ones considering the positions of max. The calculated ratios between oscillator strengths of first bright excited state in neutral and in protonated form f(PHGU)/f(PHGUH + ) are 2.6 and 1.5 for CAM-B3LYP(=0.2) and PBE0 methods, respectively.

The latter is in excellent agreement with experimentally obtained value: the area of the first experimental band is 0.399 for the neutral and 0.266 for the protonated giving the ratio of 1.5. As a next step in investigating auxxochromic effect of guanidinium group on the chromophore, we have analyzed two processes that occur upon hydrogen bond formation between guanidinium cation and anion: charge compensation and deprotonation of the guanidine subunits. We have approached this problem by computing the absorption spectrum for salts formed from complexation of PHGUH + with Fˉ, HCOOˉ, CH3COOˉ, NO3ˉ, and CF3SO3ˉ. Since both PBE0

and CAM-B3LYP(=0.2) functionals gave similar results, only the calculations performed with TD-DFT CAM-B3LYP(=0.2) in acetonitrile will be discussed further in this text. Complete set of the PBE0 data (Eexc, f, leading configurations, XYZ coordinates and comparison with experimentally available data) can be found in Supporting Information.

Two different sites of anion bonding were examined (Scheme 3). An anion can approach the guanidinium group forming a H-bonding through N7H and NH' bonds (complexation site 1) or through N10H and N9H'' bonds (complexation site 2). Moreover, different monovalent anions were considered varying their basicity: Fˉ > HCOOˉ ≈ CH3COOˉ > NO3ˉ > CF3SO3ˉ. In all cases, the structures where an anion interacts through the complexation site 1 are more stable.

Therefore, only bonding of different anions to PHGU in the site 1 was used in all subsequent TD-DFT calculations. The results of these calculations are shown in Figure 3 and Table 2. All investigated PHGU salts have the first absorption band between 210 and 250 nm and they are blue shifted with respect to the first band in parent PHGU. The shift is smaller than the corresponding shift for protonated PHGU. The value of max increases from 220 nm for CF3SO3ˉ salt to 231 nm for Fˉ salt, in the sequence CF3SO3ˉ < NO3ˉ < HCOOˉ < CH3COOˉ < Fˉ. Also the calculated oscillator strength increases from 0.244 for CF3SO3ˉ to 0.371 for Fˉ salt. It should be stressed that these shifts correlate to the anion basicity and binding energies of PHGUH + with anions (Table 2). The highest binding energy was calculated for fluoride and acetate anions, which are the most basic anions in the considered set. In spite of the relatively high basicity of these two anions, optimizations in acetonitrile do not predict full proton transfer from the guanidinium cation to the anion; although a significant N-H bond elongation takes place (Table 2). Thus, N7-H bond length assumes values of 1.064 and 1.057 Å for fluoride and acetate complexes, which are 0.050 and 0.043 Å larger than for the isolated phenylguanidinium cation. The analogous bond elongation induced by nitrate binding amounts to only 0.019 Å. This partial deprotonation is also weakest for acidic anion CF3SO3ˉand its spectrum is similar to the spectra of protonated PHGU.

For highly basic anions, as Fˉ and CH3COOˉ where partial deprotonation is more emphasized, their maxima of the first band are shifted toward neutral PHGU. a) with respect to protonated PHGU b) calculated as the difference with respect to the non-coordinated guanidinium cation PHGUH + c) pKa values in acetonitrile are given in parentheses. [START_REF] Muckerman | Toward the Accurate Calculation of pKa Values in Water and Acetonitrile[END_REF] To shed more light on the extent of UV-absorption changes caused by guanidine/ion interaction, we followed the reaction paths for three processes: (i) stepwise approach of the acetate anion, (ii) proton exchange between phenylguanidinium cation and acetate anion within ion pair, and (iii) departure of the neutral acetic acid. Process (i) encompasses partial charge compensation (shielding) with concomitant partial deprotonation of the guanidine subunit. These two processes were separated by comparison of two scans for the approach of the acetate anion toward the protonated phenylguanidine: relaxed and "semi-rigid" (partially relaxed). The "semi-rigid" scan was conducted by fixing the geometry of the guanidinium subunit that corresponds to the optimized free phenylguanidinium cation. In this case, no elongation of the NH bonds was allowed preventing thus partial deprotonation process. The scans were carried out at B3LYP/6-31G(d) level of theory, while UV spectra were calculated using CAM-B3LYP(=0.2) at seven selected partially optimized scan points (Figure SI4 in Supplemental Material). [START_REF]partial optimization of each point along the scan causes a small bathochromic shift[END_REF] The results were compared to the fully relaxed scan in which partial deprotonation takes place. Both scans were started from the optimized structure of the phenylguanidinium acetate complex and acetate anion was gradually displaced until distance between C10 and carboxylate carbon atom reached 6 Å. The acetate induced bathochromic shift was calculated as the difference in max between the first and last point of the scan. Both scans predicted similar anion induced bathochromic shift amounting to 7 and 6 nm for the relaxed and semi-rigid scan, indicating minor influence of the partial deprotonation. Since the UV spectrum of the phenylguanidinium acetate is in good agreement with the experimental spectrum, we conclude that the acetate induced batochromic shift is primarily due to the electrostatic interaction between two charged species.

The second minimum of the phenylguanidinium acetate complex, which is formed by association of two neutral species (non-ionized complex, Scheme 4), was also optimized and the calculated first max is at 249 nm, which is significantly higher than the experimental value (236 nm). Thus, this structure does not contribute significantly to the UV spectra of the phenylguanidinium acetate salt.

N N N H H H H O O CH 3 Ar H N N N H H H H O O CH 3 Ar H ionized complex (ic)
non-ionized complex (nic) Scheme 4. Two minima of the guanidinium acetate salts.

Extension of aromatic ring: naphthalene and anthracene.

The simulated spectra of 1naphthylguanidine (NGU, structure shown in Scheme 1) is characterized by three bands (Figure 4a). A very weak absorption is featured around 305 nm from HOMO-LUMO transition. There are two intense bands located at 230 and 203 nm, which are strongly overlapped. The former band is dominated by single transition from the ground state to the excited state with Eexc = 5.39

eV and f = 0.6561 (leading configuration is 47a-50a, orbitals are shown in Supporting Info).

Excitations to several excited states with calculated oscillator strengths between 0.1 and 0.25 contribute to the third band (maximum at 203 nm). In the protonated form, the first excited sate transition (HOMO-LUMO, Eexc = 4.41 eV) is shifted to higher energies by 0.35 eV (24 nm) with respect to the parent NGU. The intensity of the related band decreases, in accordance to the observed changes observed upon protonation of PHGU. In the region between 190 and 250 nm, there is only one band, mostly dominated by a transition from the ground state to the third excited state (HOMO-LUMO+1, Eexc=5.70 eV, f=1.288).

The theoretical predictions are compared with the measured UV-Vis spectra of NGU and its protonated form (Figure 4b). Clearly, there is a qualitative agreement between theoretical predictions and measured data. Experimentally observed hypsochromic shift of the lowest energy band upon protonation (24 nm) is in excellent agreement with predictions. Two overlapping bands were measured for the neutral NGU between 200 and 250 nm, as predicted by the calculation. However, in the experiment overlapping is more pronounced (max are closer) than predicted and the intensity of these two bands differs substantially. The rise of sharp band upon protonation is well defined experimentally and theoretically (max= 222 nm (measured) and 215 nm (calculated)). The absorption spectra of 1-naphthyl-guanidinium salts with Fˉ, HCOOˉ, CH3COOˉ, NO3ˉ, and CF3SO3ˉ anions were also calculated. Structural relationship among aromatic system, guanidinium group and anion, as well as complexation site, is identical to what was found in phenyl-guanidinium salts. However, inspection of the simulated UV/Vis spectra (Figure 5) reveals that the effect of complexation on the max is less pronounced than in the case of PHGU.

The position of the HOMO-LUMO band maxima for salts with CF3SO3ˉ and NO3ˉ anions is almost identical to the max calculated for protonated NGU. Calculated Eexc for the first excited state in protonated NGU and its salts with CF3SO3ˉ, and NO3ˉ are 4.41, 4.40, and 4.41 eV, respectively (Table 3). In the case of complexes with HCOOˉ, CH3COOˉ and Fˉ, the HOMO-LUMO transition is shifted toward lower energies. The shift is decreased with respect to the complexation effects found in PHGU. Further extension of the aromatic system to three condensed aromatic rings was studied on the 1-anthracenylguanidine (ACGU, Scheme 1). The simulated spectra of neutral ACGU, its protonated form and different antracenyl-guanidinium salts are shown in Figure 6. Excitation energies for the first excited state (HOMO-LUMO) are collected in Table 3. The overall pattern of the spectrum for ACGU calculated in acetonitrile (Figure 6a) is similar to the spectrum of PHGU. There are two bands: one very intense between 220 and 270 nm (max=252 nm) with the main contribution from the excitation to the sixth excited singlet state (f = 1.515), and the other very weak between 360 and 410 nm (max = 385 nm, f = 0.1442, HOMO-LUMO transition).

Variations in max values caused by protonation are smaller than in PHGU. The effect of the complexed anions on the UV spectra significantly decreases from phenyl to anthracenyl derivative. Thus, the first absorption band shift upon acetate binding amounts to 1.1 eV (phenyl), 0.12 eV (naphthyl), and 0.02 eV (anthracenyl), which corresponds to 64, 35, and 15 % of the protonation induced shift (difference between neutral and protonated form). (Table 3).

The results for the phenylguanidine overestimate the experimental results, where acetate induced bathochromic shift contributes 44% of deprotonation process. Evidently, anion induced batochromic shift is due to formation of the complex and not to the deprotonation of the guanidine. This is somewhat surprising having in mind the relatively high basicity of the acetate anion in acetonitrile (pKa(CH3COOH) = 23.5 [START_REF] Muckerman | Toward the Accurate Calculation of pKa Values in Water and Acetonitrile[END_REF] ). In contrast, the experimental data for NGU and its salts show the maximum of the first band either around 305 nm (neutral, acetate and fluoride)

or around 281 nm (protonated form, nitrate and trifluormethansulfonate).

This finding corresponds to the two-state situation: protonated and deprotonated form depending on the basicity of the anions. Measured pKa of the 1-naphthylguanidine is 19.4±0.2, which is ca 0.5 lower than that of phenylguanidine. Therefore, it is of no surprise that the former guanidine undergoes deprotonation by acetate and fluoride more readily than the latter one. The first band in UV spectra of the formate salt is slightly shifted and broadened with respect to the nitrate (protonated system) and centered at 287 nm. This position corresponds to the hydrogen bonded ion pair in which no proton transfer from the guanidinium cation to the formate anion yet happened (ionic complex). It is well known that the formate anion is less basic than the acetate 69 and the difference in their basicities (Table 2) is apparently sufficient to regulate preferential structure upon anion binding. The other investigated guanidines are even less basic than 1naphthylguanidine and, similarly, full deprotonation by fluoride, acetate and formate is expected, while other anions are significantly less basic and most likely will give the spectra of the protonated form.

Quinoline, anthraquinone and cumarine derivatives. Quinoline, anthraquinone and coumarine are heteroaromatic chromophores which have the possibility of H-bond interactions with guanidine subunit attached to C1 carbon atom of chromophore (Scheme 1). Indeed, in the optimized structures, guanidine subunit is oriented parallel to the chromophore aromatic system.

For instance, in 8-quinolinylguanidine (QGU, Scheme 1) the dihedral angle C10-N9-C1-C2 amounts to 0.5°, while corresponding dihedral angles in CGU and AQGU are 3.0° and 1.6°, respectively. Moreover, the N9-H bond is tilted toward the nitrogen lone pair on quinoline subunit, establishing classical hydrogen bond, while the imino nitrogen atom is oriented in opposite direction, allowing for non-classical C2-H...N11 H-bonding interaction. In fact, atoms in molecules (AIM) analysis of the quinolinylguanidine identifies (3,-1) bond critical point between N12 and C2-H (Figure SI5) with the values of 0.07 and 0.024 for the b and  2 b, respectively, values typical for the non-bonding component of the hydrogen bonding interaction. [START_REF] Koch | Characterization of C-H-O Hydrogen Bonds on the Basis of the Charge Density[END_REF] Surprisingly, we were unable to locate the corresponding bond critical point that would confirm C-H...N hydrogen bonding in the 1-guanidylathraquinone.The same orientation of guanidine group is found in both anthraquinone and coumarine derivative.

Protonation at the guanidine subunit removes the non-classical hydrogen bonding interaction and increases the repulsion between guanidinium and chromophore moieties, which in turn leads to the increase in the dihedral angle C10-N9-C1-C2 to 29.2, 34.6, and 60.8° in QGUH + , AQGUH + , and CGUH + , respectively. This is not surprising since the amino N9-H…O hydrogen bond in coumarine is expected to be the weakest of all considered intramolecular hydrogen bonds due to the interaction with ether-type oxygen atom. Consequently, complexation with anions in quinoline and anthraquinone goes into position 2 and, in the case of coumarine salt, which has a structure similar to what was found in PHGU, NGU and ACGU, the anion sits in position 1 (for example, see CGUHF in Scheme 5). Both anthraquinoyl-and quinolinyl-guanidines could also form "site 1" complexes at expense of hydrogen bond cleavage, but formation of these complexes is less likely. Difference between "site 1" and "site 2" complexes in anthraquinone derivatives amounts to ca 2 kcal mol -1 , as calculated at the CAM-B3LYP(=0.2) level of theory.

Since the difference corresponds to 97 : 3 ratio according to the Boltzman's distribution, we do not expect significant contribution of second isomer to the absorption spectrum. In the calculated spectrum of CGU (Figure 7e), two maxima can be observed as the main feature (max = 284 nm and max = 240 nm). Absorption associated to excitation from the ground state to the first excited state (HOMO-LUMO transition) does not appear as a separate weak band as found for all other guanidines. Instead, it contributes to a long tail (360 >  > 310 nm) connected with the stronger band positioned at 280 nm. Protonation shifts the HOMO-LUMO transition to the higher energies and it is completely hidden under the strong band with max = 277 nm. Very small difference in UV spectra of neutral and protonated CGU renders this compound practically insensitive to the presence of anions in the region 300-350 nm (Figure 7e).

As was already discussed, the most stable tautomer of the neutral form in QGU, AQGU and CGU has imino nitrogen on one of two distant nitrogen atoms within the guanidine subunit, while the -nitrogen bears the proton involved in the intramolecular hydrogen bond. The guanidine is coplanar with the chromophore, while the lone pair at the imino nitrogen atom is directed toward H(C2) hydrogen atom, additionally stabilizing the structure by N...HC unclassical hydrogen bond. Since the lone pair lies on the  plane, the dominant interaction is between the highest occupied fragment orbital of the chromophore and the  orbital of guanidine subunit (Figure 8). Fragment orbital interaction is similar as discussed for the aryl guanidines (PHGU, NGU ACGU), with destabilization being somewhat stronger in the case of anthraquinone and quinoline derivatives with respect to their naphthalene and anthracene counterparts. We note that the HOMO of 1-guanidylanthraquinone is strongly localized within the ring bearing the substituent, while the protonation induces an electron density shift toward the farthest anthraquinone ring. On the other hand, no such dramatic redistribution of the electron density was observed for the LUMO. The localization of the HOMO in one ring implies large orbital coefficients around junction bond and makes it sensitive to the change in electronic properties of guanidine subunit. Namely, neutral guanidine acts as electron donating, while the guanidinium cation acts as electron withdrawing substituent. [START_REF] Leffek | Basicity of Substituted 2-Phenyl-1,1,3,3-tetramethylguanidines and Other Bases in Acetonitrile Solvent[END_REF] Consequently, protonation of 1guanidylanthraquinone will stabilize the HOMO significantly more than the LUMO, resulting in the aforementioned largest hypsochromic effect within the investigated set of arylguanidines.

As discussed earlier, anion binding of the anthraquinone and quinoline derivatives occurs most likely at the "site 2" position and therefore does not induce deprotonation at the N9 of the guanidinium fragment. Nevertheless, the anion binding effect on the UV spectra is only slightly smaller than for the naphthylguanidine derivative. Examination of the changes in the partial charges in AQGU upon acetate and nitrate binding indicates no conformational dependence on the extent of electron density transfer. Regardless of whether guanidine is in perpendicular or almost co-planar conformation with respect to the chromophore plane, the decrease in the partial charge amounts to 0.17 and 0.10 |e| for acetate and nitrate binding, respectively, indicating a minor effect of the guanidine conformation on the UV spectra. For these derivatives, binding of fluoride, acetate and formate anions is expected to lead to deprotonation of the guanidinium cation subunit resulting in regaining UV spectra of the neutral guanidine CONCLUSIONS UV/Vis spectra of phenylguanidine (PHGU) in the gas phase and acetonitrile have been simulated by TD-DFT calculations. Several DFT hybrid and long-range corrected functionals (BMK, PBE0, M062X, B97XD, LC-BLYP, LC-PBE, and CAM-B3LYP) were tested with respect to CASPT2 calculations for the gas-phase. In addition, comparison with isoelectronic phenyl urea and related phenyltiourea were carried out. Results indicate that PBE0 and longrange functional CAM-B3LYP method with the  parameter modified to 0.2 a0 -1 provide consistently good agreement between experimental and calculated UV spectra for neutral and protonated phenyl and naphthyl guanidines, as well as for their salts with CH3COOH, HCOOH, HNO3, HF and CF3SO3H.

UV spectra were also calculated for the series of aryl guanidines, in which the influence of the guanidine moiety on the UV spectra of the selected basic chromophore substructure (anthracene, anthraquinone, quinoline, and coumarine) was analyzed. All investigated guanidine derivatives showed HOMO-LUMO transition to be sensitive to protonation/deprotonation process. The observed changes in UV spectra upon protonation were interpreted in terms of changes in orbital interactions, conformational changes and presence/absence of intramolecular hydrogen bonds.

It was shown that the preferential binding site of anion involves interaction with the most acidic NH bond, which is the one closest to the aromatic moiety unless it is involved in the intramolecular hydrogen bond. Different response of the phenyl and naphthylguanidine toward acetate anion was observed, indicating formation of the ionic complex in the former case and deprotonation in the latter. That was attributed to the differences in pKa and stability of these two modes of the salt substructure. The less basic formate anion forms an ionic complex with 1naphthylguanidine without deprotonation, as deduced from the comparison of the experimental and calculated UV spectra.

Although anion binding induces elongation of the associated NH bonds, calculations predict minimal effect of the "partial deprotonation" to the UV spectra changes. These changes are mostly associated with charge compensation upon anion binding. The second, even larger in its extent, contribution is "neutralization", where complete proton transfer to anion takes place forming nonionic acid-base associate. Although, the dynamic situation is present in solution, we have shown that both situations are distinguishable and could be recognized from the UV spectra.

ASSOCIATED CONTENT Supporting Information. Synthetic procedures and details of pKa measurements, excitation energies, oscillator strengths and leading configurations at CASPT2, CAM-B3LYP(=0.2) and PBE0 levels of theory, CASSCF active orbitals, relaxed and semi rigid scan for the formation of phenylguanidinium acetate salt, AIM calculation on NGU, Cartesian coordinates of optimized structures. This material is available free of charge via the Internet at http://pubs.acs.org.

2 NHScheme 2 .

 22 Scheme 2. Synthesis of the target aromatic guanidines.

Figure 1 .

 1 Figure 1. Orbital interactions in PHGU. The energies were taken from CAM-B3LYP(=0.2)/6-

  by using PBE0 and CAM-B3LYP(=0.2), respectively. These two DFT functionals predict a large blue shift of the first absorption band caused by protonation. The shift is 31 nm (0.75 eV) at the CAM-B3LYP(=0.2) level of theory. It was also predicted that the intensity of the band is decreased by protonation. Still at the CAM-B3LYP(=0.2) level, the calculated oscillator strength of the bright state in neutral PHGU was 0.342 and in protonated PHGUH + it was 0.131, a reduction of f = 0.21. The PBE0 values for the shifts obtained by protonation are very similar: max = 27 nm (0.62 eV) and f = 0.12.

Figure 2 .

 2 Figure 2. Comparison between UV/Vis spectra of neutral and protonated PHGU: (a) Theoretical

Scheme 3 .

 3 Scheme 3. Protonated phenylguanidine and its salts with Fˉ, COOHˉ, CH3COOHˉ, NO3ˉ and

Figure 3 .

 3 Figure 3. Comparison between simulated UV/Vis spectra (CAM-B3LYP(=0.2) method) of

69

 69

Figure 4 .

 4 Figure 4. (a) The simulated absorption spectra of NGU and NGUH + in acetonitrile at the TD-

Figure 5 .

 5 Figure 5. Comparison between simulated UV/Vis spectra by using CAM-B3LYP(=0.2)

Figure 6 .

 6 Figure 6. (a) The simulated absorption spectrum of ACGU in acetonitrile at theTD-CAM-

Scheme 5 .

 5 Scheme 5. Different complexation sites.

Figure 7 .

 7 Figure 7. (left) The simulated absorption spectra of aryl guanidines (QGU, AQGU and CGU) in

Figure 8 .

 8 Figure 8. Orbital interactions in QGU. The energies were taken from CAM-B3LYP(=0.2)/6-

Table 1 .

 1 Comparison

		of selected vertical excitation energies (Eexc/eV) calculated with different
	density functionals. CASPT2 results are also provided. Oscillator strengths (f) are given in
	parentheses.
	Method	1 st exc. state		Bright state
		Eexc/eV	f	Eexc/eV	f
	Gas phase
	BMK	5.14	0.020	5.37	0.283
	PBE0	4.91	0.019	5.16	0.267
	M062X	5.12	0.019	5.35	0.281
	B97XD	5.07	0.018	5.36	0.281
	LC-BLYP	5.25	0.016	5.54	0.252
	LC-BLYP(ω=0.2)	4.76	0.017	5.05	0.242
	LC-PBE	5.21	0.016	5.53	0.274
	LC-PBE(ω=0.2)	4.80	0.016	5.10	0.290
	CAM-B3LYP	5.06	0.018	5.32	0.275

(36a) to LUMO+1 (38a). The second excited state has HOMO-LUMO (36a-37a) character. Substitution of benzene with guanidinium group gives rise to small but nonzero oscillator strengths in contrast to the parent benzene, where f for corresponding excited states is zero. The calculated excitation energy of the first excited state is 4.67 eV using multi-state CASPT2 approach. All DFT methods give higher excitation energies up to 5.25 eV (LC-BLYP). PBE0, LC-PBE(=0.2) and CAM-B3LYP(=0.2) give excitation energy values lower than 5 eV being in better match with CASPT2 values. The first observable absorption band in the spectrum of PHGU is produced by transition from the ground state to the second excited state. Its calculated oscillator strength is predicted to be 0.505 at CASPT2 level of theory. TD-DFT gives lower intensity of the first bright state than the CASPT2. The f value varies from 0.129 to 0.290, as calculated with different density functionals (Table

Table 2 .

 2 Comparison of the anion binding energy (BE) to protonated PHGU with the shift of the first absorption band of different salts with respect to the protonated PHGU in acetonitrile. The pKa of conjugated acids HA measured in water are given as well.

	Anion	Eexc/eV a	max a	r(N7-H) b / Å	BE/kcal mol -1	pKa(HA) c
	CF3SO3ˉ	0.26	9	0.014	-12.0	-5.9 67 (2.60)
	NO3ˉ	0.35	13	0.019	-14.4	-1.64 68 (8.80)
	HCOOˉ	0.47	18	0.037	-19.9	3.75 69
	CH3COOˉ	0.50	19	0.043	-21.6	4.76 69 (23.51)
	Fˉ	0.53	20	0.050	-21.9	3.17

Table 3 .

 3 Calculated excitation energies Eexc (eV) for the lowest excited state. CAM-

	B3LYP(=0.2) density functional was used in conjunction with 6-311+G(d,p) basis set. Solvent
	= acetonitrile. pKa of the investigated arylguanidinium cations (BH + ) and complexation energies
	between acetate anion and protonated guanidines are given as well.
	Chromophore GU(	CH3COOˉ Fˉ	HCOOˉ NO3ˉ CF3SO3ˉ GU(p	pKa(	BE /	Erel c /
		neutr al)						roton ated)	BH + ) calc	kcal mol -1	kcal mol -1
	Benzene	4.87(	5.16(5.40)	5.12	5.17(5.	4.10(	5.28(5.6	(5.89	20.2	-21.6 -2.60
		5.13)	a	(5.3	42) a	5.54)	3) a	) a	(20.0
		a		7) a		a			) b

c) Relative stabilities (Erel) of two minima of the guanidinium acetate salts (Scheme 4) with respect to the proton movement within the hydrogen bond Erel = Etot(ic) -Etot(nic)

AUTHOR INFORMATION

Corresponding Author

* Electronic mail: iantol@emma.irb.hr

ACKNOWLEDGMENT

The work was supported by the Croatian Science Foundation grant No. 9310. The calculations were performed on the Isabella cluster (isabella.srce.hr) at the Zagreb University Computing Center (SRCE). Funds for Germany-Croatia bilateral collaboration from the German Academic Exchange Service (DAAD, Project 54368738) and Ministry of Science, Education and Sports (Project: Optimizing bidirectional "highway" for photoactive response through guanidinechromophore junction) are also acknowledged. MB thanks the support of the A*MIDEX grant (n° ANR-11-IDEX-0001-02) and of the project Equip@Meso (ANR-10-EQPX-29-01), both funded by the French Government "Investissements d'Avenir" program.


Wilmer Arbelo-González

Rachel Crespo-Otero

Mario Barbatti

Steady and time-resolved photoelectron spectra based on nuclear ensembles

Semiclassical methods to simulate both steady and time-resolved photoelectron spectra are presented. These approaches provide spectra with absolute band shapes and vibrational broadening beyond the Condon approximation, using ensemble of nuclear configurations built either via distribution samplings or nonadiabatic dynamics simulations. Two models to account for the electron kinetic energy modulation due to vibrational overlaps between initial and final states are discussed. As illustrative examples, the steady photoelectron spectra of imidazole and adenine and the time-and kinetic-energy-resolved photoelectron spectrum of imidazole were simulated within the frame of time-dependent density functional theory. While for steady spectra only electrons ejected with maximum allowed kinetic energy need to be considered, it is shown that to properly describe time-resolved spectra, electrons ejected with low kinetic energies must be considered in the simulations as well. The results also show that simulations based either on full computation of photoelectron cross section or on simple Dyson orbital norms provide results of similar quality.

Auger, 3, 9 and X-ray photo-scattering 3 spectroscopies. These developments have been based on a broad range of approximations and electronic structure methods, from very simple estimates of transition probabilities 3, 8c, 8d, 10 to more involved modelling with full computation of transition moments. 8b, 8e, 10c, 11 Here, we present a new implementation of the method aiming at photoelectron spectroscopy, which, although it shares a number of common features with previous implementations (especially those reported in refs. 3, 10a ), distinguishes itself by its generality; allowing simulations of steady and time-resolved spectra, use of arbitrary ensembles, investigation of general molecular systems, control of diverse parameters, choice of particular models for intensity calculation and vibrational overlap modulation, direct integration with any electronic structure method, and computation of intensities with absolute units.

To simulate steady and time-resolved photoelectron spectra we should care of three different aspects: computation of ionization energies, computation of ionization probabilities, and how to put these results together to build the spectrum.

 In view of the flexibility of the  methods to be systematically applied to a large number of points in the ensemble, we have based our current developments on them, even though we had to pay the price of dealing with non-orthogonal sets of orbitals for computation of intensities, as discussed later.

Introduction

To provide results directly comparable to experimental data is a major goal in computational theoretical chemistry. This goal represents a special challenge with the natural trend of dealing with always larger and more complex molecular systems. Whenever nonlocal quantum effects can be neglected, resorting to semiclassical simulations turns out to be a good option, as it allows closely emulating experimental techniques at relatively modest computational costs compared to full quantum simulations.

In the last years, within the development of the NEWTON-X platform, [START_REF] Barbatti | NEWTON-X: A package for Newtonian dynamics close to the crossing seam[END_REF] we have worked out and implemented diverse of such semiclassical approaches for dynamics and spectrum simulations. A central point in these developments has been the extensive use of population sampling via nuclear ensembles. [START_REF] Crespo-Otero | Spectrum Simulation and Decomposition with Nuclear Ensemble: Formal Derivation and Application to Benzene, Furan and 2-Phenylfuran[END_REF] In the present work, this approach will once more play an important role for the implementation of semiclassical methods for simulations of steady and timeresolved photoelectron (PE) spectra.

The nuclear ensemble approach is likely the simplest method to obtain absolute spectral bands. It works in three steps: i) an ensemble of nuclear geometries is built to represent the vibrational distribution in the source state; ii) spectral intensities between the source and the target states are computed for each point in the ensemble; iii) the final spectrum is obtained as an incoherent sum over all these individual transitions. Thus, the vibrational features of the spectrum are first supposed to be essentially dependent on source-state nuclear wavefunction and any property depending on the correlation between the source and target nuclear wavefunctions is neglected (vibrational structures in the electronic spectrum, for instance).

Nuclear ensembles have been chiefly used to simulate steady electronic spectra and to sample initial conditions for dynamics simulations. Going beyond its very intuitive background, we have shown in ref. 2a how the nuclear ensemble approach is correlated to more formal methods.

Recently, Bennett and co-workers, [START_REF] Bennett | Probing electronic and vibrational dynamics in molecules by time-resolved photoelectron, Auger-electron, and X-ray photon scattering spectroscopy[END_REF] presented the approach as a particular case of a unified description of time-resolved spectroscopies. Also, Petit and Subtonik [START_REF] Petit | Appraisal of Surface Hopping as a Tool for Modeling Condensed Phase Linear Absorption Spectra[END_REF] have developed ensemblebased methods for recovering source/target correlation. The impact on the spectrum and dynamics results due to different ways of building the ensemble has been discussed in refs. 2b, 5 .

The nuclear ensemble approach has been applied for simulations of different types of timeresolved spectra, including two-dimensional, [START_REF] Rivalta | Ab initio simulations of twodimensional electronic spectra: The SOS//QM/MM approach[END_REF] stimulated emission, [START_REF] Polli | Conical intersection dynamics of the primary photoisomerization event in vision[END_REF] photoelectron, [START_REF] Bennett | Probing electronic and vibrational dynamics in molecules by time-resolved photoelectron, Auger-electron, and X-ray photon scattering spectroscopy[END_REF][START_REF] Ruckenbauer | Photoelectron spectra of 2-thiouracil, 4-thiouracil, and 2,4-dithiouracil[END_REF] ultrafast Computation of ionization probabilities have a long history dating back to several decades. [START_REF] Starace | Photoionization of Atoms[END_REF] These calculations are not of our direct concern here. Instead, we have used a thirdparty program 16c, 17 to compute ionization cross-sections using standard methods, as explained below. However, for using these methods, Dyson orbitals are needed and their computation within TDDFT frame is also discussed below.

Lastly, having computed ionization energies and probabilities, we must put them together in a spectral representation. In the present work, we do this using the nuclear ensemble approach.

Formal quantum approaches for steady and time-resolved photoelectron spectrum simulations of molecules have also been available for many years 3, 12a, 16d, 18 and their success is well documented. [START_REF] Bennett | Probing electronic and vibrational dynamics in molecules by time-resolved photoelectron, Auger-electron, and X-ray photon scattering spectroscopy[END_REF][START_REF] Cederbaum | Correlation Effects in the Ionization of Molecules: Breakdown of the Molecular Orbital Picture[END_REF] Thus, the approach presented here must not be understood as a new theory aiming at replacing the previous ones. Instead, it should be taken as a routine approach, especially useful in the context of trajectory-based dynamics simulations for large molecules, where ensembles are automatically generated and must be analyzed. We, yet, emphasize that the nuclear ensemble is a low-resolution semi-classical approach. It cannot be expected to compete with formal quantum methods, which are obviously the most indicated option for problems requiring high accuracy.

We will demonstrate here the potential applications and caveats of the nuclear ensemble approach for photoelectron spectroscopy with simulations of the steady PE spectra of adenine and imidazole (Scheme 1) and the time-resolved PE spectrum of imidazole, always in the gas phase.

The choice of these systems followed a pragmatic logic: first experimental data are available for comparison; [START_REF] Devine | High resolution photofragment translational spectroscopy studies of the near ultraviolet photolysis of imidazole[END_REF] and second we have previously studied the photoelectron spectroscopy of both systems, 8d, 15, 21 which will spare us of discussing their very interesting photophysics this time to primarily focus on the method implementation itself.

Finally, we should add that all developments discussed here have been implemented in NEWTON-X and are available in the new versions of the program.

Scheme 1

Theory

Steady PE spectra

Consider the photoionization process depicted in Figure 1 Schematic representation of the photoionization for the steady case.

From a semiclassical standpoint, the probability of occurrence of such a process is proportional to the cross section per unit of electron kinetic energy:

( ) ( ) ( ) ( ) ( ) ( , ) , , , k Im IF k k IF mn F E E d E E E E V K     = -+  +    R R R R R (1)
where Im  is the probability distribution of nuclear coordinates R for the source state I,m , IF  is the photoionization cross section between I and F, and the delta function imposes the energy resonance condition involving the ionization potential

IF F I V V V  = -
and the kinetic energy

6 difference mn n m K E E  = -
between the two states. The sum runs over all target F states contributing to the process.

Supposing that the photoprocess is instantaneous and that the nuclear momentum does not change, 0 mn K  . Additionally, replacing the ( )

x  function by a sharp function, one obtains

( ) ( ) ( ) ( , ) , , , , , k Im IF k s k IF F E E d E E w E E V     -       R R R R (2)
where s w is given either as a normalized Gaussian

  ( ) ( ) ( ) ( ) ( ) 1/ 2 2 2 1 exp for , ,, 2 / 2 0 for , 2 / 2 k IF k IF s k IF k IF E E V E E V w E E V E E V      - --    -  - =      -   (3) or a normalized Lorentzian   ( ) ( ) ( ) ( ) ( ) 2 2 2 /2 1 for , /2 ,, /2 0 for . k IF s k IF k IF k IF E E V w E E V E E V E E V       -   - = -- +    -   (4)
In both cases,  is an arbitrary parameter determining the line width. It should be much smaller than the band width to not interfere with the results, usually 1 eV  is enough to satisfy this requirement.

Nuclear ensemble approach for steady spectra

The integral over R in Eq. ( 2) can be solved by a Monte-Carlo procedure, leading to

( ) ( ) 1 1 ( , ) , , = , , , p N k IF k l s k IF l Fl p E E E E w E E V N  =  -    RR (5)
where a set of In the particular case when the system is prepared before the ionization in the electronic and vibrational ground states 0 I = and 0 m = , it is fair to assume that the harmonic approximation is valid. Under these conditions, it is more natural, 2b as well as numerically efficient, to perform the sampling of the nuclear configurations in the normal mode coordinates q, where the nuclear ensemble is defined by the marginal Wigner distribution function for the quantum harmonic oscillator [START_REF] Case | Wigner functions and Weyl transforms for pedestrians[END_REF] 1/2 2 00 1 ( ) exp .

d N i i i i i i q       =   =-     q (6)
Here, i  is the angular frequency associated with the i th normal mode with reduced mass μi. Nd is the number of normal modes in the system. Once Np nuclear geometries ql are generated according to 0 ()  q , they are transformed back to Cartesian coordinates Rl.

Time-resolved PE spectra

Suppose now that the molecule is at time t = 0 in the electronic and vibrational ground state, with energy E0, when a laser of energy E1 pumps it to the excited state I'. This first excitation is considered fully vertical, so that the nuclear coordinates and conjugate momenta remain constant.

Once in the I' state, the system is allowed to evolve freely and, at t  = , the dynamics is probed by ionizing the molecule with a second laser of energy E2, exactly under the same conditions as in Section 2.1. As nonadiabatic transitions are allowed during the dynamics, the electronic state I at the moment of the ionization, with total energy 01 I E E E =+, may in general be different of I'. This is schematically depicted in Figure 2.

Given the equivalence, the analysis of this time-resolved situation parallels the development of the steady case, considering that the effects of the two laser pulses are uncorrelated. Nevertheless, there are some fundamental differences: first, the nuclear state of the molecule at t = τ is described by a wavepacket, not by a single stationary state; second, the laser pulse duration is in the femtosecond scale, impacting the energy resolution. Bennett et al. [START_REF] Bennett | Probing electronic and vibrational dynamics in molecules by time-resolved photoelectron, Auger-electron, and X-ray photon scattering spectroscopy[END_REF] have shown that the time-resolved case can still be written analogously to Eq. ( 5), but with the initial ensemble distribution given by the population

I   of state I at time .
The photoelectron spectrum is then given by ( )

( ) ( ) ( ) ( ) 2 2 2
( , , ) , , , = ,.

k I IF k s k IF F E E d E E w E E V         -     R R R R (7)
This approximation, which like in the steady case still assumes that the nuclear momentum does not change, implies that the electron is always ejected with the maximum allowed kinetic energy,

max ,2 k IF IF E E V = - . Figure 2.
Schematic representation of the photoionization for the time-resolved case.

Because the molecule at time  is described by a wavepacket rather than by a stationary eigenvector of state I, the Franck-Condon overlaps between I and F are much more complex than in the steady case. 19b Therefore, to assume that the nuclear momentum remains constant during the photo-transition (

IF k IF r k IF k IF E V E E V w E E V E E V -  -  -  - =   -   (8)
With this new assumption, which has also been applied by Fuji et al., 10a the semiclassical expression for the time-resolved spectrum is

( ) ( ) ( ) ( ) ( ) 2 2 2 = ( , , )
, , , .

k I IF k r k IF F E E d E E w E E V        -     R R R R (9)
In the remaining paper, when using sharp s w functions, we will refer to it as the peaked vibrational background (PVB) model; when using rectangular r w functions, we will refer to it as the constant vibrational background (CVB) model.

Nuclear ensemble approach for time-resolved spectra

Either with Eq. ( 7) or with Eq. ( 9), the integral over R is solved by a Monte-Carlo procedure, leading to

( ) ( ) ( ) ( ) 1 1 ( , , ) , , , , = p N k IF k l k IF l Fl p E E E E w E E V N     =   -   RR (10)
where a set of

) ( ) ( ) ( ) ( ) ( ) 1 11 1
.

nn IF k IF k nn IF k IF k k k nn kk EE E E E E EE   - -- - - = + - - (11)

Cross sections

From the light-matter interaction theory up to first order and in the electric-dipole approximation, it is possible to show that the state-resolved photoionization cross section for this process is given by the expression 24 2 0 ( , , ) ( ) ,  also describes the ejected electron with wavevector k . The remaining terms in Eq. ( 13) are the electric dipole operator μ and the unit vector ê in the direction of the electric field of the laser. Integration in Eq. ( 13) is over the electronic coordinates r.

n k IF k IF E E E D c    = RR (12)
Usually, the transition dipole matrix is computed within the Condon approximation at the nuclear equilibrium geometry 0 R . The nuclear ensemble, however, is intrinsically a post-Condon approach, as the transition moments are by construction computed for a distribution of nuclear geometries. For this reason, working equations are derived here implicitly retaining the dependence of

k IF D on R .
Now, assuming the photoelectron ejection is fast, the final electronic state can be represented by the uncorrelated product [START_REF] Oana | Dyson orbitals for ionization from the ground and electronically excited states within equation-of-motion coupled-cluster formalism: Theory, implementation, and examples[END_REF] ,

kk F F F   =  (14) k F
 is the wavefunction of the ejected electron and F  the electronic wavefunction describing the F state of the remaining 1 N -electron species. Also assuming that k F  is orthogonal to the orbitals of the initial state (strong orthogonality conditions), then

2 2 | • | . N k k d IF F IF D  =   r μe (15)
Integration in Eq. ( 15) is over only one electron coordinate, N r in this case.

()

d IF N
 r is the Dyson orbital (DO) associated with the particular I→F transition, formally defined as 25

1 {} ( ) | , N d IF N F I N  - =    r r (16)
where the integration is here over the remaining 1 N -electron coordinates. Note that

d IF
 is defined for a given nuclear configuration R . Introducing the norm of the Dyson orbital d IF  , Eq. ( 15) can be rewritten as

2 22 | • | , N k d k d IF IF F IF D    =   r μe (17)
where /

d d d IF IF IF    = is just the DO normalized to one.
Once the DOs and their norms are known, the right-hand side of Eq. ( 17) can be evaluated.

In this work, after computing the DOs as explained later in Section 2.4, we have used the EZDYSON

( ) ( ) 1/2 1/2 , k e Fk m kF  = rr (18)
where ( ) k F R is the electron wave expansion in a convenient basis and

k F
 is normalized to energy interval (i.e., it has units of (volume × energy) -1/2 ). Thus, the transition moment is

The factor 1/3 stems from the isotropic averaging, while I g accounts for spin and orbital degeneracies of state I. Replacing Eq. ( 19) into Eq. ( 12) renders

0 2 2 3 | | . 3 N dd eI IF I IF k F m kE F g c       =   r (20)
In addition to simulations based on the full computation of the cross sections, we have also simulated the spectrum based on a second approach, which consists of simply approximating Eq.

(

) to 2 , d IF IF C   ( 20
) 21
where C is an arbitrary constant. In this case, all functional dependence of the transition dipole moments on the geometries and final states are supposed to be contained in the DO norms.

In the remaining text, when using Eq. ( 20) to simulate the spectrum, we will refer to it as the cross section approach, and when using Eq. ( 21), as the DO norm approach.

Dyson orbitals

The DO associated with a particular I → F transition (Eq. ( 16)) is a single electron wavefunction containing information on where the ejected electron was removed from. According to our previous developments, once the DOs are known, the cross sections can be evaluated and the photoelectron spectrum fully computed. In the Supporting Information (SI-1), we provide a detailed discussion on how to compute DOs. Here, we will outline only few key aspects.

As shown in the SI, Eq. ( 16) can be rewritten as a linear combination of spin-orbitals q Using Eq. ( 22), the DO norm can be easily computed as

 as 1 ( ) ( ),
1/ 2 2 1 , bf N d IF s s b  =  =     (25)
which, in general, is not equal to one. In fact, d IF  may range from 0 to 1 and, as can be inferred from Eq. ( 21), it is a measure of the photoelectron ejection probability. [START_REF] Ortiz | Electron Propagator Theory of the Photoelectron Spectrum of Methanesulfenic Acid[END_REF] The closest the norm to 0 (1) the less (more) probable the ionization.

While, the Slater determinant overlaps required in Eq. ( 24) can be readily computed in terms of the overlap matrix between atomic orbitals

| uv u v S  = 
 , a standard output when computing the electronic states, the CI coefficients in the framework of linear-response TDDFT requires some additional discussion, which is done in the next section.

Dyson orbitals with TDDFT

The theory presented so far in this section to compute the DOs is general and can in principle be applied for any method used to solve the electronic problem. The only condition is the representation of the electronic wavefunctions as a linear combination of Slater determinants.

Within the frame of Hartree-Fock based methods, that introduces no problem as it is a common assumption of the methodology. Therefore, the expansion coefficients n c and m c are directly computed. In the case of TDDFT, approximated wavefunctions in the CI form should be built.

According to the Casida's Ansatz for state assignment, [START_REF] Casida | Time-dependent density functional response theory for molecules[END_REF] the electronic wavefunction corresponding to a given excited state K can be represented as

, K K ov ov ov C  =   ( 26
)
where o and v stand for occupied and virtual spin-orbitals of the same spin, respectively. Denoting The use of Eq. ( 26) for building wavefunctions out of TDDFT amplitudes has become very popular recently. [START_REF] Barbatti | Surface Hopping Dynamics with DFT Excited States[END_REF] It has been extensively used for computations of nonadiabatic couplings in dynamics simulations [START_REF] Barbatti | Non-Adiabatic Dynamics of Pyrrole: Dependence of Deactivation Mechanisms on the Excitation Energy[END_REF] and also employed to compute different types of quantities, including spinorbit couplings, [START_REF] Franco De Carvalho | Derivation of spin-orbit couplings in collinear linear-response TDDFT: A rigorous formulation[END_REF] transition dipole moments, 10b, 32 nonadiabatic coupling vectors, [START_REF] Tavernelli | Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm-Dancoff approximation and beyond[END_REF][START_REF] Tavernelli | Nonadiabatic coupling vectors within linear response time-dependent density functional theory[END_REF] and Dyson orbitals. 8b In fact, this same methodology has been generalized [START_REF] Plasser | Surface Hopping Dynamics with Correlated Single-Reference Methods: 9H-Adenine as a Case Study[END_REF] to build wavefunctions to other linear-response-based methods as well, like ADC(2) and CC2. [START_REF] Christiansen | The 2nd-Order Approximate Coupled-Cluster Singles and Doubles Model CC2[END_REF] The expansion coefficients in Eq. ( 26) can be explicitly computed as [START_REF] Casida | Time-dependent density functional response theory for molecules[END_REF] 1/ 2 ( ),

K K K vo ov K ov ov K C A X Y E   - =+    (27)
where is a normalization factor introduced to ensure electronic wavefunctions normalized to unity.

K E  is the transition energy 0 , KK VV E  =- (28)

Steady PE spectra of imidazole and adenine

According to the developments of Section 2.1, calculation of the steady spectrum at a given value of kinetic energy Ek of the photoelectron can be pictured in three main steps: (i) generation of the nuclear ensemble; (ii) calculation of the DOs and IPs for each nuclear configuration and F electronic state considered; and (iii) calculation of the individual photoelectron intensities, from which the full spectrum is statistically computed. Along this section, we will illustrate each of these steps when applied to the He(I) photoionization of imidazole 20c and adenine. 20b

Nuclear ensemble

The steady PE spectra of imidazole and adenine were computed at E = 21.21 eV, which corresponds to the energy of a He(I) source. As the ionization is assumed to occur from the electronic and vibrational ground states, the nuclear configurations were sampled according to Eq. ( 6). An ensemble of Np = 500 nuclear geometries was generated in normal modes coordinates for each molecule and then transformed back to Cartesian coordinates. The equilibrium geometries R0, normal mode frequencies ωi, and normal mode eigenvectors were computed at B3LYP/augcc-pVTZ level. DFT and TDDFT calculations here, as in rest of the paper, were all performed with GAUSSIAN 09. 36

IPs and DOs

For each nuclear configuration, the electronic ground state of the neutral molecules was computed within DFT at CAM-B3LYP/aug-cc-pVDZ level. The first 40 states of imidazole cation and the first 10 of adenine cation were considered. All those doublet excited states were computed within TDDFT with the same functional [START_REF] Yanai | A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP)[END_REF] and basis set. [START_REF] Dunning | Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen[END_REF] After building approximated electronic wavefunctions for all these states, the DOs corresponding to each particular I = 0 (neutral) → F (cation) transition were computed according to the formalism presented in Section 2.4. To illustrate this step in more details, Table 1 (imidazole) and As can be seen from the tables, theoretical IPs are in excellent agreement with experimental ones. Another interesting feature that can be appreciated is that although all 0 → F transitions are energetically allowed ( 0F VE  ), not all of them contribute to the spectrum. For imidazole, for instance, only a few transitions really do so, the rest being practicably negligible, given their small DO norms. Moreover, among the significant transitions, we can find very intense ones (norms close to one), like the 0 → 0 in both molecules, and some much less intense, like the 0 → 9 for imidazole or 0 → 7 for adenine. Thus, with DO norms alone, one can not only identify which transitions really contribute to the spectrum, but also what the relative contribution from each transition will be.

F exp 0 t F V  ( ) 0 0 F V  R ( ) 0 0 d F  R F exp 0 t F V  ( ) 0 0 F V  R ( ) 0 0 d F  R 0 8.
TD-CAM-B3LYP OVGF a F exp 0 t F V  ( ) 0 0 F V  R ( ) 0 0 d F  R ( ) 0 0 F V  R P 1/

Steady PE spectrum

Once IPs and DOs (Eq. ( 22)) are computed for each nuclear configuration Rl of the ensemble (Eq. ( 6)) and for F electronic state of the cation, the spherically averaged total cross section for the same geometries ( ( )

IF l
 R , Eq. ( 20)) are computed, and the spectrum is simulated (Eq. ( 5)).

( )

IF l
 R are computed with the EZDYSON 3.2 program. The free-electron wavefunction was expanded in Coulomb partial waves to an angular momentum of lmax = 6, which we found out to be enough to converge the ionization probabilities. The photoelectron dipole matrix elements were averaged over all molecular orientations, which is justified by the non-polarized character of the laser used in the experiments and the random orientation of the molecules before the ionization.

Alternatively, we have also computed the spectrum based on the DO norm approach, using Eq. transition moments in the cross section approach tends to deliver better balanced relative intensities.

 = eV.
Concerning intensity, only the cross section approach can provide absolute values. In intensities for these molecules have not been experimentally reported, and the same normalization procedure was applied. The good agreement between the spectra computed with cross sections and DO norms in the low binding energy region implies that the cross section approximation in Eq. ( 21) is valid for IF EV   . Thus, as long as absolute intensities are not required, the low-energy region of the photoelectron spectrum may be simulated with the DO norm approach, significantly reducing computational costs.

The bottom graphs of Figure 3 and Figure 4 show the contribution of each S0→Dn transition to the total cross section. For both molecules, only few cation states contribute to the spectrum up to 14 eV. In some cases, a single experimental band may correspond to the overlap of transitions into more than one state, as for instance transitions into D1 and D2 forming the second photoelectron band of the two molecules. Above 14 eV, the number of states needed to compute the spectrum increases substantially.

For imidazole, for instance, transitions into nine states (D3 to D12) contribute to the broad band starting at 13 eV. For the next band starting at 18 eV, even considering 24 states (D15 to D39), we have not been able to reproduce the experimental band shape. This large demand for states in the high energy region points to a major limitation of the method. Not only the simulation costs may be prohibitive, but also the computed properties of such highly excited states are not fully reliable, especially within linear-response approximation.

Time-resolved PE spectra of imidazole

As in the steady case, the calculation of the time-resolved PE spectrum at a given Ek can also be pictured in three main steps: (i) generation of the nuclear ensemble, (ii) computation of the DOs and IPs, and (iii) calculation of photoelectron intensities, from which the spectrum is statistically computed. However, according to the developments of Section 2.2, step (i) is fundamentally different now: the nuclear ensemble has to be selected from nonadiabatic dynamics. Along this section, we illustrate all these steps, when applied to simulate the time-and kinetic-energyresolved PE spectrum of imidazole.

Nuclear ensemble from surface hopping

As imidazole is initially pumped from the vibrational ground state, the sampling of the nuclear coordinates and conjugate momenta at t = 0 was performed according to a Wigner distribution function for S0. The same R0, ωi, and normal-mode eigenvectors as in Section 3.1 were used. A set of 500 nuclear configurations and conjugate momenta were generated and projected onto the adiabatic electronic states to compute the absorption spectrum shown in Figure 5.

In the experimental setup, 8d a pump laser of energy E1 = 6.18 eV (200.8 nm) was used to directly excite imidazole from the electronic ground state into the 1 ππ* state. From the computational side, using TDDFT at CAM-B3LYP/aug-cc-pVDZ level, we found out that the simulated absorption spectrum is blue-shifted by 0.2 eV compared to the experimental spectrum 20a (Figure 5). Therefore, to excite the 1 ππ* state of imidazole in the same region as done in the experiments, an energy E1 = 6.4 eV is necessary in the computational modelling. This value was used in the simulations.

Figure 5. Simulated and experimental 20a absorption spectrum of imidazole in the gas phase.

The intensity of the experimental spectrum was normalized to match the simulated one.

To initiate the dynamics, 500 phase-space points sampled for the absorption spectrum were screened to select those with excitation energy 0 I EEwithin the narrow energy interval E1 ± 0.1 eV, and resampled using their corresponding oscillator strength as transition probability. 2b A set of 100 points matching these energy-window and oscillator-strength criteria was selected to be used as initial conditions for trajectories.

The number of initial conditions per adiabatic electronic state is shown in Table 3. This distribution reflects the geometric distortions in the sampling. Although the vertical excitation into the bright * state is S3 for the equilibrium geometry, Table 3 shows that, depending on the geometry, this state may shift as down as S2 and as high as S6. The electronic energies, energy gradients, and nonadiabatic coupling terms were computed 'on-the-fly' within the frame of TDDFT at CAM-B3LYP/aug-cc-pVDZ level. Excited electronic states up to S8 were included in the dynamics. Each trajectory was propagated for a maximum of t = 500 fs, with an integration step of 0.5 fs for classical equations and 0.025 fs for quantum equations.

Nonadiabatic transitions between different electronic states were treated with the fewestswitches surface hopping [START_REF] Tully | Molecular-Dynamics with Electronic-Transitions[END_REF] including decoherence corrections ( = 0.1 Hartree). [START_REF] Granucci | Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping[END_REF] Nonadiabatic couplings with TDDFT were computed by finite differences with the method discussed in ref. 30a , which is based on the Hammes-Schiffer/Tully approach. [START_REF] Hammes-Schiffer | Proton-Transfer in Solution -Molecular-Dynamics with Quantum Transitions[END_REF] As a single-reference method, TDDFT cannot provide reliable nonadiabatic couplings for crossings with the ground state. For this reason, when a trajectory reached an S1-S0 energy gap smaller than 0.15 eV before the maximum simulation time, it was stopped. This procedure did not affect the spectrum simulations, as at this point the probe energy was already smaller than the ionization energy.

After completing the dynamics, trajectories were split in regular intervals of  = 25 fs starting from  = 0. For each time interval i between 0 and the maximum simulation time, Np = 500 nuclear geometries () i l R were randomly selected from the trajectories and used to compute the spectrum.

Initial conditions, semi-classical dynamics, absorption spectrum, and photoionization spectrum were computed with NEWTON-X interfaced with GAUSSIAN 09.

IPs and DOs

According to ref. Before proceeding with the spectrum discussion, it is illustrative to characterize the ionization process for the S0 equilibrium geometry 0 R , which approximately corresponds to the ionization at  = 0. Excitation energies, IPs, and DO norms for the lowest I → F transitions are shown in Table 4. As can be seen, only electronic states of the cation up to F = 2 (D2) need to be considered for this particular geometry. (In fact, this is also true for all remaining geometries.)

The DOs corresponding to ionizations from S0 (closed shell) and S3 (*) are shown in Figure 6. The main configuration for each singlet and doublet state is schematically shown in the figure as well. The doublet configurations differ from that for S0 by a single spin orbital. Therefore, according to the ionization rules discussed in the Supporting Information (SI-2), such ionization processes are allowed. This is corroborated by the large DO norms for S0→Dn processes reported in Table 4. In the case of S3, only ionization into D0 is allowed according to the ionization rules, which is also corroborated by the results in Table 4. As we discuss this ionization pattern from S3 will have major consequences for the simulation of the time-resolved spectrum.

Time-resolved spectrum

Once IPs and DOs are known for each nuclear geometry () i l R between the current state I and F cation states, the time-resolved PE spectrum can be computed either with the cross section approach (Eq. ( 20)) or with the DO norm approach (Eq. ( 21)). Then, if we suppose that the electron is ejected with the maximum kinetic energy, the photoelectron spectrum is computed based on peaked line shapes (PVB model), as given by Eq. ( 7). In Figure 7-top, the simulated spectrum using the PVB model is shown for the early dynamics, collecting configurations generated in the first 25 fs of dynamics simulations. The experimental spectrum from ref. 8d for zero time delay is also shown, normalized to the maximum of the simulated result. (For an analysis of the experimental results, see ref. 8d, 20d, 20e ). The spectra are plotted in terms of the binding energy The large biding energy in the experimental data shown in Figure 7-top implies that electrons are being ejected with low kinetic energies, which has been attributed by Humeniuk et al. 8b to a rearrangement of the nuclear wavepacket due to its interaction with the probe pulse. This also means that the hypothesis underlying the PVB model, i.e. that electrons are ejected at the maximum kinetic energy, does not hold in the present case and simulations based on the CVB model may be more adequate. In Figure 7-bottom, we show the spectrum simulated with this model, as given by Eq. ( 9), for times smaller than 25 fs. As before, transition dipoles were computed with Coulomb partial waves (lmax = 6) between the current neutral state and all cation states up to F = 4. 10 li n = points were used in the linear interpolation of Eq. ( 11). The experimental spectrum from ref. 8d for zero delay is also shown, but now normalized to the intensity of the S4 → D0 contribution.

The agreement of the CVB model with the experiment is still not perfect, but it is significantly better than with the peaked model. The simulation correctly predicts a series of substructures in the spectral intensity distribution. As shown in Figure 7-bottom, the trace of the experimental spectrum at  = 0 exhibits inflexion changes at 8.7, 9.2, and 10 eV. The simulation shows equivalent inflection changes at 8.7, 9.0, and 9.7 eV. The data analysis revealed that they are related to which neutral states are contributing to the ionization. Below 8.7 eV, only ionization into D0 coming from S6 contributes to the spectrum. Above this value, D0 ionization of S4 to S5 also contributes causing the spectral shift. Starting from 9.0 eV, ionization of S2 and S3 into D0 starts to contribute. S1 ionization appears at above 9.7 eV. A small contribution from ionization of S4 and S5 into D1 is observed above 10 eV.

The early dynamics spectrum computed with the DO norm approach and the CVB model is also shown in Figure 7-bottom. It is normalized to match the cross-section based spectrum at the S4 → D0 contribution. The agreement between the two approaches is very good, once more indicating that the DO norm approach can be applied as an inexpensive alternative to the full the experimental spectrum peaks at 10.5-10.7 eV and quickly drops to zero before the maximum biding energy (11.33 eV), the simulations do not show this important feature, but only a flat plateau extending up to the maximum binding energy (this is better seen in Figure 7-bottom). The reason underlying this difference can be traced back to three factors. First, the CVB model itself, which, as already discussed, provides only a very approximated guess on how vibrational overlaps modulate the distribution of electron kinetic energies. Second, the experimental setup in which electrons with small kinetic energy are not fully collected, causing the intensity drop before the maximum binding energy. Third, a TDDFT failure to describe the multiconfigurational character of the * state of imidazole. This last point is discussed next.

We know from CASPT2 calculations for imidazole that the * state has strong contributions from 1* and * configurations (see, for instance, ref. [START_REF] Serrano-Andrés | Theoretical study of the electronic spectrum of imidazole[END_REF] ). This multiconfigurational character plays a central role for the spectrum, splitting the ionization signal in two components, depending on whether the hole is created in the 1 or  orbital. The 1 hole , which is formed after ejection of the electron in the 1* orbital (approximately the DO 30 d  in

Conclusions

We have implemented semiclassical methods based on the nuclear ensemble approach to simulate steady and time-resolved photoelectron spectra. The current implementation in the NEWTON-X program works with TDDFT provided by GAUSSIAN 09, but the methods are rather general and can be easily adapted to work with other electronic structure levels and programs.

For both, steady and time-resolved photoelectron spectra, we have developed and tested two levels of approximations, one based on full computation of transition dipole moments (cross

  . A N-electrons molecule in the electronic state I and the stationary vibrational level m, with eigenvalue m I E , is excited with a monochromatic laser of fixed energy E. As a result, an electron with asymptotic kinetic energy k are the mass and wavevector of the photoelectron) is ejected, leaving the remaining N-1 electron system in the electronic state F and vibrational level n, with eigenvalue n F E .

Figure 1 .

 1 Figure 1.Schematic representation of the photoionization for the steady case.

E

  R at time t = τ is built by first running a conventional surface hopping simulation,[START_REF] Barbatti | Nonadiabatic Dynamics with Trajectory Surface Hopping Method[END_REF] and then collecting geometries R within a time window t      +  after the photoexcitation. For each R, regularly spaced between zero and 2 IF EV - in the case of Eq. (9). For evaluation of  , we search this grid for the values immediately inferior ( ) compute  with the linearly-interpolated cross section

(

where 0 

 0 photoelectron transition dipole matrix element as a function of the nuclear coordinates R, formally defined as

  terms of Slater determinant overlaps | j mn    and configuration interaction (CI) coefficients n c and m c defining the electronic wavefunctions of the N and N-1 electron systems, respectively. To reduce computational costs, q nj d terms with expansion coefficient n c or m csmaller than an arbitrarily small value cis  can be neglected. In all results discussed here, we have adopted 0.01 cis  = .

by 0 

 0 the Kohn-Sham Hartree-Fock ground state determinant, ov  are singly excited Slater determinants, where the o th occupied spin-orbital has been replaced by the v th virtual one of the same spin, analogous to a configuration interaction with single excitations. Notice that only excited  =  .



  ) is the energy of the molecular orbital o (v), and KK + XY is the linear response TDDFT vector associated with the K th electronic state. The remaining term,

Figure 3 .

 3 Figure 3. Simulated (this work) and experimental 20c steady PE spectra of imidazole for a laser energy of 21.21 eV. The intensities of the experimental spectrum and of the simulation based on DO norms are normalized to match the maximum of the first peak. Bottom: Individual contributions from each cationic state up to D39; the dotted line shows the sum over all components. 0.1  = eV.

Figure 4 .

 4 Figure 4. Top: Simulated (this work) and experimental 20b steady PE spectra of adenine for a photon energy of 21.21 eV. The intensities of the experimental spectrum and of the simulation based on DO norms are normalized to match the maximum of the first peak. Bottom: Individual contributions from each cationic state up to D9. The dotted line shows the sum over all components. 0.1

Figure 3 and

 3 Figure 3 and Figure 4, the intensities of the spectra based on DO norms are normalized to match the intensity of the first peak of the spectra based on cross sections. Unfortunately, absolute

  Photoabsorption cross section (Mb)

  8d , a laser of E2 = 4.93 eV (251.6 nm) was used in the experiment to probe the dynamics after the E1 = 6.18 eV (200.8 nm) pump excitation. In the simulation, for each () i l R , the DOs and their norms associated with each I → F transition were computed, where I is the current electronic state of the neutral molecule at the moment of the ionization and F are all cation states from 0 to 4.

Figure 6 .

 6 Figure 6. DOs for ionization from S0 (closed shell) and S3 (*) into the first three cation states (D0 to D2) computed for the equilibrium geometry of imidazole with TDDFT.

  . Transition dipoles were computed with Coulomb partial waves(lmax = 6) between the current neutral state at a certain time step and all cation states up to F = 4 (D3). The line width was assumed to be

Figure 7 .

 7 Figure 7. Photoelectron spectrum for configurations sampled within the first 25 fs of the dynamics simulations. Top: Spectrum based on the PVB model. Bottom: Spectrum based on the CVB model. In both cases, experimental data from ref. 8d are shown. Different normalizations of experimental data are applied in each panel (see text).

  It is clear from Figure7-top that the spectrum computed with the PVB model poorly compares to the experimental result. The simulation has several peaks at the resonant points , with the main contribution coming from ionization of S2 and S3 into D0. The experimental spectrum, on its turn, is much broader and it peaks at much larger binding energies than predicted by the simulation.

Figure 6 ).

 6 Figure 6), leads to an ionization signal spanning from min

Table 2

 2

	(adenine) present the values of	0F V 	and 0 d F 	for the corresponding
	equilibrium geometries of each system, together with the experimental IPs for each molecule. 20b,
	20c

Table 1

 1

	.	IPs (in eV) and DO norms 0 d F 	corresponding to each 0 → F transition for the
	equilibrium geometry of imidazole. The experimental IPs reported in ref. 20c are given as well.

Table 3

 3

	.	Number of initial conditions per electronic state for which	0 EE -I	is within the
	energy interval E1 ± 0.1 eV.
		SI = S1 S2 S3 S4 S5 S6 S7 Total
		0 26 46 21 6 1 0 100

Table 4

 4

	.		Vertical excitation energies (singlet-singlet;	0I E 	in eV) and vertical IPs (singlet-
	doublet;	IF V 	in eV) computed at the equilibrium geometry 0 R . The values of the oscillator
	strength (f ) and of the squared DO norm (	2  R ) are shown in brackets. State assignments 0 () d IF
	in terms of the main orbital contribution are given in parenthesis. cs -closed shell; 3sX -3s
	Rydberg orbital on atom X.
				0I E 	(eV) [f ]	IF V  (eV) [	2  ] d IF
		I/F			0	1	2
					D0 (1 hole )	D1 (n hole )	D2 (2 hole )
		0 S0 (cs)		0.00 8.99 [0.96] 10.30 [0.91] 10.62 [0.83]
		1 S1 (3sN)	5.59 [0.000] 3.40 [0.47] 4.71 [0.00] 5.02 [0.00]
		2 S2 (3sC)	6.36 [0.030] 2.64 [0.47] 3.95 [0.00] 4.26 [0.00]
		3 S3 (*)	6.43 [0.169] 2.56 [0.42] 3.87 [0.00] 4.19 [0.01]
		4 S4 (3sC)	6.63 [0.000] 2.36 [0.38] 3.67 [0.00] 3.99 [0.00]
		5 S5 (n*)	6.72 [0.004] 2.28 [0.00] 3.59 [0.31] 3.90 [0.00]
		6 S6 (3sN)	6.99 [0.002] 2.00 [0.00] 3.31 [0.00] 3.62 [0.36]

Acknowledgments

MB thanks the support of the A*MIDEX grant (n° ANR-11-IDEX-0001-02) funded by the French Government « Investissements d'Avenir » program. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) also within the program « Investissements d'Avenir » supervised by the Agence Nationale de la Recherche. The authors thank C. Angeli and S. Ullrich for discussions.

the slope of the spectra for large binging energies. It is caused by a small effective near-linear dependence of the transition dipole on the electron kinetic energy, which is completely neglected in the DO norm approximation.

The time evolution of the spectrum using the CVB model with the cross section approach is shown in Figure 8. The sub-picosecond time distribution of the spectrum is well predicted. In particular, the simulation clearly reproduces the dependence of the time decay on the binding energy, with a systematic increase of the lifetime between 9 and 11 eV. In spite of the qualitative agreement between experiment and simulation with CVB model, the predictions for the binding energy distribution are, however, not entirely satisfactory. While section approach) and another based on an approximation of the transition dipole moments by Dyson orbital norms (DO norm approach). Moreover, the vibrational modulation of the electron kinetic energy distribution was also modeled with two different approaches. In the first one, vibrational overlaps between N and N-1 electron systems were supposed to be significant only for the electrons ejected with the maximum allowed kinetic energy (PVB model), a common approximation resting on the sudden-ionization hypothesis. In the second model, vibrational overlaps were supposed to be constant over the whole electron kinetic energy domain (CVB model).

Applications of the methods have been done for imidazole (steady and time-resolved) and for adenine (steady). The comparison to experimental data shows that steady spectra can be nicely predicted with the PVB model, with good description of intensities and band shapes.

For time-resolved spectra, the PVB model failed and the CVB model rendered significantly better results. The CVB simulations have been able to reproduce a series of substructures in the spectrum, which were assigned to specific ionization processes. Nevertheless, the overall agreement between simulation and experiment was less satisfactory than in the steady case due to the hypotheses underlying the CVB model, the multiconfigurational character of the key state contributing to the spectrum, and the instrumental signal loss not included in the simulations.

Considering all hypotheses and approximations invoked, it is truly encouraging that the main qualitative features of the spectrum have been predicted by the nuclear ensemble modelling.

For all tested cases, the approximation of the cross sections by Dyson orbital norms delivered results of similar quality as those based on full computation of cross sections, with an enormous economy of computational effort.

Finally, all these results make us confident that photoelectron spectrum simulations based on the nuclear ensemble approach can be an effective tool to aid deconvolution and assignment of experimental data for large molecules.

Author information

Corresponding author * E-mail: mario.barbatti@univ-amu.fr.