QTAIM and Stress Tensor Characterization of ... - ACS Publications

Article pubs.acs.org/JPCA

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*

J. Phys. Chem. A 2017.121:4778-4792. Downloaded from pubs.acs.org by UNIV OF SUNDERLAND on 10/03/18. For personal use only.

Key Laboratory of Chemical Biology and Traditional Chinese Medicine Research and Key Laboratory of Resource Fine-Processing and Advanced Materials of Hunan Province of MOE, College of Chemistry and Chemical Engineering, Hunan Normal University, Changsha, Hunan 410081, China S Supporting Information *

ABSTRACT: A quantum theory of atoms in molecules (QTAIM) and stress tensor analysis was applied to analyze intramolecular interactions influencing the photoisomerization dynamics of a light-driven rotary molecular motor. For selected nonadiabatic 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.

1. INTRODUCTION

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 orientation;4,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.2,11,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 intersection.10,11,13 In the literature, several alternative

Light driven rotary molecular motors represent a novel class of molecules capable of performing unidirectional rotary motion as a consequence of photoisomerization about a double bond.1−9 The first models of molecular motors were based on overcrowded alkene framework,4−8 where the double bond connecting two blades, the rotor and the stator, was strained by the steric repulsion. The photoexcitation 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 photoisomerization proceeds further by nonadiabatic population transfer through conical intersection as was unequivocally demonstrated theoretically10,11 as well as experimentally.12 Having reached the ground state metastable conformation the © 2017 American Chemical Society

Received: March 13, 2017 Revised: June 5, 2017 Published: June 6, 2017 4778

DOI: 10.1021/acs.jpca.7b02347 J. Phys. Chem. A 2017, 121, 4778−4792

Article

The Journal of Physical Chemistry A Scheme 1a

Definition of the dihedral torsion angle θ is provided in subfigure (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 product of the electronic stress tensor response vector e1σ and the projection vector Y is presented in subfigure (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 subfigure (c). a

molecular motor designs were proposed where the necessity to pyramidalize was lifted.13−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 analogues;13 these conjectures were later corroborated by the nonadiabatic molecular dynamics (NAMD) simulations.16 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 photoisomerization 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 before.16 Particularly, two types of excited state dynamics of the 3-[(2S)-2-fluoro-2-methyl-1-indanylidene]-1-methyl-2-methylindole (F-NAIBP) motor were identified;16 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 rotation after the molecule spent prolonged period of time flexing near the Franck−Condon

(FC) point upon the photoexcitation. 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 2, 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 bond critical points (BCPs) are often formed and annihilated along torsion coordinates thus enabling the switching of bondpath connectivity within a nuclear skeleton. As suggested in recent studies,20 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 structure.19 4779


Article

The Journal of Physical Chemistry A

(RCP), respectively], and (3, 3) [the cage critical points (CCP)]. The eigenvector e3 indicates the direction of the bond-path at the BCP. The most and least preferred directions of electron accumulation are e2 and e1, respectively.29−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 ρ(r b) 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-path.32 A degree of covalent character is indicated for H(rb) < 0,33,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 in ice Ih to provide a quantitative explanation of the unusually high strength of ice Ih hydrogen bonds.35 This result was explained in terms of coupling, i.e., cooperativity between the hydrogen-bonds and adjoining (σ) O−H bondpaths that share an H NCP causing covalent bond character to “leak” from the O−H BCPs to the hydrogen-bond BCPs.34 More recently we demonstrated a two-way exchange of chemical character; covalent character was donated from a σ O−H BCP to the hydrogen bond H--O BCP that shared a proton and hydrogen bond character donated back to the σ O−H BCP from the hydrogen bond H--O BCP.36 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 number 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 BCP,37,38 F-H BCP,39 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 closedshell 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 rearrangements.32 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 π-bond.42 Previously, it was found that the stress tensor stiffness, Sσ, produced results that were in line with physical intuition.32,42 We can represent the QTAIM topologies of sets of isomers of allowed, forbidden, and unstable solutions to the Poincaré−Hopf relationship43−45 using the quantum topology phase diagram (QTPD).46−48 For molecules and clusters, the relation is expressed as

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 publication.16 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 in ref 19. 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 publication:16 fast (F) and slow (S) trajectories of dihedral torsion angle θ shown in Figure 4 in an earlier publication.16 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 formalism.24,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 stateaveraged spin-restricted ensemble-referenced Kohn−Sham (SISA-REKS or SSR, for brevity) method,26,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 other;16,28 therefore, the OM2/MRCI geometries can be used with confidence for the SSR calculations. 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 nonzero 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

n−b+r−c=1 4780

(1) DOI: 10.1021/acs.jpca.7b02347 J. Phys. Chem. A 2017, 121, 4778−4792

Article


Figure 1. Variation of the ellipticity ε with time for the axial C1−C2 BCP for the fast (F) trajectory is given in subfigure (a), and the corresponding values for the slow (S) trajectory are given in subfigure (b). The variation of the ellipticity ε with the dihedral torsion angle θ for the fast (F) and (S) trajectories are given in subfigures (c) and (d), respectively. See the C1−C2 BCP in Scheme 1.

NCPs, BCPs, and RCPs; 1-DQT contains NCPs and BCPs, and a 0DQT comprises only NCPs. 2.2. QTAIM Eigenvector and Stress Tensor BCP Descriptors; the Response β and 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 states.42 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

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 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 contains 4781


Article

The Journal of Physical Chemistry A range −90.0° ≤ β ≤ 90.0° using the additional provision that β → (β − 180.0°) if β > 90.0°, using the definition for the response β:

β = arccos(e2 ·y )

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 FNAIBP 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 photoisomerization 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 Figure 4 of ref 16. The two trajectories shown in Figure 4 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 publication,16 was caused by a delayed departure from the Franck−Condon (FC) point. In dihedral torsion angle θ in an earlier publication,16 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 40) 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 Figures S1−S3a−d of the Supporting Information. For the two trajectories shown in Figure 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 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 hoptimes within ca. 10 fs. Note however that, due to the use of the geometry obtained by the OM2/MRCI semiempirical method, the SSR S1−S0 energy gaps may not completely vanish at the time of hop. In the NAMD TSH formalism, the surface hop is initiated when the nonadiabatic 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 SA-REKS functional,26,27 see the Supporting Information for more detail. The occupation numbers of the strongly occupied active orbital

(2)

In eq 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 antiparallel to the Y-axis. Therefore, the QTAIM vector-based interpretation of a double bond and single bond, the response β, relates to the degree of 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 toward 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 motor,49 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. 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 an accurate account of both dynamic and nondynamic 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 multireference 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 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 and 51, where the SSR geometries matched the geometries from MRCISD within 0.06 A. In this 4782


Article


Figure 2. Variation of the difference (S1−S0) of the stress tensor stiffness Sσ with time for the fast (F) trajectory is given in subfigure (a), and the corresponding values for the slow (S) trajectory are given in subfigure (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 subfigures (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 given in subfigures (e) and (f), respectively. The stiffness S which does not vary in accordance with physical intuition is provided in Supporting Information S4.

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 Figures S1−S3e−h. The SSR-BH&HLYP/631G* charges are generally consistent with the charges obtained

of SSR are shown in Figures S1−S3e−h of the Supporting Information. According to Figures S1−S3e−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 4783


Article


Figure 3. Variation of the stress tensor response βσ with time and 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 subfigures (a) and (c), respectively, the corresponding values for the S trajectory are given in subfigures (b) and (d), respectively. See the caption of Figure 1 for 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 afterward.

by the OM2/MRCI method, see Figure 5 of the Supporting Information of ref 16. The QTAIM and stress tensor analysis was performed with the AIMAll58 suite. The SSR calculations were carried out using the COLOGNE2012 program.59

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 toward the product conformation. Note, however, that conical intersections are almost never reached in the NAMD simulations and the nonadiabatic population transfer usually occurs in their close proximity. As illustrated by Figures S1−S2 of the Supporting Information, 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 that hold down the progress of this type of trajectories toward the conical intersection region. The nature of these factors will be analyzed in the following sections 3.1 and 3.2. 3.1. 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 demonstrates differences in the fast (F) and slow (S) EP → ZM trajectories. In the next section, section 3.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 F-NAIBP motor. The definition of the dihedral

3. 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 Figures S1−S2 of the Supporting Information. The excited electronic state remains almost purely diradicaloid on the segment of the trajectory leading away from the FC point toward the S1 4784


Article


Figure 4. Molecular graphs with lowest and highest topologically complexities (Σbrc) of the fast (F) trajectory are shown in subfigure (a) at the hop-time = 352 fs and subfigure (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 subfigure (c) and the highest Σbrc value at time = 456 fs in subfigure (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 closedshell BCPs and bond-paths.

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. For the F trajectory the F-NAIBP motor proceeds directly toward the S1 minimum and the CI, see Figure S1a,c,e of the Supporting Information S1. Conversely, the S trajectory is characterized by an extended time period when the molecule remains in conformations near the FC point, see Figure S1b,d,f. Although this behavior is not apparent from the variation of the relative energy ΔE with time in Figure S1b,d, it is confirmed by inspecting the dihedral torsion angle θ in Figure 4 of an earlier publication.16

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 and 574 fs for the F and S trajectories, see Figure 1a,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. 4785


Article


Figure 5. Plots of the variation of the numbers of BCPs (b) with time for the fast (F) and slow (S) trajectories are presented in subfigures (a) and (b), respectively, and corresponding plots for the variation of BCPs (b) with the dihedral torsion angle θ are given in subfigures (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 2D quantum topology phase diagram (QTPD) of the solutions of the Poincaré−Hopf relation of the F-NAIBP molecular motor of the F and S trajectories and hop-times in fs are shown in subfigures (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 Poincaré−Hopf relation. The blue crosses represent the forbidden solutions of the Poincaré−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 maintenance of constant values of the ellipticity ε of the axial C1−C2 BCP for approximately the first 150 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 3.2 we will examine the role that

the intramolecular closed-shell 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 1c,d, particularly for the S trajectory in 4786


Article


Figure 6. Variation of the local total energy density H(rb) with time for the closed-shell BCPs with H(rb) < 0 for the fast (F) trajectory in subfigures (a), (c), and (e). H(rb) values for the slow (S) trajectory results are presented in subfigures (b), (d), and (f). All plots correspond to the S1 electronic state. The C--C BCPs and 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; see the caption of Figure 1 for further details. The plots for the S0 excited state are virtually indistinguishable and are provided in Supporting Information S5.

the S1 state that is stalled at θ = 20.0°, see Figure 1d. 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 CI.16

The variations of the dif ference S1−S0 of the stress tensor stiffness Sσ of the axial C1−C2 BCP with time for the F and S trajectory, see Figure 2a,b, is consistent with the variation of ellipticity ε with time from Figure 1a,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 and 400 fs, respectively, in agreement with the time period for which the 4787


Article

The Journal of Physical Chemistry A 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 that the 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 2c,d, respectively. 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 trajectories, see Figure 2e,f, respectively. For the duration of the F and S trajectories there is no discernible 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 3a−d. We notice that for the variation of βσ of the axial C1−C2 BCP with the dihedral torsion angle θ, for the S trajectory, the FNAIBP motor is stalled at βσ ≈ 20.0° for both the S0 and S1 states consistent with the ellipticity ε and stress tensor stiffness Sσ, see Figures 1d and 2d, 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 in section 2.2 for further discussion 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 the axial C1−C2 BCP bond-path breaks completely and therefore does not respond to the application of the dihedral torsion θ. The nonstep-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. 3.2. 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 3.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. 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 bondpaths and BCPs that separate the stator and the rotor blades of the F-NAIBP motor other than the axial C1−C2 BCP, see Figure

4a. 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 5a. The corresponding variation of the numbers of RCPs (r) also shows the minimum possible value of r for the F-NAIBP motor, see Supporting Information 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 5b. The corresponding variations of the numbers of BCPs (b) with dihedral torsion angle θ for the S and F trajectories are given in Figure 5c,d, 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 5e,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 5e. 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 5f. 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 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 Figures 1 and 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 6a,b. This behavior does not occur for the corresponding shared-shell H13−C9 BCP in the H-NAIBP motor, see Supporting Information 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 4788


Article


Figure 7. Variation of the local total energy density H(rb) with dihedral torsion angle θ for the closed-shell BCPs with H(rb) < 0 for the fast (F) and slow (S) trajectory results are presented in subfigures (a)−(c). All plots correspond to the S1 electronic state. See the caption of Figure 1 for further details. The corresponding values for the S0 state are provided in Supporting Information S5.

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 Figures 1d, 2d, and 3d. 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. 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 7a−c. The fact that the FNAIBP motor is stalled at θ ≈ 20.0° indicates the effect that the

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 Supporting Information S8 and S9. 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 Figures 1c and 2c. This can be compared with the 4789


Article


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, which 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 Supporting Information S8 and 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 center (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 anchored, may have on the detected intramolecular interactions; especially if the intermolecular interactions will strengthen or weaken the intramolecular ones.

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 Supporting Information 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.

4. CONCLUSIONS Here, we report on the first investigation that combines nonadiabatic 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 were 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 F-NAIBP 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 hop-times 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 F-NAIBP motor resembles the form of a dielectric constant for a damped oscillator passing through a resonance; 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) were minimal at the hop-time for both trajectories. Along the F trajectory, the only intramolecular interaction between the blades of the FNAIBP 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 hoptime would hinder the progress of the motor.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.7b02347. 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 4790


Article


■

(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.; Schapiro, I.; Valentini, A.; Frutos, L. M.; Pieri, E.; Olivucci, M.; Léonard, J.; Gindensperger, E. Directionality of DoubleBond 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. (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; Ferré, 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. (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.

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) (PDF)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: mike.fi[email protected]. *E-mail: [email protected]. ORCID

Samantha Jenkins: 0000-0002-4288-7653 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS 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.

■

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 LightDriven Rotary Molecular Motors. Org. Biomol. Chem. 2008, 6, 507−512. (8) Ferná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. (10) Kazaryan, A.; Kistemaker, J. C. M.; Schä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.; Schä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. 4791


Article

The Journal of Physical Chemistry A (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.; Herná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) Herná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-2-Butene and Its Isomerization Energy to Trans-2-Butene, Revisited. J. Comput. Chem. 2016, 37 (1), 143−154. (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.; Huan, G.; Kirk, S. R.; Jenkins, 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. (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.; MartIń ez, 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; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Springer: Heidelberg, 2015; Vol. 368, pp 445− 476. (54) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 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: Condens. Matter Mater. Phys. 1988, 37, 785−789. (56) Becke, A. D. A New Mixing of Hartree−Fock and Local Densityfunctional 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.; Gräfenstein, J.; Izotov, D.; Gauss, J.; He, Y.; Wu, A.; Polo, V.; Olsson, L.; et al. COLOGNE2012; Southern Methodist University: Dallas, TX, 2012.

4792


QTAIM and Stress Tensor Characterization of ... - ACS Publications

Recommend Documents