The term promolecule, originally used by Hirshfeld and Rzotkiewicz,
[27 ... “Electrostatic binding in the first-row AH and A_{2}
diatomic molecules”
Hirshfeld, F. L.; Rzotkiewicz, S. Mol. Phys. 1974,
27, 1319-1343.]
refers to a reference electron density consisting of overlapping,
spherically averaged, undeformed atoms prior to molecule formation.
It is equivalent to the independent atom model (IAM) in electron scattering,
and the spherical atom approximation in X-ray crystallography.
[…]
The properties of the promolecule are reasonably well understood, and
it is known to provide chemically sensible atomic radii,
[30 ... Gibbs, G. V.; Spackman, M. A.; Boisen, M. B. Am Mineral 1992, 77, 741.
31 ... Spackman, M. A.; Maslen, E. N. J Phys Chem 1986, 90, 2020.]
electrostatic binding energies,
[27 ... Hirshfeld, F. L.; Rzotkiewicz, S. Mol Phys 1974, 27, 1319.]
31 ... Spackman, M. A.; Maslen, E. N. J Phys Chem 1986, 90, 2020.
32 ... Trefry, M. G.; Maslen, E. N.; Spackman, M. A.
J Phys C: Solid State Phys 1987, 20, 19.]
atom–atom potentials,
[33 ... Spackman, M. A. J Chem Phys 1986, 85, 6579.]
diamagnetic shielding constants, and diamagnetic susceptibilities.
[34 ... Maksic, Z. B.; Eckert–Maksic, M.; Rupnik, K.
Croat Chem Acta 1984, 57, 1295.
35 ... Maksic, Z. B. J Mol Struct 1988, 170, 39.
36 ... Maksic, Z. B. In Molecules in Physics, Chemistry, and Biology;
Maruani, J., Ed.; Kluwer: Amsterdam, 1989, p. 49.]
Most recently, it is providing a novel means of partitioning space in a
crystal into molecules for visualization and other purposes.
[37 ... Spackman, M. A.; Byrom, P. G. Chem Phys Lett 1997, 267, 215.
38 ... McKinnon, J. J.; Mitchell, A. S.; Spackman, M. A.
Chem Eur J 1998, 4, 2136.
39 ... McKinnon, J. J.; Mitchell, A. S.; Spackman, M. A.
Chem Commun 1998, 2071.
xx ... “Molecular Surfaces from the Promolecule: A Comparison with
Hartree-Fock Ab Initio Electron Density Surfaces”
Mitchell, A. S.; Spackman, M. A. J. Comput. Chem. 2000,
21, 933-942.]
The use of high quality atomic wave functions provides an excellent description of the atomic densities, and our promolecule electron density surface is defined by an isosurface of the function ρ(r) = ∑_{A ∈ molecule} ρ_{A}(r) where ρ_{A}(r) is the spherically averaged atomic electron density due to atom A, and the total electron density ρ(r) is simply the sum of the individual atomic electron densities, each located at the appropriate nucleus.
A promolecule is defined to be a model of a molecule where the electron density distributions of each of its atoms have been spherically averaged and placed at their minimum energy positions (Hirshfeld and Rzotkiewicz, 1974).
The term promolecule, originally used by Hirshfeld and Rzotkiewicz, [Hirshfeld, F. L.; Rzotkiewicz, S. Mol Phys 1974, 27, 1319.] refers to a reference electron density model prior to molecule formation. It is equivalent to the IAM model in electron scattering, [Bonham, R. A,; Fink, M. High Energy Electron Scattering; Van Nostrand-Reinhold: New York, 1974.] and the spherical atom approximation in X-ray crystallography. [(a) Hirshfeld, F. L. Isr. J. Chem. 1977, 16, 87-229. (b) Becker, P., Ed. Electron and Magnetization Densities in Molecules and Crystals; Plenum: New York, 1980. (c) Coppens, P., Hall, M. B., Eds. Electron Distributions and the Chemical Bond; Plenum: New York, 1982.]
In 1970 Politzer and Harris proposed a definition of the charges associated
with the atoms in a molecule in terms of the molecular electronic density:
“if the total space of the molecule could be partitioned into regions
belonging to the individual atoms, then the electronic charge associated with
a given atom could be determined by integrating the electronic density over the
region of space belonging to that atom”[1Politzer, Peter; Harris, Roger R.
J. Am. Chem.
Soc. 1970, 92, 6451-6454.]
p_{i} = ∫ V_{i} ρ(r_{i}) d^{3}r_{i} | (1) |
where the integration is to be performed over the region of space V_{i} belonging to atom i. Based on this, the sum of all atom populations, p_{i}, in the molecule provides the total electronic charge (i.e. the number of electrons) of the molecule.
N_{atoms} ∑ i=1 p_{i} = ∫ V ρ(r) d^{3}r | (2) |
where the integration is now performed over all the region of space thet contains the molecule.
Politzer and Harris used this definition to calculate the atomic charges
in three linear molecules, i.e. acetylene, lithium acetylide and
fluoroethyne.[1Politzer, Peter; Harris,
Roger R.
J. Am. Chem. Soc. 1970, 92,
6451-6454.] Here we review their method for educational purposes.
The components of the position vector r in eq (1) are the
cartesian coordinates x, y and z, and the volume element
dV is given by the product of the differentials of the
cartesian coordinates, dV = dx dy dz
= d^{3}r.
The shape of linear molecules actually imposes the use of cylindrical
coordinates
x = R cosφ, y = R sinφ, z = z | (3) |
where z is assumed to be the molecular axis. The volume element
dV changes by the Jacobian determinant of the coordinate
change,[L2PAMoC manual. Online resource:
“Methods of Numerical Integration”
“§ 1.5 - Transformation of coordinates”, equation
(1.32)] i.e.
dV = dx dy dz = det
| (4) |
Then, the volume integral (2) can be written as a three-dimensional integral in terms of the cylindrical coordinates R, φ, and z:
∫ V ρ(r) d^{3}r = +∞ ∫ −∞ dz +∞ ∫ 0 R dR 2π ∫ 0 ρ(z,R,φ) dφ = +∞ ∫ −∞ g(z) dz | (5) |
where the radial probability density (also referred to as the radial
distribution), g(z), gives the probability of finding the
electron anywhere on the xy-plane perpendicular to the molecular axis at
z. In other words, it “is the density of electrons in this
plane”[2Brown, Richard E.; Shull,
Harrison
Int. J. Quantum Chem. 1968, 2,
663-685.] and can be designated “planar density”:[2Brown, Richard E.; Shull, Harrison
Int. J. Quantum Chem. 1968, 2, 663-685.]
g(z) = +∞ ∫ 0 R dR 2π ∫ 0 ρ(z,R,φ) dφ = +∞ ∫ 0 R f(z,R) dR | (6) |
where the function f(z,R) represents the angular integral
f(z,R) = 2π ∫ 0 ρ(z,R,φ) dφ | (7) |
The basic procedure[1Politzer, Peter;
Harris, Roger R.
J. Am. Chem. Soc. 1970, 92,
6451-6454.] consists in computing the “electron-count”
function for the linear molecule; the electron count is defined as[2Brown, Richard E.; Shull, Harrison
Int. J. Quantum Chem. 1968, 2, 663-685.]
G(z) = z ∫ −∞ g(z) dz | (8) |
It simply denotes the number of electrons on the negative side of the xy-plane intersecting the z-axis at the point z.
As the first step, we need the rules for numerical quadratures over the cylindrical coordinates R, φ, and z. We can choose the most appropriate ones, among those described in the PAMoc manual page on “Methods of Numerical Integration”. For the angular integral (7) over the interval [0, 2π], we choose the closed extended trapezoidal rule
f(z,R) = 2π ∫ 0 ρ(z,R,φ) dφ ≈ 2π n_{φ} n_{φ} ∑ i=1 ρ(z,R,φ_{i}) | (9) |
where φ_{i} = (i − 1) (2π/n_{φ}), ∀ i = 1, 2, …, n_{φ}. To solve the radial integral (6) we need a quadrature rule that spans the interval [0, +∞], like the generalized Gauss-Laguerre rule. However, the integrand R f(z,R) in eq (6) doesn't exactly correspond to the integrand of the quadrature formula, so let's rewrite it as follows
R f(z,R) = R e^{−R} e^{+R} f(z,R) = R e^{−R} F(z,R) | (10) |
where we have set F(z,R) = e^{R} f(z,R). The last expression in eq (10) has the proper form required for the integrand of the generalized Gauss-Laguerre quadrature rule for α = 1, then we obtain a solution of the radial integral (6):
g(z) = +∞ ∫ 0 R e^{−R} F(z,R) dR ≈ n_{R} ∑ i=1 w_{i} F(z,R_{i}) = n_{R} ∑ i=1 w_{i} e^{Ri} f(z,R_{i}) | (11) |
where R_{i} are the roots of the generalized Laguerre polynomial of order n_{R} and α = 1, with the associated weights w_{i}. Finally, the integral (8) is estimated by evaluating g(z) at a number n_{z} of points, z_{i}, evenly spaced by a small interval h and then applying the closed extended trapezoidal rule
G(z_{nz}) ≈ h n_{z} ∑ i=1 g(z_{i}) | (12) |
$contrl scftyp=rhf dfttyp=none runtyp=energy AIMPAC=.T. $end $system timlim=4 $end $guess guess=huckel $end $basis gbasis=N31 ngauss=6 ndfunc=1 $end $ELMOM IEMOM=3 IAMM=4 CUM=.TRUE. $END $ELFLDG IEFLD=2 $END $ELPOT IEPOT=1 $END $ELDENS IEDEN=1 $END $data CO at the experimental ground state geometry Cnv 4 C 6.0 0.0 0.0 0.0 O 8.0 0.0 0.0 1.128323 $end
The experimental bond length (r_{CO} = 1.1283 Å
or 2.1322 bohr)[3(a) Tiemann, E.;
Arnst, H.; Stieda, W. U.; Törring, T.; Hoeft, J.
Chem. Phys. 1982 67, 133-138.
(b) Frenking,G.; Loschen, C.; Krapp, A.; Fau, S.; Strauss, S. H.
J. Comput. Chem. 2007, 28, 117-126.
(c) Demaison, J.; Császár, A. G.
J. Mol. Struct. 2012, 1023, 7-14.]
is used for the calculation.
The z axis is chosen to be identical with the molecular axis, and
z increases from left to right when the molecule is written CO, with
z = 0 corresponding to the carbon atom.
Thus the electron count, G(z), is the number of electrons to
the left of a plane passing through the molecular axis at the point
z.
Thanks to the option AIMPAC=.T. specified in the $CONTROL group of the input, GAMESS provides an AIMPAC .wfn file that can be extracted from the .dat file at the end of the computation.
This wfn file can be opened in the interactive version of PAMOC to calculate the electron count function. The plot of G(z) vs. z is shown in Figure 1.
At this point, “it is necessary […] to decide upon a
reasonable and well-defined criterion in terms of which to partition
the space of the molecule. The criterion which is being proposed is as
follows. The atomic regions will be defined such that in the limiting
case of no interactions between the atoms, the electronic charge associated
with each one would be the same as for the free atom. This limiting case
can be represented by superposing free atom electronic densities, the atoms
being placed at the same relative positions as in the molecule.”[1Politzer, Peter; Harris, Roger R.
J.
Am. Chem. Soc. 1970, 92, 6451-6454.]
In the case of acetylene, lithium acetylide and fluoroethyne,
“the regions belonging to the individual atoms were determined
[…] in terms of the electron-count function for the superposed atoms.
The boundaries of the regions were established by means of hypothetical
planes which perpendicularly intersect the z axis at those points
for which G(z)_{promol} = 1, 7, and 13.
These points on the z axis then define the regions belonging to
the atoms in the molecule. The number of electrons associated with each
atom could subsequently be obtained from the molecular electron-count
function; since the number of protons is known for each atom, the net atomic
charge follows immediately.”[1Politzer,
Peter; Harris, Roger R.
J. Am. Chem. Soc. 1970, 92,
6451-6454.]
This recipe is applied to CO in Figure 2, which shows an enlarged detail
of Figure 1. The electron count function of the promolecule,
G(z)_{promol}, equals the number
of protons in the carbon nucleus (Z_{C} = 6) when the
z-coordinate has the value z_{boundary} = 1.1212 bohr.
The line perpendicular to the z-axis at z_{boundary}
intersects the electron count function of CO at
G(z_{boundary})_{CO} = 5.860 electrons, so
that the carbon is found to be positively charged by q(C) =
6 e − 5.860 e = + 0.140 e. This result is consistent with the
calculated Mulliken and Lowdin atomic populations:
TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS
ATOM MULL.POP. CHARGE LOW.POP. CHARGE
1 C 5.715758 0.284242 5.916101 0.083899
2 O 8.284242 -0.284242 8.083899 -0.083899
It is worth noting that the limiting case of non-interacting atoms used by Politzer and Harris as a reference is represented by the so-called promolecule density,[N1The promolecule.] even if the two authors have not made explicit use of this definition. In addition, the use of planes perpendicular to the internuclear axis to identify the boundaries of the atomic regions is an example of partitioning of the electron density in the real three-dimensional space into two unbounded Voronoi atomic cells.[N2Voronoi cell.] Again, the two authors do not give any explicit definition of their method of electron density partitioning, which however can be considered an extension of Voronoi's partitioning, since the boundary planes can be displaced from the bond midpoint (where, by definition, are located the faces of a Voronoi polyhedron).[N2Voronoi cell.]
Ultimately, the Politzer and Harris method consists of three steps:
q(C) = G(z_{boundary})_{promol} − G(z_{boundary})_{CO} = Z_{C} − G(z_{boundary})_{CO} = + 0.140 e | (13) |
The definition of atomic charge proposed by Politzer and Harris
“admittedly contains an element of arbitrariness in that the atomic
regions could be defined in terms of other criteria than the requirement
that the atoms have zero charges in the limit of no interaction.”.[1Politzer, Peter; Harris, Roger R.
J. Am. Chem. Soc. 1970, 92, 6451-6454.]
Indeed, by removing this constraint, imposed by the first step in the
above scheme, and allowing to choose for z_{boundary} any
appropriate value between the two bonded atoms, we can obtain partial charges
that better describe the properties of the molecule. As a consequence of this
choice, the two differences in eq (13) turn out to provide two different
models, so that in general they'll give different results.
In particular, the first difference defines the partial atomic charge
due to the deformation density, that is the density change
when going from the superposition of atomic densities to the final molecular
density.
q_{VDD} = −[ G(z_{boundary})_{mol} − G(z_{boundary})_{promol} ] = −[ z_{boundary} ∫ −∞ g(z)_{mol} dz − z_{boundary} ∫ −∞ g(z)_{promol} dz ] = − z_{boundary} ∫ −∞ g(z)_{DD} dz | (14) |
where g(z)_{DD} = g(z)_{mol} − g(z)_{promol} is the deformation density. Essentially, eq (14) is an example of application of the VDD method.
The last difference in eq (13) is the usual definition of partial atomic charge
z_{boundary}, bohr | G(z_{boundary})_{promol}, e | G(z_{boundary})_{mol}, e | q_{VDD}, e | q_{VP}, e |
---|---|---|---|---|
1.0661 (= ½r_{CO}) | 5.895 | 5.740 | 0.155 | 0.260 |
1.1212 (⇒ G_{promol} = 6) | 6.000 | 5.857 | 0.143 | 0.143 |
1.1258 | 6.009 | 5.867 | 0.142 | 0.133 |
1.1411 | 6.039 | 5.900 | 0.139 | 0.100 |
1.1863 (⇒ G_{mol} = 6) | 6.129 | 6.000 | 0.129 | 0.000 |
1.2422 | 6.245 | 6.129 | 0.116 | −0.129 |
Although in general the choice of the midpoint of the internuclear axis
would be rather unrealistic,[1Politzer, Peter; Harris, Roger R.
J. Am. Chem.
Soc. 1970, 92, 6451-6454.] and indeed
in the case of CO the criterion proposed by Politzer and Harris finds a
boundary plane placed 0.0551 bohr beyond the bond midpoint
(½ r_{CO} = 1.0661 bohr), and actually this is not
even enough, because it would be necessary to move the boundary plane at
z > bohr to find a partial negative charge on the carbon, as it was
observed experimentally.[4Muenter, J. S.
J. Mol. Spectr. 1975, 55, 490-491.]
On the other hand,
placing the boundary plane at z = ½ r_{CO}.
V_{i} = {x ∈ ℝ : d(x, P_{i}) ≤ d(x, P_{j}), j = 1, 2 , …, n}, | (N2.1) |
V_{i} = ⋂ j H_{ij}. | (N2.2) |