A Multiscale Approach - ACS Publications - American Chemical Society

Subscriber access provided by UNIVERSITY OF SOUTH CAROLINA LIBRARIES

Article

Polaron Effects on Charge Transport through Molecular Wires: A Multi-scale Approach Tomas Kubar, Marcus Elstner, Bogdan Popescu, and Ulrich Kleinekathöfer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00879 • Publication Date (Web): 01 Dec 2016 Downloaded from http://pubs.acs.org on December 6, 2016

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

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

Page 1 of 34

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

Journal of Chemical Theory and Computation

Polaron Effects on Charge Transport through Molecular Wires: A Multi-scale Approach Tomáˇs Kubaˇr,∗,† Marcus Elstner,‡ Bogdan Popescu,¶,§ and Ulrich Kleinekathöfer¶ Institute of Physical Chemistry & Center for Functional Nanostructures, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany, Institute of Physical Chemistry, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany, and Department of Physics and Earth Sciences, Jacobs University Bremen, Campus Ring 1, 28759 Bremen, Germany E-mail: [email protected]

Abstract Modeling charge transport along molecular wires immersed in polarizable environments poses a grand challenge due to the high dimensionality of the problem and the various time scales involved. A previous multi-scale non-equilibrium Green’s function simulation scheme (Phys. Rev. Lett. (2012) 109, 176802) has been extended significantly, so that the present approach provides a much more complete physical description of the process. While the previous scheme involved the environmental fluctuations and their influence of the electronic structure of the wire, several previously neglected effects were added to the formalism: the electric field between the leads, the polarization of the dielectric environment in response to the charge present on the wire, and ∗

To whom correspondence should be addressed IPC & CFN, KIT ‡ IPC, KIT ¶ Jacobs University § Current address: Max-Planck-Institute for the Physics of Complex Systems, N¨ othnitzer Str. 38, 01187, Dresden, Germany †

1

ACS Paragon Plus Environment


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

the relaxation of the electronic structure of the wire. Still, the underlying Hamiltonian of the wire is evaluated with electronic structure calculations, and the dynamics of the molecular system are described using molecular dynamics simulation, so that (i) the formalism remains free of any model parameters, and (ii) no assumption on the underlying transport mechanism are being made. All the newly introduced details prove to affect the charge transfer along the wire markedly, while interestingly their effects compensate each other partially. The new method is suitable for application to charge transport in junctions composed of well defined molecular fragments, which is the case, e.g., in typical organic electronics materials. In this work, the method has been applied to hole transport through a double-stranded DNA, which nicely displays the influence of all the newly introduced effects.

1

Introduction

Charge transport (CT) in organic materials has been in the spotlight of experimental as well as theoretical research, on length scales ranging from nanometers for molecular junctions, to micrometers for bulk polymers and organic crystals. Molecular wires have attracted considerable interest not only because of being on a frontier between molecular electronics and organic electronics, but also because a transition between different transport mechanisms is expected to occur. 1–6 In a futuristic view, the emerging field of molecular electronics – the use of single molecules as wiring elements in electronic devices – may pass beyond the physical limits of conventional field effect transistor design, opposing further miniaturization. Relevant systems where different CT mechanisms appear include self-assembled monolayers of oligomers, single oligomers, and π-stacked systems such as DNA. 5,6 While charge transfer and transport processes are ubiquitous in biological systems, 7 it is double-stranded DNA that has earned special attention in physical and chemical research of CT. Due to the π-stacked arrangement of nucleobase pairs, dsDNA has been employed in numerous studies as a material constituting one-, two- and three-dimensional nanostructures. 8 From the mech2


Page 2 of 34

Page 3 of 34

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


anistic point of view, it has recently been recognized that the picture of separated coherent and incoherent processes is incomplete for dsDNA, and complex non-adiabatic electronphonon interactions could be the leading process. 5,9 Moreover, the emergence of molecular electronics has spurred the interest in DNA molecules as active wiring elements. 10 For these reasons, model dsDNA-based systems are considered as showcase applications in this work. A broad range of experimental and theoretical studies have been performed to understand the details of DNA charge transfer. 5,8,11 Particularly the issue of length dependence and the connected question of transfer mechanism have inspired a multitude of theoretical investigations. 12–18 To explain the short-range transfer, in which the resistance increases exponentially with length, a tunneling mechanism is usually assumed. Thus, quantum effects play a major role in this regime. 19 In case of long-range transfer, for which the resistance increases algebraically with length, a hopping mechanism is the main mode of action. A recent study, however, propagates an intermediate tunneling–hopping regime to be present for CT through DNA. 20 Although many theoretical studies are based on a time-independent picture of the system, it recently became clear that fluctuations in biomolecular charge transfer systems are usually large and non-negligible; 21,22 thermal fluctuations were proven to be of huge importance also in polymers 23 and organic semi-conductors. 24 To this end, multi-scale propagation approaches were recommended for the description of charge transfer in such systems. 25 Combinations of classical molecular dynamics (MD) simulations and quantum chemical calculations of ionization energies can be linked to quantum dynamical simulations in a hybrid approach to obtain the charge transfer or transport through DNA. 15 The transfer of charges through a DNA molecule from one end to the other end of the molecule is slightly easier to model since only a single charge is moving. 26–29 In case of charge transport, however, the DNA needs to be linked to electrodes making a constant flow of charges possible. 30,31 In a multi-scale approach the dynamics of the quantum dynamics onto the classical MD simulation shows significant effects for charge transfer systems 27 but also for excitation energy

3



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

transfer complexes. 32,33 So far, such a projection from the quantum to the classical system has not been tested for transport in a molecular junction, although a major effect is expected based on previous transfer simulations. 27 Thus, the aim of the present investigation is to extend and combine the previous work on charge transfer in DNA including the projection 27 and the studies on charge transport in a molecular junction without projection. 31 No assumptions regarding the mechanism of CT itself are necessary. Moreover, no additional parameters are needed due to the employed multi-scale approach connecting the quantum transport calculations to molecular dynamics simulations and electronic structure calculations. In a previous study by the present authors, 31 a MD simulation of a solvated DNA molecule (without leads) was performed, ionization energies and charge transfer couplings were determined along this trajectory using the Density Functional Tight Binding (DFTB) method, and this information was forwarded into a quantum transport calculation. In this way, hole transport through a fluctuating DNA molecular junction was described, and this approach is extended and improved in several respects in the present study. In the beginning, the MD representation of the system shall be improved by including the contacts with leads. The electronic parameters are determined along an MD simulation performed with such a setup and fed into a non-equilibrium Green’s function (NEGF) transport calculation; the charge distribution in the system does not influence the MD simulation, yet. We denote this set of simulations as NNN (No electric field in MD and electronic structure calculations, No projection of quantum charges onto MD systems, No 2nd-order energy correction in DFTB). As a first improvement, the respective bias potentials are applied to the leads, and the corresponding electric field is included in both the MD simulation and the calculation of site energies, leading to a simulation setup to be denoted as FNN (electric Field in MD and electronic structure calculations, No projection of quantum charges onto MD systems, No 2nd-order energy correction in DFTB). In the next step, the charge distribution yielded by the quantum transport calculation is projected onto the classical MD simulations, giving the

4


Page 4 of 34

Page 5 of 34

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


FPN setup (electric Field in MD and electronic structure calculations, Projection of quantum charges onto MD systems, No 2nd-order energy correction). Finally, a second-order correction to the calculations of electronic parameters is introduced to consider the relaxation of the electronic structure, in a setup denoted as FP2. This latter approach is the most advanced one developed in the present study. The four different levels of approximation NNN, FNN, FPN, and FP2 are tested numerically for an example system of a DNA junction containing several nucleobase pairs. The results obtained are discussed in detail concerning the site energies of the nucleobases, the excess charges on the individual nucleobases as well as the net current passing through the molecular junction.

2 2.1

Methods Classical molecular dynamics part of the simulation

The molecular example system under investigation consists of a polyA DNA sequence with 2, 4, 6 or 8 A–T base pairs in aqueous solution. This DNA molecule is placed in between two parallel gold plates. In the two dimensions parallel to the gold plates, the system is assumed to be periodic in the classical simulations while it is terminated beyond the gold plates. For the MD simulations, the AMBER BSC0 force field 34 is employed for the DNA, together with the TIP3P water model. The representation of gold is more challenging. To this end, the recently introduced GolDNA-AMBER approach is applied. 35 Thus, the surface polarization effect is described through a dipole located in the position of each gold atom. One such dipole constitutes of two opposite charges, connected by a rigid rod free to rotate around one end. When a molecule lies above the surface, the partial atomic charge on its atoms forces the rods to preferentially rotate in a specific orientation, modifying the average dipole moment. This model was shown to be reasonably accurate and computationally efficient at the same time. The gold atoms are held frozen during the simulation, therefore it is not 5



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

necessary to define an interaction between them. GolDNA-AMBER was parametrized with a multi-step approach based on DFT reference data, as described in the original publication. 35 In GolDNA-AMBER, the topmost layer of the surface slab includes real gold atoms and virtual interaction sites. Virtual sites are placed in hollow positions of the real lattice and are the ones interacting with adsorbing species via the LJ terms, while the LJ parameters of real gold atoms are set to zero. Dipolar rods are located only on the gold atoms. This fact is important for the present study, because also DNA bases prefer top adsorption sites, 36 although they are characterized by a pronounced lateral mobility. 37 The GolDNA-AMBER force field reproduces both these observations. For more details concerning the GolDNAAMBER parameters the interested reader is referred to the original publication. 35 Besides the intermolecular forces, additional forces arise by virtue of the applied electric field between the parallel gold electrodes. To take this effect into account, a constant voltage is applied in the MD simulation as well as in the electronic structure calculation, as external electric field. The intensity of this field is assumed to be constant, which is a standard approach, e.g., in the modeling of ion transport across nanopores. 38 Thus the electric potential is decreasing along the DNA linearly. This FNN scheme with electric field is the first improvement over the NNN approach, which in turn is similar to the earlier work in Ref. 31 where no explicit contacts were included in the MD setup. All MD simulations were performed with a local version of the Gromacs 4.6 package. 39,40 To this end, the dynamics of the molecular system was obtained integrating Newton’s equations of motion using the leap-frog scheme with a time step of 1 fs, while any bond involving a hydrogen atom was kept at constant length using the LINCS algorithm. Furthermore, the temperature was kept constant at 300 K using the Nosé–Hoover thermostat while the simulations were performed at constant volume (NVT canonical ensemble). Periodic boundary conditions were employed in two dimensions only, and the system was terminated beyond the gold using the wall functionality in Gromacs. Correspondingly, the long-range electrostatic interations were evaluated using a two-dimensional particle–mesh Ewald implementation.

6


Page 6 of 34

Page 7 of 34

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


2.2

Coupling of the classical and quantum systems

A crucial point of the current development is that Newton’s equations of motion for the atomic coordinates and the NEGF equations for the charge density are propagated simultaneously and coupled. A step of MD simulation of 1 fs is followed by a step of electron dynamics of 1 fs. The electron dynamics calculation proceeds in sub-steps of flexible length that is adjusted automatically to guarantee a pre-set level of accuracy. Then, another step of MD simulation follows, and so on. To mediate the interaction between the excess charge in the QM region and the atomic point charges in the MM region, the MM point charges qi on atoms, which are also modeled in the QM region, are updated in every step of the MD simulation. Performed is a linear interpolation between two sets of atomic charges, determined beforehand for a neutral nucleobase (qi0 , standard force-field charges) and a nucleobase carrying a full excess charge (hole or electron, qi+ ), qi = qi0 + Qi · (qi+ − qi0 ) .

(1)

The charges Qi are described in more detail in the next subsection and vary between 0 (neutral nucleobase) and 1 (nucleobase carries full excess charge). During the MD simulation, this procedure covers the polarization response of the MM environment due to the timedependent electron density on the CT-active nucleobases. This is the second improvement introduced in the current work, leading to the FPN simulation setup.

2.3

Electronic structure calculations

Our approach to describe the electronic structure of DNA has been designed for simulations of electron transfer processes in a chemical context, 27–29,41 and will take advantage of the methods developed in that context. An efficient description of the time-dependent electronic structure of DNA is based on several elements: (i) A QM/MM partitioning is applied in 7



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

Page 8 of 34

which the QM part includes those nucleobases that are active in CT, while the remainder of the molecular system is described using a MM force field. (ii) The QM component itself is made very efficient due to, firstly, a linear scaling fragmentation technique, and secondly, the excess electron/hole wave function being built from the frontier orbitals of the CT-active nucleobases. (iii) Application of the approximative DFTB2 method (a.k.a. SCC-DFTB) 42 brings on an acceleration by a factor of 100 to 1000 with respect to standard ab initio methods without the need to make significant compromises in accuracy. The approach will be sketched in the following paragraphs while the interested reader is referred to Refs. 28 and 29 for more details. The total hole (electron) wave function Ψ is considered to be composed of the frontier orbitals of the combined system of purine nucleobases in the DNA oligonucleotide. In turn, these frontier orbitals are represented by superpositions of the HOMO orbitals φHO m of the individual nucleobases m. Thus, Ψ itself is a linear combination of φHO m ,

Ψ(t) =

X

am (t)φHO m .

(2)

m

For DNA, this is a rather good approximation since the QM system is built from weakly interacting fragments (nucleobases). Here, am (t) are time-dependent expansion coefficients, which determine the instantaneous wave function of the excess charge and thus the excess electron/hole density, via the ‘coarse-grained’ (CG) charges on the nucleobases, Qm = a∗m am . In a first step, the fragment molecular orbitals (FMOs) φHO m are obtained using DFTB2. A separate SCF calculation is performed for every nucleobase within a QM/MM formalism, so that linear scaling is achieved instead of the cubic scaling of DFTB, alleviating the computational effort for systems composed of a large number of fragments; but still the FMOs are polarized by virtue of the MM point charges representing the molecular environment (DNA backbone, other nucleobases, and the solvent). Then, the electronic Hamiltonian of ˆ 0 is represented in this FMO basis, leading to an m × m matrix if m purine the system H

8


Page 9 of 34

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


nucleobases are included in the CT pathway, D E ˆ HO 0 Hmn = φHO . H φ m 0 n

(3)

0 Note that the matrix elements Hmn are influenced or polarized by the molecular environment

due to the inclusion of QM/MM terms in the calculation of the FMOs. The fragmentation scheme was shown to work very well as long as the system naturally decomposes into molecular units, which act as individual charge carriers. We have recently analyzed the accuracy of the fragmentation scheme, comparing the energies of CT states built from the fragment orbitals to those obtained with a full atomistic basis set, see Fig. 2 in Ref. 29. The errors are sufficiently small to guarantee that the CT energetics be free of any fragmentation artifacts. Obviously, this fragmentation scheme does not work for systems where the CT states cannot be described as superpositions of FMOs. This is the case in solid state systems, or whenever the electron tunneling through individual bonds is investigated, like, e.g., in Beratan’s pathway model. 43 0 The CG molecular Hamiltonian Hmn is time-dependent through the motion of the nuclei. 0 For hole transport, the diagonal elements Hmm possess the meaning of the negative of instan0 taneous vertical ionization energies, which will be also called the site energies ε0m = −Hmm

in the following. Notably, the major part of its variation is accountable to the environmental fluctuations through the QM/MM coupling. 44 Importantly, it is these fluctuations that brings about the extreme fluctuations of instantaneous currents through DNA junctions observed in previous studies. 30,31,45 0 The Hamiltonian Hmn does not involve the excess electron/hole present in the system.

To include the effect of electronic relaxation, the total energy of the system with density ρ is expanded up to second order around the charge-neutral system with density ρ0 , and after

9



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

Page 10 of 34

several approximations, the total energy reads 27 E[ρ] ≈ E MM ±

X

0 a∗m an Hmn +C ·

mn

X

Qm Qn Γmn ,

(4)

mn

where the − sign corresponds to hole and the + sign to excess electron transport. E MM denotes the energy of the charge-neutral system represented by an MM force field, Γmm ≡ Um is the chemical hardness of nucleobase m (U = 5.666 eV/e2 for adenine), and Γmn = 1/Rmn represents the Coulomb interaction of the excess charge between the nucleobases m and n. Following Mantz et al. 46 and our previous work, 27 the constant factor C < 1 in the secondorder term corrects for the DFT self-interaction error leading to an over-delocalization of excess charge. This total energy expression provides Hamiltonian elements corrected for the effect of electronic relaxation, 0 Hmn = ±Hmn + δmn · C ·

X

Qk Γkm .

(5)

km

A straightforward interpretation of a charge-dependent diagonal element Hmm ≡ εm is that the site energy (ionization energy in case of hole transport) of a nucleobase carrying an excess charge will increase, thus a larger amount of energy is required to transport more charge onto a site that is already charged. The description of the relaxation of electronic structure described in this paragraph is the third major improvement introduced in this work, leading to the FP2 simulation setup. Note that the off-diagonal elements are charge-independent up 0 to the second order in density, Hmn = Hmn if m 6= n.

It should be also noted that the interaction of the excess charge with the MM environment is overestimated because of the missing description of electronic polarization in the used nonpolarizable MM force field. Here, the simplest way will be employed to take these effects into account, namely scaling down the electrostatic QM/MM interaction by a certain constant factor. This factor should correspond to the inverse square root of the optical dielectric constant of water, and a value of 1/1.4 will be used here. 29,47,48 10


Page 11 of 34

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


2.4

Propagation of the charge density

The aforementioned techniques yield site energies εm ≡ Hmm and couplings between site m and site n Tmn ≡ Hmn for a DNA molecule in a tight-binding model. The CT from one lead to the other through the molecule, however, needs to be determined using a quantum approach. Due to the strong and fast fluctuations triggered by the environment, only a limited set of quantum transport theories is applicable. Here we employ the same time-dependent nonequilibrium Green’s function scheme 49 used in our previous study. 31 This method is based on a parametrization of the spectral density describing the energy-dependent coupling between the molecule and the electrodes 50 as well as on a decomposition of the Fermi function in terms of poles in the complex plane. 51 These expansions enable one to express the reservoirs self-energies as a sum of weighted exponentials, which further allows for a convenient analytic treatment. 49,52 An alternative approach would use, e.g., a Chebyshev expansion technique. 53 Note that a clear advantage of the present scheme is that time-dependent effects, even strong in magnitude and/or triggered by heavily fluctuating quantities, can be treated without any further approximation. In general, the CT through the molecule coupled to electrodes can be studied by the (single-particle) reduced density operator σ(t). Its time evolution is given by an equation of the following form, 49,54

i¯h

X ∂σ(t) Πα (t) + Π†α (t) . = [H(t), σ(t)] + i ∂t

(6)

α∈{L,R}

The current matrices Πα (t) are given by another set of equations, 49 which use the above mentioned twofold decomposition of molecule–lead coupling and Fermi functions. The reader is referred to Refs. 49 and 52 for details. These current matrices Πα (t) are important quantities also to obtain the time-dependent current Iα (t) from or to the left (α = L) or right lead (α = R) using Iα (t) =

2e Re Tr {Πα (t)} . h ¯ 11


(7)


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

Usually, however, one is only interested in the net current through the molecular junction given by I(t) =

1 2

(IL (t) − IR (t)). We like to highlight that the present approach is not

assuming any transport mechanism such as pure coherent transport but is generally applicable and also includes explicit populations on the sites of the molecule. Moreover, the time dependence of the site energies leads to additional dephasing effects in the junction. 31 As a start for the NEGF theory, the Hamiltonian of the total system HT (t) has to be separated into a relevant system part H(t) describing the molecule and HR accounting for the reservoirs. Moreover, HSR is describing the wire-reservoir coupling, leading to HT (t) = H(t) + HR + HSR . Within second quantization and using the creation (annihilation) operators at site m denoted by c†m (cm ), one can describe the DNA molecule P PN † † by H(t) = N m=1 εm (t)cm cm − m6=n Tmn (t)cm cn . Here εm and Tmn are the site energies and couplings already introduced above, whereas N denotes the total number of molecular sites. Note that the explicit time dependence of the tight-binding parameters εm , Tmn stems directly from the electronic structure calculations (see previous section). Next, the fermionic reservoirs are described using the creation (annihilation) operators within the P † reservoirs b†αk (bαk ) by HR = αk εαk bαk bαk . In this expression εαk denotes the energy (†)

of a particle in state k of reservoir α. Note that the system cm and reservoirs operators (†)

bαk anticommute with each other. Finally, one needs to specify the wire–lead Hamiltonian α with the coupling constant Tkm between reservoir state k and device state m leading to P P α † α∗ † HSR = α km Tkm bαk cm + Tkm cm bαk . In particular, for the present study we made the

assumption that only the left-most site of the wire and the right-most one couple to the left and right electrode, respectively. For the present application of DNA between two gold electrodes this assumption is perfectly fine.

12


Page 12 of 34

Page 13 of 34

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


2.5

Adding electric field, projection on classical simulation and 2nd-order correction

Let us summarize the features added to the simulation setup, which bring on a more physical, more realistic description. In an earlier study we have determined in the CT dynamics through a DNA in a molecular junction geometry. 31 The leads as well as the electric field were not explicitly treated in the MD setup and, more importantly, there was no projection from the charge dynamics onto the classical dynamics, i.e., polaron effects of the environment were not included. These limitations are lifted here, step by step. NNN Different from the previous study, 31 gold contacts are now explicitly modeled in the MD setup. Other than this and a different DNA sequence, these two simulation approaches are the same. FNN The effect of the electric field between the electrodes can be taken into account in the equations of motion. This concerns two components of the simulation scheme: First, the field has to be included in the underlying MD simulation, inducing additional forces on atoms carrying non-zero charges (which often are all of the atoms of the (bio)molecule and the solvent). And second, the quantum-chemical calculation of site energies has to include an additional term describing the field as well, and this can be done easily with the DFTB calculations employed in the present investigation. FPN The electron density present on the molecular wire may be included in the dynamics in the same way as in our previous multi-scale simulations of biomolecular electron transfer: 27 The instantaneous charge population on every site can be mapped onto the MD atomic charges in every step of the leap-frog integration, by interpolating between two sets of potential-derived charges of the respective molecule (here, the nucleobases), as described in Eq. 1. Then, the molecular system will polarize according to the charge density on the wire; note that only orientational polarization is described with the

13



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

current setup because the molecular-mechanics force field does not include electronic polarization. FP2 The effect of electronic relaxation on the site energies may be described by passing formally to the second-order formalism in the CG DFTB model of the DNA wire, which results in the charge dependence of site energies as introduced in Eq. 5. Effectively, the site energies of the nucleobases carrying some amount of excess charge will increase and also, to a lesser extent, the site energies of neighboring nucleobases. While this is an integral component of our biomolecular electron transfer model, the previous NEGF formalism did not consider this electronic relaxation. Thus, two effects act on the site energies: the polarized molecular environment induced electric field that leads to a decrease of site energies (captured by including the QM/MM term in the DFTB calculation of site energies), and the relaxation of electronic structure results in their effective increase (due to Eq. 5). The contrary character of these effects are depicted in Fig. 1. We note that the NEGF ansatz would require modification in order to account for occupation-dependent site energies, leading to considerably increased computational cost. Thus, we decided to implemented a ‘poor man’s’ second-order treatment instead. To this end, every site energy is calculated according to the current charge populations on all sites prior to every step of the quantum propagation. These site energies (but not populations) are being held constant during that MD time step, which may be considerably longer than the flexible time step of the NEGF propagation. In this way, the dependence of the site energies on charge populations is taken into account in an approximate fashion keeping the computational overhead minimal.

14


Page 14 of 34

Page 15 of 34

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


Figure 1: Site energies of guanine nucleotides in a polyG dsDNA oligonucleotide; mean values from simulation described below (error bars designate the standard deviation). The nucleobase designated by (+) carries one full hole, and all others are uncharged. The electric field is switched off for the purpose of these calculations. The black line shows simulations with no effect of the hole on the MM environment and no relaxation of the electronic structure is considered (equivalent to FNN in this work). The effect of the hole on the MM environment during the dynamics is shown by the blue curve. The neighboring sites are effected by the polarization of environment as well (FPN). The red curve depicts simulations with MM polarization and, additionally, the effect of electronic relaxation considered in the calculation of the site energies according to Eq. 5 (FP2). The MM polarization and the electronic relaxation induce opposite effects with the effect of the MM polarization being larger in magnitude.

2.6

Choice of system parameters

As mentioned already, the polyA DNA sequence which is put in between two gold electrodes consists of 2, 4, 6 or 8 A–T base pairs. Adenine nucleobases are considered as hole carriers and they compose the QM region in the QM/MM calculations. Thus, the number of QM atoms is 30, 60, 90 and 120, respectively, and there are 90, 180, 270 and 360 basis functions with DFTB, respectively, in these systems. Not all different method combinations have been, however, studied for all DNA lengths. The most detailed analysis has been performed for four A–T base pairs. Certainly the transport current through the junction strongly depends on the applied voltage, at least in certain voltage ranges. We fix the difference of the electrode potentials

15



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

at 2 V, which is a value that occurs in some of the relevant experiments. 55 The magnitude of the electric field for both the force field calculation and the DFTB calculation is constant, and given as the voltage of 2 V divided by the distance between the two gold plates, which was 10.8 ˚ A, 18.0 ˚ A, 25.2 ˚ A and 31.4 ˚ A in the systems with 2, 4, 6 and 8 nucleobase pairs, respectively. The absolute Fermi energies of the electrodes deserve a second thought. The Fermi energy of gold lies at 5.5 eV, and a symmetric bias around this value would lead to a conduction window of 4.5–6.5 V. However, the HOMO energy of adenine obtained within DFTB (5.6 eV) is 2.6 eV lower than the ionization energy of adenine (8.2 eV). Therefore, the Fermi levels of the electrodes should also be decreased by the same amount to correct for the too low ionization energy of adenine obtained with the current methodology. The result is a pair of Fermi levels of 1.9 eV and 3.9 eV. Due to the rather approximative character of this estimation we will consider several pairs of electrode energies in the present investigation. The present propagation method for the charge density allows to determine the current for energy-dependent molecule–lead couplings. This energy dependence can be modeled with a sum of Lorentzian functions, which may be obtained with fitting to any energy-dependent coupling function, e.g., to represent the band structure of gold. At the same time, such a detailed modeling of the molecule–lead interaction would introduce additional parameters and would make the interpretation of results significantly more difficult. Thus, in order to study the effects of the above described enhancements of the theory, we refrained from using any energy-dependent molecule–lead coupling. Instead, a constant coupling value of 1 meV was considered following Ref. 31.

16


Page 16 of 34

Page 17 of 34

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


3

Results

3.1

Sequential coupling without explicit electric field (NNN)

In this first variant of the present scheme, a sequential coupling of the MD simulations, quantum chemistry calculations and CT dynamics is employed in order to mimic the results of Ref. 31 for the present MD setup. To this end, first the MD simulation of the molecular junction is performed without explicitly including the field effect in the setup. Subsequently, the ionization energies are determined along this MD trajectory without second-order correction and without taking the electric field into account in theses DFTB calculations. The site energies are only modulated by the fluctuations of the molecular environment with the known standard deviation of 0.30 eV. 44,56 Note that this fluctuation reduces to ca. 0.22 eV upon scaling of the electrostatic field induced by the environment by a factor of 1.4, which is intended to mimic the effect of electronic polarization. 47,48 The latter effect is not included in the present MD simulations with non-polarizable force fields. Table 1: Average currents, site populations and energies for sequential coupling without explicit electric field. leads / V I / nA Q1 Q2 Q3 Q4 ε1 / eV ε2 / eV ε3 / eV ε4 / eV

4.2–6.2 109.52 0.452 0.509 0.522 0.533 5.68 5.61 5.61 5.59

3.7–5.7 74.98 0.310 0.353 0.361 0.369 5.64 5.57 5.58 5.57

3.2–5.2 9.09 0.038 0.043 0.044 0.045 5.65 5.57 5.59 5.59

Similar simulations have been performed earlier with a different DNA sequence and without including leads in the MD part. 31 For a conduction window between 4.2–6.2 eV a large average current of 110 nA is found (see Tab. 1, averaged over a 700 ps simulation). The reason for this significant current is certainly that all site energies are well within the conduction window fluctuating around 5.6 eV. Moreover, all sites are roughly half-populated as can been seen from the NEGF calculations. Since spin is not included in the present calculations, the site populations can vary between zero and one. In case of shifting the conduction window 17



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

Figure 2: The first 200 ps of a simulation with sequential coupling without explicit electric field are shown for the case of conduction window 3.7–5.7 eV. From left: site energies (eV), occupations of sites (e), and current (nA). to 3.7–5.7 eV, the observed average current is reduced to 75 nA (see also Fig. 2). In this case the the site energies are already at the edge of the conduction window and leaving it due to the fluctuations part of the simulation time. At the same time, the average populations decrease to 0.3–0.4. The current is further reduced using a conduction window of 3.2–5.2 eV and drops drastically at 2.2–4.2 eV (not simulated). Overall, we can say that these results qualitatively agree with those from the earlier study with different system parameters. 31

3.2

Sequential coupling with explicit electric field (FNN)

The next set of simulations is identical to the above described sequential coupling of MD, quantum chemistry and dynamics but now includes the effect of the electric field also in the MD and quantum chemistry and not only in the NEGF transport calculations. There is a rather significant effect on the determination of the site energies since the electric field leads to a shift in these energies. Consequently, the averaged sites energies are not equal for all sites anymore, rather, they decrease monotonously from the site closest to the source to the one closest to the drain. This shift results in energies of the site energies fluctuating between 5.8 and 6.2 eV for the site closest to the source and between 5.1–5.4 eV for the site closest to the drain. When choosing the conduction window between 4.2 eV and 6.2 eV, all site energies fit into this range so that a reasonably large current of 42 nA is observed. The current is smaller 18


Page 18 of 34

Page 19 of 34

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


than the largest current in the previous setup indicating that the positioning of the site energies is not always optimal for CT. The average populations in the present case increase from drain to source as 0.17–0.25–0.32–0.36. Using other choices of Fermi energies for the electrodes, the observed range of site energies moves gradually out of the conduction window. This fact leads to decreased populations on the sites as well as a decreased current. The average current is only 4 nA with a conduction window 3.7–5.7 eV and it nearly vanishes for a conduction window 3.2–5.2 eV (see Tab. 2). Table 2: Average currents, site populations and energies for sequential coupling with explicit electric field. leads / V I / nA Q1 Q2 Q3 Q4 ε1 / eV ε2 / eV ε3 / eV ε4 / eV

4.2–6.2 41.63 0.173 0.251 0.323 0.362 5.05 5.37 5.79 6.17

3.7–5.7 4.42 0.018 0.027 0.032 0.032 5.08 5.44 5.81 6.16

3.2–5.2 0.07 0.000 0.000 0.000 0.001 5.04 5.37 5.80 6.18

It is noteworthy that while the course of the electron dynamics seems to achieve a certain equilibrium on a long time scale, significant fluctuations can be observed on short time scales. For instance, in the simulation with conduction window 3.7–5.7 eV (see Fig. 3), the occupations of the sites are 0.18–0.27–0.32–0.32 on average, but these occupations fluctuate in a rather broad range as do the site energies. Moreover, this dynamics leads to the rather complex time dependence of the current. Remarkably, the time series of the current do not possess this character in the above simulations without explicit electric fields or as in Ref. 31 in which the site energies fluctuate around a common average value.

3.3

Simulations including the projection from the quantum onto the classical dynamics (FPN)

These simulations are the first in the series where the MM molecular environment and its dynamics are coupled to the dynamics of the electronic system. Thus, not only does 19



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

Figure 3: The first 200 ps of a FNN simulation (sequential coupling with explicit electric field) are shown for the case of conduction window 3.7–5.7 eV. From left: site energies (eV), occupations of sites (e), and current (nA). The scale for site energies is the same as in Fig. 2, while this is not the case for the occupations and for the current. the electronic system feel the molecular structure via the time-dependent electron-transfer Hamiltonian, but in addition, the atomic charges of the adenine nucleobases in the MM model are varied according to the instantaneous population of each respective site. These charges are adopted at every MM time step of the coupled simulation. Consequently, the site energies will vary in much broader intervals than in the previous simulations above. The reason is the strong polarization of the solvent by the charge present in the system. Note that we use the quantity ‘population’ to describe the amount of charge present. In our earlier work, we have quantified the strength of the polarization of environment, in terms of the electric potential induced on the DNA molecule. We have found that when one electron hole is located on a nucleobase, then a potential of as much as −3 V is induced on that molecule, while a smaller additional potential is observed on the neighboring nucleobases, see also Fig. 1. Note further that the magnitude of this effect is somewhat smaller, as it should be corrected for the missing description of the electronic polarization of the environment. Reasonable values may be obtained using a factor of 1/1.4, so that the additional potential is approximatively −2 V. The results from simulation including the projection of charge density onto the classical MD are listed in Tab. 3. Actually the results are quite different to the previous results in that either all sites get charged completely or all get discharged completely. Either way, the

20


Page 20 of 34

Page 21 of 34

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


resulting net current is zero. In the simulations with conduction windows 4.2–6.2 eV, 3.7– 5.7 eV, and 3.2–5.2 eV all sites are completely charged after some initial transient dynamics. For conduction windows 2.3–4.2 eV and 1.2–3.2 eV all site become completely discharged. Table 3: Average currents, site populations and energies for sequential coupling with explicit electric field. leads / V I / nA Q1 Q2 Q3 Q4 ε1 / eV ε2 / eV ε3 / eV ε4 / eV

4.2–6.2 0.00 1.000 1.000 1.000 1.000 0.75 0.27 0.78 2.00

3.7–5.7 0.00 1.000 1.000 1.000 1.000 0.76 0.28 0.79 2.00

3.2–5.2 0.00 1.000 1.000 1.000 1.000 0.75 0.28 0.76 2.02

2.2–4.2 0.00 0.000 0.000 0.000 0.000 5.04 5.40 5.81 6.16

1.2–3.2 0.01 0.000 0.000 0.000 0.000 5.04 5.39 5.79 6.15

For an enhanced understanding of these results it is interesting to look at the mechanism of the complete charging or discharging of the wire. Note that all simulations were preceded by an equilibration MD run performed with the nucleobases occupied by 0.14–0.38–0.62–0.86 holes, and the actual electron dynamics simulations (FPN) are started with these occupations forming the diagonal of the density matrix. The first picoseconds of the simulations including projection with conduction windows 4.2–6.2 eV and 1.2–3.2 eV are presented in Fig. 4 in terms of site energies and site populations. In case of the 4.2–6.2 eV conduction window, the site energies lie at the lower edge of the conduction window. The site in contact with the drain electrode is even below that value. Thus, the charge flows into the system and cannot leave. As the sites are being charged gradually, the polarization of the environment becomes simultaneously stronger step by step resulting in decreasing site energies. Consequently, the site energies are lowered even further diminishing any possibility for the site energies to come closer to the Fermi energy of the drain. The exactly opposite scenario happens in case of the conduction window 1.2–3.2 eV. The site energies are at or above the upper edge of the conduction window, initially. Therefore, hardly any charge is flowing into the system, whereas the charge that is present in the system is leaving the system into the drain electrode. As the sites are being discharged, the polarization of the environment is becoming gradually weaker resulting in increasing site energies. Thus, the energy difference between 21



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

the source and the sites is increasing, making it even less likely for any charge density to flow onto the wire. Thus, for this second example with conduction window starting at a lower energy again a positive feedback is found although it now leads to a complete discharging of the sites.

Figure 4: The first 20 ps of FPN simulations (sequential coupling with explicit electric field and dynamic polarization of molecular system due to the excess charge). The Fermi energies of leads set in two different ways: left, 4.2–6.2 eV, and right, 1.2–3.2 eV. Top, site energies (eV), and bottom, occupations of sites (e).

Furthermore, the minimum of the conduction window, i.e., the Fermi level of the source electrode, has been shifted in small steps in order to identify a position that would warrant stable current. The results are listed in Tab. 4. No such Fermi level was found while keeping the bias voltage at 2 V. Thus, in all situations tested the projection of the quantum dynamics onto the classical MD dynamics led to a complete charging or discharging of all sites. In the following, we will show that including an additional correction will cure this issue.

22


Page 22 of 34

Page 23 of 34

Table 4: Average currents, site populations and energies for sequential coupling with explicit electric field. ∗ Q1 changes suddenly from 0 to 1 about half-way through the simulation. left lead / V I / nA Q1 Q2 Q3 Q4 ε1 / eV ε2 / eV ε3 / eV ε4 / eV

2.2 0.00 0.000 0.000 0.000 0.000 5.04 5.40 5.81 6.16

2.225 0.00 0.000 0.000 0.000 0.000 5.04 5.39 5.79 6.14

2.25V 0.00 0.000 0.000 0.000 0.000 5.05 5.40 5.80 6.16

2.2525 0.02 0.009 0.994 0.996 1.000 2.95 1.59 1.39 2.27

2.255 0.02 0.012 0.994 0.997 1.000 2.93 1.57 1.37 2.28

2.2575 0.00 0.000 0.000 0.000 0.000 5.04 5.39 5.78 6.15

2.26 0.00 0.000 0.000 0.000 0.000 5.06 5.40 5.79 6.15

2.2625∗ −0.13 0.578 0.997 0.999 1.000 1.69 0.84 1.04 2.13

2.275 −0.13 0.990 1.000 1.000 1.000 0.78 0.27 0.76 2.00

2.2875 0.00 1.000 1.000 1.000 1.000 0.78 0.30 0.78 2.01

2.3 −0.13 0.869 0.998 0.999 1.000 1.04 0.46 0.88 2.05

Table 5: FP2 results for DNA with 4 nucleobases; the simulation sets were performed with different initial populations: Left – from gradually increasing populations, 0.2–0.8; center – from zero populations for all sites; right – from populations of 0.5 for all sites. 23

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


leads / V I / nA Q1 Q2 Q3 Q4 ε1 / eV ε2 / eV ε3 / eV ε4 / eV

initial gradually increasing populations 4.2–6.2 3.7–5.7 3.2–5.2 2.2–4.2 1.2–3.2 0.00 87.15 0.07 0.00 0.00 1.000 0.407 0.000 0.000 0.000 1.000 0.437 0.000 0.000 0.000 1.000 0.514 0.001 0.000 0.000 1.000 0.607 0.001 0.000 0.000 0.76 3.14 5.04 5.04 5.05 0.28 2.98 5.39 5.37 5.39 0.77 3.20 5.78 5.79 5.80 2.00 3.85 6.15 6.15 6.16

4.2–6.2 0.00 1.000 1.000 1.000 1.000 0.76 0.27 0.74 1.97

initial populations of zero 3.7–5.7 3.2–5.2 2.2–4.2 53.59 0.06 0.00 0.656 0.000 0.000 0.677 0.000 0.000 0.690 0.000 0.000 0.758 0.001 0.000 2.16 5.04 5.04 1.86 5.40 5.39 2.21 5.80 5.80 3.10 6.15 6.16


1.2–3.2 0.01 0.000 0.000 0.000 0.000 5.05 5.41 5.82 6.18

4.2–6.2 0.00 1.000 1.000 1.000 1.000 0.76 0.24 0.70 1.97

initial populations of 0.5 3.7–5.7 3.2–5.2 2.2–4.2 36.49 0.06 0.00 0.755 0.000 0.000 0.771 0.000 0.000 0.780 0.000 0.000 0.834 0.001 0.000 1.76 5.03 5.03 1.41 5.39 5.38 1.80 5.81 5.79 2.78 6.16 6.14

1.2–3.2 no data


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

3.4

Simulations including the second-order term, projection of charge density as well as the effect of the electric field (FP2) – the ‘most complete’ model

In the final step, the physical model is expanded once more to involve the relaxation of the electronic structure due to the presence of excess charges in the system. This effect is described by the variation of the site energies as functions of charge present on the nucleobases, see Eq. 5. Note that the polarization of the environment (via the electric field that it induces on the nucleobases) and the relaxation of the electronic structure affect the site energies in opposite directions. Thus, these two phenomena compensate each other to a certain extent while the effect of polarization prevails as shown in Fig. 1. A possibly important technical detail is the self-interaction correction represented by the scaling factor C in Eq. 5. This extremely simplified correction possesses a rather effective but empirical character and thus it is difficult to assign proper values to the constant C. In earlier work on biomolecular electron transfer we considered C = 0.2, 27 in analogy to the original proposal for full-DFT calculations. 46 To assess the influence of the choice of C, several values were scanned in these FP2 simulations, namely 0.15, 0.2, 0.25, 0.3, 0.35 and 0.4, see the Supporting Information, Tab. S2. Qualitatively the same trends were observed in all these simulations The discussions below will only consider simulations performed with a value C = 0.2, which is the same value as used in the previous simulations of (chemical) electron transfer in DNA and proteins performed by part of the present authors. 29 Using a conduction window of 3.7–5.7 V, site occupations in the range of 0.4–0.6 are established (see Fig. 5). These occupations remain stable along the entire simulation trajectory. Larger fluctuations of site energies, howver, occur quite often driving the occupations of sites out of the steady state configuration. Nevertheless, the system always returns into that state. As a result, the site energies are 3.1–3.0–3.9 eV on average, which is at the lower edge of the conduction window leading to a relatively high average current of 87 nA.

24


Page 24 of 34

Page 25 of 34

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


Figure 5: From left: current (A), occupations and site E (eV) in FP2 simulations of 4 sites with lead Fermi energies 3.7–5.7 V; time in ps. Moreover, different initial conditions concerning the charge populations along the wire were tested and the results listed in Tab. 5. The above mentioned observations originate from simulations with gradually increasing initial occupations of 0.2–. . . –0.8. Additional simulations were performed with (i) all occupations equal to zero, and (ii) all occupations equal to 0.5. The same observations were made in these simulations: After some initial transient period of about 15 ps, occupations of the sites were established that were very similar to those in the above mentioned simulations. The obtained currents were of the same order of magnitude as well. Thus, it seems that the most complete model, FP2, is capable of supporting a selfequilibrating property of the DNA wire. In this context, the phenomena that affect the site energies have to be identified: (i) On the one hand, it is the electrostatic potential induced by the polarizable molecular environment, aqueous solution in this case. The presence of some charge population on a site leads to the polarization of the environment, which in turn induces an electric field on that site as well as its neighboring sites. This effect was not present in the previous simulations 31 and was first implemented in this work. The induced potentials lead to a decrease of site energies (for the transport of positive charge, hole transport). (ii) On the other hand, an increased charge population of a site has to be considered as a perturbation of the electronic density of that site. Consequently, the electronic structure of the molecule relaxes, and a prominent manifestation is a rather large shift of orbital energies with which the site energies in our scheme are identified. Therefore, 25



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

the site energy of a molecule carrying charge increases (for the transport of positive charge). The effects (i) and (ii) drive the site energies in opposite directions. This feature may be crucial for supporting the system to reach a steady state. Furthermore, we performed FP2 simulations for DNA oligonucleotides containing 2, 6, and 8 adenines between gold electrodes, see Tab. S1. The observations were similar in all simulations to those with four adenines presented above. Thus, with a conduction window of 3.7–5.7 V, all DNA oligonucleotides support a current at the order of tens of nA. If the Fermi energies are shifted up or down, then lower or no current were observed because the sites get charged or discharged entirely. The former leads to extremely low site energies, far below the drain. The latter leads to high site energies, far above the source. In both cases, no reasonable current is supported. Let us briefly discuss the current implementation of the relaxation of the electronic structure. The term to the second order in electron density shall describe the repulsion of excess charge present on the molecular wire. In the application to electron transfer, where this development originates, there is always exactly one excess electron or hole in the system, thus it is both necessary and meaningful to treat the second-order term with a self-interaction correction, at least, in the effective way with the scaling factor C in Eq. 5. Now, in the current study of CT, the amount of excess electron density present in the system is not constrained to correspond to exactly one electron. There may be a larger amount of density, i.e., a larger number of electrons, and then the interaction between the individual sites may not be considered to consist of self-interaction between partial densities corresponding to a single electron. Then, it may not be correct to scale down the second-order term, because two different electrons should repel each other fully, i.e., with C = 1. This missing repulsion is possibly the reason that the simulations FPN are inaccurate. This potentially too weak repulsion among different electrons composing the excess density on the wire will likely affect the balance of the occupations of individual sites. This balance will, strictly speaking, only be reached as long as two different electrons interact fully, and only the fractional occupa-

26


Page 26 of 34

Page 27 of 34

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


tions that can be ascribed to a single electron are treated with the self-interaction correction. Still, as we have tested above and shown in Tab. S2, the effect of different choices of C is too weak to change the above presented qualitative findings.

4

Conclusions

Despite the large amount of previous experimental and theoretical investigations concerning the CT in molecular junctions, a complete theoretical picture of all the dynamical effects which are involved is still missing. The present study represents a significant expansion of a previous multi-scale simulation scheme. 31 A key advantage of the current approach is that no assumption on the underlying transport mechanism are being made, and both hopping-like and coherent transport are possible, possibly depending on the length of the wire and, for example, the nucleobase sequence of the specific DNA molecule. 12–14,17 Recently, a tunnelinghopping regime of CT through DNA was proposed, 20 and all of these different transport mechanisms can be handled by the non-equilibrium Green’s function approach employed in the present investigation. Moreover, this scheme is capable of treating fluctuations in the underlying coarse-grained Hamiltonian, which are due to the dynamics of the molecular environment prominently, as established previously. 44,56 The relevant time-dependent coarsegrained Hamiltonian itself is being determined on the basis of fully atomistic molecular dynamics simulations and electronic structure calculations, making the formalism free of any model parameters. The present approach involves several enhancements with respect to previous studies. Contrary to Ref. 31, the structural effects of the gold plates are included in the MD system setup properly, although this does not affect conductances significantly. The inclusion of electric field effects in the MD simulation and, especially, in the electronic structure calculation leads to a significant variation of the site energies, and this influences conductances of the DNA junction strongly. All the above calculations used a sequential coupling of an

27



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

MD simulation and quantum chemistry calculation followed by a transport calculation. A significant step towards a more realistic description is the introduction of a feedback loop, i.e., the projection of the charge density from the transport calculation onto the MD simulation. Only then is it possible for the molecular environment to re-polarize in response to the time-dependent charge density of the molecular junction, and this re-polarization in turn induces a huge decrease of the site energies of nucleobases. This effect is actually so strong that the net current through the junction vanished for all example cases in this work, unless the electronic relaxation of nucleobases was considered. The effective description of this relaxation, which is provided by the second-order correction within DFTB, results in an increase of the site energies, thus counter-acting the effect of environmental polarization to some extent. Relatively large currents are recovered as a consequence. Evidently, there is a complex balance and counterplay of the individual effects participating in the control of the CT process. Introducing individual corrections while omitting others may thus lead to misleading results. In summary, several restrictions present in earlier studies of CT in molecular junctions are lifted by the current multi-scale approach. It allows to describe the transport without assuming any transport mechanism, hopping or coherent, a priori. This formalism can certainly also be applied to similar systems, i.e., other molecules in molecular junctions. Further developments of the method will concentrate on a proper calculation of the electronic couplings between the molecule and the leads, and on an acceleration of the transport calculation, which is the computational bottle neck presently. For instance, the performance of the transport part could be significantly improved by using the more efficient new version of the time-dependent NEGF scheme of Ref. 49, which was recently reported in Ref. 57.

Acknowledgement This work has been supported by the Deutsche Forschungsgemeinschaft (DFG), joint grant EL 206/13-1 & KL 1299/14-1. 28


Page 28 of 34

Page 29 of 34

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


Supporting Information Available Results from FP2 simulations performed for DNA oligonucleotides containing 2, 4, 6, and 8 adenines; results from FP2 simulations performed with different values of the scaling parameter C.

This material is available free of charge via the Internet at http://pubs.

acs.org/.

References (1) Yan, H.; Bergren, A. J.; McCreery, R.; Della Rocca, M. L.; Martin, P.; Lafarge, P.; Lacroix, J. C. Proc. Natl. Acad. Sci. USA 2013, 110, 5326–5330. (2) Zhao, X.; Huang, C.; Gulcur, M.; Batsanov, A. S.; Baghernejad, M.; Hong, W.; Bryce, M. R.; Wandlowski, T. Chem. Mater. 2013, 25, 4340–4347. (3) Hong, W.; Li, H.; Liu, S.-X.; Fu, Y.; Li, J.; Kaliginedi, V.; Decurtins, S.; Wandlowski, T. J. Am. Chem. Soc. 2012, 134, 19425–19431. (4) Choi, S. H.; Risko, C.; Delgado, M. C. R.; Kim, B.; Brédas, J.-L.; Frisbie, C. D. J. Am. Chem. Soc. 2010, 132, 4358–4368. (5) Genereux, J. C.; Barton, J. K. Chem. Rev. 2010, 110, 1642–1662. (6) Ratner, M. Nature 1999, 397, 480–481. (7) Blumberger, J. Chem. Rev. 2015, 115, 11191–11238. (8) Kawai, K.; Majima, T. Acc. Chem. Res. 2013, 46, 2616–2625. (9) Grozema, F. C.; Tonzani, S.; Berlin, Y. A.; Schatz, G. C.; Siebbeles, L. D. A.; Ratner, M. A. J. Am. Chem. Soc. 2008, 130, 5157–5166. (10) Cuniberti, G., Fagas, G., Richter, K., Eds. Introducing Molecular Electronics; Lecture Notes in Physics; Springer: Berlin, 2005. 29



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

(11) Wagenknecht, H.-A., Ed. Charge Transfer in DNA: From Mechanism to Application; Wiley: Weinheim, 2005; 978-3-527-31085-2. (12) Berlin, Y. A.; Burin, A. L.; Ratner, M. A. J. Am. Chem. Soc. 2001, 123, 260–268. (13) Bixon, M.; Jortner, J. Chem. Phys. 2002, 281, 393. (14) Siriwong, K.; Voityuk, A. A. WIREs Comput. Mol. Sci. 2012, 2, 780–794. (15) Kubaˇr, T.; Gutiérrez, R.; Kleinekathöfer, U.; Cuniberti, G.; Elstner, M. Phys. Status Solidi B 2013, 250, 2277–2287. (16) Li, G.; Govind, N.; Ratner, M. A.; Cramer, C. J.; Gagliardi, L. J. Phys. Chem. Lett. 2015, 6, 4889–4897. (17) Levine, A. D.; Iv, M.; Peskin, U. Chem. Sci. 2016, 7, 1535–1542. (18) Liu, C.; Beratan, D. N.; Zhang, P. J. Phys. Chem. B 2016, 120, 3624–3633. ˇ aˇc, J.; Lévy, B.; Sanders, B. C.; Salahub, D. R. (19) de la Lande, A.; Babcock, N. S.; Rez´ Phys. Chem. Chem. Phys. 2012, 14, 5902–5918. (20) Xiang, L.; Palma, J. L.; Bruot, C.; Mujica, V.; Ratner, M. A.; Tao, N. Nat. Chem. 2015, 7, 221–226. (21) Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Annu. Rev. Phys. Chem. 2010, 61, 461–485. (22) Beratan, D. N.; Liu, C.; Migliore, A.; Polizzi, N. F.; Skourtis, S. S.; Zhang, P.; Zhang, Y. Acc. Chem. Res. 2015, 48, 474–481. (23) Troisi, A.; Cheung, D. L.; Andrienko, D. Phys. Rev. Lett. 2009, 102, 116602. (24) Fogel, Y.; Zhi, L.; Rouhanipour, A.; Andrienko, D.; Räder, H. J.; M¨ ullen, K. Macromolecules 2009, 42, 6878–6884. 30


Page 30 of 34

Page 31 of 34

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


(25) Cheung, D. L.; Troisi, A. Phys. Chem. Chem. Phys. 2008, 10, 5941–5952. (26) Kubaˇr, T.; Kleinekathöfer, U.; Elstner, M. J. Phys. Chem. B 2009, 113, 13107–13117. (27) Kubaˇr, T.; Elstner, M. J. Phys. Chem. B 2010, 114, 11221–11240. (28) Kubaˇr, T.; Elstner, M. Phys. Chem. Chem. Phys. 2013, 15, 5794–5813. (29) Kubaˇr, T.; Elstner, M. J. R. Soc. Interface 2013, 10, 20130415. (30) Gutiérrez, R.; Caetano, R. A.; Woiczikowski, B. P.; Kubaˇr, T.; Elstner, M.; Cuniberti, G. Phys. Rev. Lett. 2009, 102, 208102. (31) Popescu, B.; Woiczikowski, P. B.; Elstner, M.; Kleinekathöfer, U. Phys. Rev. Lett. 2012, 109, 176802. (32) Kim, H. W.; Kelly, A.; Park, J. W.; Rhee, Y. M. J. Am. Chem. Soc. 2012, 134, 11640–11651. (33) van der Vegte, C. P.; Prajapati, J. D.; Kleinekathöfer, U.; Knoester, J.; Jansen, T. L. C. J. Phys. Chem. B 2015, 119, 1302–1313. ˇ (34) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham III, T. E.; Laughton, C. A.; Orozco, M. Biophys. J. 2007, 92, 3817–3829. (35) Rosa, M.; Corni, S.; Di Felice, R. J. Chem. Theory Comput. 2014, 10, 1707–1716. (36) Rosa, M.; Corni, S.; Di Felice, R. J. Chem. Theory Comput. 2013, 9, 4552–4561. (37) Rapino, S.; Zerbetto, F. Langmuir 2005, 21, 2512–2518. (38) Modi, N.; Winterhalter, M.; Kleinekathöfer, U. Nanoscale 2012, 4, 6166–6180. (39) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701–1718.

31



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

(40) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (41) Kubaˇr, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. J. Phys. Chem. B 2008, 112, 7937–7947. (42) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260–7268. (43) Beratan, D. N.; Onuchic, J. N.; Winkler, J. R.; Gray, H. B. Science 1992, 258, 1740– 1741. (44) Kubaˇr, T.; Elstner, M. J. Phys. Chem. B 2008, 112, 8788–8788. (45) Woiczikowski, P. B.; Kubaˇr, T.; Gutiérrez, R.; Caetano, R. A.; Cuniberti, G.; Elstner, M. J. Chem. Phys. 2009, 130, 215104. (46) Mantz, Y. A.; Gervasio, F. L.; Laino, T.; Parrinello, M. J. Phys. Chem. A 2007, 111, 105–112. (47) Blumberger, J.; Lamoureux, G. Mol. Phys. 2008, 106, 1597–1611. (48) Leontyev, I. V.; Stuchebrukhov, A. A. Phys. Chem. Chem. Phys. 2011, 13, 2613–2626. (49) Croy, A.; Saalmann, U. Phys. Rev. B 2009, 80, 245311. (50) Welack, S.; Schreiber, M.; Kleinekathöfer, U. J. Chem. Phys. 2006, 124, 044712. (51) Ozaki, T. Phys. Rev. B 2007, 75, 035123. (52) Popescu, B.; Kleinekathöfer, U. Phys. Status Solidi B 2013, 250, 2288–2297. (53) Popescu, B.; Rahman, H.; Kleinekathöfer, U. J. Chem. Phys. 2015, 142, 154103. (54) Zheng, X.; Wang, F.; Yam, C. Y.; Mo, Y.; Chen, G. Phys. Rev. B 2007, 75, 195127.

32


Page 32 of 34

Page 33 of 34

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


(55) Porath, D.; Bezryadin, A.; De Vries, S.; Dekker, C. Nature 2000, 403, 635–638. (56) Voityuk, A. A.; Siriwong, K.; Rösch, N. Angew. Chem. Int. Ed. 2004, 43, 624–627. (57) Popescu, B.; Croy, A. New J. Phys. 2016, 18, 093044.

33



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

Page 34 of 34

A Multiscale Approach - ACS Publications - American Chemical Society

Recommend Documents