A hybrid classical-quantum algorithm to simulate ECD spectra – the case of tryptophan zwitterions in water

Abstract

The vibronic structure of electronic circular dichroism (ECD) spectra is simulated in silico using a hybrid method that encompasses ab initio time-dependent density functional theory (TD-DFT) calculations and the Doktorov quantum algorithm. Four low-energy electronic excitations of the chiral tryptophan amino acid are characterized and their spectra are simulated with and without considering the Duschinsky rotation effects. The normal mode mixing shows a very small impact on the vibronic structure. Moreover, it demonstrates the versatility of the quantum algorithm to focus on individual normal mode contributions to the vibronic structure of the spectrum.

---

1. Introduction

The last decades have witnessed the development of quantum chemistry methods and algorithms to predict molecular properties, ranging from vibrational and NMR signatures to linear and nonlinear optical responses.^1–6 More recently, developments have focused on elaborating quantum computing algorithms.^7,8 These encompass quantum algorithms for predicting the molecular ground state energies,^9–18 the vibrational and vibronic signatures of molecules,^19–32 accounting for solvent effects,³³ their excited states and excitation energies,^34–41 and nonlinear optical properties.⁴² Recently, the vibronic structures of UV-vis absorption spectra have also been calculated using quantum algorithms.⁴³

The current contribution deals with the vibronic structure of electronic circular dichroism (ECD) spectra. ECD is one of the typical optical activity (OA) phenomena, routinely used for determining the absolute configuration of chiral molecules, and their enantiomeric purity.^44,45 Enantiomers of the same molecule show the same degree of OA but with opposite signs (opposite signs of the ECD peaks) so that when a solution contains an enantiomeric excess of a chiral molecule, the enantiomers of that molecule can be distinguished using ECD. While the interpretation of ECD spectra is commonly performed using time-dependent density functional theory (TD-DFT) within the Franck–Condon (FC) approximation,^46–49 this approach often neglects vibronic coupling effects which can significantly influence spectral features.^50–55 More advanced methods including Herzberg–Teller coupling and explicit vibronic simulations are necessary for capturing fine structure and intensity modulation in chiroptical spectra.^52,56

To address this, the present study aims to exploit the promise of quantum computers for simulating ECD spectra with their vibronic structure. Quantum computing offers an efficient encoding and manipulation of vibronic transitions through unitary transformations. The Doktorov formalism⁵⁷ can be implemented as a sequence of Gaussian operations,^20,23,26 or as in our recent works, as a trotterized quantum circuit,^32,43 enabling the calculation of FC profiles and other vibronic observables with potential quantum speedup. Given the limitations of near-term quantum hardware, a hybrid classical–quantum approach is employed, where classical algorithms provide the electronic and structural information, which then serve as input for the Doktorov-based quantum algorithm.

As a case study, the zwitterionic form of the chiral tryptophan molecule is chosen (Fig. 1). It is an essential amino acid of biological relevance, notable for its chiral centers and well-characterized zwitterionic forms in solution. The Rectus (R-) and Sinister (S-) enantiomers exhibit distinct ECD signals, making them suitable benchmarks for validating quantum simulations of optical activity.

Fig. 1 Sinister- and Rectus-tryptophan molecules in their zwitterion form. (top) The 2D line-dash-wedge representation. (bottom) The 3D stick-ball representation of the ground state geometry.

The proof of principle of this methodology has already been given for UV-vis absorption spectra.⁴³ In this work, the methodology is adapted to retrieve ECD spectra while studying the effects of the Duschinsky transformation. This article is organized as follows: The Theory section presents the ECD absorption equation, the Doktorov and circuit formalism for the vibronic framework and the main structural–physical factors. The Methods section presents the level of theory adopted for the TD-DFT calculations, the methodology used to obtain the vibronic structure, and the data processing complement. The Results and discussion section presents the simulated ECD spectra of tryptophan, comparing cases with and without the inclusion of Duschinsky effects. Finally, the Conclusions section highlights the key findings and insights of the study.

2. Theory

The rotatory strength R₁₂ of a vibronic transition between an initial state 1 and a final state 2, relevant for ECD spectroscopy, under the Born–Oppenheimer approximation, is expressed as the dot product of two vectors multiplied by the Franck–Condon (FC) factor:The FC factor represents the square of the vibrational overlap between the initial and final nuclear wavefunctions, |ψ^vib₁〉 and |ψ^vib₂〉, respectively. The second part corresponds to the electronic chiroptical factor, specific to ECD given in its length-form, involving the imaginary component of the scalar product between the electric transition dipole moment and the magnetic transition dipole moment vectors. This term captures the differential absorption of left- and right-circularly polarized light and is characteristic of the chiroptical activity.

The FC factors can be seen as a distribution of the probability of a given electronic transition into the vibrational states. This energy conservation is outlined by the following property:

where the sum is over all possible vibronic transitions ν of the considered state transition. One major complexity when estimating the FC integrals lies in the coordinate system used to express its components. In practice, the vibrational wavefunction is written in terms of the molecule's normal mode (NM) coordinates, since these provide a natural basis for describing its vibrational motion. In order to calculate the overlap, a transformation is therefore needed to change the initial state (1) NM coordinate basis into the final state (2) NM coordinate basis (or vice versa). In first quantization calculations, this relation between PES (and their coordinates) is usually handled via the Duschinsky transformation.⁵⁸ In the boson ladder operator formalism, another strategy is typically applied: one of the bosonic wavefunctions is transformed from the initial (ground) state into the final (excited) state via the Doktorov's operator:The primary framework employed in this work to obtain vibronic intensities associated with electronic transitions is based on the second-quantized formulation of Doktorov's expressions.⁵⁷ As detailed in our previous study,³² the Doktorov operator Û_Dok is implemented through a trotterized Doktorov ansatz within a (superconducting) quantum circuit simulator. This quantum circuit representation, combined with the SWAP test circuit^59,60 and supplemented by classical post-processing, enables the evaluation of vibronic intensities corresponding to each bosonic transition.

This theoretical framework is valid within the harmonic oscillator approximation. In turn, the Doktorov's tensor is the product (i) of a rotational operator Û_r that mixes the NM coordinate bases of the initial and final states, (ii) of the squeezing operator Û_s, which converts the vibrational frequencies of the initial state into those of the final state, and (iii) of the translation operator Û_t that accounts for the displacements of the initial and final states PES in the direction of the NM coordinates.

The rotational operator is as big as the dimension spanned by the considered normal modes to mix. Conversely, the squeezing and translation operator can be decoupled into operators acting on each NM independently:³²

where ∀ l∈ [1,L],The operator Û_t,l denotes the translational operator associated with a specific normal mode l. Its definition involves the final (excited-state) vibrational harmonic frequency (or pulsation) ω_2,l, the mass-weighted displacement d_l, and the bosonic creation and annihilation operators, and b̂_l, corresponding to the lth mode. Similarly, the operator Û_s,l represents the individual squeezing transformation for a given mode. It depends on both the initial (ground-state) and final (excited-state) vibrational frequencies, denoted ω_1,l and ω_2,l, respectively, as well as on the corresponding bosonic ladder operators. To facilitate the parameterization, it is useful to introduce dimensionless coefficients that capture the essential physical characteristics of the transformation:These coefficients are unique for a specific normal mode l and also a given electronic transition. The rotational operator admits a more complex form due to its multidimensional nature. However, due to hardware limitations (more about it in the Methods section), only a pair of normal modes are coupled together herein when accounting for the Duschinsky effects. Therefore, two-dimensional rotation operators are defined:where θ_l1,l2 is the angle of rotation of the considered normal modes. It is found by applying the natural logarithm function to the 2 × 2 -sub-Duschinsky matrix:where U_D is the Duschinsky matrix, and U_D,ij is its i-line- j-row element. The previous equality is valid provided that the sub-Duschinsky matrix is a unitary rotational matrix. Complementary information on the Doktorov operators and their implementation can be found in our previous study.³²

The Huang–Rhys factor, HR_l, is a useful quantity to describe the displacement along the l normal coordinate. It measures the ratio between the vibrational relaxation energy and the quantum of vibrational energy. It is connected to the translational displacement coefficient via the relation:

This parabolic formulation highlights its role in quantifying the squared displacement between the minima of the excited-state potential energy surface (PES) and the vertical (adiabatic) transition point from the ground state. Any displacement along a given vibrational NM can be expressed (and, therefore, calculated) in terms of the atomic Cartesian shifts arising from the geometric relaxations occurring during the electronic transition:⁶¹where ΔX_Aζ denotes the Cartesian shift of atom A along the ζ ∈ {x,y,z} direction, m_A is the atomic mass of atom A, and Q_Aζ,l corresponds to the elements of the orthonormal matrix that maps the vibrational eigenmodes onto the normal displacements. In other words, it is the matrix that diagonalizes the mass-weighted Hessian. A second way to estimate the d_l and HR_l factors is to use the Gradient method, which consists of estimating the mass-weighted displacement by considering the relation:⁶²where E^elec₂ is the electronic energy of the final state that is differentiated with respect to the normal mode coordinate q_l, and the subscript 0 refers to the equilibrium geometry of the ground state.

With the Doktorov operator fully characterized, both bosonic vibrational states can be represented within the initial NM coordinate frame. The final form used to compute a given Franck–Condon (FC) factor is:

The remainder of this study is dedicated to implementing this formalism on a superconducting quantum processor simulator to predict the vibronic features of ECD spectra, with a specific demonstration for the case of tryptophan.

3. Methods

3.1. Classical algorithms: TD-DFT calculations

The molecular geometries were constructed using the Drawmol package.⁶³ A multi-step protocol was adopted to perform the quantum chemical calculations, as outlined below:

(I) Ground-state geometry optimization was carried out employing the CAM-B3LYP^64,65 exchange–correlation (XC) functional, which belongs to the range-separated hybrid family. This XC functional incorporates 19% Hartree–Fock (HF) exchange at short range and 65% at long range, using a Coulomb-attenuating scheme with a range-separation parameter of μ = 0.33 Bohr⁻¹. The chosen atomic basis set was the triple-ζ, polarized, 6-311G(d) Pople basis. The self-consistent field (SCF) and geometry-optimization convergence criteria were set to the 'tight' thresholds implemented in Gaussian 16,⁶⁶ and the DFT calculations employed the UltraFine integration grid.

(II) Using the same level of theory as in step I, the Hessian matrix, vibrational frequencies, and normal modes were computed. This step provides the vibrational characteristics of the ground state.

(III) The first 30 vertical electronic excitation energies were calculated at the TD-DFT level, along with their corresponding oscillator and rotatory strengths.

(IV) Steps I and II were repeated at the TD-DFT level for each excited electronic state of interest. This included full geometry optimizations and vibrational frequency analyses, thereby giving the vibrational data for the final (excited) states.

This final step is crucial to account for the vibrational frequency shifts that accompany electronic excitation. Only transitions with sufficiently large rotatory strengths (those expected to contribute visibly) were considered for this step, i.e. here the 3rd, 5th, 6th, and 7th excited states.

For all steps, incorporating the solvent effects was achieved by using the integral equation formalism polarizable continuum model (IEF-PCM).⁶⁷ The solvent used was water, modeled with a static dielectric constant ε = 78.355 and a dynamic dielectric constant ε_dyn = 1.78, in accordance with reference values for water at T = 298 K,⁶⁸ as implemented in Gaussian 16.⁶⁶

As mentioned in the Introduction section, several quantum algorithms have been proposed to compute the electronic and vibrational properties. However, current quantum hardware is not yet capable of performing such calculations without errors, and the classical simulation of these quantum algorithms is not presently feasible because of the limited memory resources. Consequently, all vibronic parameters required in this work were computed using classical algorithms.

Once all relevant data were obtained, a custom Python script was developed to import the previously computed ground-state Hessian matrix. The program then calculates the geometric displacements between the ground and excited states, determining the Huang–Rhys factors viaeqn (9), (10), or eqn (11). Obtaining the Duschinsky matrix was achieved by applying the Santoro et al. algorithm^56,69 implemented in Gaussian,⁶⁶ which requires the information retrieved in steps II and IV.

3.2. Quantum algorithms: vibronic structures

Quantum hardware simulations were performed using the state-vector backend of Qiskit⁷⁰ and Tangelo,⁷¹ which ensures an idealized, error-free simulation. The construction and application of the hybrid quantum algorithm were extensively detailed in our prior work.³² As such, the following summary focuses on the essential components, while preserving the necessary rigor.

For any given electronic transition, and within the harmonic approximation, the nuclear wavefunction is expressed as a tensor product over (vibrational) normal mode wavefunctions:

where it is assumed that the system is at 0 K, so the initial state is always taken to be the vibrational ground state. Preparing the initial state in higher vibrational levels would be equivalent to introducing finite-temperature effects. The wavefunctions define the Fock basis for each bosonic mode. Since these bases are infinite-dimensional, a truncation is required. The under-sizing is governed by the number of qubits M allocated per mode in the chosen bosonic encoding.^25,28 In this work, the compact Standard Binary (SB) encoding^25,72 is employed. Within the SB encoding, initializing a specific Fock state is straightforward: a X gate is applied to each qubit corresponding to a 1 in the binary representation of the desired vibrational level. For example, with M = 3 qubits per normal mode (supporting up to 8 vibrational states), state |5〉 is encoded as |101〉, and initialized via X₂X₀|000〉. This preparation procedure is simply referred to as initialization. Given the considerable size of the molecule with its L = 75 normal modes, the full tensor product is intractable with current quantum simulation capabilities. The available memory resources restrict to consider at most two vibrational modes simultaneously. Therefore, this study focuses on computing spectra for individual modes as well as for pairs of coupled modes. This will give an estimate of the impact of Duschinsky rotations. These partial (sub-)spectra are subsequently combined to reconstruct the full vibronic response, as described further below.

To represent the Doktorov operator within the framework of quantum computing, the bosonic ladder operators are expressed in matrix form, truncated to a dimension of 2^M × 2^M. Since these operators are defined in the Fock basis, their ket-bra decomposition is straightforward and can be mapped onto qubit operators via the following bosonic encoding:

The three components of the Doktorov operator (rotation, squeezing, and translation) are then represented as complex exponentials of linear combinations of Pauli strings. As in previous works, these exponentials are decomposed using a high-order Trotterization scheme^73–77 up to the 64th order. The Trotter decomposition is implemented using a qubit-wise strategy, resulting in a product of exponentials, each involving a single Pauli string. The full decomposition procedure and examples are provided in the SI of our previous work³² (SI).

The next step involves preparing the quantum circuit to compute the Franck–Condon (FC) factors, as defined in eqn (12). Two distinct scenarios must be considered: independent normal modes and coupled pairs of normal modes. For independent modes, only translation and squeezing operations are required, as no vibrational mode mixing occurs. In contrast, coupled modes necessitate the application of the full Doktorov operator, which includes the additional Duschinsky rotation operator.

The retained method for evaluating the overlap between vibrational states is the SWAP circuit.^59,60 This method involves a sequence of CNOT gates applied between corresponding qubits of the two wavefunctions, followed by a partial layer of Hadamard gates to perform measurements in the Bell basis. For the case of independent modes, the simulation uses 8 qubits: 4 qubits to encode the ground vibrational state and 4 for the target excited state |n_l1〉. For coupled modes, a total of 16 qubits are used: 8 for the double ground state |0,0〉 and 8 for the (double) excited state |n_l1,n_l2〉. The SWAP based approach allows for an easy choice of the input wavefunction. The ansatz is depicted in Fig. 2.

Fig. 2 (Top) Specifications for the individual and the coupled normal mode circuits. (a) Doktorov circuit representation. The three operations are applied for a pair of normal modes, whereas only the two operations inside the dotted rectangle are applied for a single independent mode. (b) The complete ansatz with the Doktorov circuit and the SWAP method.

To summarize the sequence of operations: the circuit is initialized in the single (or double) ground vibrational state |0〉 using (2) M qubits. For a selected mode (or pair of modes), l₁ (and l₂) ∈ [1,L], a vibrational state |n_l1(n_l2)〉 is initialized. Each mode is associated with a specific Huang–Rhys factor HR_l1 (and HR_l2), which is transformed into a displacement parameter d_l1 (and d_l2), giving the Doktorov parameters α_l1, β_l1 (and α_l2, β_l2, and θ_l1,l2). These parameters are then used to specify the trotterized Doktorov operator. The operator transforms the initialized excited vibrational state into its translated and squeezed (and, for coupled modes, also rotated) actual expression.

Finally, after running the Doktorov-SWAP circuit and after statistically measuring over 8 (or 16) qubits, the output is processed via a dot product with the post-processing vector . The resulting scalar value corresponds to the simulated FC factor for the given vibrational state, normal mode(s), and electronic transition. This procedure is repeated for all vibrational levels n_l1 (or n_l1,n_l2) ranging from 0 to 2^M − 1, yielding the individual (or coupled) normal mode spectra. The entire process is then reiterated across all independent modes l₁ ∈ [1,L], and all relevant mode pairs (l₁,l₂) to construct the full vibronic profile.

3.3. Simulating the UV/vis absorption and ECD spectra

Once all individual normal mode spectra were obtained, a custom parallelized Fortran script was employed to combine them. One can notice that directly mixing the 16⁷⁵∼10⁹⁰ sub-FC factors is computationally infeasible with current technology. To circumvent this limitation, three pre- and post-screening strategies were implemented.

First, a lower threshold was applied to the HR factors. When the HR factor is sufficiently small, the nuclear displacement is minimal, implying that the ground and excited PESs are nearly identical (along that coordinate). In such cases, the ground-state vibrational wavefunction overlaps almost entirely with its excited-state counterpart, and thus |〈0|0′〉|² ≈ 1. The accuracy of this assumption depends on the amplitude of the squeezing (the difference in vibrational frequencies). A threshold value of HR_min = 10⁻⁶ was chosen.

Second, a minimum FC factor threshold was imposed based on previous works.⁴⁷ During the computation of FC factor products, if the contribution of an individual factor falls below FCmin = 10⁻⁷, or if the product of multiple FC factors drops below this threshold, the calculation is halted for that specific vibrational configuration. This avoids spending resources on negligible contributions.

Third, a normal mode level/class⁵⁶ restriction was applied to control the mode couplings. This pre-screening strategy limits the number of vibrational modes that can be simultaneously excited. For example, level 0 includes only the ground vibrational state |0,…,0〉; level 1 consists of states where only one mode is excited, such as |1,0,…,0〉, |2,0,…,0〉, …, up to |0,…,0,2^M − 1〉. More generally, level p includes all vibrational wavefunctions in which exactly p modes have non-zero occupation numbers. In this work, level 4 is chosen, as it is the largest level that could be calculated with the available computer resources.

Together, these three strategies drastically reduce the computational effort required to assemble the full vibrational spectrum from the individual mode contributions while their accuracy is monitored by assessing how well eqn (2) is satisfied. To simulate the final vibronic spectrum, the vibrational spectrum is scaled by its corresponding rotatory strength and shifted by the vertical electronic excitation energy. To account for phenomenological broadening, a Gaussian broadening was applied to each FC factor using a specified full width at half maximum (FWHM).

4. Results and discussion

The first step consists of selecting the targeted ECD transitions. Then, as a function of the retrieved HR factors and Duschinsky matrices, the respective vibronic spectra are simulated with or without accounting for rotational coupling.

4.1. Electronic excited states and ECD intensities

The results of the TD-DFT calculations on the tryptophan molecule are summarized in Table S1 of the SI. The resulting ECD spectra without vibronic structure are depicted in Fig. 3. The TD-DFT results show the expected well-symmetric ECD spectra: the Rectus and the Sinister are the respective mirror images of the other. The two first electronic transitions are separated by more than 1 eV from the others. Then, a more compact spectrum is observed where the remaining 28 transitions are clustered in a 3 eV range. The experimental curve extracted and adapted from Liu et al. work⁷⁸ is depicted for comparison. It shows a strong negative-to-positive progression that spans more than 1 eV and is localized in the 6–7 eV excitation energy range. Based on these results, we focused on the S₃ and S₇ electronic transitions, which exhibit large rotatory strengths (see Table S1 in the SI). The intermediate transitions were also considered, since their proximity in energy influences the overall spectral shape and enables the inclusion of both positive and negative rotatory strengths. In contrast, the S₀ → S₁, S₀ → S₂, and S₀ → S₄ transitions were not considered as their ECD intensities are low. Additional calculations using the 6-311+G(d) and the 6-311+G(d,p) basis sets were performed to assess their effects (see Fig. S1 in the SI). Although the inclusion of diffuse functions alters the electronic structure, the overall trends and progression in the region of interest remain similar. Then, given the limited resolution of the available experimental spectra and considering the illustrative purpose of this spectrum, there is no compelling reason to prefer one basis set over the others, so the 6-311G(d) basis set is retained. For completeness, the ECD spectra using the velocity-form representation were also depicted (Fig. S2 in the SI), showing negligible differences with respect to the length-form representation.

Fig. 3 TD-DFT/CAM-B3LYP/6-311G(d) ECD spectra as obtained by accounting for solvent (water) effects at the IEF-PCM level without vibronic structure. Both, the R molecule (red) and the S (blue) results are sketched. The vertical red lines are the R rotatory strengths of the different electro-magnetic transitions. The structure followed a Gaussian broadening with a FWHM = 0.2 eV. The ECD experimental spectrum in black extracted from Liu et al. work,⁷⁸ blue-shifted by 0.67 eV and scaled for better visualization and comparison.

4.2. Vibrational properties and HR factors

The vibrational frequencies and associated HR factors are listed in Table S2 in the SI for the targeted 4 transitions. With the exception of the 5th transition, all the electronic transitions have HR factors larger than 1 (∃l ∈ [1,L], HR_l > 1, red values in Table S2 in the SI). This can be problematic as several high HR values imply that the vibronic structure is very spread out and it impedes the calculation to converge (eqn (2) to be satisfied). However, the large HR factors are associated with NM with frequencies lower than 100 cm⁻¹. It must be noted that the frequencies under 200 cm⁻¹ are often inaccurately estimated,⁵ and furthermore, accounting for the low frequency modes tends to locally broaden the vibronic spectra. Due to these reasons, it seems reasonable to assume that neglecting their contributions would not affect the general shape of the spectra. For the rest of this work, the normal modes associated with vibrational frequencies below 200 cm⁻¹ of the S₀ → S₃ transition and below 100 cm⁻¹ for the other ones are neglected.

4.3. Vibronic structures of the S₀ → S₃ electronic transition

With the HR factors, the ground state- and the excited state-vibrational frequencies in handy, the Doktorov operator can be prepared to apply the translational and the squeezing operators.

4.3.1. Independent normal modes. In this sense, the reference spectrum to calculate is the vibronic structure without NM mixing. This is further called the independent normal mode (INM) spectrum. Therefore, the Duschinsky matrix is considered to be the identity, so no quantum circuit for the rotation is applied. The resulting INM spectrum is depicted in Fig. 4 for the S₀ → S₃ transition. Note that the chosen FWHM is considerably smaller to better visualize the vibronic progression.

Fig. 4 Vibronic structure of the S₀ → S₃ transition. (Top) Sums of FC factors, numbers of coupled NM pairs, and three largest HR factors. (Bottom) Simulated ECD absorption spectra with and without rotational coupling. The INM spectrum considers no Duschinsky effects (red). The Group I spectrum admits three pairs of rotationally-coupled NM (blue), while Group I & II considers two additional pairs of rotationally-coupled NM (green dashed line). The structure followed a Gaussian broadening with a FWHM = 0.05 eV.

The vibronic structure is very wide, spanning over 1.2 eV, so that it overlaps with other (purely) electronic transitions. This originates from the number of normal modes with HR factors ≥0.1, which implies that more vibrational excitations with smaller intensities have to be considered. Consequently, for the same coupling level, the sum of the FC factor is compromised, explaining that only about 70% of the vibronic total intensity was retrieved. As a check, our results were compared with those obtained with the classical FCellipse algorithm,⁴⁷ which also retrieves the vibronic structure without normal-mode coupling. Both approaches resulted in virtually identical spectra, further supporting the reliability of the INM implementation.

4.3.2. Pairs of coupled normal modes. The ground and excited state structures are used to retrieve the rotational matrix. The resulting Duschinsky matrix for the 3rd transition is depicted in Fig. 5. The rotational matrix is overall diagonal. Despite that, in principle, the whole Duschinsky matrix should be considered within the rotational operator, memory resources exponentially grow with the number of qubits used. With the current memory resources, only pairs of normal modes can be considered for rotational coupling. The couples of normal modes were selected on the basis of the amplitude of their couplings, while ensuring that the chosen (sub-)Duschinsky matrices remained close to unitary in order to avoid imaginary coefficients. The chosen pairs of normal modes for Duschinsky mixing are listed in the table in Fig. 5. This table also proposes a classification of the (sub-)rotational matrices: group I admits strong coupling values which either leave similar the normal mode intensities (strong diagonal terms, the case of the (35,36) and (69,70) tuples), or a almost complete exchange of vibronic intensity (the (64,65) tuple); group II are strong mixing of the two normal modes, giving a rather intensity sharing; group III are the last tuples with non-unitary matrices, where the mixing is not trivial. The idea is to estimate the impact on the ECD spectra of the successive groups. The reminder of the normal modes are considered as individual (or independent) normal modes with no rotational coupling.

Fig. 5 (Left) Pairs of normal modes considered for rotational coupling for the S₀ → S₃ transition; (Right) Absolute Duschinsky matrix for the S₀ → S₃ transition.

Yet, since the θ_3,4 and θ_46,47 rotational angle parameters are (found to be) imaginary, the corresponding couplings were dismissed because the current method is only valid for real parameters. This issue can be avoided by considering larger sub-Duschinsky matrices, relating more NM, and thus having a unitary transformation with real parameters. However, as previously mentioned, the available memory resources impede the implementation. On the other hand, Group I and II possess real angle parameters, and the results are depicted in Fig. 4.

The Group I spectrum shape is very similar to the non-rotated one. However, the absorption intensity has shifted to lower energies. It is difficult to assign a clear reason to this difference due to the dual effects to consider: the NM couplings redistribute the total intensity, changing the contributions, and thus, different peaks become more or less relevant, affecting the intensity. For instance, regarding the sum of FC, Group I retrieved approximately ∼0.016 more intensity than the uncoupled one while requiring slightly less FC factors. Nevertheless, upon normalization (which is tantamount to considering an equal intensity recovery), both spectra are almost identical and just the higher energies part of the spectrum is affected (see Fig. S3 in the SI). When including the remaining two extra pairs of coupled NMs, no remarkable change is noticed as the Group I spectrum is identical to the Group I & II spectrum.

4.4. The next dominant low-energy transitions

The Duschinsky matrices and the chosen coupled pairs for the S₀ → S₅, S₀ → S₆, and S₀ → S₇ transitions are depicted in Fig. S4, S6, and S8 in the SI, respectively. In summary, 5 rotational mixing pairs were chosen for the 5th transition, while 8 for the 6th and 7th transitions. They mostly correspond to strong couplings (Group I). Moreover, when comparing the INM and coupled spectra, they are found to be virtually identical (see Fig. S5, S7, and S9 in the SI). This is further appreciated in Fig. 6 where the 4 transitions have been merged. The two vibronic spectra are nearly identical. Their difference lies around the 6.4 eV energy, where the S₀→ S₃ transition is also contributing to the total spectra. Hence the difference is due to the previously discussed results in Fig. 4. The four spectra contributions to the total spectrum can be visualized in Fig. S10 in the SI. The almost exact overlap of the spectra over the remaining energies is easily explained when looking at the quantitative data of the individual results. The other transitions retrieve almost the same vibronic intensity when considering independent or (the considered) coupled NM.

Fig. 6 (Top) Summary table of the retrieved sum of FC factors, the number of considered coupled NM pairs, and the largest HR factor associated with the transition. (Bottom) Simulated ECD absorption spectra for the 4 transitions with and without rotational coupling. The INM spectrum considers no Duschinsky effects (red). The Coupled spectra considers all the coupled pairs of normal modes (black dashed line, cfr. Fig. S3). The structure followed a Gaussian broadening with FWHM = 0.05 eV.

5. Conclusions

This study has adopted a hybrid quantum-classical algorithm method to simulate the ECD absorption spectra of the tryptophan zwitterion molecule in water solvent. TD-DFT calculations have been performed for retrieving the necessary electronic and structural properties, which are subsequently used as parameters for the Doktorov quantum algorithm. After studying the structural deformations, the S₀ → S₃, S₅, S₆, and S₇ transitions were chosen for illustrating the calculations of the vibronic structure. Their Duschinsky matrices were examined, and three groups of rotational coupling pairs were defined. The Doktorov-SWAP ansatz was then employed to compute the Franck–Condon factors, and thus, the vibronic structure. It was applied to either independent or coupled pairs of vibrational normal modes, enabling the comparison between various levels of approximation. The effects of the different couplings show small impacts on the vibronic structure and the retrieved spectra. Further analyses of the vibronic structures, especially the rotationally coupled ones, are required to shed light on the spectral variations. The method is inherently limited by the power and memory computational resources, as well as by the absence of a theoretical extension capable of treating Herzberg–Teller contributions within the Doktorov second-quantization framework, which appears as a challenging area of investigation. Finally, this original work demonstrates the versatility of the method, allowing for controlled and in-depth examinations.

Author contributions

A. Soldera and B. Champagne planned and conducted the project. R. Olarte Hernandez developed the simulation programs and performed the numerical calculations. All authors discussed the results and wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: complementary Tables and Figures on the electronic properties, the vibrational frequencies, important structural parameters, and individual and merged considered vibronic structures. Electronic results calculated using TD-DFT calculations, table of Normal modes, vibrational frequencies and Huang–Rhys factors for the four considered electronic transitions. Comparisons between results obtained with different basis sets and with different representation of the dipole moment operator. Normalized spectra for the 3rd transition, and Duschinsky matrices and vibronic spectra of the 5th, 6th and 7th electronic transitions. See DOI: https://doi.org/10.1039/d5cp03805g.

Acknowledgements

R. O. H. thanks the Fonds Spécial de Recherche of the University of Namur and the University of Sherbrooke for his PhD grant. The authors thank Freddy Zutterman for fruitful scientific discussions and Vincent Liégeois for the implementation, in the DrawMol package, of the velocity-form representation in the simulation of the ECD spectra. This research used resources of the “Plateforme Technologique de Calcul Intensif (PTCI)” (https://www.ptci.unamur.be) located at the University of Namur, Belgium, which is supported by the FNRS-FRFC, the Walloon Region, and the University of Namur (Conventions No. U.G006.15, U.G018.19, U.G011.22, RW1610468, RW/GEQ2016, RW2110213, and RW2210148). All calculations were performed using the Hercules cluster. The PTCI is a member of the “Consortium des Équipements de Calcul Intensif (CÉCI)” (https://www.ceci-hpc.be). Likewise, this research was made possible in part by the support provided by Calcul Québec and the Alliance de recherche numérique du Canada (alliancecan.ca/en/). Armand Soldera gratefully acknowledges the Natural Sciences and Engineering Research Council (NSERC) of Canada (Grant No. RGPIN-2018-06126). We thank the Quantum AlgoLab at Institut Quantique for providing quantum computing resources and support.

References

1. R. G. Parr, Horizons of Quantum Chemistry: Proceedings of the Third International Congress of Quantum Chemistry Held at Kyoto, Japan, 1979, vol. 1989, pp. 5–15.
2. P. Sałek, O. Vahtras, T. Helgaker and H. Ågren, J. Chem. Phys., 2002, 117, 9630–9645.
3. T. Helgaker, P. Jorgensen and J. Olsen, Molecular electronic-structure theory, John Wiley & Sons, 2013.
4. C. J. Cramer, Essentials of computational chemistry: theories and models, John Wiley & Sons, 2013.
5. J. Bloino, A. Baiardi and M. Biczysko, Int. J. Quantum Chem., 2016, 116, 1543–1574.
6. S. Löffelsender, P. Beaujean and M. de Wergifosse, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2024, 14, e1695.
7. B. Bauer, S. Bravyi, M. Motta and G. K.-L. Chan, Chem. Rev., 2020, 120, 12685–12717.
8. M. Motta and J. E. Rice, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1580.
9. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704–1707.
10. S. Bravyi, J. M. Gambetta, A. Mezzacapo and K. Temme, arXiv, 2017, preprint, arXiv:1701.08213 DOI:10.48550/arXiv:1701.08213.
11. H. R. Grimsley, S. E. Economou, E. Barnes and N. J. Mayhall, Nat. Commun., 2019, 10, 3007.
12. S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin and X. Yuan, APS, 2020, 92, 015003.
13. K. Setia, R. Chen, J. E. Rice, A. Mezzacapo, M. Pistoia and J. D. Whitfield, J. Chem. Theory Comput., 2020, 16, 6091–6097.
14. H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes and S. E. Economou, APS, 2021, 2, 020310.
15. C. Cao, J. Hu, W. Zhang, X. Xu, D. Chen, F. Yu, J. Li, H.-S. Hu, D. Lv and M.-H. Yung, Phys. Rev. A, 2022, 105, 062452.
16. M. A. Jones, H. J. Vallury, C. D. Hill and L. C. L. Hollenberg, Sci. Rep., 2022, 12, 8985.
17. A. Eddins, M. Motta, T. P. Gujarati, S. Bravyi, A. Mezzacapo, C. Hadfield and S. Sheldon, PRX Quantum, 2022, 3, 010309.
18. F. Pavosevic, I. Tavernelli and A. Rubio, J. Phys. Chem. Lett., 2023, 14, 7876–7882.
19. S. Joshi, A. Shukla, H. Katiyar, A. Hazra and T. S. Mahesh, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 90, 022303.
20. J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean and A. Aspuru-Guzik, Nat. Photonics, 2015, 9, 615–620.
21. J. Huh and M.-H. Yung, Sci. Rep., 2017, 7, 7462.
22. Y. Shen, Y. Lu, K. Zhang, J. Zhang, S. Zhang, J. Huh and K. Kim, Chem. Sci., 2018, 9, 836–840.
23. N. P. D. Sawaya and J. Huh, J. Phys. Chem. Lett., 2019, 10, 3586–3591.
24. C. Wang, J. Curtis, B. Lester, Y. Zhang, Y. Gao, J. Freeze, V. Batista, P. Vaccaro, I. Chuang and L. Frunzio, et al. , BAPS, 2020, 65, 1.
25. S. McArdle, A. Mayorov, X. Shan, S. Benjamin and X. Yuan, Chem. Sci., 2019, 10, 5725–5735.
26. C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. Girvin and R. J. Schoelkopf, Phys. Rev. X, 2020, 10, 021060.
27. P. J. Ollitrault, A. Baiardi, M. Reiher and I. Tavernelli, Chem. Sci., 2020, 11, 6842–6855.
28. N. P. D. Sawaya, T. Menke, T. H. Kyaw, S. Johri, A. Aspuru-Guzik and G. G. Guerreschi, npj Quantum Inf., 2020, 6, 49.
29. N. P. D. Sawaya, F. Paesani and D. P. Tabor, Phys. Rev. A, 2021, 104, 062419.
30. M. Majland, P. Ettenhuber, N. T. Zinner and O. Christiansen, J. Chem. Phys, 2024, 160, 154109.
31. S. Bao, N. Raymond, T. Zeng and M. Nooijen, J. Chem. Theory Comput., 2024, 20, 5882–5900.
32. R. Olarte Hernandez, B. Champagne and A. Soldera, J. Phys. Chem. A, 2024, 128, 4369–4377.
33. E. R. Kjellgren, P. Reinholdt, A. Fitzpatrick, W. N. Talarico, P. W. K. Jensen, S. P. A. Sauer, S. Coriani, S. Knecht and J. Kongsted, J. Chem. Phys., 2024, 160, 124114.
34. P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta and I. Tavernelli, APS, 2020, 2, 043140.
35. Y. Ibe, Y. O. Nakagawa, N. Earnest, T. Yamamoto, K. Mitarai, Q. Gao and T. Kobayashi, Phys. Rev. Res., 2022, 4, 013173.
36. A. Asthana, A. Kumar, V. Abraham, H. Grimsley, Y. Zhang, L. Cincio, S. Tretiak, P. A. Dub, S. E. Economou, E. Barnes and N. J. Mayhall, Chem. Sci., 2023, 14, 2405–2418.
37. A. Kumar, A. Asthana, V. Abraham, T. D. Crawford, N. J. Mayhall, Y. Zhang, L. Cincio, S. Tretiak and P. A. Dub, J. Chem. Theory Comput., 2023, 19, 9136–9150.
38. P. W. K. Jensen, E. R. Kjellgren, P. Reinholdt, K. M. Ziems, S. Coriani, J. Kongsted and S. P. A. Sauer, J. Chem. Theory Comput., 2024, 20, 3613–3625.
39. Q.-X. Xie and Y. Zhao, arXiv, 2024, preprint, arXiv.2406.00296 DOI:10.48550/arXiv.2406.00296.
40. Y. Zheng, Z. Sun, J. Liu, Y. Fan, Z. Li and J. Yang, J. Chem. Theory Comput., 2024, 20, 9032–9040.
41. C. Cianci, L. F. Santos and V. S. Batista, J. Chem. Theory Comput., 2024, 20, 8940–8947.
42. M. Bruschi, F. Gallina and B. Fresch, J. Phys. Chem. Lett., 2024, 15, 1484–1492.
43. R. Olarte Hernandez, A. Soldera and B. Champagne, J. Phys. Chem. A, 2025, 129, 396–406.
44. C. Bannwarth and S. Grimme, Comput. Theor. Chem., 2014, 1040, 45–53.
45. M. Monti, M. Stener and M. Aschi, J. Comput. Chem., 2022, 43, 2023–2036.
46. E. Botek and B. Champagne, J. Chem. Phys., 2007, 127, 204101.
47. F. Zutterman, V. Liégeois and B. Champagne, ChemPhotoChem, 2017, 1, 281–296.
48. H. Jang, N. J. Kim and J. Heo, Comput. Theor. Chem., 2018, 1125, 63–68.
49. P. Garca-Cerezo, M. D. Codesal, A. H. David, L. Le Bras, S. Abid, X. Li, D. Miguel, M. Kazem-Rostami, B. Champagne and A. G. Campaña, et al. , Adv. Mater., 2025, 2417326.
50. P. L. Polavarapu, Chirality, 2002, 14, 768–781.
51. F. Egidi, V. Barone, J. Bloino and C. Cappelli, J. Chem. Theory Comput., 2012, 8, 585–597.
52. A. Baiardi, J. Bloino and V. Barone, J. Chem. Theory Comput., 2013, 9, 4097–4115.
53. P. L. Polavarapu and C. L. Covington, Chirality, 2014, 26, 539–552.
54. P. L. Polavarapu and E. Santoro, Nat. Prod. Rep., 2020, 37, 1661–1699.
55. J. L. Johnson and P. L. Polavarapu, J. Phys. Chem. A, 2021, 125, 8000–8013.
56. F. Santoro, A. Lami, R. Improta, J. Bloino and V. Barone, J. Chem. Phys., 2008, 128, 224311.
57. E. Doktorov, I. Malkin and V. Manko, J. Mol. Spectrosc., 1977, 64, 302–326.
58. H. Kupka and P. H. Cribb, J. Chem. Phys., 1986, 85, 1303–1315.
59. H. Buhrman, R. Cleve, J. Watrous and R. de Wolf, Phys. Rev. Lett., 2001, 87, 167902.
60. L. Cincio, Y. Subas, A. T. Sornborger and P. J. Coles, New J. Phys., 2018, 20, 113022.
61. H. Hameka, Physica, 1970, 46, 271–278.
62. J. Jortner, J. Chem. Phys., 1976, 64, 4860–4867.
63. V. Liégeois, DrawMol, UNamur, 2024, https://www.unamur.be/drawmol, Accessed 2024-01-01.
64. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
65. M. J. G. Peach, P. Benfield, T. Helgaker and D. J. Tozer, J. Chem. Phys., 2008, 128, 044118.
66. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, 2016, Gaussian Inc., Wallingford CT.
67. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
68. D. P. Fernández, Y. Mulev, A. R. H. Goodwin and J. M. H. L. Sengers, J. Phys. Chem. Ref. Data, 1995, 24, 33–70.
69. F. Santoro, R. Improta, A. Lami, J. Bloino and V. Barone, J. Chem. Phys., 2007, 126, 084509.
70. A. Javadi-Abhari, M. Treinish, K. Krsulich, C. J. Wood, J. Lishman, J. Gacon, S. Martiel, P. D. Nation, L. S. Bishop, A. W. Cross, B. R. Johnson and J. M. Gambetta, arXiv, 2024, preprint, arXiv:2405.08810 DOI:10.48550/arXiv:2405.08810.
71. V. Senicourt, J. Brown, A. Fleury, R. Day, E. Lloyd, M. P. Coons, K. Bieniasz, L. Huntington, A. J. Garza, S. Matsuura, R. Plesch, T. Yamazaki and A. Zaribafiyan, Tangelo: An Open-source Python Package for End-to-end Chemistry Workflows on Quantum Computers, 2022.
72. L. Veis, J. Visnak, H. Nishizawa, H. Nakai and J. Pittner, Int. J. Quantum Chem, 2016, 116, 1328–1336.
73. H. F. Trotter, Proc. Am. Math. Soc., 1959, 10, 545–551.
74. J. E. Cohen, S. Friedland, T. Kato and F. P. Kelly, Linear Algebra Appl., 1982, 45, 55–95.
75. L.-A. Wu, M. S. Byrd and D. A. Lidar, Phys. Rev. Lett., 2002, 89, 057904.
76. Y. Sun, J.-Y. Zhang, M. S. Byrd and L.-A. Wu, arXiv, 2018, preprint, arXiv:1805.11568 DOI:10.48550/arXiv:1805.11568.
77. A. A. Avtandilyan and W. V. Pogosov, arXiv, 2024, preprint, arXiv.2405.01131 DOI:10.48550/arXiv.2405.01131.
78. M. Liu, L. Chen, T. Tian, Z. Zhang and X. Li, Anal. Chem., 2019, 91, 13803–13809.

---