Received 26 May 2006
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
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|-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), 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].
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):
where total includes both electronic and nuclear charges:
and (r) is the electron density. For a continuous charge distribution, the potential is obtained by integration of the density over all space:
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):
while the negative of the second derivative of the ESP is the electric field gradient (EFG):
Just as the total density total(r), the ESP in the molecule can be separated into electronic Velec(r) and nuclear Vnuc(r) contributions:
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 , both terms need to be taken into account, while, at , 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; 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:
Each pseudoatom is modeled using the modified Laplace series:
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, 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), 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:
where m 0 and ylm± is the normalized linear combination of complex spherical harmonics:
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). 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). 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.
Here Vcore(r), Vval(r) and V(r) are the spherical core, spherical valence and aspherical deformation contributions, respectively, and are given by
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) 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 57Fem nuclear quadrupolar moment from experimental X-ray diffraction data.
where r< is the smaller and r> is the greater of and . 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) is preserved, i.e.
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 ylm± functions. Thus,
Essentially the same formula is given by Ghermani, Bouhmaida, Lecomte and co-workers (Lecomte et al., 1992; Ghermani et al., 1992; Ghermani, Lecomte & Bouhmaida, 1993; Ghermani, Bouhmaida & Lecomte, 1993; Ghermani et al., 1994; Bouhmaida et al., 1997). Integrals I1(n,l,r) and I2(n,l,r) have the form of the incomplete gamma function (te Velde, 1990; te Velde et al., 2001):
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 I1(n,l,r), and can be avoided by rearranging equation (22):
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.
Alternatively, integral I(n,l,r) can be replaced with its Taylor-series expansion:
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:
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:
Note that all terms in (25) 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 derivative of d41-()r4 w.r.t. xy is L41-(-6xz), where L41- = 0.474 is the density function normalization factor (Paturle & Coppens, 1988). 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
The second derivatives of dlm±()rl vanish at r = 0 for all functions except for
We define function G(r) as follows:
To simplify the notation of G, the dependence on the n and l indices is implicitly understood. Let
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:
For example, when n = l then
When n = l + 1, then
Define A(r) as follows:
Define B(r) as follows:
The terms r-1 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:
These formulae should be used for small r.
At r = 0, derivatives are much simplified:
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.
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 Å-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). 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 (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 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.
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 Å-2 are omitted for clarity. The size of the map is 6 × 6 Å with a grid spacing of 0.2 Å.
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:
where e is the elementary charge, Vzz is the largest eigenvalue of the traceless EFG tensor
and is the asymmetry parameter defined as
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 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, 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). Following the procedure of Su & Coppens (1996), the EFG tensor at the Fe nucleus is partitioned into its peripheral and central contributions:
The central component includes only the contribution of the Fe pseudoatom, while the peripheral component 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) 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). Sternheimer functions are therefore included in expression (54) as
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; Zhang et al., 2002). 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), 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) 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.
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 are discontinuous at the origin and thus non-differentiable, direct implementation of expressions containing 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 , which eliminates this problem. Special care is required when treating the integral 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 r0 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). 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). The fact that X-ray determinations of Q(57Fem) 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.
Financial support of this work by the National Science Foundation (CHE0236317) is gratefully acknowledged.
Becke, A. D. (1988). Phys. Rev. A, 38, 3098-3100.
Becke, A. D. (1993). J. Chem. Phys. 98, 5648-5652.
Bertaut, E. F. (1978). J. Phys. Chem. Solids, 39, 97-102.
Bouhmaida, N., Dutheil, M., Ghermani, N. E. & Becker, P. (2002) J. Chem. Phys. 116, 6196-6204.
Bouhmaida, N., Ghermani, N. E., Lecomte, C. & Thalal A. (1997). Acta Cryst. A53, 556-563.
Brown, A. S. & Spackman, M. A. (1994). Mol. Phys. 83, 551-566.
Clementi, E. & Roetti, C. (1974). At. Data Nucl. Data Tables, 14, 177-478.
Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.
Dufek, P., Blaha, P. & Schwarz, K. (1995). Phys. Rev. Lett. 25, 3545-3548.
Dunning, T. H. Jr (1989). J. Chem. Phys. 90, 1007-1023.
Edmonds, A. R. (1974). Angular Momentum in Quantum Mechanics. New Jersey: Princeton University Press.
Epstein, J. & Swanton, D. J. (1982). J. Chem. Phys. 77, 1048-1060.
Farrugia, L. J. & Evans, C. (2005). J. Phys. Chem. A109, 8834-8848.
Gatti, C. (1999). TOPOND98 Users' Manual, CNR-CSRSRC, Milano, Italy.
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.
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, N. E., Bouhmaida, N. & Lecomte, C. (1993). Acta Cryst A49, 781-789.
Ghermani, N. E., Bouhmaida, N., Lecomte, C., Papet, A.-L. & Marsura, A. (1994). J. Phys. Chem. 98, 6287-6292.
Ghermani, N. E., Lecomte, C. & Bouhmaida, N. (1993). Z. Naturforsch. 48, 91-98.
Halvin, R. H., Godbout, N., Salzman, R., Wojdelski, M., Arnold, W., Schulz, C. E. & Oldfield, E. (1998). J. Am. Chem. Soc. 120, 3144-3151.
Hansen, N. K. & Coppens, P. (1978). Acta Cryst. A34, 909-921.
Hariharan, P. C. & Pople, J. A. (1973). Theor. Chim Acta, 28, 213-222.
He, X. M. (1983). The MOLPOT Computer Program. Tech. Rep., Department of Crystallography, University of Pittsburgh, PA, USA.
Hohenberg, P. & Kohn, W. (1964). Phys. Rev. 136, B864.
Jackson, J. D. (1975). Classical Electrodynamics. New York: John Wiley and Sons, Inc.
Kendall, R. A., Dunning, T. H. Jr & Harrison, R. J. (1992). J. Chem. Phys. 96, 6796-6806.
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.
Lecomte, C., Ghermani, N., Pichon-Pesme, V. & Souhassou, M. (1992). J. Mol. Struct. (Theochem), 255, 241-260.
Lee, C., Yang, W. & Parr, R. G. (1988). Phys. Rev. B, 37, 785-789.
Miehlich, B., Savin, A., Stoll, H. & Preuss, H. (1989). Chem. Phys. Lett. 157, 200-206.
Paturle, A. & Coppens, P. (1988). Acta Cryst. A44, 6-8.
Perdew, J. P., Burke, K. & Ernzerhof, M. (1996). Phys. Rev. Lett. 77, 3865-3868.
Perdew, J. P., Burke, K. & Ernzerhof, M. (1997). Phys. Rev. Lett. 78, 1396-1396.
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.
Schwarzenbach, D. & Thong, N. (1979). Acta Cryst. A35, 652-658.
Slater, J. C. (1932). Phys. Rev. 42, 33-43.
Sternheimer, R. M. (1986). Z. Naturforsch. Teil A, 41, 24-36.
Stevens, E. D. (1978). Acta Cryst. B34, 544-551.
Stewart, R. F. (1979). Chem. Phys. Lett. 65, 335-342.
Stewart, R. F., Spackman, M. A. & Flensburg, C. (1998). VALRAY98, Users Manual, Carnegie Mellon University, Pittsburgh, PA, USA, and University of Copenhagen, Denmark.
Su, Z. & Coppens, P. (1992). Acta Cryst. A48, 188-197.
Su, Z. & Coppens, P. (1996). Acta Cryst. A52, 748-756.
Toh, K. K. H. (1995a). PlotMTV - Fast Multi-Purpose Plotting Program for X11-Windows.
Toh, K. K. H. (1995b). MTV Plot Data Format, Version 1.4.1, Rev. 0.
Tsirel'son, V. G., Strel'tsov, V. A., Makarov, E. F. & Ozerov, R. P. (1987). Sov. Phys. JETP, 65, 1065-1069.
Velde, B. te (1990). BAND - a Fortran Program for Band Structure Calculations. PhD thesis, Vrije Universiteit, Amsterdam, The Netherlands.
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.
Volkov, A., King, H. F. & Coppens, P. (2006). J. Chem. Theory Comput. 2, 81-89.
Volkov, A., Koritsanszky T. S. & Coppens, P. (2004). Chem. Phys. Lett. 391, 170-175.
Zhang, Y., Mao, J., Godbout, N. & Oldfield, E. (2002). J. Am. Chem. Soc. 124, 13921-13930.