Vibronic dynamics from real-time time-dependent density-functional theory coupled to the Ehrenfest scheme: the example of p-coumaric acid

We investigate the vibronic dynamics of a modified version of the p-coumaric acid using real-time time-dependent density-functional theory coupled with the Ehrenfest scheme in the adiabatic local density approximation. Due to the issues of this functional to yield a reliable starting point for the evolution of the electron-nuclear system triggered by a pulse, we start off the simulations constraining the electronic occupation of the molecule in two excited states corresponding to a bright, delocalized transition, and a dark, charge-transfer-like excitation. By monitoring the kinetic energy spectral density, we analyze the nature of the nuclear motion over a time window of 300 fs. Anharmonic effects appear at low frequencies, below 500 cm-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}, and are particularly pronounced in the charge-transfer excitation. In this case, after about 200 fs, the molecular backbone becomes largely distorted and the initially constrained occupations evolve toward a different electronic configuration. On the other hand, the dynamics initialized from the delocalized bright excitation are electronically and structurally stable, and the resulting nuclear motion is markedly harmonic. Our results provide indications to decipher the vibronic dynamics of this chromophore and related systems in view of more elaborated simulations embedding the molecules in a realistic environment.


Introduction
Understanding the fundamental mechanisms of light-matter coupling in photoactive molecules is of paramount importance for their potential application in the field of optoelectronics and photonics.[1][2][3][4] From an experimental standpoint, state-of-the-art laser techniques are able to resolve the electronic and vibrational dynamics of such systems on the sub-femtosecond time scale.[5][6][7][8][9][10] Computational approaches that are suitable to describe these physics go beyond the Born-Oppenheimer approximation, [11][12][13] probing simultaneously the dynamics of the electronic and vibrational degrees of freedom.For example, the surfacehopping scheme, which is popular in quantum chemistry, has proven to be very accurate in predicting the dynamics of photoexcited molecules, including photo-isomerization, [14][15][16] conical intersections, [17][18][19][20] and chemical reactions.[21,22] However, the range of applicability of these methods is limited to systems composed of a few tens of atoms.
On the other hand, real-time time-dependent densityfunctional theory (RT-TDDFT) coupled with Ehrenfest molecular dynamics has demonstrated its efficiency in reproducing the electron-vibrational dynamics of isolated molecules and clusters, [23][24][25][26] as well as of crystalline solids [27][28][29] and hybrid interfaces.[30][31][32] Given the good parallelization of the involved algorithms, the costs of these calculations are affordable on the currently available supercomputing architectures even for systems with several hundreds of atoms.[33][34][35] In case of even larger systems, real-time time-dependent density functional tight-binding (RT-TDDFTB) [36] also in conjunction with Ehrenfest dynamics [37] offers a good alternative.For example, this method has proven to be suitable to describe charge-transfer excitations in hybrid nanoclusters [38,39] and to assess outof-equilibrium interactions of solute-solvent dynamics.[40] However, in contrast to RT-TDDFT which is fully ab initio, RT-TDDFTB requires some material parameters, and thus, 1 3 110 Page 2 of 12 its applicability is limited to those scenarios where such parameters are available.
In a simple molecule like ethylene, RT-TDDFT+Ehrenfest is able to describe the relaxation of the electron density upon the elongation of the C=C bond in the first excited state, [25] as predicted by the Franck-Condon model.However, it fails to reproduce the torsion exhibited by the photoexcited molecule.[41] The pitfalls of this methodology can occur even on a much lower level, namely, in the description of the linear absorption spectrum or even in the electronic structure of the system.The case of strongly dipolar molecules characterized by charge-transfer excitations at lowest energy is exemplary, [42] but unfortunately, even for purely conjugated compounds such as pyrene, the order of the excitations reproduced by RT-TDDFT calculations [26] does not coincide with the one predicted by quantum-chemical references.[43] It is evident that in these problematic scenarios, where already the static starting point is flawed, the description of the dynamics provided by RT-TDDFT+Ehrenfest is unreliable.
One way to overcome this issue is to computationally decouple the electronic and vibrational excitation process, for instance, by starting off the simulation of the vibrational dynamics with the system constrained in an electronically excited state.This approach was able to reproduce the photoisomerization path of azobenzene [44] and to decipher the physics of organic photovoltaic compounds.[45,46] With the implementation of external, time-dependent fields in the most popular packages for RT-TDDFT calculations, [23,[47][48][49] this method lost most of its appeal.Yet, in those scenarios where the underlying description of the electronic structure and linear optical response is even qualitatively inaccurate due to the above-mentioned issues of the common RT-TDDFT approximations, electronically constraining the system to mimic the target excited state remains a reliable approach to monitor the evolution of its electron-nuclear dynamics.
This is precisely the path followed in this work, where we explore the vibrational dynamics of an electronically excited molecule derived from the well-known p-coumaric acid (pCMA, see Fig. 1a), the chromophore of the photoactive yellow protein (PYP) that is present in cyanobateria.[50] This molecule has become popular due to its ultrafast tras-to-cis photoisomerization, [51] which appears particularly favorable for the design, for example, of bio-inspired photo-switches.[52,53] From a fundamental perspective, this system exhibits interesting characteristics.The presence of several atomic species (hydrogen, carbon, sulfur, and oxygen) in such a relatively small molecule gives rise to different types of chemical bonds coexisting in the same moiety that may follow specific vibrational dynamics upon different electronic excitations.In order to explore these effects, it is necessary to access the spectrum of excitations in the system.
As discussed below, the adiabatic local density approximation (ALDA) is unable to correctly reproduce it.For this reason, by constraining the initial electronic occupations, we prepare the molecule in two excited states, namely, a bright transition between the frontier orbitals and an optically forbidden charge-transfer-like excitation energetically close to the absorption onset.In these settings, we evolve the coupled electronic and nuclear degrees of freedom of the molecule and discuss the resulting vibronic dynamics in terms of vibrational kinetic energy spectral density and electronic population dynamics.The non-perturbative character of the adopted approach reveals both harmonic and anharmonic vibrational effects which appear with different magnitudes in different frequency ranges.

RT-TDDFT and Ehrenfest molecular dynamics
The electron-nuclear dynamics simulations are performed from RT-TDDFT with the time-dependent Kohn-Sham (TDKS) equations coupled with the Ehrenfest scheme.[54] At each iteration, the electronic density (r, t) is calculated by solving self-consistently the TDKS equations in real-space where R n (r, t) are the TDKS states at the fixed nuclear geometry R that can be expanded in the basis of static Kohn- Sham (KS) orbitals, m,0 (r|R) , at fixed nuclear geometry R , as with the normalization condition ∑ m �c R nm (t)� 2 = 1 ∀R and with M being the dimension of the basis that includes both the occupied and unoccupied electronic states.The electronic density at a fixed geometry is expressed as where the index n runs over the total number of electrons in the system, N e , and f R m (t) = ∑ N e n �c R nm (t)� 2 is the time-dependent occupation (or population) of the m-th static KS orbital, with the constraints that is the electronic coherence between m-th and m ′ -th KS states at the nuclear geometry R .In Eqs.( 1)-(3), we assume the dependence of the KS orbitals and of the electronic density on all time-dependent nuclear coordinates grouped in the 3N-dimensional vector In Eq. ( 1), the electron-nuclear interaction is described by the time-dependent external potential (1) while VHxc [(t)] is the time-dependent Hartree plus exchange-correlation (xc) potential in the adiabatic approximation that takes into account the electron-electron coupling and that functionally depends on the electronic density.
In the Ehrenfest scheme, [55] the electronic density is effectively coupled with the nuclear dynamics determined by the equation of motion where M i and R i (t) = R eq,i + ΔR i (t) are the mass and the spatial coordinate of i-th atom, respectively, ΔR i (t) is its displacement around the equilibrium position R eq,i , and is the time-dependent nuclear potential energy.The nuclear dynamics are analyzed by means of the nuclear kinetic energy spectral density (KESD) correlation function, [26] S( , h ) , calculated as where the angular bracket ⟨⟩ av indicates the temporal average 1 T ∫ T dt over a sufficiently long propagation time T; the KESD correlation function is S( , where V n ( ) is the Fourier transform of the normal velocity v n (t) = qn (t) obtained as the time derivative of the normal coordinate q n (t) .The normal coordinate associated with the normal mode e n is obtained by projecting the mass-scaled t i m e -d e p e n d e n t n u c l e a r d i s p l a c e m e n t ve c to r � onto the n-th normal mode as q n (t) = e n ⋅ Δ R(t) .The normal mode eigenvec- tors e n satisfy the eigenvalue equation ℍ − 2 n  e n = 0 , where n are the harmonic frequencies that enter the second argument of the KESD correlation function and are the matrix elements of the mass-scaled nuclear potential energy Hessian.The eigenvectors of ℍ are the column vectors of the (orthogonal) transformation matrix that diagonalizes it such that gives the eigenvalue of e n .In practice, the KESD gives infor- mation about how each normal mode contributes to the total (average) kinetic energy by decomposing it in its Fourier components.The degree of harmonicity (anharmonicity) of each normal mode n is quantified by the diagonal (off-diagonal) components The KESD can be conveniently visualized in a map where one

Dynamics of an electronically constrained molecule
In this work, we adopt an electronically constrained approach that has already been used in other contexts to quantify nonadiabatic couplings.[56] We prepare the system in an excited state evaluated at the ground-state equilibrium geometry R eq .From the resulting non-equilibrium configuration, the system is evolved in real-time by propagating the TDKS equations coupled with the Ehrenfest nuclear dynamics.The starting point is where the initial occupation coefficients, f m , which can be summarized in the vector f = f 1 ⋅ ⋅ ⋅ f M , are constrained such that the system is in an excited state that differs from the ground state configuration where H and L are abbreviations for the highest-occupied molecular orbital (HOMO) and lowest-unoccupied molecular orbital (LUMO), respectively.For example, when one electron is promoted from the HOMO to the LUMO, the occupation is Once the constrained occupations are chosen, they are kept fixed to recalculate the TDKS orbitals KS n,0 (r|R) at the current geometry R at each time step.The time-dependent occupation coefficients are calculated as where ⟨⟩ r means integration with respect to the electronic coordinates.The rate of change of the occupation is given by where ⟨ R m,0 �∇ R R m � ,0 ⟩ and Ṙ are the non-adiabatic couplings (NAC) and the nuclear velocity vectors, respectively.Equation (10) allows evolving in real-time the electronic ( 8) occupations by projecting each instantaneous TDKS orbital � n (t)⟩ onto the static one � m,0 ⟩ evaluated at the current nuclear geometry R .We should expect that the occupations do not change when the nuclear dynamics is adiabatic, i.e., the NAC are negligible and Ṙ ≃ 0 , in analogy with the Born-Oppenheimer nuclear dynamics.On the other hand, when there is non-adiabatic nuclear motion, the occupations will change over time with a rate controlled by the nuclear velocity vector and the NAC vectors.

Computational details
The calculations presented in this work are performed with the code Octopus [57] implementing a real-space numerical grid and pseudopotentials.A minimum simulation box size of interlocking spheres of radius 5 Å and grid resolution of 0.20 Å is adopted along with the set of Hartwigsen-Goedecker-Hutter norm-conserving pseudopotentials.
[58] The local density approximation (LDA) is adopted for the ground-state calculations and its adiabatic extension (ALDA) in the time-dependent runs.The molecular structure is optimized with a threshold for the interatomic forces equal to 10 −3 eV/ Å. Optical spectra are calculated in linear response (Casida scheme [59]) to access the composition of the individual excited states in terms of single-particle transitions.These results are compared against reference data computed with Gaussian 16 [60] using the hybrid functional CAM-B3LYP [61] and the cc-pVTZ basis set.The dynamics are calculated using RT-TDDFT+Ehrenfest employing the enforced time-reversal propagator scheme with a time step of 1.28 as and a total propagation time of 300 fs.To evaluate the population dynamics, the KS orbitals KS n,0 (r|R) are recalculated every 2 fs.Finally, to obtain the KESD, the vibrational dynamics are computed in the ground state geometry where the nuclear forces are minimized below 10 −3 eV/Å.
To give an idea of the complexity of these calculations, the RT-TDDFT+Ehrenfest results presented in this work for a molecule including 23 atoms and 68 valence electrons were obtained in about 13 h on 192 processors distributed over 2 computing nodes.Given the excellent parallelization available in Octopus, [33] we expect similar performance on larger systems upon a corresponding upscale of the employed computational resources.

Electronic structure and linear optical properties
We start our analysis by discussing the electronic and optical properties of the pCMA chromophore.In cyanobacteria, this system is attached to a protein via its sulfur atom and is negatively charged in the solvent environment.[51] Since in this work, we focus on an artificially isolated chromophore in vacuo, we slightly modify its structure by bonding a methyl group to the sulfur atom and by protonating the oxygen atom on the opposite side of the molecule to make it overall charge-neutral.The resulting geometry is shown in Fig. 1a.Due to the presence of different chemical species such as carboxyl, methyl, hydroxyl groups, carbon, and sulfur, the pCMA chromophore is highly asymmetric within the molecular plane and belongs to the C s point group, which contains only a mirror symmetry with respect to the horizontal plane.Furthermore, this molecule is electronically inhomogeneous, which is manifested in a static dipole moment.Thus, the LDA will unlikely deliver an accurate description of the electronic structure of the system.Rangeseparated hybrid functionals such as CAM-B3LYP [61] are currently considered the reference approximations in DFT for these types of molecules.[62] The energies of the frontier orbitals and of the gap between them calculated with both functionals are reported in Table 1.It is evident that the LDA leads to a substantial underestimation of the HOMO-LUMO gap (2.69 eV) compared to CAM-B3LYP (6.73 eV) due to the incorrect description of the HOMO and LUMO energies, which are too shallow and too deep, respectively.
With the knowledge gained so far, we move on to the inspection of the optical properties of pCMA obtained from ALDA (Fig. 2, black curve).In this approximation, the spectrum features two maxima at the onset, which are named S b1 and S b2 .As summarized in Table 2, S b1 is domi- nated by a transition between the HOMO and the LUMO accompanied by an additional one between the HOMO-1 and the LUMO.It is polarized along the long molecular axis, in accordance with the spatial distribution of the frontier orbitals (see Fig. 1b).The second bright excitation, S b2 , is also given by a mixture of the same transitions but with reverse weights with respect to their contributions to S b1 .Hence, S b2 , is polarized along the x-axis, too.In addition to the two above-mentioned optically active excitations, the ALDA spectrum features a dark transition at lowest energy corresponding to HOMO-2→LUMO and named S d .This forbidden excitation corresponds to a charge-transfer-like transition due to the localization of the HOMO-2 around the S atom and the C=O group in the moiety (see Fig. 1b) and the delocalized character of the LUMO.
Considering now the reference spectrum obtained adopting the CAM-B3LYP functional (see Fig. 2, red curve), we notice some remarkable differences with respect to the result from ALDA.First and foremost, the onset is characterized by a single absorption peak, S b1 , dominated by the HOMO→LUMO transition.The dark transition between the HOMO-2 and the LUMO ( S d ) is present also in the CAM-B3LYP spectrum but at slightly higher energy (10 meV) with respect to S b1 .Notice that in this result, the forbidden transition HOMO-2→LUMO represents again the dominant contribution to S d but it is also accompanied by the transition from the HOMO-2 to the LUMO+2 (not shown) which makes the oscillator strength very weak but non-vanishing, in contrast with the ALDA result.The second bright excitation in the CAM-B3LYP spectrum has a different character than its counterpart obtained from ALDA.It stems from the transitions from the HOMO to the LUMO+1 and from the HOMO-3 to the LUMO.Its oscillator strength is weak (see Fig. 2 and Table 2) and its polarization is split between the x and y directions in the adopted system of reference (see Fig. 1).
The analysis of the linear absorption spectrum of pCMA reveals the limitations of ALDA to capture even qualitatively the optical features of this molecule.While the characters of the lowest excited states, S b1 and S d , agree with the CAM-B3LYP result, their energetic order is swapped.Higher-energies excitations, starting from S b2 , diverge sub- stantially both in terms of relative energies with respect to the onset as well as in terms of composition.These findings are not surprising and are in agreement with earlier results obtained for related pCMA derivatives.[63] However, these discrepancies cast a shadow for the use of the ALDA in the study of the electron-vibrational dynamics of the considered molecule using the RT-TDDFT+Ehrenfest formalism in combination with an ultrafast time-dependent pulse to excite the system.[47] In fact, due to the finite width of the applied field and the non-equilibrium charge distribution obtained during the joint time evolution of the electronic and vibrational degrees of freedom, the presence of parasitic excitations can negatively affect the quality of the results.In such a scenario, it is recommendable to use an alternative method, in which the system is initially constrained to the excited state that one wants to target during the dynamics.

Vibrational dynamics from electronically excited state
In this section, we investigate the electron-nuclear dynamics of pCMA in the excited-state configurations corresponding to S b1 and S d and realized by promoting one electron from an initially occupied state to an unoccupied one, namely, from the HOMO to the LUMO in the former case and from the HOMO-2 to the LUMO in the latter.The dominant contributions of these two transitions in the excitations computed from CAM-B3LYP (see Table 2) support this choice.Note that, being optically forbidden, S d cannot be accessed directly via absorption from the ground state.However, it can participate in more complex (de)excitation processes, which makes it an interesting case to examine here.Having prepared the system in these excited electronic configurations, we monitor their electron-nuclear dynamics in terms of KESD maps.We split this analysis in two regimes, namely, in the region below 500 cm −1 (Sec.3.2.1)and the one between 500 and 3600 cm −1 (Sec.3.2.2).

Nuclear dynamics below 500 cm −1
In Fig. 3, we display the KESD maps in the low-frequency region between 0−500 cm −1 for the electron-nuclear dynamics calculated for the systems in excited state S b1 (panel a) and in S d (panel b).We notice that the kinetic energy associated with the nuclear motion from S d is one order of magnitude larger than the one from S b1 , as it is visible from the values assigned to the color scale.The vibrational modes mostly triggered in these dynamics are marked in the KESD maps and visualized in Fig. 3c.They correspond to translations [(i) and (ii)] and to a rotation of the phenyl ring combined with slow stretching modes on right side of the molecule [(iii) and (iv)].
The kinetic energy distributions visualized in Fig. 3a, b suggest an anharmonic response of the system in both excited electronic configurations.In fact, the maxima are not spread along the diagonal, as expected in the case of harmonic nuclear dynamics when the excited frequencies coincide with those of the normal modes.Instead, for each normal frequency, the KESD maxima are on the lower part of the frequency axis (y-axis) in Fig. 3a, b.This behavior is not surprising considering the normal modes excited in this dynamics (Fig. 3c), which correspond to collective motions of the molecular backbone or of specific groups.Notably, the largest relative contribution to the KESD from S b1 is associ- ated with the mode at 70 cm −1 involving the motion of the phenyl ring and the other C atoms in the molecule, where the Table 2 Energy, components of the transition dipole moment in the adopted reference system displayed in Fig. 1a, oscillator strength (OS), and composition (only contributions larger than 10% are included) of the first three excited states of pCMA computed from adiabatic LDA (calculations performed with Octopus) and CAM-B3LYP (calculations performed with Gaussian 16) transition density is mostly localized according to the spatial distribution of the HOMO and the LUMO (see Fig. 1b).On the other hand, in the dynamics from S d , the KESD map is dominated by the mode at 222 cm −1 , which corresponds to the distortions of the C=O bond together with the C-S-C angle, and to the rotation of the phenyl ring.Again, this result is unsurprising considering the charge-transfer character of S d , with the initial state (HOMO-2) localized on the side of the molecule opposite to the phenyl ring (see Fig. 1b).Quantitatively, the kinetic energy accumulated in the dynamics from S d is almost two orders of magnitude larger than the one from S b1 [compare the values of the color scales in Fig. 3a, b].This finding suggests that anharmonic effects, while present in both dynamics, are dominant below 500 cm −1 when the system starts off in the excited-state configuration associated with S d .A better understanding of this behavior can be grasped by looking at selected snapshots of the molecular geometry at subsequent stages of nuclear dynamics, see Fig. 4.While the structure of pCMA does not undergo significant changes during the whole dynamics (300 fs) initiated from S b1 , when the propagation starts off in the electronic configuration of S d , we notice that the system is subject to a rotation of its backbone with respect to the adopted system of reference and to a substantial deformation of the C-S-C bond in the last stage of the simulation.Considering the normal modes that mostly participate in the dynamics from S d (see Fig. 3c), the information about the structural changes in real time provided in Fig. 4b enables explaining the anharmonic nature of this vibronic motion.

3.2.2
Nuclear dynamics between 500 cm −1 and 3600 cm −1   In Fig. 5, we report the KESD maps calculated in the frequency region 500-3600 cm −1 for the dynamics evolved from the excited states S b1 (panel a) and S d (panel b).In this range, we notice a different behavior with respect to the one analyzed below 500 cm −1 .The nuclear dynamics are substantially harmonic as suggested by the diagonal distribution of the maxima in the KESD maps.The normal modes excited at these frequencies mainly correspond to carbon-carbon bond stretches.In the dynamics initiated from S b1 , where one electron is promoted from the HOMO to the LUMO, the dominant mode is the one with frequency 1190 cm −1 , see Fig. 5c.Other activated modes are in its energetic proximity [between 1100 and 1600 cm −1 , labeled as (i)-(v) in Fig. 5a and c] and correspond to C-C stretches, too.An additional weaker contribution to the KESD comes from higher energy modes with frequency 2966 cm −1 , which involve oscillations of the C-H bonds in the methyl group attached to the S atom as well as in the phenyl ring.
Moving now to the analysis of the KESD obtained from S d , one notices a larger number of relevant contributions.We wish to emphasize that in this frequency region, the magnitude of the kinetic energy contributions for the dynamics from S b1 and S d are of the same order of mag- nitude as indicated by the identical color bar adopted in Fig. 5a, b.Consistent with the charge-transfer nature of S d , some of the modes that are mostly activated in its dynamics pertain to the motion of one side of the molecule only.This is the case of the mode at 939 cm −1 , involving the S-C bond, as well as of the ones at 999 cm −1 , 1068 cm −1 , and 1151 cm −1 , corresponding to C-C stretching modes of the phenyl ring [64] coupled with C-H bending modes (see Fig. 5d).In addition, the modes participating in the dynamics from S b1 in the range 1190-1600 cm −1 are acti- vated also from S d .The two modes at 2966 cm −1 contrib- ute to the KESD map from S d , too.Notably, these modes induce weak anharmonicities as can be seen by carefully inspecting the map shown in Fig. 5c, with the signal at 2966 cm −1 extending to the lower part of the vertical frequency axis.This result can be better understood with the support of the real-time snapshots of the molecule displayed in Fig. 4b.The pronounced distortion of the C-S-C angle and of the methyl group bonded to the S atom in the latest stage of the dynamics impact the kinetic energy of the modes close to 3000 cm −1 , which leads to the slight anharmonicity visualized in Fig. 5c.Evidently, this is a minor effect compared to the prominent anharmonic dynamics of the molecule in the frequency region below 500 cm −1 .

Interplay between electronic and nuclear dynamics
The behavior of the nuclear dynamics started from S b1 and S d can be rationalized by looking at the evolution of the electronic population shown in Fig. 6.When the system is in the electronic configuration corresponding to S b1 , the HOMO and the LUMO occupations, both starting off at 1, remain essentially constant during the explored time window of 300 fs with no relevant population transfer onto other KS states.This is a signature of the dominantly adiabatic nature of the electronic dynamics where the nuclei  move harmonically.On the other hand, when looking at the population dynamics initiated from S d , where one electron is promoted from the HOMO-2 to the LUMO, a significant fluctuation of the occupations is seen in Fig. 6b.In this case, the population of the HOMO-2 increases to almost 2 within a few tens of fs and fluctuates between 1 and 2 until approximately 250 fs when it drops below 0.5.The occupation of the LUMO remains slightly below 1 for about 200 fs and then it drops close to zero.Other initially occupied and unoccupied states that are visible as shaded black and red lines in Fig. 6b, respectively, also behave regularly within the first 200 fs; afterward, they undergo significant drops and raises, respectively.
These behaviors indicate the presence of NACs that control the rate of change of occupations as expressed in Eq. (10).According to Eq. ( 3), these rapid variations cause a dynamical redistribution of the charge density that, in turn, influences the nuclear motion [see Eq. ( 5)].As opposed to S b1 , the non-adiabatic dynamics initiated from S d can be considered responsible for the presence of the low frequency anharmonic kinetic features below 200 cm −1 and for the additional harmonic features at about 1000 cm −1 (see Fig. 5b).In fact, the NACs enter into the electron-nuclear potential energy in Eq. ( 6) through the electronic density that can modify dynamically the local curvature and shape of the multidimensional nuclear potential energy landscape and thus the degree of harmonicity or anharmonicity of the nuclear dynamics.

Summary and conclusions
In summary, we investigated the excited-state nuclear dynamics of a modified version of the pCMA chromophore in vacuo adopting RT-TDDFT+Ehrenfest in the adiabatic local-density approximation.The evolution of the nuclear degrees of freedom is initiated with the system constrained in an excited state.We considered two configurations, one corresponding to the optically active HOMO→LUMO transition ( S b1 ) and the other one representing a dark excitation with partial charge-transfer character dominated by the HOMO-2→LUMO transition ( S d ).We analyzed the contributions of the nuclear modes by means of kinetic energy spectral density maps which also offer information about the (an)harmonic character of the motion.In the low-frequency region, below 500 cm −1 , collective modes corresponding to translations and rotations of the molecular backbone dominate the dynamics and lead to sizeable anharmonicities.This behavior is particularly relevant for the dynamics initiated from S d and can be related to the partial charge-transfer nature of this electronic excitation.Residual anharmonic effects can be seen in this regime even in the dynamics from S b1 but their magnitude is overall substantially lower than for S d .At higher frequencies, between 500 and 3600 cm −1 , the response of the system is essentially harmonic regardless of the initial electronic configuration.In this region, the kinetic energy spectral density is dominated by the C-C stretching modes and by other oscillations of the atoms along their bonds.A careful inspection of the modes participating in the dynamics from the two excited-state configurations reveals again differences related to the distribution of the electron density in the corresponding transitions.In the case of S b1 , where the electrons distribute homogeneously over the backbone, the C-C stretching mode prevail.In contrast, by starting off the dynamics from S d , non-negligible contributions to the kinetic energy spectral density come from the modes associated with the phenyl ring, in agreement with the charge-transfer character of this excited state.In this configuration, slight anharmonicities are visible also for the C-H bond oscillations around 3000 cm −1 and can be associated to the substantial structural rearrangement of the molecule during the dynamics.These signatures can be further interpreted in light of the analysis of the timedependent occupations of the electronic levels.While the dynamics from S b1 is stable, in the sense that the occupations of the HOMO and the LUMO remain essentially constant during the explored time window of 300 fs, starting from S d , a relevant redistribution of the electronic populations occurs around 200 fs.This effect is ascribed to non-adiabatic couplings which are in turn considered responsible for the substantial distortion of the molecules in the same time frame.
The results of this study are relevant from both a methodological and a physical point of view.Mimicking excited-state configurations by constraining the electronic occupations grants access to the dynamics of the system even in those cases, such as the one examined here, in which the underlying approximations do not guarantee a reliable starting point for the dynamics triggered by a pulse.Clearly, this constrained approach enables the simulation only of the resonant condition and is blind to the transient phase in which the excitation builds up.From a physical perspective, in spite of the simplicity of the model system adopted in this work (neutral moiety in vacuo, no effects from the environment, etc.), we gained an understanding of the intrinsic connection between the electronic configuration (energy level distribution and spatial extension of the charge density) and the nature of the nuclear motion.In the example considered here, the dynamics triggered from a delocalized excitation are much more stable and lead to a predominantly harmonic nuclear response compared to the case in which the evolution is initiated from an electronic configuration mimicking a charge-transfer excitation.We expect that the insight gained from this work will be useful to simulate the vibronic dynamics of chromophores in more realistic 110 Page 10 of 12 configurations, including, for example, the effects due to their environment.

Fig. 1 a
Fig. 1 a Ball-and-stick representation of the modified pCMA considered in this work; C atoms are depicted in black, O atoms in red, S atoms in yellow, and H atoms in white.b Kohn-Sham orbitals of pCMA calculated with LDA; results obtained from CAM-B3LYP are visually identical.The isosurfaces are fixed to 10% of their maximum value

Fig. 3 Fig. 4
Fig.3KESD maps of pCMA in the region between 0-500 cm −1 obtained for a the system initiated in the excited state S b1 where one electron is promoted from the HOMO to the LUMO, and b in the (dark) excited state S d where one electron is promoted from the HOMO-2 to the LUMO.The integrals of the KESD are shown on the side panels.The frequency of each normal mode is associated with a linewidth of 20 cm −1 along the horizontal axis.In panel c, the atomic displacements of the dominant normal modes highlighted in the KESD maps are reported

Fig. 5
Fig.5KESD maps of pCMA in the region between 500-3600 cm −1 obtained for the system initiated a in the excited state S b1 where one electron is promoted from the HOMO to the LUMO, and b in the (dark) excited state S d where one electron is promoted from the

Fig. 6
Fig. 6 Electronic population dynamics of the occupied (black) and unoccupied (red) time-independent Kohn-Sham orbitals when the system is evolved a from S b1 with the initial occupations f HOMO = 1 , f LUMO = 1 and b from S d with initial occupations f HOMO−2 = 1 , f LUMO = 1 .The shaded lines indicate the populations of the other orbitals in the range (HOMO-20)-HOMO (gray) and LUMO-(LUMO+20) (pink) Page 4 of 12axis corresponds to the 3N discrete harmonic frequencies n associated with 3N − 6 vibrational, 3 rotational and 3 trans- lational modes, while the other axis is the actual frequency associated with the kinetic energy spectrum carried by n-th normal mode |V n ( )| 2 .

Table 1
Energies