Delayed photoisomerisation of the trans-PSB3 retinal toy model using on-the-fly quantum dynamics

Abstract

We explore the trans–cis photoisomerisation process in a representative retinal protonated Schiff base known as trans-PSB3, employing the quantum dynamics method direct dynamics variational multiconfigurational Gaussian -DD-vMCG- in full dimensionality, i.e., 36 degrees of freedom on potential energy surfaces computed on-the-fly using the SA(2)-CAS(6,6)SCF electronic structure method with the 6-31G basis set. Although the toy molecule has been extensively studied using trajectory methods such as Tully surface hopping and ab initio multiple spawning, our application of the on-the-fly quantum dynamics method DD-vMCG reveals a previously unreported trans–cis isomerisation pathway through the S₁ electronic state minimum that occurs hundreds of femtoseconds more slowly, despite using the same electronic structure method. This slower timescale and alternative deactivation route can be explained in terms of the accesibility to the conical intersections connecting the ground and the excited state.

---

Sandra Gómez

Dr Sandra Gómez is currently an Assistant Professor in the Chemistry Department at the Universidad Autónoma de Madrid, a position she has held since 2024. She received her Dr rer. nat. in Chemistry from the Universität Wien in 2019, where she carried out her doctoral research under the supervision of Prof. Leticia González. Following her PhD, she undertook postdoctoral research with Prof. Graham A. Worth at University College London (2019–2022), before securing an independent research fellowship at the University of Salamanca (2022–2024). Dr Gómez's research centres on nonadiabatic dynamics in photochemistry and photophysics, with applications spanning femtochemistry, optoelectronics, and bioimaging. She is an active member of the JACS Au Early Career Advisory Board and contributes to the organisation of the Virtual Winter School for Computational Chemistry, which marks its twelfth edition in 2025.

---

1 Introduction

Vision in humans begins when light enters the eye and reaches the retina, where photoreceptor cells (rods for low-light vision and cones for color vision) convert light into neural signals interpreted by the brain.^1,2 In rods, the primary photoreceptor is rhodopsin, a complex of the opsin protein and the 11-cis-retinal chromophore, linked by a protonated Schiff base (PSB).³ Upon photon absorption, the retinal chromophore undergoes a rapid photoisomerisation from the 11-cis to the all-trans configuration,^4,5 triggering a cascade of biochemical reactions that culminate in optic nerve stimulation and vision.

This cis–trans photoisomerisation occurs on an ultrafast timescale (within 200 fs) and with high quantum efficiency,^6,7 proceeding through a conical intersection between the first excited (ππ*) and ground electronic states.⁸ The mechanism has been extensively investigated through both experimental and computational studies,^9–12 focusing on the excited state dynamics of the chromophores and their role in vision.

Given the complexity of the full rhodopsin system, simplified models such as the penta-2,4-dieniminium cation (PSB3) have become invaluable for probing the fundamental photochemistry of retinal.^13–15 PSB3 preserves the essential features of the retinal protonated Schiff base, including the central double bond and isomerisation pathway, while its small size enables high-level quantum chemical calculations of excited state potential energy surfaces.¹⁶

Nonadiabatic molecular dynamics (NAMD) simulations are crucial for studying photoinduced processes, but the diversity of available methods presents challenges in accuracy and reproducibility. The need for reliable benchmarks has led to the widespread adoption of PSB3 as a model system for testing and validating NAMD approaches.¹⁷ Its well-characterized photochemistry and reduced complexity have facilitated the development of a plethora of computational results, supporting method refinement and reproducibility, and helping theoretical advances to be extended to more complex systems.

Early computational studies mapped the potential energy surfaces (PESs) of PSB3 and identified the conical intersections (CIs) involved in its cis–trans isomerisation, using high-level methods such as CASPT2 and MRPT2.^18–22 Further work compared quantum chemical methods (e.g., CASPT2, CASPT3, ADC(2)), highlighting the influence of methodological choices on predicted excitation energies and CI branching spaces.^23,24 These studies highlighted the computational challenges of accurately modeling excited state processes, even in minimal systems.

On-the-fly nonadiabatic dynamics, employing techniques such as ab initio multiple spawning (AIMS) and surface hopping, have further explored the role of dynamic nuclear effects and nonadiabatic coupling in PSB3 photoisomerisation.^25–27 These investigations revealed that relevant conical intersections often lie far from minimum energy regions, with vibrational modes driving the isomerisation. Additional approaches, including exact factorization methods^9,28,29 and machine learning-enhanced molecular dynamics,³⁰ have also provided insight into the mechanisms and global reaction coordinates of retinal models.

Despite lacking some structural features of the full retinal system, PSB3 (C₅H₆NH₂⁺, see Fig. 1) is one of the simplest minimal models for 11-cis-retinylidene,³¹ capturing key aspects of the isomerisation mechanism and enabling qualitative understanding of retinal photochemistry.^9,32,33 Its suitability for ab initio calculations involving the ground and first excited electronic states makes it an ideal benchmark system.

Fig. 1 The PSB3 retinal model system: 2-trans-penta-2,4-dieniminium cation.

In this work, we employ for the first time the direct dynamics variational multiconfigurational Gaussian (DD-vMCG) method to study the photoisomerisation of trans-PSB3. This approach will allow us to accurately capture nuclear dynamics and probe fundamental aspects of rhodopsin photochemistry, as well as to compare the performance of this relatively new nonadiabatic dynamics method in predicting the isomerisation.

2 Computational details

2.1 Electronic structure

The electronic structure calculations were performed using a two state average complete active space self-consistent field theory (SA(2)-CASSCF) with an active space consisting of six electrons in six orbitals, denoted as CAS(6,6). This active space, shown in Fig. 2, includes three pairs of π and π* orbitals, capturing the key electronic configurations responsible for the photoisomerisationprocess in PSB3. The state specific approach SS-CASSCF with the same active space was used to optimise the ground state minimum and calculate its vibrational frequencies (geometric parameters and frequencies shown in Tables S1 and S5 in the SI). The state-averaged approach (SA(2)-CASSCF(6,6)) was applied over the two lowest singlet electronic states, S₀ (ground state) and S₁ (excited state). In terms of molecular orbitals, S₀ and S₁ states are characterized by π² and ππ* electronic configurations, respectively. SA(3)-CASSCF(6,6) calculations were attempted to also include the effects of the second excited state (S₂), leading to convergence problems. As dynamics studies reveal that the S₂ state is not involved in the population transfer,²⁵ only two-state average calculations were finally performed. The 6-31G basis set was selected for its ability to significantly reduce computational cost while still accurately reproducing vertical excitation energies, conical intersection locations and energies, and reaction paths when compared with polarized basis sets.^9,32,34

Fig. 2 Orbitals included in the active space for the SA(2)-CAS(6,6)SCF/6-31G electronic structure calculations. The S₁, with ππ* character, corresponds to an electron excited from orbital π₃ to .

Previous works, such as those by Martinez et al.²⁵ Szymczak et al.³² and Filatov et al.,²⁶ have employed the same CASSCF active space and the 6-31G basis set for studying nonadiabatic dynamics in model retinal chromophores. In the study by Martínez’ group, a three-state averaging procedure was adopted,²⁵ whereas in our work, we focus on a two-state average, emphasizing the dynamics between the ground S₀ and the first excited state S₁. Szymczak et al.³² also employed SA(2)-CASSCF(6,6) for the initial reference wavefunctions in their study, followed by multireference configuration interaction (MRCI) calculations. Since the goal of this work is to use, for the first time, the direct-dynamics variational multi-configurational Gaussian (DD-vMCG) method in full dimensionality to study the excited state photoisomerization of this system, we chose the CASSCF(6,6)/6-31G level of theory as our starting point. This choice was guided not only by its computational feasibility, but also because it matches the level of theory previously used by other authors in on-the-fly trajectory-based dynamics studies of the same system. This ensures that any differences we observe can be attributed solely to the choice of nonadiabatic dynamics (NAMD) method, rather than differences in the underlying electronic structure calculations.

Table 1 shows the ab initio optimised and vertical excitation energies (in eV) of both electronic states (S₀ and S₁) for the trans, the cis-PSB3 isomer and the minimum energy conical intersection(MECI_opt) using a SA(2)-CAS(6,6)SCF/6-31G method with the Molpro electronic structure program.³⁵ MECI_db corresponds to the PSB3 geometry with the lowest energy gap encountered in the dynamics.

Table 1 Ab initio optimised and vertical excitation energies (in eV) of both electronic states (S₀ and S₁) for the most relevant geometries marked in light pink in Fig. 4: the trans, the cis PSB3 isomer and the MECI_opt point. Calculations are done with a SA(2)-CAS(6,6)SCF/6-31G method using the Molpro program.³⁵ Also the MECI_db, which is the most accesible point in the database with the lowest energy gap encountered in the dynamics, is given

	π^2/S0	ππ*/S1
	Optimised	Vertical	Optimised
Trans-PSB3	0.00	4.84	4.45
Cis-PSB3	0.13	4.69	4.53
MECI_opt (S0/S1)	2.70	—	2.70
MECI_db (S0/S1)	2.61	—	2.62

---

2.2 The vMCG quantum dynamics method

In the variational multi-configurational Gaussian (vMCG) method, the nuclear wavefunction is expressed as a linear combination of time-dependent Gaussian basis functions (GBFs), avoiding the need for precomputing potential energy surfaces on a full primitive grid. This approach belongs to the family of Gaussian wave packet (GWP) methods, where the system's degrees of freedom are described by multidimensional frozen Gaussians, making it well-suited for nonadiabatic quantum dynamics simulations.³⁶

The vMCG wavefunction ansatz takes the following form:

Here, G_i represents separable, frozen GBFs, which are products of one-dimensional Gaussian functions for each coordinate q and |s〉 is a vector representing the associated electronic state basis. These GBFs are parameterized by time-dependent quantities such as position, momentum, and width. As the GBFs evolve over time, they provide an accurate description of the system's nuclear wavefunction. The method employs a variational principle to determine the positions of these Gaussians over time, keeping the nuclear basis set size optimally small.³⁷ This distinguishes vMCG from other methods like full multiple spawning (FMS)³⁸ and coupled coherent states,³⁹ which rely on classical trajectories for Gaussian centers.

The time-dependent Schrödinger equation (TDSE) is solved using the vMCG ansatz (eqn (3)) and the Dirac–Frenkel variational principle:⁴⁰

This principle leads to coupled equations of motion for the expansion coefficients A^(s)_i and the GBF parameters, resulting in quantum “trajectories” for the GBF centers. These “trajectories” differ from classical ones, as the GBFs are variationally coupled, which allows the method to capture quantum coherence and avoid the need for large basis sets typical in trajectory-based methods. If the Hamiltonian matrix elements 〈G_i|H|G_j〉 are calculated exactly, the vMCG method can, in principle, converge to the numerically exact solution of the TDSE.³⁷ In practice, however, we do not possess information of the full potential energy surfaces and we use the local harmonic approximation (LHA) to evaluate the matrix elements. The potential V(q) is expanded around the center of each GBF (G_i(q₀)) up to second order:

as if each Gaussian were feeling its own quadratic vibronic coupling potential. The zeroth, first and second order energy terms are calculated on-the-fly as electronic energies, gradients, nonadiabatic couplings and hessians, using ab initio quantum chemistry software. This avoids the need for precomputed potential energy surfaces, while maintaining analytical potential form for the integrals over the conformational space. The DD-vMCG (where DD stands for direct dynamics, meaning on-the-fly) method integrates these features by dynamically calculating the potential energy surfaces as the nuclear wavefunction evolves, making the approach computationally feasible for complex molecular systems.⁴¹

In DD-vMCG, the wavefunction is propagated on smooth diabatic potential energy surfaces. The propagation diabatisation technique⁴² is used to construct global diabatic states from adiabatic ones along the nuclear trajectories, which is critical to avoid the geometric phase issues associated with the adiabatic picture, particularly near conical intersections.⁴³ New electronic structure points, including energies, gradients, Hessians, and nonadiabatic couplings, are calculated on-the-fly and stored in a database. If a new geometry is encountered, it is compared to stored points, and interpolation between the diabatic energies is performed using Shepard interpolation.⁴⁴ To further reduce computational cost, the Hessian, one of the most computationally expensive quantities to calculate, is calculated at the initial structure and then approximated during the dynamics using a Hessian updating scheme based on gradient information.^45,46 This practice, combined with avoiding redundant quantum chemistry calculations via database interpolation, dramatically improves the efficiency of the method while maintaining accuracy.

The DD-vMCG method is implemented within the Quantics package developed in the Worth group.⁴⁷ The simulation begins by projecting the nuclear ground state onto the first excited state, which initializes a single Gaussian with amplitude, while the other gaussians in the expansion have zero initial amplitude but are displaced in momentum space to prepare the system for propagation. Thanks to the variational nature of the equations of motion, unpopulated Gaussians evolve and contribute to the dynamics without the need for explicit sampling as in trajectory-based methods. This allows the method to optimally describe the evolving wavepacket during the dynamics.

3 Dynamic results

Initially, DD-vMCG quantum dynamics propagations were performed using four Gaussian basis functions (GBFs). These preliminary simulations generated a small database of electronic structure points, which served as the foundation for subsequent simulations with an increasing number of GBFs. As the wavepacket explored new regions of the potential energy surface, additional points were added to the database, gradually expanding its coverage.

In the final set of DD-vMCG simulations, the dynamics was propagated for up to 300 fs. The largest database, which contained 48 857 electronic structure points, was constructed incrementally by adding points from previous simulations. Specifically, the database for simulations with 8, 16, and 20 GBFs was built upon the earlier 4 GBF simulation, with each successive simulation expanding the database created by the preceding one. This strategy ensured that the wavepacket had access to an increasingly comprehensive representation of the PES as the number of GBFs grew.

The key outputs analyzed in this study are: diabatic electronic populations of the two involved electronic states, one-dimensional cuts of the PES connecting the Franck–Condon region to the conical intersection, the geometries of the database points and the time evolution of the expectation values for every degree of freedom (in the SI). Additionally, the trans–cis photoisomerisation process was characterized by monitoring the evolution of the dihedral angle and the bond length alternation parameter, focusing on the weighted average of the Gaussian centers to capture the isomerisation pathway.

3.1 Evolution of state populations and geometrical observables

The DD-vMCG dynamics simulations provide access to the full evolution of the Gaussian wave packet. This GWP is represented as a linear combination of frozen multiconfigurational GBFs. Since these GBFs are expanded in the nuclear basis of the molecule's normal modes, we can directly extract the expectation values and uncertainties of these modes. At each time step, we determine whether the Gaussian center is close to any previously calculated ab initio points stored in the database. This involves applying a unitary transformation from normal mode coordinates to cartesian coordinates. Using the cartesian geometries of the GBF centers and their weight in the frozen Gaussian expansion, we calculate the average torsional angle asand its time dependent standard deviation aswhere 〈x〉 (t) represents the time-dependent average of the torsional angles among atoms C2, C3, C4 and C5 (line marked as “average” in Fig. 3), StdDev(t) represents the time-dependent standard deviation, GGP are the Gaussian gross populations – the square of the A_i coefficients in eqn (1) –, and N is the total number of Gaussians (20 for the plots in Fig. 3 and 16 for Fig. 6).

Fig. 3 (a) Diabatic electronic population decay for the initially excited state (ππ*), (b) evolution of expectation value and uncertainty (standard deviation) of torsional normal mode q₄, (c) average torsional angle (C2–C3–C4–C5) in degrees and standard deviation of the Gaussian centers in cartesian coordinates, (d) evolution of bond length alternation (BLA) parameter in Å defined as r_C1–C2 − r_C2–C3 + r_C3–C4 − r_C4–C5 + r_C5–N. For DD-vMCG calculations with 20 GBFs.

The bond length alternation (BLA) parameter has been calculated following a similar approach. First, the BLA value for each Gaussian center was determined as x_i(t) = r_C1–C2 − r_C2–C3 + r_C3–C4 − r_C4–C5 + r_C5–N where r_ij represents the bond lengths between the respective atoms. This value was then weighted by the Gaussian gross population associated with each center.

Fig. 3(a) illustrates the ability of the DD-vMCG method to capture the population transfer from the S₁ state to the S₀ state, which decreases from 100% to 60% over 250 fs. At first glance, this might suggest that isomerisation is incomplete within this timescale. However, insets (b) and (c) of Fig. 3 show the expectation value of the torsional normal mode and the average torsional angle, respectively. These data confirm that by 250–300 fs, the PSB3 molecule has fully isomerised from the trans to the cis form, as the torsional angle transitions from approximately 180 degree to 0.

Although the torsional motion is inactive prior to 150 fs, Fig. 3(d) shows that bond lengths begin to stretch and elongate within the first 25 fs. During this period, the bond length alternation (BLA) parameter increases from values below 1.2 – characteristic of the trans-isomer – to almost 1.6, corresponding to the S₁ minimum. By 200 fs, the system stabilises, oscillating around a BLA value of 1.4, which corresponds to the cis-isomer. The narrow standard deviations indicate that the isomerisation is effectively complete within this timeframe.

When examining the time evolution of the normal modes of PSB3 (Fig. S6–S11), the system can be divided into tuning modes, which carry the gradient and move the wavepacket away from the Franck–Condon region, and coupling modes, which couple different electronic states, facilitate population transfer, and modulate the width of the wavepacket. By analyzing how the expectation value and standard deviation of each normal mode change over time, it is possible to track the most active modes during the dynamics.

During the first 30 fs, the dynamics is primarily driven by mode v27 (tuning) and mode v36 (coupling), which elongate and shorten the central CC bonds. Between 30 and 100 fs, the dynamics become more complex, involving additional tuning modes (v6, v16, and v20) and coupling modes (v16, v21, v23, and v34). These modes predominantly consist of CC stretching motions but gradually incorporate torsional motions, such as mode v16, which corresponds to a slight NH₂ torsion. The normal modes of the PSB3 cation are depicted in Fig. S5 and Table S5.

From 100 to 200 fs, the previously mentioned tuning modes continue to play a significant role, but at 115 fs, the main torsional mode, v4, emerges as a coupling mode (depicted in Fig. 3(b)). This mode regulates the width of the wavepacket and facilitates moderate population transfer. Between 200 and 210 fs, mode v4 becomes the dominant tuning mode, exhibiting the largest amplitude motion, while v23 and v19 govern the coupling between electronic states.

After 220 fs and continuing until the end of the dynamics (300 fs), modes v32 and v30 (CH stretches) display the largest amplitudes, likely driving the system toward the cis structure minimum.

3.2 Geometries of the database points

Since the DD-vMCG method stores all ab initio and interpolated points in a database containing geometries, energies, gradients, couplings, normal mode displacements, and Hessians we can use these data to represent the conformational space explored during the dynamics. In Fig. 4, we depict the torsional angles and BLA parameters of the 48 857 database structures, represented as a 2D histogram with 30 bins per dimension (approximately 6 degree of torsion per bin and 0.07 angstrom per BLA bin). The color of each bin indicates the frequency of data points within that bin, with the darkest color corresponding to bins containing 400 or more structures. The crosses in the plot represent key geometries: the initial trans-PSB3 structure (S0_trans), the final cis-PSB3 structure (S0_cis) on the ground electronic state and the minimun energy conical intersection (MECI_db), which corresponds to the structure with the lowest energy gap between the two electronic states. These three structures are found in our database during the DD-vMCG/20GBFs dynamical simulation. The figure also shows the ab initio optimised trans-PSB3 (S1_trans), cis-PSB3 (S1_cis) structures for the excited electronic states and minimum energy conical intersection (MECI_opt). Optimised structures were calculated using a SA(2)-CAS(6,6)SCF/6-31G with Molpro³⁵ and their energies are given in Table 1.

Fig. 4 Histogram of the torsional angles (C2–C3–C4–C5) and bond alternation angles (BLA) calculated as BLA = r_C1–C2 − r_C2–C3 + r_C3–C4r_C4–C5 + r_C5–N for the 48 857 database structures. S0_trans and S0_cis represent the initial and final trans and cis-PSB3 structures on the ground electronic state found in the DD-vMCG/20GBFs dynamical simulation. MECI_db is the minimum energy conical intersections also found in the dynamics. S1_trans, S1_cis are the ab initio optimised trans and cis-PSB3 on the excited state, respectively. MECI_opt is the ab initio optimised minimum energy conical intersection. MECI_Liu corresponds to the conical intersection obtained by Liu et al.²⁵Ab initio optimisation calculations were done using a SA(2)-CAS(6,6)SCF/6-31G with Molpro program.³⁵

Notably, around the top right yellow cross representing the minimum of the excited state for the trans-PSB3 geometry, we observe a high density of structures with increasing BLA values, consistent with the dynamical trends shown in Fig. 3(d). The MECI_db identified in our simulation lies near the center of the most frequently visited regions of the database. Furthermore, the high density of structures near 0 degree torsion and a BLA value of 1.4 Å indicates successful isomerisation to the cis-PSB3 geometry. Another prominent region in the conformational space corresponds to around 60 degrees of torsion and a BLA of 1.4 Å, which lies close to the optimised MECI geometry, MECI_opt. We further analysed the number of database points that lay around the yellow crosses on the map (using an interval of ±3% BLA and torsion values of the cross geometries) and found that a 3.6% of the structures lay around the S1_trans structure (1939 structures), 4.0% around the S0_cis structure (1776 geometries), 0.2% around the MECI_opt geometry (112 geometries) and 0.4% around the MECI_db structure (188 geometries), only visited slightly more frequently than the ab initio optimised structure.

The excitation energies of the cis, trans isomers and MECI structures that we report in this work are in good agreement to those reported by other authors, with MECI energies laying at 2.87,²⁵ 2.89,²⁷ and 2.87 eV.²⁶ However, it is clear to us that the deactivation pathway in our case is involving an energy barrier. From Fig. 4, the system is initially in S0_trans and gets trapped in the S1_trans minimum (almost 2000 structures of the database seem to be at the torsional and BLA values of the S1_trans minimum geometry found ab initio). After that, the system slowly manages to leave that local minimum and continues evolving towards torsional angles of 90–120 degree and BLA values of 1.4, starting to transfer population to the groundstate at 170 fs. At 300 fs, the electronic population is still 60 : 40 (S1 : S0), but the system seems to have relaxed to the S0_cis structure (1800 structures) without passing by the S1_cis minimum geometry. This suggests a similar deactivation mechanism that other authors report through a similar conical intersection, but the system is trapped for a while at the more stable S1_trans minimum and, from there, is able to relax through both the MECI_db minimum and the MECI_opt minimum.

The PES underlying the isomerization process, including the presence of a significant energy barrier, is visualised in Fig. 5. This energy barrier is a key factor in explaining the delayed isomerization observed in our DD-vMCG simulations, especially when compared to previous trajectory-based studies that report timescales as short as 60 fs.

Fig. 5 Energy profile for the two electronic states extracted from the database with 48 857 geometries. The x axis is the vector that connects the S0_trans (at 1.0 arbitrary units) and the MECI_db (at 2.0 arbitrary units) geometries.

Complementary information is shown in Fig. S13 and S14, which present the energy distributions along both the torsion and bond length alternation (BLA) coordinates. From this analysis, the minimum-energy conical intersection found via diabatic optimization (MECI_db) is clearly located at a torsional angle of 120 degree and a BLA value of around 1.4 Å. In addition, we observe a second region of near-degeneracy between the ground and excited states at a torsional angle of approximately 80 degree, which is consistent with the position of the MECI_opt obtained via adiabatic ab initio optimization, which aligns with earlier works. The existence of these two regions strongly suggests that two distinct excited state deactivation pathways are available, potentially depending on the initial wavepacket distribution and the dynamics method employed. Notably, the wavepacket in our simulations initially evolves along coordinates that do not grant immediate access to either MECI, due to the barrier along the torsional coordinate. This delayed access, coupled with our quantum treatment, leads to a slower population transfer than that observed in mixed quantum-classical approaches. It is possible that trajectory-based simulations reach the MECI_opt region sooner due to the absence of zero-point energy constraints or the use of cartesian coordinates, which could lead to artificially accelerated internal conversion.

To further support these findings, we have provided the full set of ab initio points (over 45 000 per dynamics run) in both full and reduced dimensionality, openly available through Zenodo. We hope this data will encourage future efforts, including machine learning PES fitting or comparative studies with alternative quantum dynamics methods, to clarify the interplay between PES topology, coordinate choice, and simulation methodology.

3.3 Dynamics on reduced dimensionality

To investigate the impact of including or neglecting degrees of freedom on the system's dynamics, we selectively froze several normal modes and re-ran the simulations. Initially, we performed a calculation across the full conformational space, encompassing 36 normal modes, using a nuclear basis of 16 frozen Gaussians. The results, depicted as the dark purple “36mode” solid line in Fig. 6, are based on a database containing 45 000 records. This simulation demonstrated that the system relaxed to the ground state on a similar timescale (approximately 250 fs) as the previously reported calculations with 20 Gaussians. Subsequently, reduced-dimensionality calculations were carried out, focusing on 13 (q4,q9,q15,q16,q19,q20,q22,q23,q27,q29,q30,q32,q36), 19 (q1–q4,q6–q9,q15,q16,q19–q27), and 27 modes (q1–q27). The calculations were run independently – not based on the database of the simulations in full dimensionality – and resulted in databases of 42 701 (13mode) and 44 929 (19 and 24mode), the latter being shared by both reduced-dimensionality calculations. None of these simulations successfully relaxed to locate the cis-isomer. Notably, the 27-mode calculation excluded only the CH stretch modes, highlighting the crucial role of these modes in the isomerisation mechanism. We note that the challenges in reducing the dimensionality of the system may arise from the choice of nuclear basis vectors, which do not fully isolate the torsion from other coordinates. This limitation might be alleviated by employing a different nuclear basis or by constructing more tailored collective coordinates.

Fig. 6 (a) Diabatic electronic population decay for the initially excited state (ππ*/S₁), (b) average torsional angle (C1–C2–C3–C4) in degrees and standard deviation of the Gaussian centers in cartesian coordinates, (c) evolution of bond length alternation (BLA) parameter in Å defined as r_C1–C2 − r_C2–C3 + r_C3–C4 − r_C5–N. The dashed lines correspond to DD-vMCG simulations with 16GWPs with a reduced subset of normal modes, freezing some of the degrees of freedom. The solid line corresponds to the DD-vMCG/16GWPs simulation in full dimensionality.

Conclusions

This study investigated the photochemical trans–cis isomerisation of the PSB3 retinal model using nonadiabatic quantum dynamics simulations based on the DD-vMCG method. The main findings can be summarized as follows: upon vertical excitation from the Franck–Condon region to the first electronic state (ππ*/S₁), the wavepacket initiates its dynamics, first being trapped in the minimum of the excited state with larger BLA values and then evolving toward the conical intersections (shorter BLA and 120/80 degree of torsion respectively) which facilitates population transfer back to the ground state (π²/S₀). The isomerisation is driven by coupled vibrational modes, notably torsional and CC stretching modes. These modes guide the wavepacket through key regions of the PES associated with the isomerisation pathway. Based on reduced dimensionality results, we observed that CH stretches are crucial for a succesful isomerisation.

While our simulations reveal significant S₁-to-S₀ transfer only after 200 fs, mixed quantum-classical approaches such as ab initio multiple spawning (AIMS)²⁵ and Tully surface hopping²⁷ typically report much faster internal conversion, with nearly complete transfer occurring on similar timescales. Importantly, these comparisons are made using the same electronic structure method (CASSCF(6,6)/6-31G) as employed in previous trajectory-based studies of this system. This ensures that discrepancies are not attributable to the level of electronic structure theory, but rather to the dynamics methodology itself. The final cis-PSB3 geometry is achieved by 250–300 fs, with key signatures such as a torsional angle near 0 degree and bond length alternation (BLA) values stabilizing at ≈1.4 Å, confirming successful isomerisation, as depicted in Fig. 7.

Fig. 7 Average structure of PSB3 (averaged over the 20 Gaussian centers with a weight proportional to the amplitude of each Gaussian in the linear combination – A_i in eqn (1) –) for the DD-vMCG/20GWP propagation, showing how trans-PSB3 isomerises to cis-PSB3 in 300 fs.

In conclusion, in this work we tested the DD-vMCG method on the trans-PSB3 system, which has been widely studied using mixed quantum–classical trajectory methods. The tested method is able to predict the photoisomerisation of the retinal toy model, although discrepancies in timescales and pathway exploration between quantum and mixed quantum–classical methods are observed. These results are not entirely unexpected, since the DD-vMCG method – although still in its early stages of development – is a fully variational wavepacket approach which may yield different results than independent trajectory methods. Differences between our quantum-dynamics results and those obtained with independent trajectory approaches may also arise from other factors, such as the diabatic-to-adiabatic transformation, the choice of nuclear basis (normal modes vs. cartesian coordinates), and the intrinsically quantum nature of our propagation scheme. Similar discrepancies have also been reported in previous quantum dynamics studies using grid-based MCTDH methods for molecules like ethylene^48,49 and 1,1-difluoroethylene.^50,51 These findings reinforce the importance of continued investigation into the origins of such differences and how they relate to the regions of configurational space sampled by each method. Overall, we find these results encouraging and supportive of further development and application of on-the-fly quantum dynamical approaches to complex nonadiabatic processes.

Author contributions

María Mallo: formal analysis (equal); methodology (equal); writing – original draft (lead). Susana Gómez-Carrasco: formal analysis (equal); writing – review & editing (equal). Sandra Gómez: formal analysis (equal); methodology (equal); funding acquisition (lead); methodology (equal); project administration (lead); supervision (lead); writing – review & editing (equal).

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Optimised structures for the S₀, S₁ states and the conical intersection (Section S1), vibrational frequencies (Section S2), convergence of diabatic electronic populations (Section S3), time evolution of each normal mode (Section S4), the evolution of dihedral angles for each frozen Gaussian with its weight in the expansion (Section S5), and the distribution of 48 857 database geometries along the torsional angle and the BLA coordinate and their energies (Section S6). See DOI: https://doi.org/10.1039/d5cp02117k

All data necessary to reproduce the results presented in this article are available at Zenodo at https://doi.org/10.5281/zenodo.15582165. The dataset includes the database created during the dynamics in.sql format for the 20GWP and 16GWP cases and the SI as a PDF file.

Acknowledgements

SG acknowledges NextGenerationEU funds (María Zambrano Grant for the attraction of international talent) and the grant “Programa Propio C1” from Universidad de Salamanca. The authors thank the funding by Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033) grant no. PID2020-113147GA-I00, the EPSRC under the COSMOS programme grant (EP/X026973/1) and the COST action CA21101” Confined molecular systems: from a new generation of materials to the stars (COSY)” supported by COST (European Cooperation in Science and Technology). The authors thank Graham A. Worth and Olivia Bennett for their valuable input and feedback.

Notes and references

1. G. Wald, Science, 1968, 162, 230–239.
2. A. Warshel, Nature, 1976, 260, 679–683.
3. M. Huntress, S. Gozem, K. R. Malley, A. E. Jailaubekov, C. Vasileiou, M. Vengris, J. Geiger, B. Borhan, I. Schapiro, D. Larsen and M. Olivucci, J. Phys. Chem. B, 2013, 117(35), 10053–10070.
4. K. C. Hasson, F. Gai and P. Anfinrud, Proc. Natl. Acad. Sci. U. S. A., 1996, 93(26), 15124–15129.
5. R. W. Schoenlein, L. A. Peteanu, R. A. Mathies and C. V. Shank, Science, 1991, 254(5030), 412–415.
6. K. Hayashi, M. Mizuno, H. Kandori and Y. Mizutani, Angew. Chem., 2022, 134(33), e202203149.
7. R. Liang, F. Liu and T. J. Martínez, J. Phys. Chem. Lett., 2019, 10, 2862–2868.
8. B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem., 2007, 58, 613–634.
9. M. Ruckenbauer, M. Barbatti, T. Muüller and H. Lischka, J. Phys. Chem. A, 2013, 117, 2790–2799.
10. S. Gozem, P. J. M. Johnson, A. Halpin, H. Luk, T. Morizumi, V. Prokhorenko, O. Ernst, M. Olivucci and R. Miller, J. Phys. Chem. Lett., 2020, 11(10), 3889–3896.
11. T. Vreven, F. Bernardi, M. Garavelli, M. Olivucci, M. Robb and H. Schlegel, J. Am. Chem. Soc., 1997, 119, 12687–12688.
12. Y. Liu and C. Zhu, Phys. Chem. Chem. Phys., 2021, 23, 23861–23874.
13. R. Palombo, L. Barneschi, L. Pedraza-González, D. Padula, I. Schapiro and M. Olivucci, Nat. Commun., 2022, 13, 6652.
14. O. Kobayashi and S. Nanbu, Chem. Phys., 2015, 461, 47–57.
15. A. Valentini, D. Rivero, F. Zapata, C. García-Iriepa, M. Marazzi, R. Palmeiro, I. Fdez. Galván, D. Sampedro, M. Olivucci and L. M. Frutos, Angew. Chem., Int. Ed., 2017, 56, 3842–3846.
16. J. K. Yu, R. Liang, F. Liu and T. Martínez, J. Am. Chem. Soc., 2019, 141(45), 18193–18203.
17. L. E. Cigrang, B. F. E. Curchod, R. A. Ingle, A. Kelly, J. R. Mannouch, D. Accomasso, A. Alijah, M. Barbatti, W. Chebbi, N. Došlic, E. C. Eklund, S. Fernandez-Alberti, A. Freibert, L. González, G. Granucci, F. J. Hernández, J. Hernández-Rodríguez, A. Jain, J. Janoš, I. Kassal, A. Kirrander, Z. Lan, H. R. Larsson, D. Lauvergnat, B. L. Dé, Y. Lee, N. T. Maitra, S. K. Min, D. Peláez, D. Picconi, U. Raucci, Z. Qiu, P. Robertson, E. S. Gil, M. Sapunar, P. Schürger, P. Sinnott, S. Tretiak, A. Tikku, P. Vindel-Zandbergen, G. A. Worth, F. Agostini, S. Gómez, L. M. Ibele and A. Prlj, J. Phys. Chem. A, 2025, 129(31), 7023–7050.
18. M. Robb, et al. , Faraday Discuss., 1998, 110, 51–70.
19. M. Huix-Rotllant, M. Filatov, S. Gozem, I. Schapiro, M. Olivucci and N. Ferreé’, J. Chem. Theory Comput., 2013, 9, 3917–3932.
20. P. Ertl, Collect. Czech. Chem. Commun., 1990, 55, 2874–2879.
21. S. Gozem, F. Melaccio, R. Lindh, A. Krylov, A. Granovsky, C. Angeli and M. Olivucci, J. Chem. Theory Comput., 2013, 9(10), 4495–4506.
22. P. Zhou, J. Liu, K. Han and G. He, J. Comput. Chem., 2014, 35, 109–120.
23. D. Grabarek, E. Walczak and T. Andruniów, J. Chem. Theory Comput., 2016, 12(5), 2346–2356.
24. D. Tuna, D. Lefrancois, L. Wolanski, S. Gozem, I. Schapiro, T. Andruniów, A. Dreuw and M. Olivucci, J. Chem. Theory Comput., 2015, 11(12), 5758–5781.
25. L. Liu, J. Liu and T. J. Martinez, J. Phys. Chem. B, 2016, 120, 1940–1949.
26. M. Filatov, S. K. Min and K. S. Kim, J. Chem. Theory Comput., 2018, 14, 4499–4512.
27. M. Barbatti, M. Ruckenbauer, J. J. Szymczak, A. J. Aquino and H. Lischka, Phys. Chem. Chem. Phys., 2008, 10, 482–494.
28. F. Talotta, D. Lauvergnat and F. Agostini, J. Chem. Phys., 2022, 156(18), 184104.
29. F. Schautz, F. Buda and C. Filippi, J. Chem. Phys., 2004, 121(12), 5836–5844.
30. G. Giudetti, M. Mukherjee, S. Nandi, S. Agrawal, O. V. Prezhdo and A. Nakano, J. Chem. Inf. Model., 2024, 64, 7027–7034.
31. M. Garavelli, P. Celani, F. Bernardi, M. Robb and M. Olivucci, J. Am. Chem. Soc., 1997, 119, 6891–6901.
32. J. J. Szymczak, M. Barbatti and H. Lischka, J. Chem. Theory Comput., 2008, 4, 1189–1199.
33. D. Grabarek, E. Walczak and T. Andruniow, J. Chem. Theory Comput., 2016, 12, 2346–2356.
34. J. J. Szymczak, M. Barbatti and H. Lischka, J. Phys. Chem. A, 2009, 113, 11907–11918.
35. MOLPRO is a package of ab initio programs designed by H. -J. Werner and P. J. Knowles and with contributions from, J. Almlöf and R. D. Amos and A. Berning and M. J. O. Deegan and F. Eckert and S. T. Elbert and C. Hampel and R. Lindh and W. Meyer and A. Nicklass and K. Peterson and R. Pitzer and A. J. Stone and P. R. Taylor and M. E. Mura and P. Pulay and M. Schütz and H. Stoll and T. Thorsteinsson and D. L. Cooper, version 2012.
36. G. A. Worth and B. Lasorne, Quantum Chemistry and Dynamics of Excited States: Methods and Applications, 2020, pp. 413–433.
37. G. A. Worth, M. A. Robb and I. Burghardt, Faraday Discuss., 2004, 127, 307–323.
38. T. J. Martinez, M. Ben-Nun and R. D. Levine, J. Phys. Chem., 1996, 100, 7884–7895.
39. J. A. Green, A. Grigolo, M. Ronto and D. V. Shalashilin, J. Chem. Phys., 2016, 144, 024111.
40. A. Raab, Chem. Phys. Lett., 2000, 319, 674–678.
41. G. Christopoulou, A. Freibert and G. A. Worth, J. Chem. Phys., 2021, 154, 124127.
42. G. W. Richings and G. A. Worth, Chem. Phys. Lett., 2017, 683, 606–612.
43. L. M. Ibele, E. Sangiogo Gil, B. F. E. Curchod and F. Agostini, J. Phys. Chem. Lett., 2023, 14, 11625–11631.
44. T. J. Frankcombe, M. A. Collins and G. A. Worth, Chem. Phys. Lett., 2010, 489, 242–247.
45. T. J. Frankcombe, J. Chem. Phys., 2014, 140, 114106–114108.
46. G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt and B. Lasorne, Int. Rev. Phys. Chem., 2015, 34, 269–308.
47. G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle and H.-D. Meyer, The QUANTICS Package, Version 1.1, University of Birmingham, Birmingham, U.K., 2015.
48. A. Viel, R. P. Krawczyk, U. Manthe and W. Domcke, J. Chem. Phys., 2004, 120, 11000–11010.
49. S. Gómez, E. Spinlove and G. Worth, Phys. Chem. Chem. Phys., 2024, 26, 1829–1844.
50. S. Gómez, L. M. Ibele and L. González, Phys. Chem. Chem. Phys., 2019, 21, 4871–4878.
51. S. Gómez, N. K. Singer, L. González and G. A. Worth, Can. J. Chem., 2022, 101, 745–757.

---

Footnote

† Present address: Instituto de Física Fundamental, CSIC, Calle Serrano 123, E-28006 Madrid, Spain.

---