RT-TDDFT study of hole oscillations in B-DNA monomers and dimers

We employ Real-Time Time-Dependent Density Functional Theory to study hole oscillations within a B-DNA monomer (one base pair) or dimer (two base pairs). Placing the hole initially at any of the bases which make up a base pair, results in THz oscillations, albeit of negligible amplitude. Placing the hole initially at any of the base pairs which make up a dimer is more interesting: For dimers made of identical monomers, we predict oscillations with frequencies in the range $f \approx$ 20-80 THz, with a maximum transfer percentage close to 1. For dimers made of different monomers, $f \approx$ 80-400 THz, but with very small or small maximum transfer percentage. We compare our results with those obtained recently via our Tight-Binding approaches and find that they are in good agreement.

We employ Real-Time Time-Dependent Density Functional Theory to study hole oscillations within a B-DNA monomer (one base pair) or dimer (two base pairs).Placing the hole initially at any of the bases which make up a base pair, results in THz oscillations, albeit of negligible amplitude.Placing the hole initially at any of the base pairs which make up a dimer is more interesting: For dimers made of identical monomers, we predict oscillations with frequencies in the range f ≈ 20-80 THz, with a maximum transfer percentage close to 1.For dimers made of different monomers, f ≈ 80-400 THz, but with very small or small maximum transfer percentage.We compare our results with those obtained recently via our Tight-Binding approaches and find that they are in good agreement.

I. INTRODUCTION
Carriers can be inserted in DNA via electrodes or generated by UV irradiation or by oxidation and reduction.We had better discriminate between the words transport (implying the use of bias between electrodes), transfer (a carrier moves without bias) and migration (a transfer over rather long distances).Charge transport, transfer and migration in DNA are affected by many external factors such as chemical environment, interaction with the substrate or quality of contacts [1], as well as by many intrinsic factors like the sequence of bases [2].Charge transfer in DNA is important for DNA damage and repair [3][4][5][6] and it may be used as an indicator for the separation between cancer tumor and healthy tissue [7].Although (unbiased) charge transfer in DNA nearly vanishes after 10 to 20 nm [8][9][10], DNA is still a promising component for (biased) charge transport in molecular electronics, e.g., as a short molecular wire [11].It may also be exploited in nanotechnology for nanosensors [12] and molecular wires [13,14].
The scope of this article is to study the effect of base sequence on charge transfer, in all possible monomers (i.e.within one base pair) and dimers (i.e.within two base pairs), using Real-Time Time-Dependent Density Functional Theory (RT-TDDFT) [15,16].In previous articles, we studied small and long B-DNA segments [8-10, 17, 18] with Tight Binding (TB) approaches, specifically with a wire model and an extended ladder model.We determined the spatio-temporal evolution of an extra carrier (hole or electron) along a N base-pair DNA segment and showed that the carrier movement has frequency content mainly in the THz domain, i.e., depending on the case, in the NIR, MIR and FIR range [19].This part of the electromagnetic (EM) spectrum is significant because it can be used to extract complementary to traditional spectroscopic information e.g. about hydrogen bond stretching, bond torsion in liquids and gases, low-frequency bond vibrations and so on, and because it is relatively non-invasive compared to higher-frequency regions of the EM spectrum [20].
In this article we study B-DNA monomers and dimers with RT-TDDFT, which is one of the few computationally viable techniques to model carrier dynamics in molecular systems [21][22][23][24].RT-TDDFT can simulate density fluctuations of the order of attoseconds to femtoseconds and can model the effects of an ultrafast, intense, short-pulse laser field, including the full (nonlinear) response of the electron density [25].RT-TDDFT also provides the full frequency dependence of properties of interest via Fourier transform from the time domain to the frequency domain.In addition, the determination of many-particle eigenstates is avoided, and the calculation is reduced to an independent particle problem carried out in real time [26].Recently, RT-TDDFT was applied [27] to study the charge transfer between the oxidized or reduced a-sulfur monomers and the neighbouring neutral ones.Another recent work used RT-TDDFT to study electron transfer through oligo-p-phenylenevinylene (OPV) and carbon bridged OPV (COPV) [28].
The rest of this article is organized as follows.In Section II we sketch the basic formalism of RT-TDDFT.In Section III we give some computational details and explain our notation.In Section IV we present our results.Finally, in Section V we state our conclusions.

II. RT-TDDFT FORMALISM
Density Functional Theory (DFT) [29,30] is an efficient method for treating ground state properties of many electron systems e.g.molecules or solids.Years after its development, it was extended [31] to time dependent systems (TDDFT).Specifically, the Time-Dependent Kohn-Sham (TDKS) equations with an effective potential energy υ KS (r, t), uniquely described by the time-dependent charge density, ρ(r, t), are, in atomic units, The charge density is the sum over all occupied orbitals j = 1, 2, . . .N occ , i.e., υ ext (r, t) includes external fields and nuclear potentials, υ H (r, t) is the Hartree potential energy.Exchange and correlation effects are included in υ xc [ρ](r, t).RT-TDDFT is based on a direct numerical integration of Eq. 1.This differs from the traditional linear-response approach, which is not actually a time-resolved method but instead solves Eq. 1 in the frequency domain for the excitation energies of a system subject to a small perturbation [23].Within RT-TDDFT, we solve the TDKS equations and obtain the electron density at each time step.The electron density is then used for the calculation of the Hamiltonian in the next cycle of the self-consistent process.

III. COMPUTATIONAL DETAILS AND NOTATIONS
With the notation YX we mean that the bases Y and X of two successive base pairs Y-Y compl and X-X compl are located at the same strand in the direction 5 ′ -3 ′ .X compl (Y compl ) is the complementary base of X (Y).In other words, the notation YX means that Y-Y compl is the one base pair and X-X compl is the other base pair, separated by 3.4 Å and twisted clockwise by 36 o relatively to the first base pair, along the growth axis of the nucleotide chain.For example, the notation AC denotes that one strand contains A and C in the direction 5 ′ -3 ′ and the complementary strand contains T and G in the direction 3 ′ -5 ′ .The notation (10) means that the hole is initially placed at the 1st base pair of the dimer and the notation (01) means that it is initially placed at the 2nd base pair of the dimer.
For our DFT and RT-TDDFT calculations we used the NWChem open-source computational package [32].Additionally, we performed Fourier analysis, implemented with MATLAB, to determine the frequency content of the charge oscillations.
The geometry of monomers and all dimers was produced via the BIOVIA Discovery Studio [33].The backbone was removed and H atoms were added at the corresponding atomic sites.
The range-separated functional CAM-B3LYP [34], which is appropriate for the correct estimation of the exchange energy, both at short and long ranges, was used for all monomers and dimers.The calculations were performed using the 6-31++G** [35], 3-21++G [36] and aug-cc-pVDZ [37] basis sets, which include diffuse functions, for all systems.Within a small number of time steps, we noticed that our results were similar (cf.Appendix).Hence, we chose to perform our longer simulations with the 6-31++G** basis set [35], that is, the least computationally expensive between the converging sets.
In a Gaussian basis set, it is most natural to use the single particle reduced density matrix, whose time evolution is governed by the von Neumann equation.The Magnus propagator is used in NWChem's RT-TDDFT implementation, which is both stable and conserves the density matrix idempotency [23].At the end of each time step the fragments' charge is calculated with an appropriate population analysis method, along with the total dipole moment.
For monomers, the initial state (i.e.hole localized at a specific base) was produced by CDFT, where the charge constraint was calculated with Löwdin population analysis that is also used in the subsequent RT-TDDFT simulation.For dimers, the initial state (i.e.charge localization) was produced by the following procedure: first, a ground state DFT calculation of the charged and neutral isolated monomers was performed, yielding the corresponding monomers' eigenstates.Then, the dimer's eigenstates were approximated in a non self-consistent manner (NOSCF) by the orthogonalized monomers' eigenstates.Since the monomers comprising each dimer are 3.4 Å apart, this is a fairly good approximation that also helps circumvent CDFT convergence problems.
Löwdin [38] population analysis was integrated into RT-TDDFT module of NWChem for the calculation of each monomer's charge at each time step.Löwdin population analysis is known to be much less basis-set dependent and also it does not suffer from the ultra-fast charge oscillations that Mulliken analysis (which is the default scheme in NWChem's RT-TDDFT) artificially introduces in RT-TDDFT charge simulations.As a result, Löwdin population analysis gives a more clear picture of charge oscillations convergence towards the basis-set limit, cf.Figs.7-8.The main frequencies of charge oscillations are extracted from the results via Fourier analysis.In order to increase the resolution, zero-padding with appropriate signal attenuation was used where necessary.

IV. RESULTS
Before discussing our RT-TDDFT results, we give a brief summary of the TB picture.In TB, charge transfer between the two bases of a monomer or between the two base pairs of a dimer can be easily analytically solved [8,17,18].In both cases, it is a movement between the two sites i, j of a two-site system and it is mathematically similar to Rabi oscillations in Quantum Optics [39].The maximum transfer percentage p, i.e., the maximum probability to find the carrier at the site where it was not initially placed is where ∆ ij is the difference between the two sites' on-site energies and t ij is the hopping parameter between the two sites.Since between the two bases of a monomer, ∆ ij is much larger than t ij , p is negligible.For movements between the two monomers of a dimer things are different: For dimers made of identical monomers, since ∆ ij = 0, it follows that p = 1.For dimers made of different monomers p < 1; in fact, it is usually negligible with one exception, a hole moving between the two monomers of a GA dimer, where p is of the order of 0.25.In TB, the frequency and period of charge oscillations is given by hence, increasing t ij and / or ∆ ij results in higher frequencies.In what follows, this simple picture is qualitatively obeyed by our RT-TDDFT results.Hole oscillations between the two bases of a base pair are in the THz regime, but have negligible maximum transfer percentage, in agreement with the TB results [18].
Hole oscillations between the two base pairs of a dimer, with the 6-31++G** basis set, are shown in Fig 1 .We place the hole initially at one of the monomers which make up the dimer and call p the maximum transfer percentage to the other monomer.For dimers made of identical monomers (AT, TA, GC, CG) the hole is transferred almost completely to the other monomer (p ≈ 0.9).For dimers made of different monomers (AC, CA, GA, AG), p is close to zero with the exception of the GA dimer, where p ≈ 0.3.Our calculations for the cases AT(01), TA(01), GC(01), and CG(01) are not included in Fig. 1, since they are identical with the cases AT (10), TA (10), GC (10), CG (10), respectively, as expected from base-pair symmetry and found also in Refs.[8,17,18].
In Fig. 2 we present the mean probabilities to find the extra hole at each monomer, having placed it initially either at the 1st monomer (10) or at the 2nd monomer (01).Let us call P i the mean probability at the monomer where the hole is initially placed, and P f the mean probability at the other monomer.We observe that for dimers made of identical monomers, where p ≈ 1, we have P i ≈ P f ≈ 0.5, while for dimers made of different monomers P i ≈ 1 and P f ≈ 0, with the exception of the GA dimer.This, again, agrees qualitatively with the TB picture [8,17,18].
In Fig. 3 we depict the maximum transfer percentage p as well as the electronic coupling energy (also called the electron transfer matrix element) V RP between the reactant and product states R and P [40][41][42][43], which reflects the magnitude of the interaction between the two monomers.In our case reactant state corresponds to the hole at the initial placement and product state corresponds to the hole at the other monomer.We also show the energy difference δE between reactant and product states.In analogy with a Rabi oscillation, we expect charge transfer to increase increasing V RP and to decrease decreasing δE.Additionally, for δE ≈ 0, we expect p to be close to 1.For the AT, TA, GC and CG dimers, δE ≈ 0, hence p is indeed close to 1.For dimers made of different monomers, we observe that AC, CA and AG have large δE ≈ 0.45-0.75eV and small V RP ≈ 0-1 eV, therefore their p is insignificant.An exception is the dimer GA which not only has a large δE ≈ 0.65 eV, but also the largest V RP of all dimers ≈ 5 eV which makes charge transfer significant, with p ≈ 0.3.
In Fig. 4 we present the frequencies of hole oscillations obtained by Fourier analysis.For dimers made of identical monomers (upper panel) the frequencies f ≈ 20-80 THz and their amplitudes ≈ 0.2-0.5.For dimers made of different  (10) means that the hole is initially placed at the 1st base pair of the dimer, while the notation (01) means that the hole is initially placed at the 2nd base pair of the dimer.
monomers (lower panel) f ≈ 80-400 THz and the amplitudes ≈ 0.003-0.02,except for the two GA cases which have amplitudes ≈ 0.1-0.2.Obviously, larger amplitudes reflect larger transfer.As expected by TB, for dimers made of different monomers we generally get higher frequencies than for dimers made of identical monomers.Additionally, in Fig. 5 we observe that for dimers made of identical monomers the frequencies follow |V RP |, which does not hold for dimers made of different monomers.This, again, agrees with the TB prediction (Eq.4).Finally, in Fig. 6 we compare our RT-TDDFT results with those of TB [8,18] for two different sets of TB parameters: (a) used in Ref. [18] (calculated in Ref. [44]) and (b) used in Ref. [8].We observe that the maximum transfer percentages obtained by RT-TDDFT are in good agreement with those obtained by TB.For the periods, the results, both for RT-TDDFT and TB are quite close especially when we compare them with parametrization (b) used in Ref. [8].

FIG. 1 :
FIG.1: Net charge versus time for the monomer where the hole was initially placed: (upper panel) dimers made of identical monomers, (lower panel) dimers made of different monomers.The notation(10) means that the hole is initially placed at the 1st base pair of the dimer, while the notation (01) means that the hole is initially placed at the 2nd base pair of the dimer.

FIG. 2 :FIG. 3 :
FIG.2:The mean probabilities to find the extra hole, at each base pair of a dimer, having placed it initially at the 1st(10) or at the 2nd (01) base pair.

FIG. 7 :
FIG.7: Time evolution of hole transfer in dimers made of identical monomers, for the three different basis sets.

FIG. 8 :
FIG. 8: Time evolution of hole transfer in dimers made of different monomers, for the three different basis sets.