Experimental and Theoretical Charge Density Studies at Subatomic

ARTICLE pubs.acs.org/JPCA

Experimental and Theoretical Charge Density Studies at Subatomic Resolution A. Fischer,† D. Tiana,† W. Scherer,*,† K. Batke,† G. Eickerling,*,† H. Svendsen,‡ N. Bindzus,‡ and B. B. Iversen*,‡ † ‡

Institute of Physics, University of Augsburg, Universit€atsstrasse 1, 86135 Augsburg, Germany Center for Materials Crystallography, Department of Chemistry and iNANO, Aarhus University, DK-8000 Aarhus, Denmark

bS Supporting Information ABSTRACT: Analysis of accurate experimental and theoretical structure factors of diamond and silicon reveals that the contraction of the core shell due to covalent bond formation causes significant perturbations of the total charge density that cannot be ignored in precise charge density studies. We outline that the nature and origin of core contraction/expansion and core polarization phenomena can be analyzed by experimental studies employing an extended Hansen-Coppens multipolar model. Omission or insufficient treatment of these subatomic charge density phenomena might yield erroneous thermal displacement parameters and high residual densities in multipolar refinements. Our detailed studies therefore suggest that the refinement of contraction/expansion and population parameters of all atomic shells is essential to the precise reconstruction of electron density distributions by a multipolar model. Furthermore, our results imply that also the polarization of the inner shells needs to be adopted, especially in cases where second row or even heavier elements are involved in covalent bonding. These theoretical studies are supported by direct multipolar refinements of X-ray powder diffraction data of diamond obtained from a third-generation synchrotron-radiation source (SPring-8, BL02B2).

1. INTRODUCTION The formalism suggested by Stewart in 1977,1 providing a decomposition of the total charge density into pseudoatoms, and the development of the popular Hansen-Coppens (HC) multipolar model in 19782 offer highly efficient tools to extract precise charge density distributions F(r) from experimental X-ray and electron diffraction data. However, despite the enormous successes of these pioneering concepts to analyze the complex nature of chemical bonding, the validation of the reliability and reproducibility of experimental charge density studies remains an important task.3 This is especially true in light of the fast development of diffraction techniques and the tremendous efforts in data reduction/correction methods that push the limits of experimental charge density studies beyond routine analyses of chemical bonding effects. Indeed, the study of physical properties4
7 and the understanding of chemical functionality (e.g., the activation of chemical bonds8
10 at variable conditions (temperature, pressure, or external fields)) usually demands the analysis of very tiny charge density variations, not only in the internuclear bonding domain (e.g., at the critical points), but also in the valence shell of atoms. However, recent evaluations of the standard HC multipolar model revealed its limited capability to fit theoretical structure factors. The authors concluded that the failure of the HC model is due to the limited flexibility of the employed single-ζ radial function.11,12 Volkov et al. therefore proposed the usage of a double-ζ multipolar model (denoted DZ-HC model in the following) to “faithfully retrieve a theoretical charge r 2011 American Chemical Society

density”.13 Indeed, the authors observed a “significant improvement of the topological properties such as the interatomic profiles of F(r), r2F(r), and λ3, atomic charges and volumes, and molecular dipole moments” upon usage of the DZ-HC model relative to the standard HC multipolar model. However, the practical applicability of the DZ-HC model remains to be demonstrated and will depend on the accessibility of data at uttermost resolution and precision.14,15 An even more sophisticated concept has been suggested by Koritsansky and Volkov16,17 by combining Stewart’s pseudoatom formalism1 and the stockholder partitioning scheme to derive radial distribution functions (RDF) for chemically bonded atoms in molecules and solids. The authors state that the projection of these stockholder-atom RDFs onto a small auxiliary STO basis allows in principle a highly precise representation of the total molecular density of chemically bonded atoms. Unfortunately, also this concept requires several Slater primitives per pseudoatom even in the case of first row elements. The situation is further complicated by the limited transferability of the stockholder-atom RDFs for pseudo atoms in chemically different environments. Furthermore, the elegance and simplicity of the standard HC Special Issue: Richard F. W. Bader Festschrift Received: May 30, 2011 Revised: August 2, 2011 Published: August 24, 2011 13061

dx.doi.org/10.1021/jp2050405 | J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A

ARTICLE

Table 1. Parameters Employed in the Combined Multipolar/ Rietveld Refinements of Diamond HC modela

kcore modelb

a

3.56932(2)

3.56932(2)

zero shift GU

0.163(8) 33(2)

0.163(8) 32(2)

GV

16(1)

16(1)

GW

1.37(9)

1.37(1)

LX

1.26(3)

1.25(3)

LY

0.1(2)

0.1(2)

R

0.45

0.47

Rw

0.66

0.67

Rp Rwp

2.45 4.82

2.45 4.82

a

VM radial functions of a free, C(s2p2) configured pseudoatom were employed. b SCM radial functions of an C(sp3) hybridized pseudoatom were employed.

multipolar model which hitherto allowed a chemical and physical interpretation of its parameters (and their plausibility) has been replaced by a fitting method which necessitates a priori theoretical calculations for each specific pseudoatom in its individual chemical environment. The latter point is especially critical for charge density studies of materials showing complex physical properties which depend on various control parameters (e.g., temperature, pressure, or external fields). Hence, a concept is needed that allows improvement to the flexibility of the standard HC multipolar model in a highly efficient way by a strictly limited number of additional parameters. We will demonstrate in the following that already the refinement of individual expansion/contraction parameters of the core and outer core
shells leads to a tremendous reduction of the residual densities and prevents refinements of physically unreasonable atomic displacement parameters.

2. METHODS Theoretical structure factors of diamond and R-silicon were obtained from full-potential linearized augmented plane-wave DFT calculations employing the WIEN2K (W2K) program18 and using experimental crystal structure information.19,32 In the case of diamond and R-silicon a mesh of 455 irreducible k-points in the Brillouin-Zone and the PBE GGA-functional have been used for the SCF calculations.20 Static X-ray structure factors (sin θ/λ < 1.8 Å
1) were calculated employing the LAPW3 routine of WIEN2K. Additional electronic structure calculations were performed using the LCAO approach and the PBE/6-311G(3d) method as implemented in Crystal06 (C06).21 The calculations were performed using a tight radial and angular quadrature with 455 points inside the Brillouin zone. In order to assist the SCF convergence, the eigenvalue level shifting technique was used. The structure factors were calculated using minimal (STOnG: n = 2
6) and split-valence Gaussian type orbitals (GTO): 3-21G, 4-31G, 6-21G, 6-31G, 6-311G. GTO structure factors were also calculated adding polarization functions (3-21G(d), 4-31G(d), 6-21G(d), 6-31G(d), 6-311G(d), 6-311G(2d), 6-311G(3d), 6-311G(3df)). The STOnG (n = 2
6) basis sets 3-21G, 3-21G(d), 6-21G, 6-21G(d) were used with standard coefficients as implemented in the Crystal06 code. For the 4-31G, 4-31G(d), 6-31G, 6-31G(d), 6-311G, 6-311G(d), 6-311G(2d), 6-311G(3d),

and 6-311G(3df) bases the coefficients were taken from literature22
25 and subsequently contracted to avoid any basis set linear dependence error and an unphysical conducting state for diamond. In case of the 6-21G(d) bases the structure factors were also calculated using different functionals (B3PW, BPW, BLYP, B3LYP, PBE0) and at the Hartree
Fock level (Supporting Information). For comparison reasons with earlier studies by Koritsansky et al.17 the ADF program26
28 was used to calculate the charge density distribution of R-glycine (BLYP/ QZ4P level of approximation). The DenProp code29 was employed to calculate the static X-ray structure factors using a numerical quadrature-based Fourier-transform of the molecular electron density (ED) of R-glycine in the STO representations and for the analysis of the topology of F(r). Molecular calculations on the model compound Si(SiH3)4 have been performed employing the ADF program and the BLYP/QZ4P/ZORA level of approximation. The geometry was optimized assuming Td symmetry and the wave function was analyzed with the NBO program.30 Combined Rietveld and multipolar refinements on the experimental powder diffraction data of diamond were performed with JANA2006.31 Refinements were based on F values using the damping method and unit weight settings. The background was fitted employing a 12th order Legendre polynomial in the range 8° < 2θ < 71.3°. The peak shape function was of pseudo-Voigt type introducing Gaussian type profile parameters (GU, GV, GW) and Lorenzian type parameters (LX, LY). The peak calculation cutoff was at 25 times the full width at half-maximum (fwhm) and no peak asymmetry parameters were used. The wavelength was kept fixed to the value determined by Nishibori et al. (λ = 0.40112 Å)32 assuming linearly polarized synchrotron radiation. A single linear shift parameter in combination with the unit cell constant was refined in the addition. The multipolar parameters and the isotropic thermal displacement parameter at the carbon atom were relaxed at the last steps of the refinements and optimized simultaneously with all Rietveld parameters employed (see Tables 1 and 2). The reconstruction of the experimental and theoretical charge density by multipolar refinements was performed using JANA200631 (and for comparison reason with the XD program;33 Table S3 in the Supporting Material). The Volkov and Macchi (VM)33 scattering data bank obtained from relativistic PBE/QZ4P/ZORA and the one by Su, Coppens, and Macchi (SCM)34,35 based on relativistic Dirac
Fock calculations, were used throughout. In case of diamond, the performance of different radial functions (RF) of s2p2 and sp3 hybridized carbon atoms (denoted C and Cv, respectively) were tested using the SCM scattering tables. The RF of the latter pseudoatom significantly improved the results based on experimental and theoretical structure factor fittings in agreement with earlier findings by Svendsen et al.19 Note that in case of the diamond or R-silicon we employed averaged single ζ density exponents in accord with the presumed sp3 configuration. All experimental and theoretical multipolar refinements were performed at a resolution of sin θ/λ < 1.45 Å
1 (Diamond) and sin θ/λ < 1.8 Å
1 (R-Silicon, R-glycine)in order to allow detailed studies of core expansion/contraction and polarization effects. Two independent multipolar functions (O2
, H0) were allowed by symmetry up to fourth order on the carbon and silicon atoms in diamond and R-silicium, while the hexadecapole H4+ was constrained to H0 by H4+ = (0.74045)  H0, employing the cubic symmetry. The local coordinate system was oriented such that a positive sign of the multipole term O2
signals a charge accumulation between a 13062

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A

ARTICLE

Table 2. Multipole Parameter and Refinement Results Based on Theoretical Fsta (W2K) Values, Experimental High Resolution Powder Diffraction (SPring8), and Single Crystal Data (Pendell€ osung method)a W2K/C06 wavefunction

W2K (HC)b

W2K (HC+kcore)c

SPring8 (HC) b

SPring8 (HC+ kcore) c

Pendelm. d (HC)

R/RF (%)

0.40

0.05

0.45

0.47

ΔF(r) [eÅ
3]

+0.22
0.12

+0.02
0.01

+0.06
0.14

+0.09
0.08

Uiso

0.0e

0.0e

0.00161(5)

0.00176(5)

scale

1.0024

1.0e

6.17

6.23

kvalence

0.947

0.970

0.950(8)

0.954(7)

k0 deformation

0.81

0.873

0.80(2)

0.86(2)

Pv

4.0e

4.013

4.0e

4.013f

4.0e

1.006

e

f

1.006

1.0e

e

f

kcore

1.0

e e

Pc O2


1.0

0.85

0.999(18)

2.0 0.42

1.987 0.338

2.0 0.44(2)

1.987 0.37(1)

2.0e 0.34(19) 0.0(1)


0.20


0.108


0.27(5)


0.18(3)

F(r)BCP [eÅ
3]

1.61/1.59

1.63

1.60

1.68

1.63

r2F(r)BCP [eÅ
5]


17.5/
11.2


13.6


11.5


14.9


14.4

λ3 [eÅ
5]

2.8/9.3

7.0

9.0

7.0

7.6

r2F(r)BCC [eÅ
5]


21.6/
20.4


19.9


21.6


21.0


23.3

H0

0.00177(13)

a

Additional models are described in the Supporting Information. b Standard multipolar refinements based on the radial functions of a free carbon atom (s2p2 configuration, 3P state; VM database). c Core k refinement, including the adoption of the core population Pc and employing the radial functions of an sp3 hybridized carbon atom (SCM database). d HC multipolar model (denoted B, O, H, k, R model in ref 46) using the Clementi and Roetti database (refs 47 and 48) and structure factors based on the Pendell€osung method (ref 50). e Unrefined value. f Values fixed and adopted from the theoretical (HC+kcore) model.

Figure 1. (a) Residual density, ΔF(r), map in the (110) plane of the diamond lattice based on static structure factors Fsta (W2K) after standard HC multipolar refinements (Uiso = 0.0 Å2; scale =1.0024); (b) ΔF(r) map after convoluting the Fsta (W2K) values with an unphysical Uiso parameter of
0.00012 Å2 (scale = 0.9934); (c) ΔF(r) map after refining a core kappa (kcore = 1.014) and a core monopole parameter Pc = 1.979 (Uiso = 0.0; scale =1.0); (d) resulting ΔF(r) map upon replacement of the VM radial functions of a free (s2p2 configured) carbon atom (a
c) by the SCM radial function of an sp3 hybridized carbon atom (kcore = 1.006, Pc = 1.987, Uiso = 0.0, scale = 1.0). Positive (solid) and negative (dashed) contour lines are drawn employing a step width of 0.01 eÅ
3 at a resolution of sin θ/λ < 1.45 Å
1.

carbon atom at (0, 0, 0) and its nearest bonding neighbors (for further details see Supporting Information).

3. RESULTS AND DISCUSSION 3.1. Theoretical Studies of the Core Shell Contraction in Diamond. In a recent publication, Svendsen et al. outlined that

standard multipolar refinements yield remarkably large residual densities when applied to theoretical structure factors (sin θ/λ < 1.45 Å
1) of diamond.19 The static theoretical structure factors (Fsta) of this earlier study were obtained from an all-electron density functional theory (DFT) calculation employing Bl€ochl’s PAW method36 within the local density approximation (LDA) and the Perdew
Burke
Ernzerhof (PBE) gradient correction.20 The 13063

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A observed residual density maps showed pronounced holes and peaks (ΔF(r) =
0.065
0.121 eÅ
3) in the core and valence shell of the carbon atom. This might suggest a deficiency of the standard HC multipolar model to recover precise theoretical charge density distributions from calculated structure factors.19 However, experimental multipolar refinements of diamond based on high-resolution X-ray powder diffraction data (SPring-8, BL02B2; 100 K data) yielded clearly lower residuals in the same data resolution range, as provided by the theoretical study (
0.051
0.086 eÅ
3).19 In the light of the inherent error sources of experimental data, the latter result is somewhat surprising. Apparently, the additional fitting of an isotropic thermal displacement parameter of the carbon atom in the experimental study adjusts the quality of the data fit such that it surpasses the respective fitting quality of the HC multipolar model based on the theoretical Fsta (PAW) data. Furthermore, the authors could demonstrate that convolution of an isotropic thermal displacement parameter Uiso of 0.00188 Å2 on the static theoretical data (denoted Fdyn (PAW)) and subsequent refinements also yields a significant flattening of the residual density maps. Careful inspection of the results reveals, however, that the impressive fit improvements based on the Fdyn (PAW) data (ΔF(r) =
0.018
0.016 eÅ
3) originate from a physically unreasonable (i) reduction of the Uiso value by approximately 3% and (ii) a scale factor deviating slightly from unity (0.9994(3)). This implies that refinements of a scale factor and thermal displacement parameters are a misleading approach to compensate potential deficiencies of the standard HC multipolar model in fitting theoretical structure factor data. For a systematic improvement of the standard HC model, it is therefore necessary to understand the origin of the unphysical scale and thermal displacement parameters obtained by refinements of theoretical Fdyn data. Especially, because this phenomenon might also affect and falsify experimental thermal displacement parameters. Figure 1a shows the residual density map in the (110) plane of the diamond after fitting the static theoretical structure factors Fsta (W2K) by a HC multipolar refinement employing the VM scattering factor database (for details, see Methods and Table 2). In these refinements, the theoretical data were based on periodic full-potential LAPW DFT calculations employing the WIEN2K program.18 The good agreement of the residual map (
0.12/ 0.22 eÅ
3) with the one obtained by Svendsen et al.19 (vide supra) suggests that the deficiencies of the standard multipolar model are systematic in origin and rather independent from the theoretical method employed (see Table 2 and Supporting Information for more examples). As in the case of the refinements on Fdyn structure factors (vide infra), we could significantly reduce the residual density (
0.03/0.03 eÅ
3) by scaling the Fsta (W2K) data and refining an isotropic Uiso,C-value (Figure 1b). However, the rather flat residuals are again obtained at the virtue of an unphysical (negative) Uiso value of
0.00012 Å2 and a scale factor deviating from unity (0.993). In standard X-ray analyses, unphysical small temperature factors might hint for a mistaken identity (too small Z-value) of the refined atom. In the present more subtle case, however, the refinement results (negative Uiso-value in combination with residual density accumulations in the core and valence shell region) might hint for a deficiency of the standard multipolar model to adjust the shell structure of the carbon atom in the diamond. Indeed, comparison of the theoretical radial charge density distribution of an isolated carbon atom with a

ARTICLE

Figure 2. (a) Radial charge density distribution F(r) (au) in the core region of the free carbon atom (black line) and the tetra-coordinated carbon atom in the diamond (along the CC bonding vector; PBE/6311G(3d) level of approximation, red line). Note that the minute flattening of radial F(r) profiles in the cusp region is an artifact originating from the usage of GTO basis sets of triple-ζ quality. The extension of the core shell has been marked by a dashed line that marks the first spherical node (zero-crossing) in the L(r) =
r2F(r) profile of the free carbon atom in its ground state (r = 0.173
1/Z(C); ref 38); (b) Radial F(r) profiles of the sp3-hybridized carbon atom (red line, SCM database) in comparison with the free carbon atom in its ground state (black line, VM database).

covalently bonded carbon atom in the diamond (PBE/6-311G(3d) level of approximation) identifies, in the case of the diamond, a charge depletion in the core and cusp region of the carbon atoms (F(r)cusp = 121.2 and 119.5 au for Cisolated and Cdiamond , respectively;37 Figure 2a). This core contraction is a direct consequence of the charge density accumulation in the C
C bonding domain of the diamond by involvement of the 2s atomic orbitals (AO). In contrast to the 2p orbital, the 2s AO displays a radial node and thus contributes significantly to the core density. The covalent C
C bond formation process is thus inherently linked with the core contraction phenomenon in diamond. Such a process can therefore not be described by the mere adjustment of the radial screening parameters kvalence and k0 deformation of the valence and deformation density distributions in the standard HC multipolar model. We further note that the rather isotropic core contraction of the carbon atom affects the scattering intensity with increasing diffraction angles in a similar way as a decreasing isotropic Debye
Waller factor. As a consequence, the scale factor smaller than unity and the negative Uiso value (Figure 1b) adjust, in an unphysical way, (i) the core contraction and (ii) the reduced cusp density of the carbon atoms in diamond. Hence, an extension of the standard model is needed by implementation of a core k model (1) that allows the fitting of (i) a radial screening parameter of the core density, kcore, and the simultaneous refinement of (ii) a core population parameter Pc FðHÞ ¼

∑j ½fPj, cfj, c ðH=kcore Þ þ Pj, v fj, v ðH=kvalenceÞ ge2πiH 3 rj 

ð1Þ

where Pj,c and Pj,v are the independent population parameters of the core and valence electrons of atom j. Hence, kcore and kvalence denote expansion/contraction parameters of the respective shells. In this respect, our model goes beyond earlier attempts to merely adjust the core contraction/expansion of light elements (e.g., Be) by high-order refinements of a single k-parameter for both core and valence shells.39 Indeed, the necessity of a variable core population has been addressed by a number of pioneering 13064

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A studies,40
42 and we will show later that for second row and heavier elements even individual Pj,c and kcore parameters need to be adjusted to fit the expansion/contraction and deformation of each shell of bonded atoms in molecules and solids. In case of diamond, core k refinements of Fsta(W2K) yield a kcore value of 1.014 in combination with a core population parameter Pcore of 1.979 when employing the VM database (Figure 1c). This model is in line with the original definition of the k-parameter by Coppens et al.43,44 “if k is larger than one, the valence density at r equals the free-atom valence density at a point more distant from the nucleus; thus, the atom is contracted relative to the free atom and vice versa”. Accordingly, the refinement results nicely reflect the contraction of the carbon atom in diamond. The residual density map (Figure 1c) shows a further improvement of the fit in comparison with the one in Figure 1b and can be classified as being virtually featureless (
0.02/ 0.03 eÅ
3). Simultaneous refinement of the carbon’s Uiso value yields a value close to zero (
0.000002(28) Å2). This result nicely underpins that our core k model (i) allows a reliable deconvolution of thermal smearing effects, (ii) provides chemically and physically relevant information about core expansion/ contraction phenomena, and (iii) is flexible enough to yield truly featureless ΔF(r) maps in case of our benchmark system diamond. In the next step of our analysis, we have tested whether our final model can be further improved by changing the shape of the normalized, nodeless radial density functions by employing (i) another scattering factor database and (ii) a different hybridization of the free reference atom (e.g., sp3 vs s2p2 hybridization45). A mere replacement of the VM by the SCM scattering factors does not yield any significant change of the models’ parameters (Supporting Information). Apparently, the theoretical method employed for the generation of both databases (VM: relativistic atomic ground state wave functions based on PBE/QZ4P/ ZORA calculations; SCM: nonlinear least-squares fittings of numerical relativistic atomic ground state wave functions by a linear combination of Slater-type functions) yield rather similar radial distribution curves for the carbon atom in its ground state. This situation changes significantly, as soon as radial functions of a sp3-hybridized carbon atom are employed. Now, the residual density peaks become even further reduced, both in the core and valence region (
0.01/0.02 eÅ
3). This is due to the fact that the F(r)sp3 profile is already well adopted to account for the reduction of the core density and the charge accumulation in the C
C valence bonding region (Figure 2b). As a consequence, the core k parameter of our final multipolar model (denoted W2K HC +kcore in Table 2) deviates less from unity (1.006) relative to the refinements based on the scattering factors of a free carbon atom (see above). Hence, three parameters seem to be essential for the systematic improvement of the standard HC model: (i) the shape of the radial functions, (ii) the adoption of the core population parameter (Pc), and finally, (iii) the core contraction/expansion parameter kcore (Table 2). We note that the systematic improvements by the W2K (HC +kcore) multipole model are not significantly reflected by a mere inspection of F(r) at the BCP (1.60 eÅ
3) because this density feature is already reasonably well recovered by the HC multipolar model (1.63 eÅ
3). Both F(r)BCP (HC and HC+kcore model) values are thus in excellent agreement with the benchmark value obtained by a direct analysis of the W2K and C06 wave functions (1.61/1.59 eÅ
3, respectively). The situation is, however, more evident in the case of the Laplacian, r2F(r), and its λ3

ARTICLE

Figure 3. Negative Laplacian, L(r) =
r2F(r) maps in the (110) plane of diamond based on the (a) direct analysis of the W2K wave function and (b) after smoothening and filtering by the W2K (HC+kcore) multipole model; (c) corresponding map based on the direct analysis of the C06 LCAO wave function (PBE/6-311G(3d) approximation). Note the good agreement of the bonded charge concentration (denoted r2F(r)BCC in Table 2) in the valence shell of the carbon atoms in all three models. Positive (solid) and negative (dashed) contour lines are drawn at 0, (2.0  10n, (4.0  10n, (8.0  10n eÅ
5, with n = (3, (2, (1, 0; F(r), L(r), and λ3 values at the BCP are specified in eÅ
3 and eÅ
5, respectively.

component, which characterizes the curvature of the charge density along the C
C bond path at the BCP (Table 2). In these instances, the core k model yields values (r2F(r) =
11.5 eÅ
5 and λ3 = 9.0 eÅ
5; Table 2), which appear to differ clearly from the benchmark values based on the W2K wave function (r2F(r) =
17.5 eÅ
5 and λ3 = 2.8 eÅ
5; Table 2). In that case, even the values of the standard model seem to be closer to the ones obtained from the W2K wave function analysis (Table 2). However, graphical analysis of the W2K L(r) =
r2F(r) map (Figure 3a) reveals that, in the special case of short C
C bonds, the muffin-tin region of both nuclei are separated by only a small plane wave region, a phenomenon that has been seen before in the case of transition metal carbides.4 Hence, the L(r) maps of the C
C BCP are affected by discontinuities arising at the borderline between the muffin tin and the plane wave region, where L(r) is not displaying a smooth 13065

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A

Figure 4. Wilson plot ln(F2o/F2sta)/(16π2) for the SPring-8 diamond powder diffraction data using the HC+kcore model at 300 K. Systematically absent reflections (except for the (222)-reflection) were excluded from the plot and the fit.

and continuously differentiable behavior. This effect can be (partially) compensated by adjusting the flexibility of the planewave basis set. Alternatively, comparison of the topological C
C characteristics of the W2K (HC+kcore) multipole model with those of the periodic PBE/6-311G(3d) LCAO calculations reveals that the core k model is also highly efficient in smoothening the L(r) maps of the W2K calculations. This is witnessed by an excellent agreement of all the salient parameters: (r2F(r) =
11.2/
11.5 eÅ
5 and λ3 = 9.3/9.0 eÅ
5 for the C06 wave function and the W2K (HC+kcore) multipole model, respectively; Table 2). Hence, the (HC+kcore) multipole model provides a highly robust method to recover charge density distributions from high-resolution theoretical structure factor data. In the next section, we will discuss whether core contraction/ expansion effects are relevant for experimental studies and which experimental resolution and data precision are needed to identify these tiny effects. We will also comment on earlier and alternative concepts to refine core density distributions.3,40,2,49,39 3.2. Experimental Core j Refinements: A Combined Rietveld/Multipolar Approach. In the next step of our analysis, we have tested the consequences of the core k refinements based on experimental data. Due to the simplicity of its structure, the highsite symmetry of the carbon atom (
43m) and low thermal motion that remains rather constant over a wide temperature range, the charge density distribution of the diamond has been investigated by numerous experimental charge density studies (see, for example, refs 19, 32, and 46). For our purposes, we have selected the high resolution powder diffraction data measured by Nishibori et al.32 at the third-generation synchrotron-radiation source SPring-8, BL02B2. In general, single-crystal data is considered to be more accurate than powder diffraction data due to inherent problems of extracting precise structure factors, especially in cases of severe peak overlap. However, the latter problem is not significant for high-symmetric compounds displaying small lattice parameters. On contrast, in cases of high crystal perfection (e.g., diamond) severe extinction effects might prevent precise charge density studies based on single-crystal diffraction data. Also, the ease of absorption corrections and data reduction (no interframe scaling), as well as the independency of the refinement success from crystal twinning phenomena, might render powder diffraction data more attractive for future charge density

ARTICLE

Figure 5. Rietveld refinement of the diamond powder diffraction pattern measured at SPring-8 at T = 300 K (ref 32) employing the Spring8 (HC + kcore) model, as specified in Table 2. Observed (black circles), calculated (solid red), and difference pattern (solid black, bottom). The inset shows the accuracy of the fit of the “forbidden” (222) reflection.

studies. Indeed, Svendsen et al. demonstrated recently that charge density studies on powder diffraction data of diamond (at 100 K) can be convincingly performed using an iterative approach of cyclic Rietveld and multipolar refinements.19 In our study, we will follow a direct approach and perform simultaneous Rietveld and multipolar refinements using JANA2006,31 which also allows the refinement of independent expansion/contraction parameters of individual shells. For our experimental study, we used the powder diffraction data set of the diamond measured by Nishibori et al.32 at room temperature, providing a direct comparison with the best multipolar refinements based on a single-crystal study using the Pendell€osung method.46,50 To our surprise, the combined Rietveld and multipolar refinements were robust and converged quickly because no significant correlation was observed between the multipolar parameters, the profile fitting parameters and the background function employed. The resulting multipolar parameters of our best model are listed in Table 2 and agree well with the ones obtained by Spackman based on single crystal data using the Pendell€osung method.46 However, due to the significant lower resolution of the single crystal data versus the powder diffraction data (sin θ/λ = 0.79 vs 1.45 Å
1, respectively) the standard deviations obtained in the single crystal study are significantly larger relative to the powder diffraction study. This is especially true in the case of the two independent multipolar population parameters at the carbon atom which are salient for the reconstruction of the chemical deformation density (O2
= 0.34(19) vs 0.37(1) and H0 = 0.0(1) vs
0.18(3) for the single crystal and our powder diffraction study, respectively). Furthermore, the multipolar parameters obtained by the combined Rietveld/multipolar model are close to the corresponding values based on the W2K (HC + Kcore) model (O2
= 0.338; H0 =
0.108; Table 2)). Hence, a significant population of the hexadecapole function H0 is clearly supported by our experimental study and the earlier powder diffraction study by Svendsen et al. at 100 K (O2
= 0.37(1); H0 =
0.14(2)),19 while the results of the single crystal study suggested that “the population of the hexadecapole function is insignificant”.46 Despite the low resolution of the single crystal data Spackman reported an isotropic Uiso value of 0.00177(13) Å2, which 13066

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A

Figure 6. (a) Residual density, ΔF(r), map in the (110) plane of the diamond lattice based on the experimental synchrotron data (ref 32) after multipolar refinement at a resolution of sin θ/λ < 0.79 Å
1 and kcore = 1.006, employing the Spring8 (HC + kcore) model, as specified in Table 2; (b) ΔF(r) map after multipolar refinement using the same model but a larger resolution limit of sin θ/λ < 1.45 Å
1. Positive (solid) and negative (dashed) contour lines are drawn employing a step width of 0.01 eÅ
3).

appeared to be rather robust and independent of the flexibility of the multipolar model employed.46 This value is, however, in good agreement with lattice dynamic calculations including thermal diffuse scattering51 (Uiso = 0.00176 Å2) and our own refinement results (Uiso = 0.00176(5) Å2). Note that the latter value is smaller when core k refinement is omitted and the scattering factors of a free carbon atom (VM database) is used (0.00161(5) Å2; Table 2). This finding is in excellent agreement with the theoretical predictions (vide supra). Further support is provided by a comparison of the observed structure factors Fobs with the theoretical static Fsta(W2K) values in terms of a simple Wilson plot.19 Applying a linear least-squares fit to the function ln(F2o/F2sta) = ln(k)
16π2Uiso(sin θ/λ)2 yields an Uiso value of 0.00175 Å2, in excellent agreement with our multipole refinements (Figure 4). Note that the Wilson plot shows the expected linear behavior over the full resolution range of the experimental data (sin θ/λ < 1.45 Å
1), which underpins the high quality of the SPring-8 powder diffraction data at 300 K.52 The excellent fitting of our final Rietveld/multipolar model is depicted in Figure 5, yielding flat and featureless residual density maps at full data resolution (sin θ/λ = 1.45 Å
1), which are shown in Figure 6. As expected from the comparison of Figure 1b and c, the differences in the experimental residual density maps are far too small to allow a simultaneous refinement of the Pc, kcore, and the Uiso parameter in the resolution limit of the experimental powder diffraction data (sin θ/λ = 1.45 Å
1). Hence, the kcore as well as the Pc parameter were restrained to the optimal value obtained from the refinements of the theoretical Fsta (W2K) data (Table 2). However, the Fsta (W2K) refinements showed that a simultaneous refinement of both critical parameters, kcore and Uiso, is possible (correlation coefficient = 0.995) already at the modest resolution of sin θ/λ = 1.45 Å
1. Pushing the resolution limit clearly beyond the limiting sphere of molybdenum radiation (sin θ/λ < 1.40 Å
1) enhances the ease of the kcore and Uiso deconvolution. This suggests that the phenomenon of core contraction/expansion in light atom compounds might be traceable by experimental charge density studies if the limits of data precision and resolutions can be pushed further in future studies. In the previous section, we demonstrated on the basis of theoretical Fdyn data that the thermal displacement parameters derived from standard multipolar models will be systematically wrong if the individual expansion/extraction coefficient of the inner shells are not considered in the refinements of high order data. However, despite the high data quality employed in our experimental charge density study on diamond we conclude that

ARTICLE

Figure 7. Residual density, ΔF(r), map in the (110) plane of R-silicon based on theoretical static structure factors Fsta (W2K) after standard multipolar refinement (denoted as HC model in Table 3); (b) ΔF(r) map after refinement of an additional isotropic thermal parameter (HC +Uiso); (c) ΔF(r) map after refinement of two additional expansion/ contraction parameters for the K and L shell (denoted as HC+kL+kK); two dominant ΔF(r) features opposite to the Si
Si bonds are marked by arrows; (d) ΔF(r) map after additional multipolar refinements of the aspherical density contributions of the K and L shell (denoted as HC +kL+kK+Core-MP). Positive (solid) and negative (dashed) contour lines are drawn employing a step width of 0.01 eÅ
3 (a
c) and 0.005 eÅ
3 (d); the resolution employed was sin θ/λ < 1.8 Å
1.

the systematic error in the thermal displacement parameter due to the omission of a core k refinement might not be significant in the range of the experimental errors (e.g., Uiso = 0.00161(5) and 0.00176(5) in the HC and (HC + kcore) models of diamond, respectively). For the time being, it is therefore safe to assume that there is no necessity to perform core k refinements (i) as long as the experimental resolution remains in the range of the limiting sphere of molybdenum radiation (sin θ/λ < 1.40 Å
1), (ii) the experiment is not performed at ultralow temperatures ( 1.8 Å
1) is unfortunately not easily available. However, the advent of more and more powerful X-ray sources might render experimental core polarization studies feasible in the near future. Hence, the statement from Bentley and Stewart in 1974 “the core deformation scattering factors are too small to be measured by current X-ray methods”56 should now become a realistic challenge for ambitious experimentalists. ’ ASSOCIATED CONTENT

bS

Supporting Information. Experimental and theoretical structure factors of diamond, multipolar refinements, and AIM

ARTICLE

analyses employing static theoretical structure factors (Fsta), comparison of multipolar refinements based on LCAO calculations using Crystal06, LAPW calculations employing Wien2K and PAW calculations, isocontour map of L(r) of the electron density distribution of one of the four Si
Si bonding NBOs in Si(SiH3)4. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]; bo@qchem. au.dk; [email protected].

’ ACKNOWLEDGMENT This work was supported by the DFG (SPP1178 and SCHE 487/12-2). ’ REFERENCES (1) Stewart, R. F. Isr. J. Chem. 1977, 16, 124–131. (2) Hansen, N. K.; Coppens, P. Acta Crystallogr. 1978, A34, 909–921. (3) Pillet, S.; Souhassou, M.; Lecomte, C.; Schwarz, K.; Blaha, P.; Rerat, M.; Lichanot, A.; Roversi, P. Acta Crystallogr. 2001, A57, 290–303. (4) Rohrmoser, B.; Eickerling, G.; Presnitz, M.; Scherer, W.; Eyert, V.; Hoffmann, R. D.; Rodewald, U. C.; Vogt, C.; P€ ottgen, R. J. Am. Chem. Soc. 2007, 129, 9356–9365. (5) Scherer, W.; Hauf, C.; Presnitz, M.; Scheidt, E. W.; Eickerling, G.; Eyert, V.; Hoffmann, R. D.; Rodewald, U. C.; Hammerschmidt, A.; Vogt, C.; P€ ottgen, R. Angew. Chem., Int. Ed. 2010, 49, 1578–1582. (6) Roquette, P.; Maronna, A.; Peters, A.; Kaifer, E.; Himmel, H. J.; Hauf, C.; Herz, V.; Scheidt, E. W.; Scherer, W. Chem.—Eur. J. 2010, 16, 1336–1350. (7) Poulsen, R. D.; Bentien, A.; Chevalier, M.; Iversen, B. B. J. Am. Chem. Soc. 2005, 127, 9156–9166. (8) Scherer, W.; Herz, V.; Bruck, A.; Hauf, C.; Reiner, F.; Altmannshofer, S.; Leusser, D.; Stalke, D. Angew. Chem., Int. Ed. 2011, 50, 2845–2849. (9) Scherer, W.; Wolstenholme, D. J.; Herz, V.; Eickerling, G.; Bruck, A.; Benndorf, P.; Roesky, P. W. Angew. Chem., Int. Ed. 2010, 49, 2242–2246. (10) McGrady, G. S.; Sirsch, P.; Chatterton, N. P.; Ostermann, A.; Gatti, C.; Altmannshofer, S.; Herz, V.; Eickerling, G.; Scherer, W. Inorg. Chem. 2009, 48, 1588–1598. (11) Volkov, A.; Abramov, Y.; Coppens, P.; Gatti, C. Acta Crystallogr. 2000, A56, 332–339. (12) Volkov, A.; Gatti, C.; Abramov, Y.; Coppens, P. Acta Crystallogr. 2000, A56, 252–258. (13) Volkov, A.; Coppens, P. Acta Crystallogr. 2001, A57, 395–405. (14) Figgis, B. N.; Iversen, B. B.; Larsen, F. K.; Reynolds, P. A. Acta Crystallogr. 1993, B49, 794–806. (15) We note, however, that an earlier charge density study of (ND4)2Cu(SO4)2 3 6D2O (ref 14), employing a “valence-orbital population model”, documented the need for flexible radial functions of doubleζ quality in the case of the deuterium atoms. (16) Koritsanszky, T.; Volkov, A. Chem. Phys. Lett. 2004, 385, 431–434. (17) Koritsanszky, T.; Volkov, A.; Chodkiewics, M. Struct. Bonding (Berlin, Ger.) 2011; DOI 10.1007/430_2010_32. (18) Blaha, P.; Schwarz, K.; Madsen, G.; Kvasnicka, D.; Luitz, J. WIEN2k 10.1, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties; Technische Universit€at Wien: Vienna, Austria, 2010. (19) Svendsen, H.; Overgaard, J.; Busselez, R.; Arnaud, B.; Rabiller, P.; Kurita, A.; Nishibori, E.; Sakata, M.; Takata, M.; Iversen, B. B. Acta Crystallogr. 2010, A66, 458–469. 13070

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071

The Journal of Physical Chemistry A (20) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. (21) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N.; Bush, J.; D’Arco, P.; Llunell, M. CRYSTAL06 User’s Manual; University Torino: Torino, 2006. (22) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724–728. (23) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257–2261. (24) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–654. (25) Mclean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639–5648. (26) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. J. Comput. Chem. 2001, 22, 931–967. (27) ADF2009.01, SCM, Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands, http://www.scm.com. (28) Guerra, C. Fonseca; Snijders, J. G.; Velde, G. te; Baerends, E. J. Theor. Chim. Acta 1998, 99, 391–403. (29) Volkov, A.; Koritsanszky, T.; Chodkiewicz, M.; King, H. F. J. Comput. Chem. 2009, 30, 1379–1391. (30) Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Bohmann, J. A.; Morales, C. M.; Weinhold, F. NBO 5.0, Theoretical Chemistry Institute, University of Wisconsin: Wisconsin, 2001. (31) Petricek, V.; Dusek, M.; Palatinus, L. Jana2006, The crystallographic computing system; Institute of Physics: Praha, Czech Republic, 2006. (32) Nishibori, E.; Sunaoshi, E.; Yoshida, A.; Aoyagi, S.; Kato, K.; Takata, M.; Sakata, M. Acta Crystallogr. 2007, A63, 43–52. (33) Volkov, A.; Macchi, P.; Farrugia, L. J.; Gatti, C.; Mallinson, P.; Richter, R.; Koritsanszky, T. XD2006, version 5.42, A computer program for multipole refinement, topological analysis of charge densities, and evaluation of intermolecular energies from experimental or theoretical structure factors; 2006. (34) Su, Z. W.; Coppens, P. Acta Crystallogr. 1998, A54, 646–652. (35) Macchi, P.; Coppens, P. Acta Crystallogr. 2001, A57, 656–662. (36) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953–17979. (37) Note that the artifical flatness of radial F(r) profiles in the cusp region is due to the usage of Gaussian-type basis sets of merely triple-ζ quality. (38) Bader, R. F. W.; Essen, H. J. Chem. Phys. 1984, 80, 1943–1960. (39) Hansen, N. K.; Schneider, J. R.; Yelon, W. B.; Pearson, W. H. Acta Crystallogr. 1987, A43, 763–769. (40) Chandler, G. S.; Spackman, M. A. Acta Crystallogr. 1982, A38, 225–239. (41) Bentley, J. J.; Stewart, R. F. Acta Crystallogr. 1976, A32, 910–914. (42) Coppens, P. Isr. J. Chem. 1977, 16, 159–162. (43) Griffin, J. F.; Coppens, P. J. Am. Chem. Soc. 1975, 97, 3496–3505. (44) Coppens, P.; Gururow, T. N.; Leung, P.; Stevens, E. D.; Becker, P. J.; Yang, Y. W. Acta Crystallogr. 1979, B35, 63–72. (45) In the case of radial functions derived from an sp3-hybridized carbon atom, one should, however, be aware that the original definition of the k parameter becomes somewhat blurred because in that case the reference “free” atom is no longer in its ground state. (46) Spackman, M. A. Acta Crystallogr. 1991, A47, 420–427. (47) Clementi, E.; Roetti, C. At. Data Nucl. Data Tables 1974, 14, 177–478. (48) Clementi, E.; Raimondi, D. L. J. Chem. Phys. 1963, 38, 2686–2689. (49) Hansen, N. K.; Schneider, J. R.; Larsen, F. K. Phys. Rev. B 1984, 29, 917–926. (50) Takama, T.; Tsuchiya, K.; Kobayashi, K.; Sato, S. Acta Crystallogr. 1990, A46, 514–517. (51) Stewart, R. F. Acta Crystallogr. 1973, A29, 602–605. (52) We note that the larger Uiso value determined at 100 K (0.00188 Å2; ref 19) relative to the one at 300 K (0.00175 Å2; this

ARTICLE

study) is most likely due to the fact that the 100 K data shows deviations from linearity in the Wilson plots at higher diffraction angles. (53) Each of the four bonding NBOs features a pronounced BCCM in the respective L(r)NBO maps and two additional maxima in the Lshell. The L(r) maximum cis to the Si
Si bond is more pronounced than the maximum in the trans position to the bond (for the L(r) map, see Supporting Information). (54) Destro, R.; Roversi, P.; Barzaghi, M.; Marsh, R. E. J. Phys. Chem. A 2000, 104, 1047–1054. (55) Zhurov, V. V.; Zhurova, E. A.; Stash, A. I.; Pinkerton, A. A. Acta Crystallogr. 2011, A67, 160–173. (56) Bentley, J.; Stewart, R. F. Acta Crystallogr. 1974, A30, 60–67.

13071

dx.doi.org/10.1021/jp2050405 |J. Phys. Chem. A 2011, 115, 13061–13071