Effect of helix angles on polaron motion in poly-DNA

Physics Letters A 372 (2008) 6013–6018

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Effect of helix angles on polaron motion in poly-DNA Zhen Qu a , Da-wei Kang a , Yuan Li a , Wen Liu a , De-sheng Liu a , Shi-jie Xie a,b,∗ a b

School of Physics, Shandong University, Jinan 250100, China National Key Laboratory of Crystal Materials, Shandong University, Jinan 250100, China

a r t i c l e

i n f o

Article history: Received 1 June 2008 Received in revised form 20 July 2008 Accepted 24 July 2008 Available online 7 August 2008 Communicated by A.R. Bishop

a b s t r a c t We investigate the effect of helix angles on polaron motion in poly-DNA in the framework of a tight binding model. A DNA molecule is supposed to have a uniform or a random distribution of the helix angles. It is found that a polaron moves faster along a DNA molecule with an overwinding helix than that with an unwinding one. In DNA molecule with randomly distributed helix angles, charge transport is suggested to be field-facilitated tunneling between spatial disordered potential wells. © 2008 Elsevier B.V. All rights reserved.

PACS: 72.80.Le 71.38.-k 87.14.Gg 87.15.He Keywords: Poly-DNA Polaron motion Helix angle

1. Introduction In recent years, deoxyribonucleic acid (DNA) has attracted much attention for its potential application in molecular electronics and spintronics. Direct observation of I–V characteristic in a DNA molecule has been widely investigated. However the argument whether DNA is a conductor [1–4], a semiconductor [5–8], or an insulator [9–11] is still continuing. A number of charge transport mechanisms are proposed to interpret the inconsistent experimental results, among which a polaron hopping mechanism has turned out to be a prospective candidate for modeling the strong coupling between electronic distribution and lattice configuration [12–26]. Some experimental findings also support that a polaron may be the charge carrier in DNA molecules [7,8]. With a one-dimensional SSH(Su–Schrieffer–Heeger)-like model [27–29], Conwell and Rakhmanova found that a polaron or a bipolaron may be the main electronic carrier in a DNA molecule. They obtained that a polaron has a width of 5–7 base pairs. The time required to form a polaron is of the order of picoseconds, which is apparently much longer than that in a conjugated polymer [13]. Wei et al. investigated bipolaron dynamics in a sequence-disordered

*

Corresponding author at: School of Physics, Shandong University, Jinan 250100, China. E-mail address: [email protected] (S.-j. Xie). 0375-9601/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2008.07.072

DNA molecule [18,19]. Zheng et al. showed that on-site energies of the bases have an important effect on polaron motion [20]. Recently, Zhu et al. found that small polaron effects arose from the strong electron-lattice coupling [21]. By a Peyrard–Bishop model, the group of Bishop investigated the effect of intrinsic base pair fluctuation on the dynamics of a polaron in DNA [22–24]. To describe the spatial structure of a DNA molecule, Hennig et al. introduced a three-dimensional tight binding model. They calculated the effect of thermal fluctuations on a polaron motion for a B-DNA molecule [25]. Later, it is found that a polaron in poly(dG)–poly(dC) DNA moves steadily. However, in poly(dA)– poly(dT) a polaron could not move but oscillates in a confined region [26]. Their conclusion seems to explain the experimental result that poly(dG)–poly(dC) exhibits better conductance than poly(dA)–poly(dT) [7]. However they only considered the pristine B-form of DNA molecules in the above calculations. In real situation, the helix, as well as the hydrogen bond (H-bond), may easily suffer the environmental effects such as stress, temperature or relative humidity, which may result in the nonuniform of the molecule structure. For example, some groups made extensive discussions on the helix angle effects in DNA molecules, in which the twist-stretch coupling of the molecule is observed and investigated [30–32]. In this Letter, we will investigate the effect of helix angles on a polaron dynamics by supposing that the helix angles have a departure from its pristine form in a poly-DNA molecule. One kind

6014

Z. Qu et al. / Physics Letters A 372 (2008) 6013–6018

of departure is that the helix angles are uniformly overwound or unwound. The other is a local or a random distribution of the angles. In the next section, a tight binding Hamiltonian in a threedimensional space is introduced and the time-dependent dynamic equations are given. The calculations are presented and the results are discussed in Section 3. Finally in Section 4, a summary is made.

as the base pair movement is treated classically in the Hamiltonian, we separate the electronic part from the lattice to solving Hamiltonian H = H el + H latt + H E . Then we can obtain the dynamics of the base pairs from the classical Newton equation and the electronic state evolution of an electron from the time-dependent Schrödinger equation. They are determined separately by

2. Model and formulae

M n r¨n = −β ρn,n − M n Ωr2 rn + k



A DNA molecule is a double helix structure with a backbone made of alternating sugars and phosphate groups that are combined through covalent bonds. According to Watson–Crick complementary base pair law, adenines in one strand of backbone only pair with thymines (A/T or T/A) in another strand, and cytosines with guanines (C/G or G/C). Bases connect with each other through H-bond. The two coupled backbones appear a double helix structure in a three-dimensional space, thus a double helix structure is built. Charges transport in a DNA molecule through the overlap between itinerant π orbitals in neighboring base pairs. There are two main independent variables to fix the spatial structure of a DNA molecule. One is the transversal radius within a base pair (or the H-bond) R n = R 0 + rn , where rn is the displacement of the bond from its equilibrium value. Another is the helix angle θn,n+1 = θ0 + θn ,n+1 between the neighboring base pairs, where θ0 = π /5 for a B-DNA molecule. A tight-binding Hamiltonian for the itinerant electrons is written as [26] H el =

 (n,s + β rn )C n+,s C n,s

  + D n−1,n t 0 α R 0 ρn−1,n + ρn∗−1,n + |e | E

 N −n 



ρn+k,n+k − N e / N

k=0



  + D n,n+1 t 0 α R 0 ρn,n+1 + ρn∗,n+1 + |e | E

 N −n 



ρn+k,n+k − N e / N

ih¯

∂ ψγ (n, t ) = (n + β rn )ψγ (n, t ) ∂t  − t 0 (1 − αdn,n+1 )ψγ (n + 1, t ) + (1 − αdn−1,n )ψγ (n − 1, t ) 

n −1  d j , j +1 ψγ (n, t ), + |e | E (n − 1)a +

− tn,n+1 (C n+,s C n+1,s + C n++1,s C n,s ),

(1)

where n,s is the on-site energy of the nth base-pair with spin s. For a poly-AT or -GC DNA molecule, n,s can be taken as zero in present work. Term β rn means a correction to the on-site energy due to the radial variation of the base pair. tn,n+1 is the hopping integral of π electrons between site n and n + 1, which is adjusted by the vertical distance between the neighboring base pair planes as tn,n+1 = t 0 (1 − αdn,n+1 ), where constant α measures the strength of electron–phonon interactions. dn,n+1 = {l20 − [ R n2 + R n2+1 − 2R n R n+1 cos θn,n+1 ]}1/2 − a with l0 =



a2 + 4R 20 sin2 (θ0 /2)

and a the lattice constant in a B-DNA molecule. Dependence of the electronic state on the molecule configuration means that there is a strong electron-lattice interaction, which reflects the soft characteristic of a DNA molecule. The H-bond variation will result in a lattice Hamiltonian

1 n

2

M n Ωr2 rn2 +

1 2

M n r˙n2 − k rn ,

(2)

where Ωr is the vibration frequency and M n the mass of a base pair. Term −k rn with k = 2β is introduced to prevent an endopen DNA molecule from collapsing. It has been indicated elsewhere that the helix angles are not sensitive to the itinerant π -electrons [25]. So we only include the H-bond relaxation and fix the helix angles with a given value in our calculations. To investigate the dynamic process of a charge carrier along the DNA molecule, an external electric field E is applied, HE =



(4)

(5)

j =1

n, s

H latt =

,

k=1

|e | Exn (C n+,s C n,s − N e / N ),

(3)

n, s

n−1

where e is an electronic charge and xn = (n − 1)a + j =1 d j , j +1 the position of the nth base pair. N e is the number of itinerant π electrons and N the number of the base pairs of the DNA molecule. In a recent review, the adiabatic and non-adiabatic limit for the one-dimensional Holstein model has been discussed [33]. Here



with D n,n+1 = (1 − cos θn,n+1 )/ l20 − 2R 20 (1 − cos θn,n+1 ). The density matrix ρn,n is defined as

ρn,n =



ψγ∗ (n, t ) f γ ψγ (n , t ),

(6)

γ

with f γ the time-independent Fermi distribution function. Spin index s has been omitted in the equations for clarity. In next section, we adopt the non-adiabatic evolution method as those in Refs. [19,24,26] and solve the coupled differential equations (4) and (5) with a Runge–Kutta method of order eight with step size control [34]. To describe the evolution of electronic state ψγ (n, t ) or the net charge density ρn,n (t ) − N e / N, we define the center of a function x(n) by

⎧ if cos ϕc   0 and sin ϕc   0, ⎨ N ϕc /2π , xc = N (π + ϕc )/2π , if cos ϕc  < 0, ⎩ N (2π + ϕc )/2π , otherwise,

(7)



where

sin 2nπ / N ϕc = arctan n xx((nn)) cos . 2nπ / N n

The calculations are presented with parameters referred in Refs. [25,35]: a = 3.4 Å, R 0 = 10 Å, θ0 = π /5, Ωr = 6.252 × 1012 S−1 and M n = M 0 = 4.982 × 10−25 kg, t 0 = 0.25 eV, β = 0.7 eV/Å and −1

α = 0.3 Å

are taken to fit some experimental data [7].

3. Results and discussions 3.1. Polaron in a B-DNA molecule 3.1.1. Polaron state in a B-DNA molecule Firstly, we give the static results of a pristine B-DNA molecule under a ground state and a positive charged state separately. They are obtained by solving the static eigen equation and the equilibrium conditions self-consistently,

Z. Qu et al. / Physics Letters A 372 (2008) 6013–6018

Fig. 1. H-bond distortion (solid) and net charge density (dash) of a stationary polaron in a B-DNA molecule.

6015

Fig. 2. Time evolution of the localization factors of the electronic states.

(n + β rn )φμ (n) − t 0 (1 − αdn,n+1 )φμ (n + 1) − t 0 (1 − αdn−1,n )φμ (n − 1) = εμ φμ (n), (8)       1 2 rn = k − βφμ (n) + 2t 0 α R 0 D n−1,n φμ (n − 1)φμ (n) M n Ωr2 μ μ   + D n,n+1 φμ (n)φμ (n + 1) (9) where εμ is the stationary eigen-level of π -electrons. Our calculations are based on the assumption that each base pair provides two itinerant electrons, which means a fully occupied ground state [13]. It should be mentioned that, as showed in Eq. (9), the H-bond is dependent upon the electronic occupation or the total number of the itinerant π -electrons. A discussion on this topic has been given in Ref. [36]. Here, we consider a 30-basepair-long DNA molecule, which has 60 itinerant electrons in the pristine state. It is found that there is no bond distortion and all the electronic states are extended. The width of the energy band is about 4t 0 . If one electron is taken out from the highest occupied molecular orbital (HOMO), it is found that a positive charged polaron appears. The HOMO level will shift up from the band and become isolated. By introducing the localization factor of a state  as ξμ = n |φμ (n)|4 , we obtain that the electronic state of the HOMO level is well localized with a localization factor of ξ p = 0.2. As shown in Fig. 1, the net charge density is found to be confined in the region of the H-bond distortion, which has a width of about 1/ξ p ∼ 5 base pairs. If we suppose that the pristine state is halffilled band (N e = 30 in present molecule), which means that each base pair provides one single electron [21], we can also obtain a charged polaron by introducing an electron into or extracting one from the molecule. To keep a consistency with some other investigations [13,18–20], we will consider a positive charged polaron from a full-filled band system in the following. 3.1.2. Polaron dynamics in a B-DNA molecule Firstly, the formation process of a polaron is studied. One electron is taken out from a pristine B-DNA molecule at the beginning. The evolution of the localization factors of all the electronic states are shown in Fig. 2. It is found that the localization factor of the HOMO state (dashed line) increases from 0.05 to 0.2 in 4 ps, which means that this state becomes well localized. While those of all the other states (solid lines) keep nearly unchanged. By checking the charge density distribution, we find that a localized charge appears in the center of the molecule. At the same time, a localized H-bond

Fig. 3. The instantaneous velocity of a polaron in a B-DNA molecule under the field of 0.4 mV/Å (solid), 0.8 mV/Å (dash) and 1.6 mV/Å (dot). The inset shows the center motion of a polaron under the corresponding electric field strength.

distortion forms around the localized charge. In present parameters, it takes about 4 ps to form a polaron, which is much longer than that in a conducting polymer such as polyacetylene, where it needs only about 0.1 ps to form a polaron [29]. This is due to the more massive base pairs in DNA molecules. Yamada et al. investigated the spread of the electronic wavepacket with a similar model [37]. However, the interaction between electron and lattice distortion was not taken into account. Rakhmanova et al. also reported a picosecond-time-scale for a poalron relaxation by using a one-dimensional SSH-like model [13]. Then the polaron motion is studied under an external deriving electric field. A preexisted polaron is centered on site 10 at the beginning. It has been known that, in a conducting polymer such as polyacetylene, a polaron is sensitive to external electric field [38]. Under a strong field, the electron (or hole) may get rid of the confinement of the lattice potential, which will result in the dissociation of a polaron. In a DNA molecule, we find a similar phenomenon that a polaron will dissociate if the field is strong enough. In this Letter, we only consider the situation under a weak field, which means that a polaron does not dissociate during moving. The instantaneous velocity of a polaron is shown in Fig. 3. It is found that polaron velocity does not keep constant but oscillates around an average value (average velocity) during moving, which is a detailed finding of the present investigation. We attribute the ve-

6016

Z. Qu et al. / Physics Letters A 372 (2008) 6013–6018

locity oscillation to the phonon excitation during the polaron motion. The oscillation period is obtained to be 1 ps, which is equal to the bare lattice vibration period T = 2π /Ωr ≈ 1 ps. At the same time, the polaron shape also appears an oscillation with the same time period. Ref. [39] also indicated the phonon emitting during a polaron moving in a conjugated polymer, where they obtain the conclusion from the lattice dynamics. It is also found that a polaron moves at a higher average velocity if under a stronger field. Center motion of a polaron is shown in Fig. 3 inset. It is clearly demonstrated that a polaron passes through longer distance if it is driven by a stronger electric field. It should be mentioned that, polaron velocity oscillates with small amplitude, which is generally less than 15% of the average velocity. So, the curves in Fig. 3 inset do not show obvious slope changing. 3.2. Angle effect on polaron motion 3.2.1. DNA molecule with uniformly distributed helix angles Now, we consider a uniform twist structure of a DNA molecule by setting θn,n+1 = θ0 + θ0 , where θ0 > 0 denotes an overwound structure and θ0 < 0 an unwound one. In present model, an overwound structure makes the neighboring base pair planes close to each other, which will result in an increasing of the overlap between π orbitals. So, a polaron tends to be extended with the

Fig. 4. Dependence of the polaron width on θ0 . The inset shows the dependence of the average velocity of a polaron on θ0 under electric field E = 1.0 mV/Å.

overwinding of the helix, as shown in Fig. 4. The dependence of the average velocity of a polaron on θ0 under the electric field E = 1.0 mV/Å is shown in Fig. 4 inset. It is found that a polaron moves faster in a DNA molecule with an overwinding structure than that with an unwinding one. We can also say that an extended polaron moves with a high velocity. 3.2.2. DNA molecule with randomly distributed helix angles Instead of a uniform structure, we further take a random distribution of the helix angles by supposing a canonical distribution of the helix angles with a maximum deviation of θa . At first, we consider a local twist θn ,n+1 = θ0 δn, K . When a polaron crosses the twist defect, it will suffer an apparent scattering. It is found that an overwinding defect with θ0 > 0 corresponds to a potential well. A polaron attends to be confined in the well. The dynamical process for a polaron to cross an overwinding defect (θ0 = 0.015) between site 18 and 19 (or at K = 18) is shown in Fig. 5. It is found that the polaron velocity appears to decrease when it passing through the well defect, which means that the polaron needs to accumulate energy to overcome the confinement of the well. At the same time, the polaron becomes a little thin as it passing through the well. For an unwinding defect (θ0 = −0.015), an opposite picture is obtained. The velocity and the width of the polaron appear to increase when it passing through the barrier defect, as shown in Fig. 6. Then a DNA molecule with randomly distributed helix angles is taken into account. There is a preexisted polaron in the left side of the molecule at the beginning. As the molecule is disordered, we have to adjust the strength of the field to drive a polaron. Of course, the driving field is related to the maximum deviation θa . Fig. 7 gives the instantaneous velocity (solid) and the width (dash) of a polaron during it moves through a DNA molecule with randomly distributed helix angles. The maximum deviation and the electric field strength are θa = 0.04 and E = 1.8 mV/Å respectively. Due to the random lattices, obvious well and barrier defects appear in the molecule. So, the period oscillations of the polaron velocity and polaron width are covered by a large random one. The polaron velocity, as well as the polaron width, evidently increases when it passing through barrier defects. Center motion of the polaron is shown in Fig. 7 inset. The polaron movement is apparently different from that in a structural-ordered DNA molecule. The polaron has a long-time-rest in a well defect, which means that it needs to accumulate energy to pass through the neighboring barrier. Driven by the electric field, net charge density decreases in the anterior well and increases in the posterior well with time, but keeps nearly unchanged at the barrier.

Fig. 5. The instantaneous velocity (solid) and the width (dash) of a polaron during it moves through an overwinding defect (θ0 = 0.015) at K = 18.

Z. Qu et al. / Physics Letters A 372 (2008) 6013–6018

6017

Fig. 6. The instantaneous velocity (solid) and the width (dash) of a polaron during it moves through an unwinding defect (θ0 = −0.015) at K = 18.

Fig. 7. The instantaneous velocity (solid) and the width (dash) of a polaron during it moves through a DNA molecule with randomly distributed helix angles. The inset shows the center motion of the polaron during moving.

Polaron motion in DNA molecule with a random distribution of helix angles is a discrete migration between spatial-disordered well defects. The velocity and the shape of a polaron come through an evident changing during moving. Charge transport in a structural-disordered DNA molecule is considered to be fieldfacilitated tunneling between potential wells [19]. It should be mentioned that, according to Ref. [25], random double helix structure does not significantly affect the mobility of the polaron (see Fig. 6 of Ref. [25]). The authors treated the helix angles as independent degrees of freedom and checked the response of the helix to the itinerant π -electron. However, we give a different picture for describing the effect of helix angles on polaron motion, in which helix angles are treated as invariable structural parameters. The distribution of helix angles only provide a lattice potential for the moving polaron but not change with the charge distribution. Additionally, we consider a strong longitudinal electron–phonon coupling [16–20] between neighboring base pairs. So, it is demonstrated in our calculations that helix angles significantly affect polaron motion in DNA. Some other works also support that helix angles between neighboring base pairs affect significantly on the electronic state of DNA molecules [40,41]. 4. Conclusions In summary, we emphasize the effect of helix angles on polaron motion in the framework of a tight binding model. H-bond vibrations and some details of spatial structure in DNA molecules are

taken into account. Firstly, we give the polaron picture in a pristine B-DNA molecule. A stationary polaron in a B-DNA molecule extends over about ∼ 5 base pairs. The relaxation time for forming a polaron is about ∼ 4 ps, which is much longer than that in a conducting polymer such as polyacetylene. Due to the phonon excitation, polaron velocity oscillates around an average value during moving, which is a detailed finding of the present investigation. Then we investigate the dynamics of a polaron in a DNA molecule with a uniform twist structure. It is found that polaron tends to be extended with the overwinding of the helix. Furthermore, the polaron moves faster in a DNA molecule with an overwinding helix than that with an unwinding one. We further investigate the case of randomly distributed helix angles by supposing a canonical distribution of the helix angles with a maximum deviation of θa . Polaron motion is a discrete migration which indicates a fieldfacilitated tunneling mechanism. Acknowledgements The authors would like to acknowledge the financial support from the Research Fund for the Doctoral Program of Higher Education (No. 20070422058) and the Provincial Science Foundations of Shandong (Grant No. Z2005A01). References [1] H.W. Fink, C. Schonenberger, Nature 398 (1999) 407.

6018

Z. Qu et al. / Physics Letters A 372 (2008) 6013–6018

[2] L. Cai, H. Tabata, T. Kawai, Appl. Phys. Lett. 77 (2000) 3105. [3] B. Hartzel, B. McCord, D. Asare, H. Chen, J.J. Heremans, V. Sogomonian, Appl. Phys. Lett. 82 (2003) 4800. [4] B. Hartzel, B. McCord, D. Asare, H. Chen, J.J. Heremans, V. Sogomonian, J. Appl. Phys. 94 (2003) 2764. [5] D. Porath, A.W. Ghosh, S. Datta, Nature 403 (2000) 635. [6] H. Watanabe, C. Manabe, T. Shigematsu, K. Shimotani, M. Shimizu, Appl. Phys. Lett. 79 (2001) 2462. [7] K.H. Yoo, D.H. Ha, J.O. Lee, J.W. Park, J. Kim, J.J. Kim, H.Y. Lee, T. Kawai, H.Y. Choi, Phys. Rev. Lett. 87 (2001) 198102. [8] H.-Y. Lee, H. Tanaka, Y. Otsuka, Appl. Phys. Lett. 80 (2002) 1670. [9] E. Braun, Y. Eich, U. Sivan, G. Ben-Yoseph, Nature 391 (1998) 775. [10] A.J. Storm, J. van Noort, S. de Vries, C. Dekker, Appl. Phys. Lett. 79 (2001) 3881. [11] Y. Zhang, R.H. Austin, J. Kraeft, E.C. Cox, N.P. Ong, Phys. Rev. Lett. 89 (2001) 198102. [12] E.M. Conwell, S.V. Rakhmanova, Proc. Natl. Acad. Sci. USA 97 (2000) 4556. [13] S.V. Rakhmanova, E.M. Conwell, J. Phys. Chem. B 105 (2001) 2056. [14] E.M. Conwell, D.M. Basko, Synth. Met. 137 (2003) 1381. [15] J.H. Park, H.Y. Choi, E.M. Conwell, J. Phys. Chem. B 108 (2004) 19483. [16] E.M. Conwell, J.H. Park, H.Y. Choi, J. Phys. Chem. B 109 (2005) 9760. [17] E.M. Conwell, S.M. Bloch, J. Phys. Chem. B 110 (2006) 5801. [18] J.H. Wei, L.X. Wang, K.S. Chan, Y.J. Yan, Phys. Rev. B 72 (2005) 064304. [19] J.H. Wei, X.J. Liu, J. Berakdar, Y.J. Yan, J. Chem. Phys. 128 (2008) 165101. [20] B. Zheng, J. Wu, W.Q. Sun, C.B. Liu, Chem. Phys. Lett. 425 (2006) 123. [21] J.X. Zhu, K.Ø. Rasmussen, A.V. Balatsky, A.R. Bishop, J. Phys.: Condens. Matter 19 (2007) 136203.

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41]

S. Komineas, G. Kalosakas, A.R. Bishop, Phys. Rev. E 65 (2002) 061905. G. Kalosakas, K.Ø. Rasmussen, A.R. Bishop, J. Chem. Phys. 118 (2003) 3731. G. Kalosakas, K.Ø. Rasmussen, A.R. Bishop, Synth. Met. 141 (2004) 93. D. Hennig, J.F.R. Archilla, J. Agarwal, Physica D 180 (2003) 256. D. Hennig, E.B. Strarikov, J.F.R. Archilla, F. Palmero, J. Biol. Phys. 30 (2004) 227. W.P. Su, J.R. Schrieffer, A.J. Heeger, Phys. Rev. Lett. 42 (1979) 1698. W.P. Su, J.R. Schrieffer, A.J. Heeger, Phys. Rev. B 22 (1980) 2099. W.P. Su, J.R. Schrieffer, Proc. Natl. Acad. Sci. USA 77 (1980) 5526. T.R. Strick, J.F. Allemand, D. Bensimon, A. Bensimon, V. Croquette, Science 271 (1996) 1835. J.F. Marko, Europhys. Lett. 38 (1997) 183. G. Jeff, B. Zev, N. Marcelo, U.L. Mai, R.C. Nicholas, B. Carlos, Nature 442 (2006) 836. H. Fehske, S.A. Trugman, in: A.S. Alexandrov (Ed.), Polarons in Advanced materials, Springer, 2007. R.W. Brankin, I. Gladwell, L.F. Shampine, RKSUITE: Software for ODE IVPS, http://www.netlib.org. M. Barbi, S. Cocco, M. Peyrard, Phys. Lett. A 253 (1999) 358. X.T. Gao, X. Fu, L.M. Mei, S.J. Xie, J. Chem. Phys. 124 (2006) 234702. H. Yamada, E.B. Strarikov, D. Hennig, Eur. Phys. J. B 59 (2007) 185. X.J. Liu, K. Gao, D.S. Liu, S.J. Xie, Phys. Rev. B 74 (2006) 17230. J.F. Yu, C.Q. Wu, X. Sun, K. Nasu, Phys. Rev. B 70 (2004) 064303. R.D. Felice, A. Calzolari, E. Molinari, Phys. Rev. B 65 (2001) 045104. F.C. Grozema, L.D.A. Siebbeles, Y.A. Berlin, M.A. Ratner, Chemphyschem. 6 (2002) 536.