---
abstract: 'The electronic structure of an infinite poly-guanine - poly-cytosine DNA molecule in its dry A-helix structure is studied by means of density-functional calculations. An extensive study of 30 nucleic base pairs is performed to validate the method. The electronic energy bands of DNA close to the Fermi level are then analyzed in order to clarify the electron transport properties in this particularly simple DNA realization, probably the best suited candidate for conduction. The energy scale found for the relevant band widths, as compared with the energy fluctuations of vibrational or genetic-sequence origin, makes highly implausible the coherent transport of electrons in this system. The possibility of diffusive transport with sub-nanometer mean free paths is, however, still open. Information for model Hamiltonians for conduction is provided.'
author:
- Emilio Artacho
- Maider Machado
- 'Daniel Sánchez-Portal'
- Pablo Ordejón
- 'José M. Soler'
date: 19 September 2002
title: Electrons in Dry DNA from Density Functional Calculations
---

Introduction
============

The idea of using DNA chains as molecular wires in nano-devices has been proposed and is being explored in the nano-science community [@Mirkin96; @Alivisatos96; @Braun98], given the extremely high selectivity of a DNA strand to pair with its complementary partner, which facilitates the design of complex circuitry at the nanometer scale. Even though some researchers have chosen to coat [@Braun98] or modify [@Rakitin01] the DNA wires in order to improve its conduction properties, the original idea was to use pure DNA, assuming that electric current could flow along the wires. This assumption has proven to be less trivial than expected, however, and has originated an important controversy.

Direct measurements show from low resistance DC conduction [@Fink99] (even superconduction [@Kasumov01]), to no conduction at all [@Braun98; @dePablo00; @Gomez02], including results with non-ohmic characteristics [@Porath00; @Yoo01] or appreciable ac conductance [@Gruner00]. The transport has been proposed to be band-structure-like [@Porath00], by tunneling [@Beratan97], solitonic [@BenJacob98], polaronic [@Schuster99; @Conwell00; @Yoo01], or facilitated by fluctuations [@Ladik99; @Bruinsma00]. This rich scenario is not surprising when considering the comparably small amount of clean-cut information available (mostly theoretical) about basic aspects of the electronic structure of dry DNA [@Ladik71; @Ladik83; @Pablo97; @Yang98; @Ye00], (in all the experiments mentioned above the samples were dried in one way or another). In addition to the fundamental complexity of the system (large system size, weak interactions, structure floppiness), there is a substantial lack of detailed knowledge about the environment the DNA molecules are left with after the drying process. Most likely, residual water molecules will stay adhered to the chain, together with the cations counterbalancing the negative charge of the molecule, all quite uncontrolled.

A theoretical simulation benefits from absolute control at the atomic scale. In this particular case, the system can be prepared in the best possible way for conduction. Previous simulations [@Ladik71; @Ladik83; @Pablo97; @Yang98; @Ye00] have offered very insightful information for related problems, in spite of being based on relatively drastic approximations in the electronic structure, due to the system size. They were mostly aimed at the description of DNA properties in live conditions, as were the recent first-principles studies in wet DNA [@Landman01; @Zhang02; @Parrinello02].

Recent developments in linear-scaling first-principles [@Soler02; @Bowler02; @Scuseria00] methods based on density-functional theory [@Payne92] (DFT) allow the study of an infinite DNA double helix from first principles. DFT functionals based on the Generalized Gradient Approximation (GGA) have shown to provide very satisfactory accuracy in molecules akin to DNA [@Sponer96; @Hutter96]. In this work, an extensive study is presented on a set of 30 nucleic base pairs that validates the theoretical method. Then, an infinite DNA chain of trivial sequence is studied: all guanines (G) on one strand, all cytosines (C) on the other, in the Watson-Crick pairing arrangement. We simulate it completely dry and in vacuo, and in an A-helix structure, the most relevant one for these conditions. The electronic structure is analyzed after relaxing the structure with the first-principles forces. Genetic sequence change effects are explored on a similar system where one guanine-cytosine pair is swapped every eleven base pairs, then relaxed and analyzed. Preliminary results of this study were published in a joint experiment-theory paper [@dePablo00].

Method
======

The first-principles DFT calculations reported here were performed using the [Siesta]{} method [@Soler02; @Pablo96; @Daniel97], which has already been applied to a large variety of systems [@websiesta], including biomolecules [@dePablo00; @Marivi00]. It is a fully linear-scaling method, both in the building of the DFT Hamiltonian and in its solution. GGA  [@PBE] was used for electron exchange and correlation. Core electrons were replaced by norm-conserving pseudopotentials [@TroullierMartins91] in their fully non-local formulation [@KB]. A uniform mesh with a planewave cutoff of 150 Ry was used for integrations in real space [@Soler02]. The basis sets used are made of numerical atomic orbitals of finite range and are described below.

The unit cell of the A-helix studied here includes eleven poly(G)-poly(C) (pGpC) base-pairs with a total of 715 atoms. Periodic boundary conditions are used to simulate infinitely long DNA chains well separated from each other, with negligible mutual interactions. Only the $\Gamma$ $k$-point of the 715-atom cell was included in the calculations, which corresponds to 11 $k$-points, if the helical symmetry is considered. The initial geometry was that of an experimental low-resolution A-helix structure for the pGpC sequence [@Landoldt-Bornstein], The geometry was relaxed by means of ab initio linear-scaling DFT, as described elsewhere [@ADNA-PRB]. The final structure had still well defined A-helix characteristics [@ADNA-PRB; @Calladine97]. A second DNA realization (swapped hereafter) was built by swapping one single G-C pair into C-G out of the eleven in the cell, followed by another full relaxation.

The electronic structure is analyzed in terms of the Kohn-Sham eigenvalues and eigenfunctions, based on Janak’s theorem [@Janak78]. They are obtained by a single cube-scaling diagonalization, using the final geometry and electron density obtained with the linear-scaling method. The known band-gap problem of DFT is not important here, since the particular value of the band gap is of no relevance for the conclusions of this study as long as there is a gap, which is the case. The band widths have been shown to be well described by this method, at least up to the kind of accuracy needed in this work. A much more costly calculation including explicit electronic correlation would not be justified. Mott-like correlations, that could seem to be relevant for the narrow band-widths characteristic of the problem (see below), are not important given the small density of carriers (electrons or holes, see discussion below) expected in the system. The effect of disorder and vibrations is explicitly considered below. The next section is devoted to a study of nucleic base pairs where the focus is mainly on the hydrogen bonds, which might in principle be worst described by our methodology. It will be followed by the results for A-DNA.

Study on base pairs
===================

A systematic study of 30 nucleic base pairs has been performed with the methodology presented above, in order to assess the reliability of the approximations, but also obtaining new results. For this study we have used a standard doble-$\zeta$ polarized (DZP) basis set, namely, a double (split valence) basis for each valence orbital plus polarization functions in all the atoms. The cutoff radii for the atomic orbitals of each element were obtained for an energy shift [@Emilio99] of 50 meV.

A significant range of configurations of the four bases guanine, cytosine, adenine and thymine (G, C, A, T) are considered, the same as those studied by $\rm{\check{S}poner}$ [*et al.*]{} [@Sponer96] in their MP2 study. The Watson-Crick configurations are designated WC, and the Hoogsteen, reversed Hoogsteen and reversed Watson-Crick appear as H, RH, RWC respectively. Other configurations are distinguished simply with numbers, eg. AA1, AA2, etc. following the nomenclature of Hobza and Sandorfy [@Hobza]. The structures of the bases and base-pairs studied in this work can be found in Figures 1 and 2 of Ref. . For the numbering of the atoms we followed Ref. . In the following our results are compared with those of other methods. The direct comparison with experimental geometries has been shown to be misleading because of the important charge rearrangements that happen in the crystal phases used in experiments [@Celia00].

The interaction energy $E_{int}$ in the following is defined as the energy of the base-pair minus the energy of each base with the same geometry it has in the pair. The total stabilization energy $E_t$ is defined as the difference between the energy of the pair and that of each base in its isolated optimal geometry. The difference between both is thus the deformation energy, [*i.e.*]{}, the increase in intramolecular energy due to the geometry change when the base-pair is formed.

In this study all the energies have been corrected for basis set superposition error (BSSE) using the standard Boys-Bernardi counterpoise correction [@Boys]. The correction found for $E_{int}$ is used for $E_t$ as well, given the fact that the relaxation of the isolated base molecules has been found to change the BSSE correction by less than 10%. No BSSE correction was included in the forces. The effect of this approximation was gauged for the AA1 base-pair, for which BSSE corrected and uncorrected forces were calculated for independent relaxations, showing a deviation of 0.02 Å in the final acceptor-hydrogen distances.

  ---------- ------- -------- ---------------
     Pair        MP2   SIESTA   Deviation (%)
     GCWC      -25.8    -26.8            -3.9
     GG1       -24.7    -25.1            -1.6
    GCNEW      -22.2    -21.7             2.2
      CC       -18.8    -17.5             6.9
     GG3       -17.8    -16.6             6.7
     GA1       -15.2    -15.5            -2.0
     GT1       -15.1    -15.0             0.7
     GT2       -14.7    -14.5             1.4
     AC1       -14.3    -14.0             2.1
     GC1       -14.3    -14.7            -2.8
     AC2       -14.1    -14.7            -4.2
     GA3       -13.8    -13.8             0.0
     TAH       -13.3    -13.7            -3.0
     TARH      -13.2    -13.6            -3.0
     TAWC      -12.4    -12.3             0.8
    TARWC      -12.4    -12.3             0.8
     AA1       -11.5    -11.7            -1.7
     GA4       -11.4    -11.7            -2.6
     TC2       -11.6    -10.8             7.5
     TC1       -11.4    -10.6             7.0
     AA2       -11.0    -11.4            -3.6
     TT2       -10.6     -9.9             6.6
     TT1       -10.6    -10.1             4.7
     TT3       -10.6    -10.2             3.8
     GA2       -10.3    -10.6            -2.9
     GG4       -10.0     -7.4            26.0
     AA3        -9.8     -9.8             0.0
   2aminoAT    -15.1    -15.2            -0.7
  ---------- ------- -------- ---------------

  : Base-pair Interaction Energies ($E_{int}$, in kcal/mol) at HF/6-31G\*\* geometries.[]{data-label="eint-hf"}

Our results for $E_{int}$ as compared with the MP2 results of $\rm{\check{S}poner}$ and coworkers[@Sponer96] are presented in table \[eint-hf\], as calculated for the same geometries obtained at the Hartree-Fock (HF) 6-31G\*\* level [@Sponer96]. For all the base-pairs except GG4, the agreement is quite satisfactory, with differences smaller than 8% and much less in most cases. GG4 is an exception to the general trend with a relative difference of 26% that is reduced to 16% by extending the basis range to the longer radii given by an energy shift of 10 meV. A similar correction for the other base pairs gives much smaller variations. The standard deviation of our results compared to the MP2 values is of 0.73 kcal/mol, the same value found by $\rm{\check{S}poner}$ [*et al.*]{}[@Sponer96]using the B3LYP functional [@b3lyp] on the same geometry. The standard deviation between both sets of DFT results (our PBE and their B3LYP) is of 0.85 kcal/mol.

Concerning geometries, the only pair relaxed with MP2 is cytosine-cytosine [@Sponer96]. Table \[cytocyto\] compares the length of the hydrogen bond and both $E_{int}$ and $E_t$ for various methods, showing again a very satisfactory agreement.

[lccccc]{} & HF & MP2/HF & MP2 & DFT & This work\
\
d(N4(H)-N3) & 3.050 & 3.050 & 2.980 & 2.900 & 2.872\
\
$E_{int}$ & -17.3 & -18.8 & -20.5 & -20.4 & -21.1\
$E_{t}$ & – & -17.5 & -18.7 & -18.1 & -18.5\

The results shown in this section validate the approximations used in our method for the system under study, including the sensitive hydrogen bonds. Even if the piling up of bases was not investigated, other studies show [@Hutter96] that GGA functionals can describe them with the required accuracy when connected by the sugar-phosphate chains.

             B3LYP   BP86   VWN-BP    DZP    DZ(P)
  --------- ------- ------ -------- ------- -------
  N2 – O2    2.93    2.87   2.930    2.872   2.892
  N2 – H       -      -     1.035    1.035   1.023
  N1 – N3    2.92    2.88   2.923    2.913   2.822
  N1 – H       -      -     1.051    1.056   1.043
  O6 – N4    2.78    2.73   2.785    2.770   2.715
  H – N4       -      -     1.055    1.057   1.054

  : Hydrogen-bond distances (in Å) for the guanine-cytosine pair. DZP and DZ(P) stand for two different basis sets used in this work.

\[gc\]

The whole set of base pairs as well as the isolated bases have been relaxed with our method. The results for geometry and interaction energies can be found in [@Maider00]. We extract in table \[gc\] the results for the guanine-cytosine pair, of relevance for the DNA calculations in the next section. The results of several authors [@Sponer96; @Celia00; @Santamaria] are compared, which include different GGA functionals, namely B3LYP [@b3lyp], BP86 [@p86; @b86], and VWN-BP [@p86; @b88; @vwn]. The results of the DZP basis of this study are also compared with those of the basis set used for DNA below, DZ(P), a double-$\zeta$ basis [@Emilio99; @radii], which was polarized for phosphorous and for the atoms involved in the hydrogen bonds. The Kohn-Sham eigenvalues around the Fermi level for an isolated GC pair were checked for DZ(P), which compare well with the corresponding values for the DZP basis. The convergence of the spatial range of the basis orbitals, important for band widths in the DNA molecule, was also tested [@radii].

Results for A-DNA
=================

Figure \[bands\] shows the bands of pGpC close to the Fermi level. As expected for a well saturated molecule, the Fermi level is in a clear band gap, which is of 2.0 eV in this case (this number is to be taken with caution since usual DFT functionals tend to underestimate the band gaps of insulators by a substantial amount). Electronic transport will thus be mediated by carriers in either the top of the valence band (holes, or radical cations in biochemical language) or in the bottom of the conduction band (electrons or radical anions), introduced by doping.

In the absence of disorder, the electronic structure close to the Fermi level shows well defined bands with eleven states per unit cell, one per base pair in the unit cell. The top-most valence band has a very small bandwidth of 40 meV (see figure \[bands\]), the energy separation with the next band below being ten times larger. This first band is associated to the $\pi$-like highest occupied molecular orbital (HOMO) of the guanines. The continuous line in the figure is the tight-binding band of a one-dimensional system with one orbital per unit cell. The excellent agreement of this extremely simple model with the ab initio points \[see figure \[bands\] (b)\] demonstrates that each base pair contributes to this band with one single orbital that interacts negligibly with other orbitals in the pair. This allows the modelling of the basic physics of transport in this molecule with one single orbital per base pair with an electron hopping interaction of 10 meV. Other terms for a model Hamiltonian are quantitatively justified below.

![(Color stereogram) Surfaces of constant density (1.5 $10^5$ $e^-$/Å$^3$) for the states corresponding to the lowest unoccupied band (red), and highest occupied band (blue) of the ordered pGpC structure. The graphs were produced with Molekel [@molekel]. []{data-label="homo-lumo"}](Fig2.eps){width="\columnwidth"}

Figure \[homo-lumo\] displays the charge density of the states associated to this band, which appears almost exclusively on the guanines, with weight neither in the backbones nor in the cytosines. The lowest conduction band has a width of 270 meV, similar to the width of the spectroscopic feature observed by Porath [*et al.*]{} [@Porath00]. It is separated from the next band by 0.7 eV. Similarly to the HOMO situation, this band is made of the lowest unoccupied molecular orbital (LUMO) of the cytosines, as shown in figure \[homo-lumo\], although the overlap in this case is more pronounced, as corresponds to an eight times larger bandwidth.

Note that any matrix element between HOMO and LUMO states is very much depressed by the very small spatial overlap between them. Optical absorption will thus be very weak for photon energies corresponding to the HOMO-LUMO gap (adequately corrected including correlation and excitonic effects). The unoccupied band above the LUMO band remains of a marked cytosine character. The third unoccupied band 3.9 eV above the HOMO band (see figure \[bands\]) is the first one with substantial weight on the guanines, and thus with an appreciable oscillator strength for absorption. The comparison of these results with those obtained for wet conditions [@Parrinello02] shows that the LUMO is quite sensitive to the environment, moving from the cytosines to the cations when in presence of Na$^+$ and water.

The effect of the swap on the electronic structure of the chain is dramatic, as can be seen in figure \[bars\]. The HOMO of the swapped guanine (indicated by an asterisk in figure \[bars\]) sinks 0.55 eV (14 times the HOMO bandwith) into lower valence band levels. This stabilization is of electrostatic origin  [@ADNA-PRB]. Figure \[homos\] (left) shows the electronic density associated to the eleven highest occupied states for the swapped structure, showing the cut in the HOMO-state channel produced by the swapped pair. The situation is similar for the unoccupied band, albeit less dramatic, with the LUMO of the swapped cytosine raising 0.5 eV above the LUMO band.

All the results presented so far refer to fixed, relaxed geometries. It has been proposed, however, that the electron transfer could be mediated by phonons [@Schuster99; @Conwell00; @Ladik99; @Bruinsma00]. Indeed, there are many modes in DNA that are very soft and related to the nucleic bases, since these flat, rigid molecules hang from the backbone by a single bond. It is beyond the scope of this study to calculate phonon modes and their interactions with electrons, but interesting qualitative information on the matter can be inferred from the conjugate-gradient (CG) minimizations. If close to the minimum, any snapshot during the minimization gives the atomic positions that correspond to the freezing of vibration modes with various amplitudes. The closer to the minimum, the more important the low-frequency modes versus the high-frequency ones [@valley]. We study thus the electronic structure of geometries chosen from the CG path with an energy of around 0.9 eV above the minimum, that would correspond to room temperature or lower [@temperature]. The highest occupied energy levels of such a geometry (called ragged) is presented in figure \[bars\]. If not as dramatic as the swapping, the effect of these “vibrations" on the electrons is also important, since the guanine HOMO one-particle states spread over an energy width eight times the original band width.

![Isosurface as in figure \[homo-lumo\] for the top eleven occupied states of the swapped structure (left). The same for one particular state (the 11th) for a density value 50 times smaller (right).[]{data-label="homos"}](Fig4.eps){width="\columnwidth"}

Discussion and Conclusions
==========================

Electronic transport in DNA requires the presence of carriers. In controlled biochemical experiments (in solution), doping is induced by intercalating or attaching to the DNA specific molecules of known tendency to donate or accept electrons [@Barton99]. In the dry experiments introduced earlier there is uncontrolled doping of DNA. Since the negative charge of each phosphate group along both chains has to be compensated by a cation from solution, in the drying process cations precipitate along the DNA helix decorating it. In addition to a very probable disorder in the cation positions and variability in the amount of water molecules hydrating them, there will also be deficiency or excess of counter-cations. Missing cations introduce holes in the valence band and cations in excess (or cations with higher charge) introduce electrons in the conduction band. The drying procedures of the dry experiments do not allow estimating the concentration of carriers, not even their sign, even if (in view of the ionization potential of G and the electron affinity of C) the presence of holes has been argued to be more likely. It is reasonable to expect in any case that carrier concentrations will remain low, and that for the more drastic drying procedures the likelihood is higher for the presence of holes, since cations will be pulled away with the water. pGpC has been reported to be $p$-type in at least one experimental setup [@Yoo01].

It has been shown above that the HOMO band width is one order of magnitude smaller than the energy scale associated to sequence disorder. In terms of a one-dimensional Anderson model, this leads to electronic localization over very few base pairs. The localization is observed when plotting the separated eigenstates. Figure \[homos\] (right) shows one electronic state within the HOMO energy range. The isosurface plotted is for a density 50 times smaller that that on the left, showing clearly the localization of the state within a few base pairs. Sequence disorder will therefore dramatically depress the possible band-structure-like DC conduction in DNA of arbitrary sequence ($\lambda$-DNA). The possibility [@Ratner00] of electron hopping or tunneling, like in proteins [@Beratan97] is open, but it is not expected to show the conductances close to one quantum reported for $\lambda$-DNA by some groups [@Fink99; @Kasumov01]. For electron doping, disorder fluctuations are still comparable or larger than the LUMO band width, and the localization of the electronic states is still considerable, albeit not as extreme as for the hole-doped case.

It is known [@Phillips91] that certain types of disorder do not localize electrons as in the Anderson model. They represent, however, very particular realizations that cannot be expected from $\lambda$-DNA in general. Very interesting recent results [@Carpena02] show possible electron delocalization in 1D systems with more general long-range correlated disorder, as the disorder of natural DNA. More research is needed to see how robust is this result when including other effects like off-diagonal disorder or vibrations. In any case, the very large disorder-to-bandwidth ratio presented above makes extended states for holes very unlikely, the typical state being localized as displayed in figure \[homos\] (right).

The particular experiment of Porath et al. [@Porath00] used a controlled DNA sequence, namely pGpC, for which the sequence disorder arguments do not apply. They explain their results in terms of band-like conduction. The results and analysis of Berlin and coworkers [@Ratner00] propose band-like conduction in this particular system as well. It is tempting to support this band-like condution assuming it proceeds via electrons in the LUMO, given the agreement of the widths of our LUMO band and the feature they observe spectroscopically (the sensitivity of the LUMO to the environment mentioned above has to be kept in mind, however). Such an explanation would be by no means conclusive in the light of our “vibration" results, since a deterioration of conduction should be expected with temperature, which does not seem to happen. The spread of HOMO states with frozen-in vibrations indicates the energy scale of electron-phonon interactions, again very substantial as compared with the band widths. Even if a more rigorous study is needed to quantify these interactions and their effect, our rough estimations above allow us to expect a substantial reduction of conduction with temperature even below room temperature.

Considering the floppy modes in DNA and their effect on the electron states, and the fact that a hole on a guanine represents a positive charge amidst negatively charged guanines, the idea of self trapped (small) polarons seems not at all unlikely. Polarons have already been suggested [@Schuster99; @Conwell00] and supported [@Yoo01] in the literature. Some of the models used have assumed the polarons to be of an origin similar to that found in polyacetylene, i.e. of off-diagonal (hopping) origin. Polarons of electrostatic origin are described by an electronic Hamiltonian in which the diagonal matrix elements are the ones that change with inter-base separation and stabilize the polaron, as in Holstein’s model [@Holstein59]. A purely classical electrostatic model based on the ab initio results for the geometry and ground-state electrostatic characteristics of the system [@ADNA-PRB] displays polaronic charge excitations with a spread of a few base pairs. DFT calculations on the matter will be presented elsewhere [@Simone].

[*Acknowledgments.*]{} We are grateful to J. $\rm{\check{S}poner}$ and R. Santamar[í]{}a for making their coordinates of bases and base-pairs available to us. We are also indebted to R. Weht for many useful discussions and his help during the first stages of this work. This work has been supported by the Spanish Ministerio de Ciencia y Tecnología under grant BFM2000-1312, and by the Fundación Ramón Areces. DSP acknowledges support from Spain’s MCyT under the Ramón y Cajal program. Calculations were performed at the CCCFC of the Universidad Autónoma de Madrid and at the PSMN of the École Normale de Lyon (SESC).

C. A. Mirkin, R. L. Letsinger, R. C. Mucic, and J. J. Storhorff, Nature [**382**]{}, 607 (1996).

A. B. Alivisatos, K. P. Johnsson, X. Peng, T. E. Wilson, C. J. Loweth, M. P. Bruchez and P. J. Schultz. Nature [**382**]{}, 609 (1996).

E. Braun, Y. Eichen, U. Sivan, and G. Ben-Yoseph, Nature [**391**]{}, 775 (1998).

A. Rakitin, P. Aich, C. Papadopoulos, Y. Kobzar, A. S. Vedeneev, J. S. Lee, and J. M. Xu, Phys. Rev. Lett. [**86**]{}, 3670-3673 (2001).

H. W. Fink and C. Schönenberger, Nature [**398**]{}, 407 (1999).

A. Y. Kasumov, M. Kociak, S. Gueron, B. Reulet, V. T. Volkov, D. V. Klinov, and H. Bouchiat, Science [**291**]{}, 280-282 (2001).

P. J. de Pablo, F. Moreno-Herrero, J. Colchero, J. Gómez-Herrero, P. Herrero, A. M. Baró, P. Ordejón, J. M. Soler, and E. Artacho, Phys. Rev. Lett. [**85**]{}, 4992-4995 (2000).

C. Gómez-Navarro, F. Moreno-Herrero, P. J. de Pablo, J. Colchero, J. Gómez-Herrero, and A. M. Baró, Proc. Natl. Acad. Sci. [**99**]{}, 8484-8487 (2002).

D. Porath, A. Bezryadin, S. De Vries, and C. Dekker, Nature [**403**]{}, 635 (2000).

K.-H. Yoo, D. H. Ha, J.-O. Lee, J. W. Park, J. Kim, J. J. Kim, H.-Y. Lee, T. Kawai, and H. Y. Choi, Phys. Rev. Lett. [**87**]{}, 198102 (2001).

P. Tran, B. Alavi, and G. Gruner, Phys. Rev. Lett. [**85**]{}, 1564 (2000).

D. N. Beratan, S. Priyadarshy, and S. M. Risser, Chemistry and Biology [**4**]{}, 3 (1997) and references therein.

Z. Hermon, S. Caspi, and E. Ben-Jacob, Europhys. Lett. [**43**]{}, 482 (1998).

P. T. Henderson, D. Jones, G. Hampikian, Y. Kan, and G. B. Schuster, Proc. Natl. Acad. Sci. USA [**96**]{}, 8353 (1999).

E. M. Conwell and S. V. Rakhmanova, Proc. Natl. Acad. Sci. USA [**97**]{}, 4556 (2000).

Y. J. Ye, R. S. Chen, A. Martinez, P. Otto, and J. Ladik, Solid St. Commun. [**112**]{}, 139-144 (1999).

R. Bruinsma, G. Grüner, M. R. D’Orsogna, and J. Rudnick, Phys. Rev. Lett. 85, 4393 (2000).

J. Ladik, Int. J. Quant. Chem. [**4**]{}, 307 (1971).

P. Otto, E. Clementi, and J. Ladik, J. Chem. Phys. [**78**]{}, 4547 (1983).

J. P. Lewis, P. Ordejón, and O. F. Sankey, Phys. Rev. B [**55**]{}, 6880 (1997).

D. M. York, T. S. Lee, and W. Yang, Phys. Rev. Lett. [**80**]{} , 5011 (1998).

Y. J. Ye and Y. Jiang, Int. J. Quantum Chem. [**78**]{}, 112-130 (2000).

R. N. Barnett, C. L. Cleveland, A. Joy, U. Landman, and G. B. Schuster, Science [**294**]{}, 567-571 (2001).

H. Y. Zhang, X. Q. Li, P. Han, X. Y. Yu, and Y. J. Yan, J. Chem. Phys. [**117**]{}, 4578-4584 (2002).

F. L. Gervasio, P. Carloni, and M. Parrinello, Phys. Rev. Lett. [**89**]{}, 108102-1 - 4 (2002).

J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Cond. Matter [**14**]{}, 2745-2779 (2002).

D. R. Bowler, T. Miyazaki, and M. J. Gillan, J. Phys.: Cond. Matter [**14**]{}, 2781-2798 (2002).

K. N. Kudin and G. E. Scuseria, Phys. Rev. B [**61**]{}, 16440-16453 (2000).

M. Payne, M. Teter, D. Allan, T. Arias, and J. D. Joannopoulos, Rev. Mod. Phys. [**64**]{}, 1045 (1992).

J. Sponer, J. Leszczynski, P. Hobza, J. Phys. Chem. [**100**]{}, 1965 (1996).

J. Hutter and M. Parrinello, J. Am. Chem. Soc. [**118**]{}, 8710 (1996).

P. Ordejón, E. Artacho, and J. M. Soler, Phys. Rev. B [**53**]{}, R10441 (1996).

D. Sánchez-Portal, P. Ordejón, E. Artacho, and J. M. Soler, Int. J. Quant. Chem. [**65**]{}, 453 (1997).

See refs. in P. Ordejón, Phys. Stat. Solidi (b) [**217**]{}, 335 (2000), and for more recent applications see references in ref.  and in http://ww.uam.es/siesta.

M. V. Fernández-Serra, J. Junquera, C. Jelsch, C. Lecomte, and E. Artacho, Solid St. Commun. [**116**]{}, 395 (2000).

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. [**77**]{}, 3865 (1996).

N. Troullier and J.L. Martins, Phys. Rev. B [**43**]{}, 1993 (1991).

L. Kleinman and D. M. Bylander, Phys. Rev. Lett. [**48**]{}, 1425 (1982).

R. Chandrasekaran and S. Arnott, in [*Biophysics. Nucleic Acids. Crystallographic and Structural Data II*]{}, edited by W. Saenger, Landölt-Börnstein, New Series, Group VII, Vol. I, Pt. b (Springer-Verlag, New York 1989).

E. Artacho, P. Ordejón, D. Sánchez-Portal, and J. M. Soler (to be published).

M. A. El Hassan and C. R. Calladine, Phil. Trans. R. Soc. Lond. A [**355**]{}, 43 (1997).

J. F. Janak, Phys. Rev. B [**18**]{}, 7165-7168 (1978).

E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, Phys. Stat. Sol. (b), [**215**]{}, 809 (1999).

P. Hobza and C. Sandorfy, J. Am. Chem. Soc. [**109**]{}, 1302 (1987).

J. $\rm{\check{S}poner}$ and P. Hobza, J. Phys. Chem. [**98**]{}, 3161 (1994).

C. F. Guerra, F. M. Bickelhaupt, J. G. Snijders, and E. J. Baerends, J. Am. Chem. Soc. [**122**]{}, 4117-4128 (2000).

S. F. Boys and F. Bernardi, F. Mol. Phys. [**19**]{}, 553 (1970).

A. D. Becke, J. Chem. Phys. [**98**]{}, 5648 (1994).

M. Machado, P. Ordejón, D. Sánchez-Portal, E. Artacho, and J. M. Soler, arXiv:physics/9908022.

R. Santamaría and A. Vázquez, J. Comput. Chem. [**9**]{}, 981 (1994).

J. P. Perdew, Phys. Rev. B [**33**]{}, 8822 (1986)

A. D. Becke, J. Chem. Phys. [**84**]{} 4524 (1986).

A. D. Becke, Phys. Rev. A [**38**]{}, 3098 (1988); A. D. Becke, J. Chem. Phys. [**88**]{}, 2547 (1988).

S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. [**58**]{}, 1200 (1980).

The cutoff radii for the finite range atomic orbitals used in this work are (in a.u.): 4.2 for H, 4.1 for C, 3.6 for N, 3.2 for O, and 4.7 for P. For the atoms intervening in the hydrogen bonds longer radii were used (in a.u.): 5.5 for H, 4.6 for N, and 4.2 for O. The shapes of the basis functions and the radii of the double-$\zeta$ orbitals are as in [@Emilio99] with a split-norm parameter of 0.15. The possibility that the finite range of the basis orbitals could be influencing the results has been tested by recalculating the electronic structure of the ordered structure using orbital ranges increased by 30%. The results did not change significantly (bandgap 2.14 eV, HOMO bandwidth 30 meV, LUMO bandwith 265 meV).

The conjugate-gradients method tends to define paths along the bottom of valleys in the energy landscape.

Assuming a canonical energy distribution among the low-frequency modes only (3 per nucleic basis). If the energy is distributed among more modes than the ones assumed the temperature is lower.

S. O. Kelley and J. K. Barton, Science [**283**]{}, 375 (1999).

Y. A. Berlin, A. L. Burin, and M. A. Ratner, Superlattices and Microstructures [**28**]{}, 241-252 (2000).

P. Phillips and H. L.Wu, Science [**252**]{}, 1805-1812 (1991); A. Parisini, L. Tarricone, V. Bellani, G. B. Parravicini, E. Diez, F. Dominguez-Adame, and R. Hey, Phys. Rev. B [**63**]{}, 165321 (2001).

P. Carpena, P. Bernaola-Galván, P. C. Ivanov, and H. E. Stanley, Nature [**418**]{} 955-959 (2002).

T. Holstein, Ann. Phys. [bf 8]{}, 325 (1959).

S. S. Alexandre, H. Chacham, J. M. Soler, and E. Artacho (to be published).

P. Flükiger, H. P. Lüthi, S. Portmann, and J. Weber, Swiss Center for Scientific Computing, Manno (Switzerland); S. Portmann and H. P. Lüthi, Chimia [**54**]{}, 766-770 (2000).
