text:
# The Role of the Basis Set: Assessing Density Functional Theory

## I Introduction

In the past years, Density Functional Theory (DFT) has become a very important approach for computational quantum chemistry. The Kohn-Sham implementation of DFT critically depends on the quality of the exchange-correlation functional for its success. Recently, various second-generation functionals (such as PBEPBE ,mPW91mPW91 ,VSXCVSXC , PBE0PBE0 ,PKZBPKZB ,HCTH/93 and B97-1HCTH93 ,HCTH/120HCTH120 , HCTH/407HCTH407 , OPTXOPTX ,B972B972 ,$`\tau `$-HCTH and its hybridtHCTH ,mPW1KmPW1K ,B97B97 , and B98B98 ) have been developed. These add to the numerous functionals that are already available and commonly used in standard program packages (such as B88XB88X , B3P91 and B3LYPB3P91 ,VWNVWN ,P86P86 ,LYPLYP ,P91XP91X , and P91cP91c ). All of the functionals were developed from a wide variety of considerations, with most of them focused on the exchange-correlation hole and employing different philosophies in their approximations.

Usually, after deciding upon the mathematical form of the functional, its parameters have to be obtained. The latter is a difficult process where different routes can be followed. Some of these functionals only use parameters that were determined by considering known boundary conditions that the functional or density should obey, and for properties of certain idealised systems like the uniform electron gas PBE ; mPW91 ; VWN ; P86 ; LYP ; P91X ; P91c ; hence their parameters are mainly determined by physical properties. The majority of commonly used functionals have been determined by fitting at least one of their parameters to atomic or molecular data. A number of the latter functionals have been generated using numerical DFT (i.e. basis set-free) methods B97 ; B98 ; B88X ; B3P91  in their fits (although they in principle introduce numerical noise), but most of the functionals have been fit to molecular sets with limited basis sets, usually of triple-zetaVSXC ; PBE0 ; PKZB ; HCTH93 ; HCTH120 ; HCTH407 ; OPTX ; B972 ; tHCTH  or even double-zetamPW1K  quality.

Many DFT users are thus overwhelmed by the sheer number of functionals and possibilities that can be used, whereas with ab initio methods, the choices are clearly defined and mainly determined by a trade-off between rigor and computational cost. Very often because of sheer user inertia, first-generation functionals are applied rather than the more accurate second-generation functionals. Generally, a basis set of double- or triple-zeta quality is applied. Meanwhile, systematic studies on the dependency of the basis set and functionals remain sparse basquad ; basquad2 ; basdip ; basat1 ; basblyp ; basfreq ; basrev . In addition, it is not a priori clear that basis sets optimised for wavefunction ab initio methods are the most optimal choice of DFT. Recently, new basis sets especially optimised for DFT Jensen1 ; Jensen2 ; GSAW ; DFO  have been proposed. However, we have to keep in mind that in these cases the basis sets have been optimised for specific functionals, like BLYPJensen1 ; Jensen2  or the Local Spin Density Approximation (LSDA) GSAW ; DFO , further adding to the confusion. In sum, we have a plethora of functionals developed for specific basis sets and additionally, a variety of basis sets developed for specific functionals. In the former case, the question arises which basis sets can and should be used for developing new functionals, and which basis set should then be employed when applying them. This is a major consideration in functional development, since the question remains: with which basis set the error of the basis set becomes comparable to that of the functional? Even if part of the basis set error is absorbed into the parameters of the functional, the question remains of how transferable such a functional will be to other basis sets- be they further from, or nearer to the infinite basis set limit.

In this contribution, we will evaluate and compare various functionals using basis sets of double-zeta sp to quadruple-zeta spdfg quality. In addition, we will fit a functional, in this case the Hamprecht-Cohen-Tozer-Handy (HCTH) functional HCTH93 , to slightly modified G2-1G2-1  and G3 sets of moleculesG3 . These sets are denoted as the 147 HCTH120  and 407 Sets www . All these fits will be carried out for several basis sets of double-zeta, triple- and quadruple-zeta quality. The same fitting procedure was employed for the functional HCTH93 ; HCTH120 ; HCTH407 ; HCTH407+  using the TZ2P basis set.

HCTH is a post-local spin density approximation (post- LSDA) functional, meaning that it factorises the LSDA functional forms ($`F_{LSDA}`$), which can be found elsewhere PW91 :

$`E_{xc}={\displaystyle \underset{\gamma =x,c_{\sigma \sigma },c_{\alpha \beta }}{}}E_\gamma ={\displaystyle \underset{\gamma }{}}{\displaystyle \underset{q=0}{\overset{m}{}}}c_{q,\gamma }{\displaystyle F_{LSDA,\gamma }(\rho _\alpha ,\rho _\beta )f_{\gamma ,q}(\rho _\alpha ,\rho _\beta ,x_\alpha ^2,x_\beta ^2)𝑑𝐫}`$ (1)
where, $`f_{\gamma ,q}`$ denotes the perturbation from the uniform electron gas if $`c_{0,\gamma }=1`$.

$`f_{\gamma ,q}`$ $`=`$ $`u_\gamma ^q=\left({\displaystyle \frac{\theta _{\gamma \sigma }x_\sigma ^2}{1+\theta _{X\sigma }x_\sigma ^2}}\right)^q`$ (2)
$`x_\gamma `$ is closely related to the reduced density gradient, and $`\theta `$ are fixed coefficients, which have been fit to atomic data B97 .

$`x_\sigma ^2={\displaystyle \frac{(\rho _\sigma )^2}{\rho _\sigma ^{8/3}}}`$ (3)
When employing the form in Eqn. 1 up to fourth order in $`m`$, we obtain 15 linear coefficients (because of exchange, like-spin and opposite-spin correlation), which are easily parameterised by minimising $`\mathrm{\Omega }`$:

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \underset{m}{\overset{n_E}{}}}w_m(E_m^{exact}E_m^{KS})^2+{\displaystyle \underset{l,X}{\overset{n_G}{}}}w_{l,G}\left({\displaystyle \frac{E_l^{KS}}{X}}\right)^2`$ (4)
$`+{\displaystyle \underset{j,\sigma }{\overset{n_v}{}}}w_{j,v}{\displaystyle (v_{j,\sigma }^{ZMP}+k_{j,\sigma }v_{j,\sigma }^{KS})^2\rho _{j,\sigma }^{2/3}𝑑𝐫}`$
The three summations correspond to errors of the energies, gradients and exchange-correlation potentials of each molecule, respectively. In all cases, K-S denotes the calculated property; hence we have the energy difference between the exact and calculated energy in the first sum. In the second sum the exact gradients (at equilibrium geometry) should be zero. In the final term, we fit to the exchange-correlation potentials determined by the Zhao-Morrison-Parr method ZMP  from high-level ab initio densities, which are shifted by a constant $`k`$ because of the effects of the quantum-mechanical integer discontinuity. All these contributions need to be weighted by appropriate weights $`w`$, which have been determined and reported in previous papers HCTH407 . The weights $`w`$ consist of several factorised weights making contributions for each molecule in order to ensure a balanced functional.

In the next section, we will refit the HCTH generalised gradient approximation (GGA) functional to numerous basis sets. In addition, we will assess the performance of these functionals with basis sets other than those used for the parameterisation. In the third section we will apply the same procedure to hybrid functionals, focusing on the amount of exact exchange needed for those functionals depending on the basis set. In the last section, we will assess several functionals with two triple-zeta basis sets (namely, TZ2PTZ2P  and cc-pVTZPVXZfirst ; PVXZsecond ).

## II GGA Functionals and Basis Sets

For all calculations, we used the Cadpac suite of programs Cadpac , using a standard ’high’ grid for the density functional calculations. In the basis set evaluation, we used the 3-21G321G1 ; 321G2 , 6-31G631G1 ; 631G2 , 6-311G6311G1 ; 6311G2  (the latter two with various combinations of diffuse and polarisation functions), DZP, TZ2PTZ2P , DFODFO , cc-pVDZ, aug-cc-pVDZ, cc-pVTZ, aug-cc-pVTZ and cc-pVQZ basis setsPVXZfirst ; PVXZsecond . Some of the basis sets specifically constructed for density functionals discussed in the introduction are not available (yet) for second-row atoms Jensen1 ; Jensen2 ; GSAW  and thus could not be used.

In Table 1, we compare the performance of the functionals with the basis set with which they were fit. HCTH/147@3-21G, for example, denotes an HCTH functional fitted to the 147 set with a 3-21G basis set. The first column contains the RMS energy error of the functional for the atomisation and dissociation energies, electron affinities, proton affinities, and ionisation potentials in the 147 Set. The second column displays the sum of all the gradients calculated at the equilibrium geometry in atomic units. The gradients have been shown to correlate with the accuracy of the bond distances and angles obtainedHCTH407 ; tHCTH . The third column is the square of all the errors in the exchange-correlation potential contributions. All three errors are actually fitted according to equation 4. In the fourth column, we evaluate $`Q`$, which is the sum of all three weighted errors with uniform weights rather than $`\mathrm{\Omega }`$:

$`Q`$ $`=`$ $`750\times {\displaystyle \underset{m}{\overset{n_E}{}}}(E_m^{exact}E_m^{KS})^2+500\times {\displaystyle \underset{l,X}{\overset{n_G}{}}}\left({\displaystyle \frac{E_l^{KS}}{X}}\right)^2`$ (5)
$`+{\displaystyle \underset{j,\sigma }{\overset{n_v}{}}}{\displaystyle (v_{j,\sigma }^{ZMP}+k_{j,\sigma }v_{j,\sigma }^{KS})^2\rho _{j,\sigma }^{2/3}𝑑𝐫}`$
Since all three contributions are important and the GGA functionals have been fit to a value similar to $`Q`$, it is probably the most important diagnostic. Nevertheless, as only the first two sums are directly apparent from the energy calculations, we also evaluate Q1, which excludes the potential:

$`Q_1`$ $`=`$ $`750\times {\displaystyle \underset{m}{\overset{n_E}{}}}(E_m^{exact}E_m^{KS})^2+500\times {\displaystyle \underset{l,X}{\overset{n_G}{}}}\left({\displaystyle \frac{E_l^{KS}}{X}}\right)^2`$ (6)
We have to bear in mind that both the “exact” exchange-correlation potentials and densities were calculated using the TZ2P basis set. This, however, does not necessarily imply that the HCTH/147@TZ2P functional has a distinct advantage when we fit to these quantities since basis set convergence with angular momentum is much slower in ab initio methods than for density functionals. Nevertheless, the density calculations using the Brueckner Doubles method include core correlation. This might suggest that the TZ2P basis set has a slight advantage over the other basis sets when fitting it to ZMP potentials that have been obtained from BD densitiesBD . Unlike the cc-pVTZ and 6-311G basis sets, the TZ2P basis set is of triple-zeta quality in the inner-shell orbitals.

Ignoring for the moment the Pople basis sets that just differ in the polarisation component on the hydrogen, the following ordering of HCTH/147@ in $`Q`$ is observed:

TZ2P $`<`$ 6-311+G(3df,2pd) $``$ 6-311+G(3d,2p) $``$ aug-cc-pVTZ $``$ cc-pVQZ $``$ 6-311+G(2d,2p) $`<`$ cc-pVTZ $`<`$ 6-311G(2d,p) $``$ DFO2

All basis sets are ordered by their errors (unlike Table 1), with $`<`$ and $``$ corresponding to a difference of more than 5% or 20% between the functionals, respectively. From Table 1, the HCTH/147@TZ2P functional clearly yields the best overall results. The Pople large basis sets also exhibit surprisingly low errors, although some polarisation and diffuse functions are needed to achieve this performance. Interestingly, Dunning’s correlation consistent basis sets, which were optimised at the CISD level, show high errors despite having more basis functions. Even at the quadruple-zeta level, this “basis set functional” yields higher errors than both HCTH/147@TZ2P and HCTH/147@6-311+G(3d,2p). HCTH/147@DFO2, whose basis set was developed specifically for density functionals, shows an extraordinarily poor performance, probably due to built-in constraint of DFO2 that the different angular momentum functions share the same exponents.

When excluding the potentials at the triple-zeta level, $`Q`$ changes to $`Q_1`$, and the ordering of HCTH/147@basis set becomes:

6-311+G(3df,2pd) $``$ 6-311+G(3d,2p) $`<`$ cc-pVQZ $``$ TZ2P $``$ 6-311+G(2d,2p) $`<`$ aug-cc-pVTZ $``$ cc-pVTZ $`<`$ 6-311G(2d,p) $``$ DFO2

Here, the HCTH/147@TZ2P functional exhibits a somewhat higher error compared to the HCTH/147@6-311+G(3d,2p) functionals for Q1, yet still has an error comparable to the only quadruple-zeta basis set tested, HCTH/147@cc-pVQZ. The latter basis sets yield only slightly lower errors than when fitting to a TZ2P basis. For the double-zeta quality basis sets, the energies and gradients have much larger error contributions to $`Q`$ than the potentials, hence the exclusion of the latter has no impact on the ordering in the HCTH/147@ functionals:

6-31+G(2d,p) $`<`$ 6-31+G\** $`<`$ DFO1 $``$ 6-31G(2d,p) $`<`$ aug-cc-pVDZ $`<`$ 6-31G\** $`<`$ 6-31G\* $`<`$ DZP $``$ cc-pVDZ $`<<<`$ 6-31G $`<<<`$ 3-21G

As for the double-zeta basis sets, the unpolarised basis sets are clearly not very useful for the purpose. HCTH/147@DFO1 and the functionals fitted to the Pople basis sets yield the lowest errors. The former should not be surprising as it was developed specifically for DFT whereas the Pople basis sets were developed for Hartree-Fock. Nevertheless, with the exception of the HCTH/147@6-31+G(2d,p) functional, the double-zeta “basis set functionals” yield much higher errors than those obtained by the triple-zeta quality basis sets. Hence, if computationally feasible, basis sets of triple-zeta quality are preferable over basis sets of double-zeta quality when doing calculations employing DFT. Overall, and consistent with earlier observations made by Jensen Jensen2 , Dunning’s basis sets do not seem to be an optimal choice for density functional calculations. On the other hand, diffuse functions are more important with the Pople basis sets, in that they provide a significant error reduction. A detailed analysis of the results shows that predominantly the total atom energies, atomisation energies and electron affinities of anions are affected, with the error for the latter increasing by more than a factor of two when omitting diffuse functions. In comparison, with the HCTH/147@cc-pVTZ functional the error of the atomisation energies and ionisation potentials of the anions decreases by less than 20% when diffuse functions are included in the basis set. This is a general indication that despite their limited usefulness for DFT, Dunning’s basis sets are more “balanced” than the 6-31 family of basis sets with respect to diffuse functions. The HCTH/147@TZ2P functional yields a similar error as the HCTH/147@cc-pVTZ basis set for the anions, suggesting that diffuse functions might reduce its error further. Changing the basis functions on the hydrogen atom, as has been suggested in a prior assessment of different density functionals CH ; CHcomment , worsened the performance of the HCTH/147@TZ2P “basis set functional” in all tests.

Having obtained a number of functionals, comparing their coefficients in Table 2 yields further insight. Here, only the HCTH/147@3-21G functional obtained from fitting to the split-valence 3-21G basis set differs blatantly from the other functionals. Both the HCTH/147@DZP and HCTH/147@aug-cc-pVDZ functionals have large zeroth-order like-spin correlation coefficients ($`c_{C\sigma \sigma ,0}`$) probably thus compensating for the basis set error. Note that this coefficient also increases when going to larger molecules HCTH407 . The HCTH/147@6-31+G\** coefficients are closer to the ones obtained by the triple-zeta basis sets than the other double-zeta quality “basis set functionals”, confirming its lower basis set error. No clear conclusions on this effect can be made with the triple- and quadruple-zeta functionals, since both the HCTH/147@TZ2P (including the potential) and HCTH/147@6-311+G(3d,2p) functionals yield the lowest errors and not the HCTH/147@cc-pVQZ functional. The latter basis set has supposedly the lowest basis set error compared to the basis set limit Jensen2 . The exchange coefficients ($`c_{X\sigma ,n}`$ with m going from 0 to 4) of all the functionals at this level of basis set quality seem to be similar for the different functionals reported in Table 2. In contrast, the correlation parameters exhibit a wider variation upon switching between the fitting sets with 93, 147 or 407 systems HCTH407 . The effect of the diffuse functions on the basis set going from the HCTH/147@cc-pVTZ to the HCTH/147@aug-cc-pVTZ functional is not as significant as the change between the different “basis set functionals” of triple-zeta quality. Nevertheless, the coefficients still differ by a significant amount, suggesting that convergence towards the basis set limit has not been achieved. However, if similar basis set errors are absorbed into the parameterisation of functionals fitted to different triple and quadruple-zeta basis sets, the functionals are expected to be fairly transferable. In addition, the functional error is only very slightly influenced by changing some of the functionals coefficients. We have to note that using functionals with quadruple-zeta quality or higher is not useful, since the basis set error is probably about a magnitude lower than the functional error itself.

In Tables 3 to 6 we investigate the transferability of the “basis set functionals” in more detail. Here, all the new functionals developed are evaluated with the 6-31+G\**, aug-cc-pVDZ, 6-311+G(3d,2p) and TZ2P basis sets. The results for the 6-31+G\** basis for selected functionals that are listed in Table 1 are shown in Table 3. As expected, the HCTH/147@3-21G and HCTH/147@6-31G functionals are clearly insufficient to describe the properties investigated. Reasonably low errors are obtained with all functionals being fit to the 6-31G\* basis set or larger. Interestingly, in the case of the 6-31+G\** basis set, only the HCTH/147@DZP and HCTH/147@aug-cc-pVDZ functionals show $`Q_1`$ errors below 40 a.u., affected mainly by the gradient error contribution. In comparison (Table 1), the functional optimised for this basis set has a Q1 value of 32.5 a.u.. For all the triple-zeta quality functionals, the RMS energy error is lower, but the gradients (geometries) counterbalance this effect. All of these “basis set functionals”, with the exception of HCTH/147@6-311+G(3df,3pd) and HCTH/147@aug-cc-pVTZ, give very similar $`Q_1`$ errors, with RMS energy errors of 6.7 $`\pm `$ 0.3 kcal/mol and a gradient error around 3.6 a.u.. Overall, the triple-zeta and quadruple-zeta quality functionals vary by at most 20% in both errors.

The evaluation of the functionals with the aug-cc-pVDZ basis set yields similar results (Table 4). Of course, the overall error $`Q_1`$ is higher, as is evident from the larger error obtained by the HCTH/147@aug-cc-pVDZ functional compared to the HCTH/147@6-31+G\** functional displayed in Table 1. The gradient errors are especially affected. Only the double-zeta quality “basis set functionals” yield errors that can be compared to the $`Q_1`$ value of 42.4 a.u. obtained with the HCTH/147@aug-cc-pVDZ functional (Table 1). Here, the HCTH/147@6-311+G(3d,2p), HCTH/147@6-311+G(3df,2pd), HCTH/147@TZ2P, HCTH/147@aug-cc-pVTZ and HCTH/147@cc-pVQZ functionals again yield higher errors, mainly because of the gradient (and partly due to the energy) contributions to the error $`Q_1`$. Thus, the functionals fit to the double-zeta quality basis sets give the lowest errors when applied to other basis sets of double-zeta quality.

If we evaluate the functionals with triple-zeta quality basis sets, we observe, as expected, that their errors are a lot closer to the lowest error possible obtained by the basis sets to which the functionals were fit. When using the 6-311+G(3d,2p) basis set, the double-zeta “basis set functionals” yield the lowest gradient errors, although again this is only achieved when including diffuse diffuse functions. As for the other entire basis set evaluations in Tables 3 to 6, the double-zeta basis sets give a much lower error for the gradients, but not for the energy. This is an interesting phenomenon; the same observation was made in a different contextMenconi  when GGA functionals were developed solely for the description of accurate structures and frequenciesMenconif . Hence, it is probably a lot harder to develop a single functional for both accurate geometries and energies rather than one for each individual property. However, the former approach of separating the calculations can lead to other problems, namely that energetic properties will be calculated at non-equilibrium structures. The HCTH/147@cc-pVQZ, HCTH/147@cc-pVTZ and HCTH/147@aug-cc-pVTZ functionals again yield higher errors than the other functionals developed with a triple-zeta quality basis set. Noteworthy is the low error of the HCTH/147@TZ2P functional when applied to the 6-311+G(3d,2p) basis set, resulting in a $`Q_1`$ value of 15.5 compared to the minimum value of 15.0 a.u. obtained by the HCTH/147@6-311+G(3d,2p) functional. This is even lower than the error of the HCTH/147@6-311+G(3df,2pd) functional obtained with the 6-311+G(3d,2p) basis set. This emphasises the transferability of the functionals that have been fitted to the higher basis sets, since the variance is less than 10% on $`Q_1`$ and on the RMS energy error and the gradient errors (with the exception of the HCTH/147@cc-pVTZ and HCTH/147@6-311+G(2d,p) functionals for the RMS energy error). A number of the functionals parameterised for the triple-zeta basis sets yield a lower error when evaluating $`Q_1`$ with the 6-311+G(3d,2p) basis set rather than for the basis set to which they were fit (compare to Table 1). This indicates that a large amount of the remaining basis set errors absorbed into these functionals is of a similar magnitude.

In Table 6, the performance of the different “basis set functionals” with the TZ2P basis set is shown. All errors are slightly larger (by about 10%), but the trends are again the same. The HCTH/147@DZP basis set functional now gives a lower error than the HCTH/147@6-31+G\** functional, but only by a small margin. Again, the functionals fitted to Dunning’s basis sets of double-zeta quality yield larger errors than the other functionals fitted to basis sets of double-zeta quality. Among the triple-zeta level “basis set functionals”, the errors vary by only 10%, supporting the observations made for the 6-311+G(3d,2p) basis set. When evaluating the HCTH/147@3-21G and HCTH/147@6-31G functionals in Tables 4 to 6, the former functional has a lower gradient error. Hence, polarisation functions are essential when calculating geometries, reducing the gradient errors by a factor of two and more.

In summary, the triple-zeta “basis set functionals” are transferable between each other, indicating that when fitting using basis sets of this quality, the basis set error absorbed in the parameterisation does not play a role. The variance of 10% or less when evaluating one functional fit with a certain basis set with a different one is probably not important. This has further implications in the use of density functionals. Since the errors do not change significantly when going to higher basis sets, the triple-zeta basis set level is likely to be sufficient for use in density functional calculations. With still larger basis sets, the basis set truncation error will “drown in the noise” that is the inherent error of the functional itself. In general, Dunning’s correlation-consistent basis sets developed for correlated ab initio methods yield higher errors than the various Pople basis sets or TZ2P. When we investigated this behaviour in more detail, we found that the inclusion of core-valence basis functions CWNZ ; MartinCV  only lowers the RMS energy error by 0.1 kcal/mol. Inclusion of an additional tight d-functionWilsonD , which has been shown to be important for second-row elementsSO2martin , reduces this error by a further 0.3 kcal/mol. Furthermore, these contributions lower the sum of the gradient error by another 0.1 a.u. If we assume that the same contributions were added to the HCTH/147@aug-cc-pVTZ functionals results, we would probably arrive at a $`Q`$ value close to the lowest value obtained. Thus, although several enhancements for correlation-consistent basis sets could be introduced, extensive complements of higher-angular momentum functions, in particular, do not appear to be necessary. Out of all basis sets tested, the HCTH functional fit to the 6-311+G(3df,2pd) basis set gives the lowest energies and gradients errors, and the HCTH functional fit to the TZ2P basis set yields the lowest errors for the quantity $`Q`$, which includes energies, gradients and potential points. We believe that this is due to correlation effects described by the exchange-correlation potential in the core that cannot be adequately described by the 6-311+G basis sets.

## III Hybrid Functionals and Basis Sets

The method applied in the last section to the GGA functionals is now used for hybrid functionals to assess the different basis sets. There is, however, one complication: the varying amount of exact Hartree-Fock exchange. Here, we fitted hybrid functionals to the 6-31+G\**, aug-cc-pVDZ, 6-311+G(3d,2p), TZ2P and aug-cc-pVTZ basis sets, using the abovementioned procedure. While the inclusion of the exchange-correlation potential points into the fits for hybrid functionals is possibleB972 , it is not clear if the overall performance of these functionals is generally betteranharm , as we will also see in the final section. Therefore, we restrict ourselves to fitting to gradients and energies only (similar to the determination of $`Q_1`$ in the last section, see Eqn. 6). All functionals were again fit to the 147 systems, but with m=2 in the power series in Eqn. 2. This cut-off in the power series expansion yields the B97-1 form rather than the HCTH form. Here, 9 linear coefficients are fit instead of the 15 in the GGA. The amount of exact exchange was varied over a range from 0% to 50% in order to determine the minimum.

Before discussing the outcome of the results of the hybrid functionals, it is worth comparing the values at 0% Hartree-Fock exchange to the functionals obtained in the previous section. The difference between them is the fit to exchange-correlation potentials, which will raise the $`Q_1`$ error for the GGA functionals in contrast to the “hybrid” functional at 0%. Of course, another discrepancy is the different number of coefficients. The newly obtained GGA functionals are displayed in Table 7, and can be compared (with the differences mentioned) to those in Table 1. For the TZ2P, 6-31+G\** and aug-cc-pVDZ basis sets, the higher-order coefficients seem to change very little (see table 1). While the energy error increases going from m=2 to m=4, the gradient error decreases, with $`Q_1`$ close to the errors when fitting only to gradients and energies with m=2. This has been previously reportedCH ; B97  and has led to the conclusion that it is unnecessary to include orders higher than 2 in the power series. Here, mainly the error in the potential is affected as an additional calculation with the TZ2P basis set and the variable $`Q`$ shows. In this case, $`Q`$ for the HCTH/147 with m=4 form yields 37 a.u., compared to 51 a.u. for the B97-1/147 form with m=2, hence the overall $`Q`$ value is raised by 40%. Comparing Tables 1 and 7, we realise that not only is the functional somwehat dependent on the basis set for which it was parameterised, but also the basis set dependence itself is dependent on the class of functional used.

Let us now return to the hybrid functionals. Since the error of the resulting functional also correlates slightly with the starting guess used, the points will not necessarily fit a curve. Unlike the GGA functionals where we can determine the starting guess by fitting to exchange-correlation points, with hybrid functionals the initial coefficients are, at best, educated guesses. Thus, we sometimes had to fit a curve through the points, and in Figure 1, the RMS energy error is plotted as a function of amount of exact exchange for various functionals. This plot shows a disturbing property of hybrid functionals: the amount of exact exchange in the functional depends on the basis set for which it is fit. All hybrid density functionals known have their exchange coefficient fit to a specific basis set (or are using numerical DFT), hence the variation of the exact exchange, in the range of 15% to 25% might well be due to basis set effects. In the case of the tested functionals, the optimised exact exchange-fractions for all three triple-zeta basis sets (obtaining B97-1/147@6-311+G(3d,2p), B97-1/147@TZ2P and B97-1/147@aug-cc-pVTZ) hover around 18%. However, fitting to a double-zeta basis set (B97-1/147@6-31+G\** and B97-1/147@aug-cc-pVDZ) yields minima located around 28%. In Figure 2, the value of $`Q_1`$ is displayed, exhibiting generally the same behaviour as the energy in Figure 1. Since this is the value used in the fit (with more sophisticated weights), it is more informative, albeit more abstract. Here, the differences between the basis sets, which were already exhibited with the GGA functionals, can be seen again. While the B97-1/147@6-31+G\** and the B97-1/147@aug-cc-pVDZ functionals show similar minima in the energies, with the B97-1/147@6-31+G\** curve shifted a bit towards lower values. Still, the gradient error when fitting to the 6-31+G\** basis is much lower than with the aug-cc-pVDZ basis set. Similar behaviour is observed for the basis sets of triple-zeta quality. Only when the gradients are included does the B97-1/147@TZ2P functional yield a lower error for its minimum than B97-1/147@aug-cc-pVTZ, with the values obtained by the B97-1/147@6-311+G(3d,2p) functional marginally lower still. Detailed results of the minima are given in Table 7 displaying the lowest energy points calculated. The minima obtained when fitting a separate curve through the points are at 29 and 28% for the B97-1/147@6-311+G\** and B97-1/147@aug-cc-pVDZ “basis set functionals”, and the minima for the B97-1/147@6-311+G(3d,2p), B97-1/147@TZ2P and B97-1/147@aug-cc-pVTZ “basis set functionals” are at 16, 17 and 17%, respectively.

## IV Assessment of Density Functionals

In this section we assess several density functionals with two basis sets: cc-pVTZ and TZ2P. The former basis set has been used in the construction of VSXCVSXC , and the latter in the parameterisation of the HCTH-type functionalsHCTH93 ; HCTH120 ; HCTH407 ; B972 ; tHCTH . For this evaluation, we use the large 407 Setwww , again comparing RMS energy, gradient and $`Q_1`$ errors. A similar comparison for most of the functionals evaluated has been done before, using the TZ2P basis set and the much smaller sets of 93 and 147 molecules. Thus, we can assess if the results will remain transferable between the G2-1 and G3 setsCH . Table 8 shows all results obtained.

* The simplest density functional method, LSDA in its VWN parameterisation, is already a vast improvement over the Hartree-Fock method. The LSDA geometries are surprisingly accurate, even more so when comparing the gradient error to a vast number of GGA and meta-GGA functionals.
* For molecular systems, PKZB yields no clear improvement over PBE, although it was developed as an improvement over PBE including a semi-empirical fit and an extra variable, the kinetic energy density $`\tau `$. Both functionals give RMS energy errors close to 20 kcal/mol. In addition, we would expect geometry errors similar to the ones obtained by the LSDA method. We would discourage the use of these functionals for the calculation of both thermochemical data and geometries.
* mPW1K is a hybrid functional with a large percentage of exact exchange, developed for accurate reproduction of reaction barriers, reducing the error of B3LYP for this property by about 50%mPW1Keval . However, its performance for minimum geometries and energetics of stable molecules is the worst of all hybrid functionals tested and even standard GGA’s yield lower errors.
* PW91PW91, the original GGA proposed by Perdew, yields slightly lower errors than PBE, that largely come from the evaluation of the gradients.
* When replacing the PW91 exchange functional with mPW91, the RMS energy error is reduced, yet the gradient error increases.
* The BP86 GGA functional, albeit an improvement over these methods, still gives a $`Q_1`$ error almost three times larger than the best functionals available. The RMS energy error is lowered by 5 kcal/mol when using the cc-pVTZ basis set instead of TZ2P; the gradient error however increases.
* BLYP is one of the more accurate functionals, but generally overestimates bond distances, and its gradient error is extremely large. This leads to RMS errors in bond lengths that are almost twice as high as for hybrid functionalsHCTH407 . Nevertheless, it is still one of the commonly used functionals, very often employed when calculating hydrogen bondsTBH  or in Car-Parrinello Molecular Dynamics simulations of liquid waterBLYPCP . Its RMS and gradient errors both increase when using the cc-pVTZ basis set. In contrast to the smaller sets, the BLYP functional now clearly outperforms BP86 for the 407 set, whereas this was not the case for the 93 setCH .
* The BPW91 functional uses Becke’s exchange functional in combination with PW91 correlation, and yields lower errors than the PW91 and mPW91 exchange functionals for the properties tested. Here, the functional gives a performance similar to BLYP, perhaps slightly better. Again, its gradient error is higher for the cc-pVTZ basis set.
* PBE0 does not yield better energy predictions than the GGA functionals mPW91PW91, BPW91 and BLYP. Interestingly, its gradient error is lower than the one obtained for the most commonly used functional B3LYP, thus we would expect a pretty accurate description of geometries and higher-order properties.
* OLYPOLYP  is a clear improvement over BLYP for atomisation energies and reactions, and even more so for molecular structureshoe ; baker . Its overall errors are almost comparable to those of the hybrid functionals, much better than BLYP and BPW91. Here, an improved exchange functional OPTX is used instead of Becke’s exchange functional. However, its performance in hydrogen bonds is not as good as BLYP.BHunpub
* B98 does exactly the opposite of PBE0: We would expect its geometries to be further away from the equilibrium than B3LYP, but its energy error is slightly decreased. In the overall $`Q_1`$ evaluation, the error is very similar to the one obtained by the older B3LYP.
* The B3LYP functional is probably the most widely used hybrid functional, and every new functional is compared to its accuracy. The RMS error for the 407 Set is close to 10 kcal/mol, which is considerably large compared to a “chemical accuracy” of 2 kcal/mol. The gradient error is again slightly increased for the cc-pVTZ basis set compared to TZ2P.
* B97-2 is a reparameterisation of B97-1 including the ZMP potential points into the fit of the exchange-correlation functional. Hence, we can expect its energy error to be worse than B97-1, since the fit to an extra quantity usually worsens the performance to the energies. However, its advantages over B97-1 still have to be established since even its gradient error is larger.
* The performance of VSXC for the 407 set is very similar to B3LYP. This contrasts with the bad performance of VSXC for the smaller 93 Set where it returned an error barely lower than the BLYP functionalCH . Fitted to the cc-pVTZ basis set, it is the only functional yielding similar errors for both basis sets.
* All HCTH functionals give errors which differ in $`Q_1`$ by less than 3% for the TZ2P basis set and less than 8% for the cc-pVTZ basis set. These functionals, while they are pure GGA functionals, yield errors that can be compared to hybrid functionals like B3LYP for both the TZ2P and cc-pVTZ basis sets. Here, the additional value added to the functional by reparameterising it to 147 or 407 systems is not obvious, although the HCTH/407 functional outperforms the other parameterisations in the RMS energy error. The $`Q_1`$ value of HCTH/407 is the worst for the HCTH functionals with the cc-pVTZ basis set. The justification for the reparameterisation, making HCTH/407 a better functional than HCTH/93 or HCTH/147, will become only visible when considering hydrogen bonds or inorganic moleculesTBH ; Menconi  where error cancellation plays an important role. Generally, this error cancellation cannot be expected in DFT methods. In some post-Hartree-Fock methods like MP2 it is inherent, and hence the reparameterisation of the functionals remains important to recapture such effects.
* B97-1 is probably the best choice when it comes to using density functional hybrid calculations, since it is already well tested and its calculated structures are similar to those obtained by B3LYP (or even slightly better), and it outperforms B3LYP by about a third in $`Q_1`$ when predicting energetic properties.
* $`\tau `$-HCTH is for both basis sets an improvement over the HCTH functional showing that the inclusion of the kinetic energy density can lower the error further. Unfortunately, it lacks the performance of the HCTH/407 functional for weak interactionstHCTH ; HCTH407+ .
* The $`\tau `$-HCTH hybrid, when additionally including exact exchange, clearly yields the lowest $`Q_1`$ value of all methods tested. Its error for the TZ2P basis set is 50% lower than the one obtained by B3LYP in addition to yielding a lower gradient error and structurestHCTH  than B3LYP.

Summarising the results in table 8, by ranking all GGA functionals based on their $`Q_1`$ value, we get the following order (taking into account that some functionals have been fit to one of the basis sets):

HCTH/147 $``$ HCTH/93 $``$ HCTH/407 $`<`$ OLYP $`<`$ BPW91 $`<`$ BLYP $`<`$ mPW91PW91 $`<`$ BP86 $`<`$ PW91PW91 $`<`$ PBE

For the meta-GGA’s and hybrid functionals, the ordering is:

$`\tau `$-HCTH hybrid $`<`$ B97-1 $`<`$ $`\tau `$-HCTH $`<`$ VSXC $``$ B97-2 $``$ B3LYP $``$ B98 $`<`$ PBE0 $``$ mPW1K $`<`$ PKZB

Although RMS errors can give a lot of insight, the maximum errors are also considered important. In light of this, we have examined in Table 9 the number of molecules in the 407 set for each functional that have large RMS energy errors (over 15 kcal/mol). As we have discussed in the basis set evaluation, the cc-pVTZ basis set generally exhibits more outliers (molecules with atypically large errors) than the TZ2P basis set. The results are very similar to the ones obtained in Table 8, and compared to B3LYP, the best available functional (the $`\tau `$-HCTH hybrid) cuts the number of outliers in half, yielding a considerable improvement. With regard to the accuracy of density functional theory with all the evaluations which are done with the 407 Set, the RMS error of the functionals ranges between 6 ($`\tau `$-HCTH hybrid) and 21 kcal/mol (PBE), with the most commonly used hybrid functional B3LYP yielding an error around 10 kcal/mol. The mean error of these functionals is between 3.8 kcal and 15.6 kcal mol, with B3LYP yielding 6.2 kcal/mol. This can be compared to empirical correction methods like G3G3 , which yield errors around 1 kcal/mol for the G3 set that, however, does not include some of the molecules with the largest errors in our 407 set. Extrapolation methods like W2 W2  give mean errors around 0.5 kcal/mol for a considerably smaller set similar to our 147 Set. Full MP2 yielded an RMS error of about 20 kcal/mol for the 407 Set using the TZ2P basis set. This generally places the accuracy of DFT between raw perturbation theory results and coupled-cluster theory when calculating ground state energies.

## V Conclusions

From the above, we can make a number of observations concerning the use of DFT functionals and the basis sets used with them. All results regarding the basis sets are obtained by fitting functionals to basis sets, and then evaluating their errors, thus the results are independent of the functional parameterisation. For most of the properties investigated, including energies and gradients, the Pople basis sets can be recommended. Dunning’s basis sets, fit to CISD, give considerably higher errors despite having a larger number of basis functions. Although it yields a higher error than the 6-311+G(3df,2pd) basis set for energies and gradients, the TZ2P basis set still gives the lowest overall error for the GGA functionals when including the ZMP exchange-correlation potentials into the fit. The basis set error, which might still be significant at the triple-zeta level, shows that the functionals obtained by fitting to one basis set are transferable to other basis sets. Hence, it is probably not important to reach the basis set limit when developing new density functionals, since the overall DFT error is considerably larger. Basis sets developed for DFT methods might alleviate this problem, but the problem remains as to which functional to use for their development. The same analysis for hybrid density functionals shows that the amount of exact exchange obtained is dependent on the basis set itself. Whereas basis sets of double-zeta quality yielded minima around 28%, the triple-zeta basis sets evaluated had their minima around 18%. Concerning the difference between the hybrid “basis set functionals”, the same conclusions as for the GGA functionals can be drawn. The same trends are also visible when evaluating several other published functionals to a large test set. We can deduce that several hybrid functionals, such as B97-1 and the $`\tau `$-HCTH hybrid functional render errors which are significantly lower than the ones obtained by B3LYP. As a pure GGA functional, the HCTH functional types give errors comparable to B3LYP for the investigated properties.

Nevertheless, the accuracy of modern density functional theory cannot be compared to that of ab initio extrapolation methods. The best functional tested yields an RMS energy error as large as 6.3 kcal/mol for a large set of molecules, which is still far away from the desired “chemical accuracy” of 1-2 kcal/mol.

## VI Acknowledgements

A. D. Boese is grateful for financial support by the EPSRC, the Gottlieb Daimler- und Karl Benz- Stiftung and the Feinberg Graduate School. This research was supported by the Lise Meitner Center for Computational Chemistry and the Helen and Martin A. Kimmel Center for Molecular Design.

Fig. 1:
RMS error (in kcal/mol) of the 147 Set with different hybrid “basis set functionals”.

Fig. 2:
$`Q_1`$ error (in atomic units, see eqn. 6) of the 147 Set with different hybrid “basis set functionals”.


format: text

text:
# The three dimensional structure of the giant stellar stream in Andromeda

## 1 Introduction

Stellar streams are the remnants of structures that have become perturbed or disrupted due to the action of gravitational tides. In the hierarchical formation model of galaxy haloes, smaller systems merge first and form the larger galaxies we see today (White & Rees 1978). As such, these structures are expected to be common and examination of their properties gives valuable information about the galactic formation process. In the outer regions of haloes, dynamical times are sufficently long to leave many structures detectable over much of the age of the universe (eg. Johnston, Hernquist & Bolte 1996) and so these streams can trace the left-overs of even ancient accretions. Indeed, under some favourable conditions the galaxy potential can leave massive streams essentially intact for all cosmic time.

Within our Local Group there is a lot of observational evidence that galaxy mergers are an on-going process. The discovery of the Sagitarius dwarf galaxy by Ibata, Gilmore & Irwin (1994) revealed that our own Galaxy is undergoing strong interactions with its companions. Further observations detected a tidal stream associated with this galaxy’s disruption (eg. Ibata et al. 2001a). Indeed, several of our Galaxy’s globular clusters have been found to have originally been associated with this dwarf galaxy (Ibata, Gilmore & Irwin 1994; Bellazzini, Ferraro & Ibata 2003). One of the most favoured explanations for the formation of the thick disk is that it is due to a violent merging event (eg. Quinn, Hernquist & Fullager 1993, Schwarzkopf & Dettmar 2000, Feltzing et al. 2003). Disentangling the merger history of galaxies is thus of prime importance and streams are ideal tracers of this process.

The discovery of a giant stellar stream in the outer regions of our neighbouring giant galaxy M31 was reported in two earlier contributions in this series (Ibata et al. 2001b, Ferguson et al. 2002). It is an apparently vast structure and our initial INT WFC survey did not go out far enough to find its true extent. Here we report on our findings from a deeper survey using the CFH12K camera that probes further out and to fainter magnitudes than we had previously surveyed.

## 2 Observations

During the nights of August 23, September 13-14 and September 17-18 2001, fourteen fields were oberved with the 12k$`\times `$8k CFH12K camera, at the prime focus of the Canada-France-Hawaii telescope. The observations were kindly taken by CFH Telescope staff, in service mode. Conditions were generally photometric, with good seeing of $`0\text{′′}\text{.}8`$. Three sets of Mould V-band exposures and three Mould I-band exposures, each of 545s duration, were secured per field. A negligible colour correction is involved to convert these filters to the standard Johnson-Kron-Cousins system 111http://www.cfht.hawaii.edu. The data were preprocessed by CFH staff with the CFH pipeline software, to correct for bias offset, flat-fielding and fringing. We then applied the same object detection, classification, photometry and catalog generation algorithms as we successfully used before on our INT-based survey (Irwin & Lewis 2001).

The fields were selected to lie along the previously detected region of the Andromeda stream, and we extrapolated along this stream in a straight line $`3`$ degrees further out to the south-east and to the north-west. The arrangement of the survey regions presented here is displayed in Figure 1, and the RA and Dec of fields 1–8 and 12–14 are listed in Table 1. The photometry of the three heavily-crowded fields 9–11, which lie close to the center of M31, will be presented elsewhere. Field 14 is used as a reference field, as no stream component was found there.

## 3 Stellar populations of the stream

A typical field from our survey (field 6) that illustrates the quality of the CFHT data is shown in Figure 2. This field is one of several that overlaps the INT survey region, although the CFHT data extends a full magnitude deeper than our previous survey. The overlap enables us to make an external comparison of the photometric systems and calibration. The internal overlaps between (most) CFH fields and their external overlap with some of the INT data allowed us to update the calibration of two CFH fields (2,3) that were taken under slightly non-photometric conditions. After these adjustments the photometric calibration systematic errors are at the level of $`\pm `$2$`\%`$ over the entire region surveyed. Overlaid on the colour-magnitude diagram are four well-studied globular cluster sequences (NGC6397, NGC1851, 47Tuc and NGC6553) from Da Costa & Armandroff (1990) and Sagar et al. (1999) that span a range of different metallicities (left panel), and also four evolutionary tracks (for $`\alpha `$-enhanced $`0.8\mathrm{M}_{}`$ stars) from VandenBerg et al. (2000) (right panel). A wide range of stellar populations are evident in this field spanning the full range of metallicities represented by the fiducial sequences ($`0.2[Fe/H]1.9`$). Inspection of a sequence of these diagrams shows evidence for a metal-poor halo population with a metallicity spanning roughly $`[Fe/H]2`$ to $`[Fe/H]0.7`$ and, in addition to this, a more metal-rich population is evident in the majority of the fields. Using the fiducial sequences to define the metallicity distribution function we can attribute a mean metallicity of $`[Fe/H]0.5`$ to this component, with a dispersion of $`0.5`$dex. This population is clearly seen to extend out to field 1 in our CFH survey data, where it still remains a numerous component: around 1000 stream stars are detected in field 1 to I$`=23.5`$. This is in contrast to the metal poor halo component that virtually disappears by this galactocentric distance ($`4^{}`$).

The evolutionary tracks we have selected correspond to $`[\alpha /Fe]=0.3`$ since CNO and other $`\alpha `$-process elements, which dominate the heavy element composition, are thought to be enhanced relative to solar for Galactic halo stars (eg. Wheeler, Sneden & Truran 1989, Carretta, Gratton & Sneden 2000). These tracks are seen to be in good agreement with the globular cluster sequences and our data. In particular, the I magnitude of the tip of the red giant branch (TRGB) of the tracks appear to match up well as a function of metallicity. We will make use of this property in Section 4.2 in order to determine the relative distance of the metal-rich population to the bulk of the stellar population of M31.

## 4 Three dimensional structure of the Andromeda Stream

### 4.1 Stellar profile across the stream

In the INT data, the stream appeared as a fairly straight structure pointing away to the SE from the center of the galaxy. With our extra CFH data we are now in a position to better constrain the projected profile of the stream on the sky. A reference line was chosen that follows a path parallel to the direction of the CFH survey fields and that passes through the centre of M31. The density of sources perpendicular to this line satisfying $`20.5<\mathrm{I}<22.5`$ and $`\mathrm{V}\mathrm{I}>2`$, for fields 1 to 8, is displayed as a solid line in Figure 3. This figure demonstrates that there is a rapid rise in number counts from the NE side of the stream towards the SW side, with the density of sources increasing by a factor of $`2`$ over $`0.5^{}`$. However, after the peak, the distribution is almost flat up to the SW edge of our survey area. This is consistent with the sharp edge to the stream seen in the INT survey, although at this time we cannot constrain the behaviour of the SW edge of the stream. To investigate the radial dependency we divided the CFH survey data into 4 bins (dotted lines). In Figure 3, the shallowest histogram is for fields 1 – 3, next is for fields 4–5, then 6–7 and finally 8. The distribution of sources in each of these bins can be observed to have a similar shape and, in particular, we find that the position of the peak of the distribution varies by less than $`0.2^{}`$ over the $`4^{}`$ line between fields 1 and 8. However, the spatial profile of the stream beyond $`d=2.4^{}`$ is broader, showing that the stream of width $`0.5^{}`$ appears to fan out into a wider distribution in the outer regions.

Likewise, by separating the stars in Figure 3 with $`d<0.9^{}`$, we find that there is no significant difference in the shape of the histograms with that for $`d>0.9^{}`$. We had expected that for $`d<0.9^{}`$ the contribution of the stream would be small compared to the halo, and so the histograms would be flatter in appearance. Instead, it would appear that the stream contamination is still significant in the $`d<0.9^{}`$ region.

Figure 4 shows how the number density of these stars varies along the path defined above (fields 9, 10 and 11 have not been used due to crowding problems). As expected, the star counts increase as we approach the centre of M31. However, there is an excess of stars clearly visible for $`d<1.5^{}`$, which is highlighted in Figure 4 by comparison with the profile for $`d>0^{}`$ (its reflection about $`d=0^{}`$ is shown as a dot-dashed line). The CFH fields were selected to highlight the stream which is clearly delineated in azimuth (Fergusson et al. 2002) thus we attribute this excess to the stream. Although we believe stream to be present at $`d>0^{}`$ (Section 4.3), its stellar density by this point is much smaller compared to other M31 components.

### 4.2 Stream distance relative to M31

The I band magnitude of the tip of the red giant branch (TRGB) is now recognised as a good standard candle for metal poor ($`[Fe/H]<0.7`$) old ($`>2`$ Gyrs) stellar populations (eg. Salaris & Cassisi 1997). However, as we go to redder, more metal rich populations, the luminosity decreases due to the effect of increased opacity in the stellar atmosphere. This is clearly observed in our colour magnitude diagrams. As such, in order to gain a distance estimate of the stream relative to M31, we need to calibrate the redder ($`\mathrm{V}\mathrm{I}\mathrm{}>2.2`$), metal-rich stream to the bluer ($`\mathrm{V}\mathrm{I}\mathrm{}<2.2`$) M31 component and therefore need to know how the TRGB behaves a function of colour.

Vandenberg et al. (2000) present evolutionary tracks for a range of stellar masses ($`0.5\mathrm{M}_{}1.0`$), with a range of metallicities ($`2.31[Fe/H]0.30`$). A selection of these tracks are shown in Figure 2. In the left panel of Figure 5 we show the predicted location of the TRGB for $`0.8\mathrm{M}_{}`$ stars with $`[\alpha /Fe]=0.3`$ as a function of metallicity for these tracks. This is overlaid on one of our INT WFC survey fields (#76 located at $`\xi =1.068^{},\eta =0.386^{}`$) from a central region of M31 with a similar colour spread on the red giant branch. Evidently there is good agreement between these theoretical evolutionary tracks and our M31 data. The INT passbands have been converted to V and I using $`I=i^{}0.101\times (\mathrm{V}\mathrm{I})`$ and $`\mathrm{V}=V^{}+0.005\times (\mathrm{V}\mathrm{I})`$. These transformations have been derived by comparison with observations of several Landolt standard fields 222http://www.ast.cam.ac.uk/$``$wfcsur/colours.php.

In order to calibrate this model to the data, we adjust every I magnitude for stars in the M31 field using the differential model TRGB magnitudes with respect to a fiducial colour of $`\mathrm{V}\mathrm{I}=1.6`$. We then construct luminosity functions in the ranges $`1.4<\mathrm{V}\mathrm{I}<2.2`$ and $`2.2<\mathrm{V}\mathrm{I}<2.8`$ to represent the metal poor and metal rich components respectively. Using a data adaptive least-squares technique (McConnachie et al. in preparation), we measure the relative locations of the two TRGBs. The model-corrected metal-rich TRGB is systematically $``$ 0.1 mags brighter than the metal-poor TRGB and defines an empirical correction to the model-corrected TRGB magnitude.

We now apply this metallicity correction to field 8 in our CFH data, as it is both the closest field to M31 and has the highest signal. The colour magnitude diagram of this field is shown in the right panel of Figure 5. By applying the magnitude transformation on the stars and by taking into account the 0.1 mag difference between the redder and bluer TRGBs that we would expect were the components all at one distance, we conclude that in field 8 there is no significant difference between the distance of the stream component and that of the bulk of M31.

### 4.3 Relative distance change along the stream

With a link between the stream distance and the distance of M31 in field 8, it is now possible to examine how this distance changes with respect to M31 as a function of stream position. To do this we cross-correlate the stream-dominated I-band luminosity function ($`2.5<\mathrm{V}\mathrm{I}<3.5`$) in each field with the equivalent luminosity function in field 8. Any shift measured is a direct indicator of the stream distance. Examples are shown in Figure 6. In all cases field 14 is used as a reference field and a scaled (smoothed) version of its luminosity function is subtracted from each field analysed to account for varying foreground Galactic contamination. The luminosity functions are constucted using only objects classified as stellar sources in both the V and I filters. Although for the redder objects significant incompleteness due to the V-band sets in at I $``$ 22, this is sufficiently beyond the peak of the red star luminosity function to not affect the cross-correlation, even with slight variations in field-to-field depths. Further tests using I-band only stellar detections satisfying $`\mathrm{V}\mathrm{I}>2.5`$ support this conclusion.

The results of the cross-correlation are shown in Figure 7 and are tabulated in Table 1. The main contributors to the errors in the TRGB distance determinations are: the $`rms`$ errors in the least-squares fit and cross-correlations which average about $`\pm `$0.03 magnitudes; the photometric calibration error of around $`\pm `$0.02 magnitudes; and model-dependent systematics from the various fits of order $`\pm 0.03`$. This leads to a total error of around $`\pm 20`$ kpc in the distances.

The projection on the sky and the distance variation along the stream point to an almost linear structure in the south-east region, aligned at roughly 60 degrees to the line of sight and extending out to distances well beyond 100 kpc from the centre of M31. No stream component is detected in field 14, although there is a clear signal in the inner half of field 13. The stream must therefore be wrapping around the centre of M31 toward the northeast. The distance variation between the two ends of our survey is $`140`$kpc, implying the true length of the surveyed region of the stream is $`>160`$ kpc.

## 5 Conclusions

The tidal stream discovered by our survey of M31 (Ibata et al. 2001) extends linearly over at least $`6^{}`$ of the sky. In fact, we still have not surveyed out far enough to find the far south-eastern end. By analysing the metal rich red giant branch luminosity function in the I band and comparing them to VandenBerg et al’s (2000) evolutionary tracks, we are able to measure a radial distance change of order 140 kpc between the two ends of our survey, with the stream extending from approximately 100 kpc behind to 40 kpc in front of M31. We also see evidence that the stream then proceeds to wrap itself around M31. This information confirms that the Andromeda stream is a gigantic structure, angled at approximately $`60^{}`$ to our line of sight. Additionally, by analysing the stellar profile across the stream we find that it remains a coherent structure, at least on the north-eastern edge. Its stellar distribution also appears to be fairly constant with increasing galactocentric distance. However, more observations are required in order to determine the true extent of this vast object.

## Acknowledgements

This work is based on data collected at the Canada-France-Hawaii Telescope supported by INSU at Mauna Kea, Hawaii. The research of AMNF has been supported by a Marie Curie Fellowship of the European Community under contract number HPMF-CT-2002-01758. We would like to thank Professor VandenBerg for supplying us with his evolutionary tracks.


text:
# A Monitoring System for the BaBar INFN Computing Cluster

## I INTRODUCTION

Large clusters with hundreds or thousands of nodes are very difficult to manage due to their size and the complexity of the applications they run. Computing farms are routinely used in the current generation of High Energy Physics experiments, given the huge amount of data to be processed.

An efficient monitoring system can be very helpful for profiling all components of the cluster. A *monitoring system* is a hardware or software component able to observe the activity of a system (see Jain (1991)). A monitor can observe the performances of a system, record statistics, analyze the data and display the results. Monitors are useful for many reasons, such as analyzing the resource usage of an application or identifying performance bottlenecks or usage patterns suggesting better algorithms. Monitors can be used to characterize the workload of a system, or to find the parameters to be used on a simulation of a real system. Moreover, they can be used to check for, and report hardware problems such as crashed or unresponsive hosts or broken communication links.

Continuous monitoring of computing clusters is a challenging problem for different reasons. The size (number of components) of the system to be monitored is often beyond the scalability limit of many available tools. These tools are usually tailored for a particular application; they are often closed products, so the user can’t adapt them to any variation of the requirements. Moreover, they are usually hard to configure, and provide an inconvenient user interface.

In this paper we present PerfMC (Performance Monitoring for Clusters), a monitoring system for large computing clusters. PerfMC makes heavy use of XML technologies (see W3 Consortium (2000)), and uses the Simple Network Management Protocol (SNMP) (Stallings (1999)) to gather data from any device containing an SNMP agent. Since SNMP is a standard protocol, implemented by many networked equipments, the tool is very general and can be used to monitor a wide range of devices. PerfMC embeds a WEB server which is used to generate time-series graphs from the collected data. Also, the WEB server can produce arbitrary XML pages by applying transformations to an internal XML representation of the cluster’s status. In this way it is possible to produce HTML status pages which can be displayed by ordinary WEB browsers.

The paper is organized as follows. Section II presents some previous works related to the problem of monitoring computing clusters. In Section III the design goals of PerfMC are illustrated and discussed. Section IV describes the architecture and implementation of PerfMC, and a case study is presented in Section V. Finally, Section VI reports the conclusions and future works.

## II RELATED WORKS

In recent years considerable attention has been devoted to the problem of monitoring the performances of clusters and distributed systems (see Ganglia Project ; Buyya (2000); Catania et al. (1996); King and Hunt (2000); Mansouri-Samani (1995); Saab et al. (2002); Subramanyan et al. (2000); Tierney et al. (2001); Uthayopas and Phatanapherom (2001); Wismüller et al. (1998)). In Tierney et al. (2001) the authors describe an agent-based monitor targeted primarily to GRID architectures (see Foster and Kesselman (1999)), which are wide area distributed systems where components can be connected to high-latency Wide Area Networks. The monitoring architecture is based on a producer-consumer paradigm, where individual monitors can subscribe for particular kinds of events, and receive notifications only when such events are generated by some producer.

It should be noted that a computational GRID is very different in size and complexity from a computing cluster, so implementing a monitoring system on them is different. A computational GRID is usually made of an heterogeneous collection of computing systems which are geographically distributed and connected through a WAN. A cluster is made of homogeneous machines, usually residing in the same room and connected together with a high speed LAN. These differences play an important role in defining the requirements for a monitoring system. For example, a monitor for a GRID should be built with security features, given that intrinsically insecure WAN links are used for communications. On the other hand, a cluster can be treated as a single, powerful machine. It should be protected with respect to the outside world, but communications among machines in the cluster can be unencrypted. If the monitoring infrastructure is built inside the cluster, there is no need to protect its control messages. Also, LANs are characterized by low latencies, high bandwidths and low packet loss rates. The size of a cluster is usually orders of magnitude smaller than that of a large GRID. For these reasons, monitoring systems developed for computational GRID environments, such as that described in Tierney et al. (2001), have been developed with very different requirements in mind with respect to a cluster monitor.

In Puliafito and Tomarchio (2000) it is proposed a monitoring and management architecture based on the use of mobile agents written in Java. Mobile agents allow management applications to be moved to the network devices, instead of moving the data provided by the network devices to the Management Stations. The approach based on Java mobile agents obviously requires the network devices to be equipped with a Java Virtual Machine, so that they are able to accept and execute code coming from the Management Stations. This is currently only possible when the monitored elements are general-purpose computers, as other devices are generally unable to run Java code (or any other user program).

The Ganglia monitoring system described in Ganglia Project  is a distributed monitoring system for clusters. It requires each node on the cluster to run a dæmon called gmond. It collects values from the local machine and broadcasts these values to all the other gmond processes running on the cluster. To limit the network utilization, broadcasts happen only when the changes in the observed values exceed a given threshold. gmond processes can also communicate with generic Ganglia clients by sending an XML status file over a TCP connection. Ganglia dæmons do not provide any facility to log the recorded data themselves, but rely on external programs to collect statistics, perform management actions when particular events occur, and display the status of the system.

Most monitoring tools use their own data collection protocol over TCP/IP links. One exception is SIMONE Subramanyan et al. (2000), which uses the standard SNMP protocol to build a large-scale, distributed monitoring system. Hierarchical monitoring has been employed in other systems as well, such as the one described in Uthayopas and Phatanapherom (2001). Such hierarchical, tree-based monitoring systems are particularly effective when the user is mainly interested in getting aggregate informations on the cluster’s status, such as the average load of all the machines, or the least utilized node of the cluster. This is because the information can be aggregated at each intermediate node of the hierarchy, thus avoiding the possible bottleneck of a single node getting all the data from all hosts. Unfortunately, this strategy does not help when it is necessary to continuously record the values of some parameters for every single machine, for example for producing graphs showing the variation of interesting quantities over the time.

Supermon, described in Minnich and Reid (2001), is a centralized monitoring system, yet it allows efficient and frequent data collection from the nodes of a Linux cluster. The Linux kernel has to be patched for the addition of a new system call which provides status informations. A server program running on each machine collects these informations and can pass them to requesting applications using a telnet-based network protocol. A possible drawback of this approach is the necessity to use a modified kernel on the monitored machines, and the fact that it is necessary to modify the implementation of the system call if additional parameters need to be monitored.

## III PerfMC DESIGN GOALS

The BaBar Italian reprocessing farm is in production since summer 2002. The farm, hosted at INFN Padova, is made of about 200 dual CPU Linux/Intel machines, and includes a tape library with a capacity of 70TB, and 24 TB of disk space. The farm is used to process data collected at SLAC by the BaBar detector.

We identified a number of requirements for the monitoring system of the cluster, which are now briefly discussed.

Intrusion-free. The monitor must be guaranteed not to have negative impact on the correctness of the monitored system’s results. It should be noted that this cannot be achieved if the monitored system strongly depends on hard real-time constraints to operate correctly. For such systems, even a small overhead induced by the monitoring activity could affect its results. The BaBar reprocessing farm does not have such constraints.

Low overhead. A desirable property of any monitor is that of imposing a minimal overhead on the observed system. PerfMC is a purely software monitor written in C. It uses a clean design in order to be as efficient as possible. The overhead on the network and on the monitored devices is extremely low; more details are given in Section V.

Batch operation. The monitoring system should be able to operate in batch mode, without any user interaction. At the same time, a suitable user interface should be provided. PerfMC uses a standard WEB interface to communicate with the user. This allows the user to inspect the monitor from remote locations using any WEB browser.

Generality. The monitor should be able to deal with a wide range of different networked devices, including network switches, tape libraries, uninterruptible power supplies and so on. The status of those devices needs to be monitored as well. We identified the Simple Network Management Protocol (SNMP) as a suitable candidate for the remote monitoring of a wide range of devices. More details on SNMP will be given in Section IV.1.

Easy Configuration. It was very important that the monitoring system could be configured easily using a standard, structured notation. We decided to use XML as the language in which the configuration file is written. There exist many tools able to generate, verify and transform XML documents. XML parsing can be done efficiently, as XML documents must obey strict syntactical rules (see W3 Consortium (2000) for details on XML). In Section IV.3 we will give more informations about the structure of the XML-based configuration file.

Reasonable scalability. The BaBar INFN reprocessing farm is expected to grow as the BaBar detector’s luminosity increases; this means that new machines will be added in the future to cope with the increased volume of data to be processed. The monitoring tool should be able to scale at least up to moderate cluster sizes (some hundred nodes).

## IV PerfMC ARCHITECTURE AND IMPLEMENTATION

PerfMC is a tool for medium-grained, continuous monitoring of computing clusters. It is written using the C language and currently operates under the Linux Operating System, but should be easily portable on any Unix flavor. PerfMC can monitor any networked equipment implementing an SNMP agent.

PerfMC is made of two threads: a SNMP collector and a WEB server, as depicted in Figure 1.

The SNMP collector thread periodically polls the monitored hosts using SNMP requests. The list of devices to monitor, and the list of SNMP variables to request for each one is contained in the XML configuration file. The collected data are stored on Round Robin Databases (see Oetiker (2003)) on the local disk; at the same time, an up-to-date “view” of the cluster is kept in memory. This view includes the status of each device (i.e., whether it is responding to SNMP queries), and the last value for each monitored variable. These informations are used by the web server thread to produce graphs and WEB pages. WEB pages are generated by applying user defined XSL Transformations to an internal XML representation of the cluster status. It is also possible to define an external postprocessor through which the generated pages (or graphs) are piped. The SNMP collector thread and the WEB server thread will be described in more detail in Sections IV.2 and IV.4 respectively. The format of the configuration file is described in Section IV.3.

### IV.1 The Simple Network Management Protocol

The SNMP architecture has three components:

* One or more Network Management Stations (NMS), which are responsible for monitoring and managing other devices;
* Network Nodes, which may be computing nodes or other equipments; each node hosts a software component called *SNMP-agent* which collects local data and answers requests coming from the NMS;
* A connectionless communication service; SNMP is usually implemented on top of UDP.

The SNMP agent running on each node manages a set of local *variables*; the management informations pertaining to a particular class of resources is defined in a Management Information Base (MIB). The NMS access the MIB by contacting the agent using SNMP primitives such as get or getnext to read values, and set to update values. Figure 2 illustrates the SNMP architecture.

The use of SNMP as a data collection protocol has some drawbacks. The protocol itself is very simple, and requires each NMS to periodically poll the other nodes. Polling introduces additional load on the network, due to the potentially large number of request/response packets. Also, SNMP agents have a simple structure and usually communicate only in response to get/getnext requests. They don’t perform management actions on their own, but require a NMS to take the decisions about what should be done.

However, SNMP has the advantage of being implemented in virtually every equipment having a network interface. This includes workstations, network switches/routers, tape libraries, printers and Uninterruptible Power Supplies. Many vendors only include SNMP agents on their devices, and no possibility is offered to use any custom code; SNMP is thus the only way to interact with these devices.

Recent versions of the SNMP protocol have additional features which solve some of the problems above. In particular, SNMPv2 implements *bulk requests*, which can significantly reduce the load on the network by packing several requests into a single datagram. In this way, the NMS access different MIB variables on the same host by using a single get request. The responses will be contained in a single packet as well. Our experience shows that the overhead put on the network by the SNMP request/response packets is very low. More details will be given in Section V.

PerfMC is being tested on a cluster composed of about 200 dual-processor Linux/Intel-based computing nodes, running the SNMP agent developed by the Net-SNMP (see Net-SNMP Project ).

### IV.2 The Collector

The SNMP collector thread is responsible for periodically polling the various monitored devices. For each device, the user can specify the list of MIB variables to observe and the frequency of the observations.

In order to improve the efficiency of the collector, multiple hosts are polled in parallel using non-blocking SNMP requests. The maximum number of hosts polled in parallel can be defined by the user, the only limitation being the number of simultaneous opened file descriptors supported by the underlying Operating System.

The SNMP collector stores the observations into a set of Round Robin Databases (RRD). A RRD can store time-series data (such as CPU utilization, network load, machine room temperature) in a compact way. Data must be entered into a RRD with a certain frequency. Old data are compacted by means of a consolidation function (any of Average, Minimum, Maximum and Last), or discarded. For example, the user may decide to store the average network utilization for the last week with one observation every 10 seconds, and for the last month with one observation every minute. The RRDTool package takes care of compacting observations older than a week by storing the average of six observations. Data older than one month are discarded. Round Robin Databases have constant size, which is determined when they are created. The RRD library provides the capability to plot the collected data in various ways. See Oetiker (2003) for a complete description of the RRD package.

The SNMP collector records the status of each machine in the cluster while receiving the observations. Such informations are kept in memory. Status informations include, for each host:

* the last observed value for every SNMP variable;
* whether the machine is responding to SNMP requests;
* the list of SNMP error messages generated by the machine.

In this way, the collector knows the operational status of each node in the cluster with a maximum delay equal to the time between consecutive polls. and can notify the system managers as soon as a problem arises.

### IV.3 Configuration File Format

The configuration file for PerfMC is written in XML. XML documents can be created by hand using a generic text editor, or using a specialized XML editor, or automatically generated by an application. The Gnome XML Library described in libxml  is used by PerfMC to parse, create and transform XML documents. The configuration file conforms to the structure declared in the monitor Document Type Declaration (DTD), reported in Figure 3. XML and DTD are described in detail in W3 Consortium (2000).

The `<monitor>` tag is the root element of the XML configuration file. The following optional attributes can be specified:

The maximum number of concurrent SNMP connections to use. If this value is set to $`N`$, then $`N`$ hosts are polled in parallel. Default: 50.

The log file containing messages generated by PerfMC. Default: do not write any log file.

The level of verbosity of PerfMC. Values range from 0 (high verbosity) to 3 (no verbosity). Default: 3 (no verbosity).

The directory containing the Round Robin Databases used to store the observations. Default: current directory.

The directory containing the XSL Transformation files. Default: current directory.

The directory containing the static HTML files served by the WEB server. Default: current directory.

The port on which the WEB server listens to requests. Default: 8001.

The file containing messages generated by the WEB server. Default: no log file.

The filter (postprocessor) applied to the pages generated by the WEB server. Default: no filter.

A space-separated list of file extensions; the filter defined by the http-filter attribute is applied to every file with a matching extension. Default: empty extensions list.

The configuration file consists of a sequence of `<host>...</host>` blocks, each one containing informations regarding a specific host to be monitored. The `<host>` tag has the following attributes:

The name of the device. Can be the host name, or any string used to uniquely identify it.

The IP address of the monitored device. If not given, the value of the name attribute is used as the host name of the device. Default: the value of the name attribute.

The delay in seconds between two consecutive observations of this device.

This attribute may contain any sequence of strings used to characterize this device. XSL transformations can perform selections based on the value of this attribute. Default: none.

The version of the SNMP protocol supported by the device. Recognized values are “1” for SNMPv1, and “2c” for SNMPv2. Default: “2c”.

For each device, the user must provide the list of SNMP variables to poll and the list of graphs which can be generated for that device, along with the RRD library commands used to produce the graphs. SNMP variables are described by `<mib>...</mib>` tags. The following attributes can be specified:

An identifier for the monitored quantity. The SNMP variables of a device needs to be uniquely identified by the value of this attribute. Variables belonging to different devices may have the same id.

The name of the variable in dotted-decimal (e.g. “.1.3.2.4.3.5.4”) or dotted-string (e.g. “system.sysUpTime”) notation.

Three different kind of SNMP variables are supported. “GAUGE” denotes a variable holding the absolute value of a quantity (e.g., a temperature). “DERIVE” denotes a variable holding a value over time (e.g., the number of bytes per second transmitted over a network interface). “COUNTER” denotes a variable having a value which never decreases (e.g., the total number of bytes send by a network interface). Default: “GAUGE”.

The name of the SNMP community to which the variable belongs. Default: “public”.

The minimum value the variable can assume. Default: none.

The maximum value the variable can assume. Default: none.

PerfMC creates a RRD for each monitored device. The layout of the RRD can be specified inside the `<archives>...</archives>` block. Each RRD is made of a number of Round Robin Archives (RRA), each described by a `<rra>...</rra>` tag. Details about the various options can be found in Oetiker (2003). The attributes are:

The consolidation function to use. Default: “AVERAGE”.

The XForm factor, i.e., the fraction of values which must be inserted into an interval. Default: 0.8;

The width in seconds of each interval.

The total length in seconds of the RRA. Data older than this value are discarded.

Finally, PerfMC is able to create graphs from the recorded data using the functions provided by the RRDTool package. Each graph is defined in a `<rrdgraph>...</rrdgraph>` block, and is characterized by the following attributes:

The unique identifier of the graph. Graphs from different devices can have the same id.

The width in pixels of the graph. Default: 400.

The height in pixels of the graph. Default: 180

The starting point of the graph expressed in seconds from the current time, or using the more readable notation supported by the RRDTool package. Default: “-3h” (three hours ago).

(none) The title of the graph.

The body of the graphs contain the list of instructions passed directly to the RRDTool library to produce the graph.

### IV.4 The WEB Server

PerfMC provides a WEB interface through an embedded HTTP server, implemented using the SWILL library Lampoudi and Beazley (2002). The WEB server has access to the in-core status informations about the cluster, which is kept up to date by the SNMP collector thread. Also, the WEB server has read-only access to the Round Robin Databases containing the historical data collected from the cluster. Using the graphing capabilities provided by the RRD library, the WEB server is able to dynamically generate plots from the data.

The WEB server can also produce an XML page containing the status of the whole cluster. As described in Section IV.2, the status informations include whether each host is responding to SNMP polls, the last received values of every polled MIB variable and a list of error messages reported by the host. The status document includes also the names of graphs which can be generated for each host in the cluster. An example of XML status document is reported in Figure 4.

The WEB server can apply to the XML status document XSL Transformations (XSLT) defined by the user. The resulting document is sent to the user. An XSL Transformation could produce, for example, an HTML page from the XML status file, so that WEB browsers can display it. The reader is referred to W3 Consortium (1999) for details about the XSLT language.

The user can define in the configuration file an optional postprocessor, through which every page generated by the WEB server will be filtered. The command to execute is specified with the http-filter attribute of the monitor element. The postprocessor can be any program which accepts input through its standard input and sends the result to standard output. The user can restrict the application of the postprocessor to pages with a certain extension (e.g., only pages with .html or .php extension). The list of file extensions can be specified with the http-filter-extensions attribute of the monitor element.

The WEB server recognizes the following types of Uniform Resource Identifiers (URI) specified in a HTTP GET request:

*/$`<`$hostname$`>`$/$`<`$graphname$`>`$\[?querystring\]* Returns the graph *$`<`$graphname$`>`$* for the machine whose identifier is *$`<`$hostname$`>`$*. It is possible to specify a query string for specifying some parameters of the generated graphs. Recognized query variables are: *width*, for the image width in pixels; *height*, for the image height in pixels; *start*, for the starting time of the data plotted on the graph, expressed in any notation recognized by the RRDTool. For those query variables which are not specified, the default value specified in the configuration file for that graph is assumed.

For example, the graph whose id is “cpu.png” for machine “localhost”, having width of 320 pixels, height of 200 pixels and showing data starting from 3 hours ago, can be requested with an URI like this: /localhost/cpu.png?width=320&height=200&start=-3h

*/status.html?applyTransform=$`<`$XSLT-file$`>`$* Applies the XSL Transformation specified in *$`<`$XSLT-file$`>`$* to the XML document describing the whole cluster status. The result is returned to the user;

*/$`<`$hostname$`>`$/status.html? applyTransform = $`<`$XSLTfile$`>`$* Applies the XSL Transformation specified in *$`<`$XSLT-file$`>`$* to the XML document describing the status of the single host *$`<`$hostname$`>`$*. Returns the result to the user.

Generating XML files by applying user-defined transformations to the status document is particularly useful. The user can customize the appearance of the generated HTML pages by writing a set of XSL Transformations. XSLT is a very powerful language and can be used to perform any transformations on the XML status document. Examples of HTML pages which can be produced are shown in the next section.

XML is a widely accepted standard, so it qualifies as a suitable mean for interchanging informations between PerfMC and other programs. XSL Transformations can be used to restructure and filter the data returned by PerfMC before passing them to an application.

## V A CASE STUDY

PerfMC is currently being used to monitor a Linux cluster used for High Energy Physics applications. The cluster is hosted at INFN Padova, Italy and is used to process very high volumes of data using CPU-intensive batch applications. It is made of about 200 dual processor, Linux/Intel machines, partitioned in about 150 clients and 50 servers. Both classes of machines use 1.26 GHz Pentium III processors and have 1GB of RAM. Client machines have two fast Ethernet controllers (at the moment just one Ethernet board is used). Server machines have gigabit Ethernet controllers and 1TB of local disk space as EIDE Raid arrays. The cluster is interconnected using a high-performance network switch, and is attached to a tape library.

For each machine in the cluster, the value of the quantities listed in Table 1 are monitored every 30 seconds. Note that the stock Net-SNMP agent under Linux does not provide informations on disk I/O activity. We extended the agent to report these informations by accessing the /proc/stat file.

Every quantity, with the exception of host ID, Name, Location and Uptime, is stored in a Round Robin Database. For each quantity, the average and maximum values are recorded in the RRD. Data for the previous week are kept with the granularity of one observation (average and maximum) every minute. Data for the previous month are kept with the granularity of one observation every hour. Data for the previous year are kept with the granularity of one observation every day. Finally, data more than one year old are discarded. The total size of each RRD is about 8 MB, and there is one RRD for each monitored host.

SNMP bulk get requests are used to get All the monitored variables for each machine. Network statistics collected with the tcpdump(8) utility show that the average size of an SNMP request is, in our case, about $`800`$ bytes, while the average size of an SNMP response is about $`1000`$ bytes. This gives an average network utilization of approximately 11KB/s due to the monitoring activity. The CPU overhead on the machines caused by the monitoring is negligible. The machine on which PerfMC runs is a dual Pentium III running at 1.26 GHz with 1GB or RAM. PerfMC has a very low (less than 5%) CPU utilization.

Figures 5 and 6 show two HTML pages generated by applying two different transformations to the same XML status document. Figure 5 shows a page containing part of the list of all machines in the BaBar INFN farm, with different colors indicating the CPU load of each machine. Figure 6 shows a more detailed view of a single machine, with the latest collected values of all the SNMP variables and some graphs.

## VI CONCLUSIONS AND FUTURE WORK

In this paper we have illustrated the architecture of a monitoring system for large computing clusters. A prototype written using the C language, has been implemented and is being used to monitor a cluster with 200 dual-processor Linux/Intel machines. The monitoring system uses asynchronous (non-blocking) parallel SNMP bulk requests to collect status informations from a wide variety of networked devices, and incorporates a WEB server which can generate graphs from the collected data. The WEB server can also produce an XML encoding of the current cluster status, to which an XSL Transformation can optionally be applied.

We believe that the most essential goals among those stated in Section III have been satisfied. The system is not intrusive in that all it needs is an SNMP agent running on each monitored device. At low polling rates (one observation every $`10`$ seconds) the overhead on the network and on the observed devices is very low. We do not recommend the use of PerfMC if higher, sub-second polling rates are required; in those cases a more specialized profiling system such as Supermon, should be preferable.

The scalability of PerfMC has been obtained by a clean design and an efficient implementation. Polling many hosts in parallel is a very trivial idea which indeed helped very much. The use of SNMPv2 bulk operations allowed to get the values of many variables from a single machine with just a request/response pair of packet, reducing the load on the network. Also, SNMP is a standard protocol and is implemented in virtually every device. The freely available implementation by the Net-SNMP project allows the user to extend the list of standard SNMP MIBs without the need to modify the agent’s code. In this way it is possible to monitor everything one could be interested in. We are currently using this feature to monitor the status and progress informations from the processes running on our cluster.

Finally, the heavy use of XML as the format for the configuration file, and that of XSLT to transform the status informations in arbitrary ways proved to be a good idea. XSLT transformations are used at the moment to produce a set of HTML pages showing in different ways the status of our cluster.

At the moment the prototype does not implement the alarm system. An alarm system is obviously needed to notify the system administrators in case of failures, so we are currently working on it. The alarm system will likely be implemented by listing in the configuration file a set of thresholds for each MIB variable. If a threshold is crossed in the specified direction, an alarm will be triggered. An additional alarm will be associated with each machine, and will be triggered if the machine does not reply to SNMP queries. Such a threshold-based alarm system is exactly the same implemented by the RMON protocol Stallings (1999), so alarms can be triggered directly by SNMP agents implementing RMON specifications.

The prototype implementation is performing well on our cluster, and no scalability limit has been encountered so far. However, it is obvious that a centralized monitoring system, even the most efficient one, cannot scale forever. In particular, we identified the updating of the Round Robin Databases as the most likely candidate bottleneck. As a first solution, we are currently trying to identify possible sources of inefficiencies in the RRDTool package. As a more long-term fix, we are considering the idea of partitioning the whole cluster among different monitors, each one observing a subset of the system. This would alleviate the scalability problem, as arbitrarily large clusters can be monitored by simply adding more monitors running on different machines. On top of these monitors, it is possible to build a hierarchy of *monitoring proxies* which will be used to fetch and consolidate the informations collected from the nodes behind them. The top (root) node will present a global view of the system to the user, or redirect user’s requests to the monitor responsible for observing the requested resource. Fault-tolerance can be implemented by means of standard techniques, such as electing a substitute when one of the monitors crashes. The user interface based on HTTP and XML was developed because it could also be used to exchange informations among monitors. The current PerfMC implementation can be extended to cope with a hierarchical monitoring infrastructure. What is required is the addition of a suitable client HTTP interface which can be used to contact the WEB server embedded in other monitors. In the same way it is possible to add other kinds of interfaces, such as a Lightweight Directory Access Protocol (LDAP) interface, or a text-only user interface.

###### Acknowledgements.

This work has been partially supported by the Istituto Nazionale di Fisica Nucleare, Padova. The authors thank Marco Cappellini who contributed many improvements to PerfMC, especially related to the user interface, and Roberto Stroili for providing comments and suggestions.