On the calculation of the electrostatic potential, electric field and electric field gradient from the aspherical pseudoatom model

Volkov, A.; King, H.F.; Coppens, P.; Farrugia, L.J.

doi:10.1107/S0108767306026298

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 62| Part 5| September 2006| Pages 400-408

doi:10.1107/S0108767306026298

RESPONSE
A response to this article has been published. To view the response, click here

On the calculation of the electrostatic potential, electric field and electric field gradient from the aspherical pseudoatom model

Anatoliy Volkov,^a ^* Harry F. King,^a Philip Coppens ^a and Louis J. Farrugia ^b

^aDepartment of Chemistry, University at Buffalo, State University of New York, Buffalo, NY 14260-3000, USA, and ^bWestCHEM, Department of Chemistry, University of Glasgow, Glasgow G12 8QQ, UK
^*Correspondence e-mail: volkov@chem.buffalo.edu

(Received 26 May 2006; accepted 6 July 2006)

Accurate, yet simple and efficient, formulae are presented for calculation of the electrostatic potential (ESP), electric field (EF) and electric field gradient (EFG) from the aspherical Hansen–Coppens pseudoatom model of electron density [Hansen & Coppens (1978 ). Acta Cryst. A34, 909–921]. They are based on the expansion of |r′ − r|⁻¹ in spherical harmonics and the incomplete gamma function for a Slater-type function of the form R_l(r) = rⁿ exp(−αr). The formulae are valid for 0 ≤ r ≤ ∞ and are easily extended to higher values of l. Special treatment of integrals is needed only for functions with n = l and n = l + 1 at r = 0. The method is tested using theoretical pseudoatom parameters of the formamide molecule obtained via reciprocal-space fitting of PBE/6-31G** densities and experimental X-ray data of Fe(CO)₅. The ESP, EF and EFG values at the nuclear positions in formamide are in very good agreement with those directly evaluated from density-functional PBE calculations with 6-31G**, aug-cc-pVDZ and aug-cc-pVTZ basis sets. The small observed discrepancies are attributed to the different behavior of Gaussian- and Slater-type functions near the nuclei and to imperfections of the reciprocal-space fit. An EF map is displayed which allows useful visualization of the lattice EF effects in the crystal structure of formamide. Analysis of experimental 100 K X-ray data of Fe(CO)₅ yields the value of the nuclear quadrupole moment Q(⁵⁷Fe^m) = 0.12 × 10⁻²⁸ m² after taking into account Sternheimer shielding/antishielding effects of the core. This value is in excellent agreement with that reported by Su & Coppens [Acta Cryst. (1996 ), A52, 748–756] but slightly smaller than the generally accepted value of 0.16 ± 5% × 10⁻²⁸ m² obtained from combined theoretical/spectroscopic studies [Dufek, Blaha & Schwarz (1995). Phys. Rev. Lett. 25, 3545–3548].

Keywords: aspherical pseudoatom model; electric field; electric field gradient; electrostatic potential.

1. Introduction

The electrostatic potential V(r) (ESP) is one of the most important properties in the study of molecular reactivity and the analysis of molecular bonding and packing in crystals. It is related to the physically observable property of the electron density (ED) via Poisson's equation (Coppens, 1997 ):

$[\nabla^2 V({\bf r}) = - 4\pi\rho_{\rm total}({\bf r}), \eqno(1) ]$

where ρ_total includes both electronic and nuclear charges:

$[\rho _{\rm total} ({\bf r}) = \rho {}_{\rm nuc}({\bf r}) - \rho ({\bf r})\eqno(2) ]$

and ρ(r) is the electron density. For a continuous charge distribution, the potential is obtained by integration of the density over all space:

$[V({\bf r}) = \int {{{\rho _{\rm total} ({\bf r'})} \over {| {{\bf r'} - {\bf r}} |}}} \,{\rm d}^3 {\bf r'} ,\eqno(3) ]$

where r and r′ have an arbitrary common origin. The negative of the first derivative of the ESP is the electric field E(r) (EF):

$[{\bf E}({\bf r}) = - \nabla V({\bf r}) = - {\bf i}{{\partial V({\bf r})} \over {\partial x}} - {\bf j}{{\partial V({\bf r})} \over {\partial y}} - {\bf k}{{\partial V({\bf r})} \over {\partial z}}, \eqno(4) ]$

while the negative of the second derivative of the ESP is the electric field gradient (EFG):

$[E_{\alpha \beta } ({\bf r}) = - {{\partial ^2 V({\bf r})} \over {\partial r_\alpha \partial r_\beta }}.\eqno(5) ]$

Just as the total density ρ_total(r), the ESP in the molecule can be separated into electronic V^elec(r) and nuclear V^nuc(r) contributions:

$[\eqalignno{V({\bf r}) &= V^{\rm nuc} ({\bf r}) + V^{\rm elec} ({\bf r}) = \left({\sum\limits_i {V_i^{\rm nuc} } ({\bf r})} \right) + V^{\rm elec} ({\bf r}) \cr&= \sum\limits_{i = 1}^N {{{Z_i } \over {| {{\bf R}_i - {\bf r}} |}}} - \int {{{\rho ({\bf r'})} \over {| {{\bf r'} - {\bf r}} |}}} \,{\rm d}^3 {\bf r'},&(6)} ]$

where i = 1, …, N is the index of the nucleus located at R_i carrying positive charge Z_i. The second term represents the contribution of the continuous distribution of the negatively charged electron density ρ(r). When calculating V(r) at any point in space where $[{\bf r} \ne {\bf R}_i]$ , both terms need to be taken into account, while, at $[{\bf r} = {\bf R}_i]$ , the contribution of the nuclear potential of the ith nucleus (located exactly at R_i) must be omitted. The calculation of the nuclear potential and its derivatives is trivial and need not be described here. It is the calculation of the electronic potential V^elec(r) that presents more problems.

Various methods for calculating the ESP from X-ray diffraction data have been described and consequently applied in the literature. These methods can basically be split into two very different groups: (i) directly from experimentally measured structure factors (Bertaut, 1978 ; Stewart, 1979 ; Schwarzenbach & Thong, 1979 ) and (ii) from static models of electron density. Discussion of methods belonging to group (i) lies outside the scope of this paper. More important is the second group, which allows the calculation of electrostatic properties of atoms and molecules from the electron densities deconvoluted from nuclear motions. These are relatively unaffected by Fourier-series termination and allow a direct comparison with theoretical results. These methods are based on so-called pseudoatom charge-density models.

By far the most widely used aspherical pseudoatom formalism is based on the Hansen–Coppens multipolar model of ED (Coppens, 1997; Hansen & Coppens, 1978). In this formalism, the electron density at each point in space ρ(r) is described by a superposition of atomic like densities ρ_at(r), called pseudoatoms:

$[\rho (r) = \textstyle\sum {\rho _{\rm at} (r)}. \eqno(7) ]$

Each pseudoatom is modeled using the modified Laplace series:

$[\eqalignno{\rho _{\rm at} (r) &= P_{\rm core} \rho _{\rm core} (r) + P_{\rm val} \kappa ^3 \rho _{\rm val} (\kappa r) \cr &\quad+ \textstyle\sum\limits_{l = 0}^{l_{\rm max}} {\kappa '^3 R_l (\kappa'r) }\sum\limits_{m = 0}^l P_{lm \pm } d_{lm \pm } (\Omega). &(8)} ]$

The first and second terms of expansion are the spherically averaged Hartree–Fock core and valence densities (Clementi & Roetti, 1974 ). The population of the core, P_core, is always frozen, while the population of the spherical valence shell P_val is refined together with the expansion–contraction parameter κ. The third term describes the aspherical deformation density. Coefficients P_lm± are the population parameters and κ′ are the dimensionless adjustment coefficients of the radial functions R_l. In equation (8), r is the distance from the pseudoatom center, r = |r − R|, and Ω are the corresponding angular coordinates. The angular factor d_lm± is a real density-normalized spherical harmonic:

$[d_{lm\pm}(\Omega) = D_{lm} y_{lm\pm}(\Omega), \eqno(9) ]$

where m ≥ 0 and y_lm± is the normalized linear combination of complex spherical harmonics:

$[\eqalign{y_{00}& = Y_{00}\cr y_{lm + } & = ({ - 1}){}^m ({Y_{lm} + Y_{l, - m} })/\sqrt 2, \,\,\,\,\quad {m \,\gt\, 0} \cr y_{lm - } & = ({ - 1}){}^m ({Y_{lm} - Y_{l, - m} })/\sqrt { - 2}, \quad {m \,\gt\, 0} .}\eqno(10) ]$

Renormalization factors D_lm are given by Paturle & Coppens (1988 ) and Coppens (1997).

Calculation of the electrostatic potential and its first and second derivatives, i.e. negatives of electric field and electric field gradient, respectively, from the pseudoatom model is not straightforward. It was shown that the electrostatic potential can be evaluated from the pseudoatom model in various ways, e.g. directly from ρ(r) or from its truncated Fourier-series expansion (Brown & Spackman, 1994 ).

The method of Su & Coppens (1992 , 1996) is based on the Fourier convolution theorem previously applied by Epstein & Swanton (1982 ) to a calculation of the EFG. While it does not contain approximations and is formally exact at any point in space, it involves the evaluation of fairly complicated one-electron two-center integrals. This method was encoded in the computer program MOLPROP and later included in the experimental charge-density package XD (Koritsanszky et al., 2003 $[Koritsanszky, T., Howard, S. T., Richter, T., Macchi, P., Volkov, A., Gatti, C., Mallinson, P. R., Farrugia, L. J., Su, Z. & Hansen, N. K (2003). XD - a Computer Program Package for Multipole Refinement and Topological Analysis of Charge Densities from Diffraction Data. Middle Tennessee State University, TN, USA; University of Milano, Italy; University at Buffalo, NY, USA; CNR-ISTM, Milano, Italy; University of Glasgow, UK.]$ ). However, the program does not always reproduce the correct results when theoretical structure factors are used to evaluate the physical properties. It is not clear whether this is due to programming errors or errors in the derivation of two-center integrals [labeled as A_{N,l₁, l₂, k} (Z,R) in the original paper]. Note that the method of Su & Coppens only allows calculation of the traceless EFG tensor at the nuclear positions.

Several methods for the calculation of the ESP/EFG were proposed by Brown & Spackman (1994). The direct-space evaluation of the EFG, also based on the Fourier transform theorem, is closely related to that of Su & Coppens. It was included in the experimental charge-density package VALRAY (Stewart et al., 1998 ). These authors also presented the combined Fourier/direct-space evaluation of the EFG, as well as calculation of the EFG via numerical differentiation of the ESP. While giving more or less reasonable results, these methods are either too computationally demanding or have convergence problems with Fourier sums.

Ghermani, Bouhmaida, Lecomte and co-workers (Lecomte et al., 1992 ; Ghermani, Bouhmaida & Lecomte, 1993 ; Ghermani, Lecomte & Bouhmaida, 1993 ; Ghermani et al., 1994 ; Bouhmaida et al., 1997 ) reported expressions for the ESP derived directly from the Hansen–Coppens density model [equation (8)], in which the ESP due to a pseudoatom is expanded in the same way as the density ρ_at(r) itself [equation (8)], i.e.

$[V_{\rm at} ({\bf r}) = V_{\rm core} (r) + V_{\rm val} (r) + \Delta V({\bf r}).\eqno(11) ]$

Here V_core(r), V_val(r) and ΔV(r) are the spherical core, spherical valence and aspherical deformation contributions, respectively, and are given by

$[\eqalignno{V_{\rm core} ({\bf r}) &= {Z \over {| {{\bf r} - {\bf R}} |}} - \int\limits_0^\infty {{{\rho _{\rm core} (r')} \over {| {{\bf r} - {\bf R} - {\bf r'}} |}}\,{\rm d}^3 {\bf r'}}, &(12)\cr V_{\rm val} (r) &= - \int\limits_0^\infty {{{P_{\rm val} \kappa ^{\prime 3} \rho _{\rm val} (\kappa 'r')} \over {| {{\bf r} - {\bf R} - {\bf r'}} |}}} \,{\rm d}^3 {\bf r'}, &(13)\cr\Delta V({\bf r}) &= - 4\pi \sum\limits_{lm} {{\kappa 'P_{lm} } \over {2l + 1}} \bigg [{1 \over {\kappa '^{l + 1} | {{\bf r} - {\bf R}} |{}^{l + 1} }}\int\limits_0^{\kappa ^\prime| {{\bf r} - {\bf R}} |} {t^{l + 2} R_l (t )\,{\rm d}t}\cr &\quad + \kappa ^{\prime l} | {{\bf r} - {\bf R}} |{}^l \int\limits_{\kappa ^\prime| {{\bf r} - {\bf R}} |}^\infty {{R_l (t )} \over {t^{l - 1} }}\,{\rm d}t \bigg]d_{lm} (\Omega),&(14)} ]$

where R is the position of the nucleus and r′ is the position vector relative to R [unlike that in equation (3)]. Note that, in expression (12), the contribution of the nuclear potential is explicitly combined with the core-electron contribution. Formula (14) was derived using Green's function and the property of orthogonality of the spherical harmonic functions and is valid for any point in space, since no approximations are used. These formulae were implemented in the program ELECTROS (Ghermani et al., 1992 $[Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1992). ELECTROS: Computer Program to Calculate Electrostatic Properties from High Resolution X-ray Diffraction. Internal Report URA CNRS 809. Université de Nancy I, France.]$ ) derived from the earlier program MOLPOT (He, 1983 ). Numerous studies have been published using this approach.

In the current paper, we review the derivation of equation (14), providing more detail than available in the literature. In particular, we call attention to practical problems associated with numerical instabilities, and present simple stable expressions for the ESP, EF and EFG that overcome these problems. We apply our methods in a theoretical study of the EF in the crystal structure of formamide and in a determination of the ⁵⁷Fe^m nuclear quadrupolar moment from experimental X-ray diffraction data.

2. Calculation of the electronic potential and its derivatives for Slater-type functions

Each term in the pseudoatom formalism (8) can be reduced to a linear combination of Slater-type (Slater, 1932 ) density basis functions with the general form

$[\zeta (r) = R(r)d_{lm \pm } (\Omega) = N(n,\alpha)r^n \exp({ - \alpha r}) d_{lm \pm } (\Omega),\eqno(15) ]$

where α is the effective exponent and N(n,α) is the normalization factor (Coppens, 1997). The corresponding electrostatic potential is

$[V_\zeta ({\bf r}) = - \int {{{\zeta ({\bf r'})} \over {| {{\bf r'} - {\bf r}}|}}\,{\rm d}^3 {\bf r'}}. \eqno(16) ]$

$[| {{\bf r'} - {\bf r}} |{}^{ - 1}]$ can be expanded in real spherical harmonics (Jackson, 1975 ):

$[{1 \over {| {{\bf r'} - {\bf r}}|}} = \sum\limits_l {{4\pi } \over {2l + 1}}{{r_ \lt ^l } \over {r_ \gt ^{l + 1} }}\sum\limits_m y_{lm \pm } (\Omega )y_{lm \pm } ({\Omega '} ) ,\eqno(17) ]$

where r_< is the smaller and r_> is the greater of $[| {{\bf r} - {\bf R}} |]$ and $[| {{\bf r'} - {\bf R}}|]$ . Note that the y_lm± are obtained by unitary transformation of the Y_lm, which implies that the form of the spherical harmonic addition theorem (Edmonds, 1974 ) is preserved, i.e.

$[\textstyle\sum\limits_{m = - l}^l {Y_{lm}^* (\Omega )Y_{lm} ({\Omega '})} = \sum\limits_{m = 0}^l {y_{lm \pm } (\Omega )y_{lm \pm } ({\Omega '} )}.\eqno(18) ]$

When (15) and (17) are inserted into (16), the only term in the double sum that survives after integration over Ω′ is that term for which the lm± indices match those in ζ(r). This follows immediately from orthonormality of the y_lm± functions. Thus,

$[V_\zeta ({\bf r}) = {{4\pi } \over {2l + 1}}N({n,\alpha })d_{lm \pm } (\Omega )\int\limits_0^\infty {{r_ \lt ^l } \over {r_ \gt ^{l + 1} }}({r'}){}^{n + 2} \exp({ - \alpha r' }) \,{\rm d}r'. \eqno(19) ]$

The radial integral in (15) splits into two terms:

$[\eqalignno{ I({n,l,\alpha r}) & = \int\limits_0^\infty {{{r_ \lt ^l } \over {r_ \gt ^{l + 1} }}r^{n + 2} \exp({ - \alpha r}) } \,{\rm d}r \cr & = {1 \over {| {{\bf r} - {\bf R}} |{}^{l + 1} }}\int\limits_0^{| {{\bf r} - {\bf R}} |} {r^{n + l + 2} \exp({ - \alpha r}) \,{\rm d}r} \cr&\quad+ | {{\bf r} - {\bf R}}|{}^l \int\limits_{| {{\bf r} - {\bf R}}|}^\infty {r^{n + 1 - l} \exp({ - \alpha r} )\,{\rm d}r} \cr& = I_1 ({n,l,\alpha r} ) + I_2 ({n,l,\alpha r} ), &(20)} ]$

so expression (19) becomes simply

$[V_\zeta ({\bf r}) = {{4\pi } \over {2l + 1}}N(n,\alpha)d_{lm \pm } (\Omega)I({n,l,\alpha r}).\eqno(21) ]$

Essentially the same formula is given by Ghermani, Bouhmaida, Lecomte and co-workers (Lecomte et al., 1992; Ghermani et al., 1992 $[Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1992). ELECTROS: Computer Program to Calculate Electrostatic Properties from High Resolution X-ray Diffraction. Internal Report URA CNRS 809. Université de Nancy I, France.]$ ; Ghermani, Lecomte & Bouhmaida, 1993; Ghermani, Bouhmaida & Lecomte, 1993; Ghermani et al., 1994; Bouhmaida et al., 1997). Integrals I₁(n,l,αr) and I₂(n,l,αr) have the form of the incomplete gamma function (te Velde, 1990 ; te Velde et al., 2001 ):

$[\eqalignno{I_1 ({n,l,\alpha r}) &= {{({n + l + 2} )!} \over {\alpha ^{n + 2} }}{1 \over {({\alpha r} ){}^{l + 1} }}\left[{1 - \exp({ - \alpha r}) \sum\limits_{i = 0}^{n + l + 2} {{{({\alpha r}){}^i } \over {i!}}} } \right],\cr&&(22)\cr I_2 ({n,l,\alpha r}) &= {{({n - l + 1} )!} \over {\alpha ^{n + 2} }}({\alpha r} ){}^l \exp({ - \alpha r}) \sum\limits_{i = 0}^{n - l + 1} {{{({\alpha r}){}^i } \over {i!}}}.&(23)} ]$

The exponential terms in equations (22) and (23) represent the effects of the interpenetration of the charge-density distributions. The remaining term in equation (22) is the radial factor for the potential of a point multipole. These expressions are satisfactory for large values of αr but their straightforward evaluation for small values of αr results in severe numerical round-off errors as illustrated in Table 1. Numerical errors are relatively insensitive to n but grow rapidly with increasing l and decreasing αr. The corresponding errors for first and second derivatives required for evaluation of EF and EFG are significantly greater than those for the ESP. The numerical errors are entirely due to the term I₁(n,l,αr), and can be avoided by rearranging equation (22):

$[\eqalignno{I_1 ({n,l,\alpha r}) &= {{({n + l + 2} )!} \over {\alpha ^{n + 2} }}{1 \over {(\alpha r){}^{l + 1} }}\exp({ - \alpha r}) \sum\limits_{i = n + l + 3}^\infty {{{({\alpha r}){}^i } \over {i!}}} \cr&= {{({n + l + 2} )!} \over {\alpha ^{n + 2} }}\exp({ - \alpha r}) \sum\limits_{i = n + 2}^\infty {{{({\alpha r} ){}^i } \over {(i + l + 1)!}}}.&(24)} ]$

This is a well behaved expression for any r. The rate of convergence of this infinite series is reported in Table 2. Thus, expression (24) is recommended for use for `small' values of αr while expression (22) is preferred for medium and large values of αr as it converges more rapidly in that case.

Table 1
Number of significant digits lost in direct evaluation of formulas (22) and (23)

	l
αr	0	1	2	3	4
0.25	1	4	6	9	12
0.50	1	2	5	8	9
0.75	0	3	4	6	9
1.00	0	2	4	5	7
1.25	0	1	3	5	7
1.50	1	2	1	4	6

Table 2
Number of terms in equations (24) and (25) needed to achieve given accuracy

	Number of significant digits
αr	10	15
0.25	10	13
0.50	12	16
0.75	14	18
1.00	15	20
1.25	17	21
1.50	18	23

Alternatively, integral I(n,l,αr) can be replaced with its Taylor-series expansion:

$[\eqalignno{I({n,l,\alpha r} ) &= {1 \over {\alpha ^{n + 2} }}\bigg [({n - l + 1} )!({\alpha r} ){}^l \cr&\quad- ({\alpha r}){}^{n + 2} \sum\limits_{k = 0}^\infty {{{({ - \alpha r} ){}^k ({2l + 1} )} \over {k!({k + n - l + 2} )({k + n + l + 3})}}} \bigg].\cr&&(25)} ]$

Table 2 reports the number of terms required to achieve given accuracy when computing both I(n,l,αr) and its derivatives. Errors are relatively insensitive to l and n.

Formula (21) can be directly applied to the deformation part of the pseudoatom expansion because the spherical harmonics of deformation functions are already density normalized (Coppens, 1997). The spherical core and valence densities are, however, calculated from products of wavefunction-normalized Slater functions. Nevertheless, the product of two Slater functions on the same center is still a Slater function:

$[\eqalignno{&N_1 r^{n_1 } \exp({ - \alpha _1 r}) N_2 r^{n_2 } \exp({ - \alpha _2 r}) \cr&\quad= (N_1 N_2)(r^{n_1 + n_2 })\exp[{ - (\alpha _1 + \alpha _1)r}] \cr&\quad= N_3 r^{n_3 } \exp({ - \alpha _3 r}),&(26)} ]$

therefore, formula (21) can still be applied.

The first (EF) and second (EFG) derivatives of the ESP are then obtained by straightforward differentiation of expression (21).

3. Calculation of the electronic potential and its derivatives near r = 0

Functions d_lm±(Ω) and their derivatives are poorly behaved near the origin, as they contain direction cosines, so we remove a factor r^l from I(n,l,αr) and incorporate it into the angular factor:

$[V_\zeta (r) = {{4\pi } \over {2l + 1}}N(n,\alpha) [{d_{lm \pm } (\Omega)r^l } ] [{r^{ - l} I({n,l,\alpha r} )} ].\eqno(27) ]$

Note that all terms in (25) contain powers of αr of l or higher. For l ≤ 4, derivatives of d_lm±(Ω)r^l functions are simple and well behaved at any r. For example, the second partial derivative of d₄₁₋(Ω)r⁴ w.r.t. xy is L₄₁₋(−6xz), where L₄₁₋ = 0.474 is the density function normalization factor (Paturle & Coppens, 1988). No second derivative is more complicated than that of d₄₀(Ω)r⁴ w.r.t. xx, which is simply L₄₀(36x² + 12y² − 48z²). Note that first derivatives of d_lm±(Ω)r^l vanish at r = 0 for all functions except for

$[{{\partial [{d_{11 + } r}]} \over {\partial x}} = {{\partial [{d_{11 - } r}]} \over {\partial y}} = {{ \partial[{d_{10} r}]} \over {\partial z}} = L_{1m} = 0.31831.\eqno(28) ]$

The second derivatives of d_lm±(Ω)r^l vanish at r = 0 for all functions except for

$[\eqalign{{{\partial ^2[{d_{21} r^2 }]} \over {\partial x\partial z}} & = {{\partial ^2 [{d_{21 - } r^2 } ]} \over {\partial y\partial z}} = {{\partial ^2 [{d_{22} r^2 }]} \over {\partial x\partial x}} = {{\partial ^2 [{d_{22 - } r^2 } ]} \over {\partial x\partial y}} = L_{2m},\cr {{\partial ^2 [{d_{20} r^2 }]} \over {\partial x^2 }} & = {{\partial ^2 [{d_{20} r^2 }]} \over {\partial y^2 }} = - 2 L_{20},\cr {{\partial ^2 [{d_{20} r^2 }]} \over {\partial z^2 }} &= 4 L_{20}, \quad {{\partial ^2 [{d_{22} r^2 } ]} \over {\partial y^2 }} = - L_{22}.}\eqno(29) ]$

We define function G(r) as follows:

$[G(r) = r^{ - l} I({n,l,\alpha r} ).\eqno(30) ]$

To simplify the notation of G, the dependence on the n and l indices is implicitly understood. Let

$[j = n -l \ge 0. \eqno(31) ]$

The G(r) function is finite and smooth at the origin and decays to zero as r^−2l−1 for large r. For small r, the Taylor expansion can be used. Table 2 is equally applicable to the Taylor-series expansion of G(r) and its derivatives:

$[G(r) = {1 \over {\alpha ^{j + 2} }} \left[{({j + 1})! + ({ - 1} ){}^{j + 1} \sum\limits_{k = j + 2}^\infty {{{({ - \alpha r} ){}^k ({2l + 1})} \over {k({k - j - 2} )!({2l + 1 + k} )}}} } \right].\eqno(32) ]$

Alternatively,

$[G(r) = {{({j + 1})!} \over {\alpha ^{j + 2} }} - r^{j + 2} \sum\limits_{k = 0}^\infty {{{({ - \alpha r} ){}^k ({2l + 1} )} \over {k!({k + j + 2} )({2l + 3 + j + k} )}}}. \eqno(33) ]$

For example, when n = l then

$[\eqalignno{G(r) &= {1 \over {\alpha ^2 }} - r^2 \left [{{{({2l + 1} )} \over {2(2l + 3)}} - {{({\alpha r} )({2l + 1} )} \over {3({2l + 4} )}} + {{({\alpha r}){}^2 ({2l + 1} )} \over {8({2l + 5} )}} - \ldots } \right], \cr&& j = 0.\qquad (34)} ]$

When n = l + 1, then

$[\eqalignno{G(r) &= {2 \over {\alpha ^3 }} - r^3 \left [{{({2l + 1})} \over {3(2l + 4)}} - {{({\alpha r} )({2l + 1} )} \over {4({2l + 5} )}} + {{({\alpha r} ){}^2 ({2l + 1} )} \over {10({2l + 6} )}} - \ldots \right],\cr&& j = 1.\qquad(35)} ]$

Expressions (34) and (35) both have non-zero terms in r³, which creates singularities in the second derivatives d²r³/dq², where q = x,y,z. This can be circumvented as discussed below.

Define A(r) as follows:

$[A(r) = \left({{{{\rm d}G(r)} \over {r\,{\rm d}r}}} \right) = {{({ - 1}){}^{j + 1} } \over {\alpha ^j }}\sum\limits_{k = j}^\infty {{{({ - \alpha r} ){}^k ({2l + 1})} \over {({k - j})!({2l + 3 + k})}}}. \eqno(36) ]$

Alternatively,

$[A(r) = - r^j \sum\limits_{k = 0}^\infty {{{({ - \alpha r}){}^k ({2l + 1})} \over {k!({2l + 3 + j + k} )}}}. \eqno(37) ]$

Note that A(r) is well behaved for small r for all l and n. Note also that (37) has a non-zero term linear in r when j = 0 or j = 1. These linear terms have implications for B(r).

Define B(r) as follows:

$[\eqalignno{B(r) &= \left[{{{\rm d} \over {r\,{\rm d}r}}\left({{{{\rm d}G(r)} \over {r{\rm d}r}}} \right)} \right] = \left({{{{\rm d}A(r)} \over {r\,{\rm d}r}}} \right)\cr& = {{({ - 1}){}^{j + 1} } \over {\alpha ^{j - 2} }}\sum\limits_{k = j - 2}^\infty {{{({ - \alpha r} ){}^k ({k + 2} )({2l + 1})} \over {({k + 2 - j})!({2l + 5 + k} )}}}. &(38)} ]$

Alternatively,

$[B(r) = - r^{j - 2} \sum\limits_{k = 0}^\infty {{{({ - \alpha r} ){}^k ({k + j} )({2l + 1})} \over {k!({2l + 3 + j + k} )}}}. \eqno(39) ]$

Note that the first term vanishes in the sums (38) and (39) when j = 0, giving

$[B(r) = {\alpha \over r}\sum\limits_{k = 0}^\infty {{{({ - \alpha r}){}^k ({2l + 1} )} \over {k!({2l + 4 + k} )}}}, \quad j = 0 .\eqno(40) ]$

Alternatively,

$[B(r) = {{\alpha ({2l + 1} )} \over {r({2l + 4} )}} - \alpha ^2 \sum\limits_{k = 0}^\infty {{{({ - \alpha r} ){}^k ({2l + 1} )} \over {({k + 1} )!({2l + 5 + k} )}}},\quad j = 0 .\eqno(41) ]$

If n = l + 1 then equation (39) can be expressed as follows:

$[B(r) = {{ -({2l + 1})} \over {r({2l + 4} )}} + \alpha \sum\limits_{k = 0}^\infty {{{({ - \alpha r}){}^k ({k + 2} )({2l + 1} )} \over {({k + 1})!({2l + 5 + k} )}}}, \quad j = 1. \eqno(42) ]$

The terms r⁻¹ in (41) and (42) require special care when r = 0. They become infinite at the origin but are multiplied by angularly dependent factors which are zero at the origin. It turns out in these cases that the zeros are of higher order than infinity so their product is zero.

The first and partial second derivatives are then straightforward:

$[{{\partial G(r)} \over {\partial x}} = \left({{{{\rm d}G(r)} \over {{\rm d}r}}} \right)\left({{{\partial r} \over {{\rm d}x}}} \right) = xA(r). \eqno(43) ]$

Similarly,

$[{{\partial G(r)} \over {\partial y}} = yA(r) \eqno(44) ]$

and

$[{{\partial ^2 G(r)} \over {\partial x^2 }} = \left({{{\partial [{xA(r)}]} \over {\partial x}}} \right) = A(r) + x^2 B(r). \eqno(45) ]$

Similarly,

$[{{\partial ^2 G(r)} \over {\partial x\partial y}} = \left({{{\partial [{xA(r)} ]} \over {\partial y}}} \right) = xyB(r). \eqno(46) ]$

These formulae should be used for small αr.

At r = 0, derivatives are much simplified:

$[\eqalignno{G(0) &= {{({j + 1} )!} \over {\alpha ^{j + 2} }} &(47)\cr \left({{{\partial G(r)} \over {\partial x}}} \right)_0 &= \left({{{\partial G(r)} \over {\partial y}}} \right)_0 = \left({{{\partial G(r)} \over {\partial z}}} \right)_0 = 0 &(48)\cr \left({{{\partial ^2 G(r)} \over {\partial x^2 }}} \right)_0 &= \left({{{\partial ^2 G(r)} \over {\partial y^2 }}} \right)_0 = \left({{{\partial ^2 G(r)} \over {\partial z^2 }}} \right)_0 \cr&= A(0) = \cases { {{- ({2l + 1} )} \over {(2l + 3)}}& when $j = 0$ \cr 0 & when $j \hskip2.6pt\gt\hskip2.6pt 0$ } &(49) \cr \left({{{\partial G(r)} \over {\partial x\partial y}}} \right)_0 &= \left({{{\partial G(r)} \over {\partial x\partial z}}} \right)_0 = \left({{{\partial G(r)} \over {\partial y\partial z}}} \right)_0 = 0 .&(50)} ]$

4. Applications of the method

Calculations of electrostatic properties from the Hansen–Coppens pseudoatom formalism were performed using the newly derived formulae encoded in the new version of XDPROP, part of the XD package. To test the new formulae both theoretical and experimental data were used.

4.1. ESP, EF and EFG at nuclear positions in the formamide molecule from theoretical data

In a first example, a formamide molecule with a geometry extracted from the crystal of formamide (Stevens, 1978 ) was chosen. Theoretical calculations were performed with the Gaussian03 (2004 ) suite of programs at the density-functional (Hohenberg & Kohn, 1964 ) level of theory using the 1996 exchange and correlation functionals of Perdew, Burke & Ernzerhof (1996 , 1997 ) (PBE) and 6-31G** (Hariharan & Pople, 1973 ), aug-cc-pVDZ (Dunning, 1989 ; Kendall et al., 1992 ) and aug-cc-pVTZ (Dunning, 1989; Kendall et al., 1992) series of basis sets (labeled as PBE/6-31G**, PBE/aug-cc-pVDZ and PBE/aug-cc-pVTZ, respectively).

For the PBE/6-31G** calculation, complex static valence-only structure factors in the range of 0 < sinθ/λ < 1.1 Å⁻¹ were obtained by analytic Fourier transform of the molecular charge densities for reciprocal-lattice points corresponding to a pseudocubic cell with 30 Å edges. These data were fitted in terms of pseudoatom parameters as given in the Hansen–Coppens pseudoatom model [equation (8)] using the XD program suite (Koritsanszky et al., 2003 $[Koritsanszky, T., Howard, S. T., Richter, T., Macchi, P., Volkov, A., Gatti, C., Mallinson, P. R., Farrugia, L. J., Su, Z. & Hansen, N. K (2003). XD - a Computer Program Package for Multipole Refinement and Topological Analysis of Charge Densities from Diffraction Data. Middle Tennessee State University, TN, USA; University of Milano, Italy; University at Buffalo, NY, USA; CNR-ISTM, Milano, Italy; University of Glasgow, UK.]$ ). Phases of all reflections were reset to the theoretical values after each refinement cycle. Both radial screening factors (κ,κ′) were refined independently for each atom, with the exception of the chemically equivalent H atoms which shared the same κ and κ′ parameters. The multipolar expansion was truncated at the hexadecapolar level (l_max = 4) for the non-H atoms and at the quadrupolar level (l_max = 2) for H atoms, for which only bond-directed functions of l, m = 1, 0 and 2, 0 were refined. In order to reduce the number of least-squares variables, the following local-symmetry constraints were imposed: mm2 symmetry for N and m symmetry for O and C atoms. A molecular electroneutrality constraint was also applied. Refinement of valence-only structure factors yield an R factor of 3%, with a ratio of the number of calculated structure factors to the number of refined variables of 10502.

Tables 3–5 list the ESP, EF and EFG values at the nuclear positions as obtained directly from theoretical calculations and from the pseudoatom model, the latter labeled XD/PBE/6-31G**. Electrostatic properties at the nuclear positions from theoretical data were calculated with the Gaussian03 program. Agreement in all quantities is very good, taking into account the differences between Gaussian- and Slater-type functions and the fact that the projection of the Gaussian density onto the pseudoatom model is not perfect. Especially important is the fact that the asphericity of the EFG tensor at the nuclear positions of the H atoms is reproduced rather well despite small differences in each of the values. In general, the pseudoatom model gives slightly higher values of all properties at the nuclear positions compared to the Gaussian calculations. This is attributed to the different behavior of Gaussian- and Slater-type functions near r = 0, as well as the imperfect fit of the theoretical data.

Table 3
ESP and vector components of the EF at the nuclear positions in formamide (atomic units) from different methods

All components in this and the following tables refer to a global Cartesian coordinate system (atomic coordinates listed in Table 6, Appendix A)

		Electric field
	Potential	X	Y	Z
O(1)
PBE/6-31G**	−22.34	−0.096	−0.051	0.001
PBE/aug-cc-pVDZ	−22.35	−0.065	−0.033	0.001
PBE/aug-cc-pVTZ	−22.36	−0.031	−0.016	0.000
XD/PBE/6-31G**	−22.42	−0.274	−0.174	0.004
N(2)
PBE/6-31G**	−18.30	−0.002	−0.001	−0.001
PBE/aug-cc-pVDZ	−18.32	−0.004	0.001	−0.001
PBE/aug-cc-pVTZ	−18.32	−0.005	0.002	−0.001
XD/PBE/6-31G**	−18.44	−0.019	−0.009	0.000
C(3)
PBE/6-31G**	−14.65	−0.019	0.006	0.001
PBE/aug-cc-pVDZ	−14.66	−0.015	0.010	0.002
PBE/aug-cc-pVTZ	−14.67	−0.013	−0.001	0.001
XD/PBE/6-31G**	−14.76	−0.026	0.001	0.000
H(4)
PBE/6-31G**	−1.00	0.021	0.006	0.000
PBE/aug-cc-pVDZ	−0.98	0.027	0.011	0.000
PBE/aug-cc-pVTZ	−0.99	0.018	0.004	0.000
XD/PBE/6-31G**	−1.00	0.013	−0.010	0.000
H(5)
PBE/6-31G**	−1.00	0.005	−0.017	0.000
PBE/aug-cc-pVDZ	−0.98	0.006	−0.024	0.000
PBE/aug-cc-pVTZ	−0.99	0.005	−0.013	0.000
XD/PBE/6-31G**	−1.02	0.001	0.010	0.000
H(6)
PBE/6-31G**	−1.10	0.005	0.026	−0.001
PBE/aug-cc-pVDZ	−1.07	0.005	0.031	−0.001
PBE/aug-cc-pVTZ	−1.09	0.006	0.021	−0.001
XD/PBE/6-31G**	−1.09	0.003	0.022	−0.007

Table 4
Components of the EFG tensor at the nuclear positions in formamide (atomic units) from different methods

	XX	XY	XZ	YY	YZ	ZZ
O(1)
PBE/6-31G**	−1216	−1.02	0.02	−1215	0.015	−1216
PBE/aug-cc-pVDZ	−1246	−1.00	0.02	−1245	0.014	−1246
PBE/aug-cc-pVTZ	−1248	−1.00	0.02	−1247	0.014	−1249
XD/PBE/6-31G**	−1306	−0.80	−0.03	−1305	0.047	−1306
N(2)
PBE/6-31G**	−798	−0.023	−0.032	−798	−0.004	−797
PBE/aug-cc-pVDZ	−819	−0.024	−0.031	−819	−0.004	−818
PBE/aug-cc-pVTZ	−821	−0.018	−0.031	−821	−0.003	−819
XD/PBE/6-31G**	−863	−0.015	−0.012	−864	−0.013	−863
C(3)
PBE/6-31G**	−492	0.050	0.000	−492	0.001	−493
PBE/aug-cc-pVDZ	−505	0.053	0.000	−505	0.000	−506
PBE/aug-cc-pVTZ	−507	0.055	0.000	−507	0.001	−507
XD/PBE/6-31G**	−534	0.048	−0.002	−534	−0.012	−535
H(4)
PBE/6-31G**	−1.95	−0.28	−0.012	−1.78	−0.008	−1.49
PBE/aug-cc-pVDZ	−1.83	−0.28	−0.012	−1.66	−0.008	−1.38
PBE/aug-cc-pVTZ	−2.00	−0.28	−0.012	−1.83	−0.008	−1.55
XD/PBE/6-31G**	−2.29	−0.30	−0.012	−2.09	−0.009	−1.82
H(5)
PBE/6-31G**	−1.54	0.048	−0.002	−2.11	0.002	−1.47
PBE/aug-cc-pVDZ	−1.43	0.048	−0.001	−2.00	0.002	−1.37
PBE/aug-cc-pVTZ	−1.59	0.048	−0.001	−2.16	0.002	−1.53
XD/PBE/6-31G**	−1.87	0.049	−0.001	−2.46	0.002	−1.82
H(6)
PBE/6-31G**	−1.68	0.009	0.000	−2.11	−0.008	−1.66
PBE/aug-cc-pVDZ	−1.50	0.011	0.000	−1.96	−0.008	−1.48
PBE/aug-cc-pVTZ	−1.72	0.010	0.000	−2.15	−0.007	−1.71
XD/PBE/6-31G**	−2.00	0.007	−0.001	−2.45	−0.003	−1.99

Table 5
Eigenvalues of the EFG tensor at the nuclear positions in formamide (atomic units) from different methods

	λ₁	λ₂	λ₃
O(1)
PBE/6-31G**	−1217	−1216	−1214
PBE/aug-cc-pVDZ	−1247	−1246	−1244
PBE/aug-cc-pVTZ	−1249	−1249	−1246
XD/PBE/6-31G**	−1306	−1306	−1305
N(2)
PBE/6-31G**	−798	−798	−797
PBE/aug-cc-pVDZ	−819	−819	−818
PBE/aug-cc-pVTZ	−821	−821	−819
XD/PBE/6-31G**	−864	−863	−863
C(3)
PBE/6-31G**	−493	−492	−492
PBE/aug-cc-pVDZ	−506	−505	−505
PBE/aug-cc-pVTZ	−507	−507	−507
XD/PBE/6-31G**	−535	−534	−534
H(4)
PBE/6-31G**	−2.16	−1.57	−1.49
PBE/aug-cc-pVDZ	−2.04	−1.45	−1.38
PBE/aug-cc-pVTZ	−2.21	−1.62	−1.55
XD/PBE/6-31G**	−2.50	−1.88	−1.82
H(5)
PBE/6-31G**	−2.11	−1.54	−1.47
PBE/aug-cc-pVDZ	−2.00	−1.42	−1.37
PBE/aug-cc-pVTZ	−2.17	−1.59	−1.53
XD/PBE/6-31G**	−2.46	−1.87	−1.82
H(6)
PBE/6-31G**	−2.11	−1.68	−1.66
PBE/aug-cc-pVDZ	−1.96	−1.50	−1.48
PBE/aug-cc-pVTZ	−2.15	−1.72	−1.71
XD/PBE/6-31G**	−2.45	−2.00	−1.99

4.2. EF in the crystal structure of formamide from theoretical data

Note that, for an isolated molecule, the electric field at each nuclear position in the solid state must be zero once equilibrium is established. However, the molecules `extracted' from the crystal are not at equilibrium, therefore any calculation of the isolated molecule with the crystal geometry results in significant electric forces acting on the nuclei. These forces should be `compensated' in the crystal by those of the environment generated by the packing (of course after a final relaxation of the molecular electron density).

Fig. 1 shows the projection of the electric field in the O—C—N plane of the `central' formamide molecule due to the eight nearest-neighbor molecules as found in the crystal structure of formamide. The electric field was plotted with the program PlotMTV (Toh, 1995a ) using the MTV plot data format (Toh, 1995b ). The contribution of the `central' molecule to the electric field is omitted. In general, the external electric field is directed from the `negative' part of the central molecule (i.e. the O atom) towards the `positive' end (the amino group), and has nearly the same direction as the dipole moment of formamide. Indeed, molecules in crystals tend to orient themselves so as to achieve electrostatic stabilization. By the same token, the external electric field tends to polarize the central molecule so as to increase its dipole moment and enhance this electrostatic stabilization. Our own fully periodic (Saunders et al., 1998 ; Gatti, 1999 ) calculations at the B3LYP/6-31G** level of theory (Becke, 1988 ; Lee et al., 1988 ; Miehlich et al., 1989 ; Becke, 1993 ) show that the molecular dipole moment of formamide is increased to approximately 5.5 Debye in the solid state from about 4 Debye for the free molecule in the crystal geometry.

Figure 1
EF vectors in the N2—C3—O1 plane of the `central' formamide molecule due to the eight nearest neighbors in the crystal (the contribution of the `central' molecule to EF is not included). Subscripts of atom names identify the neighboring molecules. Vectors with magnitudes larger than 0.15 e Å⁻² are omitted for clarity. The size of the map is 6 × 6 Å with a grid spacing of 0.2 Å.

4.3. Determination of the ⁵⁷Fe^m nuclear quadrupole moment from experimental X-ray data

The nuclear quadrupole moment Q(⁵⁷Fe^m) cannot be directly measured because of the short lifetime of the excited nuclear state. However, Q(⁵⁷Fe^m) is directly related to the Mössbauer quadrupole splitting ΔE_Q:

$[\Delta E_Q = \textstyle{1 \over 2}\displaystyle eQ({{}^{57}{\rm Fe}^{m} } )V_{zz} \left({1 + {{\eta ^2 } \over 3}} \right)^{{1 / 2}}, \eqno(51) ]$

where e is the elementary charge, V_zz is the largest eigenvalue of the traceless EFG tensor

$[| {V_{zz} }| \,\gt\, | {V_{yy} }| \,\gt\, | {V_{xx} } | \eqno(52) ]$

and η is the asymmetry parameter defined as

$[\eta = {{V_{xx} - V_{yy} } \over {V_{zz} }}. \eqno(53) ]$

In order to obtain the three principal components of the EFG tensor defined by expression (5), the tensor must first be converted to its traceless form and then diagonalized. Note that ΔE_Q and Q(⁵⁷Fe^m) are usually given in Doppler-speed units of mm s⁻¹ and m², respectively, while the EFG tensor components are usually reported in e Å⁻³ or atomic units. The conversion between diffraction and spectroscopic units is discussed in detail by Coppens (1997, pp. 223–224).

Experimental X-ray data for iron pentacarbonyl Fe(CO)₅ were taken from the recent low-temperature (100 K) study of Farrugia & Evans (2005 ). Following the procedure of Su & Coppens (1996), the EFG tensor at the Fe nucleus is partitioned into its peripheral $[E_{\alpha \beta }^{\rm peripheral}]$ and central $[E_{\alpha \beta }^{\rm central}]$ contributions:

$[E_{\alpha \beta } ({\bf r}) = E_{\alpha \beta }^{\rm central} ({\bf r}) + E_{\alpha \beta }^{\rm peripheral} ({\bf r}). \eqno(54) ]$

The central component $[E_{\alpha \beta }^{\rm central}]$ includes only the contribution of the Fe pseudoatom, while the peripheral component $[E_{\alpha \beta }^{\rm peripheral}]$ includes both electronic and nuclear contributions of all other atoms. Because the Hansen–Coppens pseudoatom formalism uses a flexible valence shell but assumes a frozen core configuration, it is important to include the Sternheimer functions (Sternheimer, 1986 ) $[\gamma _\infty ^{\rm core}]$ and R^core in order to properly describe the shielding/antishielding of the EFG at the nuclear position due to the polarization induced in the atomic density by the quadrupolar components of the density distribution (Su & Coppens, 1996). Sternheimer functions are therefore included in expression (54) as

$[E_{\alpha \beta } ({\bf r}) = E_{\alpha \beta }^{\rm central} ({\bf r})({1 - R^{\rm core} }) + E_{\alpha \beta }^{\rm peripheral} ({\bf r})({1 - \gamma _\infty ^{\rm core} } ). \eqno(55) ]$

Values of R^core = 0.0730 and $[\gamma _\infty ^{\rm core}]$ = −8.933 as derived by Su & Coppens (1996) for the neutral Fe atom are used in the present study.

The recent experimental value of ΔE_Q = +2.51 mm s⁻¹ in Fe(CO)₅ was taken as the reference. This value agrees very well with the theoretical value of 2.54 mm s⁻¹ obtained from density-functional (B3LYP) calculations performed on the gas-phase optimized geometry of Fe(CO)₅ (Halvin et al., 1998 ; Zhang et al., 2002 ). The reference value gives Q(⁵⁷Fe^m) = 0.11 × 10⁻²⁸ m² without Sternheimer correction and 0.12 × 10⁻²⁸ m² after taking into account shielding/antishielding effects of the core. These values agree very well with those previously reported by Su & Coppens for iron pyrite FeS₂, sodium nitroprusside Na₂[Fe(NO)(CN)₅]·2H₂O and Fe(TPP)(pyridyl)₂, but slightly smaller than the earlier value of 0.14 (2) × 10⁻²⁸ m² determined (omitting shielding/antishielding effects of the core) by Tsirel'son et al. (1987 ), based on studies of sodium nitroprusside and Fe₂O₃. By averaging over three compounds, Su & Coppens obtained values of Q(⁵⁷Fe^m) = 0.09–0.10 × 10⁻²⁸ m² and 0.11–0.12 × 10⁻²⁸ m² from uncorrected and corrected calculations, respectively. When our value is included in the average, the mean corrected value of Q(⁵⁷Fe^m) becomes 0.12 × 10⁻²⁸ m² with a standard uncertainty of 0.02. It is within three standard uncertainties of the most precise up-to-date determination of Q(⁵⁷Fe^m) = 0.16 ± 5% × 10⁻²⁸ m² reported by Dufek et al. (1995 ) based on the comparison of spectroscopic values with EFG's from linearized augmented-plane-wave (LAPW) theoretical densities on a series of solids. It is interesting to note that all four recent experimental X-ray charge-density studies consistently show lower values of Q(⁵⁷Fe^m) than combined theoretical/spectroscopic studies. This discrepancy merits further study.

5. Concluding remarks

New practical formulae for calculation of the electrostatic potential (ESP), electric field (EF) and the electric field gradient (EFG) from the aspherical pseudoatom model are presented, which allow computation in regions near the nuclear centers. As real spherical harmonic density functions $[d_{lm \pm } (\Omega)]$ are discontinuous at the origin and thus non-differentiable, direct implementation of expressions containing $[d_{lm \pm } (\Omega)I({n,l,\alpha r})]$ in the evaluation of the electrostatic potential and its derivatives in this region is not possible. Instead the expressions have been reformulated in the form $[[{d_{lm \pm } (\Omega)r^l }] [{r^{ - l} I({n,l,\alpha r} )}]]$ , which eliminates this problem. Special care is required when treating the integral $[r^{ - l} I({n,l,\alpha r} )]$ when n = l or n = l + 1.

The expressions are applied to a theoretical density of formamide and to the derivation of the ⁵⁷Fe^m nuclear quadrupolar moment from experimental X-ray diffraction data. For formamide, the ESP, EF and EFG at the nuclear positions, calculated with the new expressions and a projection of PBE/6-31G** densities onto the Hansen−Coppens pseudoatom model, agree very well with theoretical values calculated directly from the wavefunction. Small differences observed are attributed to the different behavior of Slater- and Gaussian-type functions as r→0 and to imperfections in the fitting procedure.

The new expressions have further been applied in the detailed visualization of the electric field exerted on the `central' formamide molecule by the crystal environment. This was simulated by the electric fields of the eight closest neighboring molecules and omitting the contribution of the `central' molecule. The direction of the EF in the central molecule almost exactly coincides with the direction of the molecular dipole moment of the formamide molecule before its incorporation into the crystal, demonstrating the importance of the electrostatic forces in determining the crystal packing. The coincidence of the dipole moment and electric field directions provides a direct explanation for the enhancement of molecular dipole moments in crystals in accord with results of numerous experimental and theoretical studies.

Determination of the nuclear quadrupole moment of iron from the experimental X-ray diffraction data of Fe(CO)₅ yields a value of Q(⁵⁷Fe^m) = 0.12 × 10⁻²⁸ m², after taking into account shielding/antishielding effects of the core, which is in excellent agreement with previous X-ray studies by Su & Coppens (1996). However, this value is slightly smaller than the generally accepted value of 0.16 ± 5% × 10⁻²⁸ m² obtained from combined theoretical/spectroscopic studies (Dufek et al., 1995). The fact that X-ray determinations of Q(⁵⁷Fe^m) using different crystals and data sets consistently yield slightly lower values than those obtained from theoretical and spectroscopic studies requires further examination.

Note that the formulae presented for calculation of the ESP have already been used in our previous studies on the calculation of the electrostatic interaction energies in molecular crystals (Volkov et al., 2004 , 2006 ). Numerical quadrature evaluation of the two-centered Coulomb integral in the exact potential and multipole moment (EPMM) method requires an accurate yet efficient evaluation of the ESP at any r.

Application of the new method to topological analysis of the ESP, as recently performed by Bouhmaida et al. (2002 ), is currently being pursued.

APPENDIX A

Table 6
Cartesian coordinates (Å) of the atoms in formamide used in all calculations

Atom	x	y	z
O(1)	−1.1985900	−0.2390700	0.0036000
N(2)	1.0736400	−0.1651500	−0.0015000
C(3)	−0.1305900	0.3877800	−0.0088500
H(4)	1.8732000	0.4369500	0.0209100
H(5)	1.1409300	−1.1690400	0.0014700
H(6)	−0.1574100	1.4741700	0.0124500

Acknowledgements

Financial support of this work by the National Science Foundation (CHE0236317) is gratefully acknowledged.

References

Becke, A. D. (1988). Phys. Rev. A, 38, 3098–3100. CrossRef CAS PubMed Web of Science Google Scholar
Becke, A. D. (1993). J. Chem. Phys. 98, 5648–5652. CrossRef CAS Web of Science Google Scholar
Bertaut, E. F. (1978). J. Phys. Chem. Solids, 39, 97–102. CrossRef CAS Web of Science Google Scholar
Bouhmaida, N., Dutheil, M., Ghermani, N. E. & Becker, P. (2002) J. Chem. Phys. 116, 6196–6204. Web of Science CrossRef CAS Google Scholar
Bouhmaida, N., Ghermani, N. E., Lecomte, C. & Thalal A. (1997). Acta Cryst. A53, 556–563. CrossRef CAS Web of Science IUCr Journals Google Scholar
Brown, A. S. & Spackman, M. A. (1994). Mol. Phys. 83, 551–566. CrossRef CAS Web of Science Google Scholar
Clementi, E. & Roetti, C. (1974). At. Data Nucl. Data Tables, 14, 177–478. CrossRef CAS Google Scholar
Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press. Google Scholar
Dufek, P., Blaha, P. & Schwarz, K. (1995). Phys. Rev. Lett. 25, 3545–3548. CrossRef Web of Science Google Scholar
Dunning, T. H. Jr (1989). J. Chem. Phys. 90, 1007–1023. CrossRef CAS Web of Science Google Scholar
Edmonds, A. R. (1974). Angular Momentum in Quantum Mechanics. New Jersey: Princeton University Press. Google Scholar
Epstein, J. & Swanton, D. J. (1982). J. Chem. Phys. 77, 1048–1060. CrossRef CAS Web of Science Google Scholar
Farrugia, L. J. & Evans, C. (2005). J. Phys. Chem. A109, 8834–8848. CrossRef Google Scholar
Gatti, C. (1999). TOPOND98 Users' Manual, CNR-CSRSRC, Milano, Italy. Google Scholar
Gaussian 03 (2004). Gaussian 03, Revision C.02, by M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery Jr, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci et al. Gaussian, Inc., Wallingford, CT, USA. Google Scholar
Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1992). ELECTROS: Computer Program to Calculate Electrostatic Properties from High Resolution X-ray Diffraction. Internal Report URA CNRS 809. Université de Nancy I, France. Google Scholar
Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1993). Acta Cryst A49, 781–789. CrossRef CAS Web of Science IUCr Journals Google Scholar
Ghermani, N. E., Bouhmaida, N., Lecomte, C., Papet, A.-L. & Marsura, A. (1994). J. Phys. Chem. 98, 6287–6292. CrossRef CAS Web of Science Google Scholar
Ghermani, N. E., Lecomte, C. & Bouhmaida, N. (1993). Z. Naturforsch. 48, 91–98. CAS Google Scholar
Halvin, R. H., Godbout, N., Salzman, R., Wojdelski, M., Arnold, W., Schulz, C. E. & Oldfield, E. (1998). J. Am. Chem. Soc. 120, 3144–3151. Google Scholar
Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921. CrossRef CAS IUCr Journals Web of Science Google Scholar
Hariharan, P. C. & Pople, J. A. (1973). Theor. Chim Acta, 28, 213–222. CrossRef CAS Web of Science Google Scholar
He, X. M. (1983). The MOLPOT Computer Program. Tech. Rep., Department of Crystallography, University of Pittsburgh, PA, USA. Google Scholar
Hohenberg, P. & Kohn, W. (1964). Phys. Rev. 136, B864. CrossRef Web of Science Google Scholar
Jackson, J. D. (1975). Classical Electrodynamics. New York: John Wiley and Sons, Inc. Google Scholar
Kendall, R. A., Dunning, T. H. Jr & Harrison, R. J. (1992). J. Chem. Phys. 96, 6796–6806. CrossRef CAS Web of Science Google Scholar
Koritsanszky, T., Howard, S. T., Richter, T., Macchi, P., Volkov, A., Gatti, C., Mallinson, P. R., Farrugia, L. J., Su, Z. & Hansen, N. K (2003). XD – a Computer Program Package for Multipole Refinement and Topological Analysis of Charge Densities from Diffraction Data. Middle Tennessee State University, TN, USA; University of Milano, Italy; University at Buffalo, NY, USA; CNR-ISTM, Milano, Italy; University of Glasgow, UK. Google Scholar
Lecomte, C., Ghermani, N., Pichon-Pesme, V. & Souhassou, M. (1992). J. Mol. Struct. (Theochem), 255, 241–260. CrossRef Google Scholar
Lee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785–789. CrossRef CAS Web of Science Google Scholar
Miehlich, B., Savin, A., Stoll, H. & Preuss, H. (1989). Chem. Phys. Lett. 157, 200–206. CrossRef CAS Web of Science Google Scholar
Paturle, A. & Coppens, P. (1988). Acta Cryst. A44, 6–8. CrossRef CAS Web of Science IUCr Journals Google Scholar
Perdew, J. P., Burke, K. & Ernzerhof, M. (1996). Phys. Rev. Lett. 77, 3865–3868. Web of Science CrossRef PubMed CAS Google Scholar
Perdew, J. P., Burke, K. & Ernzerhof, M. (1997). Phys. Rev. Lett. 78, 1396–1396. CrossRef CAS Web of Science Google Scholar
Saunders, V. R., Dovesi, R., Roetti, C., Causà, M., Harrison, N. M., Orlando, R. & Zicovich-Wilson, C. M. (1998). CRYSTAL98 User's Manual, University of Torino, Italy. Google Scholar
Schwarzenbach, D. & Thong, N. (1979). Acta Cryst. A35, 652–658. CrossRef CAS IUCr Journals Web of Science Google Scholar
Slater, J. C. (1932). Phys. Rev. 42, 33–43. CrossRef CAS Google Scholar
Sternheimer, R. M. (1986). Z. Naturforsch. Teil A, 41, 24–36. Google Scholar
Stevens, E. D. (1978). Acta Cryst. B34, 544–551. CSD CrossRef CAS IUCr Journals Web of Science Google Scholar
Stewart, R. F. (1979). Chem. Phys. Lett. 65, 335–342. CrossRef CAS Web of Science Google Scholar
Stewart, R. F., Spackman, M. A. & Flensburg, C. (1998). VALRAY98, Users Manual, Carnegie Mellon University, Pittsburgh, PA, USA, and University of Copenhagen, Denmark. Google Scholar
Su, Z. & Coppens, P. (1992). Acta Cryst. A48, 188–197. CrossRef CAS Web of Science IUCr Journals Google Scholar
Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748–756. CrossRef CAS Web of Science IUCr Journals Google Scholar
Toh, K. K. H. (1995a). PlotMTV – Fast Multi-Purpose Plotting Program for X11-Windows. Google Scholar
Toh, K. K. H. (1995b). MTV Plot Data Format, Version 1.4.1, Rev. 0. Google Scholar
Tsirel'son, V. G., Strel'tsov, V. A., Makarov, E. F. & Ozerov, R. P. (1987). Sov. Phys. JETP, 65, 1065–1069. Google Scholar
Velde, B. te (1990). BAND – a Fortran Program for Band Structure Calculations. PhD thesis, Vrije Universiteit, Amsterdam, The Netherlands. Google Scholar
Velde, B. te, Bickelhaupt, F. M., van Gisbergen, S. J. A., Fonseca Guerra, C., Baerends, E. J., Snijders, J. G. & Ziegler, T. J. (2001). J. Comput. Chem. 22, 931–967. Web of Science CrossRef Google Scholar
Volkov, A., King, H. F. & Coppens, P. (2006). J. Chem. Theory Comput. 2, 81–89. Web of Science CrossRef CAS Google Scholar
Volkov, A., Koritsanszky T. S. & Coppens, P. (2004). Chem. Phys. Lett. 391, 170–175. Google Scholar
Zhang, Y., Mao, J., Godbout, N. & Oldfield, E. (2002). J. Am. Chem. Soc. 124, 13921–13930. Web of Science CrossRef PubMed CAS Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 62| Part 5| September 2006| Pages 400-408

doi:10.1107/S0108767306026298

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search term		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

On the calculation of the electrostatic potential, electric field and electric field gradient from the aspherical pseudoatom model

1. Introduction

2. Calculation of the electronic potential and its derivatives for Slater-type functions

3. Calculation of the electronic potential and its derivatives near r = 0

4. Applications of the method

4.1. ESP, EF and EFG at nuclear positions in the formamide molecule from theoretical data

4.2. EF in the crystal structure of formamide from theoretical data

4.3. Determination of the 57Fem nuclear quadrupole moment from experimental X-ray data

5. Concluding remarks

APPENDIX A

Acknowledgements

References

research papers

4.3. Determination of the ⁵⁷Fe^m nuclear quadrupole moment from experimental X-ray data