Mourad Roudjane

Toncho I Ivanov

Maria Ofelia Vieitez

Kees De Lange

Lydia Tchang-Brillet

Wim Ubachs

C A De Lange

W.-Ü L Tchang-Brillet

Extreme ultraviolet laser calibration of D 2 Lyman and Werner transitions

Keywords: vacuum UV, molecular spectroscopy, hydrogen isotopologues, precision metrology

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The hydrogen molecule is a benchmark system for the investigation of molecular structure. Notwithstanding the fact that it is the smallest neutral molecule, its spectrum is complex: H 2 and its natural isotopologues HD and D 2 do not exhibit band structures normally associated with molecules in the gas phase, but rather atomic-like spectra without obvious regularities. Moreover, due to the low nuclear masses application of the Born-Oppenheimer approximation becomes somewhat questionable. Hence, adiabatic and non-adiabatic effects on the level structure are most pronounced in hydrogen, thus making the hydrogen molecule an ideal testing ground for assessing mass-dependent effects in spectra of diatomic molecules. Finally, due to the large energetic separation between 1s and 2p levels in the H-atom, and the large binding energy of the singly-bonded ground state of the hydrogen molecule, the electronic absorption and emission spectra of H 2 and its isotopologues are in the vacuum ultraviolet (VUV) domain, which for many years has been difficult to access with high-resolution methods. A detailed understanding of the hydrogen VUV-spectra is important in astrophysics due to the large abundance of hydrogen, with the strongest absorption features formed by the Lyman and Werner band systems of H 2 and HD. D 2 is of relevance for plasma science, most notably in Tokamak environments aiming at deuterium nuclear fusion.

Over the years the VUV spectrum of D 2 has been investigated by a number of spectroscopists viz. Bredohl and Herzberg 1 , Dabrowski and Herzberg 2 , Takezawa and Tanaka 3 , and Larzilliere et al. 4 , all in connecting wavelength regions. Hinnen et al. 5 thereupon performed a study of a limited part of the D 2 spectrum using a tunable laser source in the vacuum ultraviolet at a somewhat better resolution than in the classical spectroscopic studies.

At the Meudon Observatoire in the last decade techniques of high-resolution grating spectroscopy have been applied to produce a comprehensive database on the VUV emission spectrum of hydrogen, culminating in an atlas for H 2 6 . Work is in progress to extend these studies to HD and D 2 . Before the HD lines can be assigned, the D 2 spectrum has to be investigated, because Recently investigations on the D 1 Π u → X 1 Σ + g and D 1 Π u → X 1 Σ + g 7 and the B 1 Σ + u → X 1 Σ + g systems 8 have been published. In such studies several thousand spectral lines are recorded. An internal calibration standard for such comprehensive spectra may be provided by accurately calibrated wavelength positions for a limited subset. This is precisely the purpose of the present laser-based investigation of D 2 ; some 39 D 2 lines in the wavelength range 99.98 -102.07 nm have been calibrated with an accuracy of ∆λ/λ = 6 × 10 -8 . These transitions may serve as a calibration standard in the ongoing emission studies on the

B 1 Σ + u → X 1 Σ + g Lyman and C 1 Π u → X 1 Σ + g Werner band systems in D 2 .

Experiments

The methods of performing high-resolution spectroscopy with a narrowband tunable extremeultraviolet (XUV) laser source on diatomic molecules in a pulsed and skimmed molecular beam have been developed at the Laser Center (Vrije Universiteit, Amsterdam) and described previously 9 . Coherent XUV-radiation is produced by harmonic upconversion of the output of a pulsed dye amplifier, which is injection seeded by a stabilized ring-dye laser. The visible output in pulses of 5 ns duration is frequency-doubled in a KD*P-crystal and thereupon tripled in a xenon gas jet. The pulsed XUV-radiation of nearly Fourier-transform limited bandwidth (∆ν ∼ 250 MHz) then perpendicularly intersects a molecular beam of D 2 . It should be noted that the XUVradiation produced with the present laser system possesses a resolution and accuracy more than an order of magnitude better than with the system used in Ref. 5 . Possible Doppler shifts on the central resonance wavelengths are carefully addressed through probing of different velocity clusters in the beam by delaying the laser pulses with respect to the molecular beam pulses. The absolute frequency calibration is derived from a reference standard in the visible range, based on saturation spectroscopy of molecular iodine 10 ; interpolation is performed with the use of a stabilized étalon and the calibration in the XUV-domain is then derived from multiplication by the harmonic order. The full analysis of the uncertainty budget of the XUV-spectroscopic measurements has been discussed in recent studies performed by the group of the Laser Center, for limited sets of rotational lines of Lyman and Werner bands in H 2 and HD 11,12,13 and for the recently published work on freqency calibration of Lyman transitions in H 2 14 . In Fig. 1 a typical recording of a D 2 spectral line is displayed, with the I 2 and étalon traces used for the calibration.

Results and discussion

The wavelengths of the 39 accurately determined rotational lines and the corresponding frequencies in wavenumber units are listed in Table 1. From the error budget, the uncertainty in the line positions is estimated to be ± 0.000 006 nm or ± 0.006 cm -1 in wavenumbers. Note that all wavelength values pertain to vacuum. Combination differences ∆ 20 between P(J +2) and R(J ) lines can further be verified with ground-state rotational splittings from far-infrared FTS-Raman spectroscopy 15 to test the accuracy and the internal consistency of our results. The present set of Lyman and Werner band transitions contains measurements of four pairs of R(0) and P(2) lines whose combination differences are plotted in Fig. 2. This shows the consistency of the set of results involving these eight spectral lines. The deviation of the combination differences from their average value is 0.002 cm -1 . This suggests that the determination of transition frequencies is more accurate than would follow from our estimated error budget for these pairs of transitions. A statistical analysis yields a value for the combination difference of 179.066 (1) cm -1 , in agreement with the IR data 15 given by 179.068 (2) cm -1 for the ∆ 20 energy splitting. A similar procedure can be performed for the R(1)-P(3) combination differences, where an average value of 297.533 (2) cm -1 is derived from a set of four pairs of transitions, in agreement with 297.533 Roudjane, Vieitez, Ivanov, de Lange, Tchang-Brillet, Ubachs 3 Table 1. Measured wavelengths and the corresponding frequencies of rotational lines belonging to the B 1 Σ + u ← X 1 Σ + g Lyman bands (Lv for B-X(v -0)) and C 1 Πu ← X 1 Σ + g Werner bands (Wv for C-X(v -0)) for D2. The figures in parentheses represent the uncertainties in the last digit. ∆ o-h refers to the difference between the present laser-based results and the previously obtained classical values for the transition frequencies 2 . ∆o-c lists the difference between the observed and calculated frequencies obtained in the present work. (3) cm -1 given for the ∆ 31 energy splitting 15 . This result provides an independent validation of the calibration procedures followed in the present study. Note also that the present XUV results yield improved values for the ground-state combination differences.

Transition

λ ν ∆ o-h ∆ o-c (nm) (cm -1 ) ( c m -1 ) (cm -
In the spectral range 99.98 -102.07 nm the most accurate VUV wavelength measurements of D 2 were performed by Dabrowski and Herzberg 2 , using a spectrograph with concave grating. They claimed an uncertainty of 0.15 cm -1 on transition energies of unblended lines. The differences between the presently observed transitions with these previously reported data, ∆ o-h , are listed in Table 1 and are plotted in Fig. 3. From this comparison, we find that the transition energies of the B(v =9-11)← X(0) Lyman bands from Ref. 2 compared to the present laser-based results, with a standard deviation of σ = 0.17 cm -1 . The transition energies of C(0)← X(0) Werner band lines from Ref. 2 are on average smaller by 0.04 cm -1 than the laser-based results, with a standard deviation of σ = 0.12 cm -1 . In the work on the H 2 molecule by Philip et al. 12 , the authors pointed out a systematic shift of + 0.06 cm -1 with respect to previous data from the Meudon Observatory group 6 , also using a concave grating spectrograph. This shift was attributed to the calibration procedure where reference lines emitted by an auxiliary discharge, probably with a different illumination of the grating, were superimposed on the experimental spectrum. Since a similar procedure was also used by Dabrowski and Herzberg 2 for calibration of their D 2 spectrum, a systematic shift of their data with respect to the present laser-based results may have the same origin.

A second comparison is made between the observed transition frequencies and calculations based on the theoretical framework as presented in Refs. 7,8 . The upper-level energies belonging to the excited states B 1 Σ + u and C 1 Π u were calculated 7 by solving a set of coupled equations using the newest ab initio potentials by Wolniewicz and co-workers 16,17 , taking into account the nonadiabatic coupling terms between B 1 Σ + u , C 1 Π u , B 1 Σ + u and D 1 Π u states 18 . The differences between the observed and calculated transition energies, ∆ o-c , for the investigated lines belonging to the Lyman and Werner bands, are given in Table 1. For most of the transitions satisfactory agreement was obtained. The remaining discrepancies, are probably due to nonadiabatic couplings with higher electronic states, that were neglected in the calculations.

Conclusion

We present highly accurate transition wavelengths (∆λ/λ = 6×10 -8 ) and frequencies of some 39 rotational lines of D 2 in the B 1 Σ + u ← X 1 Σ + g , (v , 0) Lyman bands, for v = 9-11, and in the C 1 Π u ← X 1 Σ + g (0, 0) Werner band. The present values of the transition wavelengths are substantially more accurate than the previously published data of Dabrowski and Herzberg 2 , and will be used, together with the highly accurate XUV measurements on H 2 13,19 , as calibration lines for VUV emission studies, in particular for the D 2 and HD molecules. UV ionization spectroscopy, using a narrowband laser system tunable in the XUV domain. The absolute calibration of this line is performed by comparing the visible laser radiation before harmonic conversion to the t-hyperfine component in the I 2 B ← X(11,2) R(85) line, marked with (*) and corresponding to 16 508.79510 cm -1 .

Figure 2: Energy splitting ∆ 20 between the J = 2 and J = 0 ground-state levels derived from combination differences between four pairs of R(0) and P(2) transitions measured in the Lyman and Werner bands. Indicated are the (weighted) average value of ∆ 20 from the four pairs and the resulting uncertainty : 179.066 (1) cm -1 .

  are on average lower by 0.16 cm -1

Figure 1 :

 1 Figure1: Recording of the R(3) line in the B ← X (10,0) Lyman band in D 2 via 1 XUV + 1 UV ionization spectroscopy, using a narrowband laser system tunable in the XUV domain. The absolute calibration of this line is performed by comparing the visible laser radiation before harmonic conversion to the t-hyperfine component in the I 2 B ← X(11,2) R(85) line, marked with (*) and corresponding to 16 508.79510 cm -1 .

Figure 3 :

 3 Figure 3: Comparison between transition frequencies obatained in the present laser-based study and those obtained from the classical work of Ref. 2 .

  Lyman linesan HD emission discharge spectrum will always be contaminated by H 2 and D 2 transitions.

	March 4, 2008	16:40	Molecular Physics	D2˙ROUDJANE˙et˙al
	2			Extreme ultraviolet laser calibration of D 2
			F o
			r
			P
				e e r
				R e v i e w
				O n
				l y

Acknowledgements

This study was supported by the EC Integrated Infrastructure Initiative action (RI3-CT-2003-506350), including its Access to Research Infrastructures Program (Project-ID: lcvu001276).


J Demaison
email: jean.demaison@univ-lille1.fr

J Liévin

Equilibrium structure of the hydrogen bonded dimer H 2 O … HF

Keywords: ab initio, structure, hydrogen bond, H2O-HF ab initio, structure, hydrogen bond, H 2 O … HF

The equilibrium structure of the hydrogen bonded complex H 2 O … HF has been calculated ab initio using the CCSD(T) method with basis sets up to sextuple-ζ quality with diffuse functions and taking into account the basis set superposition error correction. The calculations carried out confirm the importance of diffuse functions and of counterpoise correction to obtain an accurate geometry. The most important point is that the basis set convergence is extremely slow and, for this reason an accurate ab initio structure requires a very large basis set. Nevertheless, the ab initio structure is significantly different from the experimental r 0 and r m structures. Analysis of the basis set convergence and of the approximations used for the determination of the experimental structures indicates that the ab initio structure is expected to be more reliable.

Introduction

The hydrogen-bonded heterodimer H 2 O … HF has been observed in gas phase for the first time in 1975 by infrared spectroscopy [START_REF] Thomas | Proc. R. Soc. Lond. A[END_REF]. The same year, its rotational spectrum was measured [START_REF] Bevan | [END_REF]. Then, its microwave spectrum was analyzed in great detail by Kisiel, Legon, Millen et al. permitting the determination of its effective (r 0 ) structure [3,4], dipole moment [5], harmonic force field [6] and dissociation energy [7]. It was established that the equilibrium configuration at the oxygen atom is pyramidal (C s symmetry) with a double minimum potential (see Figure 1). More recently, its submillimeterwave spectrum has been investigated [8] and its empirical mass-dependent r m structure determined [9]. This dimer has been subject to many ab initio calculations leading to three important conclusions [10,11,12] : i) the nondynamical electron correlation is small ii) the basis set superposition error (BSSE) is large iii) the convergence of the basis sets is slow. However, most of these calculations are devoted to the determination of dynamical properties (potential barrier, dissociation energy, …) and there is no very high level calculation of the geometrical structure.

The goal of this work is threefold : i) to determine a reliable ab initio equilibrium structure at a high level of theory using basis sets as large as possible and to compare it to the 3 experimental structures ii) to check the structural changes upon dimerization because it is often assumed that the geometry of the monomers does not change upon dimerization iii) to check the performance of density functional theory (DFT) for the calculation of the structure of hydrogen-bonded complexes because it seems that some popular functionals perform rather poorly [13]. This is important, because the procedure used here to determine the ab initio structure of H 2 O … HF cannot be used for larger molecules. Using DFT would thus be an interesting alternative, whose accuracy is still to be checked.

Methods of computation

Most correlated-level ab initio electronic structure computations of the present study have been carried out at two levels: second-order Møller-Plesset perturbation theory (MP2) [14] and coupled cluster (CC) theory with single and double excitation [15] augmented by a perturbational estimate of the effects of connected triple excitations [CCSD(T)] [16]. The Kohn-Sham density functional theory [17] was also used in this study using a variety of functionals as described in section 4.

We used correlation-consistent polarized n-tuple zeta basis sets cc-pVnZ [18] with n ∈ {D, T, Q, 5, 6}, that are abbreviated as VnZ in the text. For hydrogen-bonded complexes, the inclusion of diffuse functions in the basis set is important. For this reason, the augmented VnZ (aug-cc-pVnZ, AVnZ in short) basis sets [19] were employed. The effect of adding a second set of diffuse functions has also been investigated by using the doubly augmented VnZ (d-aug-cc-pVnZ, dAVnZ in short) basis sets. For some DFT calculations, the Jensen's polarization consistent basis set with diffuse functions apcn (n = 2, 3, 4) was used [20].

The core-core and core-valence correlation effects on the computed equilibrium geometries [21], were estimated thanks to the correlation-consistent polarized weighted corevalence quadruple zeta (cc-pwCVQZ) [22,[START_REF]Basis sets were obtained from the Extensible Computational Chemistry Environment Basis Set Database[END_REF] basis sets. In principle, for first-row atoms, it is sufficient to use the MP2 method to estimate this correction [START_REF] Margulès | [END_REF]. The frozen core approximation (hereafter denoted as fc), i.e., keeping the 1s orbitals of F, and O doubly occupied during correlated-level calculations, was used extensively. Some geometry optimizations were also carried out by correlating all electrons (hereafter denoted as ae).

The CCSD(T) calculations were performed with MOLPRO [START_REF] Molpro | package of ab initio programs[END_REF][START_REF] Hampel | [END_REF]27] electronic structure program packages while most other calculations utilized the GAUSSIAN03 (g03) program [START_REF] Frisch | Gaussian[END_REF]. The basis set superposition error (BSSE) was taken into account by the counterpoise correction (CP) [START_REF] Boys | [END_REF]. Most calculations were performed on the HP-XC4000 cluster of the ULB/VUB computing center.

Ab initio results

The first interesting result is the value of the T 1 diagnostic [30] calculated at the CCSD(T) level whose value is 0.0069 with the wCVQZ(ae) basis set and 0.0092 with the AV5Z basis set. These values are much smaller than the cut-off value of 0.020 and thus confirm that the single-reference CCSD(T) method is well appropriate for properly describing electron correlation effects.

The equilibrium geometries obtained at MP2 level with and without CP corrections are reported in Table 1. The geometrical parameters can be visualized in Figure 1 (see also the legend for the definition of α and β). The first point to be noted is the large effect of the diffuse functions on the distances between non-bonded atoms O … H (F … O), but also on the H-F bond length. At the VQZ level, addition of a single set of diffuse functions decreases the r(O … H) bond length by as much as 0.0033 Å, while at the V5Z level, the decrease is only of 0.0004 Å. On the other hand, it is found that the addition of a second set of diffuse functions is still important at the VTZ level (decrease of r(O … H) by 0.0019 Å from AVTZ to dAVTZ, but that the corresponding effect is negligible for larger n values (0.0003 Å on r(O … H) at V5Z). Another important point is the effect of the CP correction, which is quite large for the r(O … H) bond length and for small basis sets. It is still important with the AV5Z basis set, the O … H bond length increasing by 0.0069 Å. Figure 2 shows the convergence of the O … H and O … F bond lengths as a function of the cardinal number of the basis set for the MP2 calculations with and without CP correction. It appears that the convergence is more regular when the CP correction is taken into account but it remains slow. For the O … H bond length, it is obviously not yet achieved at the AV5Z level. However, if we extrapolate r(O … H) to complete basis set (CBS) using a simple polynomial of the form

r n = r CBS + an -3 (1)
as suggested by Paizs et al. [11] and only using the AVQZ and AV5Z results (i.e. n = 4 and 5), the result is almost the same with and without CP correction: r CBS = 1.7116 Å to be compared with r CBS CP = 1.7102 Å. As a check, the r(O … H) bond length was optimized at the MP2/AV6Z level of theory with the CP correction, all other parameters being fixed at their Finally, the core correction was calculated at the CCSD(T) and MP2 levels with and without CP correction. The regular wCVQZ basis set as well as the augmented AwCVQZ basis set were used. The effect of the CP correction was also investigated, see Table 3. The conclusion is that, to calculate the core correction, it is not necessary to take into account the CP correction, and it is sufficient to use the non-augmented wCVQZ basis set. Finally, the MP2 and CCSD(T) methods give similar results, the MP2 value being slightly larger (by 0.0004 Å) for the F … O distance. Likewise, for the r(O … H) bond length, the core correction is 0.0024 Å at the MP2 level and 0.0019 Å at the CCSD(T) level. It is also worth to be noted that the core correction is large for the angle β.

The best estimate of the equilibrium structure is obtained by correcting the CCSD(T)/CBS structure for the effects of core correlation calculated at the CCSD(T)/wCVQZ level. The extrapolation to CBS is not expected to be highly accurate. The equilibrium structure is given in Table 4 where it is compared to the previous calculations and to the experimental empirical structures as well as to the equilibrium structures of the monomers. The equilibrium value for the F … O distance is 2.649 Å, it is slightly shorter than the previous ab initio estimate [11] which gave 2.652 Å and it is much shorter than the experimental empirical structures. It has already been observed that the variation of the rovibrational contribution ε = I 0 -I e upon isotopic substitution is much larger in complexes than in rigid molecules [31] mainly because they are inversely proportional to vibrational frequencies, some of them being very low in complexes. The basic assumption of the r 0 structure is that ε is negligible while the basic assumption of the substitution (r s ) structure is that ε remains constant upon isotopic substitution. Obviously, neither assumption is valid in the particular case of complexes. It explains the failure of these empirical methods. On the other hand, the r m method assumes that ε varies as the square root of the moment of inertia [32]:

ε g = c g I g g = a, b, c (2)
This gives the r m (1) method. However, for XH bonds (here X = F, O), it is necessary to take into account the variation of the bond length upon deuteration, called Laurie correction.

Watson et al. [33] assume that the apparent bond X-H and X-D consist of a common part r m which is identical for both bonds (and assumed to be close to r e ) plus an additional bond elongation due to the H or D atom asymmetrically vibrating against the rest of the molecule.

This elongation is taken proportional to the respective frequency, i.e. proportional to the reciprocal square root of the reduced mass of the vibrator. In other words, the X-H bond length is given by:

r m eff (XH) = r m (XH) +δ M m H (M -m H )       1/ 2 ( 3
)
where δ is an adjustable parameter, M is the total mass of the molecule and m H the mass of hydrogen (or deuterium).

A simple way to check the validity of Eq. ( 3) is to plot ε as a function of I . As H 2 O … HF is a near-symmetric top, the A rotational constant could not be determined for its isotopologues.

Furthermore, it is possible to define a mean rotational constant B = (B + C) 2 and a mean moment of inertia I = h 2 8π 2 B . ε is then calculated using the experimental ground state rotational constants [3] and the equilibrium structure of Table 4. A nice, almost linear relationship is found for the deuterated species with an almost parallel line for the nondeuterated species, see Figure 3. It was expected to obtain different correlations for the deuterated and non-deuterated species because of Laurie correction, Eq. ( 4). This confirms the conclusion of Kisiel [9] that the r m method can be used for H-bonded complexes. However, the accuracy of the r m structure of H 2 O … HF is limited by several factors: i) there is no isotopic substitution for fluorine ii) the coefficients c b of Eq. ( 3) and δ are fully correlated.

These two factors severely limit the accuracy of the experimental r m structure and may explain the difference found with the r e structure of Table 4. The traditional assumption of unchanged monomer geometries does not apply however for the HF moiety, the HF bond length being lengthened by as much as 0.156 Å.

Since the counterpoise correction was carried out for most calculations, the dimerisation energy is easy to determine. It is listed in table 1 for the MP2 calculations. For the CCSD(T)/AV5Z level of theory, it is ∆E = -37.54 kJ/mol. This value is slightly smaller (in absolute value) than the MP2 value (Table 1). As the dimerisation energy decreases (in absolute value) with the size of the basis set, the CBS value is expected to be close to -36 kJ/mol. Indeed, Halkier et al. 10 found ∆E = -35.61 kJ/mol using a two-points extrapolation.

This result is in good agreement with the value found by Boese et al. 34

DFT results

Calculations with the CCSD(T) method and very large basis sets (up to AV5Z), like those presented in the previous section, cannot be carried out for molecules much larger than H 2 O … HF. For this reason, it would be useful to check the accuracy of more approximate methods. H 2 O … HF is indeed a fairly strongly bond complex with dominant electrostatic contribution to the binding energy and a relatively unimportant dispersive term.

First, it is interesting to note that the MP2 and CCSD(T) methods give similar values for the F … O distance.

It is also useful to make a comparative study of the results obtained with the density functional (DFT) methods. They are indeed widely used, mainly because of their relatively low cost. There is already a large amount of papers on this subject. For instance Tuma et al.

[33] compared the results of CCSD(T) calculations with those of a large variety of density functionals. Their work was recently extended to other density functionals by Boese et al. [34]. However, in both cases, the reference structure was calculated with the AVQZ basis set.

In other words, this structure may be significantly different from the equilibrium geometry because the convergence of the basis set is not achieved. One difficulty is that the world of DFT methods is becoming a jungle of acronyms corresponding to choices of parameters more • B3LYP, Becke 3-parameter [35] Lee-Yang-Parr [36] hybrid functional which is extremely popular and is known to give reliable results in many cases [START_REF] Bauschlicher | Recent Advances in Density Functional Methods[END_REF] although its accuracy for hydrogen bonds is disputed [13].

• BLYP, Becke-Lee-Yang-Parr functional which uses the Becke 88 [START_REF] Becke | [END_REF] exchange and the LYP [START_REF] Bauschlicher | Recent Advances in Density Functional Methods[END_REF] correlation functional and which is often employed when calculating hydrogen bonds [34]. Furthermore, Jensen's [21] apcn basis sets have been optimized

for this functional.

• PW91, Perdew and Wang's exchange and correlation functional [39] which was shown to perform much better than the BLYP and B3LYP functionals for the prediction of intermolecular interaction energies [13].

A number of second-generation functionals have been proposed as successors. The most successful ones are:

• B97-1, reparametrization of Becke's B97 Functional [40] by Handy et al. [41] • BMK, Boese-Martin for Kinetics [42] which was found to perform well for hydrogenbonded complexes.

A big advantage of the DFT methods is that, as shown by Boese et al. [43,44], a basis set of triple zeta quality is sufficient. They also showed that Pople's 6-311+G(3df,2pd) basis set [45] performs quite well for DFT calculations. This confirms the conclusion of Bauschlicher et al. [START_REF] Becke | [END_REF] who found that this basis set is better than the AVTZ one. However, as a check, for some calculations, the AVTZ and Jensen's apc2 and apc3 basis sets were also used. In most cases, the CP correction was taken into account. The results are given in Table 5.

As expected, the best results are obtained with the 6-311+g(3df,2pd) basis set and with the B97-1 and BMK methods which confirms the conclusions of Boese et al. [43,45].

However, they are not more accurate than with the MP2 method; this is particularly true for the angles. In particular β values spread over a range as large as 11°, while for MP2 and CCSD(T) this range does not exceed 3°. This angle which is very sensitive to the level of theory is a good indicator of the quality of the calculation. A similar scattering is also observed for the DFT F…O and H…O distances (up to 0.050 Å) while MP2 and CCSD(T) values only change by less than 0.02 Å. Finally, it has to be noted that, for the dimerisation energy, the DFT methods are less accurate than the MP2 method. For instance, for the B97-

Barrier height

The potential energy function governing the out-of-plane bending mode is of the double minimum type. The potential barrier between the pyramidal minimum and the planar top (β = 0) of the barrier was determined experimentally by Kisiel et al. [4] to be V = 126 cm - 1 (1.5 kJ/mol). The barrier was calculated at the CCSD(T)/AV5Z level of theory giving V = 115 cm -1 in good agreement with the experimental value. The effect of the core correlation is small, decreasing the barrier by only 6 cm -1 . The effect of the zero point energy is more difficult to estimate because of the double minimum nature of the potential.

Furthermore, at the AV5Z level, the CCSD(T) energy is not yet completely converged. In conclusion, the good agreement between the experimental and the ab initio values is likely to be due to a compensation of errors.

Conclusions

Geometry optimizations were carried out for the H 2 O … HF complex using the CCSD(T), MP2, and some DFT methods and with large basis sets (up to AV6Z). The bond lengths and angles were extrapolated to the CBS limit with and without the CP correction.

The core correlation was also taken into account.

For the hydrogen bond distance r(H … O), the MP2 method gives a result close to the value found at the CCSD(T) level. On the other hand, the different DFT methods give a large scattering of values, B97-1 appearing to be the best one, although inferior to MP2.

Our best estimate of the structure has been obtained from the CCSD(T)/CBS results corrected for core correlation effects. A comparison between the CBS and the AV5Z results

indicate that the accuracy is probably better than 0.003 Å for the r(F … O) and the r(O … H) bond lengths. For the other bond lengths, the accuracy is thought to be better than 0.002 Å and for the bond angles, it is of the order of 0.1°. This ab initio structure is expected to be more reliable than the published experimental ones whose accuracy is limited by several factors analyzed in detail in this work.

Upon formation of the hydrogen bond, the intramolecular OH distance is not significantly affected, the ∠(HOH) bond angle increases by 1.1° and the intramolecular HF bond length experiences a rather large elongation of 0.156 Å.

  + CP value. The result r(O … H) = 1.7112 Å confirms the validity of the extrapolation.The CCSD(T) results are given in Table2for basis sets AVnZ with n = D, T, Q, 5, and the convergence curve is shown on Figure2for comparison with the corresponding MP2 results. Each parameter is extrapolated to CBS using Eq. (1) and the AVQZ and AV5Z results. When the CP correction is taken into account at the AV5Z level, r(F … O) = 2.652 Å to be compared with the CBS value which is r(F … O) = 2.651 Å. It again confirms the reliability of the extrapolated value. It is worth noting that the CP correction is larger when calculated at the MP2 level: the CP correction increases the r(O … H) bond length by 0.0069 Å at the MP2/AV5Z level of theory while this increase is only of 0.0045 Å at the CCSD(T)/AV5Z level.

Finally, the comparison

  of the structure of the dimer to those of the monomers shows that the H 2 O moiety is not too much affected by the formation of the hydrogen bond (the ∠(HOH) bond angle increases by 1.1° only and the r(OH) bond length is almost constant).

  hoc. Nevertheless, among the numerous exchange-correlation functionals currently used, one should consider:

16 Figure 1 .Figure 2 . 18 Figure 3 .

 1612183 Figure 1. Equilibrium geometrical structure of H 2 O-HF. Parameter β is the angle between the F-O direction and the bisectrix of the HOH angle. Parameter α corresponding to the OFH angle characterizes non-linearity of the F-H-O skeleton. This angle is not visible on the Figure because of its small value (< 2 degrees).

  using the W2 theory: ∆E = -36.35 kJ/mol but it is significantly different of the experimental value: ∆E = -42.68 kJ/mol.7 However, this last value was obtained from the zero-point dissociation energy,

D 0 = 34.3(3) kJ/mol, which seems to be compatible with the ab initio values. The observed discrepancy might be explained by the difficulty to estimate a meaningful zero-point energy for a floppy molecule with a double minimum potential.

Table 1 .

 1 MP2 structures for H 2 O … HF (distances in Å, angles in degrees and dimerisation energy in kJ/mol).

	Page 11 of 21				Molecular Physics
	basis set	r(OH)	r(HF)	r(F … O)	∠α	∠(HOH)	∠β	r(O … H)	∆E a
	VTZ	0.9607	0.9344	2.6426	2.308	104.655 126.355	1.7094
	VQZ	0.9591	0.9344	2.6407	1.747	105.071 129.468	1.7070
	V5Z	0.9591	0.9354	2.6428	1.491	105.307 132.130	1.7079
	AVDZ	0.9668	0.9427	2.6575	1.398	104.684 129.522	1.7152
	AVTZ	0.9625	0.9398	2.6432	1.446	105.093 130.684	1.7039
	AVQZ	0.9599	0.9367	2.6400	1.424	105.283 132.365	1.7037
	AV5Z	0.9594	0.9361	2.6432	1.439	105.353 132.391	1.7075
	dAVDZ	0.9670	0.9424	2.6495	1.226	104.667 128.621	1.7075
	dAVTZ	0.9626	0.9401	2.6416	1.550	105.056 129.926	1.7020
	dAVTZ + CP	0.9624	0.9393	2.6610	1.466	105.070 131.841	1.7222	-38.69
	dAVQZ	0.9600	0.9368	2.6404	1.413	105.268 131.797	1.7041
	dAV5Z	0.9595	0.9361	2.6434	1.465	105.358 132.314	1.7078
	AVDZ + CP	0.9669	0.9406	2.7008	1.485	104.716 130.595	1.7606	-38.34
	AVTZ + CP	0.9623	0.9389	2.6615	1.437	105.093 131.986	1.7231	-38.32
	AVQZ + CP	0.9599	0.9359	2.6538	1.445	105.264 132.602	1.7184	-38.08
	AV5Z + CP	0.9594	0.9357	2.6497	1.458	105.369 132.627	1.7144	-37.74
	a ∆E = dimerisation energy.

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 2 .

 2 CCSD(T) structures for H 2 O … HF (distances in Å and angles in degrees).

						Molecular Physics	Page 12 of 21
	basis set	r(OH)	r(HF)	r(F … O)	∠α	∠(HOH)	∠β	r(O … H)
	AVDZ	0.9672	0.9403	2.6646	1.593	104.723	129.007	1.7249
	AVTZ	0.9625	0.9372	2.6477	1.653	105.137	130.090	1.7111
	AVQZ	0.9598	0.9341	2.6446	1.638	105.357	131.904	1.7110
	AV5Z	0.9593	0.9335	2.6477	1.632	105.420	132.030	1.7147
	CBS	0.9588	0.9329	2.6509	1.626	105.490	132.160	1.7186

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 3 .

 3 Calculation of the core correlation for H 2 O … HF (distances in Å and angles in

			Molecular Physics
	degrees).
	method a		r(OH) r(HF) r(F … O) ∠α	∠(HOH)	∠β r(O … H)
	MP2/wCVQZ	fc	0.9592 0.9344 2.6419	1.739 105.082 129.833 1.7082
		ae	0.9583 0.9339 2.6391	1.725 105.219 130.360 1.7059
	fc -ae		0.0009 0.0005 0.0028		   3 -0.137 -0.526 0.0023
	MP2/AwCVQZ fc	0.9599 0.9367 2.6422	1.410 105.300 132.594 1.7059
		ae	0.9590 0.9362 2.6393	1.410 105.438 133.104 1.7035
	fc -ae		0.0009 0.0005 0.0029		     -0.138 -0.510 0.0024
	MP2/wCVQZ+CP fc	0.9587 0.9330 2.6657	1.226 105.088 135.503 1.7330
		ae	0.9578 0.9325 2.6628	1.208 105.226 136.105 1.7306
	fc -ae		0.0009 0.0005 0.0029		   8 -0.138 -0.602 0.0024
	CCSD(T)/wCVQZ fc	0.9591 0.9318 2.6464	1.743 105.206 131.459 1.7153
		ae	0.9582 0.9313 2.6440	1.730 105.332 131.966 1.7134
	fc -ae		0.0009 0.0005 0.0024		   13 -0.126 -0.506 0.0019

a fc = frozen core approximation; ae = all electrons correlated. Page 13 of 21 URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 4 .

 4 Structures of H 2 O … HF (distances in Å and angles in degrees).

					Molecular Physics
				r(OH)	r(HF)	r(F … O) ∠α ∠(HOH)	∠β	Ref.
	r e (H 2 O) a	0.9578			104.50	46
	r e (HF) a		0.9168		47
	MP2/6-311++G(2d,2p) 0.959	0.934 2.663	2.1	105.24 130.5	48
	MP2/AVQZ			2.652	1.4	132.6	11
	r 0	b			2.676(3)	143.2(35) 9
	r m	c			2.6593(7)	131.05(44) 9
	r e	a		0.9579 0.9324 2.6485	1.61	105.62 132.67	pi d
	a Equilibrium structure.

b Empirical effective structure. c Empirical mass-dependent structure. d pi = present investigation. Page 14 of 21 URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 5 .

 5 DFT structures for H 2 O … HF (distances in Å and angles in degrees).

	Method	basis set a	r(OH)	r(HF)	r(F … O)	∠α	∠(HOH)	∠β	r(O … H)
	B3LYP	AVTZ no CP	0.9628	0.9450	2.6366	1.618	106.112	130.767	1.6922
	B3LYP	AVTZ	0.9628	0.9448	2.6385	1.609	106.121	131.008	1.6943
	B3LYP	apc3	0.9614	0.9423	2.6400	1.598	106.158	131.488	1.6983
	B3LYP	6-311+G(3df,2pd)	0.9618	0.9419	2.6454	1.686	106.103	131.873	1.7041
	PW91	6-311+G(3df,2pd)	0.9692	0.9548	2.6211	2.028	105.327	124.676	1.6672
	BMK	6-311+G(3df,2pd)	0.9582	0.9356	2.6517	1.904	106.672	135.661	1.7169
	B97-1	6-311+G(3df,2pd)	0.9606	0.9385	2.6512	1.922	105.724	129.942	1.7135
	B97-1	apc2	0.9608	0.9397	2.6456	1.800	105.795	129.723	1.7066
	B97-1	apc3	0.9605	0.9391	2.6452	1.819	105.773	129.516	1.7068
	B97-1	AVTZ	0.9617	0.9415	2.6433	1.789	105.743	129.185	1.7025
	BLYP	acp2	0.9721	0.9563	2.6622	1.634	105.492	127.416	1.7065

a With CP correction unless otherwise stated.

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

This work was sponsored by the Fonds National de la Recherche Scientifique (FNRS, contracts FRFC and IISN) and the « Action de Recherches Concertées de la Communauté française de Belgique ». It was also performed within the "LEA HiRes " and the EU project QUASAAR (MRTN-CT-2004-512202). Dr. Demaison is indebted to ULB for providing an invited international chair to him (2006). The authors thank Professor M. Herman for his support.


Florent ; Goujon

Christine Bonal

Benoit Limoges

Patrice Malfreyt
email: patrice.malfreyt@univ-bpclermont.fr

Molecular simulations of grafted metal-chelating monolayers: methodology, structure and energy

Keywords: molecular dynamics, long range electrostatic interactions, monolayers, grafted surface, methodology

come

Introduction

Monolayers of organic molecules grafted onto solid surfaces have been the focus of intense research because of their potential use in a plethora of technological applications ranging from biosensing applications to microfabrication technology [1]. One of the most investigated monomolecular layer on solid surfaces are the self-assembled monolayers formed via strong chemisorption of alkanethiols on gold surfaces [START_REF] Ulman | An Introduction to Ultrathin Organic Films[END_REF][START_REF] Bain | [END_REF][4][5].

An important aspect of such interfacial engineering is to control the organization, orientation, and functional properties of immobilized molecules at the monolayer. This is particularly important for the immobilization of proteins that are relatively sensitive to unfolding, steric hindrance and inactivation on surface. Consequently, during the last years, considerable efforts has been done to develop new immobilization strategies that circumvent the problem of biomolecules inactivation on surfaces. Among these strategies, the method based on metal-ion chelating agents and their specific affinity for proteins containing histidine residues was proved particularly useful to control protein orientation on surfaces and to preserve their function [6]. In this system, the imidazole moieties of the protein can chelate to the free coordination sites of divalent transition-metal ions such as Cu(II), which are themselves immobilized onto solid supports via appropriate chelating agents such as the tetradentate nitrilotriacetic acid (NTA) ligand. Self-assembled monolayers of thiol-functionalized-NTA on gold supports have been the most widely investigated. However, the use of gold as solid support presents a number of disadvantages such as the oxygen sensitivity of thiol compounds, the moderate stability of the self-assembled monolayer upon storage. To overcome some of these problems, carbon surfaces were recently functionalized by self-assembled monolayers of NTA-modified loaded vith Cu(II) [7] (Figure 1).

A better understanding of the general properties governing these monolayers of metal-chelating ligands attached to electrodes is required. In fact, interfacial self-assembled monolayers properties are determined by molecular structure on the surface which is strongly affected by molecular intra-and inter chain interactions and interactions with solvent as well as with the substrate [8,9]. It also appears that there is still a dire need for probes capable of providing structural information to know the details of molecular-scale surface architecture. Computer simulation has already made important contribution to the understanding of the processes involved about molecules and reactions and is now a vital adjunct to experimental studies [START_REF] Leach | Molecular modelling principles and applications[END_REF]. Molecular simulation employ force field functions that calculate the overall potential energy of a system based on a summation of individual atom-atom pair interactions. The force field equations take into account the contributions from binding interactions and non-binding interactions. Following the pioneering work of Klein's group [START_REF] Hautman | [END_REF]12], a number of simulations have employed force field to explore the structure and dynamics of self-assembled monolayer (SAM). However, the majority of the work has been limited to the study of alkanethiols of different tail-groups which self assembled on Au surface in the presence of water [START_REF] Hautman | [END_REF][12][13][14][15][16][17][18][19][20][21]. Several molecular dynamics and Monte Carlo computer simulations studies of carboxylic acid-functionalized self-assembled monolayers [22][23][24] and of oligo(ethylene oxide) (OEO) or poly(ethylene oxide) (PEO) self-assembled monolayers have also recently been reported [25][26][27]. However, very few theoretical chemical studies of monolayers on others surfaces that gold have been reported [28][29][30][31] and none of them have been focused on negatively charged monolayers of chelating metal compounds immobilized onto carbon surface via the tetradentate nitrilotriacetic acid.

The presence of a surface creates a nonuniformity of the local density along the direction normal to the surface and gives rise to important issues concerning the truncation procedures and the methodology used to account the Coulombic interactions. Additionally, this type of system is represented by a slab geometry in molecular simulations where periodic boundary conditions are used in two of the three directions and physical boundaries preventing molecules from leaving the primary cell. For this slab geometry, the calculation of the long-range Coulombic interactions cannot be performed directly using the standard Ewald summation method [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Heyes | [END_REF] which is the most widely applied technique for this kind of interactions in a three-dimensionally periodic system. Different methodologies can be thus established. Some of them [34][35][36] are designed to calculate long-range interactions in systems periodic in two directions but bound in the third one. The Hautmann-Klein Ewald method (HKE) [35] is one of these methods and corresponds to an adaptation of the Ewald method for systems which are periodic in two dimensions only. Another alternative consists of applying the standard three dimensional Ewald summation method (EW3D) to a simulation box elongated in the z direction so that a sufficiently large empty space is left between the periodic images [37,38] to dampen out the interslab interactions. A more elaborated treatment of the long-range interactions (EW3DC method) consists of modifying the EW3D method by the addition of a correction term to remove the total dipole moment of periodically repeating slabs.

It is essential from a computational viewpoint to compare the different alternatives to estimate the long-range coulombic interactions. The HKE method is in fact well adapted for slab-like systems however, it is very time consuming with a simulation time 10-20 times higher than that of EW3D. Two reference systems are then studied using either one two-dimensional method (the Hautman-Klein method) or two three dimensional methods (EW3D, EW3DC). These reference systems are built in such a way that the size of the system allows the comparison in a reasonable simulation time. We then perform a series of comparisons on the structural, dynamical and energetical properties. The first reference system is constituted by only water molecules embedded between two parallel graphite plates. We have then extended the comparison on a second reference system of our main interest which concerns the negatively charged monolayers of chelating copper compounds immobilized onto a graphite surface in an aqueous solution. Finally, we apply the EW3DC methodology to obtain a molecular description on the structure, hydration, dynamics and energetic properties of a more sophisticated system of these monolayers.

The outline of the paper is as follows. In section 2, we present the potential model with the different methods for calculating the long-range Coulombic interactions. In section 3, we discuss the results as a function of the method used by paying attention on the system formed by monolayers of charged copper complexes covalently attached to a graphite surface in water. We focus on the electrical, structural and diffusion properties on the grafted NTA-Cu(II) ligands using the EW3DC method. We conclude in section 4 by a brief summary of our main results.

U = bonds k b (r -r o ) 2 + angles k θ (θ -θ o ) 2 + dihedrals k φ [1 + cos(lφ + δ)] + N -1 i=1 N j=i+1 4 ij σ ij r ij 12 - σ ij r ij 6 + ∞ l q i q j |r ij + nL| (1)
where k b , k θ and k φ are the force constants for deformation of bonds, angles and dihedrals, respectively.

The equilibrium values of bond distances and valence angles correspond to r o and θ o , respectively.

In the dihedral angle term, l is the periodicity and δ is the phase factor. The C-H and O-H covalent bonds were kept of fixed length by use of the SHAKE algorithm [40] and the aromatic rings were kept planar using six improper torsional potentials. The intermolecular and intramolecular interactions consist of a van der Waals repulsion-dispersion term calculated using the Lennard-Jones (6-12) potential, represented by the penultimate term in Eq.( 1). In the AMBER force field, the nonbonded interactions between atoms separated by exactly three bonds (1-4 van der Waals interactions) are reduced by a factor of 0.5 [39]. The Lennard Jones potential parameters for the interactions between unlike atoms were calculated by using the Lorentz-Berthelot mixing rules (quadratic and arithmetic rules for ij and σ ij parameters, respectively). The water molecules were represented with the TIP4P/2005 model [41].

This representation (Table 1) consists of four interaction sites: three of them are placed at the oxygen and hydrogen atom positions, respectively. The fourth site is placed at the bisector of the H-O-H angle.

The metal-ligand interactions are treated as being fully covalent and are modeled using the Morse potential [42].

Long-range electrostatic interactions

The last term in Eq.(1) corresponds to the electrostatic energy (E ELEC ) of the system. Consider N molecules in a volume V = L x L y L z with center of mass r i , each molecule contains n i charges q ia at position r ia . The electrostatics interactions handled with the Ewald sum method [START_REF] Allen | Computer Simulation of Liquids[END_REF]43] orthogonal axis that is overall charge neutral are given by the following contributions

U ELEC = 1 2 o V k =0 Q(h)S(h)S(-h) + 1 8π o i a j =i q ia b q jb erfc (αr iajb /r iajb ) - α 4π 3/2 o i a q 2 ia - 1 8π o i a b =a q ia q ib r iaib erf (αr iaib ) + 1 2 o V M 2 z ( 2
)
where erfc(x) is the complementary error function and erf(x) is the error function. α is chosen so that only pair interactions in the central cell need to be considered in evaluating the second term in Eq.( 2).

In the last term of Eq.( 2), M z is the net dipole moment of the simulation cell given by N i=1 q i r i . This contribution is the correction term from Yeh and Berkowitz [37] which results by the plane-wise summation method proposed by Smith [43]. Adding this term to the total energy amount to using a z-component force for each atom given by

F i,z = - q i o V M z (3)
The EW3DC method differs from the standard EW3D method only by the presence of this dipole correction. The functions S(h) and Q(h) are defined using the Eqs.( 4) and ( 5), respectively

S(h) = i a q ia exp(ih.r ia ) (4)
Q(h) = 1 h 2 exp(- h 2 4α 2 ) ( 5
)
where the reciprocal lattice vector h is defined as h = 2π(l/L x , m/L y , n/L z ) where l, m, n take values of 0,±1,±2,• The total Coulombic interactions for a two dimensional periodic system which is finite in the z-dimension can be calculated using the Hautman and Klein method [35]. The operational expressions of the total electrostatic interactions can be found elsewhere [35] for completeness.

The atomic charges on the grafted NTA-Cu(II) complex have been calculated from the density-functional theory (DFT) [44,45] (B3LYP) [46][47][48] with effective core potential (SDDALL) Gaussian basis sets using the GAUSSIAN 03 package [49] and the CHELPG [50] procedure as a grid-based method. The calculated charges are reported on Figure 1b.

Set-up conditions

The equations of motions are integrated using the Verlet Leapfrog algorithm scheme at T = 298 K with 2 fs as timestep. The simulations are performed in the NVT ensemble using the Hoover algorithm [51] with a coupling constant of 0.5 ps. SHAKE algorithm is used to preserve the rigidity of the water molecules and the OH and CH bonds in the NTA molecule. The cutoff for the Lennard-Jones contribution is fixed to 12 Å whereas it is equal to 18 Å for the real space of the electrostatic interactions. We also use a reciprocal space cutoff radius of 1.14 Å-1 . The box dimensions, the number of water molecules, the Ewald parameters are reported in Table 1 as a function of the system. The configurations are generated using the parallel version of the modified DL POLY MD package [START_REF] Dl | POLY is a parallel molecular dynamics simulation package developed at the Daresbury Laboratory Project for Computer Simulation under the auspices of the EPSRC for the Collaborative Computational Project for Computer Simulation of Condensed Phases (CCP5) and the Advanced Research Computing Group (ARCG) at the Daresbury Laboratory[END_REF] by using up to 8 processors at a time.

The parameters of the long-range electrostatics interactions evaluated by either the Ewald summation technique or the HKE techniques have been calculated in order to satisfy a relative error of 10 -6 . In the case of water embedded between two parallel plates, a run consisted of 200 ps of equilibration followed by a production phase of additional 600 ps. To allow a possible comparison between the different methods in a reasonable simulation time, an equilibration period of 200 ps and an acquisition phase of 300 ps were performed over the NTA-Cu(II) grafted surfaces. The corresponding CPU times are reported in Table 1 for comparison. Firstly, a system formed by 2340 water molecules embedded between two parallel graphite plates was studied. This latter has the particularity to present the total dipole moment of the simulation cell close to zero. It is then expected that the EW3D and EW3DC methods give the same results. A series of molecular dynamics simulations of liquid water between two parallel graphite plates with no defects with C-C distance of 1.42 Å have been performed at room temperature (T = 298 K) and fixed densities. The graphite surface is composed by one layer of 16×18 unit cells leading to 576 carbon atoms and to an area of (39.4×38.3 Å2 ). The carbon atoms of the graphite layer are treated as fixed on the initial rigid lattice.

We have also described the graphite as a set two or three flat sheets separated by 3.35 Å and stacked perpendicularly to the z axis. We have observed no difference on the structural and energetical properties of the system when the number of layers is increased. As a result, the graphite surface is modelled as only one flat sheet in the present work. The system is thus composed of two sheets of highly oriented pyrolytic graphite (HOPG) separated by a distance d of 50 Å. The intersurface spacing d is adjusted to yield the bulk density across the simulation box. When we apply the HKE method for the calculation of the long-range Coulombic interactions, the z-dimension of the box is fixed to 50 Å. We extend this box length dimension up to 190 Å in the case of the calculation of the electrostatic interactions using the EW3D and EW3DC methods to be in line with the recommendations of Spohr [START_REF] Spohr | [END_REF]. Molecular simulations have been already reported for this system in the literature [54,55]. However, our simulations differ by a greater separation distance and a systematic comparison between HKE, EW3D and EW3DC methods.

Figure 2a shows the molecular density profiles of water along the z-direction normal to the surfaces. We observe that the profiles can be considered as identical for the three methods and that a region of bulk

k B T αβ (z) = N i=1 H(z i ) m i (v i ) α .(v i ) β N i=1 H(z i ) (6)
where k B is Boltzmann's constant, m i and (v i ) α are the molecular mass and α-component of the velocity of the center-of-mass of water molecule i, respectively. H(z i ) is a top-hat function defined such that:

H(z i ) = 1 for z -∆z 2 < z i < z + ∆z 2 0 otherwise ( 7
)
where ∆z is the slab thickness. The analysis of the temperature profile establishes that the temperature is identical for each z and equal to the input temperature for EW3D and EW3DC. For HKE, we observe a slight decrease of the local temperature in the middle of the box (deviation of 1 % from the input T ) whereas there are slight increases of T near the surfaces with the same deviation. However, the average temperature calculated over the slabs is 300 ± 3 K for HKE and 300 ± 1 K for EW3D and EW3DC. We also check that the diagonal elements of the temperature tensor are equivalent showing an identical equipartition of the kinetic energy in the three directions for all the methods. The profiles of the off-diagonal elements of the temperature tensor fluctuate about a zero mean value, as expected for a system at thermal equilibrium.

We examine the structure of water in Figure 2c by calculating the oxygen-oxygen pair correlation function in 1.0 Å slab centered in the middle of the box for the three methods. We add for comparison the three dimensional oxygen-oxygen radial distribution function calculated in a bulk liquid simulation box. The fact that the plots resulting from the two-dimensionally periodic systems are identical with that of the three-dimensionally periodic bulk water system confirm that we have a region of bulk water in the middle of the cell. We also check that the structure of water is described in the same way with HKE, EW3D and EW3DC. To complete the comparison between the methods on this system, it is very interesting to compare the degree to which the water molecules are correlated over long distances. This can be done by

G k = M 2 -M M N µ 2 (8)
where M is the total dipole moment of the simulation box and µ 2 is the mean-squared molecular dipole moment. For the TIP4P/2005 model, µ is equal to 2.305 D. Figure 3a displays the cumulative time average of the factor G k over 600 ps. We check that G k is a slowly fluctuating property and that it requires significant simulation time to converge. We obtain final values of G k of 2.2, 2.2 and 2.1 for HKE, EW3D and EW3DC, respectively. For comparison, G k is equal to 3.1 for bulk water system. We show that the correlation of water molecules over long distances where long-range electrostatic continue on playing a major role does not depend on the method used (KHE, EW3D, EW3DC) in the case of water embedded between two graphite plates. We also check that the cumulative time average of the M M quantity (Figure 3b) converges to zero at long time. This quantity measures the absolute values of the fluctuations and is a good parameter to establish the convergence of M . We find that the total dipole moment M is equal to zero in agreement with the fact that there is not net polarisation in the simulation cell. This is confirmed by the profile of the z-component of M across the simulation cell shown in the inset of Figure 3b. The inset shows that the total M z is zero for each slab and for all the methods. Table 2 reports the total energy contributions of the simulation cell for the Lennard-Jones and electrostatic contributions.

We observe that the energy contributions can be considered as similar within the errors bars indicating a equal repartition of the different contributions for the three methodologies. The fact that the values of the Coulombic interactions match perfectly for the three methods also contributes to validate the EW3D methodologies for the calculation of these long range interactions.

We confirm the good agreement between the methods by the calculation of the parallel and perpendicular components of the diffusion coefficient of water by using the Einstein expression [START_REF] Mcquarrie | Statistical Mechanics[END_REF]. The values of the diffusion coefficient are listed in Table 2 as a function of the methodology. We observe that the different values of the components match satisfactory between HKE, EW3D and EW3DC. We also check that the components calculated from the Green-Kubo expression [START_REF] Mcquarrie | Statistical Mechanics[END_REF] are in line with those reported in This explains why we find values of diffusion coefficient close to that calculated in bulk water with the TIP4P/2005 model (2.4×10 -5 cm 2 s -1 ) and to the experimental self-diffusion coefficient of water [START_REF] Mills | [END_REF] at 298 K (2.3×10 -5 cm 2 s -1 ). This assertion can be checked through the analysis of the perpendicular component D ⊥ which does not differ significantly from the parallel component.

We now focus of a system consisting in grafting randomly 10 NTA-Cu(II) complexes onto a HOPG graphite surface formed by 576 atoms, the surface coverage being then 110 pmol cm -2 [7]. Each chain is covalently attached to a carbon of the graphite surface (Figure 1). A second graphite surface bounds the simulation in the z-direction. We take the route of distributing the total positive charge (+10 e) on each carbon atom with a charge of (10/576 e) of the second graphite surface to obtain a neutral system. This method avoids both the addition of counter-ions and the increase of the box size along the z axis. We have chosen this system model in order to have a system size allowing a comparison in a reasonable simulation time between the two-dimensional (HKE) and three dimensional Ewald methods. Whereas M z is equal to 0 D in the case of water embedded between two graphite plates, M z is about -0.6 D for the grafted system. for HKE, EW3D and EW3DC indicate that the structure of the water is correctly described within the system. Finally, all structural properties studied here show that the three-dimensional methods give results in perfect accordance with the two-dimensional HKE method. Let us consider results obtained for energetical and dynamical properties (Table 2). The values of the Lennard-Jones and Coulombic interactions deviate from less than 1% between the three methods. We also observe that the diffusion The components of the diffusion coefficients calculated in zone 1 defined by 25 Å < z < 50 Å are in fact very close to those of the water embedded between graphite planes. In the case of zone 2 defined by 0 Å < z < 25 Å, water molecules strongly interact both with the grafted molecules to give preferential hydrogen bonds and with the graphite surface involving preferential orientations. This leads to a slower diffusion of the water molecules in the interfacial region.

The fact that the EW3D and EW3DC method lead to the same results for the two studied systems can be explained by the dimensions of the simulation cell that makes an empty space length five times larger than the space of the primary occupied cell. Using such box dimensions, the methods are made very close. However, from a theoretical viewpoint, it is more rigorous to use the EW3DC method.

Structure and energy

Given that the EW3DC methodology correctly reproduces the results of a two-dimensional methods, we apply then it one a more sophisiticated system of grafted NTA-Cu(II) system. In fact, we decide to improve our methodology by taking into account the presence of counter-ions in solution to counterbalance the negative charges of the grafted complexes and to better model the experimental conditions. We then add to the simulation 10 mobile Na + cations to obtain a neutral system. We extend then the primary z-dimension from 50 Å to 80 Å. It amounts to extending within the EW3DC methodology the total z-dimension from to 190 Å to 320 Å with two empty spaces of 120 Å surrounding the primary cell. The different parameters for this simulation are given in Table 1. The simulation time is increased up to 1 ns for the production phase. The previous methodology consisting in charging the opposite graphite surface leads to the creation of dielectric field across the direction normal to the surfaces. In the case of the addition of counter-ions, we observe in Figure 5a, a quasi-zero profile of the local M z (z) across the z-direction from z = 60 Å. This indicates that the z-dimension is developed enough to yield a bulklike region from the grafted surface. To further investigate this behavior, we have plotted in part b of Figure 5 the angular density of probability n(cos θ) as a function of cos θ for three different zones A, B and C. θ is the angle between the molecular dipole moment vector of water molecules and the constant unit vector normal to the graphite surface. We observe an uniform distribution of cos θ for water in zone C as expected for distributions in bulklike region.

Zone A, which is close to the ends of grafted molecules, is characterized by a preferential orientation of water molecules with respect to the graphite surfaces: the hydrogen atoms of water molecules point to the surface. This effect is much less marked for the intermediate region B and disappears for zone A. To check the influence of the counter-ions on the electrical properties of the system, we calculate the charge density profile across the simulation box from the following expression

ρ q (z) = i q i ρ i (z) (9)
where q i is the charge of the atom i and ρ i (z) is the density number profile. The profile of the electric field E(z) along the z direction can be calculated by integrating the charge density profile according to

E(z) = 1 0 z -∞ ρ q (z )dz ( 10
)
where 0 is the vacuum permittivity. As expected from the nonisotropic distribution of the water dipole moment at the surface, the charge density profile shown in Figure 5c does not vanish close to the surface with strong fluctuations at distances closer than 10 Å from the graphite surface. From a distance greater than 20 Å, the charge density profile is zero. The analysis of the electric field profile in Figure 5d shows that the electric field reaches a constant zero value from 20 Å and confirms the fact that the presence of the counter-ions does not prevent the development of a fully bulklike region in the centre of the simulation box with an isotropic distribution of the water molecules. These electrical features could be recovered in the use of oppositely charged surfaces as long as there is sufficient separation distance between the surfaces.

We have shown that the use of the EW3DC methodology and the explicit modelling of the counter-ions produce consistent description of the electrical and thermodynamic properties of the grafted NTA-Cu(II) system. These checks are essential before an extensive study of the energetical and structural properties of the system. Our conclusions concerning the use of the EW3DC method in interfacial systems which are non-periodic in one of the three directions are in line with those drawn from previous works [37,[62][63][64]. We now focus on the description of the physical insights of the grafted NTA-Cu(II) system. A typical snapshot of the system with the water molecules is shown in Figure 6 in the case of a surface coverage of 110 pmol cm -2 . A detailed description of the modified graphite surface requires to represent the density profiles of atoms of the grafted molecules as a function of the distance from the surface. We represent in Figure 7a the density profiles of several atoms such as of the oxygen atoms of the carboxylate end groups, the copper and sodium atoms. We also report the profile of the carbon atom of the amide group close to the phenyl group covalently attached to the graphite surface. The profiles of the oxygen atoms of the carboxylate groups extend over a region of 13 Å (5 Å< z < 18 Å) indicating a high mobility of the terminal group of the grafted molecules. The profile of the amide carbon shows a well-defined narrow peak indicating that the first part of the molecule is rigid with very little degrees of freedom due to the covalent bond between the aryl group and the graphite surface. As expected from the modelling of bonds between the copper cation and some oxygen atoms of the terminal carboxylate groups, the profile of the copper ion is found to be similar to that of the carboxylate groups. The counter-ions present an uniform profile in the zone of 5 Å< z < 20 Å, in line with the fact that these regions are occupied by the polar groups of the grafted molecules. We also note that a very weak proportion of counter-ions can diffuse far from the surface but they do not affect the electrical properties of these zones. Figure 7b shows the density profiles of hydrogen bonds between the carboxylate end groups and water molecules and between the amide group and the water molecules. Hydrogen bonds can exist between the grafted molecules and water molecules within the region 4 Å< z < 20 Å. The inset of Figure 7b shows the total density of water-water hydrogen bonds and water-grafted molecules hydrogen bonds. We observe in the resulting profile that the value of 0.10 gcm -3 obtained in the bulklike zone is not recovered in the zone occupied by the grafted molecules. This means that the number of hydrogen bonds between the grafted molecules and water do not compensate the loss of hydrogen bonds between water molecules due to grafting of NTA-Cu(II) complexes.

F
Figure 8a displays the mean squared displacements of the molecular centres of water molecules in the parallel and perpendicular directions for the grafted zone (0 Å < z <20 Å) and the bulklike zone (35 Å < z < 80 Å). We give for comparison the mean squared displacements of the molecular centres of water in the case of isotropic bulk water. In the case of the grafted system, we take into account only the water molecules which have never left the region of interest during the calculation. The parallel component of the diffusion coefficient in the grafted and bulklike regions is equal to 1.0×10 -5 cm 2 s -1 and to 2.4× 10 -5 cm 2 s -1 , respectively, whereas the perpendicular component decreases from 2.0×10 -5 cm 2 s -1 in the bulk region to 0.8×10 -5 cm 2 s -1 in the region close to the surface. In the bulklike region, the value of the lateral component of the diffusion coefficient is equal to that of the isotropic diffusion coefficient for the TIP4P/2005 water model at 298 K (2.4×10 -5 cm 2 s -1 ) whereas the value of the perpendicular component is slightly weaker. This indicates that the diffusive behaviour of water molecules in the z direction ranging from 35 Å < z < 80 Å is weakly retarded due to the presence of the graphite surface.

However, Figure 8b shows that the normalized velocity autocorrelation functions of water molecules in the bulklike region match very well with that of the isotropic bulk water and tend slowly to zero. The lateral and perpendicular components of the velocity correlations functions calculated in the grafted zone contrast with those calculated in the bulk region by a more negative first peak and by the fact that these curves tend to zero with damped oscillations. These oscillations mean that the water molecules are trapped in a cage formed by neighbouring water molecules and NTA-Cu(II) grafted ligands. The fact that the first oscillations appear for negative values of the correlation functions indicate that the reversal of the velocities is correlated and that the cage effect has a certain lifetime. From 0.6 ps, the correlation functions oscillate around zero indicating that uncorrelated cages are formed successively.

The energetic description of the system requires to calculate different energy contributions between grafted molecules and between grafted molecules and the water molecules. We report in Figure 9a the average values of the dispersive contribution to the Lennard-Jones potential and of the electrostatic energy between a given molecule and the rest of the grafted molecules as a function of the molecule number.

The dispersive term modelled as -1 r 6 is used to model the van der Waals interactions whereas the -1 the direction normal to the surface; these molecules cannot relax to reduce the repulsion between their partial charges. This leads to an unfavorable electrostatic energy contribution. Part b of Figure 9 shows the dispersive and electrostatic energy contributions between each grafted NTA-Cu(II) ligand and the water molecules. We also give for completeness the number of water molecules interacting for each grafted molecule with the Lennard-Jones (r c = 12 Å) and electrostatic (r c = 18 Å) cutoffs. Interestingly, we observe that each NTA-Cu(II) molecules interacts approximately with the same number of water molecules.

This uniform distribution of the electrostatic and dispersive contributions shows that the hydration of the grafted molecule does not depend on the way of grafting the molecule for this value of surface coverage.

Conclusions

In view of the modelling of monolayers of copper complexes covalently attached to graphite surface in water, we have compared one two-dimensional (HKE) and two three-dimensional methods (EW3D) and (EW3DC) for the calculation of the long range electrostatic interactions. For a slab geometry, it is necessary to extend enough the primary simulation cell in the z-axis so that there is no artificial influence between the periodic images in this direction to use a standard three-dimensional Ewald summation.

From a theoretical viewpoint, EW3DC method must be applied when there is a net polarization in the simulation box because it adds a correction term required to remove the leading dipole of the primary cell. We have firstly checked on a typical system formed of water molecules embedded between two graphite surfaces that the density profiles, hydrogen bonds, energy contributions and components of the diffusion coefficient are independent on the method used. We note that this observation is in line with the recommendations given by Yeh and Berkowitz [37] i.e that the current EW3DC method is sufficient in terms of accuracy and efficiently.

Finally, using the same methodology, we have also modelled a grafted system modified in such a way to present a larger z-dimension and to contain counter-ions to counterbalance the negative charges of the grafted molecules. The comparison between the different methods were necessary to check some fundamental aspects concerning the thermodynamic equilibria and the electrical properties of the grafted system. The fact that the methods compare very well is also very important from a computational viewpoint because the three dimensional methods are 20 times faster than two-dimensional method.

We have then reported a detailed molecular description of the grafted monolayers of copper complexes.

We have shown that the carboxylate end groups maintain degrees of freedom and can form hydrogen bonds with the water molecules. The first part of the grafted molecules is strongly constraint with very narrow density profiles due to the covalent bond between the aryl group and the grafted surface. We have shown that the hydrogen bond formed between the grafted molecules and the water cannot compensate the loss of hydrogen bonds between water molecules due to the presence of grafted molecules. The calculation of the diffusion coefficient shows that a bulklike region can also be identified with a slight retard of the diffusion along the normal to the surface. The analysis of the velocity correlation functions in the grafted zone exhibits damped oscillations suggesting a cage effect of water molecules in this zone. The energetic description of the system establishes that the hydration of the grafted molecules is identical for each grafted molecules with a non-homogeneous distribution of the molecule-molecule energy contributions.

Given that the methodology and the system size are now well adapted, we plan to study the impact of some parameters (nature of the surface (gold or graphite), length of the grafted molecule, surface coverage, pH) on the properties of the grafted monolayers.

F

1

 1 Comparison between HKE, EW3D and EW3DC methods

  liquid of approximately 25 Å is formed in the middle of the cell. Part b of Figure2shows the molecular temperature profiles of the local temperature T defined by 1/3(T xx + T yy + T zz ) components. The T αβ component along the z-direction is calculated via Eq.(6)

  finite-system g-Kirkwood factor[56] G k from the following equation

Figure 4a presents the

  Figure4apresents the density profiles of the water molecules and of the grafted NTA-Cu(II) complexes for the three methodologies. Interestingly, we observe that the different methods yield the same density profiles. The profile of the number of hydrogen bond per water molecules is reported in Figure4b. Let us recall that two water molecules are identified as being hydrogen bonded only if the oxygen-oxygen distance is less than 3.5 Å and simultaneously the OH...O distance is less than 2.45 Å. We check that respecting these two conditions leads to the oxygen-oxygen-hydrogen angle is less than 30 o[59][60][61]. The profile shows a decrease in the number of hydrogen bonds until z = 16 Å. Up to 16 Å, the number of hydrogen bonds increases again. The fact that the profiles of the number of hydrogen bond per water molecules and the oxygen-oxygen pair correlation function in the bulklike (Figure4 c) are very close

  two different locations match satisfactory between HKE, EW3D and EW3DC.

r 12 represents

 12 the repulsive interactions. We observe an inhomogeneous distribution of the dispersive term as a function of the grafted molecule number whereas the distribution of the electrostatic contribution is much more homogeneous. The dispersive and repulsive terms of the Lennard-Jones potential give important contributions at short distances; the electrostatic contribution in ( 1 r ) continue to yield contributions at larger distances. It means that the repulsive and dispersive terms are very sensitive on the way of grafting the molecules and of the local environment whereas the electrostatic contribution gives rather an average contribution for the grafted surface. Grafting NTA-Cu(II) ligands amounts to constraining molecules along

Figure 1 :

 1 Figure 1: a) Modified Nitrilotriacetic acid (NTA) complexed by Cu(II) cation b) atomic charges calculated from DFT/B3LYP method and c) snapshot of MD configuration of 10 grafted NTA-Cu(II) complexes onto the HOPG surface without the water molecules. Oxygen (red), Carbon (black), Nitrogen (blue), Cu(yellow), Hydrogen (white).

Figure 2 :

 2 Figure 2: a) Water molecular density profiles along the z-direction for HKE (open circles), EW3D (solid lines) and EW3DC (dashed lines). b) Temperature profiles of T (z) calculated from 1/3(T xx (z) + T yy (z) + T zz (z)) for the three methods. c) O-O pair correlation functions for water in 1.0 Å slab which is centered in the middle of the simulation cell. The diamond points correspond to the three dimensional radial distribution function calculated in a bulk water simulation cell.

Figure 3 :

 3 Figure 3: a) Cumulative time average of the factor G k of TIP4P/2005 water model at 298 K as a function of the different methods as indicated in the legend. b) Convergence of the square of the average total dipole moment for water. The inset shows the M z (z) profile along the z-direction for the three methods. M z is given in Debye ( 1 D = 10 -18 esu cm).

Figure 4 :

 4 Figure 4: a) Molecular density profiles of water and atomic density profiles of the grafted molecules along the z-direction for HKE (open circles), EW3D (solid lines) and EW3DC (dashed lines). The density profiles of the grafted molecules are scaled down by 2 for clarity. b) Profiles of the number of hydrogen bonds between water molecules along the direction normal to the surfaces. c) O-O pair correlation functions for water in 1.0 Å slab located at z = 30 Å.

Figure 5 :Figure 6 :

 56 Figure 5: a), z-component of the total dipole moment of the slab centered at z along the z-direction; b) angular density profile of probability n(cos θ) as a function of the cosinus of θ for three different zones A(20 Å < z < 30 Å), B(40 Å < z < 50 Å) and C(60 Å < z <70 Å). θ is formed by the dipole moment of water molecules and the z-axis direction, represented by a unit vector perpendicular to the graphite surface; c) charge density and d) electric field profiles along the direction normal to the surface.

Figure 7 :

 7 Figure 7: a) Density profiles of different atom types as a function of z and b) density profiles of the hydrogen bond between the grafted molecule and the water molecules.

Figure 8 :

 8 Figure 8: a) Parallel and perpendicular components of the mean-squared displacements and b) velocity autocorrelations functions of water molecules for two different zones: a grafted zone defined by 0 Å < z <20 Å and a bulklike zone defined by 35 Å < z < 80 Å.

Figure 9 :

 9 Figure 9: a) Dispersive and electrostatic energy contributions between the grafted molecule of interest and the remaining 9 grafted molecules; b) Moleculewater energy contributions and number of water molecules interacting with the considered grafted molecule within two spheres of radii 12 Å and 18 Å.

  Figure 1.

  • • ± ∞. The reciprocal space sum is truncated at an ellipsoidal boundary at the vector |h max |. The values of the components of the vector |h max | are given in Table 1 for each system.

	May 28, 2008 14:23 Molecular Physics pub1˙v4
	7
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 2

 2 within the statistical uncertainties. In the case of separation distance of 50 Å, a large percentage of water molecule is far from the surfaces suggesting that the influence of the surfaces is reduced.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1 .

 1 Dimensions (Lx, Ly, Lz) ( Å) of the simulation box, reciprocal lattice vectors (|h max The CPU times (hour) are reported for a production phase of 600 ps for water embedded between two parallel plates and a production phase of 300 ps and for the grafted NTA-Cu(II) ligands. The reported CPU time represents the real time of the job carried out over 8 processors.

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 2 .

 2 Lennard-Jones (U LJ ) and Coulombic (U ELEC ) energy contributions (kJ mol -1 ), overall coefficient diffusion D, parallel (D // ) and perpendicular (D ⊥ ) components of D (10 -5 cm 2 s -1 ). The diffusion coefficients are calculated over 60 time origins separated by 0.5 ps. The calculation is performed over a total time of 50 ps. In the case of the grafted NTA-Cu(II) complexes, the diffusion coefficients are calculated at two different locations (zone 1/ zone 2). The zone 1 and zone 2 are defined by 0 Å < z < 25 Å and 25 Å < z < 50 Å, respectively. The experimental self-diffusion coefficient (10 -5 cm 2 s -1 ) of water[START_REF] Mills | [END_REF] at 298 K is reported for comparison.

	Method	U LJ	U ELEC	D	D //	D ⊥
			Water between graphite plates
	HKE	1.736×10 4	-1.295×10 5	2.1	2.3	1.8
	EW3D 1.722×10 4	-1.296×10 5	2.5	2.7	1.9
	EW3DC 1.722×10 4	-1.296×10 5	2.6	3.0	1.8
	Grafted NTA-Cu(II) complexes onto graphite plates with charges distributed on the opposite surface
	HKE	1.802×10 4	-1.450×10 5	0.9/2.4 0.8/3.0	0.9/1.9
	EW3D 1.793×10 4	-1.460×10 5	1.0/2.6 1.0/3.1	0.9/1.9
	EW3DC 1.790×10 4	-1.460×10 5	1.0/2.5 1.1/2.8	0.8/2.2
	Experimental diffusion coefficient of water [58]	2.3


Anna Grochola

Wlodzimierz Jastrzebski

Pawel Kowalczyk

Spectroscopic study of the 6 1 Π u state in Li 2

Keywords: laser spectroscopy, alkali dimers, electronic states, molecular constants, potential energy curves laser spectroscopy, alkali dimers, electronic states, molecular constants, potential energy curves

come

Polarisation labelling spectroscopy (PLS), a technique first demonstrated by

Schawlow and coworkers in 1976 [1] has been widely used in studies of small molecules (see e.g. Wang and Xia [START_REF] Wang | Molecular and Laser Spectroscopy[END_REF] and references therein). In particular, the selectivity of double resonance inherent to this method allows complex molecular spectra to be resolved and assigned. For several years we have employed this powerful technique to study excited electronic states of alkali dimers, and recently have focused on the lightest of them, Li 2 . Our previous experiments furnished molecular constants and accurate potential energy curves of three excited states of lithium dimer, C 1 Π u , 2 1 Σ + u [3] and 5 1 Π u . [4] In the present work we extend our study of Li 2 to higher excitation energies (between 36800 and 39500 cm -1 above the minimum of the ground state potential well) and characterise experimentally a newly observed electronic state of 1 Π u symmetry.

Experimental details

The experimental procedure has been described in our previous papers [3,5] so only an outline of the apparatus will be given here. The Li 2 molecules were produced by heating solid lithium (2 cm 3 ) of natural isotopic composition in a linear heat-pipe oven, 80 cm long and with 3 cm inner diameter. About 4 mbar of helium buffer gas was admitted to the oven before heating it to approximately 1050 K. Under these conditions the oven did not operate in a heat pipe regime [6] but instead contained a mixture of helium and metal vapour in the central zone.

As a tuneable pump light we employed frequency doubled output of a parametric oscillator/amplifier system (OPO/OPA Sunlite EX, Continuum) pumped with the third harmonic of an injection seeded Nd:YAG laser (Powerlite 8000). The laser system delivered pulses with a spectral width of about 0. wavelength calibration was effected by recording the optogalvanic spectrum of argon [START_REF] Striganov | Tablicy Spektralnych Linii Atomov i Ionov[END_REF] and transmission peaks of a Fabry-Pérot etalon with FSR=1 cm -1 (both signals observed at fundamental frequency of tuneable OPO/OPA radiation). The rotationally resolved spectrum of the investigated 1 Π u ← X 1 Σ + g transition in Li 2 was recorded in the region 255 -280 nm.

We estimate the accuracy in absolute frequency measurements to ±0.1 cm -1 .

As a probe (labelling) laser we used a home built dye laser (Shoshan-type design, [START_REF] Shoshan | [END_REF] DCM dye) pumped synchronously with the same Nd:YAG laser. The dye laser operated on chosen wavelengths between 630 and 650 nm, controlled with HighFinesse WS-5 wavemeter.

Alternatively, fixed frequency radiation of the Ar + -Kr + laser (Carl Zeiss ILM120, λ = 476.5, 488.0, 496.5, 501.7 or 647.1 nm) was utilized. At all the employed wavelengths the probe beam excited known transitions in the

A 1 Σ + u ← X 1 Σ + g or B 1 Π u ← X 1 Σ + g band systems of
Li 2 , [9,10] thus providing labelling of individual levels in the molecular ground state which is essential for the PLS technique.

The pump and probe beams co-propagated at a small angle, intersecting in the molecular sample. During a frequency scan of the pump light, the intensity of the probe beam was monitored through a crossed polariser-analyser system (Glan-Thompson prisms placed either side of the heat-pipe) so that double resonance signals appeared on a black background.

Experiments were performed with linearly polarised pump light, with the plane of polarisation inclined by 45° against that of the probe light. Such an arrangement was chosen to enhance intensities of Q lines in the observed spectra [11] (vide infra). This part of the probe light that was transmitted through the analyser was directed onto the cathode of a photomultiplier tube coupled to the boxcar averager (Stanford SR250). The resulting signals were stored in a computer which simultaneously recorded the reference optogalvanic spectrum of argon together with the frequency markers from the etalon and also controlled tuning of the OPO/OPA system.

Measurements and analysis

In the spectral region under investigation we observed a band system corresponding to a transition from the ground state of Li 2 , X 1 Σ + g . As the progressions display strong Q lines, the excited state must be 1 Π u . Even a cursory inspection of the spectra shows perturbation of the upper state, resulting in irregular patterns of P, Q and R lines; often only Q lines could be reliably assigned. In the preliminary analysis we found that positions of P and R lines corresponding to excitation of e parity levels in the 1 Π u state behaved erratically (deviations up to 3 cm -1 ) while Q lines (leading to excitation of f parity levels) formed regular progressions. This observation suggests that the investigated state is perturbed by a neighbouring state of 1 Σ + u symmetry, containing only levels of e parity. However, such a 1 Σ + u state could not be directly seen in the spectra, presumably because of too low Franck-Condon factors for excitation from the ground molecular state of Li 2 . This made any attempt to deperturb the spectra pointless and we had to confine our analysis to the observed Q lines, i.e.

to f parity levels of the investigated 1 Π u state. Note that we have encountered similar problems in the previous experiments on high electronic states in Li 2 , [3] Na 2 [12] and K 2 . [13,14] A total of 124 Q lines were assigned to the investigated 1 Π u ← X 1 Σ + g system, spanning the range of vibrational and rotational quantum numbers in the upper state υ = 0 -11, J = 2 -47 (positions of all the recorded lines are available from the authors upon request). Although a slight breakdown of the Born-Oppenheimer approximation in Li 2 is well known, [3,4] observation of the isotope effect in 7 Li 2 and 6 Li 7 Li molecules allowed to establish unambiguously the absolute vibrational numbering in the upper state. However, even when taking into account an isotopic shift of the dissociation limit of the observed state (see [15] for the relevant atomic data) the lines of the mixed isotopic form deviated from the expected positions by more than twice the experimental uncertainty. Therefore they were excluded isotopologue. The measured wave numbers of lines were fitted to differences in term values, ν = T(υ', J') -T(υ", J"), of the upper and lower levels, which were expressed by the standard Dunham formula

∑ Λ - + + + = n m n m mn e J J υ Y T J υ T , 2 ] ) 1 ( [ ) 5 . 0 ( ) , ( . (1)
The symbols above have their usual meaning, with Λ 2 = 0 for the ground X 1 Σ + g state and 1 for the upper 1 Π u state. As the coefficients for the ground state are well known, [16] they were treated as fixed constants in the fits. A degree of the polynomial in (1) was chosen by trial and error to obtain a set of statistically significant Dunham coefficients (i.e. with errors less than 5% of each value) and, simultaneously, a good quality of the fit. The final, best-fit values of the constants are listed in Table 1. They reproduce the observed line positions with the root mean square deviation of 0.12 cm -1 , which is about the precision of our measurements. The number of digits required to represent each constant with no loss of precision in reproducing the experimental data was minimised using the sequential rounding and refitting procedure suggested by Le Roy. [17] Once the coefficients in the Dunham expansion have been determined, we applied the Rydberg-Klein-Rees (RKR) method to obtain the rotationless potential energy curve of the investigated state. Table 2 lists the vibrational energy levels and the corresponding classical turning points.

In the absence of direct experimental observations, the identity of the 1 Π u state under consideration and its assignment to an atomic asymptote can be done by comparison with theoretical calculations. The only computation covering excitation energies high enough originates from the unpublished thesis by Spies, [18] who extended earlier valence CI calculations with a core polarisation potential by Schmidt-Mink et al. [19]. Figure 1 displays the experimental potential curve together with theoretical curves of 1 Π u states placed in a similar energy range. As the 5 (2)

By using the values of D e (X 1 Σ + g ) = 8516.61 ± 0.08 cm -1 [20] and ∆E(4d -2s) = 36623.30 ± 0.01 cm -1 , [21] we arrive at value of D e (6 1 Π u ) = 8276.7 ± 0.1 cm -1 , where the errors of all three summands contribute to the uncertainty of the result.

The calculations by Spies allow also to assign the state responsible for perturbation of e parity levels in 6 1 Π u as the 5 1 Σ + u state (Figure 1). Because of its broad potential energy curve Franck-Condon factors for the 5 1 Σ + u ← X 1 Σ + g transition are expected to be low which explains why the corresponding band system did not show up in the experimental spectra. Dunham coefficients (in cm -1 ) for the 6 1 Π u state of 7 Li 2 representing energies of rovibrational levels of f parity in the range υ = 0 -11, J = 2 -47, compared with theoretical values. [17] In the bottom the equilibrium distance and dissociation energy are also given. --------------------------------------------------T

F

Constant Value Theory

-------------------------------------------------a The number in brackets that follows a coefficient is one standard deviation in units of the last digit. ------------------------------------------------

  of molecular constants which is based solely on the results for the7 Li 2

  1 Π u state has been identified at lower energies and characterised Π u is the only reasonable assignment for the state observed here. A comparison of the main molecular constants resulting from the experiment and theoretical calculations is included in Table1and shows fair agreement. It is interesting to note that the calculations underestimate the equilibrium distance, similarly as for the closelying 5 1 Π u state described in reference[4]. Spies assigns the adiabatic asymptote of the 6 1 Π u state as 2s+4d. The dissociation energy of the investigated state can be obtained then from combination of the measured T e value, the well depth of the ground state D e (X 1 Σ

	F experimentally in our previous work,[4] 6 1 + o r separation of atomic energy levels of lithium, as	g ) and the
	e D	(	6 1	Π	P ( ) e u X D =	1 Σ	g	)	+	∆	E	(	4	d	-	s 2	)	-	e T	(	6 1	Π	u	)
					e
					e r
									R
											e
													v i e
																		w
																					O n l
																								y
																									7

+

.

Table 1

 1

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 2

 2 Rotationless RKR potential energy curve for the 6 1 Π u state of7 Li 2 molecule.

	υ	G υ (cm -1 )	R -(Å)	R + (Å)
	-

  The experimental RKR potential energy curve of the 6 1 Π u state (full circles connected by a solid line) compared with theoretical potential curves for the 5 -7 1 Π u states calculated by Spies[17] (open symbols connected by dashed lines). Theoretical potential of the 5 1 Σ + u state, responsible for perturbation of e parity levels in the 6 1 Π u state, is also drawn (stars and dash-dot line).

	FIGURE CAPTION Figure 1. Page 12 of 13 Molecular Physics Page 13 of 13
	0.0 126.392 o 0 F r	2.9809	3.1639 a	3.3715	-
	1 2 3	376.338 622.716 865.589 F o r P e e 2.8599 2.7811 2.7189 r		3.5404 3.6650 3.7709
	4 5	1105.017 1341.062 P R 2.6656 2.6180 e e r R e v i 3.8656 3.9526 e w e v i e w O n l y
						O n l
						y

URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

This work received partial funding from the Polish Ministry of Science and Higher Education through Grant No. N202 103 31/0753. We thank Z. Jędrzejewski-Szmek for experimental assistance.