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

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733
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

CROSSMARK_Color_square_no_text.svg

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[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). Acta Cryst. A34, 909–921]. They are based on the expansion of |r′ − r|−1 in spherical harmonics and the incomplete gamma function for a Slater-type function of the form Rl(r) = rn 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)5. 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)5 yields the value of the nuclear quadrupole moment Q(57Fem) = 0.12 × 10−28 m2 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[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]), A52, 748–756] but slightly smaller than the generally accepted value of 0.16 ± 5% × 10−28 m2 obtained from combined theoretical/spectroscopic studies [Dufek, Blaha & Schwarz (1995). Phys. Rev. Lett. 25, 3545–3548].

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[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.]):

[\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 Velec(r) and nuclear Vnuc(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 Ri carrying positive charge Zi. 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 Ri) 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 Velec(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[Bertaut, E. F. (1978). J. Phys. Chem. Solids, 39, 97-102.]; Stewart, 1979[Stewart, R. F. (1979). Chem. Phys. Lett. 65, 335-342.]; Schwarzenbach & Thong, 1979[Schwarzenbach, D. & Thong, N. (1979). Acta Cryst. A35, 652-658.]) 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[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.]; Hansen & Coppens, 1978[Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.]). 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[Clementi, E. & Roetti, C. (1974). At. Data Nucl. Data Tables, 14, 177-478.]). The population of the core, Pcore, is always frozen, while the population of the spherical valence shell Pval is refined together with the expansion–contraction parameter κ. The third term describes the aspherical deformation density. Coefficients Plm± are the population parameters and κ′ are the dimensionless adjustment coefficients of the radial functions Rl. In equation (8)[link], r is the distance from the pseudoatom center, r = |r − R|, and Ω are the corresponding angular coordinates. The angular factor dlm± is a real density-normalized spherical harmonic:

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

where m ≥ 0 and ylm± 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 Dlm are given by Paturle & Coppens (1988[Paturle, A. & Coppens, P. (1988). Acta Cryst. A44, 6-8.]) and Coppens (1997[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.]).

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[Brown, A. S. & Spackman, M. A. (1994). Mol. Phys. 83, 551-566.]).

The method of Su & Coppens (1992[Su, Z. & Coppens, P. (1992). Acta Cryst. A48, 188-197.], 1996[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]) is based on the Fourier convolution theorem previously applied by Epstein & Swanton (1982[Epstein, J. & Swanton, D. J. (1982). J. Chem. Phys. 77, 1048-1060.]) 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 AN,l1, l2, 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[Brown, A. S. & Spackman, M. A. (1994). Mol. Phys. 83, 551-566.]). 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[Stewart, R. F., Spackman, M. A. & Flensburg, C. (1998). VALRAY98, Users Manual, Carnegie Mellon University, Pittsburgh, PA, USA, and University of Copenhagen, Denmark.]). 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[Lecomte, C., Ghermani, N., Pichon-Pesme, V. & Souhassou, M. (1992). J. Mol. Struct. (Theochem), 255, 241-260.]; Ghermani, Bouhmaida & Lecomte, 1993[Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1993). Acta Cryst A49, 781-789.]; Ghermani, Lecomte & Bouhmaida, 1993[Ghermani, N. E., Lecomte, C. & Bouhmaida, N. (1993). Z. Naturforsch. 48, 91-98.]; Ghermani et al., 1994[Ghermani, N. E., Bouhmaida, N., Lecomte, C., Papet, A.-L. & Marsura, A. (1994). J. Phys. Chem. 98, 6287-6292.]; Bouhmaida et al., 1997[Bouhmaida, N., Ghermani, N. E., Lecomte, C. & Thalal A. (1997). Acta Cryst. A53, 556-563.]) reported expressions for the ESP derived directly from the Hansen–Coppens density model [equation (8[link])], in which the ESP due to a pseudoatom is expanded in the same way as the density ρat(r) itself [equation (8[link])], i.e.

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

Here Vcore(r), Vval(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[link])]. Note that, in expression (12)[link], the contribution of the nuclear potential is explicitly combined with the core-electron contribution. Formula (14)[link] 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[He, X. M. (1983). The MOLPOT Computer Program. Tech. Rep., Department of Crystallography, University of Pittsburgh, PA, USA.]). Numerous studies have been published using this approach.

In the current paper, we review the derivation of equation (14)[link], 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 57Fem 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)[link] can be reduced to a linear combination of Slater-type (Slater, 1932[Slater, J. C. (1932). Phys. Rev. 42, 33-43.]) 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[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.]). 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[Jackson, J. D. (1975). Classical Electrodynamics. New York: John Wiley and Sons, Inc.]):

[{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 ylm± are obtained by unitary transformation of the Ylm, which implies that the form of the spherical harmonic addition theorem (Edmonds, 1974[Edmonds, A. R. (1974). Angular Momentum in Quantum Mechanics. New Jersey: Princeton University Press.]) 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)[link] and (17)[link] are inserted into (16)[link], 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 ylm± 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)[link] 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)[link] 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[Lecomte, C., Ghermani, N., Pichon-Pesme, V. & Souhassou, M. (1992). J. Mol. Struct. (Theochem), 255, 241-260.]; 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, N. E., Lecomte, C. & Bouhmaida, N. (1993). Z. Naturforsch. 48, 91-98.]; Ghermani, Bouhmaida & Lecomte, 1993[Ghermani, N. E., Bouhmaida, N. & Lecomte, C. (1993). Acta Cryst A49, 781-789.]; Ghermani et al., 1994[Ghermani, N. E., Bouhmaida, N., Lecomte, C., Papet, A.-L. & Marsura, A. (1994). J. Phys. Chem. 98, 6287-6292.]; Bouhmaida et al., 1997[Bouhmaida, N., Ghermani, N. E., Lecomte, C. & Thalal A. (1997). Acta Cryst. A53, 556-563.]). Integrals I1(n,l,αr) and I2(n,l,αr) have the form of the incomplete gamma function (te Velde, 1990[Velde, B. te (1990). BAND - a Fortran Program for Band Structure Calculations. PhD thesis, Vrije Universiteit, Amsterdam, The Netherlands.]; te Velde et al., 2001[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.]):

[\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)[link] and (23)[link] represent the effects of the interpenetration of the charge-density distributions. The remaining term in equation (22)[link] 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[link]. 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 I1(n,l,αr), and can be avoided by rearranging equation (22)[link]:

[\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[link]. Thus, expression (24)[link] is recommended for use for `small' values of αr while expression (22)[link] 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)[link] and (23)[link]

  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)[link] and (25)[link] 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[link] 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)[link] 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[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.]). 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)[link] can still be applied.

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

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

Functions dlm±(Ω) and their derivatives are poorly behaved near the origin, as they contain direction cosines, so we remove a factor rl 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)[link] contain powers of αr of l or higher. For l ≤ 4, derivatives of dlm±(Ω)rl functions are simple and well behaved at any r. For example, the second partial deriv­ative of d41−(Ω)r4 w.r.t. xy is L41−(−6xz), where L41− = 0.474 is the density function normalization factor (Paturle & Coppens, 1988[Paturle, A. & Coppens, P. (1988). Acta Cryst. A44, 6-8.]). No second derivative is more complicated than that of d40(Ω)r4 w.r.t. xx, which is simply L40(36x2 + 12y2 − 48z2). Note that first derivatives of dlm±(Ω)rl 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 dlm±(Ω)rl 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[link] 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)[link] and (35)[link] both have non-zero terms in r3, which creates singularities in the second derivatives d2r3/dq2, 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)[link] 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)[link] and (39)[link] 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)[link] 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−1 in (41)[link] and (42)[link] 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[Stevens, E. D. (1978). Acta Cryst. B34, 544-551.]) was chosen. Theoretical calculations were performed with the Gaussian03 (2004[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.]) suite of programs at the density-functional (Hohenberg & Kohn, 1964[Hohenberg, P. & Kohn, W. (1964). Phys. Rev. 136, B864.]) level of theory using the 1996 exchange and correlation functionals of Perdew, Burke & Ernzerhof (1996[Perdew, J. P., Burke, K. & Ernzerhof, M. (1996). Phys. Rev. Lett. 77, 3865-3868.], 1997[Perdew, J. P., Burke, K. & Ernzerhof, M. (1997). Phys. Rev. Lett. 78, 1396-1396.]) (PBE) and 6-31G** (Hariharan & Pople, 1973[Hariharan, P. C. & Pople, J. A. (1973). Theor. Chim Acta, 28, 213-222.]), aug-cc-pVDZ (Dunning, 1989[Dunning, T. H. Jr (1989). J. Chem. Phys. 90, 1007-1023.]; Kendall et al., 1992[Kendall, R. A., Dunning, T. H. Jr & Harrison, R. J. (1992). J. Chem. Phys. 96, 6796-6806.]) and aug-cc-pVTZ (Dunning, 1989[Dunning, T. H. Jr (1989). J. Chem. Phys. 90, 1007-1023.]; Kendall et al., 1992[Kendall, R. A., Dunning, T. H. Jr & Harrison, R. J. (1992). J. Chem. Phys. 96, 6796-6806.]) 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 Å−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[link])] 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 κ′ param­eters. The multipolar expansion was truncated at the hexadecapolar level (lmax = 4) for the non-H atoms and at the quadrupolar level (lmax = 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[link][link][link] 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 pseudo­atom 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[link], Appendix A[link])

    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

  λ1 λ2 λ3
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 equi­librium 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 en­viron­ment generated by the packing (of course after a final relaxation of the molecular electron density).

Fig. 1[link] 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[Toh, K. K. H. (1995a). PlotMTV - Fast Multi-Purpose Plotting Program for X11-Windows.]) using the MTV plot data format (Toh, 1995b[Toh, K. K. H. (1995b). MTV Plot Data Format, Version 1.4.1, Rev. 0. ]). The contribution of the `central' mol­ecule 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[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.]; Gatti, 1999[Gatti, C. (1999). TOPOND98 Users' Manual, CNR-CSRSRC, Milano, Italy.]) calculations at the B3LYP/6-31G** level of theory (Becke, 1988[Becke, A. D. (1988). Phys. Rev. A, 38, 3098-3100.]; Lee et al., 1988[Lee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785-789.]; Miehlich et al., 1989[Miehlich, B., Savin, A., Stoll, H. & Preuss, H. (1989). Chem. Phys. Lett. 157, 200-206.]; Becke, 1993[Becke, A. D. (1993). J. Chem. Phys. 98, 5648-5652.]) 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]
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 Å−2 are omitted for clarity. The size of the map is 6 × 6 Å with a grid spacing of 0.2 Å.

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

The nuclear quadrupole moment Q(57Fem) cannot be directly measured because of the short lifetime of the excited nuclear state. However, Q(57Fem) is directly related to the Mössbauer quadrupole splitting ΔEQ:

[\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, Vzz 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)[link], the tensor must first be converted to its traceless form and then diagonalized. Note that ΔEQ and Q(57Fem) are usually given in Doppler-speed units of mm s−1 and m2, respectively, while the EFG tensor components are usually reported in e Å−3 or atomic units. The conversion between diffraction and spectroscopic units is discussed in detail by Coppens (1997[Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.], pp. 223–224).

Experimental X-ray data for iron pentacarbonyl Fe(CO)5 were taken from the recent low-temperature (100 K) study of Farrugia & Evans (2005[Farrugia, L. J. & Evans, C. (2005). J. Phys. Chem. A109, 8834-8848.]). Following the procedure of Su & Coppens (1996[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]), 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[Sternheimer, R. M. (1986). Z. Naturforsch. Teil A, 41, 24-36.]) [\gamma _\infty ^{\rm core}] and Rcore 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[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]). Sternheimer functions are therefore included in expression (54)[link] 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 Rcore = 0.0730 and [\gamma _\infty ^{\rm core}] = −8.933 as derived by Su & Coppens (1996[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]) for the neutral Fe atom are used in the present study.

The recent experimental value of ΔEQ = +2.51 mm s−1 in Fe(CO)5 was taken as the reference. This value agrees very well with the theoretical value of 2.54 mm s−1 obtained from density-functional (B3LYP) calculations performed on the gas-phase optimized geometry of Fe(CO)5 (Halvin et al., 1998[Halvin, R. H., Godbout, N., Salzman, R., Wojdelski, M., Arnold, W., Schulz, C. E. & Oldfield, E. (1998). J. Am. Chem. Soc. 120, 3144-3151. ]; Zhang et al., 2002[Zhang, Y., Mao, J., Godbout, N. & Oldfield, E. (2002). J. Am. Chem. Soc. 124, 13921-13930.]). The reference value gives Q(57Fem) = 0.11 × 10−28 m2 without Sternheimer correction and 0.12 × 10−28 m2 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 FeS2, sodium nitroprusside Na2[Fe(NO)(CN)5]·2H2O and Fe(TPP)(pyridyl)2, but slightly smaller than the earlier value of 0.14 (2) × 10−28 m2 determined (omitting shielding/antishielding effects of the core) by Tsirel'son et al. (1987[Tsirel'son, V. G., Strel'tsov, V. A., Makarov, E. F. & Ozerov, R. P. (1987). Sov. Phys. JETP, 65, 1065-1069.]), based on studies of sodium nitroprusside and Fe2O3. By averaging over three compounds, Su & Coppens obtained values of Q(57Fem) = 0.09–0.10 × 10−28 m2 and 0.11–0.12 × 10−28 m2 from uncorrected and corrected calculations, respectively. When our value is included in the average, the mean corrected value of Q(57Fem) becomes 0.12 × 10−28 m2 with a standard uncertainty of 0.02. It is within three standard uncertainties of the most precise up-to-date determination of Q(57Fem) = 0.16 ± 5% × 10−28 m2 reported by Dufek et al. (1995[Dufek, P., Blaha, P. & Schwarz, K. (1995). Phys. Rev. Lett. 25, 3545-3548.]) 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(57Fem) 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 57Fem 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)5 yields a value of Q(57Fem) = 0.12 × 10−28 m2, after taking into account shielding/antishielding effects of the core, which is in excellent agreement with previous X-ray studies by Su & Coppens (1996[Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.]). However, this value is slightly smaller than the generally accepted value of 0.16 ± 5% × 10−28 m2 obtained from combined theoretical/spectroscopic studies (Dufek et al., 1995[Dufek, P., Blaha, P. & Schwarz, K. (1995). Phys. Rev. Lett. 25, 3545-3548.]). The fact that X-ray determinations of Q(57Fem) using different crystals and data sets consistently yield slightly lower values than those obtained from theor­etical 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[Volkov, A., Koritsanszky T. S. & Coppens, P. (2004). Chem. Phys. Lett. 391, 170-175.], 2006[Volkov, A., King, H. F. & Coppens, P. (2006). J. Chem. Theory Comput. 2, 81-89.]). 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[Bouhmaida, N., Dutheil, M., Ghermani, N. E. & Becker, P. (2002) J. Chem. Phys. 116, 6196-6204.]), is currently being pursued.

APPENDIX A

[link]

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

First citationBecke, A. D. (1988). Phys. Rev. A, 38, 3098–3100.  CrossRef CAS PubMed Web of Science Google Scholar
First citationBecke, A. D. (1993). J. Chem. Phys. 98, 5648–5652.  CrossRef CAS Web of Science Google Scholar
First citationBertaut, E. F. (1978). J. Phys. Chem. Solids, 39, 97–102.  CrossRef CAS Web of Science Google Scholar
First citationBouhmaida, N., Dutheil, M., Ghermani, N. E. & Becker, P. (2002) J. Chem. Phys. 116, 6196–6204.  Web of Science CrossRef CAS Google Scholar
First citationBouhmaida, N., Ghermani, N. E., Lecomte, C. & Thalal A. (1997). Acta Cryst. A53, 556–563.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationBrown, A. S. & Spackman, M. A. (1994). Mol. Phys. 83, 551–566.  CrossRef CAS Web of Science Google Scholar
First citationClementi, E. & Roetti, C. (1974). At. Data Nucl. Data Tables, 14, 177–478.  CrossRef CAS Google Scholar
First citationCoppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.  Google Scholar
First citationDufek, P., Blaha, P. & Schwarz, K. (1995). Phys. Rev. Lett. 25, 3545–3548.  CrossRef Web of Science Google Scholar
First citationDunning, T. H. Jr (1989). J. Chem. Phys. 90, 1007–1023.  CrossRef CAS Web of Science Google Scholar
First citationEdmonds, A. R. (1974). Angular Momentum in Quantum Mechanics. New Jersey: Princeton University Press.  Google Scholar
First citationEpstein, J. & Swanton, D. J. (1982). J. Chem. Phys. 77, 1048–1060.  CrossRef CAS Web of Science Google Scholar
First citationFarrugia, L. J. & Evans, C. (2005). J. Phys. Chem. A109, 8834–8848.  CrossRef Google Scholar
First citationGatti, C. (1999). TOPOND98 Users' Manual, CNR-CSRSRC, Milano, Italy.  Google Scholar
First citationGaussian 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
First citationGhermani, 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
First citationGhermani, N. E., Bouhmaida, N. & Lecomte, C. (1993). Acta Cryst A49, 781–789.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationGhermani, 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
First citationGhermani, N. E., Lecomte, C. & Bouhmaida, N. (1993). Z. Naturforsch. 48, 91–98.  CAS Google Scholar
First citationHalvin, 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
First citationHansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909–921.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationHariharan, P. C. & Pople, J. A. (1973). Theor. Chim Acta, 28, 213–222.  CrossRef CAS Web of Science Google Scholar
First citationHe, X. M. (1983). The MOLPOT Computer Program. Tech. Rep., Department of Crystallography, University of Pittsburgh, PA, USA.  Google Scholar
First citationHohenberg, P. & Kohn, W. (1964). Phys. Rev. 136, B864.  CrossRef Web of Science Google Scholar
First citationJackson, J. D. (1975). Classical Electrodynamics. New York: John Wiley and Sons, Inc.  Google Scholar
First citationKendall, R. A., Dunning, T. H. Jr & Harrison, R. J. (1992). J. Chem. Phys. 96, 6796–6806.  CrossRef CAS Web of Science Google Scholar
First citationKoritsanszky, 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
First citationLecomte, C., Ghermani, N., Pichon-Pesme, V. & Souhassou, M. (1992). J. Mol. Struct. (Theochem), 255, 241–260.  CrossRef Google Scholar
First citationLee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785–789.  CrossRef CAS Web of Science Google Scholar
First citationMiehlich, B., Savin, A., Stoll, H. & Preuss, H. (1989). Chem. Phys. Lett. 157, 200–206.  CrossRef CAS Web of Science Google Scholar
First citationPaturle, A. & Coppens, P. (1988). Acta Cryst. A44, 6–8.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationPerdew, J. P., Burke, K. & Ernzerhof, M. (1996). Phys. Rev. Lett. 77, 3865–3868.  Web of Science CrossRef PubMed CAS Google Scholar
First citationPerdew, J. P., Burke, K. & Ernzerhof, M. (1997). Phys. Rev. Lett. 78, 1396–1396.  CrossRef CAS Web of Science Google Scholar
First citationSaunders, 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
First citationSchwarzenbach, D. & Thong, N. (1979). Acta Cryst. A35, 652–658.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationSlater, J. C. (1932). Phys. Rev. 42, 33–43.  CrossRef CAS Google Scholar
First citationSternheimer, R. M. (1986). Z. Naturforsch. Teil A, 41, 24–36.  Google Scholar
First citationStevens, E. D. (1978). Acta Cryst. B34, 544–551.  CSD CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationStewart, R. F. (1979). Chem. Phys. Lett. 65, 335–342.  CrossRef CAS Web of Science Google Scholar
First citationStewart, R. F., Spackman, M. A. & Flensburg, C. (1998). VALRAY98, Users Manual, Carnegie Mellon University, Pittsburgh, PA, USA, and University of Copenhagen, Denmark.  Google Scholar
First citationSu, Z. & Coppens, P. (1992). Acta Cryst. A48, 188–197.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationSu, Z. & Coppens, P. (1996). Acta Cryst. A52, 748–756.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationToh, K. K. H. (1995a). PlotMTV – Fast Multi-Purpose Plotting Program for X11-Windows.  Google Scholar
First citationToh, K. K. H. (1995b). MTV Plot Data Format, Version 1.4.1, Rev. 0.  Google Scholar
First citationTsirel'son, V. G., Strel'tsov, V. A., Makarov, E. F. & Ozerov, R. P. (1987). Sov. Phys. JETP, 65, 1065–1069.  Google Scholar
First citationVelde, B. te (1990). BAND – a Fortran Program for Band Structure Calculations. PhD thesis, Vrije Universiteit, Amsterdam, The Netherlands.  Google Scholar
First citationVelde, 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
First citationVolkov, A., King, H. F. & Coppens, P. (2006). J. Chem. Theory Comput. 2, 81–89.  Web of Science CrossRef CAS Google Scholar
First citationVolkov, A., Koritsanszky T. S. & Coppens, P. (2004). Chem. Phys. Lett. 391, 170–175.  Google Scholar
First citationZhang, 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.

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733
Follow Acta Cryst. A
Sign up for e-alerts
Follow Acta Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds