QTAIM and Stress Tensor Characterization of Intramolecular

Subscriber access provided by Binghamton University | Libraries

Article

QTAIM and Stress Tensor Characterization of Intramolecular Interactions Along Dynamics Trajectories of a Light-Driven Rotary Molecular Motor Ling Ling Wang, Guo Huan, Roya Momen, Alireza Azizi, Tianlv Xu, Steven R. Kirk, Michael Filatov, and Samantha Jenkins J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 06 Jun 2017 Downloaded from http://pubs.acs.org on June 8, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

QTAIM and Stress Tensor Characterization of Intramolecular Interactions Along Dynamics Trajectories of a Light-Driven Rotary Molecular Motor

Lingling Wang, Guo Huan, Roya Momen, Alireza Azizi, Tianlv Xu, Steven R. Kirk*, Michael Filatov* and Samantha Jenkins*

Key Laboratory of Chemical Biology and Traditional Chinese Medicine Research and Key Laboratory of Resource FineProcessing and Advanced Materials of Hunan Province of MOE, College of Chemistry and Chemical Engineering, Hunan Normal University, Changsha, Hunan 410081, China

*email: [email protected] *email: email: [email protected] *email: [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

Abstract A QTAIM and stress tensor analysis was applied to analyze intramolecular interactions influencing the photoisomerization dynamics of a light-driven rotary molecular motor. For selected non-adiabatic molecular dynamics trajectories characterized by markedly different S1 state lifetimes, the electron densities were obtained using the ensemble density functional theory method. The analysis revealed that torsional motion of the molecular motor blades from the Franck-Condon point to the S1 energy minimum and the S1/S0 conical intersection is controlled by two factors: greater numbers of intramolecular bonds before the hop time and unusually strongly coupled bonds between the atoms of the rotor and the stator blades. This results in the effective stalling of the progress along the torsional path for an extended period of time. This finding suggests a possibility of chemical tuning of the speed of photoisomerization of molecular motors and related molecular switches by reshaping their molecular backbones to decrease or increase the degree of coupling and numbers of intramolecular bond critical points as revealed by the QTAIM/stress tensor analysis of the electron density. Additionally, the stress tensor scalar and vector analysis was found to provide new methods to follow the trajectories and from this new insight was gained into the behavior of the S1 state in the vicinity of the conical intersection.


Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Light driven rotary molecular motors represent a novel class of molecules capable of performing unidirectional rotary motion as a consequence of photo-isomerization about a double bond1–9. The first models of molecular motors were based on overcrowded alkene framework4–8, where the double bond connecting two blades, the rotor and the stator, was strained by the steric repulsion. The photo-excitation breaks the π-bond, which results in a rapid rotary motion of the rotor blade with respect to the stator in the sense dictated by folding of the molecule in its ground state equilibrium conformation. The photo-isomerization proceeds further by nonadiabatic population transfer through conical intersection as was unequivocally demonstrated theoretically10,11 as well as experimentally12. Having reached the ground state metastable conformation the motor undergoes an intermediary thermal relaxation step and ends up in a ground state conformation with the rotor blade rotated through ca. 180° relative to its original orientation4,7,10,11; absorption of another photon followed by the thermal relaxation completes the full loop. The first models of molecular motors featured fairly low quantum efficiency, successful isomerizations per photon, on the order of 0.1 to 0.211,12, which was arguably caused by the necessity to undergo large scale pyramidalization distortion of the central bond required, alongside ca. 90° torsion, to reach the conical intersection10,11,13. In the literature, several alternative molecular motor designs were proposed where the necessity to pyramidalize was lifted13–16. In particular, it was argued that replacing the torsion-pyramidalization conical intersections by torsion-bond_length_alternation (BLA) ones, as a result of chemical modification of molecular frame, should lead to molecular motors capable of performing pure axial rotation of its blades and possessing greater quantum efficiency than their overcrowded alkene analogues13; these conjectures were later corroborated by the non-adiabatic molecular dynamics (NAMD) simulations16. These new motors performed pure axial rotations and the photochemical steps of the rotary cycle were found to be dominated by the fast BLA motion that enables ultrafast access to the S1/S0 intersection. This gave a handle on rational modulation of the photo-isomerization dynamics of molecular motors by purely chemical means; a goal that could previously be attained only for thermal relaxation of molecular motors. The NAMD simulations of chemically modified molecular motors, designated as NAIBP (N-alkylated indanylidene benzopyrrole) motors, see Scheme 1, revealed certain peculiarities of their dynamic behavior, which were not observed before16. Particularly, two types of excited state dynamics of the 3-[(2S)-2-fluoro-2methyl-1-indanylidene]-1-methyl-2-methylindole (F-NAIBP) motor were identified16; the former type is characterized by rapid progress of the molecule (achieved by rotation of the rotor blade) toward the S1 potential energy surface (PES) minimum and the nearby conical intersection, while the latter type displayed delayed 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

rotation after the molecule spent prolonged period of time flexing near the Franck-Condon (FC) point upon the photo-excitation. In both cases the molecule has successfully reached the S1/S0 conical intersection, however the S1 lifetime stemming from the two types of dynamics differed by almost a factor of two, ca. 300 fs vs. ca. 600 fs. As such a difference may have considerable implications for functioning of the motors, e.g. when immersed in a liquid phase, anchored on a substrate, or when working under load, it seems desirable to investigate the origin of the difference. It is the purpose of this work to investigate the aspects of molecular electronic structure which may be responsible for the difference in the S1 dynamics of F-NAIBP motors. The quantum theory of atoms in molecules (QTAIM)17,18 methodology is capable of revealing bonding interactions between individual functional groups and atoms in a molecule by analyzing the electron density distribution. In the functioning of light-driven switches and motors19 long weak closed-shell BCPs are often formed and annihilated along torsion coordinates thus enabling the switching of bond-path connectivity within a nuclear skeleton. As suggested in recent studies20 , the significance of these attractive interactions was underestimated in discussions on repulsive steric hindrance in chemistry. In a previous work the QTAIM and the stress tensor were applied to study the bonding patterns along the torsional pathway of fluorene molecular motor, the overcrowded alkene type, where it was shown that the scalar and vector-based QTAIM and stress tensor descriptors, such as the ellipticity ε, stress tensor stiffness Sσ and the stress tensor response βσ, can be sensitive probes of the electronic structure19. In this work, the QTAIM and the stress tensor will be applied to study the bonding patterns of the F-NAIBP molecular motor along the NAMD trajectories obtained in dihedral torsion angle θ in an earlier publication16. As the aspects of the S1 decay we are addressing in this work are caused by the dynamic effects, it is natural to investigate the electronic structure along the actual NAMD trajectories, rather than minimum energy pathways, e.g. the torsional path in19. For this purpose, two representative trajectories are selected from the pool of all NAMD trajectories (ca. 400) reported in dihedral torsion angle θ in an earlier publication16: a fast (F) and a slow (S) trajectories also shown in Fig. 4 of dihedral torsion angle θ in an earlier publication16. The geometries along these trajectories were obtained in dihedral torsion angle θ in an earlier publication16 with the use of the OM2/MRCI semiempirical quantum chemical method21–23 in connection with the trajectory surface hopping (TSH) NAMD formalism24,25. As semiempirical methods are not suitable for providing the density input for QTAIM analysis, in this work we employ the ensemble DFT method, the state interaction state-averaged spinrestricted ensemble-referenced Kohn-Sham (SI-SA-REKS or SSR, for brevity) method26,27, to obtain the electron densities along the TSH-OM2/MRCI trajectories. The two methods, SSR and OM2/MRCI, were shown to yield relative energies of different S0 and S1 conformations of F-NAIBP within a few kcal/mol from each other16,28; therefore, the OM2/MRCI geometries can be used with confidence for the SSR calculations. 4 ACS Paragon Plus Environment

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The purpose of this study is to analyze what factors define the observed differences in the dynamics along the fast (F) and slow (S) trajectories of F-NAIBP motor. Why does the molecule stay near the FC point for more than 400 fs when following the slow trajectory, but rapidly slides down the S1 PES on the fast trajectory? Are intramolecular interactions responsible for these differences? Addressing these questions seems important for (a) better understanding of the S1 dynamics of molecular motors and (b) improving the design of future molecular motors by chemically modulating the aspects of their electronic structure causing different dynamics behaviors. 2. Theory and Methods In this section, we briefly recapitulate the salient features of the methodologies used in this work. For a more in depth account of the methodologies the reader is advised to visit the publications cited herein.

2.1 QTAIM eigenvalue and stress tensor BCP descriptors; ellipticity ε, the total local energy density H(rb), stiffness S and stress tensor stiffness Sσ We are concerned with the properties of the Hessian matrix of the total electronic charge density ρ(r) evaluated at each critical point that is a stationary point of the total electronic charge density, at which the gradient ∇ρ(r) vanishes. Diagonalization of the Hessian matrix of ρ(r) gives the set of ordered eigenvalues λ1 < λ2 < λ3 with |λ1| > |λ2| < λ3, with a corresponding set of eigenvectors e1, e2, e3, with the Laplacian ∇2ρ(r) of ρ(r) being the algebraic sum of these eigenvalues evaluated at one of the four types of critical point. The nature of these critical points is revealed by the ordered set of the eigenvalues (λ1, λ2 and λ3) of the corresponding Hessian matrix. These critical points are labeled using the notation (R, ω) where R is the rank of the Hessian matrix the number of distinct non-zero eigenvalues and ω is the signature (the algebraic sum of the signs of the eigenvalues); the (3, -3) [nuclear critical point (NCP), a local maximum generally corresponding to a nuclear location], (3, -1) and (3, 1) [saddle points, called bond critical points (BCP) and ring critical points (RCP), respectively], and (3, 3) [the cage critical points (CCP)]. The eigenvector e3 indicates the direction of the bondpath at the BCP. The most and least preferred directions of electron accumulation are e2 and e1, respectively29–31. The ellipticity defined as ε = |λ1|/|λ2| – 1, where λ1 and λ2 are negative eigenvalues of the corresponding Hessian matrix at a BCP, provides the relative accumulation of ρ(r) at the bond critical point denoted by ρ(rb) in the two directions perpendicular to the bond-path at the BCP. Previously, we defined a bond-path stiffness, S = λ2/λ3, as a measure of rigidity of the bond-path32. A degree of covalent character is indicated for H(rb) < 033,34, where the total local energy density H(rb) = G(rb) + V(rb) where G(rb) and V(rb) are the local kinetic and virial energy densities, respectively. Previously, one of the current authors used H(rb) < 0 to quantify the degree of covalent character present in the hydrogen-bonding 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

in ice Ih to provide a quantitative explanation of the unusually high strength of ice Ih hydrogen-bonds35. This result was explained in terms of coupling i.e., cooperativity between the hydrogen-bonds and adjoining (σ) O-H bond paths that share an H NCP causing covalent bond character to ‘leak’ from the O-H BCPs to the hydrogenbond BCPs34. More recently we demonstrated a two way exchange of chemical character; covalent character was donated from a sigma O-H BCP to the hydrogen bond H--O BCP that shared a proton and hydrogen bond character donated back to the sigma O-H BCP from the hydrogen bond H--O BCP36. Coupling between the stator and rotor BCPs may still be present when H(rb) > 0, but to a reduced extent compared with instances in which H(rb) < 0. Another factor influencing the coupling of stator and rotor is the numbers of intramolecular BCPs present before the hop time that may vary due to dynamic exchanges in ρ(rb) distributions. Throughout this work, we use the notation H--H BCP37,38, F--H BCP39 and H--C BCP40 for BCPs with H(rb) < 0, because usually H---H BCPs have an H(rb) > 0; respectively so we use ‘--‘ and ‘---‘ in the text to distinguish BCPs with H(rb) > 0 from BCPs with H(rb) < 0. The terminology ‘--‘ and ‘---‘ refer to closed-shell BCPs which by definition always possess values of the Laplacian ∇2ρ(rb) > 0 and can possess H(rb) > 0 or H(rb) < 0. Conversely ‘-‘ always refers to shared-shell BCPs which by definition always possess values of the Laplacian ∇2ρ(rb) < 0 and H(rb) < 0. The quantum stress tensor, σ(r), which is directly related to the Ehrenfest force ∇by the virial theorem therefore provides a physical explanation of the low frequency normal modes that accompany structural rearrangements32. Diagonalization of the stress tensor, σ(r), returns the principal electronic stresses. In this work we use the definition of the stress tensor proposed by Bader41 to investigate the stress tensor properties within the QTAIM partitioning scheme. Thus, analogously to the QTAIM descriptor, we calculate a stress tensor stiffness, Sσ = |λ1σ|/|λ3σ|, which has been found as a good descriptor of the ‘resistance’ of the bond-path to the twist, as well as the eigenvector, e1σ indicating the direction of the π-bond42. Previously, it was found that the stress tensor stiffness, Sσ produced results that were in line with physical intuition32,42. We can represent the QTAIM topologies of sets of isomers of allowed, forbidden and unstable solutions to the Poincaré-Hopf relationship43–45 using the quantum topology phase diagram (QTPD)46–48. For molecules and clusters, the relation is expressed as:

n-b+r–c=1,

(1)

where the four types of critical points; n, b, r, and c are given by the numbers of NCPs, BCPs, RCPs, and CCPs, respectively. The sum of the number of these critical points for a given molecular graph is referred to as the topological complexity ∑brc. For a given collection of isomeric molecular graphs the number of BCPs and RCPs 6 ACS Paragon Plus Environment

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


given by b and r respectively, contained in each molecular graph are plotted along the x-axis and y-axis. In this work we consider the structure present at each time step to be an isomer. In addition to QTPDs it is possible and convenient to create plots of the variation of the number of BCPs (b) with time as well as the variation of the number of RCPs (r) with time. In this way the variation of the complicated and changing bonding topology with time can be followed more easily than using the QTPD alone. It would be possible to create a 3-D QTPD with time as the third axis, but this would not be as straightforward to interpret. The concept of quantum topology (QT) previously defined within QTAIM46, will be used to determine whether a structure is 3-DQT, 2-DQT or 1-DQT dependent on the presence of one or more CCPs, RCPs or BCPs respectively. The subscript ‘QT’ is used to denote quantum topology to distinguish from Euclidean geometry. A molecular graph that possesses a 3-DQT quantum topology contains all four types of critical points; 2-DQT – NCPs, BCPs, and RCPs are present; 1-DQT –contains NCPs and BCPs and a 0-DQT comprises only NCPs. 2.2 QTAIM eigenvector and stress tensor BCP descriptors; the response β, the stress tensor response βσ Recently some of our authors have introduced a measure of the response β of the total electronic charge density ρ(rb) to dihedral torsion about a pivot-BCP for different electronic states42. In this work, we first establish a common reference Cartesian coordinate frame for all structures using a triplet of atoms {U,V,W} thus: the position of atom U is always the origin of the reference frame and the straight line joining atoms U and V gives the direction of the reference x-axis X′. The component of the line U-W perpendicular to the line U-V is the reference y-axis Y′ and the vector cross-product Z′ = X′× ×Y′ completes the normalized reference frame. All vector quantities are projected into this reference frame before further use, and we will drop the use of the ‘prime’ notation for clarity in all subsequent descriptions. In this work, the {U,V,W} set was chosen as atoms C1,C2,C3, see Scheme 1. The response β is defined to be the angle formed between the eigenvector e2 of the BCP and a reference direction, in this case the [0.0,-1.0,0.0], or –Y, see Figure 1. The response β of the total electronic charge density ρ(rb) to the torsion of the bond-path is mapped into the range -90.0° ≤ β ≤ 90.0° using the additional proviso that β→(β - 180.0°) if β > 90.0°, using the definition for the response β:

β = arccos(e2 · y) .

(2)

In equation (2), the e2 of the pivot-BCP determines the preferred direction of density accumulation and therefore torsion of the bond-path and y is a unit vector anti-parallel to the Y-axis. Therefore, the QTAIM vector-based interpretation of a double bond and single bond, the response β, relates to the degree of 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

detachment of the {e1,e2,e3} bond-path framework from the containing nuclear skeleton. The detachment of the {e1,e2,e3} bond-path framework from the containing nuclear skeleton will result in differences between response β and the dihedral torsion angle θ. It is also possible to define an analogous stress tensor response βσ, (90.0° ≤ βσ ≤ 90.0°), where e1σ determines the preferred direction of electronic motion:

βσ = arccos(e1σ · y)

(3)

The chemical interpretation of values of βσ < 0 is that the {e1σ,e2σ,e3σ} bond-path framework is rotating towards the nuclear skeleton defined reference axis Y. Conversely, values of βσ > 0 correspond to the {e1σ,e2σ,e3σ} bond-path framework rotating away from the nuclear skeleton. For the S0 state we have previously observed in ring-opening reactions, where there were no excited states nor conical intersection seams, that the response βσ changed sign and magnitude from -90º to +90º

19

at the transition state when following the dependence of βσ

along the reaction pathway. Therefore, any change in sign of the response βσ should not be confused with the fact that after crossing through the CI seam the electronic states swap. In this previous investigation into ring opening reactions, in the ground state, the sign and magnitude of βσ changed discontinuously from -90º to +90º at the transition state. Another system that contained excited states, but no CI seam, was an investigation into the photoisomerization of the light driven fluorene molecular rotary motor49, where both the S0 and S1 state changed discontinuously from βσ = -90.0º to βσ = +90.0º at the energy maxima and energy minima respectively. Conversely, the earlier work on the 11-cis retinal protonated Schiff base demonstrated that for S1 state the response βσ changed smoothly through the conical intersection seam from -90º to +90º 50. In this work we shall consider how the sign and magnitude of βσ for the S1 electronic state varies through the CI, i.e. continuously and so approximately passing through βσ = 0.0º or discontinuously from βσ = -90.0º to βσ = +90.0º. 2.3 The REKS methodology

The SSR method used in this work for quantum chemical calculations is based on ensemble density functional theory and is capable of providing accurate account of both dynamic and non-dynamic electron correlation in the ground and excited states of strongly correlated molecular systems.26,27,51 The method was previously benchmarked in the calculation of the ground and excited states of molecules and proved to be capable of delivering the results matching the accuracy of high-level ab initio multi-reference methodologies.28,51 Importantly, the SSR method is capable of describing both avoided crossings and real crossings (conical intersections) between the ground and excited states of molecules; the traditional spin-conserving linear 8 ACS Paragon Plus Environment

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


response TD-DFT is lacking this ability.52,53 Importantly, the SSR method was validated in the calculation of conical intersection geometries of organic molecules, see refs. 28,51, where the SSR geometries matched the geometries from MRCISD within 0.06 A. In this work, we applied the SSR method in connection with the BH&HLYP density functional54–56 and the 6-31G* basis set57 to study the density variation along the NAMD trajectories of the F-NAIBP molecular motor. Two frontier KS orbitals, the highest occupied orbital and the lowest unoccupied orbital, were taken to the active space of the SSR method. In total, six exemplary trajectories were selected from the pool of NAMD TSH-OM2/MRCI simulations;16 all trajectories represent the EP→ZM photo-isomerization step. In ref. 16, a careful analysis of the dynamics characteristics of all the NAMD trajectories was carried out, whereby it was established that all the trajectories can be separated in two classes, fast and slow trajectories, differing by their S1 lifetime as well as the manner of approaching the respective conical intersection region. As was found in ref. 16, certain individual trajectories represent the major features of the two classes of dynamics in the best way, see Fig. 4 of ref. 16. The two trajectories shown in Fig. of ref. 16, a fast trajectory (F) where the S1→S0 surface hop occurs at 352 fs and a slow trajectory (S) with the hop time of 570 fs,16 were selected in this work for the QTAIM analysis. Additionally, two S trajectories with the hop times of 556 and 557 fs and two F trajectories with the hop times of 321 and 323 fs from ref. 16 were also taken; when discussing the results however we will refer primarily to the former two trajectories. The F trajectories represent fast photoisomerization process where, upon photoexcitation, the molecule in the S1 state quickly undergoes torsion about the central double bond. The S trajectories are characterized by considerably longer S1 lifetime, which, as was mentioned in the earlier publication16, was caused by a delayed departure from the Franck-Condon (FC) point. In dihedral torsion angle θ in an earlier publication16, the NAMD TSH-OM2/MRCI simulations were carried out with the time step of 0.05 fs for the overall duration of 800 fs. Of the 16000 snapshots for each of the trajectories, 400 snapshots (one in forty) were picked up for the SSR-BH&HLYP/6-31G* calculations and the QTAIM analysis; this corresponds to an effective time step of 2 fs. The S0 and S1 energies along the trajectories F and S and the S1-S0 energy gap are plotted in Figure S1-S3(a-d) of the Supplementary Materials. For the two trajectories shown in Fig. 4 of ref. 16, the S1-S0 energy gap obtained with SSR reaches its minimum at 358 fs (F) and 574 fs (S), in a good agreement with the hop times of 352 fs and 570 fs, respectively, obtained in the NAMD TSH-OM2/MRCI simulations. Hence, it can be assumed that, in the SSR simulations, the S1→S0 surface hop would be initiated at the minimal gap magnitude, that is at 352 fs (F) and 574 fs (S); for other trajectories the time tips of reaching the minimal gap are consistent with the TSH-OM2/MRCI hop times within ca. 10 fs. Note however that, due to the use of the geometry obtained by the OM2/MRCI semi-empirical method, the SSR S1-S0 energy gaps may not completely vanish at the time of hop. In the NAMD TSH formalism, the 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

surface hop is initiated when the non-adiabatic coupling between the states becomes sufficiently large; the latter becomes large even at a finite (but small) energy gap. The electron densities along the trajectories F and S were calculated using the orbitals and the orbital fractional occupation numbers (FONs) obtained in the SSR calculations. The occupation numbers of the SSR active orbitals are obtained from diagonalization of the density matrix constructed using the state-interaction coefficients of the SSR method and the FONs of the active orbitals of the REKS and ROKS states in the SAREKS functional, 26,27 see the Supporting Materials for more detail. The occupation numbers of the strongly occupied active orbital of SSR are shown in Figure S1-S3(e-h) of the Supplementary Materials. According to Figure S1-S3(e-h), the S0 state remains closed-shell along the whole duration of time (the occupation numbers close to 2 and 0) while the S1 state remains open-shell (the occupation numbers close to 1 and 1). The overall Mulliken charges on the indole moiety are shown in the lower panels of Figure S1-S3(e-h). The SSRBH&HLYP/6-31G* charges are generally consistent with the charges obtained by the OM2/MRCI method, see Fig. 5 of the Supporting Information to ref. 16. The QTAIM and stress tensor analysis was performed with the AIMAll58 suite. The SSR calculations were carried out using the COLOGNE2012 program59.

4. Results and discussions

Before starting the discussion of the results of the electronic density analysis, a brief word about the electronic structure of the F-NAIBP molecular motor seems appropriate. In the S0 state near the equilibrium conformation, the F-NAIBP motor features positively charged indole moiety, see Scheme 1. Breaking of the central double bond connecting the rotor and the stator blades of the motor occurs either by homolytic or heterolytic mechanism, where the former leads to the occurrence of a diradicaloid electronic structure (with positively charged indole) and the latter to a charge transfer (CT) electronic structure (neutral indole, positively charged indanylidene). When the central double bond is broken by torsion, the ground electronic state acquires pronounced CT character as illustrated by the total Mulliken charges on the indole moiety shown in Figure S12 of the Supporting Materials. The excited electronic state remains almost purely diradicaloid on the segment of the trajectory leading away from the FC point towards the S1 minimum and the conical intersection seam. Near the conical intersection an abrupt change of the character of the electronic states may occur and the S1 state begins to acquire discernible CT character as the molecule moves away from the conical intersection region towards the product conformation. Note however, that conical intersections are almost never reached in the NAMD simulations and the non-adiabatic population transfer usually occurs in their close proximity. As 10 ACS Paragon Plus Environment

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


illustrated by Figure S1-2 of the Supporting Materials, the S1 state of F-NAIBP retains the diradicaloid character (positive indole) for an extended period of time before the hop at 574 fs, when moving along the slow trajectory. As was speculated in ref. 16, there appear to be certain factors which hold down the progress of this type of trajectories towards the conical intersection region. The nature of these factors will be analyzed in the following sections 4.1 and 4.2.

4.1 The QTAIM and stress tensor properties at the axial C1-C2 BCP of fast and slow EP → ZM trajectories

In this section we present the results of the axial C1-C2 BCP that demonstrate differences in the fast (F) and slow (S) EP → ZM trajectories. In the next section, section 4.2, we attempt to explain the differences in the F and S trajectories in terms of the topological factors that either may facilitate or hinder the working of the FNAIBP motor. The definition of the dihedral torsion angle θ and the stress tensor response βσ as well as a sketch of the EP → ZM isomerization step of the F-NAIBP motor are provided in Scheme 1.

(a)

(b)

(c) Scheme 1. The definition of the dihedral torsion angle θ is provide in sub-figure (a), where the Cn labels correspond to the molecular graph atom numbering scheme in Figure 1. The construction of the stress tensor response βσ angle from the dot



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

product of the electronic stress tensor response vector e1σ and the projection vector Y is presented in sub-figure (b). Positive values stress tensor response βσ > +90.0º are shown in blue and are mapped back into the range -90.0º to 0.0º. Values of the response β are constructed similarly, from the electronic response vector e2 and the projection vector Y. The schematic for the EP → ZM sequence for the F-NAIBP motor that includes the location of the C1-C2 BCP, indicated by the green marker, is provided in sub-figure (c).

For the F trajectory the F-NAIBP motor proceeds directly towards the S1 minimum and the CI, see Figure S1(a, c,e) of the Supplementary Materials S1. Conversely, the S trajectory is characterized by an extended time period when the molecule remains in conformations near the FC point, see Figure S1(b,d,f). Although this behavior is not apparent from the variation of the relative energy ∆E with time in Figure S1(b,d), it is confirmed by inspecting the dihedral torsion angle θ in Fig 4 of θ in an earlier publication16. Examination of the variation of the ellipticity ε of the C1-C2 axial BCP with time for the F and S trajectories shows contrasting behavior; with the F trajectory in the S1 state reaching the minimum value of the ellipticity ε much faster than the S trajectory in the S1 state, see Figure 1. Each of the minimum values of the ellipticity ε coincides with the respective conical intersection seam (CI) at hop-times of 352 fs and 574 fs for the F and S trajectories, see Figure 1(a) and Figure 1(b), respectively. It can also be seen that during the photoisomerization reaction for each of the F and S trajectories the S1 state initially maintains a higher value of the ellipticity ε and hence more double bond character than the S0 state.

(a)

(b)


Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(c)

(d)

Figure 1.The variation of the ellipticity ε with time for the axial C1-C2 BCP for the fast (F) trajectory are given in subfigures (a) the corresponding values for the slow (S) trajectory are given in sub-figure (b). The variation of the ellipticity ε with the dihedral torsion angle θ for the fast (F) and (S) trajectories are given in sub-figures (c) and (d) respectively. See the C1-C2 BCP in Scheme1.

The maintenance of constant values of the ellipticity ε of the axial C1-C2 BCP for approximately the first 150 fs and 400 fs for the F and S trajectories respectively suggests that the conformation of the F-NAIBP motor is being stalled during these time periods. In section 4.2 we will examine the role that the intramolecular closedshell BCP bond-paths have in effectively stalling or holding the F-NAIBP motor. The stalling of the F-NAIBP motor is particularly apparent when examining the variation of ε of the axial C1-C2 BCP with the dihedral torsion angle θ, see Figure 1(c) and Figure 1(d), particularly for the S trajectory in the S1 state that is stalled at θ = 20.0º, see Figure 1(d). The spike is apparent at the CI for the variation of the ellipticity ε of the axial C1-C2 BCP with the dihedral torsion angle θ for the F trajectory in particular, although this is not related to a pyramidalization distortion but instead occurs as a consequence of the change of the character of the S1 state, from charge transfer to covalent, at the CI16. The variations of the difference S1-S0 of the stress tensor stiffness Sσ of the axial C1-C2 BCP with time for the F and S trajectory, see Figure 2(a) and Figure 2(b) is consistent with the variation of ellipticity ε with time from Figure 1(a) and Figure 1(b). The S1-S0 difference of stress tensor stiffness Sσ of the axial C1-C2 BCP is positive for the F and S trajectories approximately for the first 150 fs and 400 fs respectively, in agreement with the time period for which the ellipticity ε values remain at an approximately constant value. We again notice the stalled behavior of the F-NABIP motor at approximately θ = 20.0º for the S but not the F trajectory indicating the that motor is stalled at this orientation for about 400 fs for the S trajectory but not for the F trajectory, for the F and S trajectories see Figure 2(c) and Figure 2(d) respectively. 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

(c)

(d)

(c)

(d)

(e)

(f)

Figure 2.The variation of the difference (S1-S0) of the stress tensor stiffness Sσ with time for the fast (F) trajectory are given in sub-figure (a), the corresponding values for the slow (S) trajectory are given in sub-figure (b). The variation of the difference (S1-S0) stress tensor stiffness Sσ with the dihedral torsion angle θ for the fast (F) and slow (S) trajectories are given in sub-figures (c) and (d) respectively. See the caption of Figure 1 for further details. The variation of stress tensor eigenvalue associated with the bond-path; λ3σ with time of the axial C1-C2 BCP for the F and S trajectories are


Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


given in sub-figures (e) and (f) respectively. The stiffness S which does not vary in accordance with physical intuition is provided in Supplementary Materials S4.

The variation of stress tensor eigenvalue λ3σ of the axial C1-C2 BCP, which is associated with the bond-path, with time shows a phase-transition like behavior at the hop-time for both the F and S trajectory, see Figure 2(ef) respectively. For the duration of the F and S trajectories there is no discernable difference between the variation with time of λ3σ for the S0 and S1 states. The stress tensor response βσ also responds consistently to the changes in the PES in as far as the hop-times coincide with the turning points of βσ both in the S1 and S0 states for the F and S trajectories, see Figure 3(a-d).

(a)

(b)

(c)

(d)

Figure 3.The variation of the stress tensor response βσ with time and the variation of βσ of the axial C1-C2 BCP with dihedral torsion angle θ for the F trajectory, S1 in red and S0 in black, are given in sub-figures (a) and (c) respectively, the corresponding values for the S trajectory are given in sub-figures (b) and (d) respectively. See the caption of Figure 1 for



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

further details. Data points for βσ are greyed out S0 before the CI seam and presented in pink S1 after the CI seam since the F-NAIBP motor follows S1 before the CI and S0 afterwards.

We notice that for the variation of βσ of the axial C1-C2 BCP with the dihedral torsion angle θ, for the S trajectory, the F-NAIBP motor is stalled at βσ ≈ 20.0º for both the S0 and S1 state consistent with the ellipticity ε and stress tensor stiffness Sσ, see Figure 1(d) and Figure 2(d) respectively. As was seen earlier50 the values of the stress tensor response βσ in the S1 state pass through βσ = 0.0º in the vicinity of the CI seam but the values of βσ for S0 do not, see the theory section 2.2 for further discussions and Figure 3. The stress tensor response βσ of the axial C1-C2 BCP for the S1 state of the F-NAIBP motor continues to react to the dihedral torsion θ in the vicinity of the CI but βσ in the S0 state does not. The reaction of the S1 stress tensor response βσ to the dihedral torsion θ resembles the form of a dielectric constant for a damped oscillator passing through a resonance. In this model the damping would be caused by the action of the intramolecular closed-shell BCPs and their associated bond-paths hindering the dihedral torsion θ of the F-NAIBP motor. From a simple scalar perspective the axial C1-C2 BCP transforms from a double bond to a single bond in both the S1 and S0 states, however from a vector-based perspective there are differences in the mechanism of the bond-path deformation and the character of the axial C1-C2 BCP in the vicinity of the CI. The step-like profile of the variation of S0 stress tensor response βσ through the CI suggests that that the axial C1-C2 BCP bond-path breaks completely and therefore does not respond to the application of the dihedral torsion θ. The non-step-like continuation of the values of S1 stress tensor response βσ of this BCP through βσ ≈ 0.0º in the vicinity of the CI suggests that the bond-path is deformed since there is still a response to the applied dihedral torsion θ in the vicinity of the CI. The nature of this deformation of the axial C1-C2 BCP bond-path in the S1 state could be a twisting that could compress both the λ1σ and λ2σ and so still give rise to a very low ellipticity ε value.

4.2 The influence of the closed-shell BCPs on the stress tensor stiffness Sσ of the F-NAIBP C1-C2 axial BCP In this section, we attempt to determine the origin of the increased stress tensor stiffness Sσ of the axial C1-C2 BCP observed in section 4.1 for the fast (F) and slow (S) trajectories. In doing so we will investigate the role that the closed-shell BCPs and their associated bond-paths have in affecting the stress tensor stiffness Sσ of the axial C1-C2 BCP.


Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

Figure 4. The molecular graphs with lowest and highest topologically complexities (Σbrc) of the fast (F) trajectory are shown in sub-figure (a) at the hop-time = 352 fs and sub-figure (b) at time = 568 fs respectively. The corresponding molecular graphs of the slow (S) trajectory show the lowest Σbrc value at the hop-time = 574 fs in sub-figure (c) and the highest Σbrc value at time = 456 fs in sub-figure (d). Rotation is about the bond-path of the C1-C2 dihedral torsional (axle). The green, red and blue undecorated spheres represent the bond critical points (BCPs), ring critical points (RCPs) and cage critical points (CCPs) respectively. The dashed thick lines represent the bond-paths of the closed-shell BCPs and the red arrows indicate the additional closed-shell BCPs and bond-paths.

The molecular graphs possessing the lowest and highest values of the topological complexity Σbrc show a considerable diversity in the intramolecular connectivity, see Figure 4. In particular, for the F trajectory at the CI seam there are no intramolecular bond-paths and BCPs that separate the stator and the rotor blades of the F17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

NAIBP motor other than the axial C1-C2 BCP, see Figure 4(a). This topology, with the lowest possible topological complexity Σbrc is unique to the F trajectory and only occurs at the conical intersection seam at the hop-time = 352 fs. This result is demonstrated by the plot of variation with time of the number of BCPs (b) for the F trajectory, where the F-NAIBP topology without any intermolecular BCPs occurs at (τ = 352 fs, b = 44), indicated by the blue arrow in Figure 5(a). The corresponding variation of the numbers of RCPs (r) also shows the minimum possible value of r for the F-NAIBP motor, see Supplementary Materials S6. Examination of the variation of the numbers of BCPs (b) with time for the S trajectory maintains significantly more intramolecular connectivity, i.e. higher values of b up to the hop-time at 574 fs, indicated by the magenta arrow, than for the F trajectory, see Figure 5(b). The corresponding variations of the numbers of BCPs (b) with dihedral torsion angle θ for the S and F trajectories given in Figure 5(c) and Figure 5(d) respectively.

(a)

(b)

(c)

(d)


Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(e)

(f)

Figure 5. The plots of the variation of the numbers of BCPs (b) with time for the fast (F) and slow (S) trajectories are presented in sub-figures (a) and (b) respectively and corresponding plots for the variation of BCPs (b) with the dihedral torsion angle θ are given in sub-figures (c) and (d). Black squares and red crosses indicate the S0 and S1 states respectively, hop-times are indicated by blue and magenta arrows for the F and S trajectories respectively. The 2-D quantum topology phase diagram (QTPD) of the solutions of the Poincaré-Hopf relation of the F-NAIBP molecular motor of the F and S trajectories, hop-times in fs are shown in sub-figures (e) and (f) respectively. The number of BCPs (b) and RCPs (r) plotted on the x- and y-axes, respectively. The upper left regions marked with orange asterisks denote the region of topological instability based on there being increasingly more CCPs. The topologically stable isomers found are represented as black squares and red circles for the S0 and S1 states respectively. The green squares denote the possible but ‘missing’ solutions of the Poincaré-Hopf relation. The blue crosses represent the forbidden solutions of the Poincaré-Hopf relation; see the main text for further details. Times (initial and hop-time) are indicated in black and red for the S0 and S1 states respectively.

The QTPD of the F trajectory is more diverse than that of the S trajectory; there are ten topologies for F compared with only eight for the S trajectory, see Figure 5(e) and Figure 5(f) respectively. For the F trajectory there are six 2-DQT and four 3-DQT quantum topologies, whereas for the S trajectory there are four 2-DQT and four 3-DQT quantum topologies. This indicates that the intramolecular BCPs and associated bond-paths are not favored in the functioning of the motor. Additionally, the F trajectory possesses the least (b = 44, r = 4) and most (b = 49, r = 10) topologically complex molecular graphs respectively, see Figure 5(e). The most topological complex molecular graph for the F trajectory (b = 49, r = 10) occurs close to the same time, at approximately 570 fs, as the least topologically complex graph of the S trajectory (b = 45, r = 5), see Figure 5(f). The longest lasting topology corresponds to (b = 47, r = 8) for the S trajectory. These findings suggest a direct link between the numbers of BCPs and RCPs and the readiness to reach the CI seam, i.e. fewer BCPs facilitate the approach to the CI seam and more BCPs inhibit the approach. It is clear that the closed-shell BCPs are responsible for differences in the results seen in Figures 1-5. These results are consistent with the values of 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

the ellipticity ε and stress tensor stiffness Sσ of the axial C1-C2 BCP in that lower values of (b, r) correlate with lower values of both the ellipticity ε and the stress tensor stiffness Sσ, see Figure 1 and Figure 2 respectively. The chemical behavior of the F13--C9 BCP is very unusual; because the F13--C9 BCP oscillates between existing as the shared-shell F13-C9 BCP and the closed-shell F13--C9 BCP where the closed-shell version has consistently more negative values for H(rb) than does the shared-shell F13-C9 BCP, see Figure 6(a-b). This behavior does not occur for the corresponding shared-shell H13-C9 BCP in the H-NAIBP motor, see Supplementary Materials S7. The changing chemical character of the F13--C9 BCP indicates that the F13--C9 BCP is actively chemically coupling with the intramolecular closed-shell BCPs in a two-way exchange of chemical character, the F13--C9 BCP acquires closed-shell BCP character donated from the F13--H36 BCP. This chemical coupling to the intramolecular closed-shell BCPs generates values of H(rb) < 0 in the intramolecular BCPs that will effectively strengthen the intramolecular closed-shell BCPs. If the chemical coupling between the F13-C9 BCP and the intramolecular closed-shell BCPs occurs for a sustained time period then the axial C1-C2 BCP will be more difficult to move as determined by the stress tensor stiffness Sσ difference S1-S0 of this BCP, see Figure 2. This suggests that a stronger coupling of the F13--C9 BCP with the closed-shell intramolecular BCPs is present for the S trajectory than for the F trajectory. These differences in coupling can account for the very different sets of significant and continuous or near continuous closed-shell BCPs and bond-paths for the F and S trajectories. It should be noted that smaller differences between the S1 and S0 ellipticity ε values for the S trajectory before the hop time result in less negative values of H(rb) for the intramolecular BCPs close to the F13 BCP, see Supplementary Materials S10 and Supplementary Materials S11. Minor hindrance for the F trajectory caused by the H---H BCPs and H--C BCPs with H(rb) < 0 is apparent by the slight clustering of data points for the ellipticity ε and the stress tensor stiffness Sσ at θ ≈ 0.0º and for θ ≈ 30.0º to θ ≈ 40.0º, at θ ≈ 60.0º and 180.0º see Figure 1(c) and Figure 2(c). This can be compared with the considerable clustering of data points of the S trajectory at θ ≈ 20.0º for the ellipticity ε, the stress tensor stiffness Sσ and the stress tensor response βσ, see Figure 1(d), Figure 2(d) and Figure 3(d). Note, that the F trajectory is discontinuous at the hop-time of 352 fs because there are no closed-shell BCPs as was mentioned earlier.


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

(e)

(f)

Figure 6. The variation of the local total energy density H(rb) with time for the closed-shell BCPs with H(rb) < 0 for the for the fast (F) trajectory in sub-figures (a), (c) and (e) the H(rb) values for the slow (S) trajectory results are presented in sub-figures (b), (d) and (f), all plots correspond to the S1 electronic state. The C--C BCPs, H--C BCPs are indicated in red, the H---H BCPs, are indicated in black, the F--C BCPs and F--H BCPs are indicated in magenta and blue respectively in the electronic version, see the caption of Figure 1 for further details. The plots for the S0 excited state are virtually indistinguishable and are provided in Supplementary Materials S5.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

Examination of the variation of H(rb) with the dihedral torsion angle θ for the closed shell BCPs also indicates that for the S trajectory the F-NAIBP motor is stalled at θ ≈ 20.0º, evident by the large concentration of data points; this behavior not being apparent for the F trajectory, see Figure 7(a), Figure 7(b) and Figure 7(c). The fact that the F-NAIBP motor is stalled at θ ≈ 20.0º indicates the effect that the increased chemical coupling between the backbone F13--C9 BCP and the intramolecular F13--H36 BCP bond-path has on holding in place the F-NAIBP motor during the S trajectory. This is due to the large negative values of H(rb) of the backbone F13--C9 BCP coupling to the intramolecular F13--H36 BCP. The coupling behavior strengthens the F13--H36 BCP that effectively binds the rotor and stator blades of the F-NAIBP motor during the first 450 fs of the S trajectory. This stalling effect does not occur to the same extent for either the F or S trajectory of the H-NAIBP version16 of the motor as seen by the total local energy density H(rb) values of the intramolecular BCPs, see Supplementary materials S7. For the F trajectory of the H-NAIBP motor there are very large fluctuations in the H(rb) value with time before the hop-time which decreases after the hop-time. For the S trajectory the fluctuations in the H(rb) values are much smaller before the hop-time than for the F trajectory. There are however, no significant instances of H(rb) < 0 for either trajectory before the hop times.

(a)

(b)


Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(c) (c) Figure 7. The variation of the local total energy density H(rb) with dihedral torsion angle θ for the closed-shell BCPs with H(rb) < 0 for the for the fast (F) and slow (S) trajectory results are presented in sub-figures (a), (b) and (c), all plots correspond to the S1 electronic state. See the caption of Figure 7 for further details. The corresponding values for the S0 state are provided in Supplementary Materials S5.

Conclusions

Here, we report on the first investigation that combines non-adiabatic molecular dynamics simulations with QTAIM and the stress tensor analysis. This was done to address the origin of differences between the fast (F) and the slow (S) dynamics trajectories of the F-NAIBP motor, which remains near the FC point for more than 400 fs when following the slow trajectory while rapidly sliding down the S1 PES when following the fast trajectory. A scalar and vector-based QTAIM and stress tensor analyses was carried out to quantify the variation of the axial C1-C2 BCP in the S0 and S1 state along both F and S trajectories during the functioning of the FNAIBP motor. In particular, we considered the ellipticity ε, the stress tensor stiffness Sσ and the stress tensor eigenvalue λ3σ associated with the bond-path and the stress tensor response βσ. We found that, along the S trajectory, the axial C1-C2 BCP remained unusually stiff and resistant to the dihedral torsion as measured by the ellipticity ε and the stress tensor stiffness Sσ. The stress tensor eigenvalue λ3σ associated with the bond-path and the stress tensor response βσ were both found to clearly demonstrate phase-transition-like behaviors at the hoptimes for the F and S trajectories. The stress tensor response βσ suggested a different mechanism of the bondpath deformation and character of the axial C1-C2 BCP for the S0 state compared with the S1 state. We observed that the form of the variation of βσ of the axial C1-C2 BCP in the S1 state to the dihedral torsion θ of the FNAIBP motor resembles the form of a dielectric constant fo r a damped oscillator passing through a resonance; 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32

the damping being attributed to the action of the intramolecular closed-shell BCPs and bond-paths hindering the dihedral torsion θ of the F-NAIBP motor. The effect of the intramolecular interactions between the rotor and the stator blades of the motor and the atoms added at the chiral center (F13 or H13) of the motor on the axial C1-C2 BCP and associated bond-path was investigated by analyzing the quantum topology phase diagrams (QTPDs) and the dependence of the number of BCPs (b) and RCPs (r) along the F and S trajectories as functions of time and dihedral torsion θ. It was found that the numbers of BCPs (b) and RCPs (r) was minimal at the hop-time for both trajectories. Along the F trajectory, the only intramolecular interaction between the blades of the F-NAIBP motor at this time was the axial C1-C2 BCP. From the results of the quantum topology analysis it was inferred that increasing number of intramolecular interactions before the hop time would hinder the progress of the motor. Examination of the bonding that revealed a very unusual phenomenon; the F13-C9/F13--C9 BCP oscillated between existing as a shared-shell BCP and a closed-shell BCP. Most unusual was the observation that the closed-shell F13--C9 BCP possessed more negative values of H(rb) than the shared-shell version. This transformation of chemical character between shared-shell BCP and closed-shell BCP occurred much more frequently for the S trajectory than for the F trajectory. Furthermore, the intramolecular F13--H36 BCP was continuously coupled with the F13-H36 BCP for the first 400 fs of the S trajectory; during which period the motor was stalled near the FC point. No such change of chemical character was observed for the shared-shell H13-C9 BCP of the H-NAIBP, that remained shared-shell in character for the duration of both the F and S trajectories. This effect for the S trajectory was observed to be dependent on the existence of sufficiently large differences in the S1 and S0 ellipticity ε values before the hop time, smaller differences resulted in less negative values of H(rb) or even positive values for the intramolecular BCPs close to the F13 BCP, see Supplementary Materials S8 and Supplementary Materials S9. The observed behavior of the C1-C2 BCP in the S1 state of the F-NAIBP motor, as well as the “sticky” intramolecular interactions between the F13 NCP of the rotor and the stator blades, not only offered an explanation for the differences in the dynamics of the motor along the F and the S NAMD trajectories, but also suggests a possible way of tuning the photochemical behavior of light-driven motors and switches. Specifically, increasing/decreasing the numbers of intramolecular BCPs and affecting the presence of chemical coupling between the moving blades of the motor by changing the atoms added to the chiral centre (i.e. F13, H13), as revealed by the QTPDs and the analysis of the total local energy density H(rb), can result in speeding up/slowing down the torsional motion occurring on the S1 PES of the motor and consequently the S1 lifetime. It remains however, to be seen how strong an effect the intermolecular interactions measurable by the presence of values of H(rb) < 0, e.g. with the solvent where the motor is dissolved or with the substrate on which it is 24 ACS Paragon Plus Environment

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


anchored, may have on the detected intramolecular interactions; especially if the intermolecular interactions will strengthen or weaken the intramolecular ones.

ASSOCIATED CONTENT Supporting Information Total energies of the S1 and S0 states, S1-S0 energy gap, fractional orbital occupations, and Mulliken charges on the indole moiety of the F-NAIBP motor along the fast (F) and slow (S) trajectories for three such pairs of trajectories. The REKS equations used to construct the S0 and the S1 densities, variation of the difference S1-S0 stiffness S and with time of the F-NAIBP motor. Local total energy density H(rb) of the F-NAIBP motor with time for the closed-shell BCPs with H(rb) < 0 for the S0 electronic state and the corresponding plots for the variation of the local total energy density H(rb) with dihedral angle θ. Numbers of RCPs (r) with dihedral angle θ for the F and S trajectories of the F-NAIBP motor. Local total energy density H(rb) of the F-NAIBP motor with time for the closed-shell BCPs with H(rb) > 0 for the S0 and S1 electronic states. Local total energy density H(rb) with time for the closed-shell BCPs with H(rb) < 0 for the S1 electronic state and the corresponding plots for the variation of the local total energy density H(rb) with dihedral angle θ of the H-NAIBP motor. The two additional data sets of F and S trajectories for ellipticity ε, stress tensor eigenvalue λ3σ variation of numbers of BCPs and total local energy density H(rb). The Supporting Information is available free of charge on the ACS Publications website at DOI: xx.xxx/ acs.jpca. xxxxxxx. Acknowledgements The National Natural Science Foundation of China is also acknowledged, project approval numbers: 21273069 and 21673071.The One Hundred Talents Foundation of Hunan Province and the aid program for the Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan Province are gratefully acknowledged for the support of S.J. and S.R.K. 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 32

References (1)

Browne, W. R.; Feringa, B. L. Making Molecular Machines Work. Nat. Nanotechnol. 2006, 1, 25–35.

(2)

Balzani, V.; Credi, A.; Venturi, M. Light Powered Molecular Machines. Chem. Soc. Rev. 2009, 38, 1542–1550.

(3)

Coskun, A.; Banaszak, M.; Astumian, R. D.; Stoddart, J. F.; Grzybowski, B. A. Great Expectations: Can Artificial Molecular Machines Deliver on Their Promise? Chem. Soc. Rev. 2011, 41, 19–30.

(4)

Koumura, N.; Zijlstra, R. W. J.; van Delden, R. A.; Harada, N.; Feringa, B. L. Light-Driven Monodirectional Molecular Rotor. Nature 1999, 401, 152–155.

(5)

Koumura, N.; Geertsema, E. M.; van Gelder, M. B.; Meetsma, A.; Feringa, B. L. Second Generation Light-Driven Molecular Motors. Unidirectional Rotation Controlled by a Single Stereogenic Center with Near-Perfect Photoequilibria and Acceleration of the Speed of Rotation by Structural Modification. J. Am. Chem. Soc. 2002, 124, 5037–5051.

(6)

Pijper, D.; van Delden, R. A.; Meetsma, A.; Feringa, B. L. Acceleration of a Nanomotor: Electronic Control of the Rotary Speed of a Light-Driven Molecular Rotor. J. Am. Chem. Soc. 2005, 127, 17612– 17613.

(7)

Pollard, M. M.; Meetsma, A.; Feringa, B. L. A Redesign of Light-Driven Rotary Molecular Motors. Org. Biomol. Chem. 2008, 6, 507–512.

(8)

Fernández Landaluce, T.; London, G.; Pollard, M. M.; Rudolf, P.; Feringa, B. L. Rotary Molecular Motors: A Large Increase in Speed through a Small Change in Design. J. Org. Chem. 2010, 75, 5323– 5325.

(9)

Feringa, B. L. The Art of Building Small: From Molecular Switches to Molecular Motors. J. Org. Chem. 2007, 72, 6635–6652.


Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Kazaryan, A.; Kistemaker, J. C. M.; Schäfer, L. V.; Browne, W. R.; Feringa, B. L.; Filatov, M. Understanding the Dynamics Behind the Photoisomerization of a Light-Driven Fluorene Molecular Rotary Motor. J. Phys. Chem. A 2010, 114, 5058–5067. (11) Kazaryan, A.; Lan, Z.; Schäfer, L. V.; Thiel, W.; Filatov, M. Surface Hopping Excited-State Dynamics Study of the Photoisomerization of a Light-Driven Fluorene Molecular Rotary Motor. J. Chem. Theory Comput. 2011, 7, 2189–2199. (12) Conyard, J.; Addison, K.; Heisler, I. A.; Cnossen, A.; Browne, W. R.; Feringa, B. L.; Meech, S. R. Ultrafast Dynamics in the Power Stroke of a Molecular Rotary Motor. Nat. Chem. 2012, 4, 547–551. (13) Filatov, M.; Olivucci, M. Designing Conical Intersections for Light-Driven Single Molecule Rotary Motors: From Precessional to Axial Motion. J. Org. Chem. 2014, 79, 3587–3600. (14) Gozem, S.; Melaccio, F.; Luk, H. L.; Rinaldi, S.; Olivucci, M. Learning from Photobiology How to Design Molecular Devices Using a Computer. Chem. Soc. Rev. 2014, 43, 4019–4036. (15) Marchand, G., J.; Schapiro, I.; Valentini, A.; Frutos, L. M.; Pieri, E.; Olivucci, M., Léonard, J.; Gindensperger, E. Directionality of Double-Bond Photoisomerization Dynamics Induced by a Single Stereogenic Center. J. Phys. Chem. Lett. 2015, 6, 599–604. (16) Nikiforov, A.; Gamez, J. A.; Thiel, W.; Filatov, M. Computational Design of a Family of Light-Driven Rotary Molecular Motors with Improved Quantum Efficiency. J. Phys. Chem. Lett. 2016, 7, 105–110. (17) Bader, R. F. Atoms in Molecules: A Quantum Theory; Clarendon press: Oxford, U. K., 1990. (18) Popelier, P. L. Atoms in Molecules, Prentice Hall: New Jersey, USA., 1999. (19) Guo, H.; Morales-Bayuelo, A.; Xu, T. L.; Momen, R.; Wang, L. L.; Yang, P.; Kirk, S. R.; Jenkins, S. Distinguishing and Quantifying the Torquoselectivity in Competitive Ring-Opening Reactions Using the Stress Tensor and QTAIM. J. Comput. Chem. 2016, 37, 2722–2733. (20) Schweighauser, L.; Strauss, M. A.; Bellotto, S.; Wegner, H. A. Attraction or Repulsion? London Dispersion Forces Control Azobenzene Switches. Angew. Chem. Int. Ed. 2015, 54, 13436–13439. 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 32

(21) Weber, W.; Thiel, W. Orthogonalization Corrections for Semiempirical Methods. Theor. Chem. Acc. 2000, 103, 495–506. (22) Koslowski, A.; Beck, M. E.; Thiel, W. Implementation of a General Multireference Configuration Interaction Procedure with Analytic Gradients in a Semiempirical Context Using the Graphical Unitary Group Approach. J. Comput. Chem. 2003, 24, 714–726. (23) Thiel, W. Semiempirical Quantum–Chemical Methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 145–157. (24) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061–1071. (25) Barbatti, M. Nonadiabatic Dynamics with Trajectory Surface Hopping Method. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 620–633. (26) Filatov, M. Spin-Restricted Ensemble-Referenced Kohn–Sham Method: Basic Principles and Application to Strongly Correlated Ground and Excited States of Molecules. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2015, 5, 146–167. (27) Filatov, M. Ensemble DFT Approach to Excited States of Strongly Correlated Molecular Systems. In Density-Functional Methods for Excited States; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Springer: Heidelberg, 2015; Vol. 368, pp 97–124. (28) Nikiforov, A.; Gamez, J. A.; Thiel, W.; Huix-Rotllant, M.; Filatov, M. Assessment of Approximate Computational Methods for Conical Intersections and Branching Plane Vectors in Organic Molecules. J. Chem. Phys. 2014, 141, 124122. (29) Ayers, P. W.; Jenkins, S. An Electron-Preceding Perspective on the Deformation of Materials. J. Chem. Phys. 2009, 130, 154104. (30) Jenkins, S.; Kirk, S. R.; Cote, A. S.; Ross, D. K.; Morrison, I. Dependence of the Normal Modes on the Electronic Structure of Various Phases of Ice as Calculated by Ab Initio Methods. Can. J. Phys. 2003, 81, 225–231. 28 ACS Paragon Plus Environment

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) Bone, R. G. A.; Bader, R. F. W. Identifying and Analyzing Intermolecular Bonding Interactions in van Der Waals Molecules†. J. Phys. Chem. 1996, 100, 10892–10911. (32) Jenkins, S.; Maza, J. R.; Xu, T. L. ; Dong, J. J.; Kirk, S. R. Biphenyl: A Stress Tensor and Vector-Based Perspective Explored within the Quantum Theory of Atoms in Molecules. Int. J. Quantum Chem. 2015, 115, 1678–1690. (33) Cremer, D.; Kraka, E. A Description of the Chemical Bond in Terms of Local Properties of Electron Density and Energy. Croat. Chem. Acta 1984, 57, 1265–1287. (34) Jenkins, S.; Morrison, I. The Chemical Character of the Intermolecular Bonds of Seven Phases of Ice as Revealed by Ab Initio Calculation of Electron Densities. Chem. Phys. Lett. 2000, 317, 97–102. (35) Isaacs, E. D.; Shukla, A.; Platzman, P. M; Hamann, D. R.; Barbiellini, B.; Tulk, C. A. Covalency of the Hydrogen Bond in Ice: A Direct .X-Ray Measurement. Phys. Rev. Lett. 1999, 82, 600–603. (36) Xu, T. L.; Farrell, J.; Momen, R.; Azizi, A.; Kirk, S. R.; Jenkins, S.; Wales, D. J. A Stress Tensor Eigenvector Projection Space for the (H2O)5 Potential Energy Surface. Chem. Phys. Lett. 2017, 667, 25– 31. (37) Matta, C. F.; Hernández-Trujillo, J.; Tang, T.-H.; Bader, R. F. W. Hydrogen–Hydrogen Bonding: A Stabilizing Interaction in Molecules and Crystals. Chem. – Eur. J. 2003, 9, 1940–1951. (38) Hernández-Trujillo, J.; Matta, C. F. Hydrogen–Hydrogen Bonding in Biphenyl Revisited. Struct. Chem. 2007, 18, 849–857. (39) Castillo, N.; Matta, C. F.; Boyd, R. J. Fluorine−Fluorine Spin−Spin Coupling Constants in Aromatic Compounds: Correlations with the Delocalization Index and with the Internuclear Separation. J. Chem. Inf. Model. 2005, 45, 354–359. (40) Matta, C. F.; Sadjadi, S. A.; Braden, D. A.; Frenking, G. The Barrier to the Methyl Rotation in Cis-2Butene and Its Isomerization Energy to Trans-2-Butene, Revisited. J. Comput. Chem. 2016, 37 (1), 143– 154. 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32

(41) Bader, R. F. W. Quantum Topology of Molecular Charge Distributions. III. The Mechanics of an Atom in a Molecule. J. Chem. Phys. 1980, 73, 2871–2883. (42) Jenkins, S.; Blancafort, L.; Kirk, S. R.; Bearpark, M. J. The Response of the Electronic Structure to Electronic Excitation and Double Bond Torsion in Fulvene: A Combined QTAIM, Stress Tensor and MO Perspective. Phys. Chem. Chem. Phys. 2014, 16, 7115–7126. (43) Collard, K.; Hall, G. G. Orthogonal Trajectories of the Electron Density. Int. J. Quantum Chem. 1977, 12 , 623–637. (44) Johnson, C. K. Peaks, Passes, Pales and Pits: A Tour through the Critical Points of Interest in Density Maps, Abstract of paper ,JQ6. American Crystallographic Association Meeting, Asilomar, CA, Feb 21 25 1977; American Crystallographic Association: New York, 1977; PMSE 30. (45) Smith, V. H.; Price, P. F.; Absar, I. XVI. Representations of the Electron Density and Its Topographical Features. Isr. J. Chem. 1977, 16, 187–197. (46) Jenkins, S.; Rong, C.; Kirk, S. R.; Yin, D.; Liu, S. Spanning Set of Silica Cluster Isomer Topologies from QTAIM. J. Phys. Chem. A 2011, 115, 12503–12511. (47) Jenkins, S.; Restrepo, A.; David, J.; Yin, D.; Kirk, S. R. Spanning QTAIM Topology Phase Diagrams of Water Isomers W4, W5 and W6. Phys. Chem. Chem. Phys. 2011, 13, 11644–11656. (48) Jenkins, S. Quantum Topology Phase Diagrams for Molecules, Clusters, and Solids. Int. J. Quantum Chem. 2013, 113, 1603–1608. (49) Hu, M. X.; Xu, T. L.; Momen, R.; Guo, H.; Kirk, S. R.; Jenkins, S., Kim, K. S.; Filatov, M. A QTAIM and Stress Tensor Investigation of the Torsion Path of a Light-Driven Fluorene Molecular Rotary Motor. J. Comput. Chem. 2016, 37, 2588–2596. (50) Maza, J. R.; Jenkins, S.; Kirk, S. R. 11-Cis Retinal Torsion: A QTAIM and Stress Tensor Analysis of the S1 Excited State. Chem. Phys. Lett. 2016, 652, 112–116.


Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(51) Filatov, M. Assessment of Density Functional Methods for Obtaining Geometries at Conical Intersections in Organic Molecules. J. Chem. Theory Comput. 2013, 9, 4526–4541. (52) Levine, B. G.; Ko, C.; Quenneville, J.; MartÍnez, T. J. Conical Intersections and Double Excitations in Time-Dependent Density Functional Theory. Mol. Phys. 2006, 104, 1039–1051. (53) Huix-Rotllant, M.; Nikiforov, A.; Thiel, W.; Filatov, M. Description of Conical Intersections with Density Functional Methods. In Density-Functional Methods for Excited States; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Springer: Springer: Heidelberg, 2015; Vol. 368, pp 445–476. (54) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098–3100. (55) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. (56) Becke, A. D. A New Mixing of Hartree–Fock and Local Density‐functional Theories. J. Chem. Phys. 1993, 98, 1372–1377. (57) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self‐Consistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72, 650–654. (58) Keith, T. A. AIMAll, Revision 12.06.21; TK Gristmill Software: Overland Park KS, USA, 2012. (59) Kraka, E.; Filatov, M.; Zou, W.; Gräfenstein, J.; Izotov, D.; Gauss, J.; He, Y.; Wu, A.; Polo, V.; Olsson, L., et al. COLOGNE2012; Southern Methodist University: Dallas, TX, 2012.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic


Page 32 of 32

QTAIM and Stress Tensor Characterization of Intramolecular

Recommend Documents