Engineering the vibrational coherence of vision into a synthetic molecular device

The light-induced double-bond isomerization of the visual pigment rhodopsin operates a molecular-level optomechanical energy transduction, which triggers a crucial protein structure change. In fact, rhodopsin isomerization occurs according to a unique, ultrafast mechanism that preserves mode-specific vibrational coherence all the way from the reactant excited state to the primary photoproduct ground state. The engineering of such an energy-funnelling function in synthetic compounds would pave the way towards biomimetic molecular machines capable of achieving optimum light-to-mechanical energy conversion. Here we use resonance and off-resonance vibrational coherence spectroscopy to demonstrate that a rhodopsin-like isomerization operates in a biomimetic molecular switch in solution. Furthermore, by using quantum chemical simulations, we show why the observed coherent nuclear motion critically depends on minor chemical modifications capable to induce specific geometric and electronic effects. This finding provides a strategy for engineering vibrationally coherent motions in other synthetic systems.

The same data for E-2. Notice that, in this case however, the computed parameters (in brackets) are compared with a previously reported [1] precursor of E-2 which is "sterically close" to E-2 (i.e. a H atom is replaced with an O atom), since E-2 could not be crystallized.

E-2
In order to theoretically describe the spectroscopic and dynamic properties of the studied molecular switches a QM/MM-based protocol oriented to generate an average configuration of the solvent environment at the desired temperature was used to optimize the ground state geometries and generate initial conditions to initiate exited states dynamics. The proposed protocol [2] is based on the idea of Herbert et al. [3,4] of combining the Average Solvent Electrostatic Configuration (ASEC) model [5,6] and the free energy gradient method proposed by M. Nagaoka et. al. to optimize solvated molecular systems. [7][8][9] In this QM/MM model the solute is considered quantum mechanically (QM) while a molecular mechanic force field is used to describe the solvent environment (MM).
As stated by M. Nagaoka et. al., [7][8][9] the free energy gradient, or the forces F(q) acting on the QM part of the system, can be computed as: where q represents the coordinates of the QM part and G is the free energy of the system, which can be computed, in a very good approximation, as the average of the potential interaction, 〈 〉, between the QM and MM subsystems. In this context, the ASEC model [6,7] is used to compute, in a very efficient way, the average potential interaction 〈 〉. ASEC model is a mathematical construct called "ASEC configuration" created from a selected sampling of configurations, obtained via extensive molecular dynamics (MD), to mimic the effect of thermodynamic equilibrium conditions at a selected temperature. The solute is kept fixed along the MD in order to sample the environment, being subsequently optimized in the presence of the average external configuration (Supplementary Figure  3a). The GROMACS code [10] was used for the MD, performing an initial thermalization of 1 ns, followed by 4 ns of production in the NPT ensemble and standard room conditions. The QM/MM geometry optimization is performed using the quasi-Newton-Raphson method implemented in the MOLCAS-TINKER interface. [11,12] The QM subsystem is described using Complete Active Space Self Consistent Field (CASSCF) method. [13] For the MM subsystem, OPLS force field parameters are used. [14] Finally, the absorption energies of optimized structures are re-computed using second order perturbation theory (CASPT2) [15] with ionization potential electronic affinity (IPEA) set to zero. An active space comprising of 12 electrons in 11 orbitals is used for all QM calculations in current study.
Since the configurations of the environment selected to generate the ASEC configuration (100 solute-solvent uncorrelated configurations in our case) are obtained from the system into thermodynamic equilibrium conditions, a Boltzmann distribution of the absorption energies is expected to be observed from those structure. Therefore, these 100 configurations can be used to initiate excited state dynamics. For the two molecular switches, absorption bands were constructed (Supplementary Figure 4). A single configuration from each band, representing the average was selected to map the potential energy surfaces. The selected two configurations were re-optimized on the ground state. The setup used for optimization and subsequent QM/MM calculations is shown in Supplementary Figure 3b. The vertical excitation energies and oscillator strengths were computed for the ground state optimized structures. According to oscillator strengths, the transition from ground state to second excited state (S0 to S2) was highly probable in both switches. Therefore, minimum energy path (MEP) calculation and Franck-Condon trajectory for each of the two switches were initiated from S2. The MEP calculation resulted in minima on S2 and S1 for Z-1 and E-2 switches respectively ( Supplementary Figures 5c and 5d). Starting from these minima, torsional scans for each switch was performed by constraining the dihedrals C9'-C1'=C4-C5 and C2'-C1'=C4-C3 (see Supplementary Figure 6). The Franck-Condon trajectories were computed at 3-root-state-average CASSCF/6-31G*/OPLS level. This was performed using the Dynamix module implemented in the MOLCAS program and more details can be found in Schapiro et al. [16] Since the CASSCF method does not account for dynamic electron correlation, the energy profiles had to be recomputed using a method that doesn't suffer from this problem. Due to this reason, for each geometry along the trajectories, the energies were recomputed at 4-root-state-average CASPT2 and XMCQDPT2 levels (see Supplementary Note 3 for more details on the latter).

Supplementary Note 2: Computation of a representative QM/MM structure in solution
Parts a. and b. of Supplementary Figure 5 illustrate the steric clash/pretwist in Z-1 and geared/planar conformation in E-2. In Z-1, the steric clash of the nearly eclipsed hydrogens at C3 with the methyl groups at C2' as well as the clash of the methyl group and hydrogen at C5 and C8' respectively result in a strongly twisted S0 equilibrium conformation. Instead structure E-2 displays a substantially planarized S0 conformation, where the small H substituent at C5 assumes a geared configuration with respect to the two relatively bulky methyl groups at C2'. The corresponding equilibrium S1 structures (see parts c. and d. in Supplementary Figure 5) are computed by optimizing on the S1 state a single equilibrium S0 configuration obtained in the following way: (i) Selecting from the ASEC configuration (composed by 100 overlapping configurations) a single configuration component featuring the absorption maximum closest to the ASEC absorption maximum and (ii) performing a QM/MM geometry optimization in S0 starting from this structure and relaxing both the QM part and the mobile solvent cavity (the MM cavity of Supplementary Figure 3b). The resulting optimized structure has then an absorption maximum very close to the ASEC absorption (these two values are indicated with green and red circles in Supplementary Figure 4 respectively).
Supplementary Figure 5: Comparison of S0 optimized and S1 optimized geometries. Bond length values are in angstrom and dihedral angles are in degrees. The S0 (ground state) structures are the geometrically optimized ASEC structures. The S1 structures (excited state) are single configurations obtained by optimizing on S1 state a specifically constructed single S0 configuration. The corresponding Cartesian coordinates are given in the Supplementary Tables 1, 2, 4 and 5. Notice that the excited state structure of E-2 has actually been optimized on the S2 state as this state is actually the spectroscopic state at the CASSCF level (see Supplementary Figures 6c and 6d). Although the CASPT2 method accounts for dynamic electron correlation and results in better energies with respect to CASSCF, it fails to describe near degeneracy regions [17][18][19]. Due to this reason, the energy profiles along the trajectories were also computed at the extended multiconfigurational quasidegenerate second order perturbation theory (XMCQDPT2) [19]. This was performed using the Firefly computer package [20]. The resulting XMCQDPT2 energy profiles are reported in the Supplementary Figure 7. By comparing these with the corresponding CASPT2 energies reported in Figure 4, one can clearly see that crossings between potential energy surfaces at CASPT2 level have now become avoided crossings.