Angel Luis Esteban
email: angel.esteban@ua.es

Ernesto Diez
email: ernesto.diez@uam.es

Maria P Galache

Jesus Sanfabian

Jorge Casanueva

Rubén H Contreras

J Angel

L Esteban

Ernesto Díez

Jesus San Fabián

Angel L Esteban Elúm

Vibrational contributions to vicinal proton-proton coupling constants 3JHH

Keywords: Vicinal coupling constants, Vibrational contributions, Numerical precision, DFT/B3LYP vicinal coupling constants, vibrational contributions, numerical precision, DFT/B3LYP

Introduction

The development of a practical procedure for accurate prediction of vicinal proton-proton coupling constants 3 HH J (written hereafter J) is an objective of great importance for investigating the conformational behaviour in solution of organic and biochemical compounds, as well as biopolymers, because of the strong dependence of J on the dihedral angle φ between the coupled protons [1][2]. After the angular dependence, the next factor that determines the magnitude of J couplings is the effect of the substituents [3][4][5][6][7]. In addition to these two important factors there is a series of secondary factors that influence the magnitude of the J couplings such as changes in local geometry and contributions from torsional vibrations [3,[8][START_REF] Díez | Contributions to the Vicinal H-H Couplings from Torsional Vibrations[END_REF][START_REF] Gil | [END_REF][11][START_REF] Esteban | Vibrational contributions to the spin-spin coupling constants[END_REF]. Predictions of J couplings have been carried out, in general, using analytical expressions empirically parameterized of limited precision [6,[START_REF] Haasnoot | [END_REF][START_REF] Altona | Encyclopedia of NMR[END_REF][START_REF] Navarro-Vázquez | [END_REF]. State-of-the-art theoretical methods are the natural alternative to the use of empirical equations for predicting accurately J couplings. A procedure to predict J from DFT/B3LYP calculations using a moderate sized basis set was developed recently [16]. For a dataset of 163 couplings the σ deviation between predicted pre J and experimental exp J couplings amounted 0.17 Hz against the larger σ deviation of 0.31 Hz reported for the best empirically parameterized extended Karplus equation [START_REF] Altona | Encyclopedia of NMR[END_REF]. J couplings were calculated at the optimized equilibrium geometries for a series of selected molecules but vibrational contributions were not estimated in Ref. [16]. In this paper, these contributions are calculated for a representative dataset of 84 J couplings in order to determine their effect on the prediction of vicinal J couplings.

In recent years, there has been a considerable interest on vibrational contributions to nuclear spin-spin coupling constants [17][18][START_REF] Ruden | Calculation of NMR and EPR Parameters[END_REF] which are typically of the 4 same order of magnitude as the differences between the coupling constants calculated theoretically and those observed experimentally. Therefore, these contributions should always be computed when the performance of computational methods is analyzed and when the temperature dependence of coupling constants or secondary isotope effects is studied [START_REF] Ruden | [END_REF].

In order to calculate vibrational contributions to coupling constants the first and diagonal second derivatives of J with respect to normal coordinates must be computed.

No analytical implementation exists for evaluating these derivatives, therefore, some numerical procedure must be used [21]. In principle, first and second derivatives of J with respect to normal coordinates k Q can be obtained from three values of J computed for the equilibrium geometry e J and for positively ( )

k Q δ +
and negatively ( )

k Q δ -
displaced geometries along the normal coordinate k Q . This procedure, that uses one δ value and requires the computation of two k J values, does not provide the computational accuracy necessary for our purposes. Therefore, the performance at the DFT/B3LYP level of procedures that use two or three δ values and require computing four or six k J values have been analyzed aiming at getting computational accuracies of 0.01 Hz for vibrational contributions vib J to J and smaller than 0.05 Hz/amu 1/2 Å and 0.10 Hz/amu Å 2 , respectively, for first and second derivatives of J with respect to k Q .

Vibrational contributions

At the second-order perturbation level, a vibrational averaged <J> coupling for a near rigid molecule with low amplitude vibrations is given by [22][23],

( )

2 e k k k k k J J J Q J Q ′ ′′ < >= + < > + < > ∑ , (1)
where k Q 's are the vibrational normal coordinates, e J is the J value for the equilibrium configuration, and k J ′ and 2 k J ′′ are the first and diagonal second derivatives of J with

respect to the normal coordinate k Q , ( ) ( ) 2 2
, 2

k k k k e e J J Q J J Q ′ ′′ = ∂ ∂ = ∂ ∂ . (2)
The average values of k Q and 2 k Q in thermal equilibrium at temperature T, using the harmonic oscillator partition function, are given by [24],

∑         Φ - >= < l B l l kll k k T k Q 2 coth 4 2 ω ω ω h h , ( 3
) 2 coth 2 2 k k k B Q k T ω ω   < >=     h h , (4)
where B k is the Boltzmann constant, k ω is the vibrational angular frequency of k-th normal mode and kll Φ 's are the cubic anharmonic normal coordinate force constants defined as 3) and ( 4), the first and second order terms in equation

( ) e l k kll Q Q E 2 3 ∂ ∂ ∂ = Φ . ( 5
(1), k k J Q ′ < > and 2 k k J Q ′′ <
> , correspond, respectively, to anharmonic and harmonic vibrational contributions to <J> from the k-th normal coordinate k Q . Therefore, a vibrationally averaged coupling can be written as

har anh e J J J J + + >= < , (6) with
, anh k k anh k k k J J Q J ′ = < > = ∑ ∑ , ( 7
) 2 , har k k har k k k J J Q J ′′ = < > = ∑ ∑ . ( 8
)
Combining contributions anh J and har J yields the vibrational contribution vib J to <J>, ( )

, , vib anh har anh k har k k J J J J J = + = + ∑ . (9)
On the other hand, when T = 0 K the coth functions in equations ( 3) and ( 4 Therefore, thermal contributions from low frequency vibrations can be important as compared with the respective zero point contributions.

T vib har anh e T J J J J J ∆ + + + >= < 0 0 , ( 10
)
T har T anh T vib J J J ∆ + ∆ = ∆ . ( 11

Computational details

All computations have been performed with the G03 system of programs [START_REF] Frisch | Gaussian 03, Revision D.02[END_REF] at the Density Functional Theory (DFT) level using the hybrid functional B3LYP [START_REF] Becke | [END_REF][27][28].

Several outputs were modified in order to increase the number of printed digits.

The harmonic frequencies k ω and the cubic force constants kll Φ in equations ( 3) and ( 4) were computed at different levels of theory as the hybrid force field B3LYP/6-31+G(d,p)//B3LYP/6-31G(d) proposed by Carbonier et al [29] to achieve the best compromise between computer time and quality of results. These authors tested the performance of B3LYP density functional for computing harmonic and anharmonic frequencies for a set of near rigid molecules. That study revealed that the relatively cheap 6-31+G(d,p) basis set performs a very good job for harmonic frequency calculations and that B3LYP anharmonicities are in close agreement with the reference values irrespective of the used basis set.

Equilibrium geometries were calculated with high numerical accuracy (optimization criterion Opt=VeryTight and numerical grid Int=Ultrafine) at the B3LYP/6-31+G(d,p) level obtaining harmonic frequencies ω k and normal coordinates k Q . Cubic force constants kll Φ were computed at the B3LYP/6-31G(d) level [29].

J couplings were computed for equilibrium geometries and at positively ( )

k Q δ + and negatively ( ) k Q δ - displaced geometries along k Q by using the DFT/B3LYP

Results and discussion

Computational precision

In a recent work [START_REF] Esteban | Vibrational contributions to the spin-spin coupling constants[END_REF], first k J ′ and second 2 k J ′′ derivatives of J with respect to k Q , equation (2), were computed by numerical differentiation using the standard three points formula,

( ) ( ) 2 k k k J J J δ δ δ ′ = + - -     , (12)
( ) ( ) being the optimal value to be used of 0.1. However, the precision of values calculated in this way for k J ′ , k J ′′ and vib J was not then analyzed.

2 2 2 k k k e J J J J δ δ δ ′′ = + + --     , (13) where (
The numerical precision of the k J ′ , k J ′′ and vib J values depends on the precision ε of calculated J couplings and on the δ value used in the finite difference procedure. For a given ε, there is an optimum δ value such that numerical errors due to ε, which decrease with increasing δ, are equal to truncation errors due to the neglect of higher derivatives, which increase with increasing δ [21]. Unfortunately, neither the ε error nor the optimum δ value to be used for each normal coordinate is known in principle.

Therefore, in order to determine the precision of k J ′ , k J ′′ and vib J more than three points must be used for calculating the k J ′ and 2 k J ′′ derivatives.

F

2 3 4 T k ek k k k k J J J J J J δ δ δ δ δ ′ ′′ ′′′ ′′′′ = + + + + , ( 14
) provide ek J ∆ values, ek ek e J J J ∆ = -, (15)
which are smaller than 10 Taking ε = 10 -7 Hz, the precision of k J ′ and k J ′′ , calculated using equations [START_REF] Esteban | Vibrational contributions to the spin-spin coupling constants[END_REF] and [START_REF] Haasnoot | [END_REF], respectively, should be better than 0.01 when the δ > 10 -5 for the former, and δ > 2×10 -3 for the latter. k J ′ values thus calculated with δ equal to 0.0001, 0.001 or 0.01 differ from the corresponding fitted ones less than 0.03, 0.01 and 0.003, respectively.

For δ = 0.1 the differences are larger, but smaller than 0.4, due to the contributions from the cubic terms which are not taken into account in equation ( 12). For symmetrical vibrations, k J ′′ calculated by using equation ( 13) for δ equal to 0.001 or 0.01 differ less than 0.3 and 0.01, respectively, from the corresponding fitted ones. For δ = 0.1 the differences are larger than for δ = 0.01, but smaller than 0.3, due to the contamination from quartic terms. For antisymmetrical vibrations there is, in addition, a spurious are eliminated, consequently, the maximum error for the averaged <J> couplings decreases below 0,04 but the maximum errors for k J ′ and k J ′′ are close to those found when only the δ value of 0.10 is used.

The maximum computational error for the averaged <J> couplings, given in tables 1 to 5 as pre J , is, probably, smaller than 0.01 Hz since, though the maximum errors for the k J ′ and k J ′′ when using the three δ values of 0.01, 0.05 and 0.10 are much smaller than when using only the two δ values of 0.01 and 0.05, differences between the averaged <J> couplings at 300 K calculated for both sets of δ values are, in general, smaller than 0.01 Hz. For bromocyclopropane such differences are a bit larger in absolute value, amounting to -0.04 Hz for coupling J 1,2' , due to the fact that when, for the lowest vibrational frequency calculated as 213. approximation which assumes the molecule exists only in the most stable rotational isomers. For mono and 1,1-disubstituted ethanes, with monovalent and methyl substituents, all three rotational isomers are equivalent and RIS J is given by ( )

1 3 RIS T G G e e e J J J J + - = + + , ( 16
)
where X e J , X=T, G + and G -, are vicinal couplings for the equilibrium conformation defined in figure 1.

Figure 1 about here

Similarly, when vibrational contributions are taken into account, such averaged coupling J is given by ( )

1 3 T G G J J J J + - = + + , (17)
with

X har X anh X e X J J J J + + = , (18)
where the anharmonic X anh J and harmonic X har J vibrational contributions are given by equations ( 7) and ( 8).

X J couplings were calculated for 6 mono-and 5 1,1-disubstituted ethanes selected from those used by Altona et al [START_REF] Altona | Encyclopedia of NMR[END_REF] to derive a substituent parameter λ scale based on the relationship between J couplings and substituent electronegativities.

Calculated couplings, cal J , are overestimated but predicted pre J couplings obtained by scaling cal J with an optimized factor f = 0.838 are very close to experimental values, exp J , i.e. J , 2,3 J and 2 ,3 J ′ ′ , are smaller than those predicted previously [16], between 0.18 and 0.60 Hz, while the corresponding differences for couplings of near type skew S J , 145

J J J J = + + ∆ . ( 20
HCCH ϕ ≈ , 1,2
J ′ and 2,3 J ′ , are smaller in magnitude, between -0.17 and 0.17 Hz. Harmonic contributions X har J , 0 in each molecule are very close for the two couplings 1,2 J and 1,2 J ′ as well as for the three couplings 2,3 J , 2,3 J ′ and 2 ,3 J ′ ′ , being the values for the latter couplings, between 0.52 and 0.84 Hz, larger than for the former ones, between 0.28 and 0.54 Hz, except for X=Li.

Table 5 about here

Results obtained for three-membered heteroatomic rings are given in table 5.

Like data reported in table 4, couplings of type cis C J , ( 2,3 J and 2 ,3 J ′ ′ ), are smaller than those predicted previously, between 0.17 and 0.38 Hz, but couplings of type skew S J , ( 2,3 J ′ ), are larger, between 0.01 and 0.18 Hz. Harmonic contributions

+ = , (21)
with a correlation coefficient R of 0.97 and a standard deviation σ of 0.10 Hz. The range for the T har J , 0 contributions, between 0.94 and 1.11 Hz, is smaller than that for G har J , 0 , which is between 0.10 and 0.52 Hz. couplings are small in three-membered rings, between -0.07 and 0.02 Hz, being a bit larger in magnitude in norbornane type molecules, between -0.06 and -0.14 Hz.

+ - = , (23)
Vibrational contributions 300, X vib J ∆ at 300 K show a rough positive correlation with X e J , 300, 0.12 0.047

X X vib e J J = + , (26)
with R = 0.68 and σ = 0.19 Hz.

Overall deviations between exp J and pre J couplings, calculated including the vibrational contributions, see last column of tables 1 to 5, are of similar magnitude than those, given in parentheses, obtained previously calculating pre J from equilibrium values X e J [16]. Fittings to equation [START_REF] Ruden | Calculation of NMR and EPR Parameters[END_REF] The accuracy of vibrational contributions calculated by means of equations ( 1) to ( 9) depends, besides on the computational accuracy, on the accuracy of the calculated k J couplings, vibrational frequencies k ω , and cubic force constants kll Φ . Equations ( 3) and ( 4) show that, the remaining factors being equal, the errors in

k Q < > and 2 k Q < >
and, therefore, in vib J increase when k ω decrease. The used equations are, in general, a good approximation for near rigid molecules with low amplitude vibrations. However, the approximation will be less good for flexible molecules with large amplitude vibrations of low frequency, such as the torsional vibrations of methyl groups. The estimation of vibrational contributions to J from these vibrations requires the determination of the energy levels and wave functions for the appropriate Hamiltonian.

In the case of a molecule with one methyl group, a Hamiltonian for hindered rotation should be applied [8][START_REF] Díez | Contributions to the Vicinal H-H Couplings from Torsional Vibrations[END_REF]. Nevertheless, the deviations for the methyl couplings predicted in this work are similar to those for couplings of near rigid molecules supporting the conclusion that the contribution to the averaged coupling of one methyl group from the torsional vibration is small [8][START_REF] Díez | Contributions to the Vicinal H-H Couplings from Torsional Vibrations[END_REF]. It should be realized that for a methyl group the three rotational isomers are equivalent with equal populations of 1/3, therefore, equation ( 17) may be a reasonable approximation. On the other hand, when populations of rotational isomers are different, such as in the case of 1,2-disubstituted ethanes, the situation is much more complicated and the contributions from the torsional vibrations can be important [START_REF] Díez | Contributions to the Vicinal H-H Couplings from Torsional Vibrations[END_REF][START_REF] Gil | [END_REF][11][START_REF] Guilleme | Contributions to the Vicinal H-H Couplings from Torsional Vibrations in Ethane Derivatives[END_REF][31]. The problem of the conformational analysis is then involved and the energy levels and wave functions for the appropriate Hamiltonian for hindered The error of a predicted coupling pre J depends, among other factors, on the error of the calculated coupling for the equilibrium configuration e J . The error of e J is determined by both the accuracy of the molecular geometry used for the equilibrium configuration and of the theoretical method used to calculate the coupling constant. In the second numerical column of tables 1 to 5 differences

X e J ∆ between the X e J
couplings predicted this work and those predicted in reference [16] are given. Except for molecules with Br atoms, these differences are due to small changes in local geometry between the molecular geometries used in these work, optimized with the 6-31+G(d,p) basis set, and those used in [16], optimized with the 6-31G(d) basis set.

X e J ∆ values
provide an idea about the importance of errors in the calculated e J values from deviations of the optimized molecular configuration with respect to the true equilibrium configuration. In table 1 the only ethane derivatives which show X e J ∆ values larger than 0.1 appear for molecules with F atoms, the most electronegative substituent, reaching 0.32 Hz for T J of 1,1-difluoroethane, and for 1,1-dibromoethane. With respect to the latter molecule it should be noted that the 6-311G(d,p) basis set was used for the Br atoms in [16] for both the geometry optimization and for calculating X e J . Also,

X e J ∆
values larger than 0.1 appear in table 4 for the two cyclopropanes with substituents Br, reaching -0.31 Hz for the coupling 2 ,3 J ′ ′ , and Li, the most electropositive substituent.

Conclusions

One of the prerequisites to develop a procedure for a reliable and accurate prediction of vicinal J couplings is the feasibility of calculating vibrational contributions T vib J with a good computational accuracy which we target in 0.01 Hz. The present work

shows that to reach this aim from equations ( 1) to ( 5 The importance of vibrational contributions T vib J computed for the J couplings of near rigid molecules and methyl protons shows that these contributions must be included for a reliable prediction of vicinal couplings and, therefore, for a more accurate application of J couplings in conformational analyses. However, equations ( 1) to (5) used in this work are not applicable, in general, when dealing with the important practical situation of flexible molecules with large amplitude vibrations involved in conformational equilibriums. In these cases the energy levels and wave functions for the large amplitude vibrations must be calculated from the appropriate Hamiltonian in order to estimate the corresponding vibrational contributions and the conformer populations.

Appendix

Two points have been analyzed after completion of this work. The obtained results are presented here for the sake of completeness but they do not affect to our conclusions.

The first point refers to the temperature dependence of <J> due to centrifugal distortions. This is small for one and two bond coupling constants [45]. The contribution J cen to <J> from this source is given by [24],

∑ = k cen k k cen Q J J ' , (27)
∑ Ι = α αα αα ω ) ( ) ( 2 2 e k k B cen k a T k Q , (28)
where ) e ( αα Ι is the moment of inertia at equilibrium geometry along the principal axis α and ) (αα k a is the derivative of Ι αα with respect to the coordinate Q k evaluated at equilibrium. In equation ( 28) a classical average is used for the rotational part.

The calculated J cen contributions are smaller in magnitude than 0.01 Hz for cyclic molecules and than 0.03 Hz for ethane derivatives.

Another interesting point refers to the method chosen to get a computational accuracy of 0.01 Hz for vibrational contributions J vib based on the exhaustive application of fittings to polynomes of second, third and fourth degrees in δ of the calculated J k values for δ equal to 0.01, 0.05 and 0.10. These fittings were carried out without including the calculated J e for the equilibrium conformation value in data because differences ∆J ek between fitted values J ek and J e are often several orders of magnitude larger than the precision ε for calculated J k values. Alternatively, differentiation methods based on using interpolation formulas or arisen from Taylor series approximations could be used, in principle. The Romberg's procedure, based on the second approach, has been used by Champagne and coworkers for numerical evaluation of hyperpolarizabilities [46] and hyper-Raman spectra [47]. Increasingly large values of the step size δ=2 i-1 λ, i=1,2,3… are used. The elimination of the higher order contamination terms is then performed by applying the iterative formula J values are taken from [38] for X = H, [39] for X = Cl, Br, and CN, and [40] for X = Li and F. J values are taken from [41] for X = NH and Y = H, and for X = O and Y = H, [42] for X = O and Y = Cl, [43] for X = PH and Y = H, and [44] for X = S and Y = H.

K , 2 , 1 , 1 4 4 1 1 1 = - - = + - - j P P P j i j i j j i j , ( 29

globally considered, 300 vibJ

 300 contributions amount to about 7% of the e J values and the correlation coefficient between the 300 vib J and e J is of only 0.68 with a σ deviation of 0.20. The first and diagonal second derivatives of J with respect to each normal coordinate k Q , required to estimate the vibrational contributions, have been obtained from six k J values computed for molecular geometries positively normal coordinate k Q and using for δ the values 0.01, 0.05 and 0.10. The computational precision of the results obtained when using one, two and three δ values is analyzed.

) are equal to 1 ,

 1 and the zero point values are obtained for equations(7) to[START_REF] Díez | Contributions to the Vicinal H-H Couplings from Torsional Vibrations[END_REF], the zero point vibrational contributions to J, the corresponding ones at the zero point can be considered as the thermal vibrational contributions

  T and decreases with ω, and at T = 300 K for vibrational frequencies k ω of 1000, 300 and 100 cm -1 becomes, respectively, 1.017, 1.62 and 4.24.

-

  are J values computed for molecular geometries obtained increasing and decreasing, respectively, the equilibrium geometry in an amount k Q δ . It was found that the 2 k J ′′ derivatives show, sometimes, an excessive variation with δ

JJ

  ′′ , which amounts to about 0.5 for δ = 0.1. ′′ can be removed by selecting two δ values, instead of only one, and fitting the sets of four T k J values calculated for each normal coordinate k Q to either a second or third degree polynome. Differences of the so derived values for k J ′ and k J ′′ with respect to those provided by the fittings of eight T k J values to the four degree polynome, equation[START_REF] Altona | Encyclopedia of NMR[END_REF], will be taken hereafter as their computational errors. When one of the two selected values is equal to 0.1 and the fourT kJ values calculated are fitted to a second degree polynome the largest errors for k J ′ and k J ′′ amount to 0.32 and 0.24, respectively, due to contamination from cubic and quartic terms. However, when one of the two selected δ values is equal to 0.01 and the other one is either 0.001 or 0.0001 the errors are smaller than 0.003 since such very small. When the two δ values are 0.001 and 0.0001, the largest error for k J ′′ increases up to 0.3 due to the low precision ε of the calculated T k J couplings. However, k J ′ errors are smaller than 0.02. On the other hand, when the fourT kJ values calculated are fitted to a third degree polynome and one of the two selected δ values is equal to 0.1 the errors for the k J ′′ do no change but the errors for the k J ′ decreases below 0.03 since contributions from the cubic terms are taken into account.The precision of the k J ′ coefficients derived from the four T k J couplings computed for two δ values can be estimated, in principle, by comparison of the four k J ′ values which can be calculated: two by means of equation (12) and the other two from the second and third degree fittings. When the two selected δ values are 0.01 and 0.001 the differences among the four k J ′ values are smaller than 0.02. In this case, clearly, the contamination from cubic and quartic terms is very small.The precision of the k J ′′ coefficients cannot be estimated, in general, by selecting only two J values due to the spurious contributions equation (13) is applied and to the fact that the second and third degree fittings provide near the same value for each k J ′′ . For these fittings the odd terms k k J ′′ precision can be estimated, in principle by selecting three instead of two δ values. When these values are 0.001, 0.01 and 0.1 and the six calculated T k J couplings are fitted to a fourth degree polynome the k J ′ , k J ′′ , k J ′′′ and k J ′′′′ errors are smaller than 0.0001, 0.002, 0.001 and 0.2, respectively. These values increase to 0.02, 0.4, 173 and 3452, when the δ values are 0.0001, 0.001 and 0.01. An inadequate selection of δ values, 0.0001 is too small, provides anomalous values for k J ′′′′ and k J ′′′ which indicate larger errors in k J ′′ and k J ′ , respectively. The J ′′ values can be estimated comparing the results of the four, third and second order fittings for the three values of δ and the results of the second and third order fittings for each of the three possible pairs of δ values. C J cis coupling in cyclopropane was also studied in detail since this is a symmetrical ring molecule with degenerate vibrations. Calculated C e J is 10.08 Hz, and C k J couplings were calculated taking δ values of 0.001, 0.01 and 0.1, which provide precise k J ′ and k J ′′ values for only a part of normal vibrations. For other normal vibrations δ = 0.001 results too small and systematic deviations appear. Additional C k J values were calculated for δ values of 0.005, 0.05, 0.075, 0.175 and 0.250 and each set of sixteen C k J was fitted to both a sixth and a fourth order polynome. The differences between the results of these fittings are smaller than 0.01 and 0.03 for k J ′ and k J ′′ , respectively, showing the small effect of the fifth and sixth order terms. On the other hand, when only three δ values are used, the respective maximum deviations for the fourth degree fittings are 0.01 and 0.03 for δ values of 0.01, 0.05 and 0.10 but increase to 0.07 and 0.44 for δ values of 0.001, 0.01 and 0.10. On the basis of the above described analysis of the T J coupling of fluoroethane and the C J coupling of cyclopropane we decided to compute the vibrational contributions of all couplings for δ values of 0.01, 0.05 and 0.10. Results, multiplied by 0.838, are given in tables 1 to 5 for the 84 considered couplings. Deviations with respect to these results of those obtained when less than three δ values are used are considered as computational errors. Excluding the 5 couplings of bromocyclopropane (see below) the following maximum errors are found. When only the δ value of 0.10 is used with equations (12) and (13) the maximum errors for k J ′ and k J ′′ derivatives and vibrational <J> couplings are 0.83, -1,02 and -0.31, respectively. When only the two δ values of 0.01 and 0.05 are used the respective values decrease to 0.20, -0.26, and 0.01 due to the fact that the spurious contributions 2 ek J δ -∆ to k J ′′ are eliminated and the contamination from quartic terms is reduced. When the two δ values are either 0.10 and 0.05, or 0.10 and 0.01, the contamination from quartic terms remains but the

  couplings are small, they are between -0.04 and -0.02 Hz, compared with zero point vibrational contributions, between -0.10 and -0.05 Hz.

Figure 2

 2 Figure 2 about here

Figure 3

 3 Figure 3 about here

  different couplings in each molecule ranging from 0.1 for (X,Y ) = (O,Cl) to 0.7 Hz for (X,Y ) = (PH,H).

Figure 5

 5 Figure 5 about here Zero point anharmonic and harmonic contributions for couplings of the trans

Figure 6

 6 Figure 6 about here

  0.86 and σ = 0.07 Hz. On the other hand, thermal contributions for cis

  calculated in order to estimate the averaged J coupling and the populations of rotational isomers.

J

  ), at the DFT/B3LYP level, and using the G03 programs, six values of k J must be computed, in principle, for each normal coordinate k Q . When k J values are calculated for molecular geometries positively ( ) equilibrium geometry along the normal coordinate k Q the three δ values of 0.01, 0.05 and 0.10 constitute an appropriate set to be used. The computational error increases, in general, less than 0.01 Hz when only the two δ values of 0.01 and 0.05 are used but, in some cases, errors are introduced by contributions to T vib J from low frequency vibrations. Therefore, three δ values should be always used for these vibrations. On the other hand, although the precision of the first k J ′ derivative of J with respect to the normal coordinate k Q can be estimated from four k J calculated for two δ values, estimation of the precision of the second 2 k J ′′ derivative requires calculating six k J for three δ values. When only one δ value of 0.10 is used, the computational error of T vib J increases up to 0.3 Hz for some couplings. The actual errors of calculated vibrational contributions T vib J are, of course, larger than the computational accuracy due to several factors discussed above. Vibrational contributions 300 vib J are positively correlated with e J values for the equilibrium configurations being the correlation coefficient equal to 0.68. This fact explains, in part, the reason why the global agreement between predicted pre J and experimental exp J couplings is near equal when the predicted couplings are calculated than when they are calculated from only the e J values. The main contributions to 300 vib J are the zero point harmonic contributions 0, X har J which show positive correlations with the X e J being the correlation coefficients larger than 0.9.

with Page 25 of 41 URL

 41

  Figure 1.

Figure 2

 2 Figure 2

Figure 3 PageFigure 1 .

 31 Figure 3

Figure 2 .

 2 Figure 2. A monosubstituted cyclohexane with the substituent X in axial position (a) and in equatorial position (b).

Figure 3 .

 3 Figure 3. Norbornane type molecules.

Figure 4 .

 4 Figure 4. Monosubstituted cyclopropanes (a) and three-membered rings with one heteroatom (b).

Figure 5 .

 5 Figure 5. Zero point harmonic ( ) and anharmonic (∆) contributions for couplings of types trans and gauche versus e J values.

Figure 6 .

 6 Figure 6. Zero point harmonic ( and ) and anhamonic (∆ ) contributions for couplings of types cis ( )and skew ( ) versus e J values.

Table 2 about here

 2 Results obtained for the chair conformations of cyclohexane, cyanocyclohexane and chlorocyclohexane are displayed in table 2. These J couplings are analogous to those for methyl protons: AA J are of type trans like T J since 180 ≈

										φ	HCCH	, and AE J
	and EE J	are of type gauche like G J	+	and G J	-since	φ	HCCH	≈	60	. Comparison of
	vibrational contributions with those of the same type given in table 1 shows two small
	but remarkable differences: anharmonic contributions F o r contributions EE T vib J , are negative.	AA anh J , 0	are positive and thermal
		P
		e
		e r
		R
		e
				v i e
					w
							O n l
											y

Table 3 about here

 3 Results obtained for norbornane-type molecules shown in figure3are given in

table 3 .

 3 The predicted couplings of type cis C H 1 and H 4 , and bridge, H 7a and H 7s , hydrogen atoms are small. On the other hand, for the small couplings between apical and endo hydrogen atoms, the larger deviations found for 3 seem to be due to low accuracy of experimental values since the predicted values are close to those predicted for1,2n

	J ,	φ	HCCH	≈	0	, 2 ,3 x x J	and 2 ,3 n n J	, are

HCCH

ϕ ≈ , between apical, J in 2, 0.66 Hz, which only deviates -0.11 Hz.

Table 4 about here

 4 Results for some monosubstituted cyclopropanes are given in table 4. Like in norbornane type molecules, predicted couplings of type cis C

	J ,	φ	HCCH	≈	0	, 1,2

Table 1 .

 1 ′′ when ∆J ek is negligible. As the magnitude of ∆J ek increases the results for k Predicted a pre J and experimental b exp J methyl couplings for mono-and 1,1disubstituted ethanes, see figure1. Vibrational contributions to pre J : zero point anharmonic Theoretical values in this table correspond to computed values multiplied by 0.836, see equation 17. b exp J values are taken from the dataset constructed by Altona et al.[4]. The value for 1,1-dibromo ethane is the weighted mean of two measures[32][33]. c Differences between X e J couplings predicted in this work and those predicted in[16]

	)

a

33 of 41 URL: http://mc.manuscriptcentral.com/tandf/tmph

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

Table 2 .

 2 Predicted pre J and experimental a exp J couplings for some monosubstituted cyclohexanes, see figure 2. For notation see table 1.

	X			J	X	J ∆	X e	J	0, X anh	J	0, X har	J ∆	, T X vib	J	300, X vib	J	pre	J	exp	J	exp	-	J	pre
	H	J	1 ,2 a a	11.97 0.00 0.08 1.10 -0.05 1.12 13.10 12.85 -0.25(-0.07)
		J	1 ,2 a e	3.48 0.04 -0.04 0.24 0.08 0.28 3.77 3.61 -0.16(-0.10)
		J	1 ,2 e e	2.66 0.00 -0.09 0.22 -0.04 0.09 2.76 2.92 0.16( 0.04)
	CN	J	1 ,2 a a	11.19 0.02 0.07 1.04 -0.07 1.04 12.23 12.66 0.43( 0.60)
		J	1 ,2 a e	3.31 0.03 -0.04 0.23 0.08 0.28 3.59 3.53 -0.06(-0.01)
		J	1 ,2 e a	4.26 0.13 -0.06 0.18 0.12 0.24 4.49 4.27 -0.22(-0.18)
		J	1 ,2 e e	2.33 -0.06 -0.08 0.20 -0.05 0.07 2.40 2.69 0.29( 0.11)
	Cl	J	1 ,2 a a	10.73 0.03 0.08 1.09 -0.06 1.11 11.84 11.88 0.04( 0.34)
		J	1 ,2 a e	3.91 0.05 -0.07 0.27 0.04 0.23 4.14 4.08 -0.06(-0.08)
		J	1 ,2 e a	2.93 0.07 -0.08 0.11 0.09 0.13 3.06 2.97 -0.09(-0.12)
		J	1 ,2 e e	2.68 -0.01 -0.08 0.22 -0.06 0.08 2.76 3.03 0.27( 0.12)

a All the exp J values are taken from

[34]

.

Table 3 .

 3 Predicted pre J and experimental a exp J couplings for some norbornane type molecules, see figure3. For notation see table1.

	C			J	X	J ∆	X e	J	0, X anh	J	0, X har	J ∆	, T X vib	J	300, X vib	J	pre	J	exp	J	exp	-	J	pre
	1	J	2 ,3 x x	11.46 -0.02 -0.16 0.72 -0.14 0.42 11.88 12.22 0.34(-0.18)
		J	2 ,3 x n	4.02 0.00 0.05 0.60 0.06 0.70 4.72 4.62 -0.10( 0.28)
		J	2 ,3 n n	8.62 0.01 -0.14 0.46 -0.10 0.21 8.84 9.05 0.21(-0.25)
	2	J	2 ,3 x x	8.61 -0.01 -0.14 0.43 -0.05 0.24 8.85 9.38 0.53( 0.07)
		J	2 ,3 x n	3.48 0.00 0.07 0.53 0.03 0.63 4.10 3.91 -0.19( 0.16)
		J	2 ,3 n n	8.39 -0.01 -0.12 0.46 -0.06 0.28 8.67 9.02 0.35(-0.05)
		J	1,2 x	3.34 0.01 -0.11 0.15 -0.02 0.02 3.37 3.66 0.29( 0.06)
		J	1,2n	0.59 0.01 -0.03 0.10 0.01 0.08 0.66 0.55 -0.11(-0.07)
		J	1,7 s	1.45 0.00 -0.06 0.12 0.00 0.06 1.50 1.51 0.01(-0.05)
		J	1,7a	1.68 0.01 -0.07 0.12 -0.01 0.05 1.72 1.83 0.11( 0.03)
	3	J	2 ,3 x x	11.21 0.01 -0.16 0.62 -0.13 0.33 11.54 12.05 0.51(-0.04)
		J	2 ,3 x n	4.16 0.03 0.00 0.59 0.03 0.62 4.78 4.70 -0.08( 0.24)
		J	2 ,3 n x	3.96 -0.03 0.10 0.57 0.05 0.72 4.69 4.59 -0.10( 0.28)
		J	2 ,3 n n	8.57 0.00 -0.13 0.44 -0.11 0.20 8.77 9.12 0.35(-0.13)
		J	1,2 x	4.47 -0.01 -0.09 0.22 -0.05 0.08 4.55 4.72 0.17(-0.12)
		J	1,2n	0.56 0.01 -0.03 0.12 0.02 0.11 0.67 0.12 -0.55(-0.47)
		J	1,7 s	1.08 0.01 -0.04 0.10 0.01 0.07 1.15 1.19 0.04( 0.03)
		J	1,7a	1.54 0.00 -0.06 0.12 -0.01 0.05 1.59 1.64 0.05(-0.02)
		J	4,7 s	1.96 -0.01 -0.06 0.15 -0.01 0.08 2.04 2.10 0.06(-0.03)
		J	4,7a	1.42 -0.06 -0.05 0.10 0.00 0.04 1.46 1.60 0.14( 0.08)
		J	4,3x	4.10 0.00 -0.11 0.19 -0.03 0.05 4.15 4.30 0.15(-0.12)
		J	4,5 x	4.41 0.00 -0.09 0.24 -0.02 0.13 4.54 4.76 0.22( 0.00)
		J	4,3n	0.48 0.01 -0.02 0.12 0.03 0.12 0.61 0.12 -0.49(-0.39)
		J	4,5n	0.42 0.00 -0.03 0.12 0.02 0.12 0.54 0.00 -0.54(-0.46)
	a exp J values are taken from [35] for 1, [36] for 2, and [37] for 3.

Table 4 .

 4 Predicted pre J and experimental a exp J couplings for some monosubstituted cyclopropanes, see figure4a. For notation see table1.

	X			J	X	J ∆	X e	J	0, X anh	J	0, X har	J ∆	, T X vib	J	300, X vib	J	pre	J	exp	J	exp	-	J	pre
	H	J	1,2	8.45 -0.02 -0.15 0.54 -0.01 0.38 8.83 8.97 0.14(-0.17)
		1,2 J ′	5.07 -0.02 0.01 0.52 0.00 0.53 5.61 5.58 -0.03( 0.08)
	F	J	1,2	5.69 0.02 -0.14 0.31 -0.01 0.17 5.85 5.89 0.04(-0.23)
		1,2 J ′	2.22 -0.03 -0.02 0.28 0.02 0.28 2.50 2.39 -0.11(-0.04)
		J	2,3	9.96 -0.08 -0.20 0.82 -0.03 0.59 10.56 10.80 0.24(-0.04)
		2,3 J ′	7.00 0.03 -0.04 0.74 -0.01 0.70 7.70 7.70 0.00( 0.17)
		2 ,3 J ′ ′	11.32 -0.06 -0.26 0.84 -0.04 0.54 11.86 12.01 0.15(-0.28)
	Cl	J	1,2	6.81 0.01 -0.14 0.36 -0.01 0.21 7.03 7.01 -0.02(-0.33)
		1,2 J ′	3.43 0.02 -0.01 0.35 0.02 0.36 3.79 3.58 -0.21(-0.10)
		J	2,3	9.48 -0.08 -0.16 0.76 -0.02 0.58 10.06 10.26 0.20(-0.06)
		2,3 J ′	6.56 -0.08 -0.03 0.70 -0.01 0.67 7.23 7.14 -0.09(-0.03)
		2 ,3 J ′ ′	10.23 -0.08 -0.19 0.76 -0.03 0.54 10.77 10.58 -0.19(-0.55)
	Br	J	1,2	7.00 0.04 -0.13 0.38 -0.01 0.23 7.23 7.13 -0.10(-0.38)
		1,2 J ′	3.70 0.01 -0.01 0.37 0.02 0.39 4.08 3.80 -0.28(-0.18)
		J	2,3	9.31 -0.18 -0.17 0.74 -0.02 0.54 9.85 10.16 0.31(-0.08)
		2,3 J ′	6.28 -0.28 -0.04 0.68 -0.01 0.63 6.91 7.01 0.10(-0.07)
		2 ,3 J ′ ′	9.92 -0.31 -0.19 0.73 -0.03 0.52 10.44 10.45 0.01(-0.59)
	Li	J	1,2	10.39 0.00 -0.04 0.67 -0.01 0.63 11.02 11.31 0.29( 0.09)
		1,2 J ′	8.13 -0.08 0.11 0.73 0.01 0.85 8.98 9.12 0.14( 0.26)
		J	2,3	7.13 -0.12 -0.15 0.31 -0.04 0.13 7.25 7.71 0.46(-0.11)
		2,3 J ′	3.03 -0.09 -0.04 0.30 -0.04 0.22 3.25 3.17 -0.08(-0.19)
		2 ,3 J ′ ′	5.79 -0.12 -0.07 0.29 0.00 0.21 6.00 6.42 0.42( 0.05)
	CN	J	1,2	8.04 0.01 -0.13 0.51 0.07 0.45 8.49 8.43 -0.06(-0.24)
		1,2 J ′	4.68 0.02 -0.01 0.48 0.05 0.52 5.20 5.12 -0.08( 0.09)
		J	2,3	8.53 -0.06 -0.14 0.59 -0.06 0.38 8.92 9.18 0.26(-0.09)
		2,3 J ′	6.41 -0.06 0.00 0.64 -0.04 0.60 7.02 7.08 0.06( 0.10)
		2 ,3 J ′ ′	9.08 -0.03 -0.15 0.60 -0.08 0.37 9.44 9.49 0.05(-0.34)

Page 37 of 41 URL: http://mc.manuscriptcentral.com/tandf/tmph

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

Table 5 .

 5 Predicted pre J and experimental a exp J couplings for some three-membered rings with one heteroatom, see figure4b. For notation see table1.

	X	Y			J	X	J ∆	X e	J	0, X anh	J	0, X har	J ∆	, T X vib	J	300, X vib	J	pre	J	exp	J	exp	-	J	pre
	NH H	J	2,3	6.09 -0.05 -0.15 0.26 -0.01 0.10 6.19
			2 ,3 J ′ ′	4.45 -0.08 -0.08 0.25 0.00 0.17 4.62
			Av	5.27 -0.07 -0.12 0.26 -0.01 0.14 5.41 6.30 0.89( 0.54)
			2 ,3 J ′	2.95 -0.03 0.00 0.27 0.00 0.27 3.22 3.80 0.58( 0.59)
	O H	J	2,3	3.45 -0.02 -0.11 0.22 0.00 0.12 3.57 4.45 0.88( 0.71)
			2,3 J ′	2.12 0.09 0.01 0.24 0.00 0.25 2.37 3.10 0.73( 0.91)
	O Cl	J	2,3	2.15 0.00 -0.09 0.07 0.00 -0.02 2.13 2.68 0.55( 0.36)
			2,3 J ′	0.92 0.06 0.01 0.09 0.01 0.11 1.03 1.42 0.39( 0.49)
	S	H	J	2,3	6.22 -0.05 -0.16 0.41 0.00 0.25 6.47 6.89 0.42( 0.12)
			2,3 J ′	5.08 -0.08 0.03 0.50 0.02 0.54 5.63 5.74 0.11( 0.17)
	PH H	J	2,3	9.55 -0.08 -0.22 0.72 -0.02 0.47 10.02 10.05 0.03(-0.35)
			2 ,3 J ′ ′	8.91 -0.07 -0.16 0.64 -0.02 0.47 9.38 9.26 -0.12(-0.44)
			2,3 J ′	7.31 -0.10 0.06 0.75 0.01 0.82 8.13 7.80 -0.33(-0.20)
	a exp

Acknowledgements

ALE and MPG gratefully acknowledge support from the Dirección General de Investigación of Spain (Grant CTQ2006-15141-C02-01). RHC gratefully acknowledges support from UBACYT (X047) and CONICET (PIP 5119/05).


Hermann Stoll
email: stoll@theochem.uni-stuttgart.de

Can incremental expansions cope with high-order coupled-cluster contributions?

Keywords: ab-initio, coupled-cluster, increments

Local-correlation incremental expansions have been used for the evaluation of high-order coupled-cluster corrections to atomization energies, in exploratory calculations for small molecules. Expansion up to (at least) 3-body terms seems to be necessary in terms of localized orbitals, within the occupied space, while truncation after 2-body inter-atomic terms appears to be sufficient in the virtual space.

An important feature of quantum-chemical wavefunction-based ab-initio methods is their ability to allow for a systematic improvement towards the exact solution of the Schrödinger equation. The price to pay for that is the steep increase of computational effort, both with respect to the size of the one-particle basis set and the level of electron correlation treatment.

For example, the computational effort for a state-of-the-art coupled-cluster calculation with single and double excitations and perturbative account of triples (CCSD(T)) rises as n 3 occ n 4 virt , where n occ and n virt are the number of occupied and virtual orbitals, respectively. Going to a higher-order coupled-cluster treatment, i.e., including triple excitations beyond perturbation theory (CCSDT) or additional quadruples (CCSDTQ) leads to computation times ∼ n 3 occ n 5 virt or ∼ n 4 occ n 6 virt . Most fruitful for overcoming the computational bottleneck and making large molecules accessible to an accurate ab-initio treatment have been the ideas of local correlation and of explicitly correlated wavefunctions. These ideas, which have been pioneered by Pulay and Saebo [1,2,3] and by Kutzelnigg and Klopper [4,5,6], have been in the focus of the scientific work of Hans-Joachim Werner for nearly two decades, see e.g. [7,8,9] and references therein. He succeeded in developing correlation methods up to CCSD(T) which linearly scale with molecular size and allow for inter-electronic terms in the wavefunction thereby greatly reducing the number of basis functions per atom.

The present work is meant as a preliminary exploration on whether similar ideas, local correlation in particular, may be useful for evaluating correlation corrections beyond CCSD(T).

Such corrections are usually quite small, with contributions in the kJ/mol regime, but have proven to be important for accurate work on atomization or excitation energies, see e.g. [10,11,12,13]. The approach used in the present work is the incremental scheme (see, e.g., [14,15]) which consists in a many-body expansion of the correlation energy in terms of local corrections. After giving some computational details in Sect. 2 this approach is applied to small molecules in the following two sections, with an expansion in terms of occupied localized orbitals in Sect. 3 and an expansion in terms of local corrections in the virtual space in Sect. All calculations of the present work have been performed using the Molpro suite of ab-initio programs developed by Werner, Knowles and co-workers [START_REF] Molpro | package of ab initio programs[END_REF][START_REF] Hampel | [END_REF]18,19], and the multireference coupled-cluster (MRCC) code by Kállay [20,[START_REF] Kállay | [END_REF]22] interfaced to Molpro. Basis sets are taken from the series of correlation-consistent polarized valence n-zeta (cc-pVnZ) basis sets of Dunning and co-workers [23], or from the corresponding augmented (aug-cc-pVnZ) sets [24]. Due to the high computational effort for the higher-order coupled-cluster treatment (CCSDT and CCSDTQ), rather small basis sets of double-zeta or triple-zeta quality had to be used for evaluating the non-perturbative triples and quadruples contributions which are at the focus of the present work. (For the sake of comparison, we also give contributions from singles/doubles and perturbative triples using the same basis sets, in Sects. 3 and 4, although for an accurate determination of atomization energies larger basis sets would be needed at the CCSD and CCSD(T) levels.)

For localizing orbitals, the Foster-Boys criterion [25] has been applied. One typical step of the incremental scheme of Sect. 3 consists in correlating a group of localized orbitals (LMOs) while keeping all the other LMOs as frozen. Before correlating such a group in, e.g., CCSD(T) calculations, the active (i.e., non-frozen) orbitals were canonicalized among each other using the multi-configuration self-consistent-field (MCSCF) facility of Molpro, i.e., pseudo-canonical orbitals were generated.

The incremental scheme involves a many-body expansion for the correlation energy of the system, of the form

E corr = E corr 0 + A ǫ A + 1 2! AB ǫ AB + 1 3! ABC ǫ ABC + ... (1)
ǫ A = E corr A ( 2
)
ǫ AB = E corr AB -ǫ A -ǫ B (3)
ǫ ABC = E corr ABC -ǫ AB -ǫ BC -ǫ CA -ǫ A -ǫ B -ǫ C (4)
Here, E corr 0 is a reference correlation energy (to be specified below, see Sects. Finally, the ǫ A , ǫ AB , ǫ ABC are one-body increments and non-additive two-body and threebody corrections, respectively. In the following, we consider differential quantities, e.g. ∆ǫ A , ∆ǫ AB , ∆ǫ ABC , which arise when enhancing the correlation level for the treatment of the system, i.e., restricting the CC calculations to singles and doubles, or considering changes due to perturbative triples, or triples beyond CCSD(T), or quadruples beyond CCSDT.

Incremental expansion in occupied space

Here, the reference level is SCF (E corr 0 = 0), The indices A, B, C denote Foster-Boys localized orbitals, and the correlated orbital space for evaluating E corr A , E corr AB , E corr ABC is restricted to the LMOs indicated, i.e., only 1, 2, or 3 occupied orbitals are correlated simultaneously for a given one-body, two-body, or three-body increment.

As a first example, we consider the CH 4 molecule (with r C-H = 1.102 Å), with four equivalent LMOs describing C-H bonds. Results with cc-pVDZ basis set are listed in Table 1. As known from previous work [14,15], CCSD correlation energies rapidly converge to the full CCSD value, within an incremental expansion. While one-body terms are clearly insufficient, yielding only ∼60% of the full value due to the lack of inter-orbital correlation, two-body terms are already quite accurate, overshooting the full value by ∼3%, and inclusion of three-body terms leads to an accuracy of better than 0.1%. As may be expected, the convergence characteristics somewhat deteriorate when considering the contributions of perturbative triples. Clearly, there are no one-body terms here at all, the two-body terms amount to ∼90% of the full value, and the accuracy reached at the three-body level is ∼2%. 1) and quadruples (Q), i.e., for increasing number of LMOs simultaneously involved, the two-body level of the incremental expansion becomes clearly insufficient, with deviations from the full results by factors of ∼2, but three-body terms lead to an efficient improvement, with deviations of up to 3% only.

For post-CCSD(T) triples contributions (T-(T) in Table

Since the absolute deviations are <0.04 kJ/mol, this seems to provide sufficient accuracy for most purposes. In order to check how sensitive the above findings are with respect to the choice of the orbitals used, let us adopt the worst case, i.e., the case of completely delocalized canonical orbitals (CMO), cf. Table 2. As may be expected, the leading terms in the incremental expansion go down in magnitude by up to factors of 2, but the results reached at the three-body level still seem to be useful: accuracies around 1% are reached for CCSD and (T) correlation energies, and of 7% and 13% (0.07 kJ/mol) for T-(T) and Q contributions, respectively.

Turning now to the discussion of the influence of the one-particle basis sets, let us consider results for CH 4 with triple-zeta basis sets (cc-pVTZ), see Table 3. Compared to cc-pVDZ, the correlation energy naturally increases in magnitude, but in a different way for different excitation classes. CCSD yields the largest absolute enhancement, of course, but percentagewise the increase is largest for the perturbative triples (in CCSD(T)). Post-CCSD(T) triples (in CCSDT) and quadruples (in CCSDTQ) lead to an only moderate increase of 8% (Q) or even a slight decrease (T-(T)). The performance of the many-body expansion turns out to be very similar to that for the smaller basis. At the three-body level, the contribution of double excitations is recovered with an accuracy of <0.1%, that of perturbative triples to <2%, while the accuracy is lower again for post-CCSD(T) triples (4%) and quadruples (6%). Still, the absolute errors in the latter cases are smaller than 0.05 kJ/mol.

A final check deals with the size of the molecule: Table 4 presents results for C 2 H 6 (cc-pVDZ basis, r C-C = 1.544 Å, r C-H = 1.102 Å). It is clear that truncation of the many-body expansion at a given order becomes more stringent with increasing number of subsystems.

Still, the correlation contribution of double excitations is accurate to ∼0.1% at the threebody level, and that of perturbative triples to ∼2%. However, post-CCSD(T) triples are overestimated by ca. 30% here! The total percentage error for the triples is not so very much larger for C 2 H 6 (4%) than for CH 4 (2%), but that of post-CCSD(T) triples is. Four-body terms are needed for reducing the error of the T-(T) contribution. Actually, four-body terms with 3 LMOs sharing a common atom are sufficient to bring the deviation down to <1%, with the largest effect coming from terms connecting the C-C bond with C-H LMOs on both sides.

Further investigations are necessary in order to monitor the performance of the incremental expansion when further increasing the size of the molecules. Interestingly, the convergence of the many-body expansion seems to be more rapid for the full triple contribution, T, than for the perturbative (T) one: the leading three-body increments (involving 3 LMOs with one atom in common) are by more than an order of magnitude smaller than the corresponding two-body ones, for both D and T contributions, while the factor is only ∼5 for the (T) ones.

Anyway, the Q contribution in C 2 H 6 is more similar in accuracy again to that in CH 4 ; at the three-body level, it is by 6% off the full value.

Summarizing, we can say that incremental expansions in terms of localized orbitals form a computationally efficient way of calculating contributions of post-CCSD(T) triples and quadruples in coupled-cluster calculations. Since individual increments rapidly decrease with increasing distance between LMOs, the computational effort eventually increases only linearly with molecular size. However, at least three-body terms have to be included into the expansion in order to get meaningful results.

A decomposition of higher-order coupled-cluster contributions to atomization energies into bond increments has recently been suggested by Bakowies [13]. His suggestion differs from that of the present work in that the increments are taken from larger entities, i.e., the bond increment for C-H is taken from a full calculation for CH

Incremental expansion in virtual space

Here, the indices A, B, C in eqs. 1 -4 denote atoms, and the incremental expansion proceeds by enlarging the virtual space atom by atom. Starting-point is the correlation energy E corr evaluated with a (small) reference basis set, and this energy is improved within E corr A , E corr AB , E corr ABC , by using a (larger) target basis set at one, two, or three atoms (A, AB, or ABC). Localization in occupied space is not necessary here, and the correlated occupied orbital space does not change within the expansion (being that of the full system).

Let us begin our examples with CH 4 again and consider the incremental expansion for extending the basis set from cc-pVDZ to cc-pVTZ. As already mentioned, the perturbative triples contribution, (T), in the coupled-cluster series is that which is mostly affected percentage-wise by the basis-set change (increasing by ca. 70%). It is seen from Table 5 that it converges quite rapidly with the order of the incremental expansion, being by 14% off the cc-pVTZ results at the one-body level, and by only ∼1% at the two-body one. Similar accuracies, of ∼0.5% for doubles (D), of <1% for the T-(T) and Q ones, are achieved at the two-body level for the other contributions. However, this is no stringent test of the incremental expansion for higher-order coupled-cluster methods, since the T-(T) and Q contributions change only very little (by <10%) when going from cc-pVDZ to cc-pVTZ.

As a simple, but more significant test, we next consider the C atom, cf. Table 6. Here, we extend the cc-pVDZ basis set shell by shell towards cc-pVTZ, i.e., the indices A,B,C are taken as shell indices s, p, d, f , and the basis functions of the smaller set are replaced step by step with the functions of the larger basis. For the C atom, the T-(T) and Q contributions increase by 78% and 45%, respectively, when going from cc-pVDZ to cc-pVTZ. Here, the incremental expansion converges even faster than in the preceding example. The T-(T) and Q contributions are off the final (cc-pVTZ) result by <4%, when enlarging the basis set for only one shell at a time, and the changes at the two-body level (enlarging two shells simultaneously) mostly improve the D, (T), and T-(T) contributions to an accuracy <0.5%.

Let us now turn to a more complex and demanding system. The O 3 molecule is among those with the largest high-order coupled-cluster corrections of the molecules considered in Refs. [11,12,13], due to its multi-reference character. Our results to be discussed in the following refer to the experimental geometry (r O-O = 1.272 Å, θ = 116.8 0 ). From Table 7, it is seen that coupled-cluster calculations are sufficiently different for the sp part of the cc-pVDZ basis set and the full cc-pVDZ set, but that the difference can be systematically and efficiently decreased with an incremental expansion. For the D and (T) values, even the one-body level is quite satisfactory, with deviations of ∼1%. For the T-(T) contribution, on the other hand, which even changes sign between cc-pVDZ(sp) and cc-pVDZ, two-body terms are necessary to reach an accuracy of ∼10%. (For the Q contribution, the starting-point is nearer to the final result, and convergence is faster again.) It is interesting to note that inclusion of diffuse functions (aug-cc-pVDZ basis) significantly change the coupled-cluster correlation contributions, by 5% and 13% for D and (T), respectively, and even change sign for T-(T); such augmented basis sets have not been considered in Refs. [11,12,13]. In the same vein, as before, we consider an incremental expansion starting with the sp part of the aug-cc-pVDZ basis and leading eventually to the full aug-cc-pVDZ set, cf. Table 8. As before, the convergence of the expansion is quite fast, with deviations from the full result of 0.1% for (T), and 13% for T-(T), at the two-body level; even one-body results are meaningful here, throughout.

As a final example, let us compare the (not fully satisfactory) incremental expansion for C 2 H 6 of the previous section to a corresponding expansion in virtual (rather than occupied) space, extending the sp part of the cc-pVDZ basis set to the full cc-pVDZ one, see Table 9. It turns out that, although the individual contributions change by up to a factor 2, in the course of the expansion, the results are fully satisfactory at the two-body level, with deviations of 0.5%, 0.7%, 0.8%, and 2.3%, for D, (T), T-(T), and Q, respectively.

Discussion and Conclusion

Incremental expansions of correlation energies in terms of local entities inherently lead to a computational effort for evaluating these energies which scales linearly with the size of the system. This is true provided that correlation effects go to zero sufficiently fast with increasing distance between the entities. It is not strictly true for the case of the expansion in terms of virtual orbitals when starting with a non-minimal basis set, since the reference energy E corr 0 of eq. 1 does not scale linearly; however, in principle, the starting point could just be a basis set consisting of the occupied SCF orbitals (n v = 0). Of course, the systems considered in this paper are too small to enter such a linearly scaling regime. For C 2 H 6 , however, we can at least check the above assumption of the decrease of correlation contributions with distance: the C-H/C-H two-body increment decreases by an order of magnitude when going from nearest-neighbour to second-nearest-neighbour LMO pairs, and this is not only true for the total increment, but also for its D, (T), T-(T), and Q parts separately.

As a disadvantage, the introduction of local entities precludes the use of symmetry which is very helpful in higher-order coupled-cluster calculations. On the other hand, symmetry is not completely lost: it could be partially retained (C s ) for O 3 , e.g., and even with (partial) canonicalization the incremental scheme can still be expected to work (cf. Table 2). It would be perfectly possible, in principle, to merge symmetry-equivalent LMOs or atoms to a single entity for the incremental expansion.

From the viewpoint of computational effort, the expansion for the virtual space seems to be preferable over that for the occupied space, since a) no localization step is needed in the former case, b) truncation at the two-body level seems to be possible, while at least threebody terms are needed in the latter case, and c) the CPU time for standard coupled-cluster calculations rises more steeply with the number of virtual than occupied orbitals (see Sect. 1). However, it would also be possible to somehow combine the two types of expansions as suggested in previous work [15], i.e., to use good-quality basis sets only for those atoms which belong to the domain of the LMOs to be correlated. This would naturally fit into the local-correlation philosophy of Werner and co-workers.

4 .

 4 Some conclusions are drawn in Sect. 5.

  4 , that for C-C is calculated from full calculations for C 2 H 6 and CH 4 . Furthermore, only effective one-body (bond) terms are defined, i.e., two-body and three-body terms are set to zero, which somewhat obscures the importance of the latter found in the present work. Finally, the very different T-(T) contribution for C-H and C-C reported in Ref. [13] (-0.021 kcal/mol and -0.221 kcal/mol, respectively) is in contrast to our observation that one-, two-, and three-body increments involving C-C LMOs change only very little when comparing to the corresponding increments with C-C replaced by C-H LMOs; examples are the T-(T) contribution to the C-C/C-H and C-H/C-H two-body increments (-1.06•10 -4 a.u. vs. -1.18•10 -4 a.u.), or the C-C/C-H/C-H and C-H/C-H/C-H three-body increments (+6.7•10 -5 a.u. vs. +7.3•10 -5 a.u.).

Table 1 :Table 2 :

 12 Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of CH 4 (r C-H =1.102 Å, cc-pVDZ basis, Foster-Boys LMO); incremental expansion in terms of occupied orbitals, in parentheses: contributions to atomization energies 1118863 -0.190762 -0.185314 -0.185472 (-0.108559) Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of CH 4 (r C-H =1.102 Å, cc-pVDZ basis, canonical MO); incremental expansion in terms of occupied orbitals, in parentheses: contributions to atomization energies

Table 7 :Table 9 :

 79 u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of the C atom ( 3 P ground state); incremental expansion in terms of cc-pVDZ → cc-pVTZ basis-set extension, Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of O 3 (r O-O =1.272 Å, θ=116.8 o ); incremental expansion in terms of cc-pVDZ(sp) → cc-437276 -0.610871 -0.610907 -0.610757 (T) -0.026033 -0.032653 -0.033080 -0.033103 T-(T) +0.000346 -0.000302 -0.000130 -0.000119 Q -0.004968 -0.005282 -0.005434 -0.005431 Table 8: Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of O 3 (r O-O =1.272 Å, θ=116.8 o ); incremental expansion in terms of aug-cc-pVDZ(sp) → aug-cc-Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of C 2 H 6 (r C-H =1.102 Å, r C-C =1.544 Å); incremental expansion in terms of cc-pVDZ(sp) → cc-

Table 3

 3

	: Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus pertur-
	bative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of CH 4
	(r C-H =1.102 Å, cc-pVTZ basis, Foster-Boys LMO); incremental expansion in terms of occu-
	pied orbitals, in parentheses: contributions to atomization energies
		1-body	2-body	3-body	full
	D	-0.137088 -0.224347 -0.218738 -0.218880 (-0.127192)
	(T)	0.	-0.005430 -0.006481 -0.006369 (-0.004217)
	T-(T)	0.	-0.000806 -0.000413 -0.000396 (+0.000182)
	Q	0.	-0.000078 -0.000172 -0.000162 (-0.000120)

Table 4 :

 4

	Page 13 of 16		Molecular Physics
		1-body	2-body	3-body	full
	D	-0.204644 -0.352865 -0.340752 -0.340352 (-0.186526)
	(T)	0.	-0.006852 -0.008595 -0.008420 (-0.006466)
	T-(T)	0.	-0.001446 -0.000830 -0.000639 (+0.000008)
	Q	0.	-0.000124 -0.000367 -0.000389 (-0.000330)
				13

Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus perturbative triples, T-(T), and quadruples to the coupled-cluster correlation energies of C 2 H 6 (r C-H =1.102 Å, r C-C =1.544 Å, cc-pVDZ basis, Foster-Boys LMO); incremental expansion in terms of occupied orbitals, in parentheses: contributions to atomization energies URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 5

 5

	: Contributions (a.u.) of doubles, D, perturbative triples, (T), full minus pertur-
	bative triples, T-(T), and quadruples, Q, to the coupled-cluster correlation energies of CH 4
	(r C-H =1.102 Å); incremental expansion in terms of cc-pVDZ → cc-pVTZ basis-set extension
		start	1-body	2-body	3-body	full
	D	-0.185472 -0.229171 -0.217986 -0.218617 -0.218880
	(T)	-0.003772 -0.007261 -0.006285 -0.006342 -0.006369
	T-(T) -0.000400 -0.000412 -0.000396 -0.000394 -0.000396
	Q	-0.000151 -0.000145 -0.000161 -0.000171 -0.000162

Table 6 :

 6 Contributions (a.


MD Anette-Marie Machata
email: anette-marie.machata@meduniwien.ac.at

MD Barbara Kabon

MD Harald Willschke

Daniela Prayer Md

MD Peter Marhofer

Upper airway size and configuration during propofol-based sedation for magnetic resonance imaging: an analysis of 138 infants and children

Background

Propofol is widely used for pediatric sedation. However, increasing depth of propofol sedation is associated with airway narrowing and obstruction. The aim of this study was to objectively assess airway patency during a low-dose propofol-based sedation regimen by measuring upper airway size and configuration with MRI in spontaneously breathing infants and children.

Methods

Magnetic resonance images of the upper airway were obtained in 138 infants and children, aged up to 6 years. Cross-sectional area, anteroposterior dimension, and transverse dimension were measured at the level of the soft palate, the base of the tongue, and the tip of the epiglottis. Sedation was induced with i.v. midazolam 0.1 mg⋅kg -1 , nalbuphine 0.1 mg⋅kg -1 and propofol 1 mg⋅kg -1 , and maintained with propofol 5 mg⋅kg -1 ⋅h -1 .

Results

Median (IQR) age was 36 [START_REF] Malviya | Depth of sedation in children undergoing computed tomography: validity and reliability of the University of Michigan Sedation Scale (UMSS)[END_REF]48) months, and mean body weight was 13,7 ± 5,6 kg. Airway patency was maintained in all infants and children. The narrowest part of the pharyngeal airway was measured at the level of the base of the tongue. Anteroposterior dimensions were narrower than transverse dimensions in all age groups at all measurement sites. Transverse dimensions increased with age at all measurement sites, while anteroposterior dimensions did not increase comparably. No patient demonstrated respiratory or cardiovascular adverse events. All MRI were completed successfully without sedation failure.

Conclusion

Airway patency was maintained in all infants and children sedated with this low-dose propofol-based sedation regimen. Keywords: Anesthesia, pediatric; Anesthetic techniques, i.v.; Anesthetics i.v., propofol; Airway, anatomy; Brain, magnetic resonance imaging; Infants and children require deep sedation for successful management of magnetic resonance imaging (MRI) examinations. The noisy environment in the MRI suite and the duration of the examinations preclude a lower sedation state in these patients. Propofol is widely used for pediatric sedation, but increasing depth of propofol anesthesia is associated with the risk of upper airway obstruction [START_REF] Mallory | Propofol vs pentobarbital for sedation of children undergoing magnetic resonance imaging: results from the Pediatric Sedation Research Consortium[END_REF][START_REF] Eastwood | Collapsibility of the upper airway at different concentrations of propofol anesthesia[END_REF]. Because maintenance of a patent airway during MRI procedures under spontaneous respiration is of utmost importance, airway management is still controversial in this setting. Sedation in spontaneously breathing children is advocated by several authors [START_REF] Mahmoud | Effect of increasing depth of dexmedetomidine anesthesia on upper airway morphology in children[END_REF][START_REF] Usher | Propofol total intravenous anesthesia for MRI in children[END_REF][START_REF] Levati | Propofol anaesthesia in spontaneously breathing paediatric patients during magnetic resonance imaging[END_REF], others prefer laryngeal mask airway insertion [START_REF] Kol | Open-label, prospective, randomized comparison of propofol and sevoflurane for laryngeal mask anesthesia for magnetic resonance imaging in pediatric patients[END_REF] or tracheal intubation to secure the airway [START_REF] Bordes | Which anaesthesia for children undergoing MRI? An Internet survey in the French university hospitals[END_REF][START_REF] Young | Anaesthesia for children and infants undergoing magnetic resonance imaging: a prospective study[END_REF]. An appropriate low propofol dosage which maintains spontaneous respiration and ensures sufficient sedation for successful MRI completion would probably minimize the risk of airway obstruction.

In a previous study we described the safety and efficacy of a propofol-based sedation regimen for MRI examinations in 500 infants and children (9). Airway patency was only deduced from normal peripheral oxygen saturation and normal end-tidal carbon dioxide values in these patients. However, MRI allows direct determination of upper airway size and configuration, even if these measurements were only performed in a small number of pediatric patients until now [START_REF] Crawford | Effect of propofol anesthesia and continuous positive airway pressure on upper airway size and configuration in infants[END_REF][START_REF] Evans | Effect of increasing depth of propofol anesthesia on upper airway configuration in children[END_REF][START_REF] Litman | Morphologic changes in the upper airway of children during awakening from propofol administration[END_REF].

Therefore, the aim of this study was to objectively assess airway patency by measuring upper airway size and configuration with MRI in spontaneously breathing infants and children sedated with the aforementioned propofol-based regimen.

Materials and Methods

Patients

With ethics hospital committee approval, 150 consecutive ASA I and II infants and children, aged up to 6 years, who required sedation for elective cranial MRI examination, were enrolled in this study. Informed consent was obtained from the parents of all infants and children.

Exclusion criteria were ASA status ≥ III, severe pulmonary or cardiovascular disease, anatomic airway abnormalities or extreme tonsillar hypertrophy, gastroesophageal reflux, a body weight of at least 30% more than ideal, a history of propofol intolerance, or known fat metabolism disorder.

Procedure

On the day of the procedure, all patients were admitted with their parents at the pediatric day care ward. An intravenous access was obtained, and then the patient was transferred to the MRI suite accompanied by his/her parents. In a preprocedure interview, all parents were instructed about patient fasting as recommended by the "American Society of Anesthesiologists Preprocedure Fasting guidelines" [START_REF] Cote | Guidelines for monitoring and management of pediatric patients during and after sedation for diagnostic and therapeutic procedures: an update[END_REF][START_REF]Practice guidelines for sedation and analgesia by non-anesthesiologists[END_REF].

After the patients received midazolam 0.1 mg⋅kg -1 intravenously, they were moved into the MRI suite.

Then, nalbuphine 0.1 mg⋅kg -1 , followed by a loading dose of propofol 1 mg⋅kg -1 was administered. Supplemental doses of propofol 0.5 mg⋅kg -1 were administered until adequate sedation was achieved. Sedation was considered adequate, when the patient slept, arousable only with significant physical stimulation (University of Michigan Sedation Scale (UMSS3)) [START_REF] Malviya | Depth of sedation in children undergoing computed tomography: validity and reliability of the University of Michigan Sedation Scale (UMSS)[END_REF]. The maintenance of spontaneous respiration was verified. The head position was standardized in a neutral position such that the angle between the horizontal plane of the MRI table and a line connecting the tragus of the ear and the lateral corner of the eye was 110° [START_REF] Crawford | Effect of propofol anesthesia and continuous positive airway pressure on upper airway size and configuration in infants[END_REF][START_REF] Evans | Effect of increasing depth of propofol anesthesia on upper airway configuration in children[END_REF]. No devices (nasal or oral airway) that assist with patency of the airway were used, and no attempt was made to open or close the mouth [START_REF] Litman | Morphologic changes in the upper airway of children during awakening from propofol administration[END_REF]. Furthermore, no attempt was made to obtain images at a certain timepoint of the respiratory cycle.

Supplemental oxygen was provided by pediatric face mask with a gas flow rate of 2 l⋅min -1 , and sedation was maintained with propofol 5 mg⋅kg -1 ⋅h -1 using syring pumps suitable for MRI.

Special monitors, not interfering with the MRI procedure, were placed in the examination room and in the observation room to monitor vital signs. Heart rate, arterial oxygen saturation (SpO 2 ) and end-tidal carbon dioxide (EtCO 2 ), sampled via special nasal cannula, were monitored continuously during the procedure and recorded at 5 min intervals by the anesthesiologist. Non-invasive arterial pressure was determined immediately before induction of sedation and at the end of the examination, but not during the MRI procedure, as stimulation from the blood pressure cuff would arouse an appropriately sedated child [START_REF]Practice guidelines for sedation and analgesia by non-anesthesiologists[END_REF].

Intervention was considered necessary when a decrease in SpO 2 to < 94%, an increase in EtCO 2 > 50 mmHg, apnea (cessation of spontaneous respiration for 20 seconds), bradycardia (> 20% decrease in heart rate from baseline), or arrhythmia occurred.

For safety reasons, a camera was placed in the MRI suite to monitor patients, and the anesthesiologist closely observed the patients while sedated and provided interventions as needed: In case of inadequate sedation, such as spontaneous movements and agitation during the imaging process, an additional bolus of propofol 0.5 mg⋅kg -1 was administered, and the propofol infusion rate was increased by 1 mg⋅kg -1 ⋅h -1 up to a total of 8 mg⋅kg -1 ⋅h -1 . In case of coughing or suspected airway obstruction, the MRI examination was interrupted, the patient was taken out of the MRI unit, and the airway patency was assessed. When partial airway F o r P e e r R e v i e w 6 occlusion was diagnosed, the neck was again extended slightly and the chin was supported with a tape. When total airway occlusion was noted, bag-valve-mask ventilation was initiated, and a laryngeal mask airway was inserted, if necessary. Tracheal intubation was performed, when these maneuvers did not relieve obstruction.

After the MRI examination was completed, the propofol infusion was terminated, and the patients were transferred from the MRI suite to the induction room.

Measurements

Magnetic resonance images were obtained using a 1.5T Philips Achieva Release 1.2.2 scanner (Philips, Medizinische Systeme GmbH, Austria) with a master/NOVA gradient strength and a circular polarized birdcage headcoil. A M2D FFE sequence slice thickness 10 mm/10 mm gap was used as a scout view in a three-plane gradient echo. Then a T1 3D FFE axial/oblique sequence was used. This sequence allows unequivocal recognition of the landmarks that were used for measurements. After the examinations, DICOM images were stored on computer and analyzed by a single investigator blinded to patient data. After image magnification (x 3), anteroposterior dimension, transverse dimension, and cross-sectional area (CSA) were determined at the level of the soft palate, the base of the tongue, and the tip of the epiglottis using image analysis software (Osirix Imaging Software, version 3.6.1) (figures 1-3). Because image acquisition at different phases of the respiratory circle was impossible in this setting, measurements were obtained from two subsequent slices at each anatomic level, and the average of both was calculated.

Data analysis

Data were descriptively analyzed and presented as median (IQR) or mean (SD) depending on their distribution. Normal distribution was assessed with q-q plot and Shapiro-Wilk test.

CSA measurements among the three age groups were compared with Kruskal-Wallis test.

Friedman test was used to compare paired nonparametric data (differences between pharyngeal dimensions at different levels). Post hoc pairwise comparisons were performed

Results

One hundred fifty patients were initially enrolled in the study. Twelve were excluded because of various reasons that included abnormal laryngeal anatomy (4), technical poor images (2), and data retrieval problems [START_REF] Kol | Open-label, prospective, randomized comparison of propofol and sevoflurane for laryngeal mask anesthesia for magnetic resonance imaging in pediatric patients[END_REF]. From the remaining 138 patients, 60 were boys and 78 were girls. Median (IQR) age was 36 [START_REF] Malviya | Depth of sedation in children undergoing computed tomography: validity and reliability of the University of Michigan Sedation Scale (UMSS)[END_REF]48) Airway patency was maintained in all patients at the 3 measurement sites (table 1).

The narrowest part of the pharyngeal airway was measured at the level of the base of the tongue in all patients. Anteroposterior dimensions, transverse dimensions, and also CSA were significantly smaller at the level of the base of the tongue compared to the soft palate and epiglottis level, which were similar. Anteroposterior dimensions were narrower than transverse dimensions in all age groups at all measurement sites (table 1).

Transverse dimensions increased with age at all measurement sites, while anteroposterior dimensions did not increase comparably. There was a significant linear increase of CSA with age at all levels. However, the interindividual variation was wide, and the correlation was only moderate (figure 4).

During either phase of the study, no patient demonstrated clinically significant upper airway obstruction as manifested by decreased oxygen saturation or increased capnometric values. None of our patients suffered from cardiovascular events such as bradycardia or arrhythmia during and after sedation. Paradoxical reaction to sedation or agitation was not seen in our patients. All scheduled MRI examinations could completed without sedation failure.

Discussion

To our knowledge, this prospective, observational study documents the largest case series of MRI measurements of the pediatric airway in the recent literature. We found that all 138 infants and children had a patent airway during sedation with this low-dose propofol-based regimen. The narrowest part of the pharyngeal airway was measured at the level of the base of the tongue. Anteroposterior dimensions were narrower than transverse dimensions in all age groups at all measurement sites. An additional finding was that transverse dimensions increased with age at all measurement sites, while anteroposterior dimensions did not increase comparably. All dimensions at all measurement sites were at least as wide than previously described during propofol sedation with comparable dosages [START_REF] Crawford | Effect of propofol anesthesia and continuous positive airway pressure on upper airway size and configuration in infants[END_REF][START_REF] Evans | Effect of increasing depth of propofol anesthesia on upper airway configuration in children[END_REF]. All MRI examinations could be completed successfully.

The narrowest part of the upper airway was measured at the base of the tongue. This finding concurs with the traditional view of several studies in adults, that airway obstruction may be caused predominantly by apposition of the base of the tongue to the posterior pharyngeal wall [START_REF] Fink | Roentgenographic studies of the oropharyngeal airway[END_REF][START_REF] Morikawa | Influence of the headjaw position upon upper airway patency[END_REF], and challenges the results of more recent studies in children [START_REF] Crawford | Effect of propofol anesthesia and continuous positive airway pressure on upper airway size and configuration in infants[END_REF][START_REF] Evans | Effect of increasing depth of propofol anesthesia on upper airway configuration in children[END_REF]. These studies concluded that anesthetic-induced attenuation of upper airway muscle activity results in an airway caliber decrease most pronounced at the level of the soft palate and the epiglottis. However, in these studies sample size was small and propofol dosage up to three times as high as the dosages administered to our patients. Furthermore, these studies were designed to determine the effects of increasing depth of propofol anesthesia on airway size, and not to evaluate an appropriate low propofol dosage for pediatric MRI sedation. We excluded a posterior displacement of the tongue due to head position as the reason for our results, because all mentioned studies investigated their patients with the head in the same neutral position. It is known that chin lift and mouth closure widens the entire upper airway mostly between the tip of the epiglottis and the posterior pharyngeal wall [START_REF] Reber | Effect of combined mouth closure and chin lift on upper airway dimensions during routine magnetic resonance[END_REF]. However, we made no attempt to optimize mouth closure and chin lift, but we confirm the observation of a previous The finding that anteroposterior dimensions were narrower than transverse dimensions reflects previously described effects of propofol on pharyngeal airway morphology. Some authors suggest that propofol is associated with a differential effect between separate pharyngeal dilatator muscles, which is represented by changes of the upper airway shape during sedation and awakening from sedation (19,20). Under their sedation regimen, the upper airway was oblong shaped, with the anteroposterior dimension larger than the transverse dimension. On awakening, the transverse dimension was larger than the anteroposterior dimension, whereas the CSA was unchanged [START_REF] Litman | Morphologic changes in the upper airway of children during awakening from propofol administration[END_REF]. The authors administered 200 mcg⋅kg -1 ⋅h -1 propofol, which is more than twice the dosage our patients received.

Therefore we suggest that airway shape in our patients better represents the morphology during the awake state, which emphasizes the safety of our sedation regimen.

The finding that transverse dimensions increased with age while anteroposterior dimensions did not increase comparably, has not been described previously. Data of pharyngeal airway measurements in a comparable large collective are missing, and only one study determined developmental changes of the laryngeal dimensions in sedated children (19). However, measurements were performed more caudad between the vocal cords and the cricoid ring, and the authors concluded that anteroposterior and transverse dimensions increased linearly with age. If we assume that the more cephalad airway may have the same relation between anteroposterior and transverse dimension, we attribute the minor anteroposterior increase of 11 our patients to an additional effect of propofol sedation on the pharyngeal airway. Increasing depth of propofol anesthesia seems to decrease airway caliber predominantly to anteroposterior narrowing [START_REF] Crawford | Effect of propofol anesthesia and continuous positive airway pressure on upper airway size and configuration in infants[END_REF].

One may argue that adding the partial agonist-antagonist opioid nalbuphine to a sedation regimen is not necessary for the inherently painless MRI examination. According to recent guidelines, sedatives or hypnotics alone have to be preferred for non-painful procedures [START_REF] Cote | Guidelines for monitoring and management of pediatric patients during and after sedation for diagnostic and therapeutic procedures: an update[END_REF].

However, nalbuphine was included in our protocol, since nalbuphine 0.1 mg⋅kg -1 , administered before anesthesia induction, is described to decrease the frequency of spontaneous, unintentional movements induced by propofol (21). Furthermore, previous publications demonstrate that nalbuphine, administered at the end of MRI examinations, prevents emergence agitation (22). Due to the assumed additional sedative effect of nalbuphine to propofol, we were able to complete all MRI examinations without interruption with one of the lowest propofol doses (9). Pure propofol regimen seem to require doses up to 200 mcg⋅kg -1 ⋅h -1 , and even then the rate of MRI pausing due to unintentional movements is 10-20 % [START_REF] Usher | Propofol total intravenous anesthesia for MRI in children[END_REF]23).

In summary, we demonstrated with objective measurements of the airway morphology via MRI in a large collective of infants and children with normal upper airway anatomy and neither obesity nor extreme tonsillar hypertrophy that our low-dose propofol-based sedation is safe regarding airway patency.

  All children wear earphones to diminish the stimulation caused by the noise during MRI acquisition phases.

7 with

 7 Wilcoxon rank sum test. To compensate for multiple comparisons the Bonferroni correction was used. Pearson´s r was used to determine the relationship between CSA and age. A P < 0.05 was considered statistically significant. SPSS software was used for all analyses.

  ), that the majority of children under propofol sedation maintain their jaw in the neutral position between opened and closed. Additionally, the inspiratory narrowing of the pharyngeal airway could have influenced our results. The only attempt to minimize this limitation was made by averaging the values of two subsequent slices at each measurement site. Image acquisition was not triggered manually to early expiration, because the relatively high respiratory rate of the investigated infants prevented us from acquiring images at a fixed point of respiratory cycle.

Figure 1 A

 1 Figure 1

Figure 2 A

 2 Figure 2

Figure 3 A

 3 Figure 3

Figure 4

 4 Figure 4

  patients sedated with propofol. Anesthesiology 1999; 90: 1617-1623. 19 Litman RS, Weissend EE, Shibata D, Westesson PL. Developmental changes of laryngeal dimensions in unparalyzed, sedated children. Anesthesiology 2003; 98: 41-45. 20 Isono S. Developmental changes of pharyngeal airway patency: implications for pediatric anesthesia. Paediatr Anaesth 2006; 16: 109-122. 21 Borgeat A, Fuchs T, Wilder-Smith O, Tassonyi E. The effect of nalbuphine and droperidol on spontaneous movements during induction of anesthesia with propofol in children. J Clin Anesth 1993; 5: 12-15. 22 Dalens BJ, Pinard AM, Letourneau DR, Albert NT, Truchon RJ. Prevention of emergence agitation after sevoflurane anesthesia for pediatric cerebral magnetic resonance imaging by small doses of ketamine or nalbuphine administered just before discontinuing anesthesia. Anesth Analg 2006; 102: 1056-1061. 23 Bryan YF, Hoke LK, Taghon TA, et al. A randomized trial comparing sevoflurane and propofol in children undergoing MRI scans. Paediatr Anaesth 2009; 19: 672-681.

3 ,

 3 Abbreviations: CSA sp -cross-sectional area at the level of the soft palate, AP sp -

table 1 .

 1 months; mean body weight was 13,7 ± 5,6 kg, median height was 95 (78, 104) cm. Patients were divided in three age groups: infants < 12 months, children between 12 and < 36 months, and > 36 months. Pertinent data are shown in

Table 1 :

 1 Demographic and morphometric dataData are presented as means ± SD, or median (IQR) Comparisons between the age groups were performed with Kruskal-Wallis test. Post hoc comparisons were performed using Bonferroni corrections.

	Age (month)	< 12	12-< 36	36-< 72	P
	Number of patients	27	40	71
	Age (month)	6 (3, 10)	21 (17, 27)	48 (42, 58)
	Weight (kg)	6.9 ± 3	11.1 ± 2	17.7 ± 4.1	< 0.001
	Height (cm)	64 (54, 73)	86 ( 80, 91)	104 (100, 108) < 0.001
	CSA sp (mm 2 )	103 (84, 127) *#	135 (103, 184)	167 (119, 196) < 0.001
	AP sp (mm)	8.2 (7.6, 9.4) #	9.2 (7.8, 11.1)	10.1(8.7, 11.5)	0.002
	T sp (mm)	13 (11.6, 14.5) *#	16.5 (13.8, 18.5) 17.7 (15.6, 19.4) < 0.001
	CSA t (mm 2 )	79 (66, 101) *#	120 (101, 153)	129 (94, 165)	< 0.001
	AP t (mm)	7.3 (6.1, 8.2) *#	9.1 (7.7, 10.4)	9 (8, 10.8)	< 0.001
	T t (mm)	12.4 (10.6, 14) *#	15.1 (12.8, 17.5) 15.8 (13.1, 18.9) < 0.001
	CSA e (mm 2 )	89 (78, 127) *#	137 (113.1, 170)	143 (111, 175) < 0.001
	AP e (mm)	8.6 (6.5, 9.9) *	9.9 (8.6, 11.5)	9.3 (7.9, 11.3)	0.048
	T e (mm)	13.2 (11.5, 14.5) *# 16.5 (14.9, 18.2) & 18.2 (16.2, 20) < 0.001
	* < 12 significantly different to 12-< 36
	# < 12 significantly different to 36-< 72
	& 12-< 36 significantly different to 36-< 72


K.-A Mardal

B F Nielsen

Efficient preconditioning of optimality systems

 we explain how the mapping properties of the operators appearing in the OS can be employed to define efficient preconditioners for finite element (FE) approximations of the involved saddle point problem also for problems that are ill-posed. More specifically, it turns out that it is possible to define a scheme such that the number of iterations needed to solve the preconditioned problem is bounded independently of the mesh parameter h, used in the FE discretization, and only increases moderately as the regularization parameter α tends towards zero. In fact, if the associated energy norm is used to define the stopping criterion for the iteration process, then the number of iterations required (in the severely illposed case) cannot grow faster than O((ln(α)) 2 ). This result is obtained by carefully analyzing the properties of the involved operators and thereby revealing the distribution of the eigenvalues of the preconditioned OS. One advantage of our approach is that it only rely on standard elliptic preconditioners based on e.g. multigrid and domain decomposition.

Our theoretical results will be illuminated by a number of numerical experiments.


Cornelis Vuik

Rohit Gupta

Kees Lemmens

A fast GPU Implementation of the Deflated Preconditioned Conjugate Gradient method

A fast GPU Implementation of the Deflated Preconditioned Conjugate Gradient method

List of authors:

C. Vuik 1 R. Gupta 2
C.W.J. Lemmens 3

Large sparse linear systems are mostly solved by Preconditioned Krylov Methods combined with some form of coarse grid acceleration. We consider the Deflated Preconditioned Conjugate Gradient Method. Most of the important building blocks of this algorithm could be optimized for execution on parallel architectures. Recently, with the advent of General Purpose Computing on the Graphical Processing Unit (GPU) it is possible to achieve 10 -100 times reduction in computing times. However, many preconditioners involved in the solution of the systems of interest do not run optimally on the GPU. A large class of preconditioners is based on incomplete Cholesky decompositions. Solving the lower and upper triangulars system is inherently sequential. So the challenge is how to make these preconditioners suitable for the GPU which can only be used with parallel methods. In this paper we show that it is possible to achieve a 20 fold speedup for a combination of suitable building blocks.

We use two levels of Preconditioning with the CG method to solve a linear system [START_REF] Tang | Efficient Deflation Methods applied to 3-D Bubbly Flow Problems Electronic[END_REF]. Our system arises from the discretization of the Pressure Correction Equation. This equation is the most time-consuming step in the solution of the Incompressible Navier-Stokes Equation using the Level Set Method. This method, as suggested in [START_REF] Van Der Pijl | A mass conserving Level-Set Method for modelling of multi-phase flows[END_REF], is of interest in modeling bubbly flows. The Partial Differential Equations describing the Pressure Correction have been discretized through the use of finite differences. The linear system is of the form

Ax = b, A ∈ R n×n , (1)
where n is the number of degrees of freedom. We assume that A is symmetric positive definite.

With the advent of the Component Unified Device Architecture (CUDA) paradigm of computing available on NVIDIA GPU devices, it has become easier to write numerical software. Characteristic features of GPU devices are: one should use vectorizable algorithms (fine scale parallelism), reuse the fast shared memory as much as possible, if statements etc. lead to severe degradation of the performance, minimize communication between CPU and GPU, and the obtained accuracy and error checking is less on the GPU as compared to the CPU.

We perform our experiments using 5 point Laplacian matrix resulting from a 2D square grid.

The same experiments have been done for a matrix originating from a two phase flow. The

Scientific And Industrial Applications, Preconditioning 2011, May 2011, Bordeaux, France. inria-

Delft University of Technology, Faculty of EEMCS, Mekelweg, Delft, The Netherlands

Delft University of Technology, Faculty of EEMCS, Mekelweg, Delft, The Netherlands 3 Delft University of Technology, Faculty of EEMCS, Mekelweg, Delft, The Netherlands 1 methods below have been implemented and tested in the given order. We start with the standard CG method, add a standard preconditioner (BIC) and combine it with a second level preconditioner (deflation). After that we replace the standard preconditioner by a parallel preconditioner and optimize the code for the GPU which leads to the fastest method: DIPCG2. We use a Q9650 Intel Quad Core CPU however we only utilize a single core. We optimize it to use SSE instructions, unrolling loops and vectorizing using compiler switches. We also use the Meschach Blas Library for the Blas routines on the CPU. The GPU we use is a Tesla C1060 from NVIDIA. We use CUDA for writing our code on the GPU. The CU BLAS and M AGM A libraries are used when Blas functions are needed in the GPU version. Table 1 shows the results for a 512 × 512 grid. Note that the combination of the IP preconditioner and the deflation approach leads to a low number of iterations. Furthermore, after optimization, the gain in computing time is around a factor 20. For further details we refer to [START_REF] Gupta | GPU Implementation of Deflated Preconditioned Conjugate Gradient submitted[END_REF].