**Supplementary UserÕs Guide**

for a version of GAUSSIAN03

(rev. D.01)

modified to allow nuclear charge
stabilization calculations for metastable states

This version of GAUSSIAN03

has been modified and tested by

Prof. Piotr Skurski[1]
while working in

The Simons Group at The

Henry Eyring Center for
Theoretical Chemistry

Department of Chemistry,
University of Utah

Salt Lake City, UT 84112, U.S.A.

**Introduction**

** **

When studying
metastable electronic states that can decay by ejection of an electron, one
cannot straightforwardly employ variation-based electronic structure tools[1].
For such cases, the lowest-energy state corresponds to that of a free electron
(infinitely distant and with zero velocity) plus a system with one fewer
electrons. Standard methods thus suffer so-called variational collapse and tend
to converge to such Òfree-electronÓ descriptions. LetÕs consider an example
offered by the formamide molecule near its equilibrium geometry. The p-symmetry molecular orbitals (MOs) are depicted qualitatively below
where their nodal characteristics are emphasized. The lowest two p MOs describe the delocalized p bonding and
non-bonding orbitals. The unoccupied MO is the anti-bonding p* orbital.

An SCF calculation using a
correlation consistent augmented polarized valence double-zeta (aug-cc-pVDZ)
basis set produces the orbitals shown below. The orbital energies for the
bonding and non-bonding p MOs (labeled HOMO-2 and HOMO) are -15.4 and -11.5 eV,
respectively. The HOMO-1 orbital is a lone pair orbital on the oxygen atom.

The
SCF orbital energy of the lowest unoccupied molecular orbital (LUMO) is 0.72
eV, which suggests that an electron of 0.72 eV might attach to produce the
formamide anion. However, as the picture shown below illustrates, the LUMO is
not even of p* symmetry, nor is the LUMO + 1 or the LUMO + 2 orbital.
In fact, these three unoccupied orbitals do not have any significant valence
character; most of their amplitude is outside the formamide moleculeÕs
molecular skeleton. They are, within the finite atomic orbital basis used,
approximations to the free-electron orbital.

The
lowest unoccupied orbital of p* character is the LUMO + 3,
which is also shown below, and this orbital has an energy of + 2.6 eV. However,
in a different atomic orbital basis, the lowest unoccupied orbital of p* symmetry would not necessarily be LUMO + 3, nor would it necessarily
have an orbital energy near 2.6
eV.

To give the reader an idea of
how these orbitals look for other atomic orbital basis sets, we show below four
of the orbitals obtained when a 6-31G** basis is employed.

We see that the HOMO-2 and
HOMO are still p bonding and non-bonding orbitals and HOMO-1 is still
an oxygen lone pair orbital, but now LUMO is the p*
anti-bonding orbital. It is important to notice that the desired p* orbital may be the LUMO in one basis but might be another orbital in
a different basis as it is in the examples shown above.

Clearly,
estimating the vertical electron affinity (EA)
by using the energy of the LUMO is wrong, but we need to make clear what is
wrong and how to properly identify the correct orbital. In fact, most of the
low-energy vacant orbitals are nothing but attempts, within the finite orbital
basis used, to represent a free electron plus a neutral formamide molecule.
This illustrates the variational collapse problem mentioned above.

Let
us now explain how we can fix this error. To focus on the problem, we consider
the effective potential that an electron approaching a formamide molecule from
afar would experience. Also keep in mind that, as is conventional in electron
scattering calculations, we may wish to decompose the wave function describing
this ÒattachedÓ electron into products of radial and angular terms:

Y(r,q,f) = S_{L,M }y_{L,M }(r) Y_{L,M}(q,f).

Substituting such an
expansion into the Schršdinger equation

{- h^{2}/2m r ^{-2 }¶/¶r(r^{2}¶/¶r)
+ h^{2}L^{2}/2mr^{2}
+ V} Y = E Y

produces the usual set of
coupled-channel equations. Multiplying the Schršdinger equation on the left by
Y*_{L,M} and integrating over the angles q and f, one obtains a radial equation for the y_{L,M}
component of the wave function:

{- h^{2}/2m r^{ -2 }¶/¶r(r^{2}¶/¶r) + h^{2}<L^{2}>/2mr^{2} + <V>}y_{L,M} = E_{ }y_{L,M}

where <V> denotes the
angular average of the electron-molecule potential

<V>
= ˜ Y*_{L,M }V Y_{L,M
}sinq dq df.

When the distance r of the
electron from the molecule is large, <V> is dominated by an attractive
Coulomb factor if the electron is interacting with a cation, by a charge-dipole
potential if the electron is interacting with a polar neutral molecule, and by
a repulsive Coulomb factor if the electron is interacting with a negative ion.
However, the centrifugal potential h^{2}<L^{2}>/2mr^{2}^{ }always varies as r^{-2}.

For the formamide
case at hand, the longest-range component of <V> is the charge-dipole
term, which varies as r^{-2}. Moreover, the nodal character of the p* orbital into which the incoming electron is to attach has dominant L
= 3 character. To see this, we view this orbital from a long distance at which
its three nuclear centers are nearly on top of one another as shown below.

When viewed as having the O,
C, and N nuclei on top of one another, this p* orbital clearly
has nodal properties like that of an f-orbital which is why L = 3 is dominant
at large-r.

An electron in an
orbital having angular momentum L experiences an effective radial potential
(i.e., the sum of <V> and the centrifugal potential) that varies as shown
qualitatively below for a neutral molecule. For an electron interacting with an
anion, the repulsive long-range part of the potential would also include a
Coulomb term e^{2}/r. In such cases, the barrier that acts to constrain
the electron in the metastable state arises from both Coulomb and centrifugal
factors. Nevertheless, the methods discussed here can still be used to
characterize the metastable stateÕs energy and wave function.

If, as in the case of
formamide, the component of the potential generated by the attractive
valence-range influences of the O, C, and N centers is not strong enough, no
bound state of the anion will exist[2].
For such a case, only the metastable so-called L = 3 shape resonance state will
occur and it will have an energy (the heavy horizontal line) and a radial wave
function as shown below.

In the valence region, the
resonance function has large amplitude, suggesting that the electron is rather
localized, it decays exponentially in the classically forbidden tunneling
region, and it has sinusoidal oscillations beyond this region with the local de
Broglie wavelength relating to the electronÕs kinetic energy.

However, at energies
both above and below that of the shape resonance state, there exists a
continuum of other states. Those lying below the resonance energy vary with r
as shown below.

They have little amplitude in
the valence region and have large amplitude at large-r and they have longer de
Broglie wavelength and thus lower energy than the resonance state. There are
also non-resonant states lying energetically above the shape resonance; they
have little amplitude in the valence region and larger amplitude at large-r but
with shorter de Broglie wavelength corresponding to higher kinetic energy.

So, how can we
identify and characterize the resonance wave function (and its energy) when it
is ÒburiedÓ within a continuum of higher- and lower- energy states? The example
discussed earlier makes it clear we cannot rely on the orbital energy nor on
choosing the LUMO to find the desired state. In the nuclear charge
stabilization method[3], we attempt
to smoothly slightly enhance the valence-region attractive character of the
potential V that the attached electron experiences to an extent that pulls the
energy of the metastable resonance state below zero thus rendering it stable[4].
We do this by altering the nuclear charges of those atoms over which the
resonance MO is expected to be localized in the valence region. For formamide,
we alter the O, C, and N nuclear charges from their actual values of 8, 6, and
7, respectively, to, for example, 8+dq, 6+dq, and 7+dq, respectively. The nuclear charges appear in the electronic Hamiltonian
in the electron-nuclear interaction potential

V_{e,n
}= S_{j=1,N }S_{K=1.M }(-Z_{K}e^{2}/|**r**_{j}-**R**_{K}|)

where the sum over j runs
over all N electrons, the sum over K runs over all M nuclei of charge Z_{K}.
So, the nuclear charge scaling in formamide introduces an additional
one-electron potential

dV_{e,n }= - dqS_{j=1,N }[_{ }(e^{2}/|**r**_{j}-**R**_{O}|) _{ }(e^{2}/|**r**_{j}-**R**_{C}|) + (e^{2}/|**r**_{j}-**R**_{N}|)

that acts to differentially
stabilize the electronÕs potential energy near these nuclear centers. Visually,
we can represent the effects of the scaled nuclear charges on the
electron-molecule potential as shown below.

Here, in black, we reproduce
the original unscaled potential and its resonance wave function. In blue and
red, we illustrate potentials in which the nuclear charges have been
incremented by small (blue) or somewhat larger (red) amounts. The key thing to
notice that is, if the nuclear charge enhancement is large enough, the
valence-range component of the potential will be lowered enough to render the
resonance state bound (notice the lack of continuum components at large-r and
notice that the energy lies below zero) rather than metastable and thus
amenable to studying using conventional variational-based quantum chemistry
tools. In this manner, one can carry out conventional calculations on the nuclear
charge enhanced species for a series of dq values (all of
which must be large enough to render the desired state bound) and then
extrapolate to dq ¨ 0 to allow the resonance state to be identified from
the finite-dq calculationsÕ data.

In practice, to use
the nuclear charge stabilization method, one (for a series of dq values)

1. Identifies those nuclei
over which the valence component of the desired resonance stateÕs orbital will
be distributed;

2. Modifies the nuclear
charges of these nuclei (one need not use equal charge increments for all the
nuclei[5])
by increasing them by an amount dq;

3. Carries out a standard,
bound-state, quantum calculation (SCF, MPn, coupled-cluster, or whatever) on
the electron-attached[6]
and non-attached states using the scaled nuclear charges to obtain attached
(E*) and non-attached (E) state energies, after which

4. One plots the energy
difference (E-E*) vs. dq but using only dq values large
enough to render E-E* > 0 (i.e., to make the attached state bound relative
to the non-attached state), and

5. One then extrapolates the
plot to dq ¨ 0 to obtain an estimate of the energy of the
electron-attached state relative to that of the parent non-attached species
(i.e., in the extrapolation, one will find E-E* negative, meaning the electron
attached species lies above its parent).

In this manner, one achieves
a nuclear charge stabilization estimate of the energy of the desired metastable
state.

In the figure shown
below, we show KoopmansÕ theorem[7],
SCF-level, MP2-level, and coupled-cluster level energy differences (D) for the
formamide and formamide p* anion obtained using the aug-cc-pVDZ basis set discussed earlier. Only for
values of dq equal to 0.18 or larger was the anion found to be
electronically stable, so only data obtained with dq = 0.18 and higher were used in constructing these plots.

The dq ¨ 0 extrapolated energies for each case are as follows:

D = –4.3 eV (KT, red)
and the fitÕs correlation coefficient (squared) is: r^{2 }=^{ }0.99957.

D = –3.6 eV (SCF,
blue), r^{2 }=^{ }0.99996

D = –3.1 eV (MP2,
green), r^{2 }=^{ }0.99994

D = –3.1 eV (CCSD(T),
magenta), r^{2 }=^{ }0.99994.

When
analogous nuclear charge stabilization calculations are performed using the
6-31G** basis, the plot shown below results.

For this basis set, the
extrapolated predictions of the D-values are as follows:

D(KT) = -5.6 eV

D(SCF) = -4.3 eV

D(MP2) = -4.4 eV

D(CCSD(T)) = -4.4 eV.

These extrapolations
represent the nuclear charge stabilizationÕs predictions for the energy of the
anionÕs metastable p* resonance state relative to that of the neutral for
each level of theory.

In
addition, below we show the p* LUMO obtained for the
aug-cc-pVDZ basis
set for dq = 0.19 and 0.26 to show how it differs significantly
from the LUMO obtained in the non-scaled calculation which, as we explained
earlier, cannot be trusted to relate to the desired resonance state.

Clearly, this orbital has the
desired valence character and nodal pattern. For dq = 0.26 it
is more compact than for dq = 0.19 because of the enhanced nuclear attraction in
the former case. For the 6-31G** basis, this same orbital is shown below also
for dq = 0.16 and dq = 0.29.

We see that the qualitative
character of the desired p* orbital does not depend on the basis set employed
although the quantitative values of the resonance state energy do. For the
larger basis and using the highest-level of theory (the CCSD(T) data), the p* attached anion is predicted to lie 3.1 eV above the neutral formamide
at the neutralÕs equilibrium geometry. Experiments using so-called electron
transmission spectroscopy methods find[8]
a resonance state to lie ca. 2 eV above the neutral, so one can see that
obtaining accurate estimates of the energies of such metastable states is
difficult even when rather good atomic basis sets are employed.

**Details for using the G03 modifications.**

** **

There are
two steps that must be taken; one must identify those nuclei whose charges are
to be altered and one must tell by how much to increase or decrease the charge
of each nucleus.

**Specifiying
the nuclear centers whose charges are to be modified**

**IOp(3/110=É)
**and** IOp(7/110=É)
**

** **

These IOps are used to
tell Gaussian the **highest numbered**
nuclear center range of atoms (centers in the geometry input list in your
Gaussian .com file) in which there are the nuclei whose charges are to be
modified.

Example #1:

*To modify only the
charge on center number 1, one needs to specify ***IOp(3/110=1) ***and*** IOp(7/110=1)**.

Example #2:

*To modify two nuclei
whose numbers (in the geometry input specification) are 1 and 2, one needs to
specify* **IOp(3/110=2)** *and*
**IOp(7/110=2).**

Example #3:

*To modify four nuclei
whose numbers (in the geometry input specification) are 13, 26, 27, and 67, one
needs to specify* **IOp(3/110=67)** *and*
**IOp(7/110=67)**.

Briefly, the values of __option
110__ in __overlay 3__ and __option 110__ in __overlay 7__ have to
be set to the highest numbered center that is to be modified, no matter how
many centers will actually be modified. Note that option 110 should have the
same value in both overlays (i.e., IOp(3/110=X_{1}) and IOp(7/110=X_{2}),
where X_{1}=X_{2}).

Warning: the largest
possible value for this option is 90 in both overlays (see the ÒLimitations and
known problemsÓ section at the end of this manual), so make sure you number
your atoms so that all of the atoms you intend to subject to charge scaling
occur within the first 90 atoms.

**Modifying
the nuclear charges on the selected centers**

**IOp(3/111É200)
**and**
IOp(7/111É200) **

These IOps are used both
to denote the centers to be modified and to specify the additional charge (dq) that is to be added to that center (or subtracted
from it).

To modify nuclear charge
on center #N, the IOptions 3/(110+N) and 7/(110+N) must be set to a proper
value. The additional charge to be
added to (or subtracted from) the N^{th} nucleus must be multiplied by
10^{5} (100000) and set as the value in both IOptions 3/(110+N) and
7/(110+N).

Example #1:

*To add the charge of
0.1 a.u. to the center number 1, one needs to specify ***IOp(3/111=10000) ***and***
IOp(7/111=10000)** because:

Ÿ the center number is 1 which corresponds to options
111 (110+1) in overlays 3 and 7 , and

Ÿ the charge to be added is 0.1 a.u. which
corresponds to values of 10000 for options 111 in overlays 3 and 7 (because
0.1×10^{5} = 10000).

Example #2:

*To modify two nuclei whose
numbers in the geometry input specification are 4 and 17, by adding an extra
positive charge of 0.03 a.u. to nucleus on center 4 and charge of 0.58 a.u. to
the nucleus on center 17, one needs to specify* **IOp(3/114=3000,3/127=58000)** **IOp(7/114=3000,7/127=58000)** because

Ÿ the number of the first center to be modified is 4
which corresponds to options 114 (110+4) in overlays 3 and 7 , and

Ÿ the number of the second center to be modified is
17 which corresponds to options 127 (110+17) in overlays 3 and 7 , and

Ÿ the charge to be added to the first modified center
is 0.03 a.u. which corresponds to values of 3000 for options 114 in overlays 3
and 7 (because 0.03×10^{5} = 3000) for this center, and

Ÿ the charge to be added to the second modified
center is 0.58 a.u. which corresponds to values of 58000 for options 127 in
overlays 3 and 7 (because 0.58×10^{5} = 58000) for this center.

Recall, that even though
you are choosing the center numbers here, you still need to set the proper
range by using IOp(3/110=É) and
IOp(7/110=É)

(see the preceding
section).

**Input
examples**

Here we present a few
examples of Gaussian03 input files showing single-point energy calculations and
geometry optimization jobs that involve modified nuclear charges.

Example #1:

%chk=testSP

#p RHF 6-31G(d)
SCF=(Tight,NoVarAcc,Maxcycle=512,IntRep)

GFInput IOp(6/7=3) IOp(3/32=1)

IOp(3/110=3,3/111=33000,3/112=34000,3/113=33000)

IOp(7/110=3,7/111=33000,7/112=34000,7/113=33000)

Test job #1

0 1

O

C 1 co2

N 2 nc3 1 nco3

H 2 hc4 1 hco4
3 dih4

H 3 hn5 2 hnc5
1 dih5

H 3 hn6 2
hnc6 1 dih6

co2
1.277811

nc3
1.496812

nco3
116.154

hc4
1.002000

hco4
135.119

dih4
172.107

hn5 1.004634

hnc5
107.935

dih5 -119.695

hn6
1.002000

hnc6
112.357

dih6 0.000

In this job, the first
three atoms are modified: charges of 0.33, 0.34, and 0.33 a.u. are added to the
oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting
charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while
those for hydrogen atoms remain unchanged. This is achieved in the 4^{th}
and 5^{th} line of the input file. The Hartree-Fock single point energy
for such a system will be calculated, with the 6-31G(d) basis sets, for the
given geometry.

Example #2:

%chk=testOPT

#p UHF 6-31G(d)
SCF=(Tight,NoVarAcc,Maxcycle=512,IntRep)

GFInput IOp(6/7=3) IOp(3/32=1)
OPT GUESS=READ

IOp(3/110=3,3/111=33000,3/112=34000,3/113=33000)

IOp(7/110=3,7/111=33000,7/112=34000,7/113=33000)

Test job #2

-1 2

O

C
1 co2

N
2 nc3 1 nco3

H
2 hc4 1 hco4
3 dih4

H
3 hn5 2 hnc5
1 dih5

H
3 hn6 2 hnc6 1 dih6

co2
1.277811

nc3
1.496812

nco3
116.154

hc4
1.002000

hco4
135.119

dih4
172.107

hn5
1.004634

hnc5
107.935

dih5 -119.695

hn6
1.002000

hnc6
112.357

dih6
0.000

In this job, the first
three atoms are modified: charges of 0.33, 0.34, and 0.33 a.u. are added to the
oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting
charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while
those for hydrogen atoms remain unchanged. This is achieved in the 4^{th}
and 5^{th} line of the input file. The Hartree-Fock geometry
optimization for such an anion will be performed, with the 6-31G(d) basis sets.

Example #3:

%chk=testOPT2

#p UHF aug-cc-pVDZ
SCF=(Tight,NoVarAcc,Maxcycle=512)

GFInput IOp(6/7=3) IOp(3/32=1)
OPT GUESS=READ

IOp(3/110=5,3/111=33000,3/112=34000,3/115=33000)

IOp(7/110=5,7/111=33000,7/112=34000,7/115=33000)

Test job #3

-1 2

O 0.000000 0.000000 0.000000

C 0.000000 0.000000 1.277811

H -0.700350 0.097094 1.987802

H 1.412538 -0.830293 2.498962

N 1.343557 0.000000 1.937584

H 2.094138 0.000000 1.273783

In this job, there are
also three modified atoms: charges of 0.33, 0.34, and 0.33 a.u. are added to
the oxygen atom, carbon atom, and nitrogen atom, respectively. So the resulting
charges for atoms O, C, N are 8.33, 6.34, and 7.33 a.u., respectively, while
those for hydrogen atoms remain unchanged. This is achieved again in the 4^{th}
and 5^{th} line of the input file. The Hartree-Fock geometry
optimization for such an anion will be performed, with the 6-31G(d) basis sets.
(Please note how the centers for modifications are specified via IOps and note
that the atoms whose charges have been changed are the first, second, and fifth
atoms).

**Some
limitations and known problems**

- The maximum number of centers whose charges
can be modified in one job is limited to 90, which is limited by the array
size used for storing the IOps.

- The maximum number of atoms in a molecule is
not limited other than being subject to the same limitations as in any
regular Gaussian03 input files.

- The highest possible number for a modified
center is 90. This means that the centers must be numbered so that all
those to be charge-modified occur in the first 90.

- The modified version of Gaussian03 has been
tested for single-point energy and geometry optimization jobs. Calculating
harmonic vibrational frequencies is also possible. Note that changing the
nuclear charges does not affect the atomic masses, so quantities that
depend on atomic masses (such as vibrational frequencies) are calculated
in a standard G03 manner.

- The value of the nuclear-nuclear repulsion
energy, as well as its contributions to the gradients and Hessians, are
calculated using the conventional
(non-modified) nuclear charges.

- Some modifications of the nuclear charges are
to be avoided because they can produce unexpected errors. In particular,
changing any nuclear charge by 1 a.u. (or more) leads to problems since,
de facto, this changes the chemical identity of the element. In such cases
errors reported by
*MinBas*appear (due to the comparison made by Gaussian03 in which the atomic number is compared with its nuclear charge).

- In some cases, both the SCF procedure and the
geometry optimization process can be more difficult to converge when
modified nuclear charges are used, so it is recommended that one use the
MaxCycle keyword to increase the number of maximum SCF iterations and
optimization steps to improve chances for successful convergence.

- Molecular properties such as populations
(e.g., Mulliken population analysis, ESP-fitted charges according to the
Merz-Singh-Kollman scheme, etc.), multipole moments, polarizability
tensors, and other quantities that depend on the nuclear charges are
evaluated with the modified nuclear charges.

Please report any problems
to:

Piotr Skurski (piotr@hydrogen.hec.utah.edu).

[1] Even methods such as M¿ller-Plesset perturbation
theory (MP) and coupled-cluster theory (CC) suffer from this problem because
they are based on a Hartree-Fock (HF) self-consistent field (SCF) initial
starting point that is intrinsically variation-based.

[2] For example, the ^{2}P state of NO is electronically stable relative to NO^{+} plus
an electron , but the ^{2}P state of the
isoelectronic species N_{2}^{-} is metastable. The increased
electronegativity of the oxygen atom in NO compared to that of the nitrogen
atom in N_{2}^{-} is enough to make the valence-range potential
attractive enough to make NO stable.

[3] This method was introduced by the Peyerimhoff group
(see, for example, B. Nestmann and S. D. Peyerimhoff, J. Phys. B. **18**, 615-626 (1985); J. Phys. B **18**, 4309-4319 (1985)) and is based on the stabilization
theory put forth by Howard Taylor (see, for example, Taylor, H. S. *Adv.
Chem. Phys.* **1970**, *18*,
91).

[4] Alternatively, it is possible to alter the large-r character of the wave function by either placing the system inside a large spherical ÒboxÓ of radius R. By then varying the radius R, one can vary the de Broglie wave lengths of the orbitals generated. This, in turn, would allow one to ultimately identify a box size R that, within the finite orbital basis set employed, can generate an orbital (not necessarily and probably not the LUMO) that has a large valence-range component and a large-r component whose de Broglie wave length (and thus energy) is appropriate. The variation of the ÒboxÓ radius R can, in practice, be achieved by scaling the orbital exponents of the more diffuse basis functions (larger orbital exponents represent a smaller box). However, we have found this procedure to be more cumbersome than that which we have decided to automate in making our changes to G03.

[5] That is, one can use a different value dq_{K} for the K^{th} nucleus, but we have found it
sufficient to increase the nuclear charges of all the nuclei over which the
resonance orbital extends by equal amounts. Because one obtains the prediction
of the resonance stateÕs energy by extrapolating to dq ¨
0, it should not matter much whether one uses the same or different values of dq for the different nuclei.

[6] In carrying out the SCF calculations (that also form
the starting points of the MP2 and CCSD(T) calculations), one must be careful
to singly-occupy the desired orbital (i.e., the p* orbital in
the formamide case). The increased nuclear charges will act to lower the energy
of not only the p* orbital but also of other orbitals that have
significant density on the atomic centers whose charges have been enhanced.
Thus, one may find that the energy ordering of the virtual orbitals changes as
one moves from one dq value to another. As a result, one may need to
employ the ÒALTERÓ keyword to maintain single occupancy of the desired orbital
as part of the process of generating data to make the energy plots discussed
here.

[7] We are defining the KoopmansÕ theorem electron binding energy to be the negative of the energy of the unoccupied orbital of the neutral molecule rather than the negative of the singly-occupied orbital of the anion.

[8] M. Seydou, A. Modelli, B. Lucas, K.
Konate, C. Desfrancois, and J.P. Schermann, Eur. Phys. J. D **35**, 199 (2005).