Cubic response functions in the random phase approximation

11 August 1995

CHEMICAL PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 242 (1995) 7-16

Cubic response functions in the random phase approximation o

Patrick Norman, Dan Jonsson, Olav Vahtras, Hans Agren Institute of Physics and Measurement Technology, Linkfping University, S-581 83 LinkSping, Sweden Received 23 March 1995; in final form 19 June 1995

Abstract

Cubic response functions in the random phase approximation have been derived and their use for computations of static and dynamic second hyperpolarizabilities is demonstrated. The performance of implemented computer strategies in terms of direct one-index transformations and of direct atomic orbital constructions of key elements is discussed. A demonstration is given by calculations of various components of the second hyperpolarizability tensor of two linked thiophene rings.

1. Introduction Experimental progress in determining non-linear molecular spectra and properties has challenged a concordant development in theory and computational techniques. Some of this challenge has indeed been met by recent advances in simulations of non-linear phenomena using microscopic quantum models. These aimed to carry out accurate calculations for small species and to obtain structure-to-property relations for somewhat larger species. Response theory methods constitute one branch of this development, and have now been diversified for calculations of different non-linear properties; electric or magnetic, time-independent or time-dependent non-linear properties [1-7]. Along with the increase in the number of applications of these methods, it has become more important to reach for new combination properties that are addressed in modern spectroscopic techniques, and to explore efficient algorithms for solutions of the response equations. In the present work we derive cubic response functions with applications to time-dependent and time-independent electric properties. We demonstrate an implementation of cubic response theory for self-consistent field wavefunctions, which, in the response terminology, coins the method as the cubic random phase approximation (CRPA). The motivation for doing this lies in the great number of new non-linear properties that then come within reach, and in that the second hyperpolarizability, which relates to the cubic response function, constitutes a key quantity in the search for non-linear optical materials. Most calculations of the second hyperpolarizability 3' have referred to finite field calculations of the polarizability o~ and the first hyperpolarizability/3, or even to differentiation of the total energy E with respect to the field. However, direct non-approximate calculations of 3' have both formal and computational advantages. Earlier ab initio calculations for second hyperpolarizabilities have utilized one-determinant self-consistent field wavefunctions and the coupled perturbed Hartree-Fock approach (CPHF) for time-independent [8,9] and time-dependent [10-12] fields. The wavefunctions utilized in correlated treatments [13,14] have been both of perturbational M011er-Plesset [15,16] and of coupled cluster [17,18] types. 0009-2614/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved 6-4

SSDI 0 0 0 9 - 2 6 1 4 ( 9 5 ) 0 0 7 1

8

P. Norman et al./Chemical Physics Letters 242 (1995) 7-16

We start out from the response formalism [3] which in the linear [19] and quadratic [20] cases has demonstrated some remarkable features for calculations of extended species [7]. We derive the working formulas in which RPA equations are solved iteratively in terms of linear transformations, having the computation of the RPA matrices driven in terms of atomic orbital integrals, and with direct computations of the so-called one-index transformations. We demonstrate the present formulation by computing various components of the second hyperpolarizability tensor of bithiophene.

2. Method

The present method is formulated in terms of general Fock matrices of linearly transformed vectors as employed in current multiconfigurational SCF (MCSCF) property calculations. The approach of this general technique does not assume canonical Hartree-Fock orbitals, and allows for solving the sets of linear equations and the eigenvalue equations by means of the same atomic orbital based direct techniques. The method is double-direct because it is driven directly by the atomic orbital integrals, and, because it uses iterative techniques based on direct linear transformations for solving the RPA equations. The latter feature involves linear transformations of Hessian-like matrices on vectors (as in direct CI) that are carried out without explicit construction of the matrices. Cubic response calculations also contain a term where the fourth energy derivative matrix is multiplied on three trial vectors. We calculate this term as a triple linear transformation. This transformation is also carried out in terms of integrals in the atomic orbital basis. The linear transformations are formulated in terms of gradient vectors of an operator with one-index transformed integrals. This is the transformation that the integrals of an operator undergo when a commutator is formed with the operator and a generator of orbital rotations, and is the basis for direct MCSCF [21,22]. The one-index transformations are also key operations for linear [23], quadratic [5] and (as shown here) cubic response theory calculations with the modification in the latter case that the gradient vectors are evaluated over triply one-index transformed integrals. The most straightforward implementation is to one-index transform the integrals and store them on file. However, both storage and time consumption can be greatly reduced if the one-index transformations are carried out in a direct fashion without pre-transforming or storing the integrals [24]. These options are included in the present method.

2.1. Third-order response function We have started out from general response theory and derived expressions for the third-order response function for closed shell Hartree-Fock systems. We have used the notation of Olsen and Jcrgensen [25], which also has been used in their analogous derivations for the first and second order response functions [23,5]. The cubic response function reads ( ( A ; B, C, D))o~,,o:o~ = -'~-NjA(0)1 -'~ 0) 2 "q'-0)3)T]4}mN:(0),)N,C(0)z)NmD(0)3) --/vjA(0)I "~ 0")2 "4- 0)3) [Tj[k3}(0)1 ' 0)2 "~ 0)3)NB(0)I)NtCD(0)2, 0)3)

0)1 + 0)3)N[(0)2)NPO(0),, 'o3)+

0)1 +

+NjA( 0)1 4- 0)2 -'b 0)3)[ B}2]N[°( 0)2, 0)3) "b C~2]NBD( 0)1, 0)3) q- D~INfC( 0)1, 0)2)] [31 C (0)2)ND(0)3) + Cj(kl)N; [31 n (0)l)Nt D(0)3) +D)(kI)N; [3] B (0)I)N/ C (0)2)] --NjA(0)I + 0)2 + 0)3)[Bj(kl)Ni +a~2)k,[ Njn( O)I)NCD( 0)2, 0)3) + NjC( O)2)N:°( 0)1, 0)3) + Nj°( 0)3)N:C( 0)1, 0)2)] _a[31 IqB[, 0),)NC( 0)2) ND(0)3), "'{j*t)~'j

(1)

P. Normanet al./ ChemicalPhysicsLetters242 (1995)7-16 NjX(oja)=(E[2]-waS[2l)~lX[kq , NjXY( C-Oa, (-Ob) XYj[1]( ¢..Oa,

= ( E [2]

X~{A,B,C,D},

(2)

-- ( oj a -~- ¢..Ob)S[2])jkl Xy[l]( O)a, COb) ,

OJb) = Tj[~)(¢.Oa,

Wb)NX(oga)NtV(Wb)_ v[Zl~Xt, jk "k ~.W,~) -- X)ktz]Arky(Wb),

(3)

XY~ { BC, BD, CD}, where we have introduced

o,1)

(4)

=

Tj[~ m ~- ( E~2lm, - (.OlSl41,m) - (.02S~;~km) -

(5)

o)3S}4m],k/)).

Indices enclosed in parentheses are to be permuted symmetrically. We assume below that A, B, C, and D refer to one-electron operators. If A, B, C, and D are dipole operators the second hyperpolarizability becomes equal to the value of the third-order response function with negative sign; YABCO(--(O~] + OJ2 + W3); Wl, W2, W3) = --((A; B, C, D)),ol,o2o,3. The evaluation of the third-order response function can be separated into two steps; first solving a set of seven linear equations (2) and (3) to obtain the corresponding seven response vectors to be followed by the matrix multiplication (1) for the response function value. Compared to the second-order response function we now encounter Hessian, overlap and property matrices with one additional index; E [41, S t4], and A [3]. For implementation purposes we need a way of multiplying these large matrices on vectors without explicitly forming them. If we, in analogy with previous implementation for lower order matrices, write /N=(~,'K),

/=1,2,3

(6)

the matrix multiplications can be expressed in terms of doubly and triply one-index transformed operators

I { (OI[q~, A(1K, 2)] [0) ) A~z]t]Nk2N'=-6 (0l[qj, A(1K, 2)]10> ' 1 {(0][qJ' H ( ' K ' 2 ' 3 ) ]

Ej[4] 1 M 2 ~ ? 3 M

(7) 10)

klm'*k,*l~*m~---~ (0 i[q~, H(, K, 2 , 3 ) ] 10) , K.[4] , M 2 M 3 M

°Y"" " " " t ' "

l

: 6

3 ~ 2 t 1 ¢ 1 t (O[[qj'[(3kq+k+Kkqk),[(2tq~+Klqt),[(Kmqm+Kmqm)]]]]l 0) 3 ? 3 t 2 ," 1 "I" 1 ' (0l[q~, [( Kkq, +Kkq,) , [(2tq ~ + Ktq,) , [(K,q, + Kmqm)ll]] 10>

(8)

,

(9)

where the transformed operators are given by A('K, 2 ) =

1 t t 2 [(K,q, + 1K,q,), [(Ktq l~- +2K~qt) , A l l ,

H(1K,~,3)=[(Kkqk+lK,kq k), [(Ktqt 2 ++ZK~qt), [(3,,q+m +3K'q,.), " ] ] ] . ' +

(10) (11)

2.2. Double direct construction of Fock matrices Since (8) is of the form of a gradient vector containing a triply one-index transformed Hamiltonian it can be expressed in terms of a Fock matrix with triply one-index transformed integrals. Forming this Hamiltonian in a non-direct fashion involves tedious integral transformations. Also the need to store integrals would restrict our

10

P. Norman et al. / Chemical Physics Letters 242 (1995) 7-16

applications to small systems in comparison to what one deals with in direct SCF methods. We therefore wish to construct the necessary Fock matrices over atomic integrals as well as using a direct linear transformation for the triply one-index transformed Fock matrix. Denoting inactive MOs with j, general MOs with p, q, r, s, and AOs with greek letters we have

Frs = hrs + C~rsjj,

(12)

.ZP~sjj = 2( rs I jj) - ( rjl js).

(13)

The one-index transformed one- and two-electron integrals with ~ take the form Frs =-hrs "~ ~ r s j j '

(14)

-hrs = lKrphps -1Kpshrp ,

(15)

"~rspj = 1Krq~qspj -- 1Kqs~rqpj + 1Kpq~rsqj -- 1Kqj~rspq"

(16)

Using (14)-(16) repeatedly we obtain after expanding _~ in atomic orbitals F(1K, 2K, 3K) = F 123 + [3K, F 12] + [2K, F 13 + [ 3 , FI]] +[11¢, F 23 + [ 3 , F 2] + [2K, F 3 + [3K, F]]], where the intermediate Fock matrices and the density matrices needed for their construction are given F123 = /')123CW u~13~rs,~/3. D123 1 2 3

(17)

by (18)

1

2

3

1

2

3

a~ = Kjp Kpq KqtCottCflj -- Kjp Kpq KtjCaqCflt -- Kjp Kqj KptCatCflq .{_1 2 3 1 2 3 ,1 2 3 Kjp KqjKtqCc~pCflt -- Kpj KjqKqtCc~tCflp "l- Kpj KjqKtpCaqC~t

1Kpj 2Kqp 3KjtCatCflq -- 1Kpj 2Kqp 3KtqCajCflt ~ mn f rmn= O2e .~,,,

(19) (20)

D~"~" = . KjpKpqCaqCflj-. . . . KjpKqjCapCf3 q. -

. . KpjKjqCaqCf3 p + KpjKqpCajC#q,

(21)

F " = n~,"~,,-~so~¢,

(22)

D ~ = mrjp C.p ct~j - mrm c.j CBp.

(23)

If we permute IK, 2K, and 3 in (17) we obtain an expression for E}~l,m) times three response vectors,

(k,.,) "k " t ",,, = - 2

-/~/s]' /

/: = p(IK, 2K, 3 K ) F ( 1 K , 2K, 3K)

= ['K, [2K, [3, F] + 3F31 + [3, [2K, F] + 3F 2] + 3F (23)1 +[2K, [1K, [3, F] + 3F 3] + [3,[1K, F] + 3F 1] + 3F (13)] + [3, [2K, ['K, F] + 3F'] + [1K, [21~, F] + 3F 2] + 3F (12)] + F (123,.

(25)

The expressions for the intermediate Fock matrices in (25) have the same structure as the ones in direct SCF.

11

P. Norman et al. /Chemical Physics Letters 242 (1995) 7-16

2.3. Computational effort For the present implementation the most CPU demanding task is to solve the set of linear equations (2) and (3) for the response vectors. To see how much more complex the evaluation of the cubic response function is in comparison to the quadratic response function we recall what the corresponding equations for the latter look like [25]. ( ( a ; B, C ) ),ol,oz = --Nit (0)1-{- (.02) (Tj[~2(o.) 1, t ° e ) g f f ( t ° l ) g F ( w 2 ) -J-u[2l~ct' " '-k'yk, t°a) + Cyn[ZlNAB (to1) ) t,l oJl)A<;,)N ~c (o~:),

~x(oJo)=(Etel-~ooStel);-kIXt~l, E [31 1N 2N = - 2 j(kl)

k

I

X~{A,B,C},

,

(26) (27) (28)

-Fiis

T= P('K, 2)T(1K, 2) = ['K, [ 2 , F ] + 2 F 2] + [ 2 , [ ~ , F ] + 2 F 1] + F '~:).

(29)

In our example with the thiophene molecule 95% of the total CPU time was needed to solve the set of linear equations for the response vectors. Since this part is so completely dominating, the scaling factor between a hyperpolarizability and a second hyperpolarizability calculation will be determined by the ratio 7/3, i.e. the ratio of the respective number of response vectors. The remaining 5% of the total time is dominated by the matrix multiplication E [31 times two vectors and E [4] times three vectors. The scaling factor between second and third order theory for this part only will in the direct implementation be determined by the number of Fock matrix constructions. Since the E [4] term involves seven Fock matrices, see (25), and the I::[3] term three, see (29), we get a scaling factor of ~ for the complete matrix multiplication.

3. Application As a preliminary demonstration of CRPA we compute the various second hyperpolarizability tensor elements for thiophene and bithiophene, and investigate their basis set and frequency dependencies. Thiophenes represent a category of polymeric compounds with non-degenerate ground states that have attracted interest because of their special doping induced conductivity properties, their luminescence and third order non-linear optical properties [26-28]. As for polyacetylene the non-linear properties exhibit strong length and frequency dependencies [28].

3.1. Basis set Hyperpolarizability calculations are known to be sensitive to the choice of basis set. A good description of the tail region of the electron distribution is important, something that is accomplished by including polarizing and diffuse functions in the basis set. We used the thiophene molecule for the purpose of finding an adequate basis set. We started with a 4-31G basis set [29,30] and added to that polarizing and diffuse functions. The exponents of the diffuse functions were 0.05 and 0.01, denoted in Table 1 as d2 and dl, respectively. Thus s-dl means we have added a diffuse function with exponent 0.05 to the block of s-type functions. We explored five different basis sets, B1-B5. We see in Table 1 that adding polarizing functions on the heavy atoms diminishes the value of y, an observation that has also been made for the hyperpolarizability/3 for other compounds [20].

12

P. Norman et al. /Chemical Physics Letters 242 (1995) 7-16

Table 1 Basis sets for second hyperpolarizability ( X 103 au) calculation of thiophene sulfur s-dl p-dl p-d2 d-dl d-d2 d-pol

B1

B2

B3

x

×

× x

×

x

x

B4

B5

X x x x

x x x x ×

x x x x

x x x x x

20.7 15.6 5.1

17.9 12.9 5.0

×

carbon s-dl p-dl p-d2 d-d1 d-d2 d-pol

x

x

× x

x

x

x

x

hydrogen s-dl yxxxx(-3to0; tOo, tOo, tOo) yxxxx(O;o, o, o) /i tOo= 0.042823 au

19.3 14.6 4.7

17.5 12.8 4.7

19.6 14.7 4.9

However, the dispersion is merely changed by an additive constant. The basis set B1 is well suited for the bithiophene calculation. 3.2. F r e q u e n c y d e p e n d e n c e

In Figs. 1 and 2 we display the frequency dependence of different components of the third harmonic generation tensor for bithiophene. Only components that contribute to the averaged value given in Eq. (30) are displayed. The coordinate system is given in Fig. 3. The choice of axes is made so that x is along the polymeric in-plane long-axis, y along the polymeric in-plane short-axis, and z as the out-of-plane axis. As expected the strongest hyperpolarizability component is the long-axis Yxxxx component with values one order of magnitude larger than the others. It is also the most frequency-dependent component and becomes therefore even more prominent for larger frequencies. An orientationally averaged scalar value for T, see Fig. 1, is obtained as 1 "Y

< ~ > = 5 ( xxxx -[- "Yyyyy "~- ~/zzzz -[- ~lxxyy -[- ")lyyxx -]- ~xxzz -[- "Yzzxx "[- ~/yyzz + "~zzyy ) ,

(30)

which can be compared to experimental results. We can also see from Fig. 2 that the so-called Kleinmann symmetry, i.e. ~/xxyy = "Yyyxx,

"Yxxzz

=

"Yz. . . .

•zzyy = ~/yyzz

(31)

is a reasonable approximation for frequencies up to about 0.02 au for the bithiophene molecule. The equalities hold trivially for zero frequency. In Fig. 4 we have computed the dispersion of the xaxx component of the Kerr, the degenerate four-wave mixing (DFWM), the second harmonic generation (SHG) and the third harmonic generation (THG) properties of bithiophene. We obtain a strong dispersion of the THG well before noticeable

13

P. Norman et al. / Chemical Physics Letters 242 (1995) 7-16

104

10s

"2'.

O

~

a)

a)

0

.~IV

e-

~ 0. 3 frequency (a.u.)

b) o -----~h)

OOl frequency (a.u.)

Fig. I. Frequency dependence of (a) the long-axis component Yxxxx(- 3 to; to, w, to) and (b) the orientationally averaged value (3') of the second hyperpolarizability for bithiophene. Fig. 2. Frequency dependence of the in-plane short-axis and out-of-plane components (a) 3"yyyy(--3w; oJ, w, to); (b) Yx,~:(-3to; to, to, w); (c)3'z~zz(-3w; to, to, oJ); (d)yzzxx(-3o~; to, to, to); (e)3,ryzz(-3to; to, to, to); ( f ) y : z r y ( - 3 t o ; to, to, to); (g)yyyxx(-3to; to, to, to) and (h) 3'xxyy(- 3to; to, to, to) of the second hyperpolarizability for bithiophene.

Fig. 3. Symmetry axes for bithiophene; the central C - C axis makes an angle of 17.7 ° with the x axis. The x axis is chosen as the in-plane long-axis of the poly-thiophene polymer.

14

P. Norman et al. / Chemical Physics Letters 242 (1995) 7-16

x 105 i

t

i

i

2.5

.=

E~

O

1.5

0

0.01

0.02 0.03 frequency (a.u.)

0.04

Fig. 4. Frequency dependence of(a) 3Zxxxx(-3to; to, to, to); (b) 7xxxx(-2to; to, to, 0); (c) 3(~xxx(- to; to, to, - to) and (d) 3'xxxx(- to; to, 0, 0) for bithiophene. dispersion of any of the other properties, which is linked to the fact that singular behaviour is obtained at one third of the energy of the H O M O - L U M O optical gap.

4. Discussion In the present work we have taken one further step in response theory with respect to order, namely from linear and quadratic response functions and considered cubic response functions with applications on time-dependent and time-independent electric properties. One motivation to go one order higher in response theory, from the quadratic to the cubic case, follows from the fact that a great number of new nonlinear phenomena then become within reach. Counting only a few electric properties, cubic response theory covers phenomena such as optical rectification, third harmonic generation, dc electric field induced second harmonic generation, intensity dependent refractive index, second hyperpolarizabilities, and several effects that can be obtained by combining different hyperpolarizability tensors and time-independent or time-dependent fields. A second reason to introduce cubic response functions is that the experimental search for new non-linear materials seems to focus on organic compounds with exceptionally large second-order susceptibilities. These

P. Norman et al. / Chemical Physics Letters 242 (1995) 7-16

15

refer to the microscopic second-order hyperpolarizabilities y of the constituent molecules. The second hyperpolarizability thus plays a key role in understanding the non-linear optical response of molecules. For instance materials with strong third harmonic generations can be used for doubling and tripling of laser frequencies introducing laser technology into the blue and near-UV wavelength regions. Other examples are given by the electro-optical Kerr effects, which describe the action of an external field on the polarization directions of light, and by strong charge transfer complexes and organic molecular crystals with potential use in optical switching. All these capabilities are directly connected to the molecular hyperpolarizability tensors. We demonstrated an implementation of cubic response theory for self-consistent field wave functions. The cost of evaluating a single non-linear molecular property scales with that of optimizing the wavefunction, and with the construction of the corresponding linear [19] and quadratic [20] response functions. These scaling factors have been derived here. The bottleneck of such calculations is the same as for direct SCF wavefunctions: the calculation, storage, retrieval, and possibly transformation of two-electron integrals. In order to alleviate this we have implemented some key features, namely direct linear transformation in the solution of the RPA equations, direct atomic orbital construction of matrices and direct so-called one-index transformations. Thus the implementation is direct both in the sense that RPA matrices and other computational elements are calculated in terms of AO integrals evaluated on the fly, as in direct SCF, and it is direct in the particular way of solving linear equations and eigenvalue equations using direct linear matrix transformations, which avoids the explicit construction of the coefficient matrices. Applications can therefore be performed on species of the same size as in direct SCF. As for direct SCF, the ultimate applicability of the method in terms of the size of the system depends on the efficiency on the Fock matrix construction over the integrals, and the implementation of density screening. Parallelization over integral distributions can be achieved as in direct SCF, but also, at a more coarse level, over equation solvers. The latter is the time-consuming part in the present CRPA method. With some extra coding of residues, three-photon excitation spectra and two-photon transition amplitudes between excited states, can be obtained. Likewise, one-open-shell wavefunctions would be fairly straightforward to implement along the lines given in the present Letter.

Acknowedgement This work was supported by CRAY Research Inc.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

J. Oddershede and P. JCrgensen, J. Chem. Phys. 66 (1977) 1541. E. Dalgaard and H.J. Monkhorst, Phys. Rev. A 28 (1983) 1217. J. Olsen and P. J0rgensen, J. Chem. Phys. 82 (1985) 3235. J. Geertsen, J. Oddershede and G.E. Scuseria, J. Chem. Phys. 87 (1987) 2138. H. Hettema, H.J.Aa. Jensen, P. Jorgensen and J. Olsen, J. Chem. Phys. 97 (1992) 1174. O. Vahtras, H. ,~gren, P. J0rgensen, H.J.Aa. Jensen, T. Helgaker and J. Olsen, J. Chem. Phys. 97 (1992) 9178. Y. Luo, H..~gren, P. Jcrgensen and K.V. Mikkelsen, Advan. Quantum Chem., in press. P. Lazzeretti and R. Zanasi, J. Chem. Phys. 74 (1981) 5216. J.K. Hurley, H. Lincshitz and A. Treinin, J. Phys. Chem. 92 (1988) 5151. R. Klingbeil, V.G. Kaveeshwar and R.P. Hurst, Phys. Rev. A 1760 (1981) 5216. H. Sekino and R.J. Bartlett, J. Chem. Phys. 85 (1986) 976. S.P. Karna and M. Dupuis, J. Comput. Chem. 12 (1991) 487. H. Sekino and R.J. Bartlett, J. Chem. Phys. 94 (1991) 3665. J.E. Rice, G.E. Scuseria, T.J. Lee, P.R. Taylor and J. Alml/Sf, Chem. Phys. Letters 191 (1992) 23. J.E. Rice and N. Handy, Intern. J. Quantum Chem. 43 (1992) 91.

16 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

P. Norman et al. / Chemical Physics Letters 242 (1995) 7-16

A. Willets, J.E. Rice, D.M. Burland and D.P. Shelton, J. Chem. Phys. 97 (1992) 7590. H. Sekino and R.J. Bartlett, J. Chem. Phys. 98 (1993) 3022. H. Sekino and R.J. Bartlett, Chem. Phys. Letters 234 (1995) 87. H. Koch, H. Agren, P. Jcrgensen, T. Helgaker and H.J.Aa. Jensen, Chem. Phys. 172 (1993) 13. H. Agren, H. Koch, O. Vahtras, P. Jcrgensen and T. Helgakar, J. Chem. Phys. 98 (1993) 6417. J. Olsen, D.L. Yeager and P. J~rgensen, Advan. Chem. Phys. 54 (1983) 1. H.J.Aa. Jensen and H..~gren, Chem. Phys. Letters 110 (1984) 140. P. J~rgensen, J. Olsen and H.J.Aa. Jensen, J. Chem. Phys. 74 (1988) 265. O. Vahtras, H. Agren and H.J.Aa. Jensen, Chem. Phys. Letters 15 (1994) 573. J. Olsen and P. J~rgensen, Time Dependent Response Theory with Applications to Self Consistent Field (SCF) and Multiconfiguration Self Consistent Field (MCSCF) Wave Functions. I.F.A. PRINT, Aarhus Universitet (1994). J.L. Bredas and R. Silbey Conjugated polymers (Kluwer Acedemic Publishers, Dordrecht, 1991). C. Taliani, R. Danieli, R. Lazzaroni, N. Periasamy, G. Ruani and R. Zamboni, Synth. Metals 55 (1993) 4714. M. Zhao, P. Bhanu and P.N. Prasad, J. Chem. Phys. 89 (1988) 5535. W.J. Hehre and W.A. Lathan, J. Chem. Phys. 56 (1972) 5255. R. Ditchfield, W.J. Hehre and J.A. Pople, J. Chem. Phys. 54 (1971) 724.