Field-Induced Conformational Change in a Single ... - ACS Publications

Nov 28, 2016 - Vincent Pohl and Jean Christophe Tremblay. Institute for Chemistry and ... *E-mail: [email protected]. Phone: +49 30 838 58150...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Field-Induced Conformational Change in a Single-MoleculeGraphene−Nanoribbon Junction: Effect of Vibrational Energy Redistribution Vincent Pohl and Jean Christophe Tremblay* Institute for Chemistry and Biochemistry, Freie Universität Berlin, Takustrasse 3, 14195 Berlin, Germany ABSTRACT: First-principles mechanistic investigations of fieldinduced switching in an oligo(phenylene ethynylene) derivative attached to graphene nanoribbon leads are presented. It was shown that torsion of the oligomer unit causes an interruption of the conjugated π-system along the nanoribbon direction, thereby drastically reducing the conductivity of the graphene wire. In this article, we investigate the dynamical aspects associated with the switching process, with particular emphasis on vibrational energy redistribution. First, a microscopic model Hamiltonian for the reaction path coupled to the phonons of the nanoribbon leads is parametrized using density functional theory calculations. The conformational change to access the energetically unfavored structure is induced by applying an external static electric field in the spirit of a traditional field effect transistor. Using the reduced density matrix formalism, we perform ground state quantum dynamics to simulate the complete switching cycle. The switching process is characterized by three distinct time scales originating from different physical phenomena and is found to be quantitative and reversible for experimentally accessible gate voltages. Analysis of the energy flow during the dynamics shows that energy is mainly dissipated to only a few transversal acoustic phonons of the graphene nanoribbon frame.



INTRODUCTION Transistors are elementary switching units forming the basis of all logical operations in modern electronic applications. With semiconductor-based electronics entering the nanoscale, devices using molecules as active components increasingly appear as promising alternatives to conventional silicon technology. A variety of switching mechanisms have been investigated over the past decade (see, for example, refs 1−3 and references therein), ranging from isomerization4−7 to ringopening8−12 and hydrogen transfer reactions13,14 via STMinduced tautomerization15−17 and redox reactions at finite bias voltage.18,19 The conductance properties of these molecules bound to metal electrodes have been thoroughly studied,20−23 while embedding between graphene nanoribbons (GNR) leads have also received a fare share of attention in recent years.24−29 The latter are very promising materials as molecular wires30 and field-effect transistors,31 since they are virtually one-dimensional building blocks and have advantageous electronic transport properties even at room temperature.32 Further, they bind covalently to organic molecules, which may affect their conducting properties. As a promising candidate for this type of molecule−GNR junction, the oligo(phenylene ethynylene)s (OPEs) can be highlighted. OPEs have received great attention over the past decade, in particular, for applications as molecular wires (see, e.g., ref 33 and references therein). Besides, STM experiments for OPE derivatives attached to gold surfaces4,34,35 have shown that stochastic conformational changes can be triggered in these molecules. In order to explain this phenomenon, several © XXXX American Chemical Society

mechanisms have been investigated (for an overview, see ref 3). The current state of knowledge indicates that the molecule switches from an upright to a tilted orientation leading to a change in the distance to the STM tip, thereby affecting the conductance.35 Alternatively, it was identified theoretically that manipulating the torsion angle between the aromatic units in OPE wires can be used to break the conjugation and reversibly control the current flowing through the device (see, e.g., ref 1). In this spirit, the conduction properties of OPE as a model molecular junction were also investigated theoretically by Agapito and Cheng24 using a combination of density functional theory, atom-centered basis functions, and nonequilibrium Green’s functions. In their work, nitro-substituted OPE (OPENO2) was bound to two 3-zigzag graphene nanoribbons (ZGNR) serving as electrodes, and the length of the central scattering region was chosen such that the ZGNR−contact length was converged for transport calculations. From timeindependent quantum mechanical calculations, two distinct conformers (logical states) differing in the angle φ between the central nitrophenyl group and the ZGNR contacts were found: a conformer with high conductance (the ON state) and a conformer with low conductance (the OFF state), see Figure 1. While in the ON conformation, the molecule lies in the plane of the ZGNR (φ ≈ 0°) and the central nitrophenyl group stands perpendicular to the ZGNR plane in the OFF Received: September 24, 2016 Revised: November 21, 2016 Published: November 28, 2016 A

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

presented and analyzed in the Results and Discussion before concluding remarks summarizing our most important findings.



METHODOLOGY System−Bath Decomposition. The model system is composed of a central nitro-substituted oligo(phenylene ethynylene) unit connected covalently to two 3-zigzag graphene nanoribbons (ZGNR), see Figure 1. The quantum chemical characterization of the system is performed using a finite section of a free-standing ZGNR and without any constraints upon geometry optimization along the reaction path. The end points of the leads are saturated using hydrogen atoms to avoid dangling bonds. For subsequent normal-mode analysis, required to obtain all phonon frequencies of the finite model, these linking hydrogen atoms have been kept frozen. This simulates that the extension of the free-standing ZGNR contact is an infinite lead deposited on a metallic substrate. More technical details are given in the Computational Details. The switching process can be rationalized by following the dynamics along the torsional coordinate of the central unit, φ. In previous theoretical investigations, the planar structure was found to have a significantly higher conductivity than the conformation in which the central unit is perpendicular to the GNR plane. Despite its apparent simplicity, a realistic treatment of the one-dimensional problem at hand requires a dissipative dynamical model in which energy can be transferred from the system mode, φ, to the remaining vibrational normal modes of the molecule {Qi} in the form of heat. We treat this problem using a system−bath decomposition, where the effective system coordinate is weakly coupled to a bath consisting of harmonic oscillators and driven by an external electric field E⃗ . The full Hamiltonian of the system and the bath reads

Figure 1. (Top) Cartoon shows a side view of the proposed device. Extended ZGNR wires lie on a substrate, represented by the shaded areas. Field-induced dynamics takes place between finite sections of free-standing ZGNR leads. Dashed box defines the finite cutout of the extended system used in the present dynamical investigations, as shown in the bottom panel. (Bottom) Conceptual sketch of the proposed switching procedure of OPE-NO2 bound to two ZGNR leads, adapted from Agapito and Cheng (cf. ref 24). Upon application of a gate bias, the nitrophenyl group aligns with the external electric field E⃗ , inducing a rotation of the central unit about the torsion angle φ. This leads to a significant decrease in conductivity.

Ĥ (φ , {Q i}, E ⃗) = Ĥ S(φ , E ⃗) + ĤB({Q i}) + Ĥ SB(φ , {Q i})

conformation (φ ≈ 90°). The rotation about the triple bonds connecting OPE-NO2 to the ZGNR leads causes the collapse of the delocalized π-system. This is reflected by a significant decrease in the system conductivity, with great potential as a molecular switching device. In the present article, we model the dynamical switching cycle between the two conformers induced by static electric fields using dissipative quantum dynamics. A conceptual sketch of the mechanism is depicted in Figure 1. The fundamental questions we aim to address are (1) can the conductivity be controlled reliably and reversibly by applying static electric fields, (2) how does vibrational energy redistribution affect the switching dynamics, and (3) what is the time scale associated with the switching process. As can be seen from Figure 1, a onedimensional reaction path can describe the conformation change adequately, while all remaining vibrations of the molecular bridge and of the nanoribbon will lead to energy relaxation. We thus propose a system−bath decomposition of the single-molecule-graphene−nanoribbon junction to enable reduced-dimensional dissipative quantum mechanical treatment of the system dynamics. For this purpose, the reaction path Hamiltonian as well as the vibrational coupling between the active mode and the bath are parametrized from first-principles using density functional theory, offering a fully microscopic description of the switching process. The paper is structured as follows. Methodology introduces the model Hamiltonian and all theoretical tools used in this work. In the Computational Details, the numerical methods and the parametrization are described. The results are

(1)

with the bath Hamiltonian Ĥ B({Qi}) and the system−bath couplings Ĥ SB(φ,{Qi}). Note that the field−bath interaction is neglected throughout. The system Hamiltonian Ĥ S(φ,E⃗ ) defines the rotational nuclear eigenstates and energies of the central unit, i.e., Ĥ S|n ⟩ = ϵn|n ⟩. It is defined as ℏ2 ∂ 2 Ĥ S(φ , E ⃗) = − + VS(φ , E ⃗) 2I ∂φ 2

(2)

The moment of inertia of the nitro-substituted oligo(phenylene ethynylene) group is defined as

I=

ℏ 4πcB

(3)

where B is the associated rotational constant and c is the speed of light. The potential VS(φ,E⃗ ) is the effective one-dimensional adiabatic reaction path along the angle φ. To treat the potential distortion due to the presence of a static electric field, a semiclassical interaction term describing the linear coupling between the field and the dipole moment of the OPE unit is added to the potential energy curve as VS(φ , E ⃗) = Vel(φ) − μ ⃗ (φ) ·E ⃗

(4)

The ground state potential energy curve Vel(φ) and the associated dipole moment μ⃗(φ) are obtained by solving the time-independent electronic Schrödinger equation parametrically along the reaction path, Ĥ elΦ(r;φ) = Vel(φ)Φ(r;φ). In the following, polarization effects are neglected since they were B

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C found to have a rather marginal effect on the potential energy curve when solving the time-independent electronic Schrödinger equation with explicit inclusion of the strongest field used in the present investigation. To derive an expression for the bath degrees of freedom and their couplings to the system, we expand the total system−bath potential V(φ,{Qi},E⃗ ) as a Taylor series expansion around each configuration along the adiabatic reaction path but close to the minimum in the bath modes {Qi}

λ i (φ , E ⃗ ) =

=

⎛ ∂V ⎞ 1 ⎛ ∂ 2V ⎞ V = V0 + ⎜ ⎟φ + ⎜ 2 ⎟φ 2 + ⋯ 2! ⎝ ∂φ ⎠ ⎝ ∂φ ⎠ +

i

+

1 2!

⎛ ∂V ⎞ 1 ⎟⎟Q i + 2! ⎝ ∂Q i ⎠

∑ ⎜⎜

⎛ ∂ 2V ⎞ 1 ⎟⎟φQ i + ∂ φ ∂ ! Q 3 ⎝ i⎠

∑ ⎜⎜ i

⎛ ∂ 2V ⎞ ⎟Q 2 + ⋯ 2⎟ i ⎝ ∂Q i ⎠ ⎛ ∂ 3V ⎞ ⎟⎟φ 2Q i + ⋯ 2 ∂ φ ∂ Q ⎝ i⎠

∑ ⎜⎜ i

⎞ ⎛ 3 1 + ∑ ⎜⎜ ∂ V ⎟⎟φQ iQ j + 1 3! i , j ⎝ ∂φ∂Q i ∂Q j ⎠ 4! φ 2Q iQ j + ⋯

⎛ ∂V (φ , {Q }, E ⃗) ⎞⎛ ∂x ⎞ i ⎟⎟⎜⎜ A ⎟⎟ x ∂ ⎝ ⎠⎝ ∂Q i ⎠ A

∂VS(φ , E ⃗) ∂Q φ

Qi =

The first line defines the effective potential of the system VS(φ) and the second line is the bath potential VB({Qi}), which is hereafter treated in the harmonic approximation. The third and fourth lines describe the system−bath coupling corresponding to the excitation of one and two bath phonons, respectively. The normal modes of the bath are found to vary only marginally as a function of φ. Consequently, they are computed only once, at the nuclear configuration corresponding to the transition state. This defines uniquely the spectral density of the bath, together with a linear coordinate transformation between the Cartesian coordinates and the normal modes. This linear transformation greatly simplifies the structure of the reaction path Hamiltonian.36 Since test calculations have shown that the two-(and more)phonon contributions are small for the system at hand, neglecting all higher order system−bath couplings yields the following effective Hamiltonian

with VBharm.({Q i}) =

1 2

xA =

∑ Oi AmA−1/2Q i

(10)

with mA being the mass of atom A. Consequently, the derivatives appearing in eq 8 can be simply evaluated and the interaction potential takes the form λ i (φ , E ⃗ ) =

∂VS(φ , E ⃗) ∂Q φ

= ci

where ci =

(

∑ OAφOi A A

∂VS(φ , E ⃗) ∂φ

∑A OAφOi A I

(11)

).

Dissipative Quantum Dynamics. The reduced density matrix formalism represents an efficient tool to investigate the dissipative dynamics of the system along the one-dimensional rotational reaction path. Tracing out the bath degrees of freedom from the dynamical equations for the total density matrix operator, the time evolution of the system can be described by the dissipative Liouville von Neumann equation of the form

(6)

⎛ ∂ 2V ⎞ ∑i ⎜ 2 ⎟Q i2 and the anharmonic ⎝ ∂Q i ⎠

∂ρ ̂(t ) l = − [Ĥ S , ρ ̂(t )] + 3 D ρ ̂(t ) ∂t ℏ

∂V (φ , {Q i}, E ⃗) ∂Q i

∑ OAimA1/2xA i

coupling functions λ i (φ , E ⃗ ) =

(9)

A

V (φ , {Q i}, E ⃗) ≃ VS(φ , E ⃗) + VBharm.({Q i}) i

(8)

where the columns of the unitary matrix define the mode vectors, {Qi}, associated with frequencies {Ωi} as diagonal entries of the matrix Ω. The unitary matrix O can be straightforwardly inverted, which yields the following relations to the Cartesian coordinates

(5)

∑ λ i (φ , E ⃗ )Q i

A

FM O = ΩO

⎛ ⎞ 4 ∑ ⎜⎜ 2∂ V ⎟⎟× i , j ⎝ ∂φ ∂Q i ∂Q j ⎠

+ ⋯three‐(and more)‐phonon processes

+

⎛ ∂Q φ ⎞⎛ ∂x ⎞ ⎟⎜⎜ A ⎟⎟ ⎝ ∂xA ⎠⎝ ∂Q i ⎠

∑⎜

where xA are the Cartesian coordinates of the atoms. The last line was obtained under the assumption that the normal mode associated with the rotation around the central unit is linearly related with the torsion angle, Qφ = √I φ, consistent with our choice for the reaction path. Here, the normal modes including the system mode are defined by diagonalization of the field-free mass-weighted Hessian at the geometry of the OFF conformer, FM, defining the unitary matrix O

∑ ⎜⎜ i

∂Q i

∑ ⎜⎜ A

=

∂V (φ , {Q i}, E ⃗)

{Q i}ref

(12)

The first term represents the coherent evolution of the system along the reaction path under the influence of the system Hamiltonian, eq 2. The last term is a dissipative Liouvillian describing the coupling of the system with the bath modes for which we choose the Lindblad form

(7)

Following Andrianov and Saalfrank37 and Sakong et al.,38 the coupling coefficient λi(φ,E⃗ ) for the ith bath mode can be written using the chain rule as C

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C 3 D ρ (̂ t ) = −

1 2



normal mode of the molecule. Consequently, this mode is excluded when summing over the bath modes. Further, in order to account for discretization errors, the delta function in eq 20 is approximated by a Lorentzian, with the full width at halfmaximum γ set to the energy difference between the vibrational ground state and the first excited state of the system mode. Other physically sound choices for this energy smearing parameter will yield different values for the energy relaxation rates. Nonetheless, the time scales and the relation between the different relaxation channels were seen to remain unaffected. In the present treatment, the system and bath are considered separable and the bath remains in thermal equilibrium at all times. This picture could be altered by including electron− phonon coupling. As will be shown in the discussion, thermalization of the phonons due to electron−phonon coupling within the leads is orders of magnitude faster than the switching dynamics, justifying the use of the Markovian approximation for the bath. To quantify the energy flow between the system and the bath, the temporal change of the system energy can be simply defined in the system eigenstate basis representation as



∑ {[Ĉk , ρ (̂ t )Ĉk ] + [Ĉk ρ (̂ t ), Ĉk ]} k

(13)

This ensures that the reduced density matrix is semipositive definite and, thereby, retains it probabilistic interpretation. The Lindblad operators Ĉ k are chosen here to mediate energy relaxation. The different channels are described by a composite index k = {m → n}, i.e.

Ĉ k →

ℏΓm → n |n⟩⟨m|

(14)

where Γm→n is the transition rate from the rotational nuclear eigenstate |m⟩ to eigenstate |n⟩, mediated by vibration−phonon coupling. In this work, the reduced density matrix operator is represented in the basis of nuclear eigenstates sustained by the system Hamiltonian, eq 2, i.e.

∑ ρmn (t )|m⟩⟨n|

ρ ̂( t ) =

(15)

mn

This choice simplifies greatly the equations of motion and the definition of the relaxation process using eq 14. The matrix elements of eq 12 in a basis of the N nuclear eigenstates take the form

∂⟨Ĥ S⟩(t ) ∂t

N

∂ρnn (t )

and for m ≠ n

=

⎛l ℏγ ⎞ = −⎜ (ϵm − ϵn) + mn ⎟ ρmn 2 ⎠ ⎝ℏ

∂ρmn (t ) ∂t

∂⟨Ĥ S⟩(t ) ∂t

∑ (Γm→ j + Γn→ j)

∑ Γ(mi)→ n



Each one-phonon contribution reads Γ(mi)→ n = π |⟨n|λi(φ , E ⃗)|m⟩|2

† |Â i + Â i |2 δ(|ϵm − ϵn| − ℏΩi) Mi Ωi

(20)

with † Â i + Â i

2

⎧ downward transition ⎪⟨ni ⟩ + 1 =⎨ ⎪ upward transition ⎩ ⟨ni⟩

(22)

∂⟨Ĥ S⟩(i)(t ) , ∂t

with individual components

∑ (ϵn − ϵj){Γ(ji→) n ρjj − Γ(ni→) j ρnn )} n 60 MV/cm), the population in the target state rapidly increases within the first 10 μs, until it reaches a plateau where more than 50% of the population can be found in the sector of φ ∈ [80°, 100°]. In the best scenario Ez = −103 MV/cm, about 99% of the wave packet can be found within the target state region, i.e., the molecule is completely switched to the nonconducting conformer. Interestingly, for the reverse step where the field has been abruptly switched off (1, blue shaded area), the delocalization dynamics proceeds almost instantaneously in all cases. The rapid delocalization of the wave packet will instantaneously lead to a measurable current in the device. This asymmetry of the rising time for switching on and off the device is due to the fact that most conformations along the reaction path belong to the ON state of the system, while the OFF state where conjugation of the orbitals is completely broken can be only found in a small domain in configuration space. In order to gain deeper mechanistic insight into the dynamics, two exemplary scenarios will be discussed in more detail: Ez = −50 MV/cm (cf. dashed red curve in Figure 3) and Ez = −103 MV/cm (cf. solid red curve in Figure 3). Contour plots of the time-dependent nuclear density for the three steps of the switching cycle are depicted in Figure 5a and 5b. Red

is time independent and can be evaluated before the propagation for computational efficiency. Figure 4 shows the time evolution of ⟨P̂ 0⟩(t) for the selected electric field strengths at T = 0 K. The red shaded area

Figure 4. Target state population ⟨P̂ 0⟩(t) as a function of time t and electric field strengths Ez for steps 0 (red shaded area) and 1 (blue shaded area) of the switching cycle at T = 0 K.

corresponds to step 0 of the switching cycle, when the electric field is first turned on. The blue shaded area corresponds to step 1, when the electric field is abruptly turned off. Let us first consider step 1. For electric field strengths below |Ez| = 30 MV/cm, no population can be found in the target state, and hence, no significant switching in the conductance of the device is expected. This can be understood by looking at Figure 3,

Figure 5. (a−d) Contour plots of the nuclear density as a function of time t and torsion angle φ (left panels) and the current at a bias voltage of U = 1 V as a function of time t (right panels) for four different setups at the three steps of the switching cycle: 0 (red shaded area), step 1 (blue shaded area), and step 2 (red shaded area). Four different setups are considered: (a) Ez = −50 MV/cm at T = 0 K (green solid curve), (b) Ez = −103 MV/cm at T = 0 K (orange solid curve), (c) Ez = −50 MV/cm at T = 300 K (green dashed curve), and (d) Ez = −103 MV/cm at T = 300 K (orange dashed curve). (e) On/off ratio (at the equilibrium distribution of steps 1 and 2) as a function of the bias voltage U. Bias voltage applied for a−d is marked by a black solid line. G

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 6. Electron density ρel(t) (gray isosurface) and difference density Y(t) (blue and orange isosurface) for the initialization step 0 with Ez = −103 MV/cm and T = 0 K at eight different snapshots representing the three distinct time scales of the dynamics: coherent regime (upper panel: t = 0, 0.4, 0.7, and 1.0 ps), delocalization (lower panel: t = 0.2 and 10 ns), and vibrational relaxation (lower panel: t = 1.0 and 3.0 μs). Blue isosurface corresponds to an increase and the orange isosurface to a decrease in probability density. Isosurface values are ±0.01 a−3 0 . For the sake of orientation, the molecular structure of the initial geometry (φ = 0°) is drawn in ball-and-stick representation.

shaded areas correspond to electric field turned on, and blue shaded areas correspond to field turned off. They are separated by black solid lines. The angular range of the target state, the OFF conformation, is denoted by black dotted lines. The time evolution of the wave packet for both electric field strengths exhibits some common features: At the beginning of each step, the localized nuclear density instantaneously spreads over a wide angular range having a maximum probability distribution at the classical turning points of the PEC. This roughly corresponds to the classical pendulum behavior of a highly excited harmonic oscillator. Due to coupling with the phonon bath, the nuclear density relaxes sequentially toward the ground vibrational state without a major change in the shape of the distribution of the wave packet. After switching off the electric field in step 1, the major fraction of the density relaxes to the thermodynamically less favored ON′ conformer due to a slight asymmetry of the PEC at the transition state. It is worth noting that for both scenarios (Ez = −50 and −103 MV/cm) the dynamics of step 1 is virtually identical. Interestingly, even if the starting point of the dynamics of steps 0 and 2 are different, the time evolution of the density is identical (see red shaded areas). Further, no memory effect from the previous ON−OFF cycles can be observed, provided the system is left to relax for a long enough period. When comparing steps 0 and 2 (cf. red shaded areas in Figure 5a and 5b) for the two field strengths, it can be recognized that for the weaker field (Ez = −50 MV/cm), cf. Figure 5a, the time evolution of the density is asymmetric, with a shoulder extending toward the ON structure. This is a consequence of the asymmetric shape of the PEC, as can be seen from Figure 3 (red dashed curve). At this lower field strength the torsional ground state that is reached at the end of the relaxation dynamics is not exactly located in the region of the OFF conformation (φ ≈ 90°) but slightly bent toward the ON′ configuration. Consequently, although the final nuclear

density is well localized at around 120°, the conjugation of the π-systemand with it the conductanceis not interrupted in this case. This can be confirmed by inspection of the frontier orbitals at φ ≈ 120° (see Figure 2, third column), which remain completely delocalized on the whole device. For Ez = −103 MV/cm (cf. Figure 5b), a different picture emerges. Here, not only is the dynamics twice as fast as in the first case but also the final nuclear density is nearly perfectly located in the target state OFF, since the potential energy curve is more symmetric between ON and ON′. When comparing the nuclear dynamics with the time evolution of ⟨P̂ 0⟩(t) at step 1, the rapid drop of ⟨P̂0⟩(t) (cf. Figure 4, blue shaded area) is seen to first originate from the nearly instantaneous spread of the wave packet and the ensuing localization of the nuclear wave packet outside of the target state region. The time-dependent current can be reconstructed from the evolution of the nuclear density, yielding a statistical weight for the current through the device in a given conformation. The current dynamics at U = 1 V can be seen in the right panels of Figure 5a and 5b. In steps 0 and 2, the conductivity drops slowly by several orders of magnitude, until it reaches a plateau, when the nuclear wave packet is equilibrated. As soon as the electric field is turned off, i.e., in step 1, the conductivity increases almost instantaneously to the original level. Hence, it is the switching off phase, steps 0 and 2, that will ultimately limit the efficiency of the device. It can be stressed that for Ez = −103 MV/cm (cf. Figure 5b) the current not only drops as twice as fast but also drops by one order of magnitude more than for Ez = −50 MV/cm (cf. Figure 5a). This is reflected in the on/off current ratio, see Figure 5e. While it is 37 for Ez = −50 MV/cm, it reaches 835 for Ez = −103 MV/cm. Increasing the temperature has a qualitatively similar effect on the conductivity as reducing the field strength. This can be seen from Figure 5c and 5d, where the dynamics of the same switching cycle as in Figure 5a and 5b are presented at room H

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. Time-integrated energy flow ⟨E ⟩i for steps 0 (red curve) and 1 (blue curve) as a function of the energy of the respective bath mode νi. Peaks have been broadened by a Lorentz distribution of width σ = 0.1 meV: (left) Ez = −50 MV/cm at T = 0 K and (right) Ez = −103 MV/cm at T = 0 K.

bridge and in the ZGNR leads. The system, which is now highly rotationally excited, relaxes during the next microsecond to its equilibrium geometry, the OFF conformer (third and fourth lower panels of Figure 6). In a classical picture, the dynamics on the latter two time scales may be understood as that of a damped harmonic oscillator, that is, the central nitrophenyl group oscillates from one classical turning point to the other, and coupling to the phonon bath causes the torsional amplitude to gradually decrease. An attractive feature of the ZGNR−OPE device is that the dynamics takes place mostly in the bridge molecule, with a slight upward shift of small parts of the ZGNR contacts. This is a very promising property for a possible application as a graphene nanoribbon-based switching device. Note that the time scales discussed above shall only be seen as an order of magnitude. By varying the length of the leads, the associated phonon spectrum becomes denser, which will affect the coupling to the system mode. We found that the key features of the dynamics and properties discussed in this work are not significantly affected by size effectsat most a factor of 2 for much smaller GNR leads. This can be rationalized on the basis that the system mode is very localized on the bridge (cf. Figure 6) and mostly couples to the local phonons of the contacts. We thus consider the lead length to be well converged for the aims of the present study. An important assumption of the present dynamical investigation is that the phonon bath, i.e., the phonons in the ZGNR leads, remains in thermal equilibrium at all times. Thermalization within the leads is due to electron−phonon coupling. Let us first divide conceivable electronic excitations into (I) local excitations of the system (the bridge) and (II) excitations of the leads. (I) In order to get a rough estimate of the transition energy which is required for local electronic excitations, let us look for the first change in the local nodal structure on the bridge within the space of unoccupied molecular orbitals. The energy difference between this particular molecular orbital and the highest occupied molecular orbital corresponds loosely to the first local excitation, associated with a local electronic transition energy Δϵ. Upon inspection of the orbitals, all unoccupied MOs of the finite system−bath model used in this study are seen to preserve the same nodal structure locally on the bridge at energies below Δϵ = 2.6 eV for ON and Δϵ = 2.8 eV for OFF (see Computational Details). These are three orders of magnitude larger than the nuclear transition energies. Consequently, local electronic excitations induced by electron−phonon coupling can be safely neglected in this device. (II) For excitations of the leads a

temperature. Since the wave packets are much more spread at all times, the projector in the target region (see eq 25) and thus the on/off ratio is decreased. Interestingly, the on/off ratio for Ez = −103 MV/cm at T = 300 K is still higher than that at Ez = −50 MV/cm and T = 0 K. This implies that a higher gate bias not only yields a better signal-to-noise ratio but also improves the thermal stability of the device. The conductive properties discussed above can be further illuminated by investigating the time dependence of the electron density ρel(t) during the switching dynamics. To this end, we also make use of the difference density, which is connected to the electron yield Y(t) by Y (t ) =

∫0

t

∂ρel (t ′) ∂t ′

dt ′

= ρel (t ) − ρel (0)

(26)

The latter is a scalar field which shows as a function of time, where electrons are lost and gained in space and hence how the electron distribution evolves in time. Figure 6 shows the electron density (gray isocontours) and the electronic yield (blue and orange isocontours for positive and negative yields, respectively) for selected snapshots of step 0 with Ez = −103 MV/cm. The corresponding nuclear dynamics is depicted in the right panel of Figure 5. Note that for the other steps of the switching cycle the picture is qualitatively the same. The electron dynamics reveals three distinct time scales associated with three different physical origins: coherent regime (t < 0.1 ns), delocalization (t < 10 ns), and vibrational relaxation (t < 10 μs). In the early stages (t ≈ 1 ps, see top four panels of Figure 6), the nuclear and electronic wave packets remain coherent and the electronic density follows adiabatically the nuclear motion. The first half vibrational period also reveals that the ZGNR leads remain unaffected by the rotational motion as long as the OPE bridge has not crossed the OFF conformation (t < 0.4 ps). As soon as the OPE group points toward the ON′ conformation (third and fourth upper panels), the ZGNR leads rise above the device plane. As will be demonstrated below, this induces stronger coupling to the transversal acoustic phonon modes. On intermediate time scales (t = 0.2 ns), it can be observed that the electron density of the OPE bridge spreads over a wide angular range (bottom left snapshot). This correlates with the spreading of the nuclear wave packet and the onset of energy relaxation, see Figure 5. Up to about 10 ns, the degree of electron delocalization increases, both in the I

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. Ball-and-stick representation of the most relevant relaxation channels (bath modes νi; cf. Figure 7). Cartesian displacements have been obtained from a frequency analysis and have been scaled by a factor of 10 for visualization. Blue corresponds to a positive and orange to a negative displacement. Energies of the bath modes are given next to each structure.

red shaded area in the right panel of Figure 7), where the main fraction of the energy is dissipated via the 3rd, 8th, and 14th modes. The first two thus couple mostly through the ZGNR leads and the last one via the upward shift of the bridge. Note that although the same modes are involved in energy transfer from the system to the contacts at both field strengths, the higher energy contribution, i.e., mode ν14, plays a more prominent role at higher fields. This is a consequence of the stronger deformation of the device due to the electric field, as also observed in the potential energy landscape.

different picture emerges. The electron−phonon coupling constant, as calculated by Zhang et al.,61 was found to be strong (λ = 0.14) in periodic 2-ZGNR and 4-ZGNR. The 1 phonon lifetimes of the associated bath modes, τi ≈ λ·ω , are i

thus in the picosecond regime, much faster than the dynamics of the switching process reported here. Consequently, the Markov approximation for the phonon bath is well justified, and the bath always remains in thermal equilibrium from the system mode’s perspective. Following eq 23, the time-integrated energy flow ⟨E ⟩(i) can be computed, yielding information about the fraction of the total energy which is dissipated to each bath mode. Figure 7 shows ⟨E ⟩(i) at the end of steps 0 (red) and 1 (blue) for the intermediate (Ez = −50 MV/cm, left) and strong field (Ez = −103 MV/cm, right) scenarios. The energy flow spectra exhibit strong similarities in both cases and for both field strengths. Six vibrational modes dominate the energy flow spectrum, and their relative contributions are similar in both steps of the switching cycle. The peak size of ⟨E ⟩(i), i.e., the energy transfer to the bath, correlates perfectly with the depth of the potential well (cf. Figure 3). The integral under the peaks of the timeintegrated energy flow gives a quantitative measure of the energy loss of the nuclear wave packet, e.g., approximately 0.5 eV is transferred from the system to the bath within the step 0 for Ez = −103 MV/cm. For stronger electric fields, the potential is subject to a larger deformation, and it becomes steeper in the region, where the dynamics is initiated. On one hand, this increases the transition energies, and on the other hand, the potential becomes more harmonic. As a consequence, the transition energies are more homogeneous, which modifies the importance of particular relaxation channels. The most important skeletal vibrations are sketched in Figure 8 and can be identified as transversal acoustic phonons along the ZGNR contacts. The associated vibrational frequencies all lie in the millielectronvolt regime, in the same energy range as the anharmonic transition frequencies between nuclear rotational eigenstates populated during the dynamics. These are found in the range of ΔE ∈ [1.6 meV, 2.1 meV] for Ez = −103 MV/cm and ΔE ∈ [0.7 meV, 1.4 meV] for Ez = −50 MV/cm. The four lowest modes depicted (ν3, ν4, ν6, ν8) are mostly localized on the leads and will allow coupling at the contact positions to the bridge, which only undergoes small amplitude motion. The number of nodes along the ZGNR leads is higher for modes ν14 and ν18, with also a larger amplitude motion on the bridge. The onset of a torsional motion in modes ν6 and ν8 can be also recognized. The largest energy loss is observed at step 0 with Ez = −103 MV/cm (cf.



CONCLUSIONS In the present work, we proposed a procedure for the fieldinduced switching of a nitro-substituted oligo(phenylene ethynylene) bound between two ZGNR electrodes. In this model system, a rotation of the central nitrophenyl group about two triple bonds leads to two conformers: one with high and one with low conductance. While the planar conformer possesses a completely delocalized π-system, the molecular orbitals of the perpendicular conformer are localized on either of the two sides of the molecule, which drastically reduces the conductivity of the graphene wire. In order to drive the switching device from the planar to the energetically less favored perpendicular conformer, we proposed using an external static electric field perpendicular to the molecular plane. The focus of this contribution is the investigation of the dynamical aspects of this switching process. First, we constructed a model Hamiltonian for the quasi-one-dimensional rotational motion of the central nitrophenyl group coupled to the phonons of the nanoribbon leads. Its microscopic parametrization was performed from first principles using density functional theory. Applying the reduced density formalism, we simulated the quantum dynamics for the full switching cycle. The dynamics induced by the static electric field revealed three distinct time scales associated with three different physical origins: (i) a coherent regime, where the bridge oscillates harmonically between the two planar conformations, (ii) a delocalization regime, where the associated wave packet spreads over a wide angular range, and (iii) a vibrational relaxation regime, where the wave packet incoherently relaxes completely to the ground state. This is achieved by transferring energy from the localized motion of the system mode to only a handful of transversal acoustic phonons of the ZGNR contacts. Since the device is reinitialized as soon as the field is turned off, all steps of the switching cycle share the same characteristic time evolution without showing any memory effect. Hence, by turning the external static electric J

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C field on and off, the conductivity of the system can be controlled at will in the spirit of a traditional field-effect transistor. To reduce the required time for the relaxation from a few microseconds to shorter time scales, one could either increase or decrease the size of the central group. The former would increase the steric effects and reduce the mobility of the central group with respect to the graphene contacts, while the latter would reduce its moment of inertia and allow for faster energy relaxation.



(12) Pathem, B. K.; Zheng, Y. B.; Morton, S.; Petersen, M. Å.; Zhao, Y.; Chung, C.-H.; Yang, Y.; Jensen, L.; Nielsen, M. B.; Weiss, P. S. Photoreaction of Matrix-Isolated Dihydroazulene-Functionalized Molecules on Au{111}. Nano Lett. 2013, 13, 337−343. (13) Benesch, C.; Rode, M. F.; Č ìžek, M.; Rubio-Pons, O.; Thoss, M.; Sobolewski, A. L. Switching the Conductance of a Single Molecule by Photoinduced Hydrogen Transfer. J. Phys. Chem. C 2009, 113, 10315−10318. (14) Pan, S.; Fu, Q.; Huang, T.; Zhao, A.; Wang, B.; Luo, Y.; Yang, J.; Hou, J. Design and Control of Electron Transport Properties of Single Molecules. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 15259−15263. (15) Liljeroth, P.; Repp, J.; Meyer, G. Current-Induced Hydrogen Tautomerization and Conductance Switching of Naphthalocyanine Molecules. Science 2007, 317, 1203−1206. (16) Kumagai, T. Direct Observation and Control of Hydrogen-Bond Dynamics Using Low-Temperature Scanning Tunneling Microscopy. Prog. Surf. Sci. 2015, 90, 239−291. (17) Kumagai, T.; Grill, L. Tautomerism; Wiley-Blackwell, 2016; pp 147−174. (18) Mendes, P.; Flood, A.; Stoddart, J. Nanoelectronic Devices from Self-Organized Molecular Switches. Appl. Phys. A: Mater. Sci. Process. 2005, 80, 1197−1209. (19) Seminario, J. M.; Zacarias, A. G.; Derosa, P. A. Theoretical Analysis of Complementary Molecular Memory Devices. J. Phys. Chem. A 2001, 105, 791−795. (20) Seminario, J. M.; Zacarias, A. G.; Tour, J. M. Theoretical Study of a Molecular Resonant Tunneling Diode. J. Am. Chem. Soc. 2000, 122, 3015−3020. (21) Huang, J.; Li, Q.; Ren, H.; Su, H.; Shi, Q. W.; Yang, J. Switching Mechanism of Photochromic Diarylethene Derivatives Molecular Junctions. J. Chem. Phys. 2007, 127, 094705−094710. (22) Zhao, P.; Fang, C.-f.; Xia, C.-j.; Liu, D.-s.; Xie, S.-j. 15, 16Dinitrile DDP/CPD as a Possible Solid-State Optical Molecular Switch. Chem. Phys. Lett. 2008, 453, 62−67. (23) Chen, W.; Chen, R.; Bian, B.; Li, X.-a.; Wang, L. First-Principles Study of the Electronic Transport Properties of a DihydroazuleneBased Molecular Optical Switch. Comput. Theor. Chem. 2015, 1067, 114−118. (24) Agapito, L. A.; Cheng, H.-P. Ab Initio Calculation of a Graphene-Ribbon-Based Molecular Switch. J. Phys. Chem. C 2007, 111, 14266−14273. (25) Wu, Q.-H.; Zhao, P.; Liu, D.-S. Electronic Transport of a Molecular Photoswitch with Graphene Nanoribbon Electrodes. Chin. Phys. Lett. 2014, 31, 057304−057307. (26) Cai, Y.; Zhang, A.; Feng, Y. P.; Zhang, C. Switching and Rectification of a Single Light-Sensitive Diarylethene Molecule Sandwiched Between Graphene Nanoribbons. J. Chem. Phys. 2011, 135, 184703−184708. (27) Xie, F.; Fan, Z.-Q.; Liu, K.; Wang, H.-Y.; Yu, J.-H.; Chen, K.-Q. Negative Differential Resistance and Stable Conductance Switching Behaviors of Salicylideneaniline Molecular Devices Sandwiched Between Armchair Graphene Nanoribbon Electrodes. Org. Electron. 2015, 27, 41−45. (28) Xia, C.-J.; Zhang, B.-Q.; Su, Y.-H.; Tu, Z.-Y.; Yan, X.-A. Electronic Transport Properties of a Single Chiroptical Molecular Switch with Graphene Nanoribbons Electrodes. Optik 2016, 127, 4774−4777. (29) Jia, C.; Migliore, A.; Xin, N.; Huang, S.; Wang, J.; Yang, Q.; Wang, S.; Chen, H.; Wang, D.; Feng, B.; et al. Covalently Bonded Single-Molecule Junctions with Stable and Reversible Photoswitched Conductivity. Science 2016, 352, 1443−1445. (30) Fischetti, M. V.; Kim, J.; Narayanan, S.; Ong, Z.-Y.; Sachs, C.; Ferry, D. K.; Aboud, S. J. Pseudopotential-Based Studies of Electron Transport in Graphene and Graphene Nanoribbons. J. Phys.: Condens. Matter 2013, 25, 473202−4732038. (31) Rodríguez-Manzo, J. A.; Qi, Z. J.; Crook, A.; Ahn, J.-H.; Johnson, A. T. C.; Drndić, M. In Situ Transmission Electron Microscopy Modulation of Transport in Graphene Nanoribbons. ACS Nano 2016, 10, 4004−4010.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 30 838 58150. Fax: +49 30 838 452051. ORCID

Jean Christophe Tremblay: 0000-0001-8021-7063 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors gratefully acknowledge Lukas Eugen Marsoner Steinkasserer for his help with the transport calculations and the Scientific Computing Services Unit of the Zentraleinrichtung für Datenverarbeitung at Freie Universtät Berlin for allocation of computer time. The funding of the Deutsche Forschungsgemeinschaft (project TR1109/2-1) and from the Elsa-Neumann foundation of the Land Berlin is also acknowledged.



REFERENCES

(1) Fuentes, N.; Martín-Lasanta, A.; de Cienfuegos, L. Á ; Ribagorda, M.; Parra, A.; Cuerva, J. M. Organic-Based Molecular Switches for Molecular Electronics. Nanoscale 2011, 3, 4003−4014. (2) Pathem, B. K.; Claridge, S. A.; Zheng, Y. B.; Weiss, P. S. Molecular Switches and Motors on Surfaces. Annu. Rev. Phys. Chem. 2013, 64, 605−630. (3) Zhang, J. L.; Zhong, J. Q.; Lin, J. D.; Hu, W. P.; Wu, K.; Xu, G. Q.; Wee, A. T. S.; Chen, W. Towards Single Molecule Switches. Chem. Soc. Rev. 2015, 44, 2998−3022. (4) Donhauser, Z. J.; Mantooth, B. A.; Kelly, K. F.; Bumm, L. A.; Monnell, J. D.; Stapleton, J. J.; Price, D. W.; Rawlett, A. M.; Allara, D. L.; Tour, J. M.; et al. Conductance Switching in Single Molecules Through Conformational Changes. Science 2001, 292, 2303−2307. (5) Pati, R.; Karna, S. P. Current Switching by Conformational Change in a π-σ-π Molecular Wire. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 155419−155423. (6) Choi, B.-Y.; Kahng, S.-J.; Kim, S.; Kim, H.; Kim, H. W.; Song, Y. J.; Ihm, J.; Kuk, Y. Conformational Molecular Switch of the Azobenzene Molecule: A Scanning Tunneling Microscopy Study. Phys. Rev. Lett. 2006, 96, 156106−156109. (7) del Valle, M.; Gutiérrez, R.; Tejedor, C.; Cuniberti, G. Tuning the Conductance of a Molecular Switch. Nat. Nanotechnol. 2007, 2, 176− 179. (8) Mitchell, R. H. The Metacyclophanediene-Dihydropyrene Photochromic π Switch. Eur. J. Org. Chem. 1999, 1999, 2695−2703. (9) Dulić, D.; van der Molen, S. J.; Kudernac, T.; Jonkman, H. T.; de Jong, J. J. D.; Bowden, T. N.; van Esch, J.; Feringa, B. L.; van Wees, B. J. One-Way Optoelectronic Switching of Photochromic Molecules on Gold. Phys. Rev. Lett. 2003, 91, 207402−207405. (10) Li, J.; Speyer, G.; Sankey, O. F. Conduction Switching of Photochromic Molecules. Phys. Rev. Lett. 2004, 93, 248302−248305. (11) Roldan, D.; Kaliginedi, V.; Cobo, S.; Kolivoska, V.; Bucher, C.; Hong, W.; Royal, G.; Wandlowski, T. Charge Transport in Photoswitchable Dimethyldihydropyrene-Type Single-Molecule Junctions. J. Am. Chem. Soc. 2013, 135, 5974−5977. K

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (32) Schwierz, F. Nat. Nanotechnol. 2010, 5, 487−496. (33) James, D. K.; Tour, J. M. Molecular Wires and Electronics; Springer Science + Business Media, 2005; pp 33−62. (34) Ramachandran, G. K. A Bond-Fluctuation Mechanism for Stochastic Switching in Wired Molecules. Science 2003, 300, 1413− 1416. (35) Moore, A. M.; Dameron, A. A.; Mantooth, B. A.; Smith, R. K.; Fuchs, D. J.; Ciszek, J. W.; Maya, F.; Yao, Y.; Tour, J. M.; Weiss, P. S. Molecular Engineering and Measurements To Test Hypothesized Mechanisms in Single Molecule Conductance Switching. J. Am. Chem. Soc. 2006, 128, 1959−1967. (36) Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction Path Hamiltonian for Polyatomic Molecules. J. Chem. Phys. 1980, 72, 99− 112. (37) Andrianov, I.; Saalfrank, P. Theoretical Study of VibrationPhonon Coupling of H Adsorbed on a Si(100) Surface. J. Chem. Phys. 2006, 124, 034710−034719. (38) Sakong, S.; Kratzer, P.; Han, X.; Laß, K.; Weingart, O.; Hasselbrink, E. Density-Functional Theory Study of Vibrational Relaxation of CO Stretching Excitation on Si(100). J. Chem. Phys. 2008, 129, 174702−174712. (39) TURBOMOLE, V6.5 2013, A Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since2007; available from http://www. turbomole.com (accessed Nov 20, 2016). (40) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158−6170. (41) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (42) Jónsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. Classical and Quantum Dynamics in Condensed Phase Simulations 1997, 385−404. (43) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978−9985. (44) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (45) Smidstrup, S.; Pedersen, A.; Stokbro, K.; Jónsson, H. Improved Initial Guess for Minimum Energy Path Calculations. J. Chem. Phys. 2014, 140, 214106−214111. (46) Bahn, S. R.; Jacobsen, K. W. An Object-Oriented Scripting Interface to a Legacy Electronic Structure Code. Comput. Sci. Eng. 2002, 4, 56−66. (47) Schmidt, B.; Lorenz, U. WavePacket 4.9: A Matlab Program Package for Quantum-Mechanical Wavepacket Propagation and TimeDependent Spectroscopy; 2013; see http://wavepacket.sourceforge.net (accessed Nov 20, 2016). (48) Thygesen, K. S.; Bollinger, M. V.; Jacobsen, K. W. Conductance Calculations with a Wavelet Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys 2003, 67, 115404. (49) Thygesen, K.; Jacobsen, K. Molecular Transport Calculations with Wannier Functions. Chem. Phys. 2005, 319, 111−125. (50) Strange, M.; Kristensen, I. S.; Thygesen, K. S.; Jacobsen, K. W. Benchmark Density Functional Theory Calculations for Nanoscale Conductance. J. Chem. Phys. 2008, 128, 114714. (51) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-Space Grid Implementation of the Projector Augmented Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035109−035119. (52) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dułak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; et al. Electronic Structure Calculations with GPAW: a Real-Space Implementation of the Projector Augmented-Wave Method. J. Phys.: Condens. Matter 2010, 22, 253202−253225. (53) Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Localized Atomic Basis Set in the Projector

Augmented Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 195112−195121. (54) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (55) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396−1396. (56) Tremblay, J. C.; Klamroth, T.; Saalfrank, P. Time-Dependent Configuration-Interaction Calculations of Laser-Driven Dynamics in Presence of Dissipation. J. Chem. Phys. 2008, 129, 084302−084309. (57) Tremblay, J. C.; Carrington, T., Jr. Using Preconditioned Adaptive Step Size Runge-Kutta methods for Solving the TimeDependent Schrödinger equation. J. Chem. Phys. 2004, 121, 11535− 11541. (58) Hermann, G.; Pohl, V.; Tremblay, J. C.; Paulus, B.; Hege, H.-C.; Schild, A. ORBKIT: A Modular Python Toolbox for Cross-Platform Postprocessing of Quantum Chemical Wavefunction Data. J. Comput. Chem. 2016, 37, 1511−1520. (59) Stalling, D.; Westerhoff, M.; Hege, H.-C. The Visualization Handbook; Elsevier, 2005; pp 749−767. (60) Kokalj, A. XCrySDen−a New Program for Displaying Crystalline Structures and Electron Densities. J. Mol. Graphics Modell. 1999, 17, 176−179. (61) Zhang, T.; Heid, R.; Bohnen, K.-P.; Sheng, P.; Chan, C. T. Phonon Spectrum and Electron-Phonon Coupling in Zigzag Graphene Nanoribbons. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 62− 67.

L

DOI: 10.1021/acs.jpcc.6b09682 J. Phys. Chem. C XXXX, XXX, XXX−XXX