Volume 70 Received 3 May 2013  Statistical analysis of multipolemodelderived structural parameters and chargedensity properties from highresolution Xray diffraction experiments Radoslaw Kaminski,^{a}^{*}^{+}^{++} Slawomir Domagala,^{a}^{*}^{+++} Katarzyna N. Jarzembska,^{a} Anna A. Hoser,^{a} W. Fabiola SanjuanSzklarz,^{a} Matthias J. Gutmann,^{b} Anna Makal,^{a} Maura Malinska,^{a} Joanna M. Bak^{a} and Krzysztof Wozniak^{a} ^{a}Department of Chemistry, University of Warsaw, Pasteura 1, 02093 Warszawa, Poland, and ^{b}ISIS Neutron and Muon Source, Science and Technology Facilities Council, Rutherford Appleton Laboratory, Harwell Oxford, Didcot, Oxfordshire OX11 0QX, UK A comprehensive analysis of various properties derived from multiple highresolution Xray diffraction experiments is reported. A total of 13 chargedensityquality data sets of oxalic acid dihydrate (C_{2}H_{2}O_{4}·2H_{2}O) were subject to HansenCoppensbased modelling of electron density. The obtained parameters and properties were then statistically analysed yielding a clear picture of their variability across the different measurements. Additionally, a computational approach (CRYSTAL and PIXEL programs) was utilized to support and examine the experimental findings. The aim of the study was to show the real accuracy and interpretation limits of the chargedensityderived data. An investigation of raw intensities showed that most of the reflections (6070%) fulfil the normality test and the lowest ratio is observed for weak reflections. It appeared that unitcell parameters are determined to the order of 10^{3} Å (for cell edges) and 10^{2} ° (for angles), and compare well with the older studies of the same compound and with the new 100 K neutron diffraction data set. Fit discrepancy factors are determined within a 0.5% range, while the residual density extrema are about ±0.16 (3) e Å^{3}. The geometry is very well reproducible between different data sets. Regarding the multipole model, the largest errors are present on the valence shell chargetransfer parameters. In addition, symmetry restrictions of multipolar parameters, with respect to local coordinate systems, are well preserved. Standard deviations for electron density are lowest at bond critical points, being especially small for the hydrogenbonded contacts. The same is true for kinetic and potential energy densities. This is also the case for the electrostatic potential distribution, which is statistically most significant in the hydrogenbonded regions. Standard deviations for the integrated atomic charges are equal to about 0.1 e. Dipole moments for the water molecule are comparable with the ones presented in various earlier studies. The electrostatic energies should be treated rather qualitatively. However, they are quite well correlated with the PIXEL results. Keywords: multipole model; charge density; statistical analysis; oxalic acid. 
The experimental chargedensity field has undoubtedly reached its mature stage (Coppens, 2005). Having crystalline samples of an appropriate quality and modern Xray equipment, it is possible to collect very good highresolution diffraction data sets. Such data are indispensable to quantitatively model electrondensity distribution within crystals. Nowadays, because of the advanced development of instruments, techniques and software, chargedensity studies of simple organic molecules have become almost routine, similar to structural analysis. Nevertheless, there are still a number of challenging scientific problems in the field, such as chargedensity distribution modelling of crystal structures containing heavy atoms (Zhurov, Zhurova & Pinkerton, 2011; Zhurov, Zhurova, Stash et al., 2011b; Schmøkel et al., 2013), crystals exhibiting a high degree of ionicity (e.g. minerals) (Schmøkel et al., 2012), treatment of disorder (Bak et al., 2009; Meindl et al., 2009), description of inner atomic electronic shells (Fischer et al., 2011), proper treatment of hydrogen atoms (Hoser et al., 2009; Munshi et al., 2008) and thermal motion deconvolution (Zhurov, Zhurova, Stash et al., 2011a), among others. More complex crystal cases may exhibit a variety of interesting and useful properties (mechanical, optical etc.) and thus are worth exploring. Extraction of reliable electrondensity information for such samples often causes severe difficulties. Consequently, new approaches are constantly being developed, as for example the most recent Xray constrained wavefunction fitting technique (Jayatilaka & Grimwood, 2001). There are also various ideas for improving the most widely used HansenCoppens model (Hansen & Coppens, 1978). Modelling procedures should, however, be designed with care, as they could also become a source of significant bias affecting the final results (Hoser et al., 2009; Spackman & Byrom, 1997a; El Haouzi et al., 1996; Howard et al., 1995). Therefore, we do need to be aware of the limitations of the applied tools and methods in the widest context possible.
Recently, a number of studies explored the deficiencies of the multipole model (Bak et al., 2009, 2011, 2012; PoulainPaul et al., 2012; Dittrich et al., 2012). Nevertheless, it should be noted that recorded diffractedbeam intensities also constitute a source of experimental errors, being further propagated across all stages of electrondensity elucidation. Thus, to obtain an overall view of the method's real shortcomings all of these factors should be taken into account. This could be done by applying statistical methods to an experimental sample. Quite recently, Grabowsky et al. (2009) have investigated statistical aspects of the transferability principles regarding chargedensity distribution properties. Additionally, the authors briefly compared their results with the literature examples of standard deviations obtained from several (from two to four) chargedensity measurements of the same compound. However, the mentioned analyses are not fully reliable because of the few data sets used in the study. Hence, in this paper, we indicate the precision of the determined structural and chargedensity parameters on the basis of multiple highresolution data collections of oxalic acid. The main idea of our study can be summarized by the following question: `What are the differences between the data sets, geometrical parameters and resulting chargedensity models of the same chemical system derived from different highresolution Xray measurements which we could classify as acceptable?'
We have collected a total of 13 chargedensityquality data sets for oxalic acid dihydrate (C_{2}H_{2}O_{4}·2H_{2}O) crystals. This constitutes a good statistical sample, especially considering the very timeconsuming experiments. There are several reasons for the choice of oxalic acid. Firstly, the crystals of this compound can be easily grown and, importantly, with good reproducibility of their quality. The other reason is that oxalic acid is probably the best chargedensityinvestigated compound so far. Since the very beginning, the crystals of this compound were used as a reference. The earliest data collections were conducted by Stevens et al. (1979), followed by Coppens et al. (1984) who carried out an excellent survey of both Xray and neutron data sets collected across the globe (four Xray and five neutron data sets are reported therein). The most uptodate chargedensity study of the kind was performed by Martin & Pinkerton (1998). This was also the first application of the CCD detector for experimental electrondensity elucidation. Our investigation is supplemented by brief comparison with these older results. This allows us to gain a deeper insight into the experimental reproducibility and various factors that influence the chargedensity Xray data quality.
The studied oxalic acid dihydrate was purchased from the L.P.PH `OCH' Lublin chemical company as a crystalline powder. The compound quickly dissolves in water and, upon evaporation at room temperature under a microscope, forms crystals suitable for Xray examination. Crystals of the size appropriate for a neutron diffraction experiment were grown by a very slow evaporation of water at room temperature from a mediumsize beaker.
A singlecrystal neutron diffraction experiment was performed at the SXD beamline (Keen et al., 2006), at the neutron spallation source in the ISIS facility (Oxfordshire, UK). For the purpose of the experiment a single crystal of dimensions 4 × 2 × 1 mm was attached to an Al pin with an adhesive Al tape and cooled to 100 K in a closedcycle refrigerator with He as an exchange gas. A total of six crystal orientations were recorded with an exposure time of 1 h per orientation. Data reduction was accomplished using the locally available SXD program (Keen et al., 2006). The crystal structure (with anisotropic hydrogenatom thermal parameters) was then refined in the JANA package (Petrícek et al., 2006) using the timeofflight data. It is also worth noting that the severe extinction was successfully modelled by the BeckerCoppens approach (Becker & Coppens, 1974a,b, 1975). A summary of refinement information follows: total No. of reflections: 2630; No. of reflections with : 1583; index ranges: , , ; No. of parameters: 65; discrepancy factors: R[F] [] = 6.71%, wR[F^{2}] [] = 6.45%, R[F] (all data) = 10.83%, wR[F^{2}] (all data) = 6.83%, goodness of fit (GoF) [] = 1.64; nuclear density largest extrema: 4.75/+4.40 fm Å^{3}; extinction parameter (Gaussian isotropic): 2.36 (6). The CIF file is available in the supplementary material,^{1} or can be retrieved from the Cambridge Structural Database (Allen, 2002) (deposition number: CCDC 930687). Hereafter the model is abbreviated as sxd1.
Highresolution singlecrystal Xray diffraction experiments were performed on three different experimental setups. Details are given in Table 1. In each case a new crystal was grown and mounted on the goniometer head, installed on the fourcircle goniometer, and cooled to a temperature of 100 K (temperature instability is an additional factor affecting derived chargedensity parameters). Depending on the crystal overall quality and its orientation, the data sets were collected using different strategies optimized for each case (only scans; scan width 0.5°; exposure time from 3 to 10 s on low angles, 45 to 70 s on high angles). Data sets were integrated with the respective diffractometer software: APEX2 (Bruker AXS Inc., 2013) or CRYSALIS (Agilent Technologies, 2010). Further data reduction, correction application and merging were done with the SORTAV program (Blessing, 1987, 1989, 1995a).

All multipole refinements, within the HansenCoppens formalism, were performed in the MOPRO suite (Jelsch et al., 2005) combined with the new version of the University at Buffalo Data Bank (UBDB) (Jarzembska & Dominiak, 2012) and neutron diffraction data. All deformation parameters are defined with respect to their local Cartesian coordinate systems, which were automatically assigned by the LSDB program (Jarzembska & Dominiak, 2012; Volkov et al., 2004). All refinements were based on F and only the reflections fulfilling the condition were taken into account (see below). In all cases statistical weights were used (i.e., for the ith reflection ). Initial atomic coordinates (x,y,z) and anisotropic displacement parameters (U_{ij}'s) for each atom were taken from the spherical refinement stage, whereas initial multipolar and contractionexpansion parameters were transferred from the UBDB with the aid of the LSDB program. The charges of both moieties (i.e., acid and water molecules) were initially set neutral. The multipole expansion was truncated at the hexadecapole (l_{max} = 4) and quadrupole (l_{max} = 2) levels for all nonhydrogen and hydrogen atoms, respectively. Each atom was assigned with core and sphericalvalence scattering factors derived from Su and Coppens atomic wavefunctions (Su & Coppens, 1998). The refinement of each oxalic acid model was performed in two stages, (a) and (b), which are summarized below.
(a) During the first stage of the refinement procedure the SHADEserverderived temperature factors for hydrogen atoms were used (Madsen, 2006). The oxalic acid and water molecules were entered into the SHADE server separately, so as to estimate the atomic displacement parameters (ADPs) using solely the riding approximation approach [a reliable TLS fit (Schomaker & Trueblood, 1968; Cruickshank, 1956; Sands, 1982) is impossible for such small molecules]. At this initial refinement stage the local symmetry restrictions were strictly followed. Such a procedure was found to be essential as it prevented falling into the false minimum during the leastsquare refinement procedure [the false minimum problem is discussed e.g. by Watkin (2008) or in the JANA program `cookbook' (Petrícek et al., 2006)]. XH bond lengths (X = nonhydrogen atom) were standardized to neutronnormalized distances taken from the literature and restrained using . This approach has recently been successfully tested on a variety of uracil derivatives, providing geometrical results comparable with the corresponding theoretical periodic computation and neutron studies, and seems to be the best choice in the absence of neutron diffraction data (Hoser et al., 2009; Jarzembska et al., 2012). Additionally, hydrogenatom U_{iso} parameters (i.e., isotropic thermal parameters) were restrained to the value of 1.5 ×U_{eq}^{O} using (where an appropriate restraint weight is equal to ). The refinement was carried out according to the following scheme: (i) scale factor (which was also refined in all other stages); (ii) atomic coordinates; (iii) atomic coordinates and ADPs; (iv) a few cycles of SHADE estimation and point (iii) until convergence is reached; (v) multipole parameters in an incremental manner; (vi) all population, positional and thermal parameters simultaneously; (vii) a few cycles of SHADE estimation and point (vi); (viii) 50 cycles of block refinement of (vi) and parameters for nonhydrogen atoms; (ix) all parameters simultaneously; (x) a few cycles of SHADE estimation and point (ix).
(b) In the second stage of the refinement procedure the neutron diffraction data were included. XH bond lengths are restrained to the values derived from the sxd1 model. The hydrogenatom ADPs from neutron studies are scaled using the Blessing method (Blessing, 1995b) and subjected as fixed to the multipolar refinement performed as follows: (xi) refinement of all population, positional and thermal parameters simultaneously; (xii) 50 cycles of block refinement of point (xi) and parameters for nonhydrogen atoms; (xiii) all possible parameters simultaneously. Beginning from point (xi) the local symmetry is not preserved and all populations are allowed to vary.
The final parameters characterizing refinements for all data sets are summarized in Table 1. CIF files for each refinement are present in the supplementary material, or can be retrieved from the Cambridge Structural Database (deposition numbers: CCDC 925544925554, CCDC 925558, CCDC 925559). All identical refinements were performed smoothly and reached an acceptable convergence of maximal shift over the estimated standard deviation equal to around 10^{3} or less (except for the oxa7 model, where the value of 10^{2} was obtained, though all parameters are reasonable). This ensures that all parameters derived from multipole models are comparable one to another in the sense of statistical analysis of various factors treated as a sample. In other words, the differences between models are not `hidden' by the deficiencies of the multipole model itself. It is also worth mentioning that in all cases the Hirshfeld rigidbond test (Hirshfeld, 1976) is fulfilled, ensuring a proper deconvolution of electron density and thermal motion features.
Finally, the use of the cutoff during multipole refinement needs a brief discussion. This issue was described some time ago by Watkin (2008) in his excellent review paper; therefore here we shall focus solely on the cutoff effect related directly to the subject of the current contribution. Our definition of `weak data' is that the ratio is below 3.0. Such a cutoff criterion is based on statistical considerations related to the normal distribution (the normality of raw intensity data is discussed later on in this article). However, in general, inclusion or omission of weak data in the refinement may bias the results. In fact, the correct procedure whether to reject or include such reflections should depend on their reciprocal vector length. If reflections are located in the lowangle region they should be included as they are related to the atomic valence populations. Rejection of such reflections obviously may bias to some extent the result of the chargedensity refinement. When a reflection belongs to the highangle shell, in our opinion, it is a matter of choice. Such a reflection shall mostly affect the positional and thermal motion parameters. Following the findings of Hirshfeld & Rabinovich (1973) and Seiler et al. (1984) regarding the data cutoff, it can be concluded that in general differences on the refined parameters should be very small. Further considerations on a resolutiondependent reflection rejection were presented by Schwarzenbach et al. (1989).
Keeping in mind all the above factors, we collected our data sets in a careful way so as to obtain a minimum number of weak lowangle reflections. In fact they do not exceed 20% (out of several hundreds of recorded reflections). This number is further reduced during the outlier rejection stage in SORTAV. To ensure that our refinement results are not substantially biased by the applied cutoff, we performed all chargedensity refinements both with and without application of this condition. The results of the statistical analysis of multipole populations are presented in the supplementary material. In general, there are no statistically significant differences (being within 1 margin) in terms of the multipole populations between the two sets of results.
The experimental chargedensity distributions obtained by the multipolar approach were analysed by means of Bader's Quantum Theory of Atoms in Molecules (QTAIM) (Bader, 1994). The VMOPRO module (part of the MOPRO package) was used for Fourier syntheses and bond critical point (BCP) evaluation. The same program was used for Fourier map visualization and computation of dipole and quadrupole moments. Visualization of bond paths and surfaces was accomplished with the MOPROVIEWER program (Guillot, 2012). Neighbouring molecules, used for the bond path search, were built up using the CLUSTERGEN program (Kaminski et al., 2013). Electric dipoles and quadrupoles were derived from the multipole populations using the formulas given by Spackman (1992) and Coppens (1997). First and second moments of the charge distribution were computed with respect to the centre of mass of the considered molecule. The traceless convention was used for the quadrupole moment and eigenvalues (Q_{xx}, Q_{yy}, Q_{zz}) were determined. Evaluation and integration of atomic basins were performed with the WINXPRO program (Stash & Tsirelson, 2002, 2005).
The experimental electrostatic interaction energy values (E_{elst}^{exp}) between various bimolecular complexes composed of oxalic acid or water molecules were computed with the VMOPRO module as an integration over the whole space of the electron density of molecule A multiplied by the electrostatic potential generated by molecule B (or reciprocally):
Such integrals were computed using a numerical integration method based on a spherical grid around selected atoms. The GaussChebyshev (Becke, 1988a) and Lebedev & Laikov quadratures (Lebedev & Laikov, 1999) were used for the radial and angular parts, respectively. Radial coordinates and weights were remapped using the formula of Treutler & Ahlrichs (1995). The integrations involved 100 radial and 434 angular quadrature points. Interaction energies were calculated between pairs of neighbouring molecules in contact in the crystal for which two atoms were separated by a distance smaller than, or equal to, the sum of their van der Waals radii.
Lattice and bimolecular motif energy computation was performed using the OPIX program (Gavezzotti, 2002, 2003). All the 13 chargedensityquality structures of oxalic acid were used to calculate the corresponding molecular electrondensity grids by standard quantum chemical methods using GAUSSIAN (Frisch et al., 2009) at the MP2/631G** level of theory (Møller & Plesset, 1934; Krishnan et al., 1980). The electrondensity model of each molecule was analysed in the OPIX package. It was possible to extract the interaction energy values, and their breakdown into electrostatic, polarization, dispersion and repulsion components, for bimolecular motifs present in the crystal structure. The electrostatic and polarization energy terms were compared directly with the results obtained on the basis of the multipolar model.
Cohesive energy computations were performed with the CRYSTAL program package (Dovesi et al., 2005, 2009) at the DFT(B3LYP) level of theory (Lee et al., 1988; Becke, 1988b). Two types of molecular allelectron basis sets were used, i.e. either 631G** (Krishnan et al., 1980) or pVTZ (Dunning, 1989), which is specified in the text. Both the Grimme dispersion correction (Grimme, 2004, 2006) and a correction for basis set superposition error (BSSE) were applied (Civalleri et al., 2008). Ghost atoms were selected up to 5 Å distance from the considered molecule in a crystal lattice and were used for the BSSE estimation. The evaluation of Coulomb and exchange series was controlled by five thresholds, set arbitrarily to the values of 10^{7}, 10^{7}, 10^{7}, 10^{7}, 10^{25}. The shrinking factor was equal to 8, which refers to 170 kpoints in the irreducible Brillouin zone and assures the full convergence of the total energy. The cohesive energy (E_{coh}) was calculated as described below (Civalleri et al., 2008):
where E_{bulk} is the total energy of a system (calculated per unit cell), E_{mol} is the energy of a molecule extracted from the bulk and Z stands for the number of molecules in the unit cell. For comparative purposes two geometry optimizations at the DFT(B3LYP) level of theory with the pVTZ and 631G** basis sets were carried out. During the optimization procedure atomic coordinates were varied while cell parameters were kept fixed.
In this section the most common and useful statistical descriptors are summarized (Brandt, 1998). In the case of a set of measurements { x_{i} } () a sample mean is defined as
which stands for the estimator of the expectation value. The sample standard deviation is expressed as
It is worth noting here that, whereas s_{x}^{2} constitutes the unbiased estimator of the sample variance, s_{x} itself in general cannot be treated as the unbiased estimator of the standard deviation. Under the assumption that the sample obeys the normal distribution it is, though, possible to obtain the unbiased estimator by multiplying s_{x} by a certain correction factor [which can be derived by the use of Cochran's theorem (Cochran, 1934)]. However, this approach was not utilized here as the correction factor is very close to unity ( 0.979 for 13 data sets) so its application does not affect the final conclusions.
Numerous advanced statistical methods concern the socalled normal distribution, which means that a given statistical population can be described by a Gaussiantype function. However, to apply the appropriate statistical tests, one needs to check whether the examined statistical sample can be treated as normally distributed. The well known ShapiroWilk test (Shapiro & Wilk, 1965) is the best method for such verification. In this approach one may reject the null hypothesis saying that a given property comes from a normally distributed population, or state that there is not enough information to reject it. The test algorithm is as follows: (i) sort the values in an increasing order; (ii) compute the test statistic as
where K equals n/2 or (n  1 )/2 (for even or oddpopulated samples, respectively), whereas the a_{k,n} values are tabulated in standard statistical tables; (iii) the null hypothesis is rejected when , where the critical values are taken from statistical tables (where is a given significance level). In our case of 13 data sets, the critical values equal 0.814, 0.866 and 0.889 for equal to 1, 5 and 10%, respectively.
In the beginning of the analysis and discussion it should be noted that we were not able to perform a fully systematic study of various factors (crystal, Xray radiation source, detector etc.) affecting the final set of obtained chargedensity quantities. All 13 data sets were measured within several years mainly for the purpose of diffractometer alignment for chargedensity studies. Therefore, we limited ourselves to investigating the overall variability of different experimentally derived parameters. Nevertheless, the study can be extended in the future when more data sets, preferably measured using other Xray sources, are available.
The abovedescribed statistical methods were first applied to analyse, compare and contrast a set of unitcell parameters determined from different measurements. The corresponding numerical values are given in Table 2. It is worth noting here that the unitcell setting used in this study (P2_{1}/c) is different to that used in previous works (P2_{1}/n). The unitcell setting was changed, following Dunitz (1996) and Feast et al. (2009), to minimize parameter correlations in obliqueangle crystal systems.

It is a well known fact that the precision of CCD detectors is much worse than that obtained with point detectors (Herbstein, 2000). Our results for estimated standard deviations are, though, comparable with the early estimates given in the literature (Paciorek et al., 1999). They are even larger than the ones reported by Guzei et al. (2008) for 2dimethylsulfuranylidene1,3indanedione, where the sample standard deviation value (averaged over 223 data sets) amounts to about 0.0005 Å. It is also clearly seen that the calculated sample standard deviations are at least an order of magnitude larger than the values resulting from the leastsquare fits of orientation matrices. In addition, they are rather large, especially in the case of the unitcell angle.
It is noticeable that the obtained standard deviation values are quite similar to the corresponding parameters resulting from the orientation matrix leastsquares fit to the neutron diffraction data. It should be stressed that neutron data are processed quite differently when compared to the Xray diffraction case. The sxd1 model unitcell parameters presented in Table 2 are computed as averaged values over all detector settings using the centreofmass position of peak centroids on the detector pixels and along timeofflight as determined by the initial peak searching routine. The lattice parameters denoted as sxd2 in Table 2 use the fitted position of the peaks as determined during the integration. They are very close to the sample means computed from the Xray diffraction experiments. We found out that for the neutron data the influence of the unitcell parameters (sxd1 or sxd2) on various structural parameters (e.g. hydrogenatom ADPs) is rather negligible. These results indicate a very good agreement of the two methods, keeping in mind that the neutron diffraction technique is also significantly affected by factors such as extinction, absorption, neutronbeam wavelength distribution or flux inaccuracies.
In order to analyse the nature of the experimental errors, the ShapiroWilk test was applied. This yielded the conclusion that for all of the unitcell parameters the normal distribution hypothesis cannot be rejected at any of the considered significance levels. It is then justifiable to treat the errors as a result of random fluctuations of environmental conditions during each experiment (e.g. slightly unstable temperature, goniometer angle inaccuracy, detector distance etc.). It is worth stressing that the unitcell constants used in the study were obtained for precise and accurate chargedensityquality measurements. The highangle data allow a very precise description of the atomic positions and thus also of the unitcell parameters. Consequently, in the case of many standard lowresolution structural data, the experimental standard deviations may be even larger. Importantly, the socalled unconstrained unitcell refinements (i.e., with the and angles not being fixed at 90°) led to values of right angles in the monoclinic system not significantly different from 90°. Finally, comparison of the above results with the ones presented by Coppens et al. (1984) and Martin & Pinkerton (1998) (see Table 2 in the original paper) shows a very good agreement. The unitcell parameters in the previous works deviate by about 4 at worst (for the a parameter) from our values. This confirms a good reproducibility of unitcell parameters, even recorded with the older experimental setups.
On the whole, in the case of a single Xray diffraction experiment we suggest using the leastsquares estimates of unitcell parameters, even though they are clearly lower than the real ones [this conclusion matches the statement made by Farrugia & Evans (2005) and Guzei et al. (2008)]. We believe it is better to show the results that are at least mathematically sound than the ones that are based on rather vague reasoning. In addition, our results for unitcell parameters have serious implications for various analyses of structural phase transitions or timeresolved processes. In such cases particular attention has to be paid to unitcell changes.
Careful analysis of the raw intensity data is a desirable procedure in every chargedensity study (Kuntzinger et al., 1999). For the purpose of this study, we decided to provide a statistical analysis of the collected raw data. Thus, in order to apply the statistical methods to raw intensities we first brought them to the common (absolute) scale by applying the respective scale factors from the final multipole refinements. Subsequently, we extracted reflection indices being common for all 13 data sets. This yielded us a set of 1662 reflections with resolution up to about 1.0 Å^{1} (Table 1).
At first, we analysed raw intensities (I) and their estimated standard deviations [] across all data sets. Typical scatter plots are shown in Fig. 1, where the second data set is compared to the 12th one (all other comparisons are very much the same). For all data sets the I versus I relation is almost exactly linear (with the linear determination coefficient, R^{2}, being in the range from 0.918 to 0.999). In the case of standard deviations discrepancies from linearity are more pronounced (R^{2} = 0.3740.847); however, these are less representative as they are only estimated during the integration procedures.
 Figure 1 Example plots of (a) I_{set 2} versus I_{set 12} (R^{2} = 0.997) and (b) versus (R^{2} = 0.831) relationships for the raw (scaled) data sets. The black solid line is the reference line. 
The most important issue is the normality of such distributions. To answer this question we utilized the ShapiroWilk test based on the raw intensities. The resultant plots are shown in Fig. 2. A plot of the intensitybased ShapiroWilk test with respect to the resolution shows that at least more than 50% of reflections pass the normality test. The lowest acceptance ratio is observed for the most highresolution reflections, and for the rest this ratio varies from 60 to 70%. Interesting features are also present in the plot illustrating the results of the mentioned test with respect to the intensity values. Clearly, a maximum is observed. The largest deviations from normality are observed for the very weak reflections (25% pass the test), then the ratio of acceptance increases to 80%, but then it drops down to about 62% for the most intense reflections. The low acceptance ratio for the very weak reflections can be understood as resulting from the increased measurement problems. In the case of strong reflections the background estimation and detector characteristics could play an important role.
 Figure 2 Histograms showing the acceptance ratios ( = 10%) for the ShapiroWilk test on the reflection raw intensity values versus (a) resolution and (b) raw (scaled and averaged) intensity. 
The very first quantities derived for all the chargedensity models are the R factors and residual density distributions. About 40 years ago, Abrahams et al. (1970) investigated a set of 17 independently measured structures of D(+)tartaric acid. They concluded that the average Rfactor error, with respect to different measurements, was about 5.8%. In our case, because of the great advances made in the development of datacollection techniques and devices, it is much improved (Table 3). The standard deviations reach values of about 0.5% [i.e., an order of magnitude better than in the Abrahams et al. (1970) study]. Similar results are observed for the goodnessoffit (GoF) parameter, where the mean value is equal to 1.2 (2). In the case of two data sets out of 13, the GoF is lower than 1.0, which may suggest slight overfitting of the model. However, on the basis of the above statistics, such values of GoF are still within the 3 significance level.

Some important information may also be extracted from the extremal residual density values computed for the whole unit cell. The average minimum value equals 0.17 (3) e Å^{3}, whereas the average maximum value amounts to +0.15 (3) e Å^{3}. This suggests that the overall mean residual density is very flat. The flatness criterion seems to be best fulfilled by the oxa12 model ( = 0.11/+0.12 e Å^{3}). The largest extreme values are observed for oxa6 ( = 0.22/+0.20 e Å^{3}); however, they are still commonly acceptable.
Additional insight into the residual density distribution is achieved by analysing the twodimensional plots in the oxalic acid plane (Fig. 3). The map of the mean residual density values is very flat and almost featureless. On the other hand, the corresponding standard deviation map shows some nonuniform features, located mostly in the vicinity of the molecular fragments. It has to be noted that while summing up all maps some systematic errors start to appear in the regions not well modelled by the atomic multipole expansion (i.e., far from atomic positions) (see supplementary material).
 Figure 3 Twodimensional maps of (a) mean residual density distribution ( = 0.10/+0.06 e Å^{3}) and (b) computed residual density standard deviation ( = 0.07 e Å^{3}) [C1O2C1_{(1x, 1y, 1z)} plane; contours at 0.05 e Å^{3}, blue lines  positive values, red  negative]. 
Finally, the socalled residual density analysis introduced by Meindl & Henn (2008) was performed and the results are summarized in the supplementary material. Most of the plots generated for the analysed models indicate Gaussian noise for the corresponding residual density maps. Only in a few cases are some asymmetric features observed, suggesting some contribution of systematic errors. These are, though, still acceptable by means of the common criteria applied in the experimental chargedensity studies.
The discussion would be incomplete without a brief mention of the static deformation density maps. Such maps constitute a useful tool to analyse the quality of the multipole model fit, and in some cases they also provide additional chemical information. The twodimensional maps of the mean static deformation density and its standard deviation are presented in Fig. 4. There is a good reproducibility of the qualitative features in the density accumulation (i.e., lone pairs and bonding density) and reduction regions. However, in the standard deviation map it is clearly seen that the largest errors are present in the vicinity of atomic positions. At the nuclei they reach values up to about 1.7 e Å^{3}. This results from the well known problem regarding the multipole model inflexibility and hence its inability to describe the nucleus region properly (Abramov et al., 2000; Spackman & Byrom, 1996; Madsen et al., 2004; Volkov et al., 2001). On the other hand, much smaller maxima are observed in the chemically important areas such as the oxygenatom lone electron pairs, suggesting that in such cases quantitative information can be obtained and binding conclusions drawn. Importantly, the bonding regions (e.g. hydrogenbond interaction lines) exhibit error values smaller than 0.05 e Å^{3}.
 Figure 4 Twodimensional maps of (a) mean static deformation density distribution ( = 0.39/+1.10 e Å^{3}) and (b) its computed standard deviation ( = 1.73 e Å^{3}) [C1O2C1_{(1x, 1y, 1z)} plane; contours at 0.05 e Å^{3}, blue lines  positive values, red  negative]. The largest maxima in the standard deviation map are located at the nuclei positions. 
It is also interesting and valuable to point out the ShapiroWilk test performed at each point of the grid. Such a 01 grid (1 means that the test is not fulfilled, 0 otherwise) is then multiplied by the original mean grid (Figs. 3a and 4a, respectively). Such plots are presented in Fig. 5. Clearly, the normality test is well fulfilled for the residual density in the plane of the molecule. This is not the case for the deformation density. Plot contours indicate that the C1O2 bond density and lone pair regions are those most affected by systematic errors. This agrees well with the already mentioned discussion related to Fig. 4b. However, the nuclei positions are not much affected. Such results may be related to the different source of errors affecting the bonding and nuclei density regions.
 Figure 5 Twodimensional maps of the ShapiroWilk test performed for (a) residual density and (b) deformation density maps [C1O2C1_{(1x, 1y, 1z)} plane; contours at 0.05 e Å^{3}, blue lines  positive values, red  negative]. (c) Twodimensional map of the mean residual density after application of the Student's ttest. 
To test the significance of the residual density peaks on the residual map we used the Student's ttest, with an assumption that the normal distribution is fulfilled for the residual density. The onesample Student's ttest, with a null hypothesis that the residual density should be equal to zero, was performed on every point of the residual density map (Fig. 5c). As expected, most of the peaks were removed and the map is clearer. However, even at the significance level of , some residual density peaks are still present. Therefore, we can conclude that, at these peak positions, the residual density is significantly different than zero, which suggests that there are some systematic errors in the derived chargedensity models.
Structural parameters which can be readily analysed are bond lengths and bond angles. The corresponding mean values and the standard deviations characterizing the investigated sample are presented in Table 4. Interestingly, the calculated standard deviations are usually very similar to the ones computed by the leastsquare procedure after the last refinement cycle. This shows a rather narrow distribution of these parameters and, in consequence, a very good reproducibility of geometrical results. For covalently bound atoms, variations of bond lengths are of the order of 10^{4}, whereas for angles they are of the order of 10^{2}. CC and CO bond lengths agree with the neutron diffraction results within 3 (taking into account neutron data inaccuracies as well). The calculated mean values seem to be slightly smaller, i.e. of some parts of 10^{3}, than the sxd1 model values. Nevertheless, the agreement is quite good. The obtained results imply that the leastsquares estimates of bond lengths and bond angles may be considered as good estimators of the structure determination precision. It may also be concluded that for the standard lowresolution data of good quality such estimates would be an order of magnitude larger. It should be recalled here that the angles computed for bonds with hydrogen atoms are determined with lower precision than those for heavyatom bonds. Here, the H2O2H3 and all OHO angles show variations of several degrees (i.e., the standard deviation ranges from 2 to 7°). This is rather discouraging in terms of analysis of the very fine details of hydrogen bonding. However, the multipole model still provides a significant improvement when compared to the independentatom model (IAM) results (Hoser et al., 2009). In addition, such bondangle values and their sample standard deviations compare well with the neutron diffraction results within the 3 limit. The same is true for a comparison of our geometrical results with the ones obtained by Martin & Pinkerton (1998) (the 3 limit is generally fulfilled). It might be also estimated that the bondlength precision error in the case of the old data is roughly an order of magnitude larger than that resulting from our experiments (i.e., the positional parameter discrepancies are up to 0.002 Å, which may indicate that the bond standard deviation is even larger).

In order to compare thermal motion parameters we computed the rootmeansquare displacement factors along the principal axes for all nonhydrogen atoms. These are obtained from U_{ij}'s by well known transformations (Waser, 1955; Busing & Levy, 1958; Sands, 1982). Respective values for all models are summarized in Table 5. Such an approach enables comparison between the probability ellipsoid `sizes'. Almost all principal components exhibit rootmeansquare displacement sample deviations of around 0.002 Å (corresponding to 4 × 10^{6} Å of mean displacement) in each direction. The relative error when estimating the physical atomic displacements is then about 2% or less which is due to the thermal motion. This is a very encouraging result regarding the analysis of temperature dependence of various structural parameters (e.g. during the phase transition). For instance, in the case of the Lauemethodbased photocrystallographic techniques, the observed structural changes are much larger (of the order of 10^{1} Å) (Benedict et al., 2011; Makal et al., 2011, 2012) and, thus, in view of the above results, they cannot be associated solely with the machine precision (i.e., they seem to be a real physical effect). Similarly, well described ADPs would constitute a firm basis for reliable solidstate entropy determination from the Xray diffraction data (Madsen & Larsen, 2007). In addition, the ShapiroWilk test is well fulfilled by all principal components of oxygen atoms, whereas for the carbon atom the normality is questionable. This is also reflected in the overall ellipsoid orientation variability. The estimated deviations from the average value do not exceed 2° (Table 5), except for the C1 atom where the deviation reaches about 4°. It is worth noting that all these values are not considered significantly different from zero, which suggests that there is almost no orientation change for vibrating heavy atoms. This result can be compared to the change in the ellipsoid orientation due to the temperature increase, which reached about 8° in the case of the previously studied Zn complex (Schmøkel et al., 2010). Thus, it may be concluded that the ellipsoid orientation change is rather a subtle effect and should be studied with particular care.

The difference between the Xray and neutron data also seems important here. Clearly, the corresponding ellipsoid sizes are different (Table 5). This observation justifies the scaling of hydrogenatom ADPs, mentioned while describing the multipole refinements. Comparison of these results with the ones previously published is hampered due to the unavailability of a full set of thermal parameters in the literature.
In a standard procedure of the data treatment unitcell parameters are kept fixed while various other parameters are refined. These are usually a scale factor (k_{S}), atomic coordinates (x,y,z), thermal parameters (U_{ij}'s), and, finally, multipolar terms describing the charge redistribution (P_{v} and ) and atomic deformations (P_{lm}'s). Obviously, the accuracy of these parameters is mostly affected by errors coming from two main sources: (a) cell parameters and (b) structurefactor values. It is, however, convenient to treat the refined parameters as independent and random variables in the statistical analyses. In this section, particular attention is paid to the comparative analysis of multipolar terms. The corresponding numerical values are summarized in Table 6.

Standard deviation magnitudes of the P_{v} values seem to be quite large. However, their relative values are not that significant (circa 2% on average) when compared to the P_{v} population magnitudes. On the other hand, the same standard deviations are much more substantial when related to the net charges (i.e., Z  P_{v}) as they amount to about 42% of their values. As a consequence, it is questionable to draw any further conclusions regarding, for example, chargetransfer concepts on the basis of net charges only. On the other hand, from our experience, net charge values may constitute valuable parameters when selecting a proper form factor so as to describe a heavy atom (Domagala et al., 2009; Kaminski et al., 2011).
The hexadecapole populations proved to be an interesting subject of study. Almost all of them are statistically insignificant and equal to zero, except for the P_{40} of the C1 atom (the parameter value is greater than 4). This might suggest that in most of the refinements these parameters could be left unrefined. However, having the proper datatoparameter ratio extending the number of refined parameters with the nonstatistically supported ones (or simply taking the whole expansion up to hexadecapoles) should not be a severe problem. The R_{free} approach, recently introduced into chargedensity analysis by Zarychta et al. (2011) and Paul et al. (2011), would be a very useful tool to be utilized here. However, the idea of rejecting some part of the reflections during the refinement of charge density is questionable. This is due to the significant influence of various reflections on the fine details of the electrondensity distribution. A more appropriate approach would be to determine a `rigid' set of reflections which cannot be taken as a test set in the R_{free} refinement. Nevertheless, the original R_{free} procedure seems to be working quite well in the determination of how loose the chargedensity model is in respect to the data set given.
It is, though, worth noting that the above result is not fully representative and, thus, cannot serve as a basis for some general refinement strategy for other structures. Oxalic acid and water molecules present in the analysed crystal structure are both quite `flat' entities. Therefore, the hexadecapole deformation functions should not be significantly populated. However, in other structures these fine details can be obviously pronounced, especially in the case of heavier atoms located in a more complex environment (e.g. transition metals).
Analysis of localsymmetryforbidden multipole populations constitutes another important aspect. Possible local symmetries, assigned initially by the LSBD program, are given in the second column of Table 6. Consequently, the following constraints would be imposed on the multipolar parameters (when considered up to octupoles): (i) P_{11 + } = P_{11  } = P_{21 + } = P_{21  } = P_{22  } = P_{31 + } = P_{31  } = P_{32  } = P_{33 + } = P_{33  } = 0 (mm2 symmetry), (ii) P_{10} = P_{21 + } = P_{31  } = P_{32  } = 0 (m symmetry) and (iii) P_{11 + } = P_{11  } = P_{21 + } = P_{22 + } = P_{22  } = 0 ( symmetry for hydrogen atoms). As indicated by Table 6, most of these values are statistically insignificant and can be considered as equal to zero (most values do not even exceed the 1 cutoff).
Once again, comparison of our results with those presented by Coppens et al. (1984) and Martin & Pinkerton (1998) sheds some light on electrondensity distribution. In the latter paper, Martin & Pinkerton (1998) found a striking difference between the P_{v} parameters for the C1 atom (P_{v} = 4.31 for the pointdetector study and P_{v} = 3.87 for the CCD). Our results indicate that within the range of 3 (i.e., 4.054 ± 0.279) both values are acceptable. Indeed, as stated by Martin & Pinkerton (1998), such a discrepancy does not affect the total electrondensity distribution significantly.
Analysis of topological properties of the electrondensity distribution is one of the most common undertaken in chargedensity studies. In most cases, various descriptors are computed at bond critical points and provide a firm basis for analysis of chemical bonding properties. It is then extremely important to analyse both the variability of such properties and the stability of a molecular graph. Here, the term `molecular graph' is understood as a set of all nonredundant bond paths present in the crystalline framework (this includes both strong bonds and weak intermolecular contacts).
During the analysis of bond paths for all 13 models we found that the covalent electrondensity topology is highly reproducible and all the bond paths assigned to strong covalent interactions are always present. However, this is not the case for some of the weak intermolecular contacts. Fig. 6 shows the molecular graph of the oxalic acid and water molecules in the crystal. Golden lines indicate the bond paths which are present in all cases and thus constituted the subject for further statistical analysis.
 Figure 6 Bond paths (thick golden lines) and bond critical points (red spheres) present in all models of oxalic acid dihydrate: oxalic acid (a) and water (b) molecule environments. For oxalic acid moieties symmetry transformations are provided only for half of the molecule. 
All bond critical point properties are summarized in Table 7. The most important ones are the electron density and its Laplacian. In addition, the kinetic and potential energy densities, and ellipticity values are also listed. The electrondensity standard deviations constitute an important quantity in the analysis of interactions. Here, for strong covalent bonds they are in the order of 10^{2}, whereas for weaker contacts they are even an order of magnitude smaller. Judging from Table 7 one can see that the relative error for weak contacts is also smaller. Importantly, all electrondensity values show a small spread and are well determined. A similar observation is made for the Laplacian, although some more exceptions are present in this case. For example, the very strong O3H1 hydrogen bond exhibits the Laplacian value of 1 (1) e Å^{5}. This indicates the statistical insignificance of this quantity (i.e., it is not straightforward if it is positive or negative). However, a value close to zero might indicate that the electrostatic character of this contact is counterbalanced by the covalent contribution (the hydrogen bond is very short, Table 4). Similar behaviour is observed for the two remaining hydrogen bonds as well.

The sample mean and sample standard deviation of electron density and its Laplacian can also be computed along the bond path. Fig. 7 shows an exemplary case  the C1O1 bond. As expected, the smallest error is observed at the BCP position and increases as we move away from the BCP. The same behaviour is observed for the relative error values, although it is much less pronounced (see supplementary material).
 Figure 7 Sample means and sample standard deviations evaluated at each point of the C1O1 bond path for (a) electron density (green solid line) and (b) its Laplacian (blue solid line). Error bars (red vertical lines) are shown at the 3 and 1 levels for (a) and (b) panels, respectively. 
Other properties derived from the electron density and its Laplacian values are the kinetic and potential energy densities, which can be computed using the Abramov equation (Abramov, 1997). The Abramov assumption is relevant for the internuclear region analysis (closedshell and sharedshell interatomic interactions) and is particularly important for the chemical bond studies. The potential energy density value at the hydrogenbond BCP obtained via this approach may be utilized further for the hydrogenbond energy estimation by means of the EspinosaLecomte approach (Espinosa et al., 1998, 1999). Therefore, we checked the variability of the kinetic and potential energy density values for the covalent bonds and selected interatomic contacts across the 13 chargedensity models. It appeared that all the energy density values are statistically significant with the average relative error reaching about 10% (with the highest value for the kinetic energy density characterizing the O1H1 bond BCP equal to 41%). The relative sample error is generally noticeably lower for the potential energy density than that for the kinetic energy density for both covalent bonds and intermolecular BCPs. Keeping in mind the above considerations, it can be concluded that the energy density results are consistent, among which especially the potential energy density values are essential for hydrogenbond energy estimation. Furthermore, based on the EspinosaLecomte formula and our results (Table 7), the estimated hydrogenbond energy uncertainty (denoted by sample standard deviation) does not exceed 9 kJ mol^{1} in the case of the selected intermolecular contacts, whereas the relative error does not constitute more than 11% of the mean energy value, which reflects the precision of the derived energies. This makes the analysis of interaction energies, by the use of the EspinosaLecomte approach, feasible.
The last parameter, the ellipticity, seems to be quite sensible. For stronger interactions its values are statistically significant at the 1 level, though at 3 it might be questionable in some cases (e.g., the C1O1 bond). For weak intermolecular contacts, ellipticity is, as might be expected, rather poorly defined.
Integrated atomic charges and atomic volumes are all very important quantities which can be derived directly from the electrondensity distribution. Integrated atomic charges are very frequently used in the analysis of various effects associated with charge redistribution under weak interaction formation, molecule reactivity etc. The numerical values are shown in Table 8. The estimated standard deviations reach values of 0.1 e. The average ratio of standard deviation over the corresponding mean charge amounts to about 5%. In our opinion, this is in the range of possibly observed charge differences between different atoms in the same molecule. Thus, one needs to be very careful to draw meaningful conclusions on the basis of atomic charges obtained experimentally. For instance, in the case of oxalic acid data, the O1 and O2 atoms have practically the same charge (the difference is statistically insignificant), yet they are obviously chemically different. The same thing is true for atomic volumes. In this case the variations are in the range from 1 to 10% even.

At this point it is worth having a more careful look at the total molecular charges. To simplify the discussion the attention is focused solely on the water molecule. The results are summarized in Table 9. The total charge of the water molecule, for the ith data set, is obviously computed as
and the formula stands also for the sample means. The estimator of the standard deviation is derived from the whole population of the computed water charge values (using the basic equation for the standard deviation):
The total charge of the water molecule is statistically equal to zero, even at the 1 level. What is important is the standard deviation regarding the whole water moiety, which is slightly larger than the corresponding standard deviations for each individual atom, i.e., O3, H2 and H3. This suggests some error accumulation while summing up the atomic charges. On the other hand, the atomic charges are not independent variables (their sum must be obviously equal to zero over all atoms), and so the correlation factors may play a significant role (though here not really well observed). This result implies the importance of careful analysis of such charges when it is applied to investigate problems such as `charge transfer' or crystallographic equivalence of molecules in the asymmetric unit. It seems that an approximation of two molecules of the same compound, despite being symmetry nonequivalent, having the same charge (e.g. neutral) should hold rather well. Further modelling of such fine details is within the experimental errors and should not be carried out in most cases. This is in contrast with the recent literature contributions (Nelyubina et al., 2010, 2011), where such an approach is arbitrarily proposed without any deeper insight into the method limitations.

The average value of the dipole moment for the water molecule in oxalic acid dihydrate structure equals 2.1 (3) D, being close to the HartreeFock theoretical calculations (2.2 D) reported by Spackman (1992). The average dipole moment magnitude is also very close to the value reported by Stevens & Coppens (1980) [2.1 (2) D] on the basis of a highresolution chargedensity study. Importantly, the ShapiroWilk test is fulfilled showing the normal distribution of the observed dipole moment magnitudes across different measurements. The computed dipole moments show some spatial spread around their mean position [ = 12 (6)°]. This clearly indicates that the statistical variation is of the same range to the one observed from different multipolar refinements with respect to the same data set (PoulainPaul et al., 2012; Bak et al., 2011).
Second moments of charge distribution were also calculated. The average eigenvalues of the quadrupole moment in the traceless convention are as follows: Q_{xx} = 0.7 (1) e Å^{2}, Q_{yy} = 4.6 (4) e Å^{2} and Q_{zz} = +5.2 (5) e Å^{2}. All mean values are well above the estimated standard deviations. In addition, as might be expected, the first component is very small. It is due to the molecule planarity. Similarly to dipole moments, the ShapiroWilk test always indicates the normal distribution of quadrupole moment values.
In general, the multipolar model provides some information concerning the strength of the electrostatic interactions in the crystal lattice. The electrostatic interaction energies can be simply derived from the electrondensity distribution. However, to draw any binding conclusions on the basis of such results, it is important to check their reliability. Table 10 contains the electrostatic energy values of selected interactions, averaged over the 13 analysed multipolar refinements. The estimated standard deviations are given in parentheses. Additionally, the theoretical results from the PIXEL calculations are listed as a reference point. The calculations were performed for the model oxa8, considered to be the one least deviating in geometry from the mean value. In the case of the theoretical results, the electrostatic energy was summed up together with the polarization component, as polarization should be present in the experimental data, and, thus, to some extent accounted for by the multipolar model (Bak et al., 2012).

The estimated standard deviations are rather high, varying in the range of 517 kJ mol^{1}. Such a result is comparable to the rootmeansquare deviation values obtained for different databases of the transferable aspherical atom models with respect to the corresponding theoretical calculations (Jarzembska & Dominiak, 2012; Bak et al., 2011; Volkov et al., 2006). In the case of oxalic acid, only two electrostatic interactions are significant, i.e., with the electrostatic interaction energy values above 3. This suggests that the electrostatic interactions derived from the multipolar model should be treated more qualitatively than quantitatively, maybe except for the considerably stronger ones.
However, the general trend of the strength of interactions is preserved when compared to the sum of the electrostatic and polarization components from PIXEL calculations. The theoretical values are systematically lower in magnitude, yet within the 2 margin from the multipolar model results. The largest difference is observed for the F interaction being of the type. The correlation between the theoretical model and the experimental one is high (with the linear determination coefficient, R^{2}, being equal to 0.978). Nevertheless, the semiempirical results of PIXEL are not fully consistent with the CRYSTAL results when the total energy is considered. The balance sheet of all the methods perfectly illustrates the accuracy and credibility of the obtained results.
Electrostatic potential (ESP, ) is very frequently used in the analysis of crystal packing and interactions in molecular crystals (Spackman et al., 2008). For this purpose it is usually mapped onto molecular surfaces such as van der Waals or Hirshfeld surfaces. The latter ones are especially convenient as they do not interpenetrate each other, but rather adjoin one another (Spackman & Byrom, 1997b; Spackman & Jayatilaka, 2009). Fig. 8(a) shows the mean electrostatic potential mapped on the Hirshfeld surface of the oxalic acid moiety. The electrostatic potential distribution follows the expected behaviour. The hydrogenbonded regions (i.e., areas closest to the H1 and O2 atoms) show the extreme electrostatic potential values. The side of the molecule (i.e., the area responsible for interaction) seems to be slightly negatively charged. Additional information about the electrostatic potential distribution is gained by the analysis of its standard deviation. It appears that the best way to visualize the standard deviation distribution is to map the quantity onto the same Hirshfeld surface (Fig. 8b). It provides information about the statistical significance of the electrostatic potential values at given points on the surface. It is obvious and confirmed by the data that the hydrogenbonded regions are most significant. The vicinity of the H1 atom shows very large significance (> 10) of the ESP values. Similarly, the negative values are quite significant in the region associated with the O2 atom (> 4.5). It is worth noting that the analysis of the ESP agrees well with the BCP standard deviation summary. For hydrogenbonded regions the standard deviations are very small and, thus, the evaluated properties are well defined. Interestingly, for the region the electrostatic potential is statistically not significant. In this case, the standard deviations and the mean ESP values are of the same range. Thus, the analysis of weaker dispersive interactions should be undertaken with extra care. This statement agrees well with the experimental electrostatic energy values computed for the F dimer and with the results discussed above.
 Figure 8 (a) Mean electrostatic potential () mapped onto the Hirshfeld surface of the oxalic acid molecule. (b) quantity mapped onto the same Hirshfeld surface. 
In order to characterize the ESP distribution, the surface quantities, first introduced by Murray & Politzer (Murray et al., 2000; Murray & Politzer, 1998), were computed. Statistical analysis of the ESP distribution was performed on threedimensional grids with a 0.1 Å sampling step around the whole molecule of oxalic acid with a 3 Å margin. The surface quantities were evaluated on the Hirshfeld surface and around the van der Waals (vdW) surface within a shell of ±0.15 Å thickness. The average values together with sample standard deviations are summarized in Table 11 (for a detailed definition of all parameters see the supplementary material). A comprehensive interpretation of all the ESP parameters is given in the original papers by Murray & Politzer, whereas for the purpose of our study we shall focus on the parameter reliability and overall variability.

Most of the ESP surface quantities are generally well defined at the 3 threshold level. The mean ESP value () is equal to 0.03 (4) e Å^{1}, which correlates well with the slightly negative charge of the oxalic acid molecule, as calculated using the QTAIM approach (Table 9). However, one needs to keep in mind that the value is close to zero and its standard deviation suggests statistical insignificance (this also agrees with the QTAIM result). The average positive (), negative () and total () electrostatic potential is comparable, no matter whether it is calculated on the Hirshfeld or vdW surfaces. Generally speaking, the positive potential region is naturally well defined in the vicinity of the acidic protons and described with the relatively high amplitude of the positive potential values [V_{S,max} = +0.88 (5) e Å^{1}; value for the Hirshfeld surface]. In turn, the negative potential is much more uniformly spread over the remaining parts of the molecule with the pronounced maxima at the oxygen atoms [V_{S,min} = 0.19 (5) e Å^{1}; value for the Hirshfeld surface], but with a larger relative error than that of the positive potential. One can observe that the maximum value of the positive potential and the variance of the positive potential () is higher for the Hirshfeld surface. A larger spread of the positive and negative variances for the Hirshfeld surface leads to the significant difference in the degree of balance parameter (). In the case of the Hirshfeld surface amounts to 0.11 (2), whereas for the vdW it is equal to 0.22 (3), which approaches the theoretical limit of 0.25. The difference between the results obtained for the two surfaces might be associated with the fact that vdW surfaces interpenetrate one another in the crystal frame, and thus may not be an appropriate choice for the purpose of molecular ESP distribution investigations in the solid state.
Finally, we decided to investigate the energetic consequences of slightly different molecular geometries resulting from the obtained 13 chargedensity models. Additionally, two geometry optimizations were carried out in the CRYSTAL program using the periodic DFT(B3LYP) method with two allelectron basis set types, i.e., pVTZ (th1) and 631G** (th2). In the case of all 13 geometries, the total unitcell, molecule and cohesive energies were computed with the pVTZ basis set and analysed. However, it should be stressed that the socalled `molecule energy' stands for the energy of a molecular unit built of one acid moiety and two water molecules (it consists of two asymmetric unit fragments related through the centre of symmetry). As a result, the socalled `cohesive energy' was calculated with respect to the mentioned threemolecule fragment. Such a choice was dictated by technical facility and is sufficient for the purpose of this analysis.
All the results are grouped together in Table 12. The absolute energy values are given in the first three columns. As can be seen, the energy results are very much consistent. Such observation is supported by the corresponding calculated mean energies and standard deviations describing the sample. Standard deviation magnitudes are quite low, as they do not exceed 1 kJ mol^{1} for each of the energy types (unitcell, molecule and cohesive energy). What is more, even a 3 margin is certainly comparable to the error bars of the computation method itself.

The other columns of Table 12 contain the relative energy values with respect to those obtained for the two optimized geometries. In all cases the unitcell energy (E_{bulk}) is lowest in the case of the optimized models, especially for the one with the more advanced basis set used during the optimization procedure (th1). The differences in the unitcell energies range from 5 to almost 26 kJ mol^{1}. On the other hand, both cohesive and molecule energies contribute to the overall stability of the crystal, and so the lowest cohesive energy may not mean the optimal molecular geometry at the same time. Most balanced and energetically closest to the results obtained for th1 and th2 is the oxa10 chargedensity model, while oxa4 and oxa9 seem to be furthest away from the optimal results. However, the energetic differences constitute just a small percentage of the parent energy values.
All in all, on the basis of chargedensity features it was not possible to indicate the optimal model, which is not easily justifiable on the basis of the energetic results either. Nevertheless, it should be concluded that all the energy results are very much consistent and within the error of the methods, which makes them all equally good.
The results presented in this contribution can be summarized as follows:
(i) Raw intensity analysis indicates that the normality distribution is well preserved for about 6070% of measured reflections (up to the resolution of 1.0 Å^{1}). In the case of weak reflections the corresponding percentage is lower (around 50%).
(ii) Unitcell parameters are determined to the order of 10^{3} (for cell edges) and 10^{2} (for angles). The distribution of unitcell parameters with respect to different measurements fulfils the normality test. Values compare well with the older studies of the same compound and with the neutron diffraction data set.
(iii) Fit discrepancy factors are determined within a 0.5% range. Similarly, the standard deviation of the GoF parameter equals about 0.2. The residual density extrema are ±0.16 (3) e Å^{3}, which indicates a rather flat residual chargedensity distribution across the whole unit cell.
(iv) Twodimensional plots show that the mean residual density is very flat and featureless; however, its largest standard deviation values are located near molecular fragments and so it is not completely uniform. The ShapiroWilk test confirms the normal distribution of the experimental errors on Fourier maps.
(v) Deformation density is well reproducible between different measurements. However, the largest errors are located at and near nuclei positions. Moreover, some systematic errors (induced presumably by the multipolar model) are present, as indicated by the normality test.
(vi) Geometry is extremely well reproducible between different data sets. This is confirmed by both direct comparison and periodic computations.
(vii) The largest errors are on the valenceshell chargetransfer parameters (i.e., P_{v}'s), which are further propagated onto the net charges. The symmetry of multipolar expansions, with respect to local coordinate systems, is well preserved.
(viii) The `alwayspresent' bond paths are the ones related to the strong covalent bonds and well defined intermolecular contacts (such as relatively strong hydrogen bonds).
(ix) Standard deviations for electron density at bond critical points are well defined. They are especially small for well defined hydrogenbonded contacts. The same is true for kinetic and potential energy densities, directly related to the interaction strength.
(x) Standard deviations for the integrated atomic charges are equal to about 0.1 e or similar. The charge transfer seems to be a fine detail, which needs to be taken carefully into chemical analyses.
(xi) Dipole moments are comparable with the ones presented in different, older studies. The standard deviation is small (0.2 D); however, the dipole moment direction is less well defined [12 (6)° deviation from the mean dipole moment].
(xii) The electrostatic energies should be treated qualitatively (as suggested in the literature) as large standard deviations occur. On the other hand, they are well correlated with the PIXEL results.
(xiii) Similarly as for the point properties, the ESP distribution is best defined in the hydrogenbonded regions.
(xiv) The cohesive energy values obtained for crystal geometries resulting from the studied data sets show no significant differences. Therefore, the molecular geometry is not much affected by the final model.
On the whole, experimental chargedensity studies are useful in several respects. Firstly, the geometry is well determined after multipolar refinement, so one could use it as a subject for further theoretical computations. The derived properties, such as bond critical point electron densities, are extremely useful in the analysis of moderatestrength intermolecular contacts, occurring in the solid state. This conclusion may be extended to a weak nonclassical bonding occurring in coordination complexes, assuring the proper description of the metal centre. Particular attention should be paid to the analysis of integrated charges and, in the case of dynamical systems (phase transitions etc.), unitcell constants. Hopefully, these quantities can be much improved by the recent advances in pixel detectors. Integrated charges should be analysed with great care, especially while considering chargetransfer issues. As predicted, the very useful properties such as dipole moments, ESP distribution and electrostatic energies should be treated rather qualitatively. Of course, some quantities are affected by the multipolar model itself. Thus, it would be recommended to conduct a similar analysis using different electrondensity reconstruction methods in the future, as the awareness of the model and method limitations is of prime importance. Nevertheless, we clearly show here that the source of errors is Gaussiannoiselike and it is particularly related to the experimental setup and methodology inaccuracies. It should also be noted that the variations of the derived properties will likely be larger when measured on a greater variety of instruments. We do believe that the extension of our study to different types of detectors [e.g. complementary metaloxide superconductor (CMOS) or other pixelarray devices] and Xray sources (e.g. synchrotron radiation) will be available in the future and the presented results can be readily updated.
The authors gratefully acknowledge the Wroclaw Centre for Networking and Supercomputing for providing computer facilities. RK would like to thank the Foundation for Polish Science for financial support within the `International PhD Projects' programme. KNJ would like to acknowledge the Foundation for Polish Science for financial support within the `START' programme. KW acknowledges the Ministry of Science and Higher Education/National Science Centre (grant No. N204 135138) for financial support. The authors thank Krzysztof Durka (Warszawa, Poland) for growing excellent quality oxalic acid crystals for the neutron diffraction experiment.
Abrahams, S. C., Hamilton, W. C. & Mathieson, A. McL. (1970). Acta Cryst. A26, 118.
Abramov, Yu. A. (1997). Acta Cryst. A53, 264272.
Abramov, Y. A., Volkov, A., Wu, G. & Coppens, P. (2000). J. Phys. Chem. B, 104, 21832188.
Agilent Technologies (2010). CRYSALIS PRO. Agilent Technologies, Yarnton, Oxfordshire, UK.
Allen, F. H. (2002). Acta Cryst. B58, 380388.
Bader, R. F. W. (1994). Atoms in Molecules: A Quantum Theory. New York: Oxford University Press.
Bak, J. M., Czyznikowska, Z. & Dominiak, P. M. (2012). Acta Cryst. A68, 705714.
Bak, J. M., Domagala, S., Hübschle, C., Jelsch, C., Dittrich, B. & Dominiak, P. M. (2011). Acta Cryst. A67, 141153.
Bak, J. M., Dominiak, P. M., Wilson, C. C. & Wozniak, K. (2009). Acta Cryst. A65, 490500.
Becke, A. D. (1988a). J. Chem. Phys. 88, 25472553.
Becke, A. D. (1988b). Phys. Rev. A, 38, 30983100.
Becker, P. J. & Coppens, P. (1974a). Acta Cryst. A30, 148153.
Becker, P. J. & Coppens, P. (1974b). Acta Cryst. A30, 129147.
Becker, P. J. & Coppens, P. (1975). Acta Cryst. A31, 417425.
Benedict, J. B., Makal, A., Sokolow, J. D., Trzop, E., Scheins, S., Henning, R., Graber, T. & Coppens, P. (2011). Chem. Commun. 47, 17041706.
Blessing, R. H. (1987). Crystallogr. Rev. 1, 358.
Blessing, R. H. (1989). J. Appl. Cryst. 22, 396397.
Blessing, R. H. (1995a). Acta Cryst. A51, 3338.
Blessing, R. H. (1995b). Acta Cryst. B51, 816823.
Brandt, S. (1998). Data Analysis: Statistical and Computational Methods for Scientists and Engineers. London: Springer.
Bruker AXS Inc. (2013). APEX2. Madison, Wisconsin, USA.
Busing, W. R. & Levy, H. A. (1958). Acta Cryst. 11, 450451.
Civalleri, B., ZicovichWilson, C. M., Valenzano, L. & Ugliengo, P. (2008). CrystEngComm, 10, 405410.
Cochran, W. G. (1934). Math. Proc. Camb. Philos. Soc. 30, 178191.
Coppens, P. (1997). Xray Charge Densities and Chemical Bonding. New York: Oxford University Press, Inc.
Coppens, P. (2005). Angew. Chem. Int. Ed. 44, 68106811.
Coppens, P., Dam, J., Harkema, S., Feil, D., Feld, R., Lehmann, M. S., Goddard, R., Krüger, C., Hellner, E., Johansen, H., Larsen, F. K., Koetzle, T. F., McMullan, R. K., Maslen, E. N. & Stevens, E. D. (1984). Acta Cryst. A40, 184195.
Cruickshank, D. W. J. (1956). Acta Cryst. 9, 754756.
Dittrich, B., Sze, E., Holstein, J. J., Hübschle, C. B. & Jayatilaka, D. (2012). Acta Cryst. A68, 435442.
Domagala, S., KorybutDaszkiewicz, B., Straver, L. & Wozniak, K. (2009). Inorg. Chem. 48, 40104020.
Dovesi, R., Orlando, R., Civalleri, B., Roetti, C., Saunders, V. R. & ZicovichWilson, C. M. (2005). Z. Kristallogr. 220, 571573.
Dovesi, R., Saunders, V. R., Roetti, C., Orlando, R., ZicovichWilson, C. M., Pascale, F., Civalleri, B., Doll, K., Harrison, N. M., Bush, I. J., D'Arco, P. & Llunell, M. (2009). CRYSTAL09 User's Manual. University of Turin, Turin, Italy.
Dunitz, J. D. (1996). Xray Analysis and the Structure of Organic Molecules. Basel, Weinheim, New York: WileyVCH.
Dunning, T. H. (1989). J. Chem. Phys. 90, 10071023.
El Haouzi, A., Hansen, N. K., Le Hénaff, C. & Protas, J. (1996). Acta Cryst. A52, 291301.
Espinosa, E., Lecomte, C. & Molins, E. (1999). Chem. Phys. Lett. 300, 745748.
Espinosa, E., Molins, E. & Lecomte, C. (1998). Chem. Phys. Lett. 285, 170173.
Farrugia, L. J. & Evans, C. (2005). J. Phys. Chem. A, 109, 88348848.
Feast, G. C., Haestier, J., Page, L. W., Robertson, J., Thompson, A. L. & Watkin, D. J. (2009). Acta Cryst. C65, o635o638.
Fischer, A., Tiana, D., Scherer, W., Batke, K., Eickerling, G., Svendsen, H., Bindzus, N. & Iversen, B. B. (2011). J. Phys. Chem. A, 115, 1306113071.
Frisch, M. J. et al. (2009). GAUSSIAN09. Gaussian Inc., Pittsburgh, Pennsylvania, USA.
Gavezzotti, A. (2002). J. Phys. Chem. B, 106, 41454154.
Gavezzotti, A. (2003). J. Phys. Chem. B, 107, 23442353.
Grabowsky, S., Kalinowski, R., Weber, M., Förster, D., Paulmann, C. & Luger, P. (2009). Acta Cryst. B65, 488501.
Grimme, S. (2004). J. Comput. Chem. 25, 14631473.
Grimme, S. (2006). J. Comput. Chem. 27, 17871799.
Guillot, B. (2012). Acta Cryst. A68, s204.
Guzei, I. A., Bikzhanova, G. A., Spencer, L. C., Timofeeva, T. V., Kinnibrugh, T. L. & Campana, C. F. (2008). Cryst. Growth Des. 8, 24112418.
Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909921.
Herbstein, F. H. (2000). Acta Cryst. B56, 547557.
Hirshfeld, F. L. (1976). Acta Cryst. A32, 239244.
Hirshfeld, F. L. & Rabinovich, D. (1973). Acta Cryst. A29, 510513.
Hoser, A. A., Dominiak, P. M. & Wozniak, K. (2009). Acta Cryst. A65, 300311.
Howard, S. T., Hursthouse, M. B., Lehmann, C. W. & Poyner, E. A. (1995). Acta Cryst. B51, 328337.
Jarzembska, K. N. & Dominiak, P. M. (2012). Acta Cryst. A68, 139147.
Jarzembska, K. N., Kubsik, M., Kaminski, R., Wozniak, K. & Dominiak, P. M. (2012). Cryst. Growth Des. 12, 25082524.
Jayatilaka, D. & Grimwood, D. J. (2001). Acta Cryst. A57, 7686.
Jelsch, C., Guillot, B., Lagoutte, A. & Lecomte, C. (2005). J. Appl. Cryst. 38, 3854.
Kaminski, R., Herbaczynska, B., Srebro, M., Pietrzykowski, A., Michalak, A., Jerzykiewicz, L. B. & Wozniak, K. (2011). Phys. Chem. Chem. Phys. 13, 1028010284.
Kaminski, R., Jarzembska, K. N. & Domagala, S. (2013). J. Appl. Cryst. 46, 540543.
Keen, D. A., Gutmann, M. J. & Wilson, C. C. (2006). J. Appl. Cryst. 39, 714722.
Krishnan, R., Binkley, J. S., Seeger, R. & Pople, J. A. (1980). J. Chem. Phys. 72, 650654.
Kuntzinger, S., Dahaoui, S., Ghermani, N. E., Lecomte, C. & Howard, J. A. K. (1999). Acta Cryst. B55, 867881.
Lebedev, V. I. & Laikov, D. N. (1999). Dokl. Math. 59, 477481.
Lee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785789.
Madsen, A. Ø. (2006). J. Appl. Cryst. 39, 757758.
Madsen, A. & Larsen, S. (2007). Angew. Chem. Int. Ed. 46, 86098613.
Madsen, A. Ø, Sørensen, H. O., Flensburg, C., Stewart, R. F. & Larsen, S. (2004). Acta Cryst. A60, 550561.
Makal, A., Benedict, J., Trzop, E., Sokolow, J., Fournier, B., Chen, Y., Kalinowski, J. A., Graber, T., Henning, R. & Coppens, P. (2012). J. Phys. Chem. A, 116, 33593365.
Makal, A., Trzop, E., Sokolow, J., Kalinowski, J., Benedict, J. & Coppens, P. (2011). Acta Cryst. A67, 319326.
Martin, A. & Pinkerton, A. A. (1998). Acta Cryst. B54, 471477.
Meindl, K. & Henn, J. (2008). Acta Cryst. A64, 404418.
Meindl, K., Henn, J., Kocher, N., Leusser, D., Zachariasse, K. A., Sheldrick, G. M., Koritsanszky, T. & Stalke, D. (2009). J. Phys. Chem. A, 113, 96849691.
Møller, C. & Plesset, M. S. (1934). Phys. Rev. 46, 618622.
Munshi, P., Madsen, A. Ø., Spackman, M. A., Larsen, S. & Destro, R. (2008). Acta Cryst. A64, 465475.
Murray, J. S., PeraltaInga, Z. & Politzer, P. (2000). Int. J. Quantum Chem. 80, 12161223.
Murray, J. S. & Politzer, P. (1998). J. Mol. Struct. (Theochem.), 425, 107114.
Nelyubina, Y. V., Antipin, M. Y., Cherepanov, I. A. & Lyssenko, K. A. (2010). CrystEngComm, 12, 7781.
Nelyubina, Y. V., Dalinger, I. L. & Lyssenko, K. A. (2011). Angew. Chem. Int. Ed. 50, 28922894.
Paciorek, W. A., Meyer, M. & Chapuis, G. (1999). Acta Cryst. A55, 543557.
Paul, A., Kubicki, M., Jelsch, C., Durand, P. & Lecomte, C. (2011). Acta Cryst. B67, 365378.
Petrícek, V., Dusek, M. & Palatinus, L. (2006). JANA2006. Institute of Physics, Czech Academy of Sciences, Prague, Czech Republic.
PoulainPaul, A., Nassour, A., Jelsch, C., Guillot, B., Kubicki, M. & Lecomte, C. (2012). Acta Cryst. A68, 715728.
Sands, D. E. (1982). Vectors and Tensors in Crystallography. London, Amsterdam, Toronto, Sydney, Tokyo: AddisonWesley Publishing Company.
Schmøkel, M. S., Bjerg, L., Overgaard, J., Larsen, F. K., Madsen, G. K. H., Sugimoto, K., Takata, M. & Iversen, B. B. (2013). Angew. Chem. Int. Ed. 52, 15031506.
Schmøkel, M. S., Cenedese, S., Overgaard, J., Jørgensen, M. R., Chen, Y. S., Gatti, C., Stalke, D. & Iversen, B. B. (2012). Inorg. Chem. 51, 86078616.
Schmøkel, M. S., Kaminski, R., Benedict, J. B. & Coppens, P. (2010). Acta Cryst. A66, 632636.
Schomaker, V. & Trueblood, K. N. (1968). Acta Cryst. B24, 6376.
Schwarzenbach, D., Abrahams, S. C., Flack, H. D., Gonschorek, W., Hahn, Th., Huml, K., Marsh, R. E., Prince, E., Robertson, B. E., Rollett, J. S. & Wilson, A. J. C. (1989). Acta Cryst. A45, 6375.
Seiler, P., Schweizer, W. B. & Dunitz, J. D. (1984). Acta Cryst. B40, 319327.
Shapiro, S. S. & Wilk, M. B. (1965). Biometrika, 52, 591611.
Spackman, M. A. (1992). Chem. Rev. 92, 17691797.
Spackman, M. A. & Byrom, P. G. (1996). Acta Cryst. B52, 10231035.
Spackman, M. A. & Byrom, P. G. (1997a). Acta Cryst. B53, 553564.
Spackman, M. A. & Byrom, P. G. (1997b). Chem. Phys. Lett. 267, 215220.
Spackman, M. A. & Jayatilaka, D. (2009). CrystEngComm, 11, 1932.
Spackman, M. A., McKinnon, J. J. & Jayatilaka, D. (2008). CrystEngComm, 10, 377388.
Stash, A. & Tsirelson, V. (2002). J. Appl. Cryst. 35, 371373.
Stash, A. I. & Tsirelson, V. G. (2005). Crystallogr. Rep. 50, 177184.
Stevens, E. D. & Coppens, P. (1980). Acta Cryst. B36, 18641876.
Stevens, E., Coppens, P., Feld, R. & Lehmann, M. (1979). Chem. Phys. Lett. 67, 541543.
Su, Z. & Coppens, P. (1998). Acta Cryst. A54, 646652.
Treutler, O. & Ahlrichs, R. (1995). J. Chem. Phys. 102, 346.
Volkov, A., Abramov, Y. A. & Coppens, P. (2001). Acta Cryst. A57, 272282.
Volkov, A., King, H. F. & Coppens, P. (2006). J. Chem. Theory Comput. 2, 8189.
Volkov, A., Li, X., Koritsanszky, T. & Coppens, P. (2004). J. Phys. Chem. A, 108, 42834300.
Waser, J. (1955). Acta Cryst. 8, 731.
Watkin, D. (2008). J. Appl. Cryst. 41, 491522.
Zarychta, B., Zaleski, J., Kyziol, J., Daszkiewicz, Z. & Jelsch, C. (2011). Acta Cryst. B67, 250262.
Zhurov, V. V., Zhurova, E. A. & Pinkerton, A. A. (2011). Inorg. Chem. 50, 63306333.
Zhurov, V. V., Zhurova, E. A., Stash, A. I. & Pinkerton, A. A. (2011a). Acta Cryst. A67, 160173.
Zhurov, V. V., Zhurova, E. A., Stash, A. I. & Pinkerton, A. A. (2011b). J. Phys. Chem. A, 115, 1301613023.