Structural fluctuations and quantum transport through DNA molecular wires: a combined molecular dynamics and model Hamiltonian approach

Charge transport through a short DNA oligomer (Dickerson dodecamer) in presence of structural fluctuations is investigated using a hybrid computational methodology based on a combination of quantum mechanical electronic structure calculations and classical molecular dynamics simulations with a model Hamiltonian approach. Based on a fragment orbital description, the DNA electronic structure can be coarse-grained in a very efficient way. The influence of dynamical fluctuations arising either from the solvent fluctuations or from base-pair vibrational modes can be taken into account in a straightforward way through time series of the effective DNA electronic parameters, evaluated at snapshots along the MD trajectory. We show that charge transport can be promoted through the coupling to solvent fluctuations, which gate the onsite energies along the DNA wire.


I. INTRODUCTION
The electrical response of DNA oligomers to applied voltages is a highly topical issue which has attracted the attention of scientists belongig to different research communities. The variability of experimental results is reflected in DNA being predicted to be an insulator, 1 a semiconductor, 2,3 or a metallic-like system. 4,5 This fact hints not only at the difficulties encountered in carrying out well-controlled single-molecule experiments, but also at the dramatic sensitivity of charge migration through DNA molecules to intrinsic, system-related or extrinsic, set-up mediated factors: the specific base-pair sequence, internal vibrational excitations, solvent fluctuations, and the electrode-molecule interface topology, among others.
As a result, the theoretical modelling of DNA quantum transport remains a very challenging issue that has been approached from many different sides, see e.g., Refs. 6,7,8,9 for recent reviews. While most of the models originally used started from a static picture of the DNA structure, 10,11,12,13 it has become meanwhile clearer that charge migration through DNA oligomers attached to electrodes may only be understood in the context of a dynamical approach. 14,15,16,17,18,19,20,21 Hole transfer experiments 22,23,24,25,26 had already hinted at the strong influence of DNA structural fluctuations in supporting or hindering charge propagation. Hence, it seems natural to expect that dynamical effects would also play a determining role in charge transport processes for molecules contacted by electrodes. A realistic inclusion of the influence of dynamical effects onto the transport properties can however only be achieved via hybrid methodologies combining a reliable description of the biomolecular dynamics and electronic structure with quantum transport calculations. Ab initio calculations for static biomoecular structures can provide a very valuable starting point for the parametrization of model Hamiltonians; 27,28,29,30,31,32,33,34,35 however, a full firstprinciple treatment of both dynamics and electronic structure lies outside the capabilities of state-of-the-art methodologies.
In this paper, we will ellaborate on a recent study on homogeneous DNA sequences 20 by addressing in detail some methodological issues. Our focus will be on the so called Dickerson dodecamer 36 The dynamical information provided in this way builds the starting point of our treatment of quantum transport through biomolecular wires. In the next Section, we briefly describe the computational methodology used to obtain the effective electronic parameters of the model Hamiltonian, which will be then introduced in Sec. III, where we also illustrate how to relate the Hamiltonian of Eq. (1) to a different model describing the coupling of a timeindependent electronic system to a bosonic bath. Further, expressions for the electrical current as well as the relation between the auto-correlation function of the onsite energy fluctuations and the spectral density of the bosonic bath will be derived. Finally, we discuss in Sec. IV the transport properties of the Dickerson dodecamer in vacuum and in presence of a solvent. We stress that in contrast to other models which explicitly contain the coupling to vibrational excitations or to an environment 41,42,43,44,45 at the price of introducing several free parameters, our methodology potentially contains the full dynamical complexity of the biomolecule as obtained from the MD simulations. One main advantage of our approach is the possibilty to progressively improve the degree of coarse-graining by an appropriate re-definition of the molecular fragments.

A. Computational methodology
We will first give an overview of the fragment-orbital approach used in our computations; further details can be found elsewhere. 38,39,46 The core of our method is based on a combination of a charge self-consistent densityfunctional parametrized tight-binding approach (SCC-DFTB) and the fragment orbital concept; 46 both are used to compute the electronic parameters ǫ j (t) and V j,j+1 (t) for the effective tight binding model introduced in Eq. (1) in a very efficient way, see Fig. 1 for a schematic representation. The V j,j+1 are calculated using the highest occupied molecular orbital Φ i computed for isolated bases as V j,j+1 = Φ j |H|Φ j+1 , where the Φ j 's can be expanded in a valence atomic orbital basis η µ on a given fragment: Φ j = µ c j µ η µ . The c j µ are obtained from calculations on isolated bases and stored for subsequent use to calculate V j,j+1 . Therefore, one can write: The Hamilton matrix in the atomic orbital basis H µν = η µ |H|η ν evaluated using the SCC-DFTB Hamiltonian matrix is pre-calculated and stored, thus making this step extremely efficient, i.e., it can be calculated for geometry snapshots generated by a classical molecular dynamics simulation even for several nanoseconds. Additionally, the minimal LCAO basis set used in the standard SCC-DFTB code has been optimized for the calculation of the hopping matrix elements, 38,39 and the results are in very good agreement with other approaches. 28,46 Concerning now the coupling to the solvent, a hybrid quantum mechanics/molecular mechanics (QM/MM) approach has been used, implemented in the SCC-DFTB code, 47 leading to the following Hamiltonian matrix in the valence atomic orbital basis η ν : Here, Q δ are the Mulliken charges of the quantum-mechanical region and the Q A are the charges in the MM region (backbones, counter-ions, and water), S µν is the atomic orbital Thus, we have found that the apparently most important modes are located around 1600 cm 1 corresponding to a base skeleton mode, and at 800 cm −1 related to the water modes. 39 Both contributions modulate the onsite energies significantly on a short time scale and a long time scale of about 1 ps, respectively.

B. Model Hamiltonian for electronic transport
Using Eq. (1) directly for quantum transport calculations may mask to some degree different contributions (solvent, base dynamics) to charge propagation through the DNA π-stack.
Moreover, since Eq. (1) contains random variables through the time series, we are confronted with the problem of dealing with charge transport in an stochastic Hamiltonian. This is a more complex task which has been addressed e.g., in the context of exciton transport 48,49,50,51 but also to some degree in electron transfer theories. 52,53,54 Here, we adopt a different point of view and formulate a model Hamiltonian, where the relevant electronic system, in this case the fragment orbital-derived effective DNA electronic system, is coupled explicitly to a bosonic bath. The latter will encode through its spectral density the dynamical information drawn from the MD simulations on internal base dynamics as well as solvent fluctuations.
The Hamiltonian can be written in the following way: The time averages (over the corresponding time series) of the electronic parameters ǫ j t and V j,j+1 t have been split off to provide a zero-order Hamiltonian which contains dynamical effects in a mean-field-like level. The effect of the fluctuations around these averages is hidden in the vibrational bath, which is assumed to be a collection of a large (N → ∞) number of harmonic oscillators in thermal equilibrium at temperature k B T . The bath will be characterized by a spectral density J(ω) which can also be extracted from the MD simulations, as shown below and in Ref. 55 . Since we are interested in calculating the electrical response of the system, the charge-bath model has to also include the coupling of the system to electronic reservoirs (electrodes). The coupling to the electrodes will be treated in a standard way, using a tunneling Hamiltonian H tunnel which describes the coupling to the s-lead with s= left (L) or right (R). Later on, the so called wide-band limit will be introduced (the corresponding electrode self-energies are purely imaginary and energy-independent), thus reducing the electrode-DNA coupling to a single parameter.
The previous model relies on some basic assumptions that can be substantiated by the results of the MD simulations: 39 (i) The complex DNA dynamics can be well mimic within the harmonic approximation by using a continuous vibrational spectrum; (ii) The simulations show that the local onsite energy fluctuations are much stronger in presence of a solvent than those of the electronic hopping integrals (see also the end of the previous section), so that we assume that the bath is coupled only diagonally to the charge density fluctuations; (iii) Fluctuations on different sites display rather similar statistical properties, so that the chargebath coupling λ α is taken to be independent of the site j. This latter approximation can be lifted by introducing additional site-nonlocal spectral densities J j,j+1 (ω); this however would make the theory more involved and less transparent. In (not shown) display even shorter relaxation times, so that the approximation V j,j+1 (t) = V j,j+1 (t) t is enough for our purposes (the hopping integrals are self-averaging).
In order to deal with the previous model, we first perform a polaron transformation of the Hamiltonian Eq. (4), using the generator , which is nothing else as a shift operator of the harmonic oscillators equilibrium positions. The parameter g α = λ α /Ω α gives an effective measure of the electron-vibron coupling strength. Since we will work in the wide-band limit for the electrode self-energies, the renormalization of the tunneling Hamiltonian by a vibronic operator will be neglected. As a result, we obtain a Hamiltonian with decoupled electronic and vibronic parts and where the onsite energies are shifted as ǫ j t → ǫ j t − ∞ 0 dωJ(ω)/ω. However, as it is well-known, 56 the retarded Green function of the system is now an entangled electronic-vibronic object that can not be treated exactly; we thus decouple it in the approximate way: 41,57 In this equation, θ(t − t ′ ) is the Heaviside function and the pure bosonic operator In the last row of Eq. (5) we can pass to the continuum limit and express φ(t) in terms of the bath spectral density J(ω): 58 C. The electrical current We derive in this section the expression we are going to use to calculate the electrical current through the DNA oligomer under study. Starting point is the well-known Meir-Wingreen expression for the current from lead s: 59 Now, we can exploit the decoupling approximation used in Eq. (5) together with the wideband limit in the electrode-molecule coupling to write e.g., for the left electrode: Hereby we have used the explicit expressions for the greater-and lesser-Green functions: the index 0 indicating that the vibrational degrees of freedom have already been decoupled. The transmission-like function t(E) is given by t(E) = Tr G 0 (E)Γ L G † 0 (E)Γ R . The retarded matrix Green function G 0 (E) is calculated without electron-bath coupling but including the interaction with the electrodes: G −1 0 (E) = E + i η − H 0 + i Γ L + i Γ R , H 0 being the electronic part of the Hamiltonian of Eq. (4). Using these results, the right-going current can be written as can be performed, yielding: In the former expression √ κ therm is related to an inverse decoherence time. In the high temperature limit κ therm ∼ k B T E reorg . The Fourier transform of the previous expression can be calculated straightforward and gives: Thus, the current calculated using Eq. (8) becomes a convolution of t(E) with a Gaussian function. Notice the similarity of this expression with a Frank-Condon factor appearing in the Marcus electron transfer theory.

D. Spectral density from molecular dynamics
The next issue at stake is how to relate the bath spectral density J(ω) to the time series of the electronic parameters as obtained through the molecular dynamic simulations. To illustrate this point we will consider the simple case of a single level, whose site energy is a Gaussian random variable, coupled to left and right electrodes according to the Hamiltonian (we assume for simplicity that ǫ(t) = 0): To deal with this problem we may use equation of motion techniques for the retarded dot . A similar equation of motion for the latter function allows to introduce the electrode self-energies: Note that the Green function thus obtained is still a random function, since no average over the random variable δǫ(t) has been yet performed. Using the wide-band limit in the lead Introducing the time-exponential U(t, t ′ ) = exp(−(i / ) t t ′ ds(δǫ(t) − i Γ) allows to write a closed solution for the dot's Green function: The average Green function over the random variable δǫ(t) yields now in energy-space: Performing now a cumulant expansion 60 of the averaged exponential and taking into account that cumulants higher than the second one exactly vanish due to the Gaussian nature of the fluctuations and to the fact that δǫ(t) is a classical variable, we get (with t ′ = 0): which is the formally exact solution of the problem. We may now look at the same problem from a different point of view by considering explicitly the coupling of a single site with time-independent onsite energy to a continuum of vibrational excitations. Using the polaron transformation together with the approximations introduced at the beginning of Sec. II B, we can write the retarded Green function as: where φ(t) has been already defined in Eq. (6). By comparison of Eqs. ( 15) and (16), there must exist a relation between the (real) correlation function and the real part of φ(t). The latter can be re-written in the following way: so that from this we conclude that From here it follows by inversion: The previous result was obtained in Ref. 55 in the weak coupling, perturbative regime to the vibrational system; we have shown here that it can be also extended to the case of arbitrary coupling. Similar expressions could be derived for (spatially) non-local spectral densities.

E. Results
To illustrate our methodology we have focused on the Dickerson dodecamer (DD) which has a non-homogeneous base sequence. In contrast to our previous study 20  (1/T MD ) j T (E, t j ). This quantity is expected to provide some qualitative insight into different factors affecting transport (solvent vs. vacuum) but its use is restricted by the fact that for longer chains inelastic, fluctuation mediated channels will play an increasing role and then a pure elastic transmission can not catch all the transport physics. In this latter case, the model Eq. (5) seems to us more appropriate since it includes the dressing of the electronic degrees of freedom by the structural fluctuations. In Fig. 3 we show the averaged transmission function (upper panel) for the cases where solvent effects are included or neglected, respectively. We clearly see that solvent-mediated fluctuations considerably increase and broaden the transmission spectrum. Its fragmented structure is simply due to the fluctuations in the onsite energies, which make the system effectively highly disordered.
Another way of looking at the effect of fluctuations is via the introduction of a coherence parameter 20 which is defined as C(E) = T (E, t) 2 t / T (E, t) 2 t . If the fluctuations are weak C(E goes to one, while a strong fluctuating system will lead to a considerably reduction of C(E). Obviously, this parameter has only a clear meaning over the spectral support of the transmission function. C(E) is shown in the lower panel of Fig. 3; the coherence parameter in vacuum conditions can become larger than in solvent for some small energy regions due to the reduction of the fluctuations. However, this does not turn to be enough to increase the transport efficiency for the special base sequence of the Dickerson DNA, since this system has in the static limit a distribution of energetic barriers due to the base pair sequence. To further illustrate the influence of the solvent, we have plotted in Fig. 4 the time averaged onsite energies along the model tight-binding chain. Remarkably, the presence of the solvent "smoothes" the averaged energy profile (though the amplitude of the fluctuations clearly becomes stronger). To compute the current through the system, we use

III. CONCLUSIONS
We have presented a hybrid methodology which allows for a very efficient coarse-graining of the electronic structure of complex biomolecular systems as well as to include the influence of structural fluctuations into the electronic parameters. The time series obtained in this way allow to map the problem onto an effective low dimensional model Hamiltonian describing the interaction of charges with a bosonic bath, which comprises the dynamical fluctuations.
The possibility to parametrize the bath spectral density J(ω) using the information obtained from the time series makes our approach very efficient, since we do not need to use the typical phenomenological Ansätze (ohmic, Lorentzian, etc) to describe the bath dynamics.
The example presented here, the Dickerson dodecamer, shows in a very clear way that solvent-mediated gating may be a very efficient mechanism in supporting charge transport if the static DNA reference structure already posseses a disordered energy profile (due to the base sequence). Our method is obviously not limited to the treatment of DNA but it can equally well be applied to deal with charge migration in other complex systems like molecular organic crystals or polymers, where charge dynamics and coupling to fluctuating environments plays an important role. 52,53,54,64,65,66,67,68 IV. ACKNOWLEDGMENTS