International
Tables for
Crystallography
Volume C
Mathematical, physical and chemical tables
Edited by E. Prince

International Tables for Crystallography (2006). Vol. C, ch. 8.7, pp. 714-725

Section 8.7.3. Charge densities

P. Coppens,a Z. Sub and P. J. Beckerc

a732 NSM Building, Department of Chemistry, State University of New York at Buffalo, Buffalo, NY 14260-3000, USA,bDigital Equipment Co., 129 Parker Street, PKO1/C22, Maynard, MA 01754-2122, USA, and cEcole Centrale Paris, Centre de Recherche, Grand Voie des Vignes, F-92295 Châtenay Malabry CEDEX, France

8.7.3. Charge densities

| top | pdf |

8.7.3.1. Introduction

| top | pdf |

The charge density is related to the elastic X-ray scattering amplitude F(S) by the expression [ \left \langle \rho ({\bf r})\right \rangle =\textstyle\int F({\bf S})\exp (-2\pi i{\bf S\cdot {\bf r}}){\,{\rm d}}{\bf S},\eqno (8.7.3.1)]or, for scattering by a periodic lattice, [ \left \langle \rho ({\bf r})\right \rangle = {1\over V}\sum_{{\bf h}}F({\bf h})\exp (-2\pi i{\bf h\cdot {\bf r}}). ]As F(h) is in general complex, the Fourier transform (8.7.3.1)[link] requires calculation of the phases from a model for the charge distribution. In the centrosymmetric case, the free-atom model is in general adequate for calculation of the signs of F(h). However, for non-centrosymmetric structures in which phases are continuously variable, it is necessary to incorporate deviations from the free-atom density in the model to obtain estimates of the experimental phases.

Since the total density is dominated by the core distribution, differences between the total density ρ(r) and reference densities are important. The reference densities represent hypothetical states without chemical bonding or with only partial chemical bonding. Deviations from spherical, free-atom symmetry are obtained when the reference state is the promolecule, the superposition of free-space spherical atoms centred at the nuclear positions. This difference function is referred to as the deformation density (or standard deformation density) [ \Delta\rho ({\bf r})=\rho ({\bf r}) -\rho _p({\bf r}), \eqno (8.7.3.2)]where [\rho _p({\bf r})] is the promolecule density. In analogy to (8.7.3.1)[link], the deformation density may be obtained from [ \left \langle \Delta\rho \left ({\bf r}\right) \right \rangle = {1\over V}\sum \left [F_{{\rm obs}}\left ({\bf h}\right) -F_{{\rm calc,\,free\ atom}}({\bf h})\right] \exp \left (-2\pi i{\bf h\cdot r}\right), \eqno (8.7.3.3)]where [F_{{\rm obs}}] and [F_{{\rm calc}}] are in general complex.

Several other different density functions analogous to (8.7.3.3)[link], summarized in Table 8.7.3.1[link], may be defined. Particularly useful for the analysis of effects of chemical bonding is the fragment deformation density, in which a chemical fragment is subtracted from the total density of a molecule. The fragment density is calculated theoretically and thermally smeared before subtraction from an experimental density. `Prepared' atoms rather than spherical atoms may be used as a reference state to emphasize the electron-density shift due to covalent bond formation.

Table 8.7.3.1| top | pdf |
Definition of difference density functions

[\Delta \rho = {2 \over V} \left\{\sum\limits_0^{1/2} (A_1 - A_2)\cos2\pi{\bf h} \cdot {\bf r} + \sum\limits_0^{1/2} (B_1 -B_2)\sin 2\pi {\bf h} \cdot {\bf r} \right\}] with F = A + iB.

(a) Residual map A1, B1 from observations calculated with model phases.
A2, B2 from refinement model.
(b) X − X deformation map A1, B1 from observation, with model phases.
A2, B2 from high-order refinement, free-atom model, or other reference state.
(c) X − N deformation map [as (b) but] A2, B2 calculated with neutron parameters.
(d) X − (X + N) deformation map [as (b) but] A2, B2 calculated with parameters from joint refinement of X-ray neutron data.
(e) X − X, X − N, X − (X + N) valence map As (b), (c), (d) with A1, B1 calculated with core-electron contribution only.
(f) Dynamic model map A1, B1 from model.
A2, B2 with parameters from model refinement and free-atom functions.
(g) Static model map ρmodel − ρfree atom, where ρmodel is sum of static model density functions.

8.7.3.2. Modelling of the charge density

| top | pdf |

The electron density ρ(r) in the structure-factor expression [ F_{\rm calc}\left ({\bf h}\right)\textstyle \int\limits _{\rm unit\ cell} \rho \left ({\bf r}\right) \exp \left (2\pi i{\bf h\cdot r}\right) {\,{\rm d}}{\bf r} \eqno (8.7.3.4a) ]can be approximated by a sum of non-normalized density functions [g_{i}({\bf r})] with scattering factor [f_{i}({\bf h})] centred at [{\bf r}_{i},] [ \rho({\bf r}) =\textstyle\sum\limits_{i}g_{i}({\bf r}) *\delta ({\bf r}-{\bf r}_{i}). \eqno (8.7.3.5)]Substitution in (8.7.3.4a[link]) gives [ F({\bf h}) =\textstyle\sum\limits_{i}\, f_{i}({\bf h}) \exp(2\pi i{\bf h\cdot r}_{i}). \eqno (8.7.3.4b)]

When [g_{i}({\bf r})] is the spherically averaged, free-atom density, (8.7.3.4b[link]) represents the free-atom model. A distinction is often made between atom-centred models, in which all functions g(r) are centred at the nuclear positions, and models in which additional functions are centred at other locations, such as in bonds or lone-pair regions.

A simple, atom-centred model with spherical functions g(r) is defined by [ \rho _{{\rm atom}}({\bf r})=P_{{\rm core}}\rho _{{\rm core}}({\bf r})+\kappa ^3P_{{\rm valence}}\rho _{{\rm valence}}(\kappa {\bf r}).\eqno (8.7.3.6)]This `kappa model' allows for charge transfer between atomic valence shells through the population parameter [P_{{\rm valence}}], and for a change in nuclear screening with electron population, through the parameter κ, which represents an expansion [(\kappa \, \lt \, 1)], or a contraction [(\kappa \, \gt \, 1)] of the radial density distribution.

The atom-centred, spherical harmonic expansion of the electronic part of the charge distribution is defined by [\eqalignno{ \rho _{{\rm atom}}({\bf r}) &=P_{{\rm c}}\rho _{{\rm core}}({\bf r})+P_{{\rm v}}\kappa ^3\rho _{{\rm valence}}(\kappa r) \cr &\quad +\textstyle\sum \limits _{l=0}^{l(\max)}\kappa ^{\prime\,3}R_l(\kappa ^{\prime}\zeta r)\textstyle\sum \limits _{m=0}^l \, \textstyle\sum\limits _pP_{lmp}d_{lmp}(\theta, \varphi), & (8.7.3.7)}]where p = ± when m is larger than 0, and [R_l(\kappa ^{\prime }\zeta r) ] is a radial function.

The real spherical harmonic functions [{d}_{lmp}] and their Fourier transforms have been described in International Tables for Crystallography, Volume B, Chapter 1.2[link] (Coppens, 2001[link]). They differ from the functions [{y}_{lmp}] by the normalization condition, defined as [\int | d_{lmp} | {\,{\rm d}}\Omega=2-\delta _{l0}]. The real spherical harmonic functions are often referred to as multipoles, since each represents the component of the charge distribution ρ(r) that gives a non-zero contribution to the integral for the electrostatic multipole moment [{q}_{lmp},] [ {q}_{lmp}=-\textstyle\int \rho _{{\rm atom}}({\bf r}) {r}^l{c}_{lmp}{\,{\rm d}}{\bf r},\eqno (8.7.3.8)]where the functions [{c}_{lmp}] are the Cartesian representations of the real spherical harmonics (Coppens, 2001[link]).

More general models include non-atom centred functions. If the wavefunction ψ in (8.7.2.4)[link] is an antisymmetrized product of molecular orbitals ψi, expressed in terms of a linear combination of atomic orbitals [\chi _{\mu},\ \psi _{i}=\sum _\mu {c}_{{i}\mu }\chi _\mu ] (LCAO formalism), the integration (8.7.2.4)[link] leads to [ \rho ({\bf r}) =\textstyle\sum \limits _\mu \textstyle\sum \limits _\nu P_{\mu \nu }\chi _\mu ({\bf r}-{\bf R}_\mu) \chi _\nu ({\bf r}-{\bf R}_\nu), \eqno (8.7.3.9)]with Rμ and Rν defining the centres of [\chi _\mu] and [\chi _\nu], respectively, [P_{\mu \nu }=\sum _{i}n_{i}{c}_{i\mu }{c}_{{i}\nu}], and the sum is over all molecular orbitals with occupancy [n_{i}]. Expression (8.7.3.9)[link] contains products of atomic orbitals, which may have significant values for orbitals centred on adjacent atoms. In the `charge-cloud' model (Hellner, 1977[link]), these products are approximated by bond-centred, Gaussian-shaped density functions. Such functions can often be projected efficiently into the one-centre terms of the spherical harmonic multipole model, so that large correlations occur if both spherical harmonics and bond-centred functions are adjusted independently in a least-squares refinement.

According to (8.7.2.4)[link] and (8.7.3.9)[link], the population of the two-centre terms is related to the one-centre occupancies. A molecular-orbital based model, which implicitly incorporates such relations, has been used to describe local bonding between transition-metal and ligand atoms (Becker & Coppens, 1985[link]).

8.7.3.3. Physical constraints

| top | pdf |

There are several physical constraints that an electron-density model must satisfy. With the exception of the electroneutrality constraint, they depend strongly on the electron density close to the nucleus, which is poorly determined by the diffraction experiment.

8.7.3.3.1. Electroneutrality constraint

| top | pdf |

Since a crystal is neutral, the total electron population must equal the sum of the nuclear charges of the constituent atoms. A constraint procedure for linear least squares that does not increase the size of the least-squares matrix has been described by Hamilton (1964[link]). If the starting point is a neutral crystal, the constraint equation becomes [\textstyle\sum \Delta P_{i}S_{i}=0,\eqno (8.7.3.10)]where [S_{i}=\textstyle\int g({\bf r})\,{\rm d}{\bf r}], g being a general density function, and the [\Delta P_{i}] are the shifts in the population parameters. For the multipole model, only the monopolar functions integrate to a non-zero value. For normalized monopole functions, this gives [ \textstyle\sum \limits _{\rm monopoles} \Delta P_{i}=0.\eqno (8.7.3.11)]If the shifts without constraints are given by the vector y and the constrained shifts by yc, the Hamilton constraint is expressed as [ {\bf y}_{{c}}^T={\bf y}^T-{\bf y}^T{\bi Q}^T({\bi QA}^{-1}{\bi Q}^T)^{-1}{\bi QA}^{-1},\eqno (8.7.3.12)]where the superscript T indicates transposition, A is the least-squares matrix of the products of derivatives, and Q is a row vector of the values of [S_{i}] for elements representing density functions and zeros otherwise.

Expression (8.7.3.12)[link] cannot be applied if the unconstrained refinement corresponds to a singular matrix. This would be the case if all population parameters, including those of the core functions, were to be refined together with the scale factor. In this case, a new set of independent parameters must be defined, as described in Chapter 8.1[link] on least-squares refinements. Alternatively, one may set the scale factor to one and rescale the population parameters to neutrality after completion of the refinement. This will in general give a non-integral electron population for the core functions. The proper interpretation of such a result is that a core-like function is an appropriate component of the density basis set representing the valence electrons.

8.7.3.3.2. Cusp constraint

| top | pdf |

The electron density at a nucleus i with nuclear charge [Z_{i}] must satisfy the electron–nuclear cusp condition given by [ \lim _{{\bf r}_{i}\rightarrow 0}\left ({\partial \over \partial {\bf r}}_{ i}+2Z_{i}\right) \rho _{0{i}} ({\bf r}_{i}) =0,\eqno (8.7.3.13)]where [\rho _{0{i}}({\bf r}_{i})=(1/4\pi)\textstyle\int\rho({\bf r})\,{\rm d}\Omega_i] is the spherical component of the expansion of the density around nucleus i.

Only 1s-type functions have non-zero electron density at the nucleus and contribute to (8.7.3.13)[link]. For the hydrogen-like atom or ion described by a single exponent radial function R(r) = N exp (−ζr), (8.7.3.13)[link] gives [\zeta =2Z/a_0], where Z is the nuclear charge, and [a_0] is the Bohr unit. Thus, a modification of ζ for 1s functions, as implied by (8.7.3.6)[link] and (8.7.3.7)[link] if applied to H atoms, leads to a violation of the cusp constraint. In practice, the electron density at the nucleus is not determined by a limited resolution diffraction experiment; the single exponent function R(r) is fitted to the electron density away from, rather than at the nucleus.

8.7.3.3.3. Radial constraint

| top | pdf |

Poisson's electrostatic equation gives a relation between the gradient of the electric field ∇2Φ(r) and the electron density at r. [ \nabla ^2\Phi ({\bf r}) =-4\pi \rho ({\bf r}). \eqno (8.7.3.14)]As noted by Stewart (1977[link]), this equation imposes a constraint on the radial functions R(r). For [R_l(r)=N_l r^{n(l)}\exp (-\zeta _lr)], the condition [n(l)\geq l] must be obeyed for [R_l r^{-l}] to be finite at r = 0, which satisfies the requirement of the non-divergence of the electric field ∇V, its gradient ∇2V, the gradient of the field gradient ∇3V, etc.

8.7.3.3.4. Hellmann–Feynman constraint

| top | pdf |

According to the electrostatic Hellmann–Feynman theorem, which follows from the Born–Oppenheimer approximation and the condition that the forces on the nuclei must vanish when the nuclear configuration is in equilibrium, the nuclear repulsions are balanced by the electron–nucleus attractions (Levine, 1983[link]). The balance of forces is often achieved by a very sharp polarization of the electron density very close to the nuclei (Hirshfeld & Rzotkiewicz, 1974[link]), which may be represented in the X-ray model by the introduction of dipolar functions with large values of ζ. The Hellmann–Feynman constraint offers the possibility for obtaining information on such functions even though they may contribute only marginally to the observed X-ray scattering (Hirshfeld, 1984[link]).

As the Hellmann–Feynman constraint applies to the static density, its application presumes a proper deconvolution of the thermal motion and the electron density in the scattering model.

8.7.3.4. Electrostatic moments and the potential due to a charge distribution

| top | pdf |

8.7.3.4.1. Moments of a charge distribution1

| top | pdf |

Use of the expectation value expression [ \langle {O}\rangle =\textstyle\int \widehat {{O}}\rho({\bf r}) {\,{\rm d}}{\bf r},\eqno (8.7.3.15)]with the operator [\widehat {{O}}=\widehat {\gamma }r_{\alpha _1}r_{\alpha _2}r_{\alpha _3}\ldots r_{\alpha _l}=r_{\alpha _1}r_{\alpha _2}r_{\alpha _3}\ldots r_{\alpha _l}] gives for the electrostatic moments of a charge distribution ρ(r) [\mu _{\alpha _1\alpha _2\alpha _3\ldots \alpha _l}=\textstyle\int \rho({\bf r}) r_{\alpha _1}r_{\alpha _2}r_{\alpha _3}\ldots r_{\alpha _l}{\,{\rm d}}{\bf r},\eqno (8.7.3.16)]in which the [r_\alpha ] are the three components of the vector ri = 1, 2, 3), and the integral is over the complete volume of the distribution.

For l = 0, (8.7.3.16)[link] represents the integral over the charge distribution, which is the total charge, a scalar function described as the monopole. The higher moments are, in ascending order of l, the dipole, a vector, the quadrupole, a second-rank tensor, and the octupole, a third-rank tensor. Successively higher moments are named the hexadecapole (l = 4), the tricontadipole (l = 5), and the hexacontatetrapole (l = 6). An alternative, traceless, definition is often used for moments with [l\geq 2]. In the traceless definition, the quadrupole moment, [\Theta _{\alpha \beta }], is given by [ \Theta_{\alpha \beta }= \textstyle{1\over2}\int \rho ({\bf r}) \left [3r_\alpha r_\beta -r^2\delta _{\alpha \beta }\right] {\,{\rm d}}r,\eqno (8.7.3.17)]where δαβ is the Kronecker delta function. The term [\textstyle\int\rho({\bf r})r^2\,{\rm d}r], which is subtracted from the diagonal elements of the tensor, corresponds to the spherically averaged second moment of the distribution.

Expression (8.7.3.17)[link] is a special case of the following general expression for the lth-rank traceless tensor elements. [ { M}_{\alpha _1\alpha _2\ldots \alpha _l}^{\left (l\right) }= {\left (-1\right) ^l \over l!} \int \rho ({\bf r}) r^{2l+1} \displaystyle{\partial ^l \over \partial r_{\alpha _1}\partial r_{\alpha _2}\ldots \partial r_{\alpha _l}}\left (1\over r\right) {\,{\rm d}}{\bf r}.\eqno (8.7.3.18)]

Though the traceless moments can be derived from the unabridged moments, the converse is not the case because the information on the spherically averaged moments is no longer present in the traceless moments. The general relations between the traceless moments and the unabridged moments follow from (8.7.3.18)[link]. For the quadrupole moments, we obtain with (8.7.3.17)[link] [\eqalignno{ \Theta_{xx} &= \textstyle{3\over2}\mu _{xx}- {1\over 2} \left (\mu _{xx}+\mu _{yy}+\mu _{zz}\right) \cr &=\mu _{xx}- \textstyle{1\over2} \left (\mu _{yy}+\mu _{zz}\right), }]and [ \Theta _{xy}= \textstyle{3\over2}\mu _{xy}.\eqno (8.7.3.19)]Expressions for the other elements are obtained by simple permutation of the indices.

For a site of point symmetry 1, the electrostatic moment [\mu _{\alpha _1\alpha _2\alpha _3\ldots \alpha _l}] of order l has (l + 1)(l + 2)/2 unique elements. In the traceless definition, not all elements are independent. Because the trace of the tensor has been set to zero, only 2l + 1 independent components remain. For the quadrupole there are 5 independent components of the form (8.7.3.19)[link].

In a different form, the traceless moment operators can be written as the Cartesian spherical harmonics [{c}_{lmp}] (ITB, 2001[link]) multiplied by [r^l], which defines the spherical electrostatic moments [ \Theta _{lmp}=\textstyle\int \rho({\bf r}) {c}_{lmp}r^l{\,{\rm d}}r. \eqno (8.7.3.20)]

The expressions for [{c}_{lmp}] are listed in Volume B of International Tables for Crystallography (ITB, 2001[link]); for the l = 2 moment, the [{c}_{lmp}] have the well known form 3z2 − 1, xz, yz, (x2y2)/2, and xy, where x, y and z are the components of a unit vector from the origin to the point being described. The spherical electrostatic moments have (2l + 1) components, which equals the number of independent components in the traceless definition (8.7.3.18)[link], as it should. The linear relationships are [\eqalignno{ \Theta _{zz} &=(1/2) \Theta _{20}, \cr \Theta _{xx} &=(1/2) [3\Theta _{22+}-(1/2) \Theta _{20}] , \cr \Theta _{yy} &=(1/2) [-3\Theta _{22+}-(1/2) \Theta _{20}] , \ \cr \Theta _{xz} &=(3/2) \Theta _{21+}, \cr \Theta _{yz} &=(3/2) \Theta _{21-}, \cr \Theta _{xy} &=(3/2) \Theta _{22-}. & (8.7.3.21)}]

8.7.3.4.1.1. Moments as a function of the atomic multipole expansion

| top | pdf |

In the multipole model [expression (8.7.3.7)[link]], the charge density is a sum of atom-centred density functions, and the moments of a whole distribution can be written as a sum over the atomic moments plus a contribution due to the shift to a common origin. An atomic moment is obtained by integration over the charge distributions ρtotal,i(r) = ρnuclear,i − ρe,i of atom i, [ \mu _{\alpha _1\alpha _2\alpha _3\ldots \alpha _l}=\textstyle\int \rho _{{\rm total},{i}}({\bf r}) r_{\alpha _1}r_{\alpha_2}r_{\alpha _3}\ldots r_{\alpha _l}{\,{\rm d}}r,\eqno (8.7.3.22)]where the electronic part of the atomic charge distribution is defined by the multipole expansion [\eqalignno{ \rho _{{e,i}}({\bf r}) &={P}_{{i,c}}\rho _{{\rm core}} (r) +{P}_{{i,v}}\kappa _{i}^3\rho _{ i,{\rm valence}}(\kappa _{i}r) \cr &\quad+\textstyle\sum \limits _{l=0}^{l_{\max }}\kappa _{i}^{\prime 3}{R}_{{i},l}(\kappa _{i}^{\prime }r) \textstyle\sum \limits _{m=0}^l\, \textstyle\sum \limits _p{P}_{{i},lmp}d_{lmp}(\theta, \varphi), & (8.7.3.23)}]where p = ± when m [\gt \, 0], and [{R}_l(\kappa _{i}^{\prime }r)] is a radial function.

We get for the jth moment of the valence density [\eqalignno{ \mu ^j &=\mu _{\alpha _1\alpha _2\alpha _3\ldots \alpha _j} \cr &=-\int \left [{P}_{{i,v}}\kappa _{i}^3\rho _{i,{\rm valence}}(\kappa _{i}r)+\textstyle\sum \limits _{l=0}^{l_{\max }}\kappa _{{i}}^{\prime 3}{R}_{{i},l}(\kappa _{i}'r) \right. \cr &\quad \times \left. \textstyle\sum \limits _{m=0}^l\textstyle\sum \limits _p{P}_{{i},lmp}{d}_{lmp}(\theta, \varphi ) \right] r_{\alpha _1}r_{\alpha _2}\ldots r_{\alpha _j}{\,{\rm d}}{\bf r}, & (8.7.3.24)}]in which the minus sign arises because of the negative charge of the electrons.

We will use the symbol [\widehat {{O}}] for the moment operators. We get [ \mu ^j=-\kappa _{i}^{\prime 3}\int \widehat {{O}}_j\sum \limits _{l=1}^{l_{\max }}\left [\textstyle\sum \limits _{m=0}^l\textstyle\sum \limits _p{P}_{lmp}d_{lmp}{R}_l\right] {\,{\rm d}}{\bf r},\eqno (8.7.3.25)]where, as before, p = ±. The requirement that the integrand be totally symmetric means that only the dipolar terms in the multipole expansion contribute to the dipole moment. If we use the traceless definition of the higher moments, or the equivalent definition of the moments in terms of the spherical harmonic functions, only the quadrupolar terms of the multipole expansion will contribute to the quadrupole moment; more generally, in the traceless definition the lth-order multipoles are the sole contributors to the lth moments. In terms of the spherical moments, we get [ \Theta _{lmp}=-{P}_{lmp}\textstyle\int \widehat {{O}}_{lmp}\left [d_{lmp}{R}_l\right] {\,{\rm d}}{\bf r}.\eqno (8.7.3.26)]

Substitution with Rl = {(κ′ζ)n(l)+3/[n(l) +2] !}rn(l)exp(−ζr) and [\widehat {{O}}_{lmp}={c}_{lmp}r^l] and subsequent integration over r gives [ \Theta _{lmp}=-{P}_{lmp} {1\over(\kappa ^{\prime }\zeta ) ^l}\, {[n(l)+l+2] ! \over [n(l) +2] !}\, {1\over{D}_{lm}{M}_{lm}}\int {y}_{lmp}^2\sin \theta \,{\,{\rm d}}\theta {\,{\rm d}}\varphi, \eqno (8.7.3.27)]where the definitions [ d_{lmp}={L}_{lm}{c}_{lmp}=\left ({{L}_{lm} \over {M}_{lm}}\right) {y}_{lmp}\quad {\rm and}\quad {c}_{lmp}=\left ({1\over {M}_{lm}}\right) {y}_{lmp}\eqno (8.7.3.28)]have been used (ITB, 2001[link]). Since the [{y}_{lmp}] functions are wavefunction normalized, we obtain [ \Theta _{lmp}=-{P}_{lmp} {1\over(\kappa ^{\prime }\zeta ) ^l}\, { [n(l) +l+2] ! \over [n(l) +2] !}\, {{L}_{lm} \over \left ({M}_{lm}\right) ^2}.\eqno (8.7.3.29)]Application to dipolar terms with n(l) = 2, Llm = 1/π and Mlm = (3/4π)1/2 gives the x component of the atomic dipole moment as [ \mu _x=-\int {P}_{11+}{\,d}_{11+}{R}_1x{\,{\rm d}}{\bf r}=- \displaystyle{20 \over 3\kappa ^{\prime }\zeta }{P}_{11+}.\eqno (8.7.3.30)]For the atomic quadrupole moments in the spherical definition, we obtain directly, using n(l) = 2, l = 2 in (8.7.3.29)[link], [ \Theta _{20}=- {30 \over \left (\kappa ^{\prime }\zeta \right) ^2}\, {{L}_{20}\over \left ({M}_{20}\right) ^2}\,{P}_{20}=- {36\sqrt {3} \over \left (\kappa ^{\prime }\zeta \right) ^2}{P}_{20},\eqno (8.7.3.31)]and, for the other elements,[\Theta _{2mp}=-{30\over \left (\kappa ^{\prime }\zeta \right) ^2}\, {{L}_{2m}\over \left ({M}_{2m}\right) ^2}{P}_{2mp}=- {6\pi \over\left (\kappa ^{\prime }\zeta \right) ^2}{P}_{2mp}. \eqno (8.7.3.32)]

As the traceless quadrupole moments are linear combinations of the spherical quadrupole moments, the corresponding expressions follow directly from (8.7.3.31)[link], (8.7.3.32)[link] and (8.7.3.21)[link]. We obtain with n(2) = 2 [\eqalignno{ \Theta _{zz} &=- {18\sqrt {3}\over \left (\kappa ^{\prime }\zeta \right) ^2}{P}_{20}, \cr \Theta _{yy} &=+{9\over\left (\kappa ^{\prime }\zeta \right) ^2}\left (\sqrt {3}{P}_{20}+\pi {P}_{22+}\right), \cr \Theta _{xx} &= {9\over \left (\kappa ^{\prime }\zeta \right) ^2}\left (\sqrt {3}{P}_{20}-\pi {P}_{22+}\right),}]and [ \Theta _{xz}=-{9\pi\over{\left (\kappa ^{\prime }\zeta \right) ^2}}{P}_{21+}, \eqno (8.7.3.33)]and analogously for the other off-diagonal elements.

8.7.3.4.1.2. Molecular moments based on the deformation density

| top | pdf |

The moments derived from the total density ρ(r) and from the deformation density Δρ(r) are not identical. To illustrate the relation for the diagonal elements of the second-moment tensor, we rewrite the xx element as [\eqalignno{ \mu _{xx}\left (\rho _{{\rm total}}\right) &=\textstyle\int \rho x^2{\,{\rm d}}{\bf r} \cr &=\textstyle\int \rho _{{\rm promolecule}}x^2{\,{\rm d}}{\bf r}+\int \Delta\rho x^2{\,{\rm d}}{\bf r.} & (8.7.3.34)}]The promolecule is the sum over spherical atom densities, or [\eqalignno{ \textstyle\int \rho _{{\rm promolecule}}x^2{\,{\rm d}}{\bf r} &= \textstyle\int  \textstyle\sum \limits _{{i}}\rho _{{\rm spherical\ atom,}i}x^2{\,{\rm d}}{\bf r} \cr &= \textstyle\sum \limits _{i} \textstyle\int \rho _{{\rm spherical\ atom,}i}x^2{\,{\rm d}}{\bf r.} & (8.7.3.35)}]If Ri = (Xi, Yi, Zi) is the position vector for atom i, each single-atom contribution can be rewritten as [\eqalignno{ \mu _{{i},\,xx,\,{\rm spherical\ atom}} &= \textstyle\int \rho _{i{\rm,\, spherical\ atom}}x^2{\,{\rm d}}{\bf r} \cr &= \textstyle\int \rho _{i{\rm,\, spherical\ atom}} (x-X_{i}) ^2{\,{\rm d}}{\bf r} \cr &\quad+X_{i} \textstyle\int \rho _{i{\rm,\, spherical\ atom}}2(x-X_{{i}}) {\,{\rm d}}{\bf r} \cr &\quad+X_{i}^2 \textstyle\int \rho _{i{\rm,\, spherical\ atom}}{\,{\rm d}}{\bf r.}\ & (8.7.3.36)}]Since the last two integrals are proportional to the atomic dipole moment and its net charge, respectively, they will be zero for neutral spherical atoms. Substitution in (8.7.3.35)[link] gives, with [\langle(x-X) {_{i}^2}\rangle ={1\over3}\langle r_{{i}}^2\rangle], and [\langle r_{i}^2\rangle =\int \rho _{i}(r) r^2{\,{\rm d}}{\bf r},] [ \textstyle\int \rho _{{\rm promolecule}}x^2{\,{\rm d}}{\bf r}= \textstyle{1\over3}\textstyle\sum \limits _{{\rm atoms}}\langle r^2\rangle _{{\rm spherical\ atom}},\eqno (8.7.3.37)]and, by substitution in (8.7.3.34)[link], [ \mu _{xx}(\rho _{{\rm tot}}) =\mu _{xx} (\Delta\rho) + \textstyle{1\over3}\textstyle\sum \limits _{{\rm atoms}}\langle r^2\rangle _{{\rm spherical\ atom}},\ \eqno (8.7.3.38a) ]with [ \mu _{xx}(\Delta\rho) =\textstyle\sum \limits _{i}\left(\textstyle\int \Delta\rho _{i}x^2{\,{\rm d}}{\bf r}+2X_{i}\mu _{i}+X_{{i}}^2q_{i}\right), \eqno (8.7.3.38b) ]in which [\mu _{i}] and [q_{i}] are the atomic dipole moment and the charge on atom i, respectively.

The last term in (8.7.3.38a[link]) can be derived rapidly from analytical expressions for the atomic wavefunctions. Results for Hartree–Fock wavefunctions have been tabulated by Boyd (1977[link]). Since the off-diagonal elements of the second-moment tensor vanish for the spherical atom, the second term in (8.7.3.38a[link]) disappears, and the off-diagonal elements are identical for the total and deformation densities.

The relation between the second moments μαβ and the traceless moments [\Theta]αβ of the deformation density can be illustrated as follows. From (8.7.3.17)[link], we may write [ \Theta _{\alpha \beta }\left (\Delta\rho \right) = \textstyle{3\over2}\mu _{\alpha \beta }\left (\Delta\rho \right) - {1\over2}\delta _{\alpha \beta }\int \Delta\rho \,r^2{\,{\rm d}}{\bf r}.\eqno (8.7.3.39)]Only the spherical density terms contribute to the integral on the right. Assuming for the moment that the spherical deformation is represented by the valence-shell distortion (i.e. neglect of the second monopole in the aspherical atom expansion), we have, with density functions ρ normalized to 1, for each atom [ (\Delta\rho) _{{\rm spherical}} = \kappa ^3{P}_{{\rm valence}}\rho _{{\rm valence}}(\kappa r) -{P}_{{\rm valence}}^0\rho _{{\rm valence}}(r) \eqno (8.7.3.40)]and [\eqalignno{ \textstyle\int \Delta\rho \,r^2{\,{\rm d}}{\bf r} &=\textstyle\int \textstyle\sum \limits _{i}\left [\kappa _{i}^3{P}_{{\rm valence,}i}\rho _{{\rm valence,}i}\left (\kappa _{i}r\right) \right.\cr &\left. \quad -{P}_{{\rm valence,}i}^0\rho _{{\rm valence,}i}\left (r\right) \right] r^2{\,{\rm d}}{\bf r} \cr &=\textstyle\sum \limits _{i}\left ({P}_{{\rm valence,}i}/\kappa _{{i}}^2-{P}_{{\rm valence}}^0\right) \left \langle r_{i}^2\right \rangle _{{\rm spherical\ valence\ shell}} \cr &\quad+{R}_{i}^2\left ({P}_{{\rm valence,}i}-{P}_{{\rm valence,}i}^0\right), & (8.7.3.41)}]which, on substitution in (8.7.3.39)[link], gives the required relation.

8.7.3.4.1.3. The effect of an origin shift on the outer moments

| top | pdf |

In general, the multipole moments depend on the choice of origin. This can be seen as follows. Substitution of [{\bf r}_\alpha ^{\prime }={\bf r}_{\alpha} - {\bf R}_{\alpha}] in (8.7.3.16)[link] corresponds to a shift of origin by Rα, or X, Y, Z in the original coordinate system. In three dimensions, we get, for the first moment, the charge q, [ q^{\prime }=q,\eqno (8.7.3.42)]and for the transformed first and second moments [\eqalignno{ \mu _x^{\prime } &=\mu _x-q{\bf X}\semi \ \mu _y^{\prime }=\mu _y-q{\bf Y}\semi \ \mu _z^{\prime }=\mu _{z}-q{\bf Z}\semi \cr \mu _{\alpha \alpha }^{\prime } &=\mu _{\alpha \alpha }-2\mu _\alpha {R}_\alpha +q{R}_\alpha ^2\semi \cr \mu _{\alpha \beta }^{\prime } &=\mu _{\alpha \beta }-\mu _\alpha {R}_\beta -\mu _\beta {R}_\alpha +q{R}_\alpha {R}_\beta. & (8.7.3.43)}]

For the traceless quadrupole moments, the corresponding equations are obtained by substitution of [{\bf r}_\alpha ^{\prime }={\bf r}_{\alpha}-{\bf R}_{\alpha}] and r′ = rR into (8.7.3.17)[link], which gives [\eqalignno{ \Theta _{\alpha \beta }^{\prime } &=\Theta _{\alpha \beta }+ \textstyle{1\over2}\left (3{\bf R}_\alpha {\bf R}_\beta -{\bf R}^2\delta _{\alpha \beta }\right) q \cr &\quad - \textstyle{3\over2}\left ({\bf R}_\beta \mu _\alpha +{\bf R}_\alpha \mu _\beta \right) +\sum \limits _\gamma \left ({\bf R}_\gamma \mu _\gamma \right) \delta _{\alpha \beta }. & (8.7.3.44)}]Similar expressions for the higher moments are reported in the literature (Buckingham, 1970[link]).

We note that the first non-vanishing moment is origin-independent. Thus, the dipole moment of a neutral molecule, but not that of an ion, is independent of origin; the quadrupole moment of a molecule without charge and dipole moment is not dependent on the choice of origin and so on. The molecular electric moments are commonly reported with respect to the centre of mass.

8.7.3.4.1.4. Total moments as a sum over the pseudoatom moments

| top | pdf |

The moments of a molecule or of a molecular fragment are obtained from the sum over the atomic moments, plus a contribution due to the shift to a common origin for all but the monopoles. If individual atomic coordinate systems are used, as is common if chemical constraints are applied in the least-squares refinement, they must be rotated to have a common orientation. Expressions for coordinate system rotations have been given by Cromer, Larson & Stewart (1976[link]) and by Su & Coppens (1994a[link]).

The transformation to a common coordinate origin requires use of the origin-shift expressions (8.7.3.42)[link][link]–(8.7.3.44)[link], with, for an atom at [{\bf r}_{{i}}], [{\bf R}=-{\bf r}_i]. The first three moments summed over the atoms i located at [{\bf r}_{i}] become [\eqalignno{ q_{{\rm total}} &=\textstyle\sum q_{i}, &(8.7.3.45)\cr \mu _{{\rm total}} &=\textstyle\sum \limits _{i}\mu _{i}+\textstyle\sum {\bf r}_{{i}}q_{i}, &(8.7.3.46)}%fd8.7.3.46]and [ \mu _{\alpha \beta, {\rm total}}=\textstyle\sum \limits _{i}\big (\mu _{\alpha \beta {i}}+r_{\beta {i}}^{\prime }\mu _\alpha +r_{\alpha {i}}^{\prime }\mu _\beta +r_{\alpha {i}}^{\prime }r_{\beta {i}}^{\prime }q_{i}\big) \eqno (8.7.3.47)]with α, β = x, y, z; and expressions equivalent to (8.7.3.44)[link] for the traceless components [\Theta]αβ.

8.7.3.4.1.5. Electrostatic moments of a subvolume of space by Fourier summation

| top | pdf |

Expression (8.7.3.16)[link] for the outer moment of a distribution within a volume element [V_T] may be written as [ \mu _{\alpha _1\alpha _2\ldots \alpha _l}=\textstyle\int\limits _{V_T}\rho({\bf r}) \widehat {\gamma }_{\alpha _1\alpha _2\alpha _3\ldots \alpha _l}{\,{\rm d}}{\bf r}, \eqno (8.7.3.16)]with [\widehat {\gamma }_{\alpha _1\alpha _2\alpha _3\ldots \alpha _l}=r_{\alpha _1}r_{\alpha _2}r_{\alpha _3}\ldots r_{\alpha _l}], and integration over the volume [V_T].

Replacement of ρ(r) by the Fourier summation over the structure factors gives [\eqalignno{ \mu _{(V_T) }^l &= {1\over V}\int\limits _{V_T}\widehat {\gamma}_l\sum F({\bf h}) \exp (-2\pi i{\bf h\cdot r}) \,{\,{\rm d}}{\bf r} \cr &={1\over V}\sum F({\bf h}) \int\limits _{V_T}\widehat {\gamma }_l\exp(-2\pi i{\bf h\cdot r}) \,{\,{\rm d}}{\bf r},& (8.7.3.48)}]where [\widehat\gamma_l] is the product of l coordinates according to (8.7.3.16)[link], and μl represents the moment of the static distribution if the F(h) are the structure factors on an absolute scale after deconvolution of thermal motion. Otherwise, the moment of the thermally averaged density is obtained.

The integral [\int _{V_T}\widehat {\gamma }_l] [\exp (-2\pi i{\bf h}\cdot {\bf r}){\,{\rm d}}{\bf r}] is defined as the shape transform S of the volume [V_T.] [ \mu ^l(V_T) = {1\over V}\sum F({\bf h}) S_{V_T}(\widehat {\gamma }_l,{\bf h}). ]For regularly shaped volumes, the integral can be evaluated analytically. A volume of complex shape may be subdivided into integrable subvolumes such as parallelepipeds. By choosing the subvolumes sufficiently small, a desired boundary surface can be closely approximated.

If the origin of each subvolume is located at [{\bf r}_{i}], relative to a coordinate system origin at P, the total electronic moment relative to this origin is given by [ \mu ^l\left (\sum \limits _{i}V_{T,{i}}\right) = {1\over V}\sum_{{\bf h}}F({\bf h}) S_{V_T}(\widehat {\gamma }_l,{\bf h}) \sum _{i}\exp (-2\pi i{\bf h}\cdot {\bf r}_{i}) . \eqno (8.7.3.49)]

Expressions for [S_{V_T}] for [l\leq 2] and a subvolume parallelepipedal shape are given in Table 8.7.3.2[link]. Since the spherical order Bessel functions [j_n(x)] that appear in the expressions generally decrease with increasing x, the moments are strongly dependent on the low-order reflections in a data set. An example is the shape transform for the dipole moment. Relative to an origin O, [ S(\widehat {\gamma }_1,{\bf h}) =\textstyle\int\limits_{V_t}{\bf r}_0\exp (-2\pi i{\bf h}\cdot {\bf r}_0){\,{\rm d}}{\bf r}.]A shift of origin by [-{\bf r}_{i}] leads to [\eqalignno{ \mu ^{1^{\prime }}(V_T) &= {1\over V}\sum F({\bf h}) \left[S_{V_T}(\widehat {\gamma}_1,{\bf h}) +{\bf r}_1\int \exp (-2\pi i{{\bf h}\cdot r}_0) {\,{\rm d}}{\bf r}\right] \cr &\quad\times \exp (-2\pi i{\bf h}\cdot {\bf r}_1) \cr &=\mu ^1(V_T) +{\bf r}_1q, }]in agreement with (8.7.3.46)[link].

Table 8.7.3.2| top | pdf |
Expressions for the shape factors S for a parallelepiped with edges δx, δy, and δz (from Moss & Coppens, 1981[link])

j0 and j1 are the zero- and first-order spherical Bessel functions: j0(x) = sin x/x, j1(x) = sin x/x2 − cos x/x; VT is volume of integration.

[\hat y]Property[S[\hat y({\bf r}), {\bf h}]]
1 Charge [V_T\,j_0(2\pi h_x\delta_x)\, j_0(2\pi h_y\delta_y)\, j_0(2\pi h_z\delta_z)]
rα Dipole μα [\eqalign { -i&V_T\delta_\alpha\,j_1(2\pi h_\alpha\delta_\alpha) \cr & \times\, j_0(2\pi h_\beta\delta_\beta)\,j_0(2\pi h_y\delta_y)}]
rαrβ Second moment
μαβ off-diagonal
[\eqalign { -V_T&\delta_\alpha\delta_\beta\,j_1 (2\pi h_\alpha\delta_\alpha) \cr \times\,& j_1(2\pi h_\beta\delta_\beta)\,j_0(2\pi h_\gamma \delta_\gamma) }]
rαrα Second moment
μαα diagonal
[\eqalign { -V_T&\delta^2_\alpha\bigg\{\displaystyle{j_1(2\pi h_\alpha\delta_\alpha) \over \pi h_\alpha \delta_\alpha} - j_0(2\pi h_\alpha\delta_\alpha)\bigg\} \cr \quad \times\, &j_0(2\pi h_\beta\delta_\beta)\,j_0(2\pi h_\gamma\delta_\gamma) }]

8.7.3.4.2. The electrostatic potential

| top | pdf |

8.7.3.4.2.1. The electrostatic potential and its derivatives

| top | pdf |

The electrostatic potential Φ(r′) due to the electronic charge distribution is given by the Coulomb equation, [ \Phi({\bf r}^{\prime }) =-k\int \displaystyle{\rho({\bf r}) \over | {\bf r}-{\bf r}^{\prime }| }{\,{\rm d}}{\bf r},\eqno (8.7.3.50)]where the constant k is dependent on the units selected, and will here be taken equal to 1. For an assembly of positive point nuclei and a continuous distribution of negative electronic charge, we obtain [ \Phi({\bf r}^{\prime }) =\sum _{{M}}\displaystyle {Z_{{M}}\over | {\bf R}_{{M}}-{\bf r}^{\prime } | }-\int \displaystyle {\rho ({\bf r}) \over | {\bf r}-{\bf r}^{\prime } | }{\,{\rm d}}{\bf r},\eqno (8.7.3.51)]in which ZM is the charge of nucleus M located at RM.

The electric field E at a point in space is the gradient of the electrostatic potential at that point. [ {\bf E}({\bf r}) =-\nabla \Phi ({\bf r}) =-{\bf i}{\partial \Phi ({\bf r}) \over\partial x}-{\bf j} {\partial \Phi ({\bf r}) \over \partial y}-{\bf k} {\partial \Phi({\bf r}) \over \partial z}.\eqno (8.7.3.52)]As E is the negative gradient vector of the potential, the electric force is directed `downhill' and proportional to the slope of the potential function. The explicit expression for E is obtained by differentiation of the operator |rr′|−1 in (8.7.3.50)[link] towards x, y, z and subsequent addition of the vector components. For the negative slope of the potential in the x direction, one obtains [ {\bf E}_x({\bf r}^{\prime }) =\int \displaystyle {\rho _{{\rm total}}({\bf r}) \over | {\bf r}^{\prime }-{\bf r} | ^2}\, {({\bf r}^{\prime }-{\bf r}) _x \over | {\bf r}^{\prime }-{\bf r} | }{\,{\rm d}}{\bf r}=\int \displaystyle {\rho _{{\rm total}}({\bf r}) \over | {\bf r}^{\prime }-{\bf r}| ^3}({\bf r}^{\prime }-{\bf r}) _x{\,{\rm d}}{\bf r},\eqno (8.7.3.53)]which gives, after addition of the components, [ {\bf E}({\bf r'}) =-\nabla \Phi({\bf r}^{\prime }) =\int \displaystyle {\rho _{t}({\bf r}) ({\bf r}^{\prime }-{\bf r}) \over \left | {\bf r}-{\bf r}^{\prime }\right | ^3}{\rm d}{\bf r}. \eqno (8.7.3.54)]

The electric field gradient (EFG) is the tensor product of the gradient operator [\nabla={\bf i}{\partial\over {\partial x}}+{\bf j}{\partial \over{\partial y}}+{\bf k}{\partial\over {\partial z}}] and the electric field vector E; [ \nabla {\bf E}=\nabla: {\bf E}=-\nabla: \nabla \Phi. \eqno (8.7.3.55)]It follows that in a Cartesian system the EFG tensor is a symmetric tensor with elements [ \nabla {\bf E}_{\alpha \beta }=- {\partial ^2\Phi \over \partial r_\alpha \partial r_\beta }. \eqno (8.7.3.56)]

The EFG tensor elements can be obtained by differentiation of the operator in (8.7.3.53)[link] for Eα to each of the three directions β. In this way, the traceless result [\eqalignno{ \nabla {\bf E}_{\alpha \beta }({\bf r}^{\prime }) &= {\partial {\bf E}_\alpha \over\partial (r_\beta -r_\beta ^{\prime }) } \cr &=-\int \displaystyle{1\over | {\bf r}-{\bf r}^{\prime } | ^5}\Big \{ 3(r_\alpha -r_\alpha ^{\prime }) (r_\beta -r_\beta ^{\prime }) \cr &\quad {\,}- | {\bf r}-{\bf r}^{\prime }| ^2\delta _{\alpha \beta }\Big\} \rho _{{\rm total}}({\bf r}) {\,{\rm d}}{\bf r} & (8.7.3.57)}]is obtained. We note that according to (8.7.3.57)[link] the electric field gradient can equally well be interpreted as the tensor of the traceless quadrupole moments of the distribution −2ρtotal(r)/|rr′|5 [see equation (8.7.3.17)[link]].

Definition (8.7.3.56)[link] and result (8.7.3.57)[link] differ in that (8.7.3.56)[link] does not correspond to a zero-trace tensor. The situation is analogous to the two definitions of the second moments, discussed above, and is illustrated as follows. The trace of the tensor defined by (8.7.3.56)[link] is given by [ -\nabla ^2\Phi =-\nabla \cdot \nabla \Phi =-\left ({\partial ^2\Phi \over \partial x^2}+ {\partial ^2\Phi \over \partial y^2}+ {\partial ^2\Phi \over \partial z^2}\right). \eqno (8.7.3.58)]Poisson's equation relates the divergence of the gradient of the potential [\Phi ({\bf r})] to the electron density at that point: [ \nabla ^2\Phi({\bf r}) =-4\pi \left [-\rho _{{e}}({\bf r}) \right] =4\pi \rho _{{e}}({\bf r}). \eqno (8.7.3.59)]Thus, the EFG as defined by (8.7.3.56)[link] is not traceless, unless the electron density at r is zero.

The potential and its derivatives are sometimes referred to as inner moments of the charge distribution, since the operators in (8.7.3.50)[link], (8.7.3.52)[link] and (8.7.3.54)[link] contain the negative power of the position vector. In the same terminology, the electrostatic moments discussed in §8.7.3.4.1[link] are described as the outer moments.

It is of interest to evaluate the electric field gradient at the atomic nuclei, which for several types of nuclei can be measured accurately by nuclear quadrupole resonance and Mössbauer spectroscopy. The contribution of the atomic valence shell centred on the nucleus can be obtained by substitution of the multipolar expansion (8.7.3.7)[link] in (8.7.3.57)[link]. The quadrupolar (l = 2) terms in the expansion contribute to the integral. For the radial function [R_l=\{{\zeta ^{n(l)+3}}/{ [n(l) +2] !}\}r^{n(l) }\exp(-\zeta r)] with n(l) = 2, the following expressions are obtained: [\eqalign{ \nabla {\bf E}_{11} &=+(3/5)\left(\pi P_{22+}-\sqrt {3}P_{20}\right) Q_r, \cr \nabla {\bf E}_{22} &=-(3/5) \left (\pi P_{22+}+\sqrt {3}P_{20}\right) Q_r, \cr \nabla {\bf E}_{33} &=+(6/5) \left (\sqrt {3}P_{20}\right) Q_r, \cr \nabla {\bf E}_{12} &=+(3/5) \left (\pi P_{22-}\right) Q_r, \cr \nabla {\bf E}_{13} &=+(3/5) \left (\pi P_{21+}\right) Q_r, \cr \nabla {\bf E}_{23} &=+(3/5) \left (\pi P_{21-}\right) Q_r,} \eqno(8.7.3.60)]with [\eqalignno{ Q_r &=\langle r^3 \rangle _{3{\,{\rm d}}} \cr &=\textstyle\int\limits _0^\infty [{R}(r) /r] {\,{\rm d}}r \cr &=(\kappa ^{\prime }\zeta) ^3/ [n_2(n_2+1) (n_2+2)] \cr &=(\kappa ^{\prime }\zeta) ^3/120,} ]in the case that n2 = 4 (Stevens, DeLucia & Coppens, 1980[link]).

The contributions of neighbouring atoms can be sub­divided into point-charge, point-multipole, and penetration terms, as discussed by Epstein & Swanton (1982[link]) and Su & Coppens (1992[link], 1994b[link]), where appropriate expressions are given. Such contributions are in particular important when short interatomic distances are involved. For transition-metal atoms in coordination complexes, the contribution of neighbouring atoms is typically much smaller than the valence contribution.

8.7.3.4.2.2. Electrostatic potential outside a charge distribution

| top | pdf |

Hirshfelder, Curtis & Bird (1954[link]) and Buckingham (1959[link]) have given an expression for the potential at a point [{\bf r}_{i}] outside a charge distribution: [\eqalignno{ \Phi ({\bf r}_{i}) &= {q\over r_{i}}+ {\mu _\alpha r_\alpha \over r_{i}^3}+ \textstyle{1\over2}\left [3r_\alpha r_\beta -r^2\delta _{\alpha \beta }\right] \displaystyle {\Theta _{\alpha \beta }\over r_{{i}}^5} \cr &\quad +\left [5r_\alpha r_\beta r_\gamma -r^2\left (r_\alpha \delta _{\beta \gamma }+r_\beta \delta _{\gamma \alpha }+r_\gamma \delta _{\alpha \beta }\right) \right] {\Omega_{\alpha \beta \gamma } \over 5r_{i}^7}+\ldots, \cr& & (8.7.3.61)}]where summation over repeated indices is implied. The outer moments q, μα, Φαβ and Ωαβγ in (8.7.3.61)[link] must include the nuclear contributions, but, for a point outside the distribution, the spherical neutral-atom densities and the nuclear contributions cancel, so that the potential outside the charge distribution can be calculated from the deformation density.

The summation in (8.7.3.61)[link] is slowly converging if the charge distribution is represented by a single set of moments. When dealing with experimental charge densities, a multicentre expansion is available from the analysis, and (8.7.3.61)[link] can be replaced by a summation over the distributed moments centred at the nuclear positions, in which case [r_{i}] measures the distance from a centre of the expansion to the field point. The result is equivalent to more general expressions given by Su & Coppens (1992[link]), which, for very large values of [r_{i}], reduce to the sum over atomic terms, each expressed as (8.7.3.61)[link]. The interaction between two charge distributions, A and B, is given by [ E_{AB}=\textstyle\int \Phi _A({\bf r}) \rho _B({\bf r}) {\,{\rm d}}{\bf r},]where ρB includes the nuclear charge distribution.

8.7.3.4.2.3. Evaluation of the electrostatic functions in direct space

| top | pdf |

The electrostatic properties of a well defined group of atoms can be derived directly from the multipole population coefficients. This method allows the `lifting' of a molecule out of the crystal, and therefore the examination of the electrostatic quantities at the periphery of the molecule, the region of interest for intermolecular interactions. The difficulty related to the origin term, encountered in the reciprocal-space methods, is absent in the direct-space analysis.

In order to express the functions as a sum over atomic contributions, we rewrite (8.7.3.51)[link], (8.7.3.54)[link] and (8.7.3.57)[link] for the electrostatic properties at point P as a sum over atomic contributions. [\eqalignno{ \Phi ({\bf R}_p) &=\sum _{{{M}}\neq {P}}\displaystyle {Z_{{M}}\over | {\bf R}_{{MP}} | }- \sum_{{M}}\int \displaystyle {\rho _{{e,M}} ({\bf r}_{{M}}) \over | {\bf r}_p | }{\,{\rm d}}{\bf r}_{{M}}, &(8.7.3.62)\cr {\bf E} ({\bf R}_p) &=-\sum _{ {{M}}\neq {P}}\displaystyle {Z_{{MP}}{\bf R}_{{MP}}\over | {\bf R}_{{MP}} | ^3}+\sum_{{M}}\int \displaystyle {{\bf r}_p\rho _{{e,M}} ({\bf r}_{{M}}) \over | {\bf r}_p | ^3}{\,{\rm d}}{\bf r}_{{M}}, &(8.7.3.63)\cr \nabla {\bf E}_{\alpha \beta } ({\bf R}_p) &=-\sum_{{M}\neq {P}}\displaystyle {Z_{{M}} (3{R}_\alpha {R}_\beta -\delta _{\alpha \beta } | {R}_{{MP}} | ^2) \over | {R}_{{MP}} | ^5} \cr &\quad+\sum_{{M}}\int \displaystyle {\rho _{{e,M}} ({\bf r}_{{M}}) (3r_\alpha r_\beta -\delta _{\alpha \beta } | {\bf r}_p | ^2) \over | {\bf r}_p | ^5}{\,{\rm d}}{\bf r}_{{M}}, & (8.7.3.64)}%fd8.7.3.63fd8.7.3.64]in which the exclusion of M = P only applies when the point P coincides with a nucleus, and therefore only occurs for the central contributions. ZM and RM are the nuclear charge and the position vector of atom M, respectively, while rP and rM are, respectively, the vectors from P and from the nucleus M to a point r, such that rP = rRP, and rM = rRM = rP + RPRM = rPRMP. The subscript M in the second, electronic part of the expressions refers to density functions centred on atom M.

Expressions for the evaluation of (8.7.3.62)–(8.7.3.64)[link] from the charge-density parameters of the multipole expansion have been given by Su & Coppens (1992[link]). They employ the Fourier convolution theorem, used by Epstein & Swanton (1982[link]) to evaluate the electric field gradient at the atomic nuclei. A direct-space method based on the Laplace expansion of 1/|Rpr| was reported by Bentley (1981[link]).

8.7.3.4.3. Electrostatic functions of crystals by modified Fourier summation

| top | pdf |

Expression (8.7.3.49)[link] is an example of derivation of electrostatic properties by direct Fourier summations of the structure factors. The electrostatic potential and its derivatives may be obtained in an analogous manner.

In order to obtain the electrostatic properties of the total charge distribution, it is convenient to define the `total' structure factor [F_{{\rm total}}({\bf h})] including the nuclear contribution, [ F_{{\rm total}}({\bf h})=F_N({\bf h}) -F({\bf h}), ]where [F_N({\bf h})=\textstyle\sum _jZ_j\exp (2\pi i{\bf h}\cdot {\bf R}_j)], the summation being over all atoms j with nuclear charge [Z_j], located at [R_j ]. If Φ(h) is defined as the Fourier transform of the direct-space potential, we have [ \Phi({\bf r}) =\textstyle\int \Phi({\bf h}) \exp (-2\pi i{\bf h}\cdot {\bf r}) {\,{\rm d}}{\bf h} ]and [ \nabla ^2\Phi ({\bf r}) =-4\pi ^2\textstyle\int h^2\Phi ({\bf h}) \exp(-2\pi i{\bf h}\cdot {\bf r}) {\,{\rm d}}{\bf h}, ]which equals −4πρtotal according to the Poisson equation (8.7.3.59)[link]. One obtains with [\eqalignno{ \rho _{{\rm total}}({\bf r}) &=\textstyle\int F_{{\rm total}}({\bf h}) \exp(-2\pi i{\bf h}\cdot {\bf r}) {\,{\rm d}}{\bf h}, \cr \Phi ({\bf h}) &= +\,F_{{\rm total}}({\bf h}) /\pi h^2, & (8.7.3.65)}]and, by inverse Fourier transformation of (8.7.3.65)[link], [ \Phi({\bf r}) = {1\over\pi V}\sum F_{{\rm total}}({\bf h}) /h^2\exp (-2\pi i{\bf h}\cdot {\bf r})\eqno (8.7.3.66)](Bertaut, 1978[link]; Stewart, 1979[link]). Furthermore, the electric field due to the electrons is given by [\eqalignno{ {\bf E}({\bf r}) &=-\nabla _r\Phi({\bf r}) =-\nabla \textstyle\int \Phi ({\bf h}) \exp(-2\pi i{\bf h}\cdot {\bf r}) {\,{\rm d}}{\bf h} \cr &=2\pi i\textstyle\int {\bf h}\Phi ({\bf h}) \exp  (-2\pi i{\bf h}\cdot {\bf r}) {\,{\rm d}}{\bf h}. }]Thus, with (8.7.3.65)[link], [ {\bf E}({\bf r}) = {{2i}\over V} \sum \left [F_{{\rm total}}({\bf h}) /h^2\right] {\bf h}\exp(-2\pi i{\bf h}\cdot {\bf r}), \eqno(8.7.3.67a)]which implies [ {\bf E}({\bf h}) = {{2i}\over {h^2}}{\bf h}F_{{\rm total}}({\bf h}). ]Similarly, the h Fourier component of the electric field gradient tensor with trace [4\pi \rho ({\bf r})] is [ [\nabla: {\bf E}]({\bf h}) =+4\pi ^2{\bf h}:{\bf h}\Phi ({\bf h}) =+4\pi {\bf h}:{\bf h}F_{{\rm total}}({\bf h}) /h^2, \eqno(8.7.3.67b)]where i : k represents the tensor product of two vectors. This leads to the expression for the electric field gradient in direct space, [ [\nabla: {\bf E}] ({\bf r}) = {4\pi \over V}\sum {\bf h}:{\bf h}F_{{\rm total}}({\bf h}) /h^2\exp(-2\pi i{\bf h}\cdot {\bf r}).\eqno (8.7.3.68)](The elements of h : h are the products [h_ih_j.])

The components of E and the elements of the electric field gradient defined by (8.7.3.67a[link]) and (8.7.3.68)[link] are with respect to the reciprocal-lattice coordinate system. Proper transforms are required to get the values in other coordinate systems. Furthermore, to get the traceless [\nabla {\bf E}] tensor, the quantity −(4π/3)ρe(r) = −(4π/3V)[\sum]F(h)exp(−2π ih · r) must be subtracted from each of the diagonal elements [\nabla {\bf E}_{{ii}}.]

The Coulombic self-electronic energy of the crystal can be obtained from [\eqalignno{ E_{{\rm Coulombic,\ electronic}} &= {\textstyle{1\over 2}}\int\int \displaystyle {\rho _{{e}}({\bf r}) \rho _{{e}}({\bf r}^{\prime }) \over | {\bf r}-{\bf r}^{\prime } | }{\,{\rm d}}{\bf r}{\,{\rm d}}{\bf r}^{\prime }, \cr &= \textstyle{1\over2}\int \Phi _{{e}}({\bf r}) \rho ({\bf r}) {\,{\rm d}}{\bf r}.}]Since [\textstyle\int \Phi({\bf r})\rho({\bf r})\,{\rm d}{\bf r}= \textstyle\int\Phi({\bf h})F({\bf h})\,{\rm d}{\bf h}] (Parseval's rule), the summation can be performed in reciprocal space, [ E_{{\rm Coulombic,\ electronic}}= {1\over2\pi V}\sum F^2({\bf h}) /h^2, \eqno (8.7.3.69a) ]and, for the total Coulombic energy, [ E_{{\rm Coulombic,\ total}}= {1\over2\pi V}\sum F_{{\rm total}}^2({\bf h}) /h^2, \eqno (8.7.3.69b) ]where the integral has been replaced by a summation.

The summations are rapidly convergent, but suffer from having a singularity at h = 0 (Dahl & Avery, 1984[link]; Becker & Coppens, 1990[link]). The contribution from this term to the potential cannot be ignored if different structures are compared. The term at h = 0 gives a constant contribution to the potential, which, however, has no effect on the energy of a neutral system. For polar crystals, an additional term occurs in (8.7.3.69a[link],b[link]), which is a function of the dipole moment D of the unit cell (Becker, 1990[link]), [ E_{{\rm Coulombic,\ total}}= {1\over 2\pi V}\sum F_{{\rm total}}^2({\bf h}) /h^2+ \displaystyle{2\pi \over3V}{D}^2. \eqno (8.7.3.69c) ]To obtain the total energy of the static crystal, electron exchange and correlation as well as electron kinetic energy contributions must be added.

8.7.3.4.4. The total energy of a crystal as a function of the electron density

| top | pdf |

One can write the total energy of a system as [ E=e_c\left [\rho \right] +T+E_{xc},\eqno (8.7.3.70)]where T is the kinetic energy, [E_{xc}] represents the exchange and electron correlation contributions, and [E_c[\rho] ], the Coulomb energy, discussed in the previous section, is given by [ E_c= {\textstyle{1\over2}}\sum_{{i}\neq j} \displaystyle{Z_{i}Z_j \over R_{{i}j}}-\sum_{i}Z_{i}\Phi({\bf R}_{i}) + {\textstyle {1\over2}}\int\int \displaystyle {\rho ({\bf r}) \rho ({\bf r}^{\prime }) \over  | {\bf r}-{\bf r}^{\prime }| }{\,{\rm d}}{\bf r}{\,{\rm d}}{\bf r}^{\prime },\eqno (8.7.3.71)]where [Z_{i}] is the nuclear charge for an atom at position Ri, and Rij = RjRi.

Because of the theorem of Hohenberg & Kohn (1964[link]), E is a unique functional of the electron density ρ, so that [T+E_{xc}] must be a functional of ρ. Approximate density functionals are discussed extensively in the literature (Dahl & Avery, 1984[link]) and are at the centre of active research in the study of electronic structure of various materials. Given an approximate functional, one can estimate non-Coulombic contributions to the energy from the charge density ρ(r).

In the simplest example, the functionals are those applicable to an electronic gas with slow spatial variations (the `nearly free electron gas'). In this approximation, the kinetic energy T is given by [ T=c_k\textstyle\int \rho t[\rho] {\,{\rm d}}^3{\bf r},\eqno (8.7.3.72)]with [c_k=(3/{10})(2\pi {^2}){^{2/3}}]; and the function [t[\rho] =\rho ^{2/3} ]. The exchange-correlation energy is also a functional of ρ,[ E_{xc}=-c_x\textstyle\int \rho e_{xc}\left [\rho \right] {\,{\rm d}}^3{\bf r},]with [c_x=(3/4)(3/\pi)^{1/3}] and [\ e_{xc}[\rho] =\rho ^{1/3}].

Any attempt to minimize the energy with respect to ρ in this framework leads to very poor results. However, cohesive energies can be described quite well, assuming that the change in electron density due to cohesive forces is slowly varying in space.

An example is the system AB, with closed-shell subsystems A and B. Let ρA and ρB be the densities for individual A and B subsystems. The interaction energy is written as [\eqalignno{ \Delta E &=E_c [\rho] -E_c [\rho _A] -E_c [\rho _B] \cr &\quad +c_k\textstyle\int {\,{\rm d}}{\bf r}\left \{ \rho t [\rho] -\rho _At [\rho _A] -\rho _Bt [\rho _B] \right \} \cr &\quad-c_x\textstyle\int {\,{\rm d}}{\bf r}\left \{ \rho e_{xc} [\rho] -\rho _Ae_{xc} [\rho _A] -\rho _Be_{xc} [\rho _B] \right \}. & (8.7.3.73)}]This model is known as the Gordon–Kim (1972[link]) model and leads to a qualitatively valid description of potential energy surfaces between closed-shell subsystems. Unlike pure Coulombic models, this density functional model can lead to an equilibrium geometry. It has the advantage of depending only on the charge density ρ.

8.7.3.5. Quantitative comparison with theory

| top | pdf |

Frequently, the purpose of a charge density analysis is comparison with theory at various levels of sophistication. Though the charge density is a detailed function, the features of which can be compared at several points of interest in space, it is by no means the only level at which comparison can be made. The following sequence represents a progression of functions that are increasingly related to the experimental measurement. [\matrix{ &{\rm electrostatic} &\cr &{\rm properties} &\cr &\uparrow &\cr \psi(1,2,\ldots, n) \rightarrow \Gamma ^1(1,1)\hskip -3pt&\rightarrow \rho({\bf r}) \rightarrow &\hskip -5pt\left \langle \rho({\bf r}) \right \rangle \rightarrow F({\bf h}) \rightarrow I({\bf h}) \cr &\downarrow \cr &\Delta\rho({\bf r}) \rightarrow &\hskip -5pt\left \langle \Delta\rho ({\bf r}) \right \rangle, \hfill } \eqno (8.7.3.74)]where the angle brackets refer to the thermally averaged functions.

The experimental information may be reduced in the opposite sequence: [\matrix{ {\rm electrostatic} &\cr {\rm properties} &\cr \uparrow &\cr \rho ({\bf r})&\leftarrow & \left \langle \rho({\bf r}) \right \rangle &\leftarrow &F({\bf h}) &\leftarrow &I({\bf h}) \cr \downarrow &&\downarrow \cr \Delta\rho ({\bf r}) &\leftarrow & \left \langle \Delta\rho ({\bf r}) \right \rangle.} \eqno (8.7.3.75)]

The crucial step in each sequence is the thermal averaging (top sequence) or the deconvolution of thermal motion (bottom sequence), which in principle requires detailed knowledge of both the internal (molecular) and the external (lattice) modes of the crystal. Even within the generally accepted Born–Oppenheimer approximation, this is a formidable task, which can be simplified for molecular crystals by the assumption that the thermal smearing is dominated by the larger-amplitude external modes. The procedure (8.7.3.65)[link] requires an adequate thermal-motion model in the structure-factor formalism applied to the experimental structure amplitudes. The commonly used models may include anharmonicity, as described in Volume B, Chapter 1.2[link] (ITB, 2001[link]), but assume that a density function centred on an atom can be assigned the thermal motion of that atom, which may be a poor approximation for the more diffuse functions.

The missing link in scheme (8.7.3.75)[link] is the sequence [\psi \leftarrow \Gamma ^1(1,1)\leftarrow \rho ({\bf r})]. In order to describe the wavefunction analytically, a basis set is required. The number of coefficients in the wavefunction is minimized by the use of a minimal basis for the molecular orbitals, but calculations in general lead to poor-quality electron densities.

If the additional approximation is made that the wavefunction is a single Slater determinant, the idempotency condition can be used in the derivation of the wavefunction from the electron density. A simplified two-valence-electron two-orbital system has been treated in this manner (Massa, Goldberg, Frishberg, Boehme & La Placa, 1985[link]), and further developments may be expected.

A special case occurs if the overlap between the orbitals on an atom and its neighbours is very small. In this case, a direct relation can be derived between the populations of a minimal basis set of valence orbitals and the multipole coefficients, as described in the following sections.

8.7.3.6. Occupancies of transition-metal valence orbitals from multipole coefficients

| top | pdf |

In general, the atom-centred density model functions describe both the valence and the two-centre overlap density. In the case of transition metals, the latter is often small, so that to a good approximation the atomic density can be expressed in terms of an atomic orbital basis set [d_{{i}} ], as well as in terms of the multipolar expansion. Thus, [\eqalignno{ \rho _d &=\textstyle\sum \limits _{{i}=1}^5\textstyle\sum \limits _{j\ge {i}}^5P_{{i}j} d_{i} d_j, \cr &=\textstyle\sum \limits _{l=0}^4\kappa ^{\prime 3}\left \{ R_l(\kappa ^{\prime }r) \textstyle\sum \limits _{m=0}^l\textstyle\sum \limits _pP_{lmp}d_{lmp}\right \}, & (8.7.3.76)}]in which [d_{lmp}] are the density functions.

The orbital products [d_{i}d_j] can be expressed as linear combinations of spherical harmonic functions, with coefficients listed in Volume B, Chapter 1.2[link] (ITB, 2001[link]), which leads to relations between the [P_{{i}j}] and [P_{lmp}]. In matrix notation, [{\bf P}_{lmp}={\bi M}{\bf P}_{{i}j},\eqno (8.7.3.77)]where [{\bf P}_{lmp}] is a vector containing the coefficients of the 15 spherical harmonic functions with l = 0, 2, or 4 that are generated by the products of d orbitals. The matrix M is also a function of the ratio of orbital and density-function normalization coefficients, given in Volume B, Chapter 1.2[link] (ITB, 2001[link]).

The d-orbital occupancies can be derived from the experimental multipole populations by the inverse expression, [ {\bf P}_{{i}j}={\bi M}^{-1}{\bf P}_{lmp}\eqno (8.7.3.78)](Holladay, Leung & Coppens, 1983[link]).

The matrix M−1 is given in Table 8.7.3.3[link]. Point-group-specific expressions can be derived by omission of symmetry-forbidden terms. Matrices for point group 4/mmm (square planar) and for trigonal point groups (3, [\bar{3}], 32, 3m, [\bar3m]) are listed in Tables 8.7.3.4[link] and 8.7.3.5[link], respectively. Point groups with and without vertical mirror planes are distinguished by the occurrence of both [d_{lm+}] and [d_{lm-}] functions in the latter case, and only [d_{lm+}] in the former, with m being restricted to n, the order of the rotation axis. The [d_{{ln}-}] functions can be eliminated by rotation of the coordinate system around a vertical axis through an angle [\psi_0] given by [\psi_0=(1/n)\arctan (P_{{ln}-}/P_{{ln}+})].

Table 8.7.3.3| top | pdf |
The matrix M−1 relating d-orbital occupancies Pij to multipole populations Plm (from Holladay, Leung & Coppens, 1983[link])

d-orbital populationsMultipole populations
P00P20P22+P40P42+P44+
[P_z^2] 0.200 1.039 0.00 1.396 0.00 0.00
Pxz 0.200 0.520 0.942 −0.931 1.108 0.00
Pyz 0.200 0.520 −0.942 −0.931 −1.108 0.00
[P_{x^2-y^2}] 0.200 −0.039 0.00 0.233 0.00 1.571
Pxy 0.200 −1.039 0.00 0.233 0.00 −1.571

Mixing terms.

 P21P21−P22+P22−P41+P41−P42+P42−P43+P43−P44−
[P_{z^{2}/xz}] 1.088 0.00 0.00 0.00 2.751 0.00 0.00 0.00 0.00 0.00 0.00
[P_{z^{2}/yz}] 0.00 1.088 0.00 0.00 0.00 2.751 0.00 0.00 0.00 0.00 0.00
[P_{z^2/x^2-y^2}] 0.00 0.00 −2.177 0.00 0.00 0.00 1.919 0.00 0.00 0.00 0.00
[P_{z^{2}/xy}] 0.00 0.00 0.00 −2.177 0.00 0.00 0.00 1.919 0.00 0.00 0.00
[P_{xz/yz}] 0.00 0.00 0.00 1.885 0.00 0.00 0.00 2.216 0.00 0.00 0.00
[P_{xz/x^{2}-y^{2}}] 1.885 0.00 0.00 0.00 −0.794 0.00 0.00 0.00 2.094 0.00 0.00
[P_{xz/xy}] 0.00 1.885 0.00 0.00 0.00 −0.794 0.00 0.00 0.00 2.094 0.00
[P_{yz/x^{2}-y^{2}}] 0.00 −1.885 0.00 0.00 0.00 0.794 0.00 0.00 0.00 2.094 0.00
[P_{yz/xy}] 1.885 0.00 0.00 0.00 −0.794 0.00 0.00 0.00 −2.094 0.00 0.00
[P_{x^{2}-y^{2}/xy}] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.142

Table 8.7.3.4| top | pdf |
Orbital–multipole relations for square-planar complexes (point group D4h)

 P00P20P40P44+
P20 a1g 0.200 1.039 1.396 0.00
P21+ eg 0.200 0.520 −0.931 0.00
P21−   0.200 0.520 −0.931 0.00
P22+ b1g 0.200 −1.039 0.233 1.570
P22 b2g 0.200 −1.039 0.233 −1.570

Table 8.7.3.5| top | pdf |
Orbital–multipole relations for trigonal complexes

 P00P2−P40P43+P43−
(a) In terms of d orbitals
P20 0.200 1.039 1.396 0.00 0.00
P21+ 0.200 0.520 −0.931 0.00 0.00
P21− 0.200 0.520 −0.931 0.00 0.00
P22+ 0.200 −1.039 0.233 0.00 0.00
P22− 0.200 −1.039 0.233 0.00 0.00
P21+/22+ 0.00 0.00 0.00 2.094 0.00
P21+/22− 0.00 0.00 0.00 0.00 2.094
P21−/22+ 0.00 0.00 0.00 0.00 2.094
P21−/22− 0.00 0.00 0.00 −2.094 0.00
(b) In terms of symmetry-adapted orbitals
P1(a1g) 0.200 1.039 1.396 0.00  
P2(eg) 0.400 −1.039 −0.310 −1.975  
[P_3(e_g^\prime)] 0.400 0.00 −1.087 1.975  
[P_4(e_{g^+}e'\hskip-3pt_g+e_{g^-}e'_{g^-})] 0.00 −2.942 2.193 1.397  
The electron density in terms of the symmetry-adapted orbitals is given by: [\rho_{3d}=P_1a^2_{1g}+\textstyle{1\over2}P_2(e^2_{g^+}-e^2_{g^-})+{1\over 2}P_4(e_{g^+}e'\hskip-3pt_g+e_{g^-}e'\hskip-3pt_{g^-}),]with: [a_{1g}=d_{z^2}]; [e_{g^+}=\sqrt(2/3)d_{x^2-y^2} - \sqrt(1/3)d_{xz}]; [e_{g^-}= \sqrt(2/3)d_{xy}+\sqrt(1/3)d_{yz}]; [e'\hskip-3pt_{g^+} = \sqrt(1/3)d_{x^2-y^2}+\sqrt(2/3)d_{xz}]; and [e'_{g^-}=\sqrt(1/3)d_{xy} -\sqrt(2/3)d_{yz}].
The signs given here imply a positive [e'\hskip-3pt_g] lobe in the positive xz quadrant. Care should be exercised in defining the coordinate system if this lobe is to point towards a ligand atom.

8.7.3.7. Thermal smearing of theoretical densities

| top | pdf |

8.7.3.7.1. General considerations

| top | pdf |

In the Born–Oppenheimer approximation, the electrons rearrange instantaneously to the minimum-energy state for each nuclear configuration. This approximation is generally valid, except when very low lying excited electronic states exist. The thermally smeared electron density is then given by[\left \langle \rho ({\bf r}) \right \rangle =\textstyle\int \rho ({\bf r},{\bf R}) P({\bf R}) {\,\rm d} ({\bf R}), \eqno(8.7.3.79)]where R represents the 3N nuclear space coordinates and P(R) is the probability of the configuration R. Evaluation of (8.7.3.79)[link] is possible if the vibrational spectrum is known, but requires a large number of quantum-mechanical calculations at points along the vibrational path. A further approximation is the convolution approximation, which assumes that the charge density near each nucleus can be convoluted with the vibrational motion of that nucleus, [ \left \langle \rho({\bf r}) \right \rangle =\textstyle\sum \limits _n\textstyle\int \rho _n ({\bf r}-{\bf u}-{\bf R}_n) P_n ({\bf u}) {\,{\rm d}}{\bf u}, \eqno (8.7.3.80)]where ρn stands for the density of the nth pseudo-atom. The convolution approximation thus requires decomposition of the density into atomic fragments. It is related to the thermal-motion formalisms commonly used, and requires that two-centre terms in the theoretical electron density be either projected into the atom-centred density functions, or assigned the thermal motion of a point between the two centres. In the LCAO approximation (8.7.3.9)[link], the two-centre terms are represented by [ \rho _{\mu \nu }({\bf r}) =P_{\mu \nu }\chi _\mu ({\bf r}-{\bf R}_\mu) \chi _\nu ({\bf r}-{\bf R}_\nu ), \eqno (8.7.3.81)]where χμ and χν are basis functions centred at Rμ and Rν, respectively. As the motion of a point between the two vibrating atoms depends on their relative phase, further assumptions must be made. The simplest is to assume a gradual variation of the thermal motion along the bond, which gives at a point ri on the internuclear vector of length Rμν [ U_{{i}j}({\bf r}_{i}) =\left [U_{{i}j} ({\bf R}_\mu) | {\bf R}_\mu -{\bf r}_{i} | +U_{{i}j} ({\bf R}_\nu) | {\bf R}_\nu -{\bf r}_{i}| \right] /R_{\mu \nu }. \eqno (8.7.3.82)]This expression may be used to assign thermal parameters to a bond-centred function.

8.7.3.7.2. Reciprocal-space averaging over external vibrations

| top | pdf |

Thermal averaging of the electron density is considerably simplified for modes in which adjacent atoms move in phase. In molecular crystals, such modes correspond to rigid-body vibrations and librations of the molecule as a whole. Their frequencies are low because of the weakness of intermolecular interactions. Rigid-body motions therefore tend to dominate thermal motion, in particular at temperatures for which kT (k = 0.7 cm−1) is large compared with the spacing of the vibrational energy levels of the external modes (internal modes are typically not excited to any extent at or below room temperature).

For a translational displacement (u), the dynamic density is given by [ \rho _{\rm dyn} ({\bf r}) =\textstyle\int \rho ({\bf r}-{\bf u}) P({\bf u}) {\,{\rm d}}{\bf u}, \eqno (8.7.3.83)]with ρ(r) defined by (8.7.3.81)[link] (Stevens, Rees & Coppens, 1977[link]). In the harmonic approximation, P(u) is a normalized three-dimensional Gaussian probability function, the exponents of which may be obtained by rigid-body analysis of the experimental data. In general, for a translational displacement (u) and a librational oscillation (ω), [ \rho _{\rm dyn}({\bf r}) =\left (\textstyle\sum P_{\mu \nu }\chi _\mu \chi _\nu \right) *P({\bf u},\omega ). \eqno (8.7.3.84)]If correlation between u and ω can be ignored (neglect of the screw tensor S), P(u, ω) = P(u)P(ω), and both types of modes can be treated independently. For the translations [\eqalignno{ \left \langle \rho \right \rangle _{{\rm trans}} &=\textstyle\sum P_{\mu \nu } \{ \chi _\mu ({\bf r}) \chi _\nu ({\bf r}) *P({\bf u})  \} \cr &=\textstyle\sum P_{\mu \nu }F^{-1} \{\, f_{\mu \nu }({\bf h}) \cdot T_{{\rm tr}}({\bf h})\}, & (8.7.3.85)}]where F−1 is the inverse Fourier transform operator, and Ttr(h) is the translational temperature factor.

If R is an orthogonal rotation matrix corresponding to a rotation ω, we obtain for the librations [\eqalignno{ \left \langle \rho \right \rangle _{{\rm libr}} &=\textstyle\sum P_{\mu \nu } \{ \chi _\mu({\bi R}{\bf r}) \chi _v({\bi R} {\bf r}) *P (\omega) \} \cr &=\textstyle\sum P_{\mu \nu }F^{-1}\left \{\left\langle \,f_{\mu \nu }({\bi R}{\bf h}) \right \rangle \right \}, & (8.7.3.86)}]in which f(h) has been averaged over the distribution of orientations of h with respect to the molecule; [ \left \langle \,f_{\mu \nu }({\bf h}) \right \rangle =\textstyle\int f_{\mu \nu } ({\bi R}{\bf h}) P (\omega ) {\,{\rm d}}\omega. \eqno (8.7.3.87)]

Evaluation of (8.7.3.85)[link] and (8.7.3.86)[link] is most readily performed if the basis functions ψ have a Gaussian-type radial dependence, or are expressed as a linear combination of Gaussian radial functions.

For Gaussian products of s orbitals, the molecular scattering factor of the product ψμψν = Nμexp[−αμ(rrA)2] × Nνexp [−αν (rrB)2], where Nμ and Nν are the normalization factors of the orbitals μ and ν centred on atoms A and B, is given by [\eqalignno{ f_{{\rm stat}}^{s,s}\left ({\bf h}\right) &=N_\mu N_\nu \exp \left (- {\alpha _\mu \alpha _\nu \over\alpha _\mu +\alpha _\nu }\left | {\bf r}_A-{\bf r}_B\right | ^2\right) \left ({\pi \over\alpha _\mu +\alpha _\nu }\right) ^{3/2} \cr &\quad \times \exp \left ({-\pi ^2\left | h\right | ^2 \over\alpha _\mu +\alpha _v}\right) \exp \left (2\pi i{\bf h}\cdot {\bf r}_c\right), & (8.7.3.88)}]where the centre of density [{\bf r}_c] is defined by rc = (αμ rA + αν rB)/(αμ + αν).

For the translational modes, the temperature-factor exponent [-2\pi ^2\sum _{i}\sum _jU_{ij}h_{i}h_j] is simply added to the Gaussian exponent in (8.7.3.88)[link] to give [ \exp \left (-{\pi ^2 | h | ^2 \over \alpha _\mu +\alpha _\nu }\right) -2\pi ^2\sum \limits _{i}\sum \limits _jU_{{i}j}h_{{i}}h_j. ]For librations, we may write [ {\bi R}{\bf r}={\bf r}+{\bf u}_{{\rm lib}} ]As [({\bi R}{\bf h})\cdot {\bf r}={\bf h}\cdot{\bi R}^T{\bf r}={\bf h}\cdot {\bf r}-({\bi R}{\bf h})\cdot {\bf u}_{{\rm lib}}], for a function centred at r, [\eqalignno{ \left \langle \exp (2\pi {\bf h}\cdot {\bf r}_c) \right \rangle &=\textstyle\sum \exp [2\pi i({\bi R}{\bf h}) \cdot {\bf r}_c] P (\delta) {\,{\rm d}}\omega \cr &=\exp (2\pi {\bf h}\cdot {\bf r}_c) \textstyle\int \exp [-2\pi i({\bi R}{\bf h}) \cdot {\bf u}_{{\rm lib}}^c] P(\omega ) {\,{\rm d}}\omega, \cr&& (8.7.3.89)}]which shows that for ss orbital products the librational temperature factor can be factored out, or [ f_{\,{\rm dyn}}^{s,s}=f_{\rm stat}^{s,s}\textstyle\int \exp [-2\pi i({\bi R}{\bf h}) \cdot{\bf u}_{\rm lib}^c] P(\omega ) \,{\rm d}\omega. \eqno (8.7.3.90)]Expressions for P(ω) are described elsewhere (Pawley & Willis, 1970[link]).

For general Cartesian Gaussian basis functions of the type [ \psi({\bf r}) = (x-x_A) ^m (y-y_A) ^n (z-z_A) ^p\exp (-\alpha | {\bf r}-{\bf r}_A| ^2),\eqno (8.7.3.91)]the scattering factors are more complicated (Miller & Krauss, 1967[link]; Stevens, Rees & Coppens, 1977[link]), and the librational temperature factor can no longer be factored out. However, it may be shown that, to a first approximation, (8.7.3.90)[link] can again be used. This `pseudotranslation' approximation corresponds to a neglect of the change in orientation (but not of position) of the two-centre density function and is adequate for moderate vibrational amplitudes.

Thermally smeared density functions are obtained from the averaged reciprocal-space function by performing the inverse Fourier transform with phase factors depending on the position coordinates of each orbital product [\left \langle \rho \right \rangle = {1\over V}\sum P_{\mu \nu }\sum \limits _{{\bf h}}f_{\mu \nu }\left ({\bf h}\right) \exp \left [-2\pi i{\bf h}\cdot \left ({\bf r}-{\bf r}_c\right) \right] , \eqno (8.7.3.92)]where the orbital product χμχν is centred at [{\bf r}_c]. If the summation is truncated at the experimental limit of (sin θ)/λ, both thermal vibrations and truncation effects are properly introduced in the theoretical densities.

8.7.3.8. Uncertainties in experimental electron densities

| top | pdf |

It is often important to obtain an estimate of the uncertainty in the deformation densities in Table 8.7.3.1[link]. If it is assumed that the density of the static atoms or fragments that are subtracted out are precisely known, three sources of uncertainty affect the deformation densities: (1) the uncertainties in the experimental structure factors; (2) the uncertainties in the positional and thermal parameters of the density functions that influence ρcalc; and (3) the uncertainty in the scale factor k.

If we assume that the uncertainties in the observed structure factors are not correlated with the uncertainties in the refined parameters, the variance of the electron density is given by [\eqalignno{ \sigma ^2 (\Delta\rho) &=\sigma ^2 (\rho _{{\rm obs}}^{\prime }) +\sigma ^2(\rho _{{\rm calc}}) +(\rho _{{\rm obs}}^{\prime }) ^2 \, {\sigma ^2(k) \over k^2} \cr &\quad+\sum \limits _p\displaystyle {\delta \rho _{{\rm calc}} \over \delta u_p}\sigma(u_p) \rho _{{\rm obs}}^{\prime } {\sigma(k) \over k}\gamma (u_p,k), & (8.7.3.93)}]where [\rho_{\rm obs}^\prime =\rho_{\rm obs}/k], up is a positional or thermal parameter, and the γ(up, k) are correlation coefficients between the scale factor and the other parameters (Rees, 1976[link], 1978[link]; Stevens & Coppens, 1976[link]).

Similarly, for the covariance between the deformation densities at points A and B, [\eqalignno{ {\rm cov}(\Delta\rho _A,\Delta\rho _B) &={\rm cov} (\rho _{{\rm obs},A},\rho _{{\rm obs},B}) /k^2 \cr &\quad+{\rm cov}(\rho _{{\rm calc},A},\rho _{{\rm calc},B}) \cr &\quad+\rho _{{\rm obs},A}\rho _{{\rm obs},B} [\sigma (k) /k] ^2, & (8.7.3.94)}]where it is implied that the second term includes the effect of the scale factor/parameter correlation.

Following Rees (1976[link]), we may derive a simplified expression for the covariance valid for the space group [P\bar 1]. Since the structure factors are not correlated with each other, [\eqalignno{ {\rm cov} (\rho _{{\rm obs},A},\rho _{{\rm obs},B}) &\simeq \sum \displaystyle{\partial \rho _{{\rm obs},A} \over\partial F_{{\rm obs}}({\bf h}) }\, {\partial \rho _{{\rm obs},B} \over \partial F_{{\rm obs}}({\bf h}) }\sigma ^2\left [F_{{\rm obs}}({\bf h}) \right] \cr &\simeq { 2\over V^2}\sum _{1/2}\sigma ^2 (F_{{\rm obs}}) \left [\cos 2\pi({\bf r}_A+{\bf r}_B) \cdot {\bf h}\right. \cr &\left. \quad+\cos 2\pi ({\bf r}_A-{\bf r}_B) \cdot {\bf h}\right] , & (8.7.3.95)}]where the latter equality is specific for [P\bar 1], and [\sum _{1/2}] indicates summation over a hemisphere in reciprocal space. In general, the second term rapidly averages to zero as [h_{{\it \max }}] increases, while the first term may be replaced by its average [ \left \langle \cos 2\pi \left ({\bf r}_A-{\bf r}_B\right) \cdot {\bf h}\right \rangle =3\left (\sin u-u\cos u\right) /u^3\equiv {\rm C}\left (u\right), \eqno (8.7.3.96)]with u = 2π |rArB|hmax, or cov(ρobs, Aρobs, B) [\simeq] (2/V2)C(u) × [\sum _{1/2}\sigma ^2(F_0)] and σ2obs) [\simeq (2/{V^2})\sum _{1/2}\sigma^2(F_0)], a relation derived earlier by Cruickshank (1949[link]). A discussion of the applicability of this expression in other centrosymmetric space groups is given by Rees (1976[link]).

8.7.3.9. Uncertainties in derived functions

| top | pdf |

The electrostatic moments are functions of the scale factor, the positional parameters x, y and z of the atoms, and their charge-density parameters κ and [P_{lmp}]. The standard uncertainties in the derived moments are therefore dependent on the variances and covariances of these parameters.

If Mp represents the m × m variance–covariance matrix of the parameters [p_j], and T is an n × m matrix defined by [\partial \mu ^l/\partial p_j] for the lth moment with n independent elements, the variances and covariances of the elements of ml are obtained from [ {\bi M}_\mu ={\bi TM}_p{\bi T}^T. \eqno (8.7.3.96a) ]If a moment of an assembly of pseudo-atoms is evaluated, the elements of T include the effects of coordinate rotations required to transfer atomic moments into a common coordinate system.

A frequently occurring case of interest is the evaluation of the magnitude of a molecular dipole moment and its standard deviation. Defining [ y=\mu ^2={\boldmu}{\bi G}{\boldmu}^T, \eqno (8.7.3.96b)]where [{\boldmu}] is the dipole-moment vector and G is the direct-space metric tensor of the appropriate coordinate system. If Y is the 1 × 3 matrix of the derivatives [\partial y/\partial \mu _{{i}}], [ \sigma ^2(y) ={\bi YM}_\mu {\bi Y}^T, \eqno (8.7.3.96c)]where Mμ is defined by (8.7.3.96a[link]). The standard uncertainty in μ may be obtained from σ(μ) = σ(y)/2μ. Significant contributions often result from uncertainties in the positional and charge-density parameters of the H atoms.

References

Becker, P. (1990). Electrostatic properties from X-ray structure factors. Unpublished results.
Becker, P. & Coppens, P. (1985). About the simultaneous interpretation of charge and spin density data. Acta Cryst. A41, 177–182.
Becker, P. & Coppens, P. (1990). About the Coulombic potential in crystals. Acta Cryst. A46, 254–258.
Bentley, J. (1981). In Chemical Applications of Atomic and Molecular Electrostatic Potentials, edited by P. Politzer & D. G. Truhlar. New York/London: Plenum Press.
Bertaut, E. F. (1978). Electrostatic potentials, fields and gradients. J. Phys. Chem. Solids, 39, 97–102.
Boyd, R. J. (1977). The radial density function for the neutral atoms from helium to xenon. Can. J. Phys. 55, 452–455.
Buckingham, A. D. (1959). Molecular quadrupole moments. Q. Rev. Chem. Soc. 13, 183–214.
Buckingham, A. D. (1970). Physical chemistry. An advanced treatise, Vol. 4. Molecular properties, edited by D. Henderson, pp. 349–386. New York: Academic Press.
Coppens, P. (2001). The structure factor. International tables for crystallography, Vol. B, edited by U. Shmueli, Chap. 1.2. Dordrecht: Kluwer Academic Publishers.
Cromer, D. T., Larson, A. C. & Stewart, R. F. (1976). Crystal structure refinements with generalized scattering factors. J. Chem. Phys. 65, 336–349.
Cruickshank, D. W. (1949). The accuracy of electron density maps in X-ray analysis with special reference to dibenzyl. Acta Cryst. 2, 65–82.
Dahl, J. P. & Avery, J. (1984). Local density approximations in quantum chemistry and solid state physics. New York/London: Plenum.
Epstein, J. & Swanton, D. J. (1982). Electric field gradients in imidazole at 103 K from X-ray diffraction. J. Chem. Phys. 77, 1048–1060.
Gordon, R. G. & Kim, Y. S. (1972). Theory for the forces between closed-shell atoms and molecules. J. Chem. Phys. 56, 3122–3133.
Hamilton, W. C. (1964). Statistical methods in physical sciences. New York: Ronald Press.
Hellner, E. (1977). A simple refinement of density distributions of bonding electrons. Acta Cryst. B33, 3813–3816.
Hirshfeld, F. L. (1984). Hellmann–Feynman constraint on charge densities, an experimental test. Acta Cryst. B40, 613–615.
Hirshfeld, F. L. & Rzotkiewicz, S. (1974). Electrostatic binding in the first-row AH and A2 diatomic molecules. Mol. Phys. 27, 1319–1343.
Hirshfelder, J. O., Curtis, C. F. & Bird, R. B. (1954). Molecular theory of gases and liquids. New York: John Wiley.
Hohenberg, P. & Kohn, W. (1964). Inhomogeneous electron gas. Phys. Rev. B, 136, 864–867.
Holladay, A., Leung, P. C. & Coppens, P. (1983). Generalized relation between d-orbital occupancies of transition-metal atoms and electron-density multipole population parameters from X-ray diffraction data. Acta Cryst. A39, 377–387.
International Tables for Crystallography (2001). Vol. B, edited by U. Shmueli. Dordrecht: Kluwer Academic Publishers.
Levine, I. L. (1983). Quantum chemistry, 3rd edition. Boston/London/Sydney: Allyn and Bacon Inc.
Massa, L., Goldberg, M., Frishberg, C., Boehme, R. F. & La Placa, S. J. (1985). Wave functions derived by quantum modeling of the electron density from coherent X-ray diffraction: beryllium metal. Phys. Rev. Lett. 55, 622–625.
Miller, K. J. & Krauss, M. (1967). Born inelastic differential cross sections in H2. J. Chem. Phys. 47, 3754–3762.
Moss, G. & Coppens, P. (1981). Pseudomolecular electrostatic potentials from X-ray diffraction data. In Molecular electrostatic potentials in chemistry and biochemistry, edited by P. Politzer & D. Truhlar. New York: Plenum.
Pawley, G. S. & Willis, B. T. M. (1970). Temperature factor of an atom in a rigid vibrating molecule. II. Anisotropic thermal motion. Acta Cryst. A36, 260–262.
Rees, B. (1976). Variance and covariance in experimental electron density studies, and the use of chemical equivalence. Acta Cryst. A32, 483–488.
Rees, B. (1978). Errors in deformation-density and valence-density maps: the scale-factor contribution. Acta Cryst. A34, 254–256.
Stevens, E. D. & Coppens, P. (1976). A priori estimates of the errors in experimental electron densities. Acta Cryst. A32, 915–917.
Stevens, E. D., DeLucia, M. L. & Coppens, P. (1980). Experimental observation of the effect of crystal field splitting on the electron density distribution of iron pyrite. Inorg. Chem. 19, 813–820.
Stevens, E. D., Rees, B. & Coppens, P. (1977). Calculation of dynamic electron distributions from static molecular wave functions. Acta Cryst. A33, 333–338.
Stewart, R. F. (1977). One-electron density functions and many-centered finite multipole expansions. Isr. J. Chem. 16, 124–131.
Stewart, R. F. (1979) On the mapping of electrostatic properties from Bragg diffraction data. Chem. Phys. Lett. 65, 335–342.
Su, Z. & Coppens, P. (1992). On the mapping of electrostatic properties from the multipole description of the charge density. Acta Cryst. A48, 188–197.
Su, Z. & Coppens, P. (1994a). Rotation of real spherical harmonics. Acta Cryst. A50, 636–643.
Su, Z. & Coppens, P. (1994b). On the evaluation of integrals useful for calculating the electrostatic potential and its derivatives from pseudo-atoms. J. Appl. Cryst. 27, 89–91.








































to end of page
to top of page