On multidimensional avoided crossings of potential energy surfaces

On multidimensional avoided crossings of potential energy surfaces

151. number 3 Volume CHEMICAL ON MULTIDIMENSIONAL PHYSICS LETTERS 14October 1988 AVOIDED CROSSINGS OF POTENTIAL ENERGY SURFACES E. OHRENDORF, ...

739KB Sizes 6 Downloads 61 Views

151. number 3

Volume

CHEMICAL

ON MULTIDIMENSIONAL

PHYSICS LETTERS

14October

1988

AVOIDED CROSSINGS OF POTENTIAL ENERGY SURFACES

E. OHRENDORF, L.S. CEDERBAUM AND H. KtiPPEL Theoreiische Chemie, institut ftir Physikalrsche Chemle. Universrfdt Heldelbers,

D6YVO Heidelberg

Federal Republic

of tiermany

Received 15 April 1988; in final form 27 August 1988

In polyatomic molecules, the Wigner-van Neumann non-crossmg rule is generally relaxed because of the cxistcncc of more than a single nuclear coordinate. Potential energy surfaces may or may not intersect. There is an interesting class of states where an intersection cannot occur if only the dominant electronic configurations are considered. A general way to analyze the influence of all other electronic configurations is presented. A multidimensional avoided crossing in the ethylene dication is discussed as an illustrative example.

1. Introduction The conditions for intersections of potential energy surfaces of polyatomic molecules as well as the topology of the intersecting surfaces are of general interest [ l-41. For g surfaces to intersect at a single nuclear configuration, i.e. for a g-fold degeneracy to occur, a number of conditions must be satisfied. To simplify the discussion it is generally assumed that the electronic energies of the system are given by the eigenvalues of a real Hermitean matrix H = {FI,} of finite dimension, say n. As many as f(/=?l)(g+2)

(1)

conditions have to be fulfilled to obtain a g-fold degeneracy (g< n ). This number of conditions results from the general formula derived by von Neumann and Wigner [ 5 1. In the special case g= n the conditions can be given explicitly (see, for instance, ref. [6] ). They comprise n - 1 conditions for the diagonal elements of H, H, , =H22 = ...=H.,

,

and 1n (n - 1) ones for the non-diagonal UC,= 0,

for

if

j

(2a) elements, (2b)

Thus, all coupling elements must vanish and the uncoupled states must be degenerate at a single nuclear configuration. For two states of the same symmetry to intersect it is sufficient to fulfil only two inde0 009-2614/88/$ ( North-Holland

03.50 0 Elsevier Science Publishers Physics Publishing Division )

pendent conditions (see eq. (2)). This requires the existence of at least two independent nuclear coordinates (accidental degeneracies are rare and need not to bc considered). In a diatomic molecule the interatomic distance is the only available coordinate so that the non-crossing rule of von Neumann and Wigner [ 51 follows immediately. For polyatomic molecules additional degrees of freedom are present and potential energy surfaces belonging to states of the same symmetry can intersect. This has been shown, for example, to occur for the lower electronic states of CH,C [ 61. Here, we point out that there is an interesting class of systems where two potential energy surfaces are likely not to intersect although plenty of nuclear degrees of freedom exist, in principle, to allow for an intersection to occur. In this class of states, the leading electronic configurations (e.g. Slater determinants built of Hartree-Fock or another kind of suitable orbitals) are doubly excited with respect to each other, where the excitation involves only two spatial orbitals, say a and b. We call this excitation “pair excitation”. The matrix element H,, describing the interaction between the two electronic configurations is the common exchange integral Knb. It can be shown that such an integral Kab cannot vanish for finite internuclear distances [ 71. Consequently, considering only the two leading electronic configurations, the resulting potential energy surfaces cannot intersect. For a polyatomic system they exhibit B.V.

273

Volume 15 I, number

3

CHEMICAL

PHYSICS LETTERS

a nzultidimensional avoided crossing provided that the diagonal elements H,, and H,, cross as a function of the nuclear coordinates. In reality the exact electronic states are not fully determined from two electronic configurations only. Rather, they contain an admixture of many additional configurations, i.e. they are solutions of an electronic Hamiltonian with a much larger dimension than the number g of intersecting or nearly intersecting potential energy surfaces. Therefore, the simple conditions (2) are not directly applicable. To find out whether the multidimensional avoided crossing still exists under the influence of the additional n-g electronic configurations, one may, of course, just compute the surfaces by diagonalizing H for many nuclear distances and decide by inspection. This approach provides, however, little insight into the problem. To be able to use the simple conditions (2) we have to find a representation of the Hamiltonian in which the g states of interest are decoupled from all other states and in which the resulting matrix elements behave smoothly as a function of the nuclear coordinates to allow for a simple interpretation of the states. This goal is achieved in the present work with the aid of a block-diagonalization procedure which has been proposed earlier [ 8,9 ] and been shown to lead to the required smooth dependence of the matrix elements under very general conditions [ lo]. Its suitability for the present purposes is demonstrated by analyzing a multidimensional avoided crossing occurring in the ethylene dication.

2. Numerical

results

In this section we present results of configuration interaction (CT) calculations on two low-lying electronic states of the ethylene dication. In the electronic ground state the neutral molecule has the electronic configuration (symmetry group D2,.,, carbon-carbon bond on the z axis) la,2 lb:,2a,22b:,

lb&3ai lb:, lb:, .

By annihilating two electrons from the outermost orbitals lbJ, and lb?, two dicationic double hole configurations of A, symmetry are created which we denote by @,(lb?:) and &( lbyg2). The corresponding wavefunctions !P, and U: can be written as linear 274

14 October

1988

combinations of mainly o, and & and of many additional configurations with lower weight. At planar symmetry I Y,) is the ground state of the dication and 1rU,> the first excited state of A, symmetry which we denote by S, and Sz, respectively. The potential energy surfaces of these two states are the point of interest. The relevant orbitals 1bi, and lb,, are the carbon-carbon bonding n orbital and the highest otype orbital which is carbon-carbon antibonding and carbon-hydrogen bonding. In the CI calculations a double-zeta plus polarization basis set was used consisting of 4s 2p, Id Cartesian Gaussians on carbon and 2s, lp on hydrogen [ I I], The exponential scale factor for the s functions on hydrogen is 1.2. The exponent for the d functions on carbon and for the p functions on hydrogen are chosen to be 0.78 and 1.0, respectively. For the ground state of the dication the configuration space comprises all single and double cxcitations (SD) relative to the two closed-shell HartreeFock configurations @, and #r. i.e. WCperformed tworeference (2R) SD CI calculations using o, and g2 as reference configurations. The total energy of the dicationic ground state (denoted by E, ) is given by the first root of the corresponding secular equation, that of the doubly excited state (denoted by El) by the second one. Note that & results from 9, after a HOMO-LUMO pair excitation in the dication. In ethylene there are twelve nuclear degrees of freedom and ten parameters of free choice are still available after two have been committed to the task of satisfying the two double-degeneracy conditions (2). We look for geometries for which the potential energy surfaces S, and S2 are likely to approach each other closely in energy. We perform that by varying the coordinates to which the (adiabatic) energy difference AE=E, - E2 is most sensitive, like the carbon-carbon distance, the carbon-hydrogen distance and the torsion angle between the methylene groups which we denote by Rcc, RCH and 0, respectively. At first the influence of Rcr on E, and El was investigated for all other parameters fixed at the experimentally determined geometry of the neutral ground state [ 121. In fig. 1 the dependence of the total energies of S, and Sr on Rcc is depicted. The potential energy curves resulting from the interaction of the two configurations @, and & alone are depicted in fig. 1A. The diagonal elements of the matrix

Volume 15 I, number 3

CHEMICAL

PHYSICS LETTERS

14 October

1988

-77.02

Fig. I. Energies E, and Ez of the states S, and S1 as a function of the carbon-carbon distance. (A) Resulting from the interaction of @, and oz alone (solid lines); the broken lines are the diagonal matrix elements H, , and HZr. (B) Resulting from the 2R SD CI calculation (solid lines); the corresponding diabatic curves are shown as dashed lines (see sectron 3).

l-l, i.e. H,,= ($~lHl#,> (i= 1, 2), are shown as broken lines. Because of the different carbon-carbon bonding properties of the relevant lbs, and lb3, orbitals, the curves have different slopes and are seen to intersect at Rcc= 1.153 A. The interaction element HL2= (9, IEZl&) is equal to the exchange in1A exhibits a well-defined tefd &bTulb3,e Fig. avoided crossing situation. In comparison, fig. 1B shows the full results of the CI calculations for E, and E2 as a function of Rcc, i.e. the influence of the additional configurations is now taken into account. Although the relevant parameters have clearly changed, e.g. the value of R,, at the avoided crossing point and the split AE at this point are now smaller, we still find the avoided crossing. At an Rcc of la1197 8, the CI energy curves (solid lines) of the states S, and S2 come close to each other. Here AE is smallest and amounts to 0.2 eV. The dashed lines shown in the figure are so-called diabatic states [ 1315] and are discussed in section 3. At the crossing point of the diabatic states obviously the character of the adiabatic one changes completely. For an Rcc of 1.10 A, P, consists of 89.6% of &, in which two electrons are removed from the orbital 1bX,. At the degeneracy point of the diabatic states the weights of @, and $I are nearly equal in both the adiabatic states !P, and ‘Pz. For Rcc=l.14 A, !P, is mainly (85%) characterized by $, resulting from the annihilation of two electrons in the orbital lbs,. In order to investigate a multidimensional avoided

crossing additional parameters are varied. Clearly, if we restrict ourselves to the two configurations $, and & (as in fig. 1A) the avoided crossing should persist because the exchange integral is positive definite. It is, however, not a priori clear what happens if the interaction with all other possible configurations is included. We first study the behavior of the potential energy surfaces depending on the two totally symmetric coordinates R,-- and I&,. In fig. 2A a contour plot is presented which reflects the dependence of AE on these two parameters. It is seen that AR is always positive, i.e. the surfaces do not intersect. Three minima of AE occur along the diagonal of the selccted portion of the (Rcc, RCH) plane. On the employed CI level the first minimum (marked by a square) amounts to 0.25 eV, the second (see cross) to 0.20 eV and the third (see triangle) to 0.21 eV so that the central minimum is the absolute one. On both sides of the “minimum channel” along the diagonal AR increases from inside to outside in a similar fashion. Fig. 2B shows a contour plot of AE in the ( Rcc, p) plane. Here a totally symmetric (a,) and a nontotally symmetric coordinate (a,) are varied. For the sake of comparison again Rcc was chosen as the a, mode the torsion angle fi was taken as one of the non-totally symmetric modes. In the (Rcc, /I) plane shown in fig. 2B the minimum of AR (indicated by a cross in fig. 2A) is identical to the absolute one of the ( Rcc, RCH ) plane, i.e. in both planes the absolute 275

Volume 15 1,number

CHEMICAL

3

14 October

PHYSICS LETTERS

1988

8.0

1.08

1.12

RCC

1.16

1.20

(A)

1.10

1.11

1.12

RCC

1.1 3

1.14

(A)

Fig. 2. Contour diagrams of the difference AE=E, -El of the adiabatic potential energy surfaces: (A) in the (R,,, R,.,,) plane; (B) in the R(,(,, /I) plane. In both figures the energetically lowest contour line corresponds to a value of 0.27 eV (0.01 au). The values of the contour lines increase from inside to outside m steps of 0.27 eV (0.01 au). For the construction of the contour plot in (A) 45 Cl data points regularly distributed over the considered section are subjected to a two-dimensional standard interpolation providing 693 points in total. To get (B), the same interpolation was used but starting from a regular network of 78 Cl points. We obtained 855 points serving for the construction of the contour plot.

minimum occurs at the same set of geometrical parameters (Rcc=1.1197..& R,,=l.O86~andP=O). A/Z is symmetric with respect to an angle of 0 and increases in all directions leaving the minimal (adiabatic) energy gap. For a torsion angle different from zero the two highest occupied molecular orbitals (HOMO) of ethylene both belong to the symmetry bj within the symmetry group DL. Therefore, already in the neutral ground state an orbital repulsion exists caused by symmetry which increases fast with increasing torsion angle j3. The potential energy surfaces of the derived dicationic double-hole states S, and S2 repel each other strongly already at small torsion angles. The repulsion due to the configuration interaction already existing in DZh is reinforced by the orbital repulsion which additionally appears within Dz. As may be expected, the minimum of the energy gap occurs at a torsion angle of zero where the molecule belongs to DZh so that the HOMOs are of different symmetries (bju and bsg) and cannot influence each other.

3. Analysis of the results The topological character of the avoided crossing is better understood if we are able to characterize it 276

in terms of a few relevant parameters. In the simple cast, where the interactions of @, and & to all other configurations is neglected, the relevant parameters are obviously H,,, Hz2 and HL2. Such a reduction is indeed also possible in the general case, where the abovementioned configuration interactions are included. This task is achicvcd by a unitary transformation which block-diagonalizes the Hamiltonian matrix H of the system [ 8,9]. A block-diagonal matrix consists of square matrices (blocks) along its diagonal and is zero elsewhere. The eigenvalues of the transformed matrix are, by construction, identical to the eigcnvalues of H. Let us denote by B the block which possesses as eigenvalues the energies exhibiting the avoided crossing. The block-diagonalizing transformation is uniquely determined by the requirement that the eigenvectors of B are as similar as possible to those of H. The same result follows from the very appealing condition that the transformation “does not do anything but block-diagonalization” (for mathematical details, see ref. [ 91). In particular, the transformation guarantees that the elements B, (i, j= 1, 2) reduce to H, when the interaction of @, and & with the remaining configurations is neglected. B12 can be viewed as the effective exchange integral and the B,, (i= 1, 2) as the energies of effective configurations Q, and ~1~.Explicitly,

Volume 15 I,

number3

B can be computed

CHEMICAL

via the formula

B= (S,S8)1’2(S,‘)+E,S,‘(S,SB)1’2.

[ 8,9] (3)

The 2x2 matrix S, is the upper left corner of the unitary eigenvector matrix S of H where the eigenvectors are ordered in columns with increasing energy from left to right. The 2 x 2 diagonal matrix E, contains the corresponding eigenvalues E, and E, along the diagonal. The effective configurations 9, and 92 follow directly from the transformation (3). We may write Bii=<9#19,),

k,l=lJ,

(4)

where fi is the Hamilton operator. Under well-delined conditions the states 19, ) and 19Z) obtained from eq. ( 3 ) are diabatic states [ 13-I 5 ] (for details see refs. [8,10]). Using this block-diagonalization procedure the interaction of all the configurations occurring in the Hamiltonian matrix H can be replaced by the socalled effective interaction B,Z of only two states which are relevant for the problem under consideration, here the avoided crossing problem. The two conditions for doubly degenerate eigenvalues can now be explicitly expressed in terms of the matrix elements of B: B,, =Bzz B,2 =O

for

AB=B,,

-Bz2=0, (5)

meaning that both the energy gap AL? and the coupling B, z between the effective configurations vanish independently at the same time. In our specific example of the ethylene dication we may choose the matrix H to be the usual CI matrix discussed in section 2. The basis configurations $r ( lb3;2 ) and &( lbrz ) introduced there do not depend strongly on R,, and R,, because no orbital coupling is possible in DZh symmetry. The resulting states 19,) and I 9r> are diabatic states. Once the torsional mode B is activated, the symmetry is reduced to Dz. Now orbital coupling occurs and a third configuration characterized by one hole in lb?, and one in lb3, couples to di and &_ Then, all the above equations are still valid, but the simple interpretation of 9, and 92 is somewhat lost. A model of a threedimensional effective Hamiltonian B comprising &,

PHYSICS LETTERS

14 October

1988

9* and the additional configuration may restore the situation and lead to diabatic states. To gain insight into how the two conditions ( 5) are violated in our avoided-crossing situation, we investigate the behavior of AB and B,, separately in dependence on R,-, and R,,. The corresponding contour plots are given in figs. 3A and 3B, respectively. In fig. 3A the bold line represents m= 0 and comparison with fig. 2A shows that this line nicely reflects the minimum channel of hE (indicated by the symbols in figs. 2A and 3A). Consequently, the first condition can easily be satisfied and cannot cause the multidimensional avoided crossing. This behavior is in line with that of the difference between the matrix elements H,, -HZ, themselves which describes the situation in first order. In fig. 1B we have included the quantities B, , (marked by crosses) and Bzz (marked by rhombi) and the close analogy with H, , and H,, in fig. IA is apparent. In particular, and as expected for diabatic states, B, , and Bal cross at the geometry of minimal adiabatic energy gap, at z&-c= 1.1197 A. Now it is straightforward to analyze the second condition in (5) associated with the coupling of the diabatic states. If a crossing of the adiabatic states exists it can only occur where the “zero” lines of m and B,, cross. As expected for diabatic states, the contour lines of AB and B,> are smooth functions of the nuclear coordinates. In particular, the contour lines are nearly straight lines. As shown in fig. 3B, B,> decreases with increasing RCH but in the coordinate space considered no contour lint corresponds to a value of zero. This means that the multidimensional avoided crossing is caused by the non-vanishing effective coupling which behaves similarly to the exchange matrix element itself. In a last step we compare explicitly the first-order coupling Klb,.,,_ and the effective coupling B12 for a fixed R,,. Fig. 4 shows the dependence of both RI2 and N, 2=K,b3.,bYs on R,, using all the same geometrical parameters as chosen in fig. 1B. Obviously, the difference between the curves demonstrates the influence of all other remaining configurations on the coupling. These configurations weaken the interaction by a shift of 0.1 eV to lower energy, but do not change the slope of the coupling curve; both coupling terms are different from zero and increase slightly with growing Rec. Nearly parallel curves result in277

Volume ISI, number 3

CHEMICAL PHYSICS LETTERS

1.08

1.12 RCC

1.16

1.20

(A)

1.08

1.12 RCC

14 October I988

l.i6

1.30

(A)

-Bzz and (B) the effective couplingB,Z, in the (Rc.c, R,,) plane Fig. 3. Contour diagrams of (A) the diahntirenergy differenceAB=B,, (compare panel (A) with fig. 2A for the adiabatic M). The bold line in (A) is the zero line, i.e. AB=O. The first contour line in the lower right corner corresponds to a value of 5.44 eV (0.20). Towards the upper left corner the values of the contour lines decrease in steps of 0.27 eV (0.01 au). In (B) the first contour line on the top corresponds to a value of 0.027 eV (0.001 au). The values of the contour lines increase from top to bottom in steps of 0.0068 eV (0.00025 au). For the constructton points regularly distributed in total.

over the considered

section are subjected to a two-dimensional

dicating that the effect of all the other configurations on the coupling is essentially independent of the nuclear coordinate. The first-order and full coupling terms reveal the same physics, i.e. apart from an overall shift the coupling can already be interpreted within a simple two-state model neglecting all the non-relevant configurations. For the geometries presented in this article the sum

standard

of the contour plots, 45 Cl data interpolation providing 693 points

of the squared eigenvector components of Y, and YZ over #, and & lies in the range 0.89-0.92, indicating that the condition that Y, and y/2 are mainly of twohole character is fulfilled in some neighborhood of the avoided crossing points. This makes the above findings plausible and leads to the expectation that the multidimensional avoided crossing persists even if the remaining nine internal coordinales, which have been neglected in the present work, are considered.

0.25

n20

! +-

4. Summary and conclusions

0.15

7 ?I. >

_____-_-C--*

0.10

LL

! 0.05 0.00

1.10

1.11

1.12

1.13

1.14

RCC (A) Fig. 4. Coupling P’between the diabatic states as a function ofthe carbon-carbon distance. On the level of first-order perturbation theory the couphng is given by the exchange matrix element K Ib,.,b,S (solid line), which is compared with the effective coupling RI2 (dashed line).

278

The potential energy surfaces of states of the same symmetry and multiplicity may cross in the presence of sufficiently many degrees of freedom to fulfil the degeneracy conditions. The conditions for intersections of potential energy surfaces are discussed in the context of a block-diagonalization procedure. This procedure collects the states in question into a single block matrix and thus decouples them from all other states. If some care is taken in constructing the Cl matrix, the block-diagonalization yields diabatic states which simplify the discussion. In the specific case of the ethylene dication the

Volume

I 51, number 3

CHEMICAL

PHYSICS LETTERS

analysis shows that the diabatic energy surfaces intersect and the non-vanishing interaction between the diabatic states prevents the intersection of the adiabatic surfaces, giving rise to the multidimensional avoided crossing. The comparison of the effective coupling and that resulting from first-order perturbation theory indicates that the coupling mechanism can be understood qualitatively in the framework of a two-configuration problem. The dependence on the nuclear coordinates is governed by these two configurations and all the remaining configurations only introduce an overall (though substantial) constant shift of the coupling (see fig. 4). The relevant two configurations are doubly excited with respect to one another, where the excitation involves only two spatial orbitals. Hence, the coupling is given to first order by an exchange integral. Exchange integrals of this type are, for finite nuclear coordinates, always positive. Consequently, the interaction between diabatic states which arc dominated by such configurations differing by the occupancy of two orbitals is likely to be non-zero. The potential energy surfaces will then exhibit multidimensional avoided crossings in the vicinity of nuclear geometries where the diabatic surfaces intersect. The ground state and low-lying excited states of highly symmetric doubly ionized molecules are natural candidates for this type of avoided crossings. Of course, it may also appear in neutral and singly ionized species. Other, closely related cases, are also of interest. In NzOf [ 161 and CH2+ [ 171, for example, there are two states of the same symmetry lying within the same energy range. As in the dication of ethylene these electronic states are doubly excited with respect to each other, but the excitation is of a different type ($,(a-‘) and &(bP2), i.e. the double excitation from @, to q2 is described by removing two electrons from the spatial orbital b, placing one in c and one in a. The matrix element coupling of these configurations involves the three orbitals a, b and c and we could not show that it is necessarily definite. Also in these cases there are indications that no conical intersection of the adiabatic potential energy surfaces takes place. Here, the reason lies in the finding that the corresponding Coulomb matrix elements describing the coupling depend only weakly on the nuclear coordinates. This finding has been

14 October

1988

relevant for NzO+, simplifying the nuclear dynamics calculation which has been carried out to explain the experimental photoelectron spectrum [ 16 1. The construction of diabatic states is itself a vivid subject of current research. Several methods to construct diabatic states have been proposed, the discussion of which is beyond the scope of this work. One major method consists of solving differential equations to eliminate a part of the derivative coupling between the states. This method requires as input data the derivative coupling of the adiabatic states as a function of nuclear coordinates. It has been successfully used to compute diabatic states of method used CH: [ 17 1. The block-diagonalization in the present work has the advantage that it does not require the knowledge of the coupling of the adiabatic states, which are often cumbersome to calculate as a function of several nuclear coordinates. The only input data are the adiabatic potential energy surfaces themselves and a few components of the corresponding eigenvectors (see eq. (3)). The general implications of the block-diagonalization procedure on the quality of the diabatic states are discussed in detail in ref. [lo].

Acknowledgement The authors thank F. Tarantelli for many useful discussions and help throughout the work. Financial support by the Deutsche Forschungsgemeinschaft is acknowledged.

References [I ] G. Herzberg and H.C. Longuet-Higgins, Discussions FaradaySoc.35(1963) 77. [2]C.M. Meerman-van Renthem, A-H. Huizer and J.J.C. Mulder, Chem. Phys. Letters 5 1 ( 1977) 93. [3] G.J. Hatton, W.L. Lichten and N. Ostrove, Chem. Phys. Letters 40 (1976) 437. [4 ] T. Carrington, Accounts Chem. Res. 7 ( 1974) 20. [5] J. van Neumann and E. Wigner, Physik. Z. 30 (1929) 467. [6] J. Katriel and E.R. Davidson, Chem. Phys. Letters 76 (1980) 259. [ 7) H.-D. Meyer, private communication. [g] LX Cederbaum, H. Koppel and W. Domcke, QuantumChem. 15 (1981) 251.

Intern.

J.

279

Volume 15 I, number 3

CHEMlCAL

PHYSICS LETTERS

[9] L.S. Cederbaum, .I. Schirmer and H.-D. Meyer, to be published. [ 101 T. Pacher, L.S. Cederbaum and H. Koppel, to be published. [ 111 T.H. Dunning, J. Chem. Phys. 53 ( 1970) 2823. [ 121 G. Herzbcrg, Electronic spectra of polyatomic molecules (Van Nostrand, Princeton, 1966 ). [ 131L.D. Landau, Physik. 2. Sowjetunion 2 (1932) 46.

280

14 October

1988

[ 141 W. Lichten, Phys. Rev. 164 ( 1967 ) 13I [ 151F.T. Smith, Phys. Rev. 179 (1969) 111. [ 161 H. Kisppel, L.S. Cederbaum and W. Domcke, Chem. Phys. 69 (1982) 175. [ 171 M. Desouter-Lecomte, B. Leyh-Nihant, M.T. Praet and J.C. Lorquet, J. Chem. Phys. 86 (1987) 7025.