Stefanie Becke

Stefanie Haller

Michael A Osipov

Frank Giesselmann

Correlation between the molecular chirality index and the spontaneous polarization in series of smectic C* liquid crystals

Keywords: ferroelectric liquid crystals, chirality index, molecular structure ferroelectric liquid crystals, chirality index, molecular structure

come

Correlation between the molecular chirality index and

the spontaneous polarization in series of smectic C* liquid crystals

INTRODUCTION

Ferroelectric liquid crystals continue to attract significant attention because of the unique origin of the spontaneous polarization, and because of various potential applications of these materials which are already used, for example, in microdysplays.. It is known since the pioneering work of R.B.Meyer [1] that the molecular chirality is a necessary condition for the appearance of the spontaneous polarization in the smectic C* liquid crystal phase. The spontaneous polarization in this system can thus be considered as an experimentally observable macroscopic manifestation of chirality which should be related to microscopic chirality of the constituent molecules. It is well known that even small changes in the molecular structure may result in dramatic changes of both magnitude and sign of the spontaneous polarization. This raises a general problem of the relationship between the macroscopic spontaneous polarization and microscopic molecular parameters which represent its chirality.

Recently a detailed molecular-statistical theory of ferroelectric ordering in smectics C* has been developed [2,3] and used to obtain approximate expressions for the spontaneous polarization in terms of molecular model parameters. The theory enables one to describe various properties of ferroelectric liquid crystals and to understand the general role of molecular chirality. On the other hand, a molecular theory is inevitably based on some approximations, and explicit expressions for the polarization can only be obtained using very simple models for chiral molecules. Thus the theory can hardly be used to explain how small changes of a particular molecular structure effect the spontaneous polarization.

The relationship between the spontaneous polarization and the molecular chirality can, in principle, be studied experimentally provided their exists a "measure" of molecular chirality which can be calculated using the molecular structure and composition. Chirality measures of the two different types have been proposed to characterize the chirality of idealized objects using mainly the geometrical approach. All chirality measures vanish identically for nonchiral objects, but the measures of the first type are scalars, and thus they do not distinguish enantiomers. The primary example of a scalar chirality measure is the Hausdorff measure proposed by Mislow et.al. [4]. The Hausdorff measure is defined as the Hausdorff distance [4] between the geometrical object and its mirror image minimized over all mutual translations and orientations of the two objects. The advantage of the Hausdorff measure is related to the possibility to distinguish between chiral and nonchiral objects on a rigorous generally more difficult to establish a relationship between a scalar chirality measure and macroscopic physical chiral properties which are expressed as pseudoscalars (such as optical activity, helical twisting power in nematic liquid crystals, etc.).

A relatively general chirality measure of the second type (the pseudoscalar geometrical chirality index) has been proposed in ref. [5]. It has later been scaled in ref. [6] in order to enable a comparison of chirality for molecules of different size. The derivation of the geometrical chirality index has been inspired by the Kirkwood theory of molecular optical activity [7]. Another pseudoscalar chirality index, which characterizes the chirality of the effective molecular shape, and which has been designed to describe the helical twisting power of chiral dopants has been proposed by Nordio and coworkers [8][9][10]. Pseudoscalar indices are specially designed to characterize a relationship between the molecular chirality and measurable pseudoscalar physical properties. It should be noted, however, that the disadvantage of all pseudoscalar chirality measures is related to the problem of chiral zeros [4]. It can be shown [4,11] that any chiral set of points can be continuously transformed into its mirror image without passing a nonchiral state. It follows then that any pseudoscalar chirality measure of the set must vanish for at least one chiral object, and thus a pseudoscalar measure cannot be universal. On the other hand, pseudoscalar measures are more closely related to chiral physical properties and may be used to study the effects of relatively small changes of the molecular structure (including those occuring in series of chiral liquid crystals) on macroscopic chiral properties.

Both the scaled chirality index G OS [5,6] and the Nordio index Q can readily be calculated numerically for molecules of arbitrary size, and are widely used as chirality measures showing a good correlation with the helical twisting power of various chiral dopants in nematic liquid crystal solvents [6,[8][9][10]12]. Very recently the scaled chirality index has also been successfully used to predict chirality of proteins which opens a new field for its application [13]. The Hausdorff chirality measure has also been applied to real molecules [14,15] but the application technique is still to be optimized.

The scaled chirality index has also been used to study a relationship between the molecular chirality and the spontaneous polarization in smectic s C* [16], but only a weak correlation has been found. One notes that there are two main reasons why the correlation has been found to be weak. Firstly, the spontaneous polarization is known to be proportional to the magnitude of the transverse molecular dipole which may vary significantly. Thus one expects a better correlation between the spontaneous polarization and the product of the chirality index and the transverse dipole, which has not been taken into consideration in [16]. Secondly, only five smectic C* mesogens with significantly different molecular structure have been considered in [16]. In such an ensemble the correlation between the polarization and the molecular chirality can also be shadowed by the orientational order parameter of the transverse dipole which is expected to be different for molecules of very different structure.

In this paper we attempt to find out whether the scaled chirality index is a useful measure to describe the variation of the spontaneous polarization in series of ferroelectric smectic C* liquid crystals. We selected four series of Smectic C* mesogens where systematic minor changes of the molecular structure within each series result in significant variation of the spontaneous polarization. For all mesogens the scaled chirality index has been calculated numerically using atomic coordinates of the minimum energy conformation. The paper is arranged as follows. In Section 2 the general relationship between molecular chirality and the spontaneous polarization in the Smectic C * phase is discussed, and in Section 3 we consider correlations between the experimental values of the spontaneous polarization and the chirality index for constituent molecules for five series of ferroelectric liquid crystals.

II. SPONTANEOUS POLARIZATION AND MOLECULAR CHIRALITY IN

SMECTIC C* LIQUID CRYSTALS.

In the smectic C* phase the director, specified by the unit vector n which is tilted with respect to the smectic layer normal k by the angle Θ, and the spontaneous polarization P s is parallel to the C 2 symmetry axis of the layer which is perpendicular to the tilt plane (see Fig. 1).

According to the general phenomenological theory of ferroelectric liquid crystals [START_REF] Pikin | Ferroelectric Liquid Crystals. Principles, Properties and Applications[END_REF] the spontaneous polarization can be expressed as:

P s = χµ p (n • k)[n × k], (1)
where χ is the transverse dielectric susceptibility, (n • k) = cos Θ, and the coefficient µ p vanishes if the molecules are nonchiral. One notes that in eq.( 1) polarization P s is a polar vector while the cross product [n × k] is a pseudovector which has different transformation properties. Taking into account that the susceptibility χ is a scalar, one concludes that the equality sign in eq.( 1) is only possible if the coefficient µ p is a pseudoscalar, i.e. a quantity which changes sign under space inversion, and which possesses opposite signs for any two enantiomeric systems.

A more detailed general expression for the spontaneous polarization can be obtained combining eq.( 1) with a statistical mechanical expression for the polarization. Indeed, in the molecular-statistical theory the polarization in the Smectic C* phase can be defined as an average transverse molecular dipole per unit volume:

P s = ρ < µ ⊥ >, ( 2
)
where ρ is the molecular number density, µ ⊥ is the molecular transverse dipole and < ... > denotes the statistical average. The orientation of a rigid anisotropic molecule can be specified by three mutually orthogonal unit vectors a, b and c where a is the primary molecular axis which is ordered along the director n on average. The unit vectors b and c are in the direction of the two orthogonal short molecular axes. One of these short axes, i.e. the axis c, can be defined as the axis which orders along the C 2 symmetry axis in the Smectic C* phase(see Fig. 1), i.e. in the direction of the spontaneous polarization. The particular orientation of this axis in the molecular frame depends on the molecular structure, but the transverse molecular dipole is not necessarily parallel to it. In the general case there is an angle β between the transverse dipole and the is orthogonal to c and thus, on average, it is parallel to the tilt plane which is perpendicular to the C 2 axis. As a result the statistical average < b >= 0 because there is no polar ordering in the tilt plane. Substituting the expression for the transverse dipole into eq.( 2) one obtains

P s = ρµ ⊥ cos β < c > . ( 3
)
comparing eq.( 3) with eq.( 1) one concludes that the vector < c > must be proportional to

[n × k]. Taking also into account that (n • k)[n × k] = (1/2) sin(2Θ)w
, where w is the unit pseudovector in the direction of the spontaneous polarization, one obtains:

P s = 1 2 ρµ ⊥ cos βC sin(2Θ)w, (4)
where the pseudoscalar coefficient C =< (c • w) >. In eq.( 4) the pseudoscalar coefficient C depends on molecular chirality and also contains contributions from all intermolecular interactions which determine the spontaneous polarization. It has been shown in [2] that C = ∆Γ where ∆ s a molecular chirality coefficient and Γ is the biaxial orientational order parameter which mainly characterizes the orientational order of the short molecular axes (i.e. the transverse dipole) [2].

Finally the ratio of the polarization and tilt can be expressed as:

P 0 = P s sin 2Θ/2 = ρµ ⊥ cos β∆Γ (5)
where the pseudoscalar ∆ represents the molecular chirality, and one may assume that it is directly related to the pseudoscalar chirality index discussed above. In the next section we investigate a correlation between the polarization/tilt ratio P 0 and the scaled molecular chirality index G OS . It should be noted, however, that such a correlation can only be detected if the angle β and the orientational order parameter Γ do not vary too much within the selected group of ferroelectric SmC * materials.

The effect of the angle β between the transverse dipole and the "main" short axis c on the spontaneous polarization is illustrated in Fig. 1. For example, if the axis c is parallel to the "molecular plane" of the biaxial molecule (see Fig. 1), the molecular plane appears, on average, to be perpendicular to the tilt plane. In this case the maximum polarization is achieved if the transverse dipole is also parallel to the molecular plane. In contrast, the polarization is very low if the dipole is normal to the plane. Thus, in the general case, the spontaneous polarization strongly depends not only on the magnitude of the transverse molecular dipole but also on its orientation in the molecular frame. The orientational order parameter Γ also depends on the molecular structure. It is well known that some smectic C* mesogens with comparable transverse dipoles exhibit spontaneous polarization which may differ by an order of magnitude or more.

Thus one arrives at the conclusion that a correlation between the spontaneous polarization and molecular chirality can be best found in special series of smectics C* in which small changes of the details of the molecular structure result in a significant variation of the spontaneous polarization while the parameters β and Γ are not changing strongly.

III. CORRELATION BETWEEN THE SPONTANEOUS POLARIZATION AND THE CHIRALITY INDEX IN FIVE SERIES OF SMECTICS C*

The so-called scaled chirality index G OS is given by the following expression [6]:

G OS = 1 3 4! N 4 N i,j,k,l=1 w i w j w k w l [(r ij × r kl ) • r il ](r ij • r jk )(r jk • r kl ) (r ij r jk r kl ) 2 r il , ( 6
)
where the summation is performed over all sets of four atoms i, j, k, l with coordinates r α , α = i, j, k, l in the molecule, and where r ij = r i -r j and N is the total number of atoms in the molecule. Here w i , w j , ... are the weight factors that allow for an adjustment of the chirality index in order to describe different chirality effects by inserting e.g. atomic weights, partial atomic changes or van der Waals radii. In the simplest case the factors are simply set to be equal to one which yields a purely geometrical chirality index. In this paper such a geometrical index has been calculated for real molecules using the atomic coordinates of the minimum energy conformation which is obtained by optimization of the molecular geometry by means of semi-empirical quantum chemical calculations using MOPAC and the AM1

Hamiltonian.

The data for the spontaneous polarization and the tilt angle for the first four series of ferroelectric smectics C* have been taken from [START_REF] Furukawa | [END_REF]. In each of these series (see Tables 1,2) As an example we present the molecular models for minimum energy conformations for all molecules in series D (see Fig. 2). One can readily see that the changes in the molecular structure are mainly restricted to some rotations of the substitution groups in the chiral chain and to the changes of the dihedral angle of the adjacent benzene ring.

F
The scaled chirality index given by eq.( 6) has been calculated by a MatLab routine (for w i = w j .. = 1) for minimum energy conformations of all molecules in series A-D.

The transverse dipole moments of these molecules have been obtained using the following procedure. Firstly the molecular inertia tensor has been calculated and diagonalized, and Then the transverse components of the total dipole moment have been obtained using the MOPAC/AM1 calculations for the minimum energy conformation. The results of these calculations are presented in Table 3. The correlation between the molecular chirality index multiplied by the molecular transverse dipole and the experimentally observed saturated ratio of the spontaneous polarization FIG. 3: Linear correlation between the product of the chirality index and the transverse molecular dipole and the spontaneous polarization for mesogens taken from [START_REF] Furukawa | [END_REF].

One can readily see that the correlation with the product involving the transverse molecular dipole, presented in Fig. 3, is much better then the one involving the full molecular dipole which is illustrated in Fig. 4. The correlation between the reduced polarization and the scaled chirality index illustrated in Fig. 3 is reasonable except for the compound A-H from the series A. For this compound an abnormally large value of the scaled chirality index has been obtained which doe not go in line with values for similar molecules, and which also possesses the wrong sign. The reason for this may be related to the fact that evaluation of small molecular dipoles with a semiempirical method is generally not reliable.. Apart from this there is a good correlation within each series of ferroelectric smectics C* which and the spontaneous polarization for mesogens taken from [START_REF] Furukawa | [END_REF].

experimentally for one member of the series. One notes, however, that the slopes of the regression lines are different for different series, which means that one cannot predict the value of the spontaneous polarization using the data for another series. The difference in slopes can be explained taking into account that the molecules of different series possess significantly different structure. As a result, the angle β between the transverse dipole and the primary molecular short axis and the orientational order parameter Γ for short molecular axis, which enter eq.( 5) for the spontaneous polarization, may differ significantly for different series of smectic C* mesogens because they are expected to depend strongly on the molecular structure. terminal chain, derived from natural (S)-amino acids, varies (see Fig. 5).

The variation of the chiral chain leads to a strong negative polarization in the case of LC2

, and a moderate negative polarization in LC1. Further variation of the chiral chain results in a sign inversion of the spontaneous polarization, and the mesogen LC3 possesses a moderate positive polarization. The correlation between the product of the scaled chirality index multiplied by the molecular dipole and the experimental value of the saturated spontaneous polarization, taken from ref. [START_REF] Kelly | Handbook of Liquid Crystals[END_REF] is illustrated in Fig. 3. One can readily see that there is a good correlation with both magnitude and sign of the spontaneous polarization. Moreover, the calculated chirality index undergoes a sign inversion together with the experimentally observed spontaneous polarization. This particular example illustrates one of the practical advantages of a pseudoscalar chirality measure as the chirality sense inversion cannot be detected by a scalar chirality measure like the Hausdorff one.

IV. CONCLUSIONS

In this paper we have considered a correlation between the scaled chirality index multiplied by the transverse molecular dipole and the spontaneous polarization/tilt ratio for five series of smectic C* mesogens. In general a good correlation has been found within each series, and in the last series, with variable chiral chain, the change of the scaled chirality index has even reproduced the experimentally observed sense inversion of the spontaneous polarization. This indicates that the scaled chirality index proposed in refs. [5,6], is a useful measure of the spontaneous polarization normalised by the tilt angle and the molecular transverse dipole. In Figs. 1 and 3 the correlation between molecular chirality index and the normalized polarization is illustrated by the distribution of the corresponding points along the regression lines. At the same time, regression lines for different series of smectics C* possess different slopes. This difference can be explained using the general eq. ( 4) for the reduced spontaneous polarization. According to eq.( 4), the polarization, normalized by the tilt and the transverse dipole is not only a function of the molecular chirality, but it also proportional to the orientational order parameter of the short axes Γ and the factor of cos β where β is the angle between the primary molecular short axis and the transverse dipole.

Both parameters are expected to depend strongly on the molecular structure which may explain the difference in slopes for different series. In fact, these parameters may vary to a certain extent also within a particular series thus reducing the correlation between molecular chirality and the spontaneous polarization. Thus the scaled chirality index may be used to estimate the spontaneous polarization in a series of smectic C* mesogens with small systematic changes of the molecular structure provided the spontaneous polarization has been measured experimentally for at least one member of the series. At the same time one cannot expect any significant correlation between the polarization and the chirality index in a group of smectic C* mesogens of sufficiently different structure. This explains why only a weak correlation has been observed in ref. [16].

Perhaps the most drastic simplification in the present approach is related to the use of a single low energy molecular conformation in all calculations. It is well known that the conformation distribution for some molecules may be very degenerate, and indeed we have seen in some cases that the final configuration my depend on the initial conformation used as a starting point. In more realistic approach one should take into account a distribution of several low energy conformation in thermal equilibrium [12].

  the measure never vanishes for a chiral set of points. At the same time it is

FIG. 1 :

 1 FIG. 1: Direction and magnitude of the spontaneous polarization in the Smectic C phase. The polarization is very small if the equilibrium orientation of the molecular plane is perpendicular to the tilt plane and the transverse molecular dipole is normal to the molecular plane

  and the transverse dipole can be expressed asµ ⊥ = µ c c + µ b b,where µ c = µ ⊥ cos β, µ b = µ ⊥ sin β are the projections of the transverse dipole on the unit vectors c and b, respectively. One notes that the short axis b

10 FIG. 2 :

 102 FIG. 2: Minimum energy conformations for mesogens from the series B'

  axis (which corresponds to the largest eigenvalue) has been identified.

Page 11 of 18 URL

 18

FIG. 4 :

 4 FIG. 4: No correlation between the product of the chirality index and the full molecular dipole

  So far we have only considered the series of smectic C* mesogens where the systematic change of the molecular structure was restricted to the ring system while the chiral chain itself was not changed. In all these cases the change of the molecular chirality is caused by the change of the molecular conformation induced by variation of the substitution group in the nonchiral part of the molecule. Now we will show the results for three compounds which have the same biphenyl core and the same nonchiral n-heptyloxy-chain , whereas the chiral

FIG. 5 : 2 FIG. 6 :

 526 FIG.5: Minimum enenrgy molecular conformation, dipole moment, saturated spontaneous polarization and the chirality index for mesogens taken from[START_REF] Kelly | Handbook of Liquid Crystals[END_REF].Page 15 of 18

TABLE I :

 I Structural formulae, notation, spontaneous polarization P s , tilt angle θ and the polarization/tilt ratio P 0 for mesogens of the series A and B taken from[START_REF] Furukawa | [END_REF]. P s and θ are measured at [T = T c -10C].

					Molecular Physics			Page 8 of 18
												8
		Structural formulae				No-	P s	θ	P 0
									tation nC/cm 2 [degree] nC/cm 2
					H
			O
					O
						* *
	C8H17	o r F o r	O	O	F	O	C6H13 *	A-H	-43	42	-64
	C 8 H 17			O				C 6 H 13
		P P	CN				A-F	-101	34	-180
	C8H17	e e r O O e H e r	O	*	C6H13	A-CN	-127	22	-345
				O
	C8H17	O		O	R e O * C6H13 R F e	B-H	-45	45	-64
	C8H17	O		O O		O	v i e * C6H13 B-F v i e	-78	36	-134
	C8H17			O	Cl	O			w w
							*
	C8H17	O O		O O O	Br	O	*	C6H13 C6H13	B-Cl	O n l -123 40 O n l	-191
									B-Br	-131	y 41 y	-202
					CN
				O
	C8H17
						O
							*
	O			O				C 6 H 13
									B-CN	-160	25	-386

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

TABLE II :

 II Structural formulae, notation, spontaneous polarization P s , tilt angle θ and the polarization/tilt ratio P 0 for mesogens of the series C and D taken from[START_REF] Furukawa | [END_REF]. P s and θ are measured at [T = T c -10C].

		Structural formulae	No-	P s	θ	P 0
			tation nC/cm 2 [degree] nC/cm 2
		H
		O
		O
		*
	C 8 H 17	O
		C 6 H 13	C-H	-52	26	-120
		F
		O
		O
		*
	C 8 H 17	O
		C 6 H 13	C-F	-91	25	-211
		Cl
		O
		O
		*
	C 8 H 17	O
		C 6 H 13	C-Cl	-151	34	-270
		H
		O
	C 8 H 17	O
		*
	O	O
		C 6 H 13	D-H	-76	37	-123
		F
		O
	C 8 H 17	O
		*
	O	O
		C 6 H 13	D-F	-85	30	-170
		Cl
		O
	C 8 H 17	O
		*
	O	O
		C 6 H 13	D-Cl	-161	38	-263
		Br
		O
	C 8 H 17	O
		*
	O	O
		C 6 H 13	D-Br	-166	39	-262
		CN
		O
		O
		*
	O	O
		C 6 H 13

TABLE III :

 III Full and transverse dipole moments of the moecules from sewries A-D and the values of their chirality index G OS .

	Notation µ	µ ⊥ G OS • 10 4
	A-H 2,6610 0,4291	196
	A-F 2,3730 1,1637	-50
	A-CN 3,2660 2,8302	-70
	B-H 2,5610 1,1018	-116
	B-F 2,9570 2,4816	-121
	B-Cl 2,8030 2,1270	-126
	B-Br 2,8890 2,2671	-125
	B-CN 4,2570 4,1278	-123
	C-H 2,3480 0,4684	167
	C-F 2,0280 1,7509	135
	C-Cl 2,0690 1,7223	96
	D-H 2,4240 1,3522	80
	D-F 1,7220 1,6199	24
	D-Cl 1,6360 1,4467	12
	D-Br 1,7940 1,6284	29
	D-CN 3,0790 3,0501	18
	A1	7,4380 3,3508	65
	A2	6,3540 2,6557	219
	A4	4,3740 3,5597	13
	A6	5,1440 3,1123	23

Acknowledgments

This work has been undertaken during the stay of M.A.Osipov at the University of Stuttgart as a Mercator Guest Professor supported by DFG.


Stefan Loibl

Frederick R Manby

Martin Schütz
email: martin.schuetz@chemie.uni-regensburg.de

Density fitted, local Hartree-Fock treatment of NMR chemical shifts using London atomic orbitals

Keywords: GIAOs, NMR chemical shifts, density fitting, Hartree-Fock

An efficient program for calculating nuclear magnetic shielding tensors at the level of density fitted Hartree-Fock theory is presented. London atomic orbitals are used as AO basis functions, while ordinary Gaussians are employed as auxiliary fitting functions for the density fitting. The errors due to density fitting turn out to be very small. Accuracy and efficiency of the program are demonstrated in calculations on several test molecules.

I. INTRODUCTION

During the last decade the calculation of NMR chemical shifts has become routinely applicable at the level of ab initio electronic structure theory (for two more recent reviews see Refs. 1 and 2). One major obstacle to be overcome in such calculations is the gauge-origin problem arising from the incompleteness of the AO basis set. The Individual Gauge for Localized Orbitals (IGLOs) approach of Schindler and Kutzelnigg [START_REF] Schindler | [END_REF]4] is one way to deal with that problem, however, nowadays, the preferred solution is the use of Gauge-Including Atomic Orbitals (GIAOs, or London atomic orbitals) [5,6]. A discussion about the pros and cons of IGLOs vs. GIAOs can be found in Refs. 7, 8, and, in the context of a DFT implementation, in Ref. 9. GIAOs depend explicitly on the external magnetic field and the gauge origin. In the expressions for the atomic integrals the latter dependence cancels, rendering the desired molecular properties independent of the gauge origin [10]. Efficient programs to calculate nuclear magnetic shielding tensors on the basis of GIAOs are presently available for Hartree-Fock (HF) and Density Functional Theory (DFT) [7,[11][12][13] and linear scaling methods have also been reported [14,15]. Yet, somewhat surprisingly, to our knowledge hardly any GIAO-HF or DFT implementation for the computation of NMR chemical shifts based on the Density Fitting (DF) [16][17][18][19][20][21][22] approximation for the Electron Repulsion Integrals (ERIs) has been reported so far. An exception may be the GIAO-DFT program by Schreckenbach and Ziegler (which is part of the Amsterdam density functional package ADF), who mention in one sentence the use of auxiliary Slater type orbitals to fit the molecular density [12].

For a broad range of applications NMR chemical shifts computed at the HF or DFT level are sufficiently accurate. However, there are also many cases where a proper treatment of electron correlation, i.e., employing wavefunction based methods rather than DFT, is necessary (e.g. see sections 3.7, 3.8 in Ref. 2 and references therein). Correlated methods for calculating nuclear magnetic shielding constants with inclusion of dynamical electron correlation effects have been mainly reported for the GIAO approach. Gauss et al. presented Coupled Cluster (CC) GIAO implementations with CC models including up to triples and quadruples substitutions (see Refs. 23, 24 and references therein). Yet the application range of these methods is restricted to tiny molecules. Only the simplest electron correlation method, i.e., second-order Møller-Plesset perturbation theory (MP2) appears to be applicable to NMR chemical shift calculations of larger molecules. By virtue of integral-direct techniques [25] and full exploitation of non-Abelian point group symmetry [26] GIAO-MP2 chemical shift calculations of molecular systems comprising more than 600 basis functions could be performed [26]. However, due to the unfavorably high scaling of the computational cost with molecular size such calculations still are limited to molecules of relatively modest size or high symmetry.

From that perspective the development of correlated GIAO methods for chemical shift calculations based on local correlation approaches appears to be attractive. Efficient local correlation implementations ranging from MP2 to CC are available for the calculation of ground state energies [27][28][29][30][31][32][33] and gradients [34], as well as for properties of excited states [35][36][37]. A pilot implementation of GIAO local MP2 for chemical shifts has been presented by Gauss and Werner [38]. With this program, which simulates GIAO-LMP2 on top of a canonical GIAO-MP2 program, the authors could asses the accuracy of GIAO-LMP2. They report that the effect of the local approximation on the resulting shifts are small (smaller than other errors inherent to MP2 itself), e.g., for 13 C smaller than 1 ppm, and conclude that "GIAO-LMP2 appears to be a powerful tool for chemical applications as soon as a computationally efficient code is available".

In this contribution we report a first step towards this goal. We present an efficient GIAO program for calculating nuclear magnetic shielding tensors based on density fitted Hartree-Fock (DF-HF). DF-HF is much more efficient than conventional HF when localized occupied orbitals are used in the construction of the exchange matrix. Furthermore, linear scaling can be achieved when restricting the fitting functions to local fit domains [39]. Also the efficient DF-LMP2 analytic energy gradient (with respect to nuclear displacements) of the MOLPRO program package is based on the DF-HF energy functional (see section IIIB in Ref. 34). The GIAO-DF-HF method for chemical shifts employs Localized Molecular Orbitals (LMOs) and Projected Atomic Orbitals (PAOs) to span occupied and virtual space, respectively, in order to have sparse three-index ERIs and ERI derivatives in the Fock matrix constructions of the Coupled Perturbed Hartree-Fock step. LMOs are mutually orthogonal, while the PAOs are non-orthogonal and form an over-complete set. Ordinary Gaussians are employed as fitting functions for the orbital product densities that appear in the ERIs (since the use of GIAOs as fitting functions would trivially violate gauge origin independence of the ERIs). The error introduced by the fitting turns out to be negligible (see section III). The GIAO-DF-HF method for nuclear shieldings presented here is highly efficient, and, as for DF-HF, linear scaling could be achieved by introducing local fit domains. This work was originally considered as a necessary step towards an efficient GIAO-DF-LMP2 code. However, the GIAO-DF-HF program is interesting by itself because it makes it possible to calculate HF NMR chemical shifts efficiently within the MOLPRO program package.

II. THEORY A. The NMR chemical shielding tensor

The NMR chemical shielding tensor is defined as the mixed second derivative of the electronic energy with respect to the components of the external magnetic field B and the magnetic moment m Z of nucleus Z:

σ Z βα ≡ d 2 E HF dB α dm Zβ B,mZ=0
.

For calculations in the context of the GIAO ansatz explicitly field-dependent basis functions ω µ (r M , A M ) are used, i.e.,

ω µ (r, A M ) ≡ exp (- i c A M • r)χ µ (r M ), (2)
where A M is the vector potential representing the magnetic field at nucleus M

A M = 1 2 B × (R M -R O ) = 1 2 B × R MO , (3)
and χ µ (r M ) the field-independent basis-functions. R M and R O represent the position vector of nucleus M and the gauge origin, R MO is their difference vector. Respectively, r denotes the position vector of an electron, and r M the vector pointing from nucleus M to this electron.

Greek letters µ, ν, . . . label atomic orbitals; i, j, k, . . . occupied localized molecular orbitals (LMOs); a, b, c, . . . projected atomic orbitals (PAOs); and p general localized MOs. Orbital indices decorated with a bar, e.g. ī, j, k, . . . and ā, b, c, . . . denote canonical occupied and virtual orbitals, respectively.

The shielding tensor can be written more explicitly as

σ Z βα = µν D 0 µν ∂ 2 h µν ∂B α ∂m Zβ + µν D Bα µν ∂h µν ∂m Zβ B=0 , (4)
D 0 µν ≡ 2 i L µi L νi , (5)
D Bα µν ≡ 2 i L Bα µi L νi + L µi L Bα νi * . (6)
The derivatives of the LMO coefficients in eq. ( 6) can be expressed via the transformation matrix U Bα , i.e.,

L Bα µi = k L µk U Bα ki + c P µc U Bα ci . (7)
Here, the matrices L and P represent LMO coefficient matrix and the projection of the original AOs onto the virtual space, respectively (see e.g. eqs. [START_REF] Gauss | Molecular Properties[END_REF][START_REF] Schindler | [END_REF] in Ref. 40), and

P = 1 -LL † S AO . (8)
The one-electron derivative integrals appearing in eq. ( 4) can be written as

∂ 2 h µν ∂B α ∂m Zβ = 1 2c 2 µ (R MN × r) α (r Z × ∇) β r 3 Z ν + µ (r N • r Z ) δ αβ -(r N ) α (r Z ) β r 3 Z ν , (9)
∂h µν ∂m Zβ = -i c µ (r Z × ∇) β r 3 Z ν . (10)
The explicit dependence of the integrals on the gauge origin R O cancels: both through the vector potential in the kinetic energy operator π = -i∇ + A(r), as well as through the GIAOs cancels (for details, see e.g. Ref. 10). The derivative integrals in eqs. (9,10) hence are independent of the choice of R O . Furthermore, the derivative of the two-electron part does not occur in eq. ( 4) because it is independent of the perturbation m Z . Since σ Z βα is evaluated at zero field strength the GIAOs of the three one-electron integrals of eqs. (9,10) reduce to ordinary Gaussian basis functions.

In order to construct the perturbed density matrix according to eq. ( 6) the general-occupied block of the transformation matrix U Bα has to be calculated. The requirement for preserving orthonormality of the orbitals spanning the occupied space, as well as strong orthogonality to the virtual space yields the condition

U Bα pi + U Bα ip * = -S Bα pi , (11)
with

S Bα µν = i 2c µ |(R MN × r) α | ν . (12)
Again, the R O dependence of the integral S Bα pq (here through the GIAOs only) trivially cancels. Eq. ( 11) determines the occupied-occupied block of U Bα to be

U Bα ij = - 1 2 S Bα ij . (13)
The virtual-occupied block of U Bα is determined by the requirement for preserving the Brillouin theorem and is obtained by solving the Coupled Perturbed Hartree-Fock (CPHF) equations

0 ! = R ai = c f ac U Bα ci - k f ik U Bα ak - k f ik S Bα ak + G(D Bα ) ai + ∂h ai ∂B α + ∂G(D 0 ) ai ∂B α , (14)
∂h µν ∂B α = i 2c µ |(R MN × r) α h| ν -µ |(r N × ∇) α | ν , (15)
∂ µρ|νσ ∂B α = i 2c µρ 1 r 12 (R MN × r 1 ) α νσ + µρ 1 r 12 (R RS × r 2 ) α νσ . ( 16
)
Once more, the R O dependence cancels in these integrals, for the electron repulsion integral again trivially, since there R O occurs in all of the GIAOs and not in the operator.

To ensure fast convergence for the iterative solution of the CPHF equations the solution U Bα is not updated in the local LMO/PAO basis but in the canonical basis. The residual R is therefore transformed from the local to the canonical basis,

R āī = aci Q āc V † ca R ai W † i ī, (17)
which is then used to generate the update for the new solution in each iteration step. The matrix Q is defined as

Q = C v † S AO , (18)
with C v and S AO representing the virtual block of the canonical MO coefficient matrix and the AO overlap matrix. W is the usual unitary transformation matrix specified by the localization criterion transforming canonical occupied orbitals φ CAN i to LMOs φ LMO i . The matrix V is a pseudo-inverse of S PAO (see Ref. 41). R āī and V are sparse, however the second transformation step in eq. ( 17) involving W scales cubically, yet with a very low pre-factor. Therefore, the transformation of R to canonical basis does not present itself as a bottleneck, so far.

B. Density fitting approximation for electron repulsion integrals

For GIAO-DF-HF all ERIs are approximated by DF. Hence, for the calculation of the NMR chemical shielding tensor the expressions for the G(D Bα ) and ∂G(D 0 ) µν /∂B α in the CPHF equations (14) have to be modified accordingly. The former is originally defined as

G(D Bα ) ai = ρσ D Bα ρσ (ai|ρσ) - 1 2 (aσ|ρi) = - 1 2 ρσ D Bα ρσ (aσ|ρi), (19)
where (ai|ρσ), (aσ|ρi) represent the four-index ERIs (half-transformed to LMO/PAO basis), which are to be approximated by DF. For the last equality sign in the previous equation we have used the fact that D Bα is a purely imaginary hermitian matrix and as a consequence only the exchange part survives. Inserting eq. ( 6) yields

G(D Bα ) ai = - k (ak| ki) -(a k|ki) , (20)
where k denotes a perturbed LMO χk (r), transformed with the coefficient matrix L Bα defined in eq. ( 7).

The DF-approximation implies the expansion of the orbital product 'densities' χ k (r)χ a (r), χk (r)χ i (r), etc. appearing in the ERIs in an auxiliary basis of fitting functions Ξ P (r). This leads to a factorization of the four-index ERIs in terms of three-index objects. From eq. ( 20) one so obtains

G(D Bα ) ai = - k P c P ak (P | ki) -(a k|P )c P ki , (21)
(a k|P ) = µν P µa L Bα νk (µν|P ), (23) and the fitting coefficients

c P pi = Q J -1 P Q (Q|pi). (24)
In eqs. (21)(22)(23)(24) the indices P, Q stand for auxiliary Fitting Functions (FFs), and J P Q = (P |Q) is the Coulomb metric of these functions. Observe that the fitting coefficients c P pi are independent of the CPHF iterations and thus just computed once. Furthermore, the three-index ERIs of eqs. (22)(23) split into two parts according to eq. ( 7), with the first term being independent of the CPHF iterations, as well. The quantity G(D Bα ) ai of eqs. (20,21) is, as any other term appearing in the CPHF equations ( 14), purely imaginary. The CPHF equations thus can be considered as real valued and solved for the imaginary part of the (purely imaginary) U Bα .

An important note concerns the FFs. Using GIAOs as FFs would inevitably violate gauge invariance. Since for a given FF there is naturally no complex conjugate corresponding to the same electron the origin dependence (on R O ) would not cancel in the three-index ERIs. An alternative natural choice are ordinary Gaussian basis functions, implying that the GIAO orbital product densities are fitted at zero field, i.e., at B = 0. In the present implementation no local restrictions to the fitting basis (fit domains [31,33]) were yet introduced. However, since the orbital product 'densities' χ k (r)χ a (r), χk (r)χ i (r), etc., are intrinsically local (provided that the solution of the CPHF equations U Bα entering eq. ( 7) is local), this could be accomplished. Substantial computational savings and eventually linear scaling of the computational cost can be expected for such a local fitting scheme. Since the individual local fitting basis sets then are different for the two orbital products in each ERI of eq. ( 20) the robust three-term formula has to be applied in the context of local fitting [31,33]), leading to four, rather than two contractions in eq. ( 21). The possibilities of local fitting in the context of GIAO-DF-HF will be explored in forthcoming work.

The ∂G(D 0 ) ai /∂B α part of the CPHF equations is also independent of the iteration and can be computed just once. Originally, it involves the contraction of the four-index derivative integrals defined in eq. ( 16) with the unperturbed density defined in eq. ( 5), i.e.,

∂G(D

0 ) ai ∂B α = 2 k ∂(ai|kk) ∂B α - 1 2
∂(ak|ki) ∂B α . ( 25
)
Invoking the DF-approximation yields

∂G(D 0 ) ai ∂B α = P (ai|P ) Bα c P - P k (ak|P ) Bα c P ki + c P ak (ki|P ) Bα , (26)
with

c P = Qρσ D 0 ρσ J -1 P Q (ρσ|Q), (27) and (µν|P )
Bα = i 2c µν 1 r 12 (R MN × r 1 ) α P . (28)
In the three-index derivative ERIs defined in eq. ( 28) the R O dependence of the GIAOs of the orbital product density cancels, while the FF does not carry any R O dependence by construction. Furthermore, due to the fact that these ERIs are needed at B = 0 the GIAOs reduce to ordinary Gaussian basis functions. The new GIAO-DF-HF method for calculating nuclear magnetic shielding tensors is implemented in the MOLPRO package [START_REF] Werner | MOLPRO, version 2009.2, a package of ab initio programs Cardiff[END_REF]. Most of the time-critical subroutines are parallelized, yet the scratch files containing the transformed three-index integrals and fitting coefficients reside on a common shared file system. In this section we present results from some test calculations, demonstrating accuracy and efficiency of the new program.

Tables I, II, and III compile the 13 C, 17 O, and 15 N chemical NMR shifts, calculated for a subset of the molecules used by Gauss in Ref. [START_REF] Gauss | [END_REF]. The geometries of these molecules were optimized at the level of MP2 with the cc-pVTZ AO basis. GIAO-DF-HF calculations of the nuclear magnetic shielding tensors then were performed at these geometries, employing the cc-pVXZ AO basis, for X=D, T, and Q, respectively. As fitting basis sets the JK-fitting basis sets of Weigend [START_REF] Weigend | A fully direct RI-HF algorithm: optimized auxiliary basis sets, demonstration of accuracy and efficiencyPhys[END_REF] related to cc-pVXZ, cc-pV(X+1)Z, and cc-pV(X+2)Z (if available) were used. Reference calculations were performed by using the conventional GIAO-HF implementation in TURBOMOLE [11]. From these results it is evident that the use of DF has virtually no effect on the resulting chemical shifts. Even without enlarging the cardinal number of the fitting basis relative to that of the AO basis the deviation relative to the reference value is in the far sub-ppm range, i.e., about 0.1 ppm or less. The simple approach of employing ordinary Gaussians as fitting functions, or, in other words, to fit the orbital product densities in the ERIs at zero field (see discussion in section II B) hence appears to be the right way of applying DF in the context of property calculations involving GIAOs.

The efficiency of the GIAO-DF-HF program was explored by carrying out calculations on three somewhat more extended molecular systems, i.e. on coronene, on a tweezer host-guest complex 1@2 "clinching" the 1,4-dicyanobenzene guest molecule 2, and on a phenothiazine-isoalloxazine-pyrene triad (see Fig. 1). For the coronene molecule and the 1@2 tweezer experimental and theoretical results on NMR shifts are available in the literature (cf. Refs. 45, 46 and references therein). For the calculations presented in this work the geometries of these molecules were optimized at the HF level using the SVP (coronene) and the 6-31G * (1@2) basis sets. Since the structure of the 1@2 tweezer so obtained deviates from that reported in Ref. 46, additional convential GIAO-HF calculations at the present structure were carried out for reference, rather than adopting the 1 H NMR shifts provided in Ref. 46. The triad has been used as a test molecule already in previous work [35,37] and its structure has not been re-optimized for the present purpose.

Table IV contains the 1 H NMR shifts calculated with GIAO-DF-HF and a convential GIAO-HF in TURBOMOLE for coronene and 1@2 at their optimized geometries. The results of the GIAO-DF-HF calcaulation are essentially identical to the results of the GIAO-HF calculations. For the triad we just present the GIAO-DF-HF 17 O, and 15 N chemical shifts in Table V, a conventional reference calculation was not attempted for this case.

Tables IV and V also compile the CPU times (per processor) and the elapsed times, measured for the essential steps of the calculation, i.e., for the Fock build involving the ERI derivative integrals for the right-hand side of the CPHF equations according to eq. ( 26), and for an individual CPHF iteration. For comparison, the corresponding timings for an SCF iteration of the initial DF-HF calculation is also provided. The calculations were run on four cores, Intel(R) Xeon(R) CPU E5462 @ 2.80GHz, without exploiting any point group symmetry. For the calculation of the tweezer host-guest complex in the larger cc-pVTZ basis, solving the CPHF equations takes about 65 minutes elapsed to build up the right-hand side, plus about 30 minutes per iteration to converge the equations. An individual DF-HF SCF cycle, on the other hand, takes about 35 minutes. For the triad in the larger cc-pVTZ basis, a CPHF iteration costs about 140 minutes, which compares to 100 minutes for an SCF iteration. By comparing CPU and elapsed times in Table V, it can be seen that I/O becomes a real bottleneck in our present parallel implementation. Yet there is certainly room for further optimization of the code in that respect.

IV. CONCLUSIONS

In this contribution we present an efficient program for calculating nuclear magnetic shielding tensors at the level of density fitted (DF) Hartree-Fock for extended molecular systems. As basis functions Gauge-Including Atomic Orbitals (GIAOs, London atomic orbitals) are used. The orbital product densities appearing in the Electron Repulsion Integrals (ERIs) are expanded in a fitting basis consisting of ordinary Gaussians. GIAOs cannot be used here, since the dependence on the gauge origin would then not cancel in the ERIs. Using ordinary Gaussians as fitting functions corresponds to a fit of the orbital product densities at zero field strength of the external magnetic field. This simple approach works remarkably well: the dependence of the resulting chemical shifts on the size of the fitting basis set is virtually negligible and the agreement with results from conventional calculations without DF is typically within IV: GIAO-DF-HF 1 H chemical shifts δ (in ppm) calculated for the tweezer host-guest complex 1@2 and the coronene molecule, respectively, along with the corresponding conventional reference values. The shifts are given relative to TMS reference ab . In addition, the CPU times (per processor) and the elapsed times (in minutes) measured for the individual key steps of the GIAO-DF-HF calculation, i.e., the timings for the setup of the rhs of the CPHF equations, and the averaged CPHF iteration times are given. For comparison, the averaged SCF iteration times of the initial DF-HF calculation are also provided.

Basis cc-pVDZ cc-pVTZ

Tweezer host-guest complex 1@2 a GIAO-DF-HF ref. V: GIAO-DF-HF relative chemical shifts (in ppm) for the heteroatoms of the phenothiazine-isoalloxazine-pyrene triad. 15 N chemical shifts are relative to NH3 a , 17 O chemical shifts are relative to H2O b . In addition, the CPU times (per processor) and the elapsed times (in minutes) measured for the individual key steps of the GIAO-DF-HF calculation, i.e., the timings for the setup of the rhs of the CPHF equations, and the averaged CPHF iteration times are given. For comparison, the averaged SCF iteration times of the initial DF-HF calculation are also provided.

cc-pVDZ cc-pVTZ

Triad

  and (purely imaginary) perturbed density matrices

  one-and two-electron derivative integrals

  -index ERIs in LMO/PAO basis (P | ki) = µν (P |µν)L Bα µk L νi ,

TABLE I :

 I 13 C chemical shifts (δ, ppm) obtained with the new GIAO-DF-HF program, and compared to conventional GIAO-HF reference values. All shifts relative to CH4 a as in Ref.43.

	Basis		cc-pVDZ			cc-pVTZ		cc-pVQZ
	Fit basis	VDZ	VTZ	VQZ	ref. b	VTZ	VQZ	V5Z	ref. b	VQZ	V5Z	ref. b
	C6H6	134.14 134.12 134.13 134.14	138.48 138.49 138.49 138.50	140.45 140.45 140.45
	CF4	109.92 109.88 109.88 109.89	115.42 115.42 115.41 115.41	116.11 116.11 116.11
	CH3F	64.51	64.49	64.49 64.50	69.22	69.22	69.22 69.22	69.90	69.90 69.90
	C2H2	74.51	74.50	74.50 74.51	79.26	79.26	79.26 79.26	80.18	80.18 80.18
	C2H4	128.37 128.36 128.37 128.38	133.09 133.09 133.09 133.10	135.25 135.25 135.26
	CH2O	193.30 193.27 193.27 193.27	200.50 200.48 200.48 200.48	203.52 203.51 203.52
	(CH3)2CO	211.25 211.19 211.19 211.20	214.74 214.72 214.72 214.73	219.02 219.01 219.02
	(CH3)2CO	29.70	29.69	29.70 29.71	31.21	31.21	31.21 31.22	31.90	31.89 31.91
	CH3CH3 CH3CN CH3CN CH3CHO CH3CHO	10.95 125.03 124.99 125.00 125.01 10.93 10.93 10.94 4.18 4.17 4.17 4.18 30.95 30.94 30.95 30.96 201.07 201.03 201.03 201.03 F o r	11.38 131.96 131.95 131.95 131.96 11.38 11.38 11.38 4.10 4.09 4.09 4.11 32.49 32.49 32.49 32.49 206.89 206.87 206.87 206.88	11.91 134.61 134.61 134.62 11.91 11.91 4.49 4.49 4.50 33.41 33.41 33.42 210.48 210.47 210.49
	CH3NH2	29.77	29.76	29.76 29.77	31.27	31.27	31.27 31.28	31.99	31.99 32.00
	CH3OH CO2	47.46 136.35 136.30 136.30 136.31 47.44 47.44 47.45 P	50.69 145.40 145.38 145.38 145.38 50.69 50.69 50.70	51.68 146.59 146.59 146.59 51.68 51.69
	CO a) Absolute chemical shifts for CH4: 205.36 (cc-pVDZ/cc-pVDZ), 205.37 (cc-pVDZ/cc-pVTZ), 205.38 (cc-pVDZ/cc-pVQZ), 214.29 214.25 214.25 214.28 222.24 222.21 222.20 222.23 225.16 225.15 225.18 e 205.38 (ref., cc-pVDZ), 196.84 (cc-pVTZ/cc-pVTZ), 196.84 (cc-pVTZ/cc-pVQZ), 196.84 (cc-pVTZ/cc-pV5Z), 196.84 (ref., cc-pVTZ), 195.66 (cc-pVQZ/cc-pVQZ), 195.66 (cc-pVQZ/cc-pV5Z), 195.66 (ref., cc-pVQZ). The basis sets are specified as (AO basis / corresponding JKfit basis). e r
	b) GIAO HF values calculated with TURBOMOLE.	R
						e
						v i e
							w
								O n l
										y

TABLE II :

 II 17 O chemical shifts (δ, ppm) obtained with the new GIAO-DF-HF program, and compared to conventional GIAO-HF reference values. All shifts relative to H2O a as in Ref.43.

	Basis		cc-pVDZ			cc-pVTZ			cc-pVQZ
	Fit basis	VDZ	VTZ	VQZ	ref. b	VTZ	VQZ	V5Z	ref. b	VQZ	V5Z	ref. b
	CH3OH	-5.37	-5.41	-5.42	-5.42	-8.67	-8.68	-8.68	-8.67	-11.07 -11.07 -11.06
	H2O2	189.58 189.57 189.55 189.58	188.75 188.74 188.74 188.75	187.13 187.13 187.14
	CO	417.66 417.59 417.57 417.56	429.26 429.21 429.21 429.20	426.69 426.68 426.67
	CO2	120.44 120.36 120.34 120.35	120.84 120.80 120.81 120.82	116.10 116.15 116.16
	CH2O	784.19 784.04 784.00 784.00	785.14 785.16 785.13 785.13	781.00 781.14 781.10
	CH3CHO	700.69 700.58 700.58 700.57	700.96 700.93 700.93 700.92	697.56 697.50 697.44
	(CH3)2CO	669.37 669.30 669.31 669.32	669.73 669.69 669.70 669.69	667.06 666.99 666.95
	a) Absolute chemical shifts for H2O: 347.16 (cc-pVDZ/cc-pVDZ), 347.18 (cc-pVDZ/cc-pVTZ), 347.18 (cc-pVDZ/cc-pVQZ),
	347.19 (ref., cc-pVDZ), 334.90 (cc-pVTZ/cc-pVTZ), 334.89 (cc-pVTZ/cc-pVQZ), 334.90 (cc-pVTZ/cc-pV5Z), 334.91 (ref., cc-
	pVTZ), 330.74 (cc-pVQZ/cc-pVQZ), 330.74 (cc-pVQZ/cc-pV5Z), 330.75 (ref., cc-pVQZ). The basis sets are specified as (AO
	basis / corresponding JKfit basis).
	b) GIAO HF values calculated with TURBOMOLE.

TABLE III :

 III 15 N chemical shifts (δ, ppm) obtained with the new GIAO-DF-HF program, and compared to conventional GIAO-HF reference values. All shifts relative to NH3 a as in Ref.43.

			F o r
	Basis Fit basis	VDZ	cc-pVDZ P VTZ VQZ ref. b e	VTZ	cc-pVTZ VQZ V5Z	ref. b	VQZ	cc-pVQZ V5Z	ref. b
	HCN CH3NH2 CH3CN	309.85 309.80 309.77 309.77 14.29 14.27 14.26 14.26 302.10 302.06 302.05 302.00 e 322.86 322.85 322.85 322.85 14.82 14.81 14.82 14.83 313.92 313.90 313.90 313.85 r	322.57 322.56 322.56 13.47 13.47 13.47 313.38 313.38 313.34
	N2 pVTZ), 265.51 (cc-pVQZ/cc-pVQZ), 265.51 (cc-pVQZ/cc-pV5Z), 265.51 (ref., cc-pVQZ). The basis sets are specified as (AO 383.64 383.54 383.52 383.51 393.16 393.12 393.13 393.12 393.78 393.78 393.78 R cc-e basis / corresponding JKfit basis). b) GIAO HF values calculated with TURBOMOLE. v i e
					w
						O n l
							y

a) Absolute chemical shifts for NH3: 278.27 (cc-pVDZ/cc-pVDZ), 278.28 (cc-pVDZ/cc-pVTZ), 278.28 (cc-pVDZ/cc-pVQZ), 278.

28 (ref.

, cc-pVDZ), 269.18 (cc-pVTZ/cc-pVTZ), 269.17 (cc-pVTZ/cc-pVQZ), 269.17 (cc-pVTZ/cc-pV5Z), 269.

18 (ref.,

TABLE

  level, absolute chemical shielding 32.2 ppm (cc-pVDZ), 32.1 ppm (cc-pVTZ). b) TMS optimized at HF/SVP level, absolute chemical shieldings 32.0 ppm (cc-pVDZ), 31.9 (cc-pVTZ). c) GIAO HF values calculated with TURBOMOLE.

			c	GIAO-DF-HF	ref. c
	H a (guest)	5.8	5.8	5.8	5.8
	H b (guest)	3.7	3.7	3.5	3.5
	H 2,3,14,15	6.7	6.7	6.6	6.6
	Harom(host)	7.2 -7.8	7.2 -7.8	7.3 -7.8	7.3 -7.8
	H bridgehead H 25,28 H 26,27 Building eq. (26) CPU (elapsed) time / min CPHF iteration step CPU (elapsed) time / min DF-HF iteration step CPU (elapsed) time / min Coronene b H Building eq. (26) CPU (elapsed) time / min CPHF iteration step CPU (elapsed) time / min DF-HF iteration step CPU (elapsed) time / min o r P e e r R GIAO-DF-HF 3.7 -4.0 2.2 1.9 -2.1 13.4 (13.4) 3.7 -4.0 2.2 1.9 -2.1 10.8 (11.0) 9.4 (9.5) ref. c 9.0 9.0 0.9 (0.9) 0.6 (0.6) 0.4 (0.4) e a) TMS optimized at HF/6-31G F v i e	3.9 -4.1 2.2 2.0 -2.1 65.0 (65.2) 3.9 -4.1 2.2 2.0 -2.1 28.7 (29.0) 33.2 (33.4) GIAO-DF-HF ref. c 9.2 9.2 3.3 (3.3) 2.1 (2.1) 1.1 (1.1)
		w
			O n l
				y

*

TABLE

  NH3 optimized at MP2/cc-pVTZ level, absolute chemical shifts 278.3 (cc-pVDZ), 269.2 (cc-pVTZ). b) H2O optimized at MP2/cc-pVTZ level, absolute chemical shifts 347.2 (cc-pVDZ), 334.9 (cc-pVTZ).

	N 1 N 2 N 3 N 4 N 5 N 6 O 1 O 2 O 3 O 4 O 5 Building eq. (26) CPU (elapsed) time / min 30.1 (30.7) 129.9 166.9 217.1 457.4 164.3 91.4 180.8 400.3 404.3 364.2 444.4 CPHF iteration step CPU (elapsed) time / min 26.6 (27.2) DF-HF iteration step CPU (elapsed) time / min 28.1 (28.4) a) Page 12 of 12 134.6 171.8 224.4 474.7 170.5 95.6 187.5 406.3 406.8 366.9 446.0 99.3 (130.2) 71.0 (141.6) 97.6 (98.9)

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

Acknowledgments

S.L. gratefully acknowledges a scholarhip from the German National Academic Foundation.

0.1 ppm or less. The efficiency of the new program is competitive, though some aspects could be further improved. This new program now provides the platform to develop a correlated GIAO method for nuclear magnetic shielding tensors based on DF Local MP2 [34], which is the primary target of this development.


Timothy G Wright
email: btim.wright@nottingham.ac.uk

Larry Alan Viehland
email: l.viehland@chatham.edu

Benjamin R Gray

Interactions of Rare Gas Cations with Lighter Rare Gas Atoms

Keywords: rare gas, complex, spectroscopy, cation, mobility

High-level ab initio calculations are employed to generate potential energy curves for rare gas cations, RG + , interacting with neutral atoms of other rare gases, RG′, that are lighter (RG + /RG′ = Ne-Rn). The calculations employ the RCCSD(T) method, with doubly-augmented basis sets of quintuple-ζ quality. The interaction potential curves, with the full counterpoise correction applied, are calculated point-by-point. Spin-orbit coupling is applied analytically in an atom-based model.

The potentials are used to calculate spectroscopic parameters, which are then compared to recent experimental work, and the very limited previous theoretical work. In addition, the potentials are used to calculate ion transport properties, and the ion mobilities are compared to the few experimental data available.

I. Introduction

In recent papers 1,2 , we have reported high-level, ab initio potential energy curves for the interactions of RG + ions with He. Those calculations employed the restricted coupled cluster method with singles, doubles and non-iterative triples, RCCSD(T), with large basis sets that, when necessary, are extrapolated to the complete basis set limit. The full counterpoise correction was performed at each separation and spin-orbit coupling was taken into account. The calculated spectroscopic parameters were generally in very good agreement with those obtained from UV absorption and microwave spectra. In addition, the calculated potentials were in good agreement with empirically-fitted potentials obtained from microwave spectra of energy levels close to the dissociation limit. Finally, ion mobilities in helium calculated from the interaction potential energy curves were also in good agreement with experimental results. 1,3,4,5,6 The purpose of the present paper is to extend this work to other cases in which a rare gas cation, RG + , interacts with a lighter rare gas atom (RG′).

The first emission spectrum from a RG/RG′ mixture was obtained by Druyvesteyn 7 in Ne/He mixtures. These studies were confirmed and extended to Kr/Ar, Xe/Ar and Xe/Kr mixtures by Oskam and Jongerius in 1958/9. 8,9 Excitation by electron beams produced similar emissions. 10,11,12 Tanaka et al. 13 were the first to suggest that these emissions might be due to electronic transitions within RG + -RG′ ions, and they produced emission spectra of the ten possible pairs of species for RG + , RG′ = He-Xe by electric discharge.

We have outlined the previous work on the RG + -He species in our recent work. 1,2 For the heavier species under consideration here, we note that there are emission studies on Kr + -Ar by Holland et al. 14 and Tsuji et al. 15 , on Xe + -Ar by Huber and Lipson 16 and on Xe + -Kr by Tsuji et al. 17 Photoionization and photoelectron spectroscopic studies have been carried out on RG-Ne (RG = Ar-Xe); 18 Xe-RG′ (RG′ = Ne-Kr); 19 Kr-Ar, Xe-Ar and Xe-Kr; 20,21 Xe-Ar; 22 and Xe-Kr, 23 allowing information on the corresponding cations to be obtained. The combined results from these studies are usefully summarized in Table I of ref. 21. Threshold photoelectron-photoion coincidence (TPEPICO) spectra have also yielded information on the lower cationic levels for Kr-Ar 24 , Xe-Kr 25 , and Xe-Ar. 26 By far the most detailed information, however, has come from recent zero-kinetic energy (ZEKE) spectroscopic studies by Merkt's group, who have tackled Xe + -Kr, 27 Kr + -Ar, 28 and Xe + -Ar. 29 Reference and comparison to the results of these experimental results will be made below. Ion mobilities have been measured for a number of these systems. The mobilities of Ar + in Ne at 77 K have been reported by Helm and Elford, 30 Kr + in Ar at 300 K by Kaneko et al. 31 and Xe + in Ne at 300 K by Johnsen and Biondi, 32 and in Ar at 294 K by Jones et al. 33 For the purposes of the present work, we have used the smoothed versions of the experimental mobility data, as reported in refs. 4 and 5.

Interestingly, in contrast to the RG + -He species, 2 and particularly Ne + -He, 1,2 there have been very few theoretical studies on the heavier RG + -RG′ species. There are model potential calculations on all ten heterodiatomic species by Hausamann and Morgner 34 ; empirically-fitted results for Ar + -Ne by Brunetti et al., 35 and also HF and MP2 optimized geometries for Ar + -Ne, Kr + -Ne and Kr + -Ar, together with D 0 values calculated using G2 energetics. 36 There appear to be no previous reports of ion transport properties calculated from potential energy curves for the titular systems.

II. Calculational Details and Results

A. Potential Energy Curves

We have performed calculations with the MOLPRO suite of programs 37 at the RCCSD(T) level of theory. It has been argued 38 that single-determinant methods such as this might not be sufficient to be able to describe contributions from charge transfer in RG + -RG′ systems. However, our previous calculations 1 with Ne + -He found little multireference character in the wave functions for the lowlying states of interest here. Moreover, the RCCSD(T) method has the advantages of recovering a large fraction of the dynamic correlation energy and of being size consistent. Hence we have again chosen to use this method.

Our calculations employed aug-cc-pV5Z basis sets 39 , augmented with additional diffuse functions of all angular moment types in the original set. These additional functions were obtained in an even-tempered fashion from the two most diffuse basis functions in the original set, and the augmented set is referred to as d-aug-cc-pV5Z. For the heavier atoms, small-core effective core potentials (ECP10MDF, ECP28MDF and ECP60MDF for Kr, Xe and Rn, respectively) and were supplemented with the d-aug-cc-pV5Z-PP basis sets for the valence electrons. 40 The non-spinorbit (non-SO) interaction potentials of the 2 Σ + and 2 Π states were constructed by making the full counterpoise correction at each point. Since the interaction between RG + and RG′ is weak, the spin-orbit coupling is localized on RG + . Therefore the spin-orbit (SO) potentials for the X 2 Σ 1/2 + , A 1 2 Π 3/2 and A 2 2 Π 1/2 states were calculated directly from the non-SO potentials for the 2 Σ + and 2 Π states, using equations for the so-called atomic model that are well known. 2,41 The values of ∆, the SO splitting between the RG + atomic 2 P 3/2 and 2 P 1/2 states from ref. 42 were employed; these are 1431.41, 5371.0, 10537.01 and 30895.1 cm -1 for RG = Ar-Rn, respectively. The final potentials for both the non-SO and SO states are given as numerical tables as Electronic Supplementary Information.

B. Spectroscopic Values

Our potential energy curves were used as input to the LEVEL program. [START_REF] Leroy | Level 7.2 -A computer program for solving the radial Schrödinger equation for bound and quasibound levels, and calculating various values and matrix elements[END_REF] Least-squares fitting of calculated rovibrational energies to standard Morse expressions was employed to determine spectroscopic constants from the lowest few levels. The results are given in Tables 1 and2 for the SO and non-SO potentials, respectively.

C. Gaseous Ion Mobilities

We calculated the transport cross sections for RG + in RG′ from the SO ab initio interaction potential energy curves as functions of the ion-neutral collision energy using the classical-mechanical program QVALUES. [START_REF] Viehland | [END_REF]45 The transport coefficients depend upon weighted averages of these values. Since we do not know the ratios that were present in the previous experiments, or those that may be present in future work, we have calculated the mobilities based upon the following cross sections: a 1:1 average of the 2 Σ 1/2 + and 2 Π 3/2 molecular states that arise from RG + ( 2 Π 3/2 ) +RG′ interactions, with SO included (denoted "Ground" hereafter); pure 2 Π 1/2 , arising from RG + ( 2 Π 1/2 ) +RG (denoted "Excited" hereafter); and, a 1:1:1 average of the three SO states (denoted "Statistical" hereafter). It is expected that the Ground results are most likely to be applicable to the previous experiments, given the large spin-orbit coupling. If the method for generating RG + is such that the excited SO component is populated, then the experimental results will lie between the Excited and Ground values, depending on the method of production, as well as spin-orbit relaxation processes present in the experiment. The Statistical calculations are expected to be most applicable at large ratios, E/n 0 , of the electrostatic field strength to the gas number density, where the magnitude of the field is sufficient to scramble the atomic ion states.

The weighted cross sections as a function of collision energy were used in the program GC 46,47 to determine the gaseous ion mobility and the other gaseous ion transport coefficients as functions of E/n 0 at gas temperatures of 100, 200, 300, 400 and 500 K, as well as at any gas temperatures for which experimental values are known. The calculated mobilities are generally precise within 0.1%, which means that the numerical procedures within programs QVALUES and GC have converged within 0.1% for the ion-neutral interaction potential. Calculations were performed for every isotope of the ion that is naturally occurring, but in all cases the RG′ gas was assumed to have the natural abundance of isotopes. The numerous results have been placed in the database maintained at Chatham University. 48 Figure 2 compares the calculated mobilities with the limited data available for these systems. The differences can also be more easily understood in terms of the statistical quantities, δ and χ, given in Table 3. The dimensionless quantity δ is a measure of the relative difference compared to the combined experimental and calculational errors, while χ is a measure of the relative standard deviation compared to the sum of the squared error estimates. 45

Discussion

A. Spectroscopic Quantities

A.1 Non-SO Potentials

The spectroscopic quantities obtained from the non-SO potentials are presented in Table 1. Since it is only possible to compare the experimental results to the SO potentials, the values in Table 1 can only be compared to previous theoretical results where the spin-orbit interaction was not considered. The only such case appears to be that of Ma et al., 36 where good agreement is seen in the case of Ne-Ar + for both R e (at the MP2 level) and D 0 (calculated employing the G2 approach).

Good agreement is also seen for R e in the case of Ne-Kr + , but surprisingly this species was calculated to be unbound in that work. Similarly, for Ar-Kr + , the R e value is in agreement with the present work, whereas the D 0 value is about half of the present calculated value; this suggests that these species demand the much higher levels of theory and basis sets used here. The spectroscopic quantities obtained from the SO potentials are presented in Table 2. Note that with the states considered in the present work, the 2 Π 3/2 state parameters are identical to those for the non-SO 2 Π potential. We show the curves for Ne-Ar + in Figure 1 as an example.

For the X 2 Σ 1/2 + state of Ne-Ar + there is very good agreement between the present values, the calculated model potential of Hausamann and Morgner 34 and the CI values of Brunetti et al. 35 agreement with the present one. The value of ≤4585 cm -1 obtained from emission spectroscopy 16 is also consistent with the value we obtain; better agreement is seen with the PEPICO value of 4313±12 cm -1 and the recent ZEKE value 28 of 4326.9±5 cm -1 , although these values are still somewhat higher than ours. On the other hand, very good agreement is seen for ω e , where our value is in excellent agreement with the PEPICO 24 and ZEKE 28 studies, with the agreement for ω e x e also being good; the ω e value from the emission study 16 is somewhat too high. Similar comments apply to the 2 Π 3/2 state, where again the agreement with the ZEKE 28 work is extremely good for all values, but the D e value from the model potential 34 is significantly higher, as is the PEPICO 24 value. This is in stark contrast with the comparison for the 2 Π 1/2 state where the ZEKE value 28 for D e is significantly higher than the value obtained in the present work (with the model potential 34 again also giving a value higher than both of these); although the PEPICO 24 value is also higher, it is much closer to the present value. The ZEKE 28 ω e value is also significantly larger than the present value; in contrast to the excellent agreement seen for the other two states; on the other hand, the PEPICO 24 26 and ZEKE studies 29 , is extremely good, and similar comments apply to the 2 Π 3/2 state. As with the case of Ar-Kr + , however, the experimental D e and ω e values for the 2 Π 1/2 state are somewhat higher than the present ones in contrast to the agreement seen for the other states. No previous results have been reported for Ar + -Rn.

For Kr-Xe + , in contrast to the other species, the binding energy for the X 2 Σ 1/2 + state from the model potential 34 is smaller than the present value, while for the other states it is larger. There have been a

small number of experimental studies, with the two most precise being the recent PEPICO 25 In summary, good agreement is generally seen between the experimental results and the theoretical ones reported here, at least for the X 2 Σ 1/2 + and 2 Π 3/2 states; for the 2 Π 1/2 state, on the other hand, significant differences are seen with the ZEKE values; the agreement is sometimes slightly better between the present results and the PEPICO results.

B. Ion Mobilities

In the cases where experimental data exist, the calculated mobilities compared to experiment are given in Figure 2, and the results of the statistical comparison with the available experimental data are given in Table 3.

Since δ > 0 for all three weights of the interactions of 40 Ar + in Ne, the experimental values lie above the calculated data, as can be seen in Figure 2 The situation is essentially reversed for 84 Kr + in Ar. For 84 Kr + in Ar, the experimental values all lie below the calculated data by amounts that are statistically significant. This may indicate that the experimental error was greater than the 2% that was estimated.

For 132 Xe + in Ne, calculations for all three weights are virtually identical, and very similar for 132 Xe in Ar. The calculated results for Ne lie above the experimental values, which can be seen in Figure 2 for the ground state, but they lie slightly below for Ar. The differences are not quite significant for Ne and they are certainly not significant for Ar. There is little scatter in the data for either gas, and the disagreements do not depend on E/n 0 . Consequently, without more precise experimental data, it is not possible to ascertain what the distribution of Xe + SO states was in these experiments from the present comparison.

The precise isotope used was not identified in any of the experiments reported in Table 3. The theoretical prediction (from the first approximation of the two-temperature kinetic theory 49 ) is that the mobility is inversely proportional to the square root of the RG + -RG′ reduced mass. This indicates that 36 Ar + has a mobility in Ne that is 1.8 % higher than that of 40 High-level ab initio calculations have been employed to generate interaction potentials for RG + -RG′ complexes, where RG′ is lighter than RG. For the lowest 2 P state of RG + , 2 Σ + and 2 Π states result, and these two curves were calculated. We combined these, using relatively standard expressions that assume that the spin-orbit splitting is the same as in the isolated atom and does not vary with R. For the most part, very good agreement is seen with the most precise experimental data, suggesting that the potentials are reliable. That the spin-orbit potentials are derived from ab initio potentials using the same approach, suggests that the reliability of the curves should be similar for all of the states. However, it is seen that we obtain extremely good agreement with the experimental spectroscopic results for some of the states, but not all of them; even when all three states are studied in the same experiment. Additionally, even though there is occasionally better agreement with one study compared to another, this statement does not tend to cover all three spinorbit states, and hence it is difficult to come to a definitive conclusion regarding the reliability of the experimental data, or indeed the present results. That said, the fact that the present data ought to be internally consistent and of a comparable reliability for each state, coupled with the fact that very good agreement is seen with experimental data for some of the states of most of the species where data exist, we believe leads to the conclusion that the experimental data in some cases needs to be confirmed.

We note that, in principle, charge transfer states and their coupling to the lower states ought to be included in a fuller treatment of these states (see for example, refs 27,28 and 29). In the present case, as before, 1,2 we did not find evidence of any deviation from single-reference behaviour for the lowest 2 Σ + and 2 Π states. Inclusion of the excited states, necessitating a multireference, multiroot approach with spin-orbit coupling, might be an area for future exploration.

On this occasion, the comparison with mobility data has proven to be less useful than usual. This is primarily due to the limited systems and temperatures for which they were available. In addition, the experimental errors were usually too large to allow us to be conclusive about the population of the spin-orbit states; although the evidence appeared to suggest it was the ground state ( 2 P 3/2 ) that was the most likely -a conclusion that would be consistent with our earlier work on Ne + in He. 20 Ne, 40 Ar, 84 Kr, 131 Xe and 222 Rn. Here, the symbols represent the usual spectroscopic quantities: R e is the equlibrium bond length, D e is the depth of the potential, D 0 is the energy between the zero-point and the asymptote, ω e is the harmonic vibrational frequency, ω e x e is the anharmonicity constant, B 0 is the equilibrium rotational constant of the v = 0 level, D 0J is the first distortion constant for the v = 0 level, σ is the crossing point of the potential, when the potential = 0, and ∆G 1/2 and ∆G 3/2 are the first two vibrational separations.

Figure 1 :

 1 Figure 1: Interaction energy (E int ) curves for Ne-Ar + , relative to the energy of the lower spin-orbit asymptote Ne + Ar + ( 2 P 3/2 ).

Figure 2 :

 2 Figure 2: Comparison of calculated mobilities for RG + ( 2 P 3/2 ) to available experimental data. The lines represent the calculated data, and the data points with error bars represent the experimental data. The latter are the smoothed versions taken from refs. 4 and 5; see text for original references.E/n 0 is given in Td, and K 0 is given in cm 2 V -1 s -1 .

Figure 2 F

 2 Figure 2

  for R e and D e. There are no experimental data to which to compare. The R e values from refs. 34 and 35 are also in good agreement with each other and those obtained herein, but the D e values differ markedly, with those of Brunetti et al. 35 being in very good agreement with the present value, whereas the value of ref. 34 is significantly higher. For the 2 Π 1/2 state, similar comments apply.

For Ne-Kr + and Ne-Xe + , for the R e values there is again good agreement between the present work and the model potential results

34

of all three states, but the D e values are significantly higher than the present ones. There are no previous values for Ne-Rn + . For Ar-Kr + (X 2 Σ 1/2 + ), good agreement with ref. 34 is again obtained, but the D e value is somewhat higher in that work. The photoionization value of 4260±120 cm -1 from ref. 21 is in reasonable

  value is lower.There have been a few studies on Ar-Xe + , with the present D e value being in good agreement with most of the experimental values for the X 2 Σ 1/2 + state. Again, the agreement with the most precise experiments, the PEPICO

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

  and ZEKE 29 studies. For the X 2 Σ 1/2 + state, we see that there is good agreement between the experimental values and the present calculated ones. For the 2 Π 3/2 state, similar comments hold. It is noticeable that the D 0 value from the PEPICO study appears to be in slightly better agreement with the present value for the 2 Π 3/2 state than is the ZEKE study, while both the PEPICO and ZEKE studies yield virtually identical values for the X 2 Σ 1/2

+ state. For the 2 Π 1/2 state, we see that there is very good agreement between the PEPICO study and the present values for both ω e and D 0 ; the agreement is not so good with the ZEKE study, where significantly larger D 0 and ω e values are obtained. It is unlikely that this effect is attributable to isotope effects.

  and statistical states are not significantly different. If it is assumed that the experimental error30 of 2% is a true indication of the precision of the data, then the suggestion here is that there may have been a mixture of Ar + ( 2 P 3/2 ) and Ar + ( 2 P 1/2 ) present. For all three states, χ is only slightly above |δ|, which means that the disagreements are not significantly different at various values of E/n 0 and that the experimental values have little scatter.

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

for the ground state calculations. It can be seen that |δ| > 1 for the ground state, suggesting that there is a statistically-significant difference between the experimental data and those calculated for this state. The values of |δ| < 1 for the other two weightings indicate that the differences between the data and the values calculated for the excited

  84 Kr + (56.987 %). Although the abundances make this highly unlikely, it would increase the value of δ in Table III by 0.59 and bring the experimental and ground-state values into nearly perfect agreement (on average). The worst-case predictions for Xe + are too small to bring about disagreement between the calculated and measured mobilities in Ne or Ar, since the experimental uncertainties are already so large.

	IV. Concluding Remarks
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Ar + , which (in the unlikely event that this isotope only was monitored) would lead to even larger disagreements between the calculated and measured values. The best-case situation would be for Kr + in Ar if experiments measured the mobility of 78 Kr + (0.355 % natural abundance), which lies 1.2 % above that for

Table 1 .

 1 Parameters of the RG + -RG′ non-SO interaction potentials, with numbers in parentheses denoting the power of 10. Distances are in Å, and all other values are in cm -1 . Isotopes employed are:

	1

Table 2 .

 2 Parameters of the RG + -RG′ SO interaction potentials, with numbers in parentheses denoting the power of 10. Distances are in Å, and all other values are in cm -1 . See Table1for definition of terms, and isotopes used. Two values obtained in that work, arising from inclusion of two different sets of previous data in fit to potential.

	Page 11 of 20			Molecular Physics
	System State Ref. Ne-Ar + 2 Σ + This Work 2.10 2.42 917.2 848.7 σ R e D e D 0	ω e 139.8 6.10 0.211 2.16(-06) 127.6 115.3 ω e x e B 0 D J0 ∆G 1/2 ∆G 3/2
		36	2.41 a	810 b
	2 Π	This Work 2.67 3.06 262.7 232.7	61.4	4.25 .131	2.82(-06) 53.0 44.2
	Ne-Kr + 2 Σ +	This Work 2.35 2.71 608.5 562.5	93.9	4.12 0.139 1.37(-06) 85.7 77.4
		36	2.56 a	< 0 b,c
	2 Π	This Work 2.86 3.27 222.0 197.7	50.1	3.35 0.095 1.65(-06) 43.4 36.5
	Ne-Xe + 2 Σ +	This Work 2.66 3.05 412.4 378.8	68.8	3.26 0.102 1.02(-06) 62.3 55.7
	2 Π	This Work 3.11 3.53 177.6 157.5	41.2	2.81 0.076 1.22(-06) 35.6 29.8
	Ne-Rn + 2 Σ +	This Work 2.76 3.16 373.7 343.2	62.2	2.95 0.090 8.61(-07) 56.4 50.4
	2 Π	This Work 3.20 3.63 164.0 145.7	37.8	2.56 0.068 1.04(-06) 32.7 27.4
	Ar-Kr + 2 Σ +	This Work 2.20 2.59 5459.2 5347.0 226.2 2.09 0.093 6.54(-08) 222.0 217.9
		36	2.62 a	3000 b
	2 Π	This Work 3.03 3.50 630.3 602.1	57.2	1.56 0.050 1.68(-07) 54.1 50.9
	Ar-Xe + 2 Σ +	This Work 2.53 2.91 2452.9 2386.7 133.6 2.03 0.064 6.22(-08) 129.6 125.5
	2 Π	This Work 3.27 3.75 518.8 495.4	47.4	1.30 0.039 1.10(-07) 44.8 42.2
	Ar-Rn + 2 Σ +	This Work 2.67 3.06 1904.3 1851.9 106.0 1.69 0.053 5.43(-08) 102.6 99.2
	2 Π	This Work 3.36 3.84 489.6 468.3	43.3	1.15 0.033 8.48(-08) 41.0 38.7
	Kr-Xe + 2 Σ +	This Work 2.52 2.93 4960.4 4767.3 128.4 0.78 0.032 8.15(-09) 126.9 125.3
	2 Π	This Work 3.34 3.85 696.7 678.5	36.9	0.58 0.019 1.95(-08) 35.8 34.6
	Kr-Rn + 2 Σ +	This Work 2.65 3.06 3615.8 3561.1 109.9 0.84 0.029 8.63(-09) 108.2 106.5
	2 Π	This Work 3.43 3.94 659.6 641.9	35.7	0.58 0.018 1.82(-08) 34.6 33.4
	Xe-Rn + 2 Σ +	This Work 2.69 3.17 7240.7 7183.2 115.3 0.38 0.020 2.52(-09) 114.5 113.8
	2 Π	This Work 3.53 4.08 902.5 885.8	33.5	0.36 0.012 6.64(-09) 32.7 32.0
	a MP2/6-31G*
	b G2 value.

c Calculated to be unbound. URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph a Model potentials fitted to experimental data. b Photoionization Spectroscopy c Emission Spectroscopy d ZEKE Spectroscopy e PEPICO Spectroscopy f REMPI-PES Spectroscopy g URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 3 .

 3 Statistical comparison of calculated and experimental mobilities of RG + ions in RG′. Here T is the gas temperature (in K) used in the experiments whose smoothed values are reported in the reference cited, E/n 0 is the ratio of the electrostatic field strength to the gas number density (in Td, where 1 Td=10 21 V m 2 ), N is the number of data points, %Exp is the maximum inaccuracy claimed by the experimenters, %Calc is the maximum error in the calculations report here, Weight indicates the weighting used for the cross sections (see text), and δ and χ are the statistical quantities discussed in the text.

	RG +	RG′ T	Ref. E/n 0	N %Exp %Calc Weight	δ	χ
	40 Ar + Ne 77 5	6-20	5 2.0	0.1	Ground	1.15 1.16
							Excited	0.29 0.32
							Statistical 0.86 0.86
	84 Kr + Ar 300 4	25-150 9 2.0	0.1	Ground	-1.80 1.81
							Excited	-1.72 2.61
							Statistical -1.53 1.69
	132 Xe + Ne 300 5	15-140 10 5.0	0.2	Ground	-0.98 0.98
							Excited	-0.94 0.95
							Statistical -0.97 0.97
	132 Xe + Ar 294 5	3-90	6 7.0	0.2	Ground	0.15 0.17
							Excited	0.13 0.21
							Statistical 0.11 0.16
				90-220 3 7.0	0.2	Ground	-0.13 0.14
							Excited	-0.47 0.51
							Statistical -0.26 0.27

Acknowledgment

The work of LAV was supported by the US National Science Foundation under grant CHE-0718024. TGW is grateful for computing time at the NSCCS, funded by the EPSRC (UK).