Soliton dynamics in alternating trans-polyacetylene and in stacked systems

Journal of Molec~~r Structure ~T~eoc~m}, 188 (1989) 231-260 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

231

SOLITON DYNAMICS IN ALTERNATING ~~A~~-~OLYACETYLENK AND IN STACKED SYSTEMS

W. FORNER, J. LADIK, D. HOFMANN,

M. SEEL and A. GODZIK*

Chair for T~oret~~l Che~~t~ and ~aborato~ of the Nutiona~ Fou~atio~ for Cancer Research at the Friedrich-A~x~~er-U~~vers~ty Er~ngen-N~rnberg, D-8520 Rr~a~gen, E~er~un~tr. 3 (F.R.G.) F. MARTIN0 Physics ~epurt~ent,

Portkcnd State university, Portland, OR 972~7-~7~1 (U.S.A.)

(Received 24 May 1988)

ABSTRACT The influence of su~titution~ disorder, of random forces and of energy dissipation on solitons in alternating bran-polyacetylene using the Su-Sch~effer-Hager H~il~ni~ is discussed. Double ionization in polyene chains was also investigated. Further sol&on dynamics within a full Pariser-Parr-Pople Hamiltonian for trans-polyacetylene are presented. The presence of conformational solitons in stacked systems is shown numerically in a polyformamide model system. Possible relations to DNA and carcinogenesis are briefly reviewed.

INTRODUCTION

Solitons are represented by a special class of solutions of non-linear partialdifferential equations. Soliton solutions are characterized by the remarkable property that they are excitations with an infinite life time. They do not change their shape when travelling through the system, upon reflection at system boundaries, or upon collisions with other solitons f 1 J. Non-linearity of the forces occurring is a necessary, but not sufficient, condition for the existence of soliton solutions. However, if a potential has two (@*potential [ 2 ] as in polyacetylene) or more (Sine-Gordon potential [ 31) degenerate energy minima, the existence of topological soliton solutions is guaranteed. These solitons are zone boundaries between phases belonging to different energy minima (vacuum states in field theory) [4]. Once excited such a soliton cannot be removed from the system. Its movement may be stopped but it cannot vanish, because it connects two physically different phases of the same energy. Topo*Present address: Institute for Experimental Physics, University of Warszawa ul. Zwirki i Wigury 93,02-089 Warszawa, Poland.

0166-1280/89/$03.50

Q 1989 Elsevier Science Publishers B.V.

logical solitons are often compared with a knot in an infinite string, which cannot be removed without cutting the string [5]. In contrast to topological solitons, conformational solitons are not bound to exist forever once excited, because the same phase exists on both sides of conformational solitons. The famous water wave observed by Scott-Russel [6] is a typical example of a conformational soliton. In this respect it is necessary to distinguish between solitons and solitary waves. A conformational soliton is an excited state of a system (e.g. geometrical distortion) which, in principle, can relax, but travels around without change of shape for infinite time due to the non-linear forces acting in the system. A solitary wave is a slowly dispersing soliton. This means it is a perturbation which survives much longer than would be expected from linear theory (e.g. dispersing phonon packets) due to the action of non-linear forces. Soliton concepts are involved in many fields of physics, chemistry, and biochemistry. To mention only some examples, solitons have been introduced for the explanation of phase changes [ 71, rotation of CHz groups in polyethylene chains [8], energy transport and storage in proteins [9], coupling of base movements to the sugar puckering in DNA [lo], phase changes between B-, A- and Z-DNA [ll]. In this review we will deal with two distinct examples. One is a topological soliton proposed for al~rnating ~~~~s-polyacetylene (PA) [ 121, the other is a conformational soliton or probably a solitary wave in stacked systems [ 13,14]. Let us look first at trans-polyacetylene which is, because of the papers of Su et al. [ 121, one of the most investigated cases of a soliton in solid-state physics. trans-Polyacetylene (CH),, contains zig-zag chains of CH units (Fig. 1) with alternating long and short bonds. The latter was confirmed only recently [ 151 by X;ray and ab initio studies on oriented material leading to a difference of 0.08 A between the two bond lengths (corresponding to a projected dimerization parameter of 0.04 A). Thus in polyacetylene there are two different dimerization phases A and B (Fig. 2) with the same energy. The soliton, or phasekink, represents the zone boundary between two such phases and is therefore a topological soliton. Due to the chemical structure of the chain the neutral soliton is accompanied by an unpaired spin (Fig. 2 ). Therefore, just a single soliton can occur only in a chain with an odd number of CH units (in fact any odd number of solitons is topologieally possible in such a chain). An evennumbered chain can carry only pairs of solitons. Singlet coupled pairs occur not as pairs of neutral solitons but as pairs of charged solitons, one positively and one negatively charged [ 161. Since theory predicts that an unpaired spin tends to localize in longer PA chains and does not lead to an equidistant geometry with a complete delocalization of the unpaired spin, there is at least one soliton present in every oddnumbered PA chain. Indeed, electron parama~eti~ resonance (EPR) measurements indicate in pristine t-PA the presence of unpaired spins (400 ppm >

233

(a)

Fig. 1. ~e~nition of the coordinates used in Grady-polyace~lene in (a) an evident (b) an ideally dimerized chain. [ 171.

chain, and

Therefore, there are odd-numbered chains in the material and thus neutral solitons are present. Su et al. [ 121 introduced the soliton concept for PA to explain three experimental features. First there is a mid-gap absorption observed in PA [ 18 1. However, only in the Hiickel-type SSH model are the neutral soliton levels at mid-gap, due to symmetry. In any more-realistic model which inco~orates electron~lectron interactions, the soliton levels are shifted towards the band edges [19,20]. Indeed Weinberger et al. [Zl] found experimentally that the mid-gap absorption cannot be due to neutral solitons; its origin is still an open question. Weinberger et al. [ 211 discussed remaining catalysts as well as charged solitons as its origin. The second most remarkable feature of PA is that its electric conductivity changes over 12 orders of ma~it~de upon doping 1221 and that in the lightly doped regime the charge transport occurs via spinless carriers [23]. An additional electron added to a chain with a neutral soliton would occupy the soliton level and a spinless negatively charged quasi-particle is obtained, which is free to move along the chain [ 121. However, this is insufficient to explain the conductivity of t-PA, since the charge must get from one chain to the other to finally reach the electrode [24]. A soliton cannot hop from one chain to the other because this would require dimerization phase changes over a wide range of both chains f24]. Kivelson [25] pointed out that the electron of a negatively charged soliton on one chain can hop to a neutral soliton on a neighboring chain. However, on one hand even in the lightly doped

F’ig.2. Sk&h of several different configurations in t~Rs-poi~a~tyiene: (a) dimerization phase A; [X-r) d~~e~iz~tion phase B, (c) soliton localized on one site within the chain; (d) solitm locatized at the chain end; and (e) polaron bcalized on two sites.

regime all initially neutral solitons should already be charged, as pointed out by Brddas et al. [ 26 1. On the other hand, since the initially charged soliton will carry the spin of the initially neutraf soliton after the hopping, this process would be accompanied by a spin tra~spo~ in the opposite direction of the charge t~~spo~~ Another model for the spinless charge transport is provided by the polaron [27 ] and bipola~n [28 ] concepts. The third feature Su et al. [I2 1 wanted to explain was the observed narrowing of the EPR lines [111 which can be explained by invoking a moving spin and, therefore, a sol&on. However, this line-narrowing was also measured for co-polymers of acetylene and CO [ 29 1 where soliton migration is not possible, as our previous calculations have shown [ 30 3. However, at least there are solitons present in t-PA, but to what extent they can really explain the three above-mentioned properties of t-PA still remains, in our opinion, an open question. meanwhile other properties of t-PA have been explained using the soliton model, However, for this purpose the simple Hiickell-type SSH model is not enough and one has to in~o~or~t~ e~eGtron~~e~tro~ infractions. Using the ~~~~R-te~hni~ue~ the alternating negative and positive spin densities have

been measured in t-PA [31]. Such a spin polarization naturally cannot be explained by a Hiickel model, but requires electron-electron interactions [32]. In addition, 13C-NMR line shapes could only be explained on the basis of Pariser-Parr-Pople studies [ 33 1. Further photo-induced low- and high-energy absorptions have been found [ 341. The low-energy peak is assigned to photogenerated charged soliton pairs [34 J which would absorb mid-gap in the SSH model. The high-energy peak is a matter of controversy. Bishop et al. [35] assign it to a breather vibration left behind between a separating pair of charged solitons, Wang and Martin0 [36] found the oscillating bound state of the soliton pair responsible (exciton breather), while Su (371 and Kivelson and Wu [38] attributed it to a triplet coupled pair of neutral solitons. Recent experimental work [ 34b] reveals that the high-energy photoabsorption could be due to neutral solitons created by an intrachain excitation. This leads primarily to charged solitons which recombine to pairs of neutral solitons after 0.5 ps. However, all models require inclusion of electron-electron interactions in order to yield reliable results. Our interest was first of all in the study of substitutional impurities 1301 and environment, simulated by random forces and a dissipation term [391 on soliton dynamics in the SSH model. Another interesting point was the investigation of double ionization [40]. Recently, we have implemented the full PPP model [41] and reported the first soliton dynamics within this model [ 201, where we observed a considerable change in the soliton properties upon inclusion of electron-electron interactions. The photodynamic study by Wang and Martin0 [ 361 dealt with an extended Hubbard Hamiltonian which is not necessarily enough to describe t-PA sufficiently well [ 421. Another important field in which solitons or solitary waves most probably play an important role is the as yet unsolved problem of chemical carcinogenesis, which is the reason for 80% of all tumors [43 ] . Simple s~tistical considerations show that the probability of removing a nucleo-histone molecule from DNA to activate an oncogene by local effects of carcinogens is far too small to explain carcinogenesis [ 441. Long-range effects of carcinogens have been discussed in the literature. One possibility is the change of the tertiary structure of DNA and the surrounding environment (water, ions) through binding of a bulky metabolized carcinogen [ 451. Via internal charge transfer between DNA and the histone both biopolymers can become conducting [ 46 1. If the binding of a carcinogen causes a filling of the valence band of one or both of the two chains then the dispersion and polarization forces between them will decrease strongly [ 47,481. This, together with the broken periodicity of the DNA backbone [ 451 can cause a release of the histone from DNA and thus deblocking of an oncogene. Recently, a third mechanism was proposed [ 13,14 J which is based on solitary waves in the nucleotide base stacks of DNA. A bulky ultimate carcinogen bound to DNA will certainly cause a structural change which is coupled to the

236

change of the electron-electron interaction between the stacked bases. These changes in the base stack are localized around the bonding site. Through the action of repair enzymes, the carcinogen can be released from DNA within a comparatively short time in vivo. However, due to the occurrence of non-linear forces in DNA (caused by these two coupled changes) it is possible that the geometrical change will not relax, but travel as a solitary wave through a long range in the base stack. This solitary wave can influence the interactions between DNA and a histone [ 13,141. Ladik and &ek [ 131 have generalized the SSH Hamiltonian to three geometrical variables in order to investigate the probability of this mechanism. To study the problem numerically an extension of their Hamiltonian (one orbital per site) was performed to be able to treat the constituent molecules as individual entities [491. However, electron-electron interactions have not been explicitly included as in ref. 13, but implicitly via the parameters. These parameters have been taken from ab initio Hartree-Fock calculations. As a first simplified model stack, we considered formamide as the constituent unit, because the same heavy atoms as in formamide occur also in the nucleotide bases. METHODS

Time-simulation procedure

In all cases considered in this paper o-a separability is assumed, i.e. the x electrons are treated explicitly while the o electronic effects are represented by a potential. It is assumed that due to their delocalized nature the x electrons may give rise to long-range effects, while the o electrons are assumed to be localized so that the potential energy connected with them can be written as a sum of pair or bond contributions. Thus the total Hamiltonian is given by jk&+V,+V’+T

(1)

where & is the a-electron Hamiltonian, V, the potential due to the aelectrons (including also the 0-x coupling effects), T the kinetic energy, and V’ a potential term which must be added in all cases in order to stabilize the system in its presumed equilibrium geometry. Furthermore, it is assumed that the Born-Oppenheimer approximation is valid. Therefore, in each time step of the simulation the electronic state is determined as time independent leading to an energy term E,. Thus the total potential in which the nuclei are moving is given as V=E,+V,+V

(2)

The nuclear motion is treated in classical approximation using Newton’s equations of motion. If qi denotes some set of coordinates of the moving nuclei or units the classical equations of motion are

237

Fi=fi

&z-g

(3) z

where fi is a factor depending on the nature of qi, namely a mass mi if qi is a simple linear coordinate or a moment of inertia 0, if qi is an angular coordinate. The dots denote the total time derivative. To solve eqn. (3) one can introduce a suitably small time step, z, and proceed in a simple one-step procedure. The value of z must be small enough such that all conservation laws (especially on total energy) hold for the whole time period within a prescribed numerical accuracy. Given a set of coordinates qi (t) and velocities qi( t) at a time t one can write approximately b’i=(~i(t+z>-~i(t))/~

4itt)=

(44

(Qi(t+z)-qi(t)

)/7

From eqn. (4) the new sets of coordinates at t+ ;li(t+z)=~i(t)+b’i(t)‘z Qi(t+7)

=qi(E)

(4b)

7

can be computed as (5a)

+ii(t).7

(5b)

The accelerations ii(t) are computed using Newton’s law [ eqn. (3) ] by differentiation of V(t) with respect to qi. In the first simulations on t-PA published by Su and Shrieffer [ 161 and also in those performed in our laboratory, instead of eqn. (5b), the equation qi(t+7)=qi(t)+ii(t+7)*7

(6)

was used. In the Appendix we show, using the case of a simple exactly soluble problem (harmonic oscillator), that eqn. (6) leads to more reliable results than eqn. (5b). However, if 7 is small enough, eqns. (5b) and (6) should lead to the same solution of eqn. (3). Another possibility is to assume the forces Fi to be constant in 7 and to integrate directly the equations of motion. This leads to ii(t+7)

=~i(t)+~i(t)‘7

(74

(7b) However, simulations on t-PA using eqn. (5a) together with eqns. (6) and (7) show that with 7= 1.25 x lo-” s the total energy remains sufficiently constant and oscillates around its initial value in the former case. In the latter case, however, the total energy rises monotonically leading to completely meaningless results after 40 7 [ 201. This implies that to use eqn. (7) 7 must be chosen much smaller than 1.25 x lo-l5 s. Therefore, the use of eqn. (5a) together with eqn. (6) must be preferred. In the Appendix a comparison of the results of

238

these three approximations with the exact solution in the simple case of an isolated classical harmonic oscillator is shown for illustration. The electronicHamiltonianand the potential terms of trans-polyucetylene Most of the studies on t-PA reported here were performed using the SuShrieffer-~eeger (SSH) Hamiltonian for t-PA which is essentially of Hiickel type. %Xitself can be written easily in second quantization language. Since in practical calculations we must deal with its matrix representation in the basis space of thep, atomic orbitals at the carbon atoms, we prefer this motion. First the geometric variables must be defined. Following SSH we introduce the physical-displacement coordinates, tii, of each CH unit from its position in the equidistant chain along the chain axes. This is shown in Fig. 1. From these the staggered coor~nates, vi, are de~ned by ly;+-l)“+“ui

(3)

Obviously the choice of v/i= u. for all i corresponds to the ideally dimerized chain. The experimental value for u, is 0.04 A [ 151. However, in many studies, in order to reduce the necessary chainlen~h, u, = 0.01 A is chosen which leads to a soliton width of ca. 3 lattice sites. Using the v/i coordinates, a soliton is indicated by a change in sign. Dimerization phase A is indicated by v/i= u. and phase B by vi= - h. The position of the solution is at the zero of v (Fig. 2 ). In a chain with an odd number of CH units and vi = u. the soliton is located at the right-hand end of the chain. This is called an end-kink geometry, The end-kink is usually used as the starting geometry for a simulation because it carries potential energy. The equilibrium position for a soliton is in the middle of the chain [ 121. The SSH-Hamiltonian matrix elements can be written as

for an open chain. The diagonal parameter I, is set to zero for carbon atoms and the resonance (sometimes called “hopping” [ 121) integrals are expanded in first order about the equidistant chain bond length [ 121 Br,r+l = --PO-+ t--1)‘(rv,+w,+I)~ P

l=-Po-(-1)’

(wr+vr-l)~

(lOa) (lob)

j3:;: 2.5 eV and the electron-phonon coupling constant cy is 4.1 eV A-” [ 121. The energy levels and molecular orbital (MO) coefiei&ts are obtained by diagonalizing the Hiickel matrix [eqn. (9) ] HSSHC.-E.C. c- t E

(11)

239

From the MO coefficients, constructed.

the charge-density

bond-order

matrix can be

where the o p are the occupation numbers of the MOs for the two possible spin orientations. From eqn. (12) EzSH can finally be computed EzsH =IC P,., H;:H F.8

(13)

which is equivalent to the sum of all the occupied eigenvalues (energies), r p. To include electron-electron interactions explicitly at the semiempirical level, the Pariser-Parr-Pople (PPP) H~iltonian 1411 has been successfully applied for many different purposes [ 50 J. In this case the two-electron integrals, ‘y+g, are introduced. According to our previous experience 1201 the Ohno approximation [ 511 was used for these integrals in this work. y,,=e2 [4e*/(y,,+ys~)2+R~sl-‘/2

04)

where e is the elementary charge and R, the distance between units r and s. l/rr is the so-called on-site Hubbard constant, given by [ 521 Y,,=&--A,

(15)

where 1, is the ionization potential and A, the electron affinity of unit r in the appropriate valence state of those atoms which contribute an orbital to the unit. The Fock matrix elements for spin o (CT=cy, P; ;1= a! if CT=/? and vice versa) in the unrestricted Hartree-Fock (UHF) app~ximation is given as

P$=F

OFjCRCQQi;

P,ZPp”,+P&

(17)

and z, being the charges of the ionic cores after removal of the n-electrons. The one-electron part is then HE=

[

-4-T

ztyrt(1-&l

1

4s +H:fH (1--J,)

(18)

The MO coefficients are obtained from the self-consistent solution of the UHF equations

F’@‘=E;cP

(19)

and the n-electron energy is given by

(20) where the second sum represents the repulsion between the ionic cores. Following SSH [ 121 the o-electron energy is approximated by a harmonic potential for each o bond

The constant K is determined by the condition that v/i= u, should be the librium geometry and thus an energy minimum. However, it turned out that the system described by & + V is unstable respect to lattice shrinking [ 531. This is true at both the SSH and the levels. Therefore the term v’ must be introduced to stabilize the chain. is most easily done by using a linear ansatz for V’ leading to [ 531

equi-

V’=-A(yl+(-l)N~N)

(22)

with PPP This

To determine K and A let us consider a geometry &=u+

(-l)‘(i-l)ca

(23)

where u is the dimerization parameter and the second part leads to a shrinking of each bond by a. Requiring that E,+ V,+ V’ should be minimal for u= U, and a = 0 leads to the formulae for K and A given in refs. 20 and 30 (note that for K in ref. 30 a printing error occurred which is corrected in ref. 20). A probably more convenient way to avoid lattice shrinking is to keep some units fixed at both chain ends. Before describing the gradient computation, which is straightforward for V, and V’ but not for E,, let us turn first to the description of the formamide stack system. This seems to be convenient because the gradients of E, can be described for both models on an equal footing. Formamide stack The geometrical degrees of freedom for each molecule in each stack are shown in Fig. 3 in the case of a dimer [491. The nth molecule in a stack is described by {z,, 0,, qn} where z, is the shift of the molecule from its equilibrium position (n- l)z, along the axis of the stack (z), and 8, is the tilting angle of the molecular plane, As can be seen in Fig. 3 0, is the angle between the CH bond and its projection onto a plane perpendicular to the z axis intersecting the axes at ( n - 1 )zo + z,. (n - 1) Qiois the rotation angle of the nth molecule around the 7

241

Fig. 3. Geometry of a formamide dimer (used for the determination of the parameters).

Fig. 4. Projection of a formamide dimer in its presumed equilibrium geometry (all angles EO”, except pO)along the helix axes (through the H atoms at the CO groups).

axis in equilibrium where the CH bond in molecule 1 coincides with the x axis. v)~is the distortion out of equilibrium. z. and &, were chosen as 3.36 A and 36’) respectively, resembling the situation in B-DNA [ 541. The internal (rigid) geometry of the molecules is shown in Fig. 4. The kinetic energy of the system is given approximately [ 491 by (24)

242

Equation (24) is an approximation [as is eqn. (3) ] in this case because the coordinates are not ~~en~cular [ 551. However, since only small deviations from equilibrium in @and 9 are considered, eqn. (24) should be a reasonable approximation. M is the mass of a molecule (45 g mol-’ ), e,= 115.6 g A2 and 13p)= 162.4 g A” are the moments of inertia with respect to the two angles. The x electrons were treated in the Hiickel approximation and the necessary parameters were taken from the minimal basis set (STO-3G [ 561) HartreeFock calculations using the GAUSSIAN 74 [ 56 ] program package. In order to include implicitly the electron~lectron interactions, the converged Fock-matrix elements between pr AOs were chosen as the parameters. The x-electron Hamiltonian can be divided into a geometry-independent part %&which contains the diagonal parameters as well as resonance integrals between atoms inside one molecule, and &‘, the interaction part which is restricted to ~rst-neighbor interactions.

(25) In matrix representation

&J is given by

where g and h run over the Kcenters of a molecule. From ab initio calculations on the monomer, the parameter values are 1491 cll,= -0.139 H (hartree), C+= - 0.340 H, LYE=-0.184 H, j&= -0.195 H, and flco= -0.362 H. The interaction part of the Hamiltonian is (27) The resonance integrals &r are geometry dependent. With the definitions ~2, = (z,+~ -2,)

; de, = (en+, -4)

; 4h

= (P~+~ -tfd

(33)

the resonance integrals are expanded in a Taylor series of order p

To obtain the coefficients tijkBh,4,OOformamide dimers in different geometries chosen randomly in ( 1AZ I< 0.5 A, 1A81 < 15’) 1Ay, I < 15 ’ ) were computed at the ab initio level and the Taylor series [eqn. (29) ] was fitted to the corresponding Fock-matrix elements by a least-squares procedure [491. Usingp= 6, an error smaller than ca. 0.3 meV for ail points could be obtained. The Hiickel matrix can now be constructed for each time step and diagonaiized to obtain the P matrix [eqn. (12) 3 and the n-electron energy [ eqn. (13)] in the same way as described above. The potential energy associated with the aelectrons is expanded into a Taylor series in the same way as the resonance integrals.

243

(30) For the calculation of the coefficients the same 400 dimer geometries as above were used. For each dimer the z-electron energy using %z was calculated and substracted from the total ab initio SCF energy. Then the series [eqn. (30) ] was fitted to the difference using the same least-squares procedure. The accuracy of the fit applying p=6 was better than 0.009 mH. Thus En+ V, represents the HF energy of the system. The differentiation of V, to obtain the gradients is straightforward. However, at the HF level, the formamide dimer in our presumed “equilibrium” geometry is unstable. Therefore it was again necessary to introduce a linear potential term V

v’=A,(z,

-~~)+AB(~N-~~)+A~(Q)N-Q)~)

The constants A,, A@,A, are defined by the requirement that for a given stack the gradients in the presumed equilibrium geometry should vanish, just as in t-PA. In a 30-unit stack this leads to A ,=1.31 mH/B, Ae=3.45 mH and A,=Z.OO mH [49]. The explicit introduction of electron-electron interactions at the PPP level is under programming in our laboratory. Gradients

of the

z-electron e’nergy

In the first numerical study on soliton dynamics in t-PA, Su and Shrieffer [ 16 ] used a numerical procedure to obtain the gradients of the n-electron energy. They introduced a small shift Swi in one of the coordinates, constructed the new Hiickel matrix and recalculated the energy E,(&) by diagonalization. The gradient is thus approximated by (31) To our knowledge, in all other dynamical studies on t-PA this method, or variants of it, has been used (see, for example, refs. 35, 57, 58). However, this procedure requires (N+ 1) matrix diagonalizations (or SCF iterations in the PPP case) for each time step considered and is, therefore, very time consuming, In our simulations we have computed directly, within the SSH model, the energy shift due to the shift &vi with the help of matrix pe~urbation theory of first order [30]. In a recent paper [59] we showed that this method besides being much faster than that of Su and Shrieffer [16] in fact yields exact gradients.

244

In ref. 59 it was shown that if the x-electron energy in the Hiickel case is given by the expression (32)

E, = C P, Hr, J.8

where P is built from the eigenvectors of H, the exact gradient is given by (33) All terms containing derivatives of P, cancel out exactly. This leads to the very simple formula ~=2~(-l)i[Pi,i+l(l-~~~)-Pi,i_l(l-~~l)] Wi

(34)

In the PPP case the energy is given by (35)

Again in ref. 59 it was shown that in aE,/ayi also, all the terms containing derivatives of PFs cancel out exactly. This leads again to the very simple and exact formula dEPPP

1

a=2 Wi

C (l-S&) 1 t

-2( (P$)*+ (Pi)‘)]

2

[Pii(Ptt-2

Zt) +Ptt (Pii-

Zi)+2

ZiZj

&

+4a(-l)i

[Pj,i+l (l-aa)-P&i-,

(l_Sil)]

(36)

Equation (36) makes soliton simulations within the PPP model computationally possible on computers like CYBER 855 where for a 79-unit chain one calculation of the gradient needs ca. 16 000 CPU when applying Su and Shrieffer’s procedure, while it needs ca. 4 CPU when using eqn. (36). Equation (33) also holds for the case of the formamide stack model system, giving the final result [491

with P ,=2

5 CirCk i=l

(36)

245

for the closed-shell case, where n* is the number of occupied orbitals. The formulae for the derivatives with respect to 0, and pmare similar. RESULTS AND DISCUSSION

Alternating tram-polyacetylene

As mentioned in the Introduction, one property of t-PA which has been explained by the concept of moving solitons [12] is the motional line narrowing of EPR lines in pristine samples. However, Chien and Babu 1291 reported the same line narrowing for copolymers of acetylene and carbon monoxide. In these copolymers the distance between two CO groups was the equivalent of almost 14 CH groups, being the soliton width within the SSH model. Chien and Babu concluded that if CO represents a barrier to the soliton movement then the observed line narrowing in both their copolymers and in t-PA cannot be explained by soliton movement. In order to study this problem we concluded in our program the option of putting a CO group at any desired site of the chain 1301. The oxygen atom is treated explicitly as a z center (1, = - 2.5 eV for oxygen and - 0.25 eV for the carbon bound to oxygen). In this case the starting geometry was Yi =

u, ,iiV, i -4

(39)

where No is the site of the CO group which already exhibits a kink structure at equilibrium. In this study we used a 30-unit chain with CO at site 15, a time step of z= 1.25 lo-l5 s, u ,,=O.l A implying a soliton halfwidth of ca. 3 lattice sites, and the SSH model. The time evolutions of the staggered coordinates are shown in Fig. 5 [ 301. Figure 5 (a) shows a neutral soliton approaching CO from the right and being trapped by the CO group. The same happens to a negatively charged soliton [Fig. 5 (b) 1. In this case the soliton can even change the C=O to a C-O- group and thus creates a fully conjugated chain. However, in our model, probably due to the rigidity of the CO group, this was not observed. Finally the positively charged soliton [Fig. 5(c) f is repelled by CO and rests after some oscillations midway between the CO group and the right chain end. Therefore, CO is, at least in the SSH model, a barrier for solitons. Thus the line narrowing in the copolymers casts some doubts on the soliton model for the same effect in pure PA. However, in the more realistic PPP model [ 201 as well as in MNDO calculations [ 601, the soliton halfwidth turns out to be reduced from ca. 7 (SSH) to ca. 3 for u,w 0.04 A. Therefore, there may still be enough space in the copolymers for solitons to move. This question can be only settled by numerical investigations which are in progress in our laboratory. Similar investigations using other impurities have also been performed [30 1.

Fig. 5. Time evoiution of the staggered-order parameter A = ur,/u,, (ug= 0.1A ) , ill~trati~ the dynamics of a neutral (A), negatively charged (B) and positively charged (C) soliton in a 30-unit chain approaching a CO group at site 15 in the SSH model (TX 1.25 10-'5s).(---)Initial end: kink structure.

247

An NH group which has a kink geometry just like CO repels neutral and neg-

atively charged solitons. In addition we have considered three substituents which are isoelectronic with CH, namely NH+, N, and O+. The substituent NH+ is a trap for negatively charged solitons, independent of the position of the group, relative to the kink. A neutral soliton is trapped if the position of NH+ is such that the resulting kink structure NH + has the same topology as the approaching soliton. If NH + is in a position where the resulting NH + structure is an anti-kink, the soliton can pass the group. The N substituent slows down So to some extent and S- is again trapped, while O+ reduces the velocity of So drastically and traps S- again. Therefore each impurity could trap or repel charged solitons, while neutral solitons though affected by the groups could pass most of the substituents isoelectronic to CH. Currently the effects of continuous parameter changes on solitons (which has been reported recently for the SSH case 1581) are being studied within the PPP Hamiltonian to find the parameter range in which solitons can still pass the impurity. However, it should be pointed out that a kink which does not move is still a soliton in PA due to the topological character of the phase boundary. Besides impurities, the effects of environment on solitons are also of interest. To perform a first model study we simulated the environment by adding two terms to the equations of motion [391. One term consists of random forces, R,, at each unit n, the other term is a dissipation term - ( - l)“+lyM( ayJ a$). For the randomly fluctuating forces R, a Gaussian distribution was assumed. W(R,)=~:(2n(R~))-~‘~exp

(-R2,/2(Rz))

(40)

where ( ) denotes the equilibrium average. In this study we were interested in the case of a chain in equilibrium with its environment. Therefore the fluctuation-~ssipation theorem was applied 1611 to ensure that, on average, no energy is gained from or lost to the environment (R; >=2hfkBTy/z

(41)

The same model parameters were used as described above and y was varied. Up to y= 10 l3 s-l no changes in the behavior of the kink were seen. At y= 1013 s-l the random forces were of the same order of magnitude as the internal lattice forces and in this case the kink was trapped in the middle of the chain. Applying eqns. (39) and (40) to a stationary kink in the middle of the chain lead to a random walk as shown in Fig. 6 for y= lOI s-’ and y= 1013s-‘. Presently the effects of random forces and dissipation are investigated within the PPP model removing also the requirement of equilibrium between the chain and environment. In addition we have studied doubly ionized chains within this model [40]. In the case of kink-free chains double ionization leads to the formation of a soliton-antisoliton pair each carrying one charge. Chains already containing

246

lo

N

”

3

10

,,,

20

3C

Fig. 6. Time evolution of the normalized staggered coordinate & = vr~/h (u,=O.l A) in a 31unit chain illustrating the dynamics of a neutral middle kink (initially) in unite of T= 1.25 lO_” s for y= 10” s-l (A) and y= lOI3 5-l (B) in the SSH model.

a neutral soliton tend to desorb an acetylene molecule at the chain end after ionization. Substitution of an NH group into the kink-free chain blocked the bipolaron formation since the electrons have been removed from the N lone pair. The NH group in a soliton-carrying chain could not prevent acetylene desorption. The relevance of the results described so far is limited by the fact that they

249

have been obtained within a simple Hiickel-type model. Therefore we have implemented the PPP model into our program as described above. In our first calculations [ 201 we varied the parameter yv for carbon describing the strength of the electron-electron interaction. The most striking result was that, upon inclusion of electron-electron interaction (4% 0.04 A), the soliton width decreased from ca. 7 sites (SSH) to ea. 3 sites, while its kinetic mass increased from ca. 6 m, (SSH) to ca. 15 m,. The width of three lattice sites, consistent with MNDO results, was obtained f511 and the usually adop~dv~ue of yrr= 11.08 using the Ohno p~~etri~tion eV [52 1. The Ma~ga-Nishimoto p~ametrization [62 ] lead to spin-density waves and to a preference of equidistant geometries. The soliton velocity, ca. 3.1 X lo* m s-i, was only slightly smaller than in the SSH model. The averaged spin densities of 0.06 and -0.02 on alternating sites measured by ENDOR could be reproduced only to a qualitative degree (0.10 and - 0.06). Figure 7 gives the time evolution of the staggered coordinates for the SSH and the PPP models (the first and last four CH units were kept fixed in these calculations). The neutral soliton levels are shifted to the band edges in PPP, consistent with the experimental finding of Weinberger et al, [21] that the mid-gap absorption cannot be due to neutral solitons. Since the soliton properties change considerably upon inclusion of electronelectron interactions, the studies reported above aE currently being repeated

-j,mh -

lb

site

210

2$

Fig. ‘7.The normalized staggered coordinates vi/% for a (CH), chain (first and last four sites

kept fixed) for different timee (r=l.% lo-l6 s) obtained with the PPP model using (-) Ohno’s p~met~~tion

- 11.08eV), and (---) the SSH model. (pm-

250

using the PPP model. In addition, the photodynamics currently under investigation in our laborato~.

of PA chains are also

Formamide stacks In our first time simulation on conformational solitons in formamide stacks [49 1, the initial excitation was considered as a shift of the first molecule 0.5 A away from the second one (N= 30). In Fig. 8 all three geometrical coordinates are shown as a function of time; each molecule is symbolized by a triangle. The backbone corresponds to the dashed line. The rotations are displayed such that the B axis is perpendicular to the paper plane for each molecule, and the 0 axis is identical to the z axis. The bottom line of each triangle gives 9. If q is positive the triangle stands on the bottom line, if (pis negative it stands on its top. The variables have been m~tiplied by six for better visualization. A narrow (ea. 6 sites) solitary wave can be seen to pass through the stack in 3.6 ps. Its velocity is 2.7 km s-l and the mean value of the kinetic energy is 8.4 meV, Therefore its kinetic mass is ca. 400 me (electron mass). Figure 9 shows the kinetic and total energies as functions of time. Apparently the total-energy fluctuations are tolerably small as compared with the kinetic energy. Figure 10 shows the collision of two such waves starting from opposite chain ends and Fig. 11 shows the corresponding energies. The two waves pass through each other without any visible perturbation. The picture for the one starting at the left is completely identical to Fig. 8 where only one soliton is present. The mean kinetic energy is exactly double that of a single soliton. Therefore

cs1

s

p”

s

8

Fig. 8. Time evolution of ali the geometrical variables of a 30-unit stack after initial excitation of the first unit by z1= - 0.5 A: propagation of a solitary wave (time step T= 0.01 ps).

251

t*

1 8 I

!

3 I

0.5

L

s I 1 f 2 1 ! 3 1 1 ! 1 1 1 4 1 ! 1 1 1 ’

1.5

1

2

2.5

’ 3

’ 1 1 1 ! ’ 1 ’ ’

3.5 tc p’j 1

’ 4

1 ’ ’ 1 ! 1 ’ ’ 1

4.5

! 5

’ ’ 1 ’ 1 ’ 1 ’

J

5.5

Fig. 9. Kinetic (solid line) and total energy (dotted line) as functions of time for the propagation of a solitary wave in a 30-unit formamide stack.

0

A7

d-”

s

s

Fig. 10. Time evolution of all the g~metrical variables of a 30-unit stack after an initial excitation of the first and last unit by f 0.5 A: collision of two solitary waves (time step ~0.01 ps).

no perturbation occurs during the collision of the waves which is a typical soliton property. Finally, we have studied the time evolution after a more complex excitation of the first three coordinates (z= -0.8 A, - 0.4 L&,- 0.2 A; e=p= -w, -6”, - 3 o ). For this case a shorter chain (I?= 15) was chosen.

252

22

20

10

16 214 .z iz 12

10

8

6

4

2

0

-2

Fig. 11. Kinetic energy (solid line) and total energy (dotted line) as functions of time for the collision of two solitary waves.

As shown in Fig. 12, the soliton passing twice along the chain is visible in this case also. The energies are shown in Fig. 13. However, in this case some vibrations are excited in addition to the solitary wave. Nevertheless, the existence of solitary waves in stacked systems is numerically proven by computer experiments.

Fig. 12. Time evolution of all the geometrical variables of a 15-unit stack after a complex excitation in all variables of the first two units (see text; time step t = 0.01 ps ) .

30

=_ 25 L z 20

15

10

5

Fig. 13. Kinetic energy (solid line) and total energy (dotted line) as functions of time for the complex initial excitation shown in Fig. 12.

254 CONCLUSION

As has been discussed above, both substitution and environment may well influence the properties of solitons in PA; they are even able to stop its movement. Since pristine PA is known to contain a rather high concentration of impurities (such as cis segments, cross links [ 631, etc.) the soliton mobility is expected to be rather limited, this is also indicated by ENDOR experiments [ 311. Electron-electron interactions have been shown to change both the width and kinetic mass of the solitons as well as the relative positions of the soliton levels rather drastically as compared to the case of the Hiickel-type SSH Hamiltonian. The reduction in soliton width is consistent with the MNDO results. Therefore, influences of impurities and environment must be recalculated at the PPP level. The kinetic mass, although doubled in the PPP model with ca. 15 m,, is still small enough to suspect that quantum dynamics may lead to considerable changes in the results [ 391. Although the large mass of the moving CH units and the results of previous studies [64] seem to oppose this notion, the question can be settled only by numerical investigations which are currently in progress. The existence of solitary waves in stacked systems has been proven numerically. The energy carried by them (8.4 meV) seems to be too small to interfere with DNA-protein interactions. However, a rather crude model system was chosen (formamide) and, due to the limited range of our potential fit, only low-energy excitations could be considered. The binding of a bulky carcinogen would introduce far larger geometry changes and a study with a simplified ansatz for the potential has shown that upon variation of the initial excitation energy the solitary wave always carries half the energy away [ 551. Therefore, in more realistic cases (nucleotide base stacks with higher energy excitations ) solitary waves carrying a much larger energy are expected to be found. The high kinetic mass of the solitary wave (400 m,) together with the large masses of the moving entities suggests that quantum dynamics should have a minimal effect in this case. The velocity of these conformational solitons is independent of the energy carried, while the kinetic mass increases with increasing energy [ 551. Inclusion of the dispersion energy into the model via London’s formula is in progress, as is the extension to real nucleotide base stacks. APPENDIX

To illustrate the performance of the three time-simulation methods outlined in the text let us consider a simple one-dimensional harmonic oscillator

w

(AlI

255

where k is the force constant, m the mass of the particle and IV’ the potential energy. Introduction of dimensionless coordinates gives

x=(~+x~)-l’zx ; T+)“‘t

(A21

where x0 and u. are the position and velocity, respectively at time to=O, and eqn. (Al) is transformed to d2X -= dT2

-X

+kx;)X2

; w’ +mu;

(A3)

The kinetic energy is given by =f (mu;+kx;) Introduction

(A4)

of dimensionless energy variables, W and K, gives

W=2(mu;+kx;)-1

W

K =2(mu;+kx&‘K

645)

E =2(mu;-kx$-1E’

;

E’=W’+K’

which leads to w=x2

046)

The exact solution of eqn. (A3) is then X =sin (T+arcsin

X0)

dX ~=COS

X0)

(T+arcsin

(A7)

and E=l With an initial velocity uo= 0 for arbitrary -l
648) x0, X0=1 holds. For uo#O,

where the minus sign applies for X,, < 0. The three simulation methods are then defined by [ V=g,

with z being the

time step and using eqn. (A3)] V(Zs-1) = V(Z) -X(Z)7

(AlO)

and for method A (SSH, this work) x~Z+l)~x~Z)+v~Z+~)~

(All)

for method B (difference ansatz ) X(Zf 1) =X(Z) + V(Z)r

fAl2)

and metbod C (direct integration with constant force during z) X(z+1)=X(Z)+V(z)r-~X(z)r2

(A13)

We have performed simulations over a full period 0 < T< 2n starting with a time step of z/2 which is reduced step by step by a factor of l/2.In Table Al, for X0 = 1 the maximum deviations of X, V, and E (AX, AV and dE) from the exact solutions are given for each simulation performed. Obviously method B, which is derived directly from the difference equations [eqn (4) ] performs worst of the three methods, method C works somewhat better. It should be noted that method C can also be derived from the trapezoidal rule. method A performs best and can be rewritten as a direct integral of the equations of motion, as for method C, but using 112 [V(l) + V(l) - 1) f as the constant velocity during time r. Thus, in the sequence of improving numerical accuracy, the three methods can be rewritten as

method B: X(Z+l)

=X(Z) +X(Z)z

(A14)

.

methodC:X(Z+l)=X(Z)+~[X(Z)+8(1+1)]~

(A15)

Table Al also shows that ail three methods converge to the correct solution for small values of r. Elowever, they are all one-step procedures and thus the accuracy improves only by a factor of ca. lf2 upon decreasing z by the same factor.

257

TABLEAl Maximum deviations dX and AE of the dimensionless coordinate X and dimensionless total energy E from the exact solution for a harmonic oscillator obtained with the simulation methods A, B and C (see text) for different dimensionless time-step sizes r (in unite of x, each calculation carried out in an interval of 2x) for the initial condition X0= 1 r(n)

112

u4 l/8 l/16 l/32 l/64 l/128 l/256 l/512 l/1024

AX” (X oscillates in + 1)

AE” (E is constant with value 1)

A

B

c

A

B

c

1.6142 0.5534 0.2310 0.1063 0.0510 0.0250 0.0124 0.0062 0.0031 0.0015

-8.7163 2.9363 2.0128 0.8260 0.3590 0.1665 0.0801 0.0393 0.0195 0.0097

3.3726 1.9425 0.8204 0.3586 0.1665 0.0801 0.0393 0.0195 0.0097 0.0048

3.6207 0.5820 0.2390 0.1083 0.0516 0.0251 0.0124 0.0062 0.0031 0.0015

143.5494 45.7045 8.9214 2.3552 0.8476 0.3608 0.1667 0.0802 0.0393 0.0195

39.9343 8.8531 2.3508 0.8472 0.3607 0.1667 0.0802 0.0393 0.0195 0.0097

“The numbers obtained for &=0.8 are very similar and, therefore, are not shown.

TABLEA Values of AX and AE (see Table Al) for the fourth-order Runge-Kutta (D) and the sixth order Mime (E) methods (see text)

t(n)

112

114 l/8 l/16 l/32 l/64

AX (X oscillates between k 1)

AE (E is constant with value 1)

D

E

D

E

-0.1046 -0.0063 -0.0004 -0.2x10-4 -0.1x10-~ -0.9x10-7

-0.1309 0.0028 0.0002 0.1x10-4 0.8~10-~ 0.8x10-'

-0.1923 -0.0065 -0.0002 0.1x10-" 0.9x10-6 0.6~10-~

0.3285 0.0045 0.0~3 0.1x10-4 0.8x10-" 0.6x10-'

Finally, we would like to show for this simple example the performance of the Runge-Kutta method (D ) which is correct up to fourth order, and of Milne’s predictor-corrector method (E) which is correct up to sixth order [65]. The formulae are given explicitly in ref. 65 and the two methods are currently being implemented in our programs. As Table A2 shows, these methods are definitely superior over one-step procedures and can be used to save a large amount of computer time because they allow larger time steps to be used The very large deviations from the exact solutions in the first two rows in Table Al are due

258

to the large time-step sizes of n/2 and ~14, respectively. These data are only shown for comparison with the corresponding values in Table A2 which are much smaller. ACKNOWLEDGEMENT

The financial support of the NATO Scientific Affairs Division (grant D.736/ 86) is gratefully acknowledged by two of us (J.L. and F.M. ).

10 11 12

13 14

15

16 17

18 19

M.A. Collins, Adv. Chem. Phys., 53 (1983) 225. J.A. Knzmhansl and J.R. Shrieffer, Phys. Rev. B, 11 (1975) 3535. H.J. Mikeska, J. Phys. C, 11 (1978) L29; K. Maki, J. Low Temp. Phys., 41 (1980) 327. R. Rqjaraman, Phys. Rep., 21(c) (1975) 227,313. D.W. Mc~~~in, in E. Clementi and R.H. Sarma (Eda.), Structure and Dynamics: Nucleic Acids and Proteins, Adenine Press, New York, 1983, p. 55. J. Scott-Russel, Proc. R. Sot. Edinburgh, (1844) 319. S. Aubry, J. Chem. Phys., 64 (1976) 3392; M.A. Collins, A. Blumen, J.F. Currie and J. Ross, Phys. Rev. B, 19 (1979) 3630. K.J. Wanted, J. Chem. Phys., 82 (1985) 5247; K.J. Wahlstrand and P.G. Woylnes, 6. Chem. Phys., 64 (1976) 3392. AS. Davydov, Phys. Script., 20 (1979) 387; SOV. Phys. Usp., 25 (1982) 898; AS. Davydov and H. Kislushka, Phys. Stat. Sol. B, (1973) 465; A.C. Scott, Phil. Trans. R. Sot. London, Ser. A, 315 (1985) 423; A.C. Scott, Phys. Script., 29 (1984) 284; Phys. Rev. A, 26 (1982) 578. J.A. K~rnh~~ and D.M.Alexander, in E. Clementi and R.H. Sarma (Eds.), Structure and Dynamics of Nucleic Acids and Proteins, Adenine Press, New York, 1983, p. 61. S.W. Englander, N.R. Kallenbach, A.J. Weeper, J.A. Krumhansl and S. Litwin, Proe. Nati. Aced. Sci. U.S.A., 77 (1980) 7222. W.P. Su, J.R. Schrieffer and Ad. Heeger, Phys. Rev. L&t., 42 ( 1979) 1698; W.P. Su, Solid State Commun., 35 (19801899; W.P. Su, J.R. Scbrieffer and A.J. Heeger, Phys. Rev. B, 22 (1980) 2099. J. Ladik and J. &%k, Int. J. Quantum Chem., 26 (1984) 955. J. Ladik, in Molecular Basis of Cancer, Part A, Alan Lise, Inc., New York, 1985, p. 343; J. Lad& in E. Clementi, G. Corongiu, MR. Sarma and R.H. Sarma (Eds.), Structure and Motion: Membranes, Nucleic Acids and Proteins, Adenine Press, New York, 1985, p. 553; J. Ladik, Quantum Theory of Polymers as Solids, Plenum Press, New York, 1988, p. 394. C.R. Fincher, Jr., C.-E. Chen, A.J. Heeger, A.G. Mac Diarmid and J.B. Hastings, Phys. Rev. Lett., 48 (1982) 100; S. Suhai, Phys. Rev., Ser. B, 27 (1983) 3506; J. Ladik, QuantumTheory of Polymers and Solids, Plenum Press, New York, 1988, p. 207. W.P. Su and J.R. Shrieffer, Proc. N&l. Acad. Sci. U.&A., 77 ( 1980) 5626. H. Shirakawa, T. It0 and S. Ikeda, Die Macromol. Chem., 179 (1978) 1565; B. Goldberg, H.R. Growe, P.R. Newman, Ad. Heeger and A.G. Mac Diarmid, J. Chem, Phys., 70 (1979) 1132. C.R. Fincher, Jr., M. Ozaki, M. Tanaka, D. Peebles, L. Lauchlan, A.J. Heeger and A.G. Mac Diarmid, Phys. Rev., Ser. B, 20 (1979) 1589. Z.G. Sooa and L.R. Ducasse, J. Chem. Phys., 78 (1983) 4092.

259

20 21 22

23 24

25 26 27 28 29 30 31

32 33 34

35 36 37 38 39 40 41

42 43 44 45 46

W. Forner, C.L. Wang, F. Martin0 and J. Ladik, Phys. Rev., Ser. B, 37 (1988) 4567. B.R. Weinberger, C.B. Roxio, S. Etemad, G.L. Baker and J. Orenstein, Phys. Rev. L&t., 53 (1984) 86. C.K. Chiang, C.R. Fincher Jr., Y.W. Park, A.J. Heeger, H. Shirakawa, E.J. Louis, S.C. Gau and A.G. Mac Diarmid, Phys. Rev. I.&t., 39 (1977) 1098; K.W. Park, M.A. Druy, C.K. Chiang, A.J. Heeger, A.G. Mac Diarmid, H. Shirakawa and S. Ikeda, J. Polym. Sci. Polym. Chem. Ed., 17 (1979) 195;C.K. Chiang, N. Park, A.J. Heeger, H. Shirakawa, E.J. Louis and A.G. Mac Diarmid, J. Chem. Phys., 69 (1978) 5098. Y.W. Park, A. Denenstein, C.K. Chiang, A.J. Heeger and A.G. Mac Diarmid, Solid State Commun., 29 (1979) 747. R.L. E~enbaumer, J.E. Frommer, J.-L. Bredas and R. SiIbey, in J. Ladik, J.-M. Andre and M. Seel (Eds.), Quanta Chemistry of Polymers; Solid State Aspects, Reidel, Dordrecht, 1984, p. 221. S. Kivelson, Phys. Rev. L&t., 46 (1981) 1344; Phys. Rev., Ser. B, 25 (1982) 3798. J.-L. B&as, B. Th&nans, J.-M. Andre, R.R. Chance and R. Siibey, Synth. Methods, 9 (1984) 265. M. Seel, Solid State Commun., 52 (1984) 467. J.-L. Bredas and G.B. Street, Act. Chem. Res., 18 (1985) 309. J.L.W. Chien, G.N. Babu and J.A. Hirsch, Nature, 314 (1985) 723. J.L.W. Chien and G.N. Babu, J. Chem. Phys., 82 (1985) 441; Macromolecules, 18 (1985) 622. W. Forner, M. See1 and J. Ladik, Solid State Commun., 57 (1986) 463; J. Chem. Phys., 84 (1986) 5910. H. ‘I’homann, L.R. Dalton, Y. Tomkiewicz, N.S. Shiren and T.C. Clarke, Phys. Rev. Lett., 50 (1983) 533; H. Thomson, H. Kim, A. Morobel-Sol, L-R. Dalton, M.T. Jones, B.H. Robinson, T.C. Clarke and Y. Tomkiewicz, Synth. Method, 9 (1984) 255. H. Thoman, J.E. Cline, B.M. Hofmann, H. Kim, A. Morobel-Sosa, B.H. Robinson and L.R. Dalton, J. Phys. Chem., 89 (1985) 1994. Z.G. Soos and S. Ramashesha, Phys. Rev. Lett., 51 (1983) 2374; A.J. Heeger and J.R. Shrieffer, Solid State Commun., 48 (1983) 207. M. Sasai and H. Fukutome, Synth. Methods, 9 (1984) 295. (a) J. Orenstein and G.L. Baker, Phys. Rev. Lett., 49 (1982) 1043. B.R. Weinberger, Phys. Rev. Lett., 50 (1983) 1693; G.B. Bianchet, CR. Fincher and A.J. Heeger, Phys. Rev. Lett., 51 (1983) 2132. G.B. Blanchet, CR. Fincher, T.C. Chungand A.J. Heeger, Phys. Rev. Lett., 50 (1983) 1938. (b) L. Rothberg, T.M. Jedju, S. Etemad and G.L. Baker, Phys. Rev., Ser. B, 36 (1987) 7529. A.R. Bishop, D.K. Campbell, P.S. Lormdahl, B. Horovitz and SR. Phiilpot, Phys. Rev. I.&t., 52 (1984) 671. C.L. Wang and F. Ma~ino, Phys. Rev., Ser. B, 34 (1986) 5540. W.P. Su, Phys. Rev., Ser. B, 34 (1986) 2988. S. Kivelson and Wei-Kong Wu, Phys. Rev., Ser. B, 34 (1986) 5423. A. Godzik, M. Seel, W. Farner and J. Ladik, Solid State Commun., 60 (1986) 609. C.-M. Liegener, W. Forner and J. Ladik, Solid State Commun., 61 (1987) 203. R. Pariser and R.G. Parr, J. Chem. Phys., 21 (1953) 704,466O. J.A. Pople, Trans. Faraday Sot., 49 (1953) 1375; J. Ladik, Acts Phys. Acad. Sci. Hung., 18 (1965) 185; J. Ladik, D.K. Rai and K. Appel, J. Mol. Spectrosc., 27 (1968) 72. M. Takahashi, J. Paldus and J. Cfiek, Int. J. Quantum Chem., 24 (1983) 707. E. Boyland, in Proceedings of the Israel Academy of Sciences, Symposium on Carcinogenesis, Jerusalem, 1968. J. Ladik, Int. J. Quantum Chem. QBS, 13 (1986) 307. J. Ladik, S. Suhai and M. Seel, Int. J. Quantum Chem. QBS, 5 (1978) 35. A.K. Bakhshi, J. Ladik, M. See1 and P. Otto, Chem. Phys., 108 (1986) 233.

47 K. Laki and J. Ladik, Int. J. Quantum Chem. QBS, 3 (1976) 51. 48 F. Belexnay, S. Suhai and J. Ladik, Int. J. Quantum Chem., 20 (1981) 683. 49 D. Hofmann, W. Fisrner and J. Ladik, Phys. Rev., Ser. A, 37 (1988) 4429. 60 J. Avery, J. Packer, J. Ladik and G. Bi&o, J. Mol. Spectrosc., 29 (1969) 194; M. Kert&z, S. Suhai and J. Ladik, Acta Phys. Acad. Sci. Hung., 36 (1974) 77. M. Kert&z, Chem. Phys., 44 (1979) 349; A. Warshel and M. Karplus, J. Am. Chem. Sec., 94 (1972) 5612; A.C. Lasaga and M. Karplus, Phys. Rev., Ser. A, 16 (1977) 807; J. Chem. Phys., 71 (1979) 1218. AC. Lasaga, R.J. Aerni and M. Karplus, J. Chem. Phys., 73 (1980) 5230. J. Paldus and E. Chin, Int. J, Quantum Chem., 24 (1983) 373; J. PaIdus, E, Chin and M.G. Gray, Int. J. Quantum Chem., 24 (1983) 395; J.-L. B&as, M. Dory and J.-M. Andre, J. Chem. Phys., 83 (1985) 5242. 51 K. Ohno, Theor. Chim. A&e, 2 (1964) 219. 52 J. Iadik, Quantenchemie, F. Enke Verlag, Stuttgart, 1973. 53 S. Kivelson and D.E. Heim, Phys. Rev., Ser. B, 26 (1982) 4278. 54 S. Arnott, SD. Dover and A.J. Wonacott, Acta Crystallogr., Sect. B, 25 (1969) 2192, 55 W. F&net, Phys. Rev., Ser. A, 38 (1988) 939. 56 J.A. Pople and W. Hehre, Q. C. P, E., (1974) 236. 57 C.L. Wang, ZB. Su and F. Martino, Phys. Rev., Ser. B, 33 (1986) 1512. F. Guinea, Phys. Rev., Ser. B, 30 (1984) 1884. 58 S.R. Phillpot, D. Baeriswyl, A.R. Bishop and P.S. Lomdahl, Phys. Rev., Ser. B, 35 (1987) 7633. 59 W. Farner, Solid State Commun., 63 (1987) 941. 69 D.S. Boudreaux, R.R. Chance, J.-L. Br&les and R. Silbey, Phye. Rev., Ser. B, 28 ( 1983 ) 6927. 61 M.A. Collins, Adv. Chem. Phys., 53 (1983) 225. M.A. Cohins, A. Blumen, 59. Currie and J. Roes, Phys. I&v., Ser. B, 19 (1979) 3630; R, Kubo, Rep. Prog. Phys., 29 (1966) 255. D. Forster, Hydrodynamic Fluctuations, Broken Symmetry and Correlation Functions, Benjamin, Reading, 1975. 62 N. Mataga and K. Nishimoto, Z. Phys. Chem., 13 (1957) 140. 63 H.W. Gibson, R.J. Weagley, R.A. Masher, S. Kaplan, W&I. Prest, Jr. and J.A. Epstein, Phys. Rev., Ser. B, 31 (1985) 2338. 64 J.E. Hirsch and M. Grabowski, Phys. Rev. Lett., 52 (1984) 1713. D.K. Campbell, T.E. De Grand and S. Mazum~r, Phys. Rev. Lett., 52 (1984) 1717. Functions, Dover, New York, 65 A. Abramowitx and P.A. Stegung, Handbook of Mathe~tical 1972, p. 897.