The total density, γ(r), of a charge distribution can be expressed as a sum over atomic contributions, γA(r - RA), each of which includes the nuclear charge ZA and the electron density ρA(r - RA):
| γ(r) = ∑ A γA(r - RA) = ∑ A [ZAδ(r - RA) - ρ(r - RA)] | (1) |
where δ(r - RA) is the Dirac delta function,
RA is the position vector of nucleus A, and the vector
rA = r - RA
marks the position of an infinitesimal density element with respect to the nuclear origin
RA.
Partitioning of the electron density into atomic contributions is widely used in molecular electronic structure calculations employing density functional theories (DFT), in order to avoid difficulties in numerical integration due to the presence of a cusp singularity at each nucleus.
[1Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.]
The functions γA(rA) and
ρA(rA),
which are associated to a nuclear center at RA, are called the
atoms-in-molecules (AIM) charge density (CD) and electron density (ED) of atom A, respectively.
The integral of
ρA(rA)
over all space gives the atomic population pA, i.e. the zero-th rank moment
of atom A,
| pA = ∫ρA(rA) d3rA. | (2) |
Higher rank electrostatic moments are obtained by evaluation of single-center three-dimensional integrals of the type
| mA(n) = ∫ ρA(rA) rAn d3rA | (3) |
where mA(n) is the n-th rank unabridged cartesian multipole tensor of atom A, and the expression rAn denotes the direct or polyadic product of n factors rA. Of course, pA = mA(0). The set of multipole moments mA(n) ∀ A ∈ {molecule or unit cell} is called a distributed multipole analysis (DMA) of the ED. The DMA of the CD can be derived in a similar way, yielding
| qA = ZA - pA and MA(n) = - mA(n) for n > 0. |
The quantity qA is called atomic charge.
Defining the atomic contributions to the total ED (and, in turn, to the total CD) is not a trivial task. It can be accomplished by using a set of nuclear weight functions, wA(r):
| ρA(rA) = wA(r)ρ(r) | (4) |
The weight functions, wA(r), measure to what degree the given point of space can be considered to belong to atom A. They are required to fulfill the conditions
| 0 ≥ wA(r) ≤ 1 and ∑A wA(r) = 1 ∀ r ∈ ℝ3, | (5) |
which illustrate the concept of atomic decomposition of the identity in 3D space.[2"Atomic Decomposition of Identity:
General Formalism for Population Analysis and Energy Decomposition"
Mayer I.; Hamza A.
Int. J. Quantum Chem. 2005, 103, 798-807.]
Here is the list of the partitioning schemes which are available in PAMoC for experimental and/or theoretical densities:
Space can be partitioned with well-defined discrete boundaries or with overlapping functions, often referred to as fuzzy boundaries. Among the above partitioning schemes, only QTAIM provides non-everlapping atoms with discrete boundaries, whereas the remaining ones give fuzzy atoms with boundaries extending to infinite.
Stewart, Becke, Hirshfeld and QTAIM atoms are defined in real 3D space. On the other hand, Mulliken and Stone atoms are defined in Hilbert space and, in addition, they are restricted to theoretical wavefunctions given in analytical form (.wfn files).
Partitioning schemes in the Hilbert space are known to be unstable with respect to changes of the basis set. Indeed, if the basis set includes very diffuse functions, overlap densities involving them extend over many atoms. Therefore, partitioning of the electron density in the real 3D space should be preferred whenever it is possible.
However, DMA's in Hilbert space often give a satisfactory representation of the molecular charge distribution. In addition, they can be calculated very fast because they use an exact and very efficient Gauss-Hermite quadrature. For these reasons, they are provided by PAMoC as a default procedure for all wave functions given in analytical form (.wfn files).
"Debye suggested that the one-electron density function for the molecule or unit cell at
the equilibrium nuclear configuration be represented by a superposition of isolated
atoms [3
"Röntgeninterferenzen und Atomgrösse"
Debye, P.
Physikalische Zeitschrift 1930, 31, 419-428.]
that rigidly follow the nuclei on which they are centered. ... An extension of Debye's
suggestion is to partition the molecular charge density into pseudoatom charge densities.
Each pseudoatom density is defined by a finite surface harmonic (multipole) expansion about
its nucleus and the coefficient radial functions are determined by a least squares fit to the
given molecular density" (Stewart, R.F. [4
"Generalized X-Ray scattering factors for pseudoatoms in molecules and crystals"
Stewart, R. F.
Acta Cryst. A (Supplement) 1975, 31, s218.,
5
Stewart, R. F.
Acta Cryst. 1976, A32, 564-574.]):
| ρA(rA) = fA(rA)†pA | (6) |
The elements of vector pA are the multipole populations of the pseudoatom A and the elements of vector fA(rA) are nuclear-centered solid harmonics functions which account explicitly for the dependence of the ED on the spacial coordinates rA = |xA, yA, zA|†. The nuclear weight functions, wA(r), are implicitly assumed to be equal 1 everywhere if the multipole expansion is centered on atom A and zero otherwise.
In Fourier space the procedure is equivalent to approximating the molecular
form factor in terms of generalized X-ray scattering factors for the
pseudoatoms.[6
(a) Stewart, R. F.
J. Chem. Phys. 1969, 51, 4569-4577;
(b) Stewart, R. F.; Bentley, J.; Goodman, B.
J. Chem. Phys. 1975, 63, 3786-3793;
(c) Hansen, N. K.; Coppens, P.
Acta Cryst. 1978, A34, 909-921;
(d) Coppens, P.
X-ray Charge Densities and Chemical Bonding;
Oxford University Press: New York, 1997.
]
Atomic populations pA are the main outcome of
least-squares multipole analysis programs, like
VALRAY[7
Stewart, R. F.; Spackman, M. A.; Flensburg, C.
VALRAY User's Manual (Version 2.1),Carnegie-Mellon University, Pittsburg,
and University of Copenhagen, Copenhagen, 2000.]
and XD.[8
Koritsanszky, T.; Mallison, P.; Macchi, P.; Volkov, A.; Gatti, C.; Richter, T.;
Farrugia, L.:
XD2006. A Computer Program Package for Multipole Refinement, Topological
Analysis of Charge Densities and Evaluation of Intermolecular Enerigies from
Experimental or Theoretical Structure Factors (2007).]
PAMoC is interfaced to both VALRAY and XD and, as a consequence, it can use their multipole analysis
of experimental ED's to evaluate atomic and molecular properties, including, for instance, inner and
outer moments and electrostatic interaction energies. The variance-covariance matrix of the best values of the least-squares populations is used to estimate standard uncertainties on all the derived quantities, as a measure of the reliability of the pseudoatom model.
Though Gill has shown that Stewart atoms can be extracted directly from theoretical
ED's,[9
"Extraction of Stewart Atoms from Electron Densities"
Gill, P. M. W.
J. Phys. Chem. 1996, 100, 15421-15427.]
most applications of pseudoatom partitioning of theoretical ED's are based on the least-squares
multipole refinement of theoretical structure factors.[6b
Stewart, R. F.; Bentley, J.; Goodman, B.
J. Chem. Phys. 1975, 63, 3786-3793.,
10
(a) Spackman, M. A.
J. Chem. Phys. 1986, 85, 6587-6601.
(b) Spackman, M. A.; Byrom, P. G.; Alfredsson, M.; Hermansson, K.
Acta Cryst. 1999, A55, 30-47.
(c) Koritsanszky, T.; Volkov, A.; Coppens, P.
Acta Cryst. 2002, A58, 464-472.
(d) Volkov, A.; Li, X.; Koritsanszky, T.; Coppens, P.
J. Chem. Phys. A 2004, 108, 4283-4300.]
PAMoC is able to perform least-squares multipole refinement of theoretical ED's in 3D space.
This is the default procedure, whenever the IDF is a 3D grid file. Otherwise, when the IDF is an
AIMPAC wavefunction file, the keyword lsq is needed, in order to instruct PAMoC to build a 3D
grid of ED values and then to perform a least-squares refinement.
Least-squares refinement of experimental and theoretical ED's provide a set of pseudoatom populations in the form of spherical harmonics tensors, that can be easily converted into traceless cartesian moments by a linear transformation. However a volume integration is needed to evaluated unabridged cartesian moments. This can be achieved by means of keyword stw.
Becke's weight functions[1Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.] are defined as
| wA(r) = VA(r)/∑BVB(r), | (7) |
where
| VA(r) = ∏B≠ask(μAB) | (8) |
is the Voronoi polyhedron function constructed from the smoothing function
| sk(μAB) = ½[1 - fk(μAB)]. | (9) |
The polynomial fk is obtained by the iteration formula
| fk = ½fk-1(3 - f2k-1), | (10) |
starting from f0 = μAB and iterating k-times. As suggested by Becke, the value assigned by PAMoC to the stiffness parameter k is 3, which provides a soft transition between atoms. The dependence on spacial coordinates is described by the confocal elliptical coordinate
| μAB = (rA - rB)/RAB, | (11) |
where rA, rB and RAB denote distance to
nucleus A, distance to nucleus B, and the internuclear separation, respectively.
Becke further proposed an atomic size adjustment[1Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.]
for heteronuclear molecules by adding to the argument μAB of the smoothing
function sk the correcting term
| νAB = ¼(1 - μAB2)(1 - χAB2)/χAB | (12) |
where χAB denotes the ratio of the Bragg-Slater covalent radii[11Slater, J. C.
J. Chem. Phys. 1964, 41, 3199-3204.] of atoms A and B.
For hydrogen, however, Becke suggested a radius of 0.35 Å, somewhat larger than Slater's value
of 0.25 Å. Stone argued that using the Bragg-Slater radii in the Becke partitioning
"leads to rather implausible multipole moments (rather large and not in keeping with chemical
intuition) though they lead to accurate electrostatic potentials".[12Stone, A. J.
J. Chem. Theory Comput. 2005, 1, 1128-1132.]
For this reason, the default
procedure in Stone's GDMA program is to set all the atom radii equal
to a value of 0.65 Å, except for hydrogen (to which a value of 0.325 Å is assigned).
According to user's preferences, PAMoC can either use Becke's recipe for atomic size
adjustments, or ignore them so that the boundary between neighboring atoms comes halfway between
them.
Though Becke's method has been developed for handling the integrals of density functional theory,
it has been efficiently applied to the distribute multipole analysis of theoretical
electron densities,[12Stone, A. J.
J. Chem. Theory Comput. 2005, 1, 1128-1132.]
to estimation of overlap populations, bond orders and valences of fuzzy atoms[13Mayer, I.; Salvador, P.
Chem. Phys. Lett. 2004, 383, 368-375.]
and to energy partitioning in one- and two-atomic contributions.[14(a) Salvador, P.; Mayer, I.
J. Chem. Phys. 2004, 120, 5046-5052.
(b) Salvador, P.; Mayer, I.
J. Chem. Phys. 2007, 126, 234113-10.
(c) Mayer, I.
Phys. Chem. Chem. Phys. 2006, 8, 4630-4646.]
It can be applied to experimental electron densities as well.
Hirshfeld defined the nuclear weight function wA(r) according to the
common sense stockholder principle,[ (a) Hirshfeld, F.L. Theoret. Chim. Acta 1977, 44, 129-138. (b) Hirshfeld, F.L. Isr. J. Chem. 1977, 16, 226-229]
wA(r) = ρ°A(r - RA)/ρ°(r),
by which each atom A partecipates in the molecular profit ρ(r) in proportion
to its share wA(r) in the promolecular investment ρ°(r).
Here, ρ°A(r - RA) is the electron density of the
isolated spherical atom A (proatom), placed at its position in the molecule, and
ρ°(r) = ∑A ρ°
It has been shown that Hirshfeld atoms exhibit the minimum information distance to their free atom reference[Nalewajski, R. F.; Par, R. G. PNAS 2000, 97, 8879-8882] and satisfy the requirement of maximum transferability,[Ayers, P. W. J. Chem. Phys. 2000, 113, 10886-10898] i.e. they retain certain characteristic properties in a wide variety of environments. Hirshfeld's stockholder partitioning scheme has become to prominence in recent years with the widespread use of Hirshfeld surface analysis and fingerprint plots for visualizing and analyzing intermolecular interactions in crystals.[Spackman, M. A.; Byrom, P. G. Chem. Phys. Lett. 1997, 267, 215-220. Mc Kinnon, J. J.; Mitchell, A. S.; Spackman, M. A. Chem. Eur. J. 1998, 4, 2136-2141. Mc Kinnon, J. J.; Spackman, M. A.; Mitchell, A. S. Acta Crystallogr. 2004, B60, 627-668. Spackman, M. A.; Mc Kinnon, J. J. CrystEngComm 2002, 4, 378-392. Mc Kinnon, J. J.; Fabbiani, F. P. A.; Spackman, M. A. Crys. Growth Des. 2007, 7, 755-769. Spackman, M. A.; Jayatilaka, D. CrystEngComm 2009, 11, 19-32.]
One of the best known exhaustive partitions of 3D space is provided by the quantum theory of atoms in molecules (QTAIM) of Bader and coworkers,[Bader, R.F.W. Atoms in Molecules: A Quantum Theory; Clarendon Press: Oxford, 1990. Bader, R.F.W. Chem. Rev. 1991, 91, 893-928. (a) Popelier, P. L. A. Atoms in Molecules. An Introduction; Pearson: London, Great Britain, 2000. (b) Popelier, P. L. A. In Structure and Bonding: Intermolecular Forces and Clusters; Wales, D. J., Ed.; Springer: Heidelberg, Germany, 2005; Vol. 115; p. 1. (c) Popelier, P. L. A.; Aicken, F. M.; O'Brien, S. E. In Chemical Modelling: Applications and Theory; Royal Society of Chemistry Specialist Periodical Report, Vol. 1; Hinchliffe, A., Ed.; The Royal Society of Chemistry: Cambridge, Great Britain, 2000; Chap. 3; p. 143. (d) Popelier, P. L. A.; Smith, P. J. In Chemical Modelling: Applications and Theory; Royal Society of Chemistry Specialist Periodical Report, Vol. 2; Hinchliffe, A., Ed.; The Royal Society of Chemistry: Cambridge, Great Britain, 2002; Chap. 8; pp. 391-448.] which introduces a new concept of atom, based on the quantum chemical topology (QCT), of the electron density.[Popelier, P. L. A.; Aicken, F. M. ChemPhysChem 2003, 4, 824-829. (a) Popelier, P. L. A.; Rafat, M. Chem. Phys. Lett. 2003, 376, 148-153. (b) Rafat, M.; Popelier, P. L. A. J. Chem. Phys. 2005, 123, 204103-7. (c) Rafat, M.; Popelier, P. L. A. J. Chem. Phys. 2006, 124, 144102-7. (d) Rafat, M.; Popelier, P. L. A. J. Comput. Chem. 2007, 28, 832-838.]
According to the QTAIM, the boundary condition for a quantum subsystem leads to a unique, exhaustive and disjoint partitioning of the real space of a molecule or crystal into a set of mononuclear regions. QTAIM performs a topological analysis of ρ(r) and divides all space into non overlapping regions, the atomic basins Ω, separated by interatomic surfaces S(r) that satisfy the boundary condition of zero flux of the gradient vector field of the electron charge density, i.e. ∇ρ(r)˙n(r) = 0 ∀ r ∈ S(r), where n(r) is a unit vector normal to the surface at r. Here the nuclear weight function wA(r) equals 1 for each r ∈ ΩA and zero otherwise.
QTAIM partitioning can be applied either to the molecule in the crystal form, where the atomic basins are defined by the interatomic surfaces (IAS) in the crystal (QCT atoms), or to the molecule extracted from the crystal, whose atoms are separated by interatomic surfaces inside the molecule, but have boundaries artificially moved to infinity outside the molecule. The latter partitioning (hereafter referred as QCT*) gives molecules whose electron densities overlap with one another, so that their electrostatic energies need to be corrected for the penetration energy as in the case of any other type of fuzzy boundary molecules.
Within the LCAO MO theory, first order density functions can be expressed in a very general manner using the diagonal form:[(a) Carbó-Dorca, R.; Bultinck, P. "Quantum mechanical basis for Mulliken population analysis" J. Math. Chem. 2004, 36, 231-239. (b) Carbó-Dorca, R.; Bultinck, P. "A general procedure to obtain quantum mechanical charge and bond order molecular parameters" J. Math. Chem. 2004, 36, 201-210. , (a) Carbó-Dorca, R. "Mathematical aspects of the LCAO MO first order density function (1): atomic partition, metric structure and practical applications" J. Math. Chem. 2008, 43, 1076-1101. (b) Carbó-Dorca, R. "Mathematical aspects of the LCAO MO first order density function (2): Relationships between density functions" J. Math. Chem. 2008, 43, 1102-1118.]
| ρ(r) = φ†(r) Ω φ(r) = ∑i ωi φi*(r) φi(r) = ∑i ωi ρi(r), | (13+0) |
where the (column) vector φ(r) = {φi(r)} contains the set of molecular orbital (MO) functions and the diagonal matrix Ω = {ωi} contains the MO occupation numbers. In addition, ρi(r) = φi*(r) φi(r) defines the first-order density function of the ith MO. The first order density function has a norm equal to (i.e. integrates to) the total number of electrons, N:
| <ρ> = ∫ ρ(r)d3r = N. | (13+1) |
Within the LCAO (Linear Combination of Atomic Orbitals) formalism, the MO set is expressed as a linear combination of one-electron basis functions {χμ(r)}, generally centered on the various nuclei of the system, i.e.:
| φi(r) = ∑μcμiχμ(r) = χ(r) ci, | (13+2) |
where the column vector ci = {cμi} contains the coordinates of the ith MO with respect to the basis set χμ(r) and the row vector χ(r) = {χμ(r)} contains the AO basis functions. The LCAO expansion leads to the well-known expression of the first order electron density function ρ(r) as the sum of products of two basis functions, each one centered on a given nuclear site, with coefficients determined from the density matrix:
| ρ(r) = ∑i ωi [ ∑μ c*μi χ*μ(r) ] [ ∑ν cνi χν(r) ] = ∑μν ( ∑i ωi c*μi cνi ) χ*μ(r) χν(r) = ∑μν Dμν χ*μ(r) χν(r) = ∑μν Dμν |χμ><χν| | (13+3) |
where the (symmetric) density matrix D, also (and preferably[Carbó-Dorca, R. "LCAO MO first order density functions: Partition in monocentric and bicentric terms, reciprocal MO spaces, invariant transformations and Euclidian atomic populations" J. Mol. Struct.: THEOCHEM 2010, 943, 32-41.]) called charge and bond order matrix, is defined as:
| D = {Dμν} = ∑i ωi ci ci† = CΩC† | (13+4) |
where the diagonal matrix Ω = {ωi} contains the MO occupation numbers and the rectangular matrix C = {cμi} contains the MO coefficients as columns. The norm of the first order density function can be written as:
| N = <ρ> = ∑μν Dμν ∫ χ*μ(r) χν(r) d3r = ∑μν Dμν <χμ|χν> = ∑μν Dμν Sμν = Tr(DS) = ∑μν Pμν | (13+5) |
where the symmetric square matrix:
| S = {Sμν} = ∫ χ(r) χ†(r) d3r | (13+6) |
is a positive definite metric or overlap matrix for the AO basis set. We have also introduced the (symmetric) population matrix P with elements Pμν = Dμν Sμν as the inward product (also known as the Hadamard product or the Schur product or the entrywise product) of the two matrices D and S, i.e. P = D∘S = {Dμν Sμν}.
In any molecular case, the LCAO MO density function can be expressed as the sum of two functions:[Carbó-Dorca, R. "LCAO MO first order density functions: Partition in monocentric and bicentric terms, reciprocal MO spaces, invariant transformations and Euclidian atomic populations" J. Mol. Struct.: THEOCHEM 2010, 943, 32-41.]
| ρ(r) = ρAA(r) + ρAB(r) | (13+7) |
containing monocentric (A) and bicentric (B) elements, namely:
| ρAA(r) = ∑A ∑μ∊A ∑ν∊A DμνAA |χμA><χνA| | (13+8) |
and
| ρAB(r) = 2 ∑A ∑B>A ∑μ∊A ∑ν∊B DμνAB |χμA><χνB| | (13+9) |
The norm of the density function can be written in terms of the partition above:
| <ρ(r)> = <ρAA(r)> + <ρAB(r)> = N | (13+10) |
and, in any circumstance, it is equal to the number of electrons N. It can be proved that the norm <ρAA(r)> is positive definite in any case, but can be lesser or greater than the number of electrons, while the pseudonorm <ρAB(r)> is a real number, yielding even negative values in circumstances where <ρAA(r)> becomes greater than the number of electrons.
The sum over the μth row of the AO population matrix P leads to the definition of a vector GOP whose elements are the so-called gross orbital populations:
| GOPμ = ∑ν Pμν | (13+11) |
Finally, summing the populations of all the atomic orbitals (basis functions) centered on a particular atom A leads to the gross atom population GAPA on that atom (hereafter referred to as Mulliken population, pA):
| pA = GAPA = ∑μ∊A GOPμA = ∑μ∊A ∑ν Dμν ∫ χ*μ(r) χν(r) d3r = ∑μ∊A∑B ∑ν∊B DμνAB SμνAB = ∑μ∊A∑B ∑ν∊B PμνAB, | (13+12) |
where the superscripts on vector and matrix elements stress the centers of the involved basis functions.
Mulliken atomic population can be considered as the summation of AO population contributions pμA, i.e.
| pA = ∑μ∊A pμA | (13+13) |
where
| pμA = ∑B ∑ν∊B PμνAB = ∑ν∊A PμνAA + ∑B≠A ∑ν∊B PμνAB. | (13+14) |
This procedure, known as Mulliken Populaton Analysis (MPA),[(a) Mulliken, R. S. J. Chem. Phys. 1935, 3, 573-585. (b) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833-1840.] can be generalized to the evaluation of higher multipoles of atom A, yielding a nuclear-centered Distributed Multipole Analysis (DMA):
| mA(n) = ∑μ∊A ∑ν Dμν ∫ rAn χ*μ(r) χν(r) d3r = ∑μ∊A ( ∑ν∊A DμνAA <χμA|rAn|χνA> + ∑B≠A ∑ν∊B DμνAB <χμA|rAn|χνB> ). | (13+15) |
In the Mulliken partitioning scheme half of the multipole moments associated with the overlap densities of the bicentric term are allocated to each of the sites A and B from which the basis functions came.
Two types of primitive AO functions are considered here, namely:
a cartesian Gaussian type orbital (GTO)
| χμA(r) = NG (x - Ax)l (y - Ay)m (z - Az)n exp[-ζ(r - A)2] ≡ NG xAl yAm zAn exp(-ζ rA2) | (13+16) |
and a cartesian Slater type orbital (STO)
| χμA(r) = NS (x - Ax)l (y - Ay)m (z - Az)n exp[-ζ|r - A|] ≡ NS xAl yAm zAn exp(-ζ rA) | (13+17) |
where rA = r - A is the electron position relative to the position of atom A on which the AO is centered, ζ is the exponent and NG and NS are normalization constants. When l = m = n = 0, we have an s-type AO, when l + m + n = 1 we have a p-type AO, and so on. The multipole integrals between two GTO's, χμA(r) and χνB(r), are evaluated through a Gauss-Hermite quadrature.
As a default procedure, PAMoC provides a Mulliken DMA for all wave functions given in analytical form (.wfn files). The multipoles are given in the form of unabridged Cartesian tensors.
At variance from the Mulliken partitioning scheme, the Stone nearest-site allocation algorithm involves transferring all of the multipoles to the nearest site. Let's consider the product of two GTO's, one centered at A and the other centered at B:
| χμA(r) = N1 xAn1x yAn1y zAn1z exp(-ζ1 rA2) and χνB(r) = N2 xBn2x yBn2y zBn2z exp(-ζ2 rB2) | (14+0) |
where rA = r - A and rB = r - B are the electron positions relative to the positions of atoms A and B on which the GTO's are centered, ζ1 and ζ2 are the orbital exponents, n1 = (n1x, n1y, n1z) and n2 = (n2x, n2y, n2z) are the Cartesian angular momenta, and N1 and N2 are normalization constants. The exponential term represents the radial part of the GTO, whereas the polynomial represents the angular part, in that the sum of the Cartesian angular momenta nx + ny + nz = 0, 1, 2, 3, ... corresponds to orbital functions of type s, p, d, f, ... We have:
| χμA(r) χνB(r) = N1 N2 xAn1x yAn1y zAn1z xBn2x yBn2y zBn2z exp[-(ζ1 rA2 + ζ2 rB2)] | (14+1) |
where the exponential term satisfies the equality:
| ζ1 rA2 + ζ2 rB2 = [ζ1 ζ2 / (ζ1 + ζ2)] (A − B)2 + (ζ1 + ζ2)rP2. | (14+2) |
It constitutes the Gaussian Product Theorem for Two 1s-type GTO's,["Electronic Wave Functions. I. A General Method of Calculation for the Stationary States of Any Molecular System" Boys, S. F. Proc. Roy. Soc. A 1950, 200, 542-554.] which in other words states that the product of two 1s GTO's, one centered at A and the other centered at B, is another Gaussian centered at P = (ζ1A + ζ2B)/(ζ1 + ζ2) with exponent ζ3 = ζ1 + ζ2:
| exp(-ζ1 rA2) exp(-ζ2 rB2) = exp[- (ζ1ζ2/ζ3) (A - B)2] exp(- ζ3 rP2) ≡ KAB exp(- ζ3 rP2), | (14+3) |
where KAB = exp[- (ζ1ζ2/ζ3) (A - B)2] and rP = r - P. Since 1s GTO's are separable along the three Cartesian coordinates:
| exp(- ζ3 rP2) = exp(- ζ3 xP2) × exp(- ζ3 yP2) × exp(- ζ3 zP2) = ∏ s=x,y,z exp(−ζ3 sP2) , | (14+4) |
the product of any two GTO's can be written as equation:
| χμA(r) χνB(r) = N1 N2 KAB ∏ s=x,y,z [ sAn1s sBn2s exp(−ζ3 sP2) ] | (14+5) |
The products of the cartesian angular parts, sAn1s sBn2s for s = x, y, z, can also be expressed conveniently in terms of the components xP, yP, and zP of vector rP. For instance:
|
sAn1s = (s − As)n1s = (s − Ps + Ps − As)n1s = (sP + Ps − As)n1s =
n1s
∑
i=0
(
n1s i ) (Ps − As)n1s−i sPi | (14+6) |
where a standard binomial espansion has been used. Likewise,
|
sBn2s = (s − Bs)n2s = (sP + Ps − Bs)n2s =
n2s
∑
j=0
(
n2s j ) (Ps − Bs)n2s−j sPj. | (14+7) |
Using these, the product sAn1s sBn2s can be written as a summation of sP to various powers:
| sAn1s sBn2s = n1s ∑ i=0 n2s ∑ j=0 fij(n1s, n2s, Ps−As, Ps−Bs) sPi+j | (14+8) |
where fij(n1s, n2s, Ps−As, Ps−Bs) is the polynomial:
|
fij(n1s, n2s, Ps−As, Ps−Bs) =
(
n1s i ) ( n2s j ) (Ps − As)n1s−i (Ps − Bs)n2s−j | (14+9) |
Whence the Gaussian Product Theorem for Two GTO's of any type can be written as equation
| χμA(r) χνB(r) = N1 N2 KAB ∏ s=x,y,z [ n1s ∑ i=0 n2s ∑ j=0 fij(n1s, n2s, Ps−As, Ps−Bs) sPi+j exp(−ζ3 sP2) ] | (14+10) |
It is worth noting that the Gaussian Product Theorem holds for the product of any number of GTO's on different centers.["The general Gaussian product theorem" Besalú, E.; Carbó-Dorca, R. J. Math. Chem. 2011, 49, 1769-1784.]
Now we are able to evaluate the multipole expansion of the AO overlap density ρμν(r) = χμA(r) χνB(r) around the atomic center P = (ζ1A + ζ2B)/(ζ1 + ζ2):
|
<xPnx yPny zPnz>AB =
+∞
∫
−∞
xPnx yPny zPnz
χ*μA(r) χνB(r) d3r =
<χμA| xPnx yPny zPnz |χνB> = N1 N2 KAB ∏ s=x,y,z [ n1s ∑ i=0 n2s ∑ j=0 fij(n1s, n2s, Ps−As, Ps−Bs) +∞ ∫ −∞ sPi+j+ns exp(−ζ3 sP2) dsP ]. | (14+12) |
Odd values of i+j+ns result in odd functions whose integrals vanish. For even values of i+j+ns, a solution for the integrals is given by
| +∞ ∫ −∞ sPi+j+ns exp(−ζ3 sP2) dsP = ( π ζ3 )½ (i+j+ns−1)!! (2 ζ3)(i+j+ns)/2 , | (14+13) |
and, in those cases we have the final expression to evaluate the multipole moments of the AO overlap density of two Gaussians with arbitrary Cartesian angular factors:
| <xPnx yPny zPnz>AB = N1 N2 KAB ( π ζ3 )3⁄2 ∏ s=x,y,z [ n1s ∑ i=0 n2s ∑ j=0 fij(n1s, n2s, Ps−As, Ps−Bs) (i+j+ns−1)!! (2 ζ3)(i+j+ns)/2 ]. | (14+14) |
In the case of monocentric integrals (B = A and P = A), Eq. (14+14) simplifies to:
| <xPnx yPny zPnz>AA = N1 N2 ( π ζ3 )3⁄2 ∏ s=x,y,z (ns+n1s+n2s−1)!! (2 ζ3)(ns+n1s+n2s)/2 | (14+15) |
Of course, we must keep in mind that for the summations in Eqs (14+14) and (14+15) only terms of even values of i+j+ns and ns+n1s+n2s survive. For nx = ny = nz = 0, Eqs (14+14) and (14+15) give the overlap integrals of two Gaussian primitives.
As a default procedure, PAMoC provides a Stone DMA for all wave functions given in analytical form (.wfn files). The multipoles are given in the form of unabridged Cartesian tensors.
Stone DMA's in the form of unabridged Cartesian tensors are also calculated by GAMESS-US through the keyword "STONE". On the other hand, DMA's in the form of wave function normalized spherical harmonic multipoles can be obtained through Stone's GDMA computer program in conjunction with the formatted checkpoint file produced by Gaussian.