Physics Medicine
Chemistry,1970-2014-Chemistry-2003.jsonl.gz.gz
Chemistry
38ba1334f3bce6b18b47303feabf6cd68266c5c9
Rules for Minimal Atomic Multipole Expansion of Molecular Fields

A non-empirical minimal atomic multipole expansion (MAME) defines atomic charges or higher multipoles that reproduce electrostatic potential outside molecules. MAME eliminates problems associated with redundancy and with statistical sampling, and produces atomic multipoles in line with chemical intuition.

The problem of representing the electrostatic potential outside a molecule using atomic charges or higher atomic multipoles is very important for understanding intermolecular forces. Atomic partial charges, an important part of chemical intuition, are defined in many different ways for different purposes. Chemically-derived (CD) charges, such as Mulliken 1 or Löwdin 2 , often describe molecular fields poorly. 3 More recent schemes partition molecular density into atomic regions, which may or may not overlap. 4 Similar approaches have been developed for solids. 5 Most attractive for our purposes are potentialderived (PD) charges, which avoid representation of the density by producing the 'best' fit to the molecular potential directly. 6,7 Atomic dipoles and quadrupoles 7 are often used to increase accuracy in solvation problems 8 and force field calculations. 9 Induced atomic dipoles appear naturally in electronic polarization of molecular solids 10 to account for the small part of molecular polarization that is due to the deformation of atomic orbitals and is not captured by redistribution of charges.
Computational schemes for PD multipoles such as Merz-Kollman (MK), 11 CHelp 12 , or CHelpG 13 differ mainly in the sampling domain and the resulting atomic charges are strongly method-dependent. 14 Worse still, PD methods often yield atomic charges that are counterintuitive, such as negative charges on hydrogens in alkanes. 15 Higher multipoles only increase the redundancies inherent in distributed multipole analysis, improving on the accuracy of the field at the expense of instability in the multipole values. The severity of the problem can be somewhat reduced with SVD techniques, 16,14 or by introducing restraints. 17 Our approach does not use sampling and eliminates redundancies before they appear.
We approximate the true molecular potential, φ(r), as a sum of multipoles of strength q k centered at nuclear positions r i , where φ k (r) is the potential due to the kth multipole of unit strength: φ k (r) = 1/r for charge, (n · r)/r 3 for a dipole in the direction n, and so on. Since ∇ 2 φ approx = 0 everywhere except at r i , but ∇ 2 φ = 4πρ(r), the atomic multipole expansion can only be accurate in regions where ρ(r) ≈ 0. Furthermore, φ on any closed surface S on which ρ = 0, determines φ(r) everywhere outside S. We therefore choose S to be an isodensity surface, ρ(r) = f , where f is sufficiently small to ensure negligible charge beyond S, but with sufficient potential on S for a determining fit (Fig. 1). We chose φ approx to minimize over S which leads to a system of linear equations Atomic multipoles defined in this way are fully rotationally invariant, which is an issue with some PD schemes. 13 The error σ can be compared to φ, The crucial issue remaining is the choice of a set of multipoles. We choose a minimal set, usually one scalar value per atom, and add additional multipoles to describe lone pairs when necessary, based on the Lewis structure. This carefully chosen minimal atomic multipole expansion (MAME) set avoids redundancies but is within ∼ 1 mH everywhere beyond S.
We illustrate MAME with three molecules: n-pentane, which is a classic example of difficulties encountered in PD schemes; glycine (standard and zwitterion), as a typical application in biochemistry; and water, to see how general MAME rules apply to a small polar molecule. All densities and potentials are produced on a cubic mesh by the Gaussian 98 program, 18 at the B3LYP/aug-cc-pVTZ level (6-311++G** for pentane). Surface integrals (3) are computed by triangulation of S. The program runs within a few seconds, and is available on request. Figure 1 shows φ(r) on S for n-pentane. Red spots (φ > 0) show an excess of positive charge near each hydrogen, but all PD schemes tested yield some or all hydrogens negative. Closer inspection of Fig. 1 reveals that the positive regions occupy less solid angle around hydrogens than would be produced by a positive charge. Such a potential is consistent with a dipole with a negative charge pointing inwards.
Our first rule is therefore to assign a charge to all nuclei but protons, to which we assign a dipole moment instead. The hydrogen atom is special as its sole electron participates in the bond, leaving no electron density centered on the proton. This unique property of hydrogens is wellknown in X-ray structure analysis, which systematically underestimates the C-H bond lengths for this reason. Table 1 Partial atomic charges in n-pentane. f = 5×10 −4 au., error σ as in (2), (%)= (σ/φ), eφ = 3.6 mHartree (1mH = 27meV ∼ kT at 300K). "µ" indicates atomic dipoles, "µr" -dipoles restricted along H-C bonds.
Mulliken charges are intuitively meaningful but produce large errors in the potential (Table 1). PD charges are negative on some hydrogens and still give significant errors. Adding dipoles reduces the potential error, but at the cost of producing meaningless multipoles. 15 Our scheme with charges on all atoms produces similar (though better) results, but we do far better (line 2) when the charges on hydrogens are replaced with dipoles. All dipoles come out similar in magnitude (numbers in brackets, in au.) and point toward C within 20 • of the H-C bond. The hydrogen dipoles can be safely restricted to lie along the H-C bonds (last line) with the accuracy still better than that of charges. All multipoles have reasonable values, including small charges on carbons. Note that we have now described the field outside the molecule more accurately than any existing scheme, with only one parameter per nuclues (a charge on each carbon and a bond-directed dipole on each hydrogen).
The same choice of multipoles yields a 1.05 mH error (=2%) in the glycine zwitterion, (NH 3 ) + -CH 2 -COO − , down from 4% with charges alone and 4%-6% with standard PD schemes. The glycine zwitterion is highly polar with dipole µ = 10.3 D, which MAME recovers within 0.1% accuracy.   Table 2 lists MAME results for glycine in its standard form, NH 2 -CH 2 -COOH, and illustrates the need for special treatment of lone pairs. In the zwitterion, the NH 3 group is well-described by a charge on N and three dipoles on hydrogens, similar to methyls in pentane. The NH 2 group in glycine lacks one site, but has extra electron density associated with the lone pair. We thus assign a dipole moment to N, in addition to its charge, restricted along the sp 3 direction of the lone pair.
Similarly, each oxygen has two lone pairs. Two dipoles for the two lone pairs sum to just one dipole along the symmetry axis, leading to only one variational parameter. The potential of this single dipole, however, is axially symmetric, whereas the potential around the oxygen deviates from axial symmetry due to the particular orientation of the lone pairs. Such a deviation can be accounted for with a quadrupole moment on the oxygen. The finite system of charges sketched in the inset in Fig. 2 shows what is needed. Computing a multipole expansion of three charges we describe two lone pairs with two scalar parameters, a dipole µ r restricted along the symmetry axis, and a restricted quadrupole θ r which has angle β as a fixed parameter (β = 120 • for O= and 109.7 • for O-, due to sp 2 and sp 3 hybridization respectively). θ r and µ r are chosen negative with "−" pointing outside. Table  2 shows a clear advantage of such a multipole set.
To make the definition more transparent, θ r can be expressed in the conventional form as a carefully crafted combination of θ zz = θ r (3 cos 2 β/2 − 1), θ yy = θ r (3 sin 2 β/2 − 1), θ xx = −θ r , and θ xy = θ xz = θ yz = 0, which depend on a single scalar parameter θ r . Here z is along the symmetry axis, and the lone pairs are in the yz-plane. A single restricted quadrupole of strength θ r creates the potential φ quad (r) = θ r 3(r · n 1 ) 2 + 3(r · n 2 ) 2 − 2r 2 2r 5 , where n 1 and n 2 are the directions of the lone pairs. For the zwitterion, MAME does not require θ r on the oxygens, because of the resonance. Lone pairs in the sp 2 and sp 3 configurations lie in perpendicular planes, virtually destroying any asymmetry. MAME accuracy improves away from the molecule (Fig. 2). If f is too large (10 −3 ), there is a net charge inside S which strongly affects the asymptotic behavior. This can be repaired by fixing the total charge using a Lagrange multiplier. 16 Figure 2 demonstrates MAME's insensitivity to choice of f , provided the total charge is correct. Ratio of (φapprox − φ) to φ, both square-averaged over points with a given density ρ, which is an inverse measure of distance. The inset illustrates two lone pairs on oxygen represented with a combination of restricted dipole µr and quadrupole θr.
Last, we analyze the water molecule. The first rule leads to two dipoles on hydrogens pointing along the bonds, sketched in Fig. 3. The charge on oxygen is zero because the molecule is neutral. The dipoles are equal due to symmetry and require no calculation since their vector sum must yield the dipole moment of water, 1.847D (B3LYP/aug-cc-pVTZ value). This already reduces the error to 21%, from 45% with a single dipole on oxygen. The two distributed dipoles yield Θ xx − Θ yy = 4.06 DÅ for the quadrupole moment of water, whereas experiment 19 gives 5.12 DÅ. We note that Θ xx − Θ yy is the only invariant combination of quadrupole components, since the finite dipole makes them dependent on the center of coordinates. We next add µ r and θ r multipoles on the oxygen to describe the correction due to lone pairs, and we find excellent accuracy eσ = 0.59 mH (< 3%) on and beyond S (1.56Å from H and 2.11Å from O, f = 10 −4 ). The sp 3 choice of β = 109.5 • in θ r is crucial: accuracy deteriorates dramatically (to 12%) when β is changed to e.g. 180 • (θ r replaced with θ yy ). The sp 3 description of the oxygen lone pairs is appropriate due to invariance under unitary rotations of occupied orbitals. In conclusion, molecular fields are represented to chemical accuracy with a minimal set of atomic multipoles carefully chosen based on the Lewis structure of the molecule. All H atoms are represented as dipoles. Lone pairs are treated with extra multipoles, avoiding additional off-nuclear expansion sites. 20 The scheme yields multipole values that conform to chemical intuition, are unique, fully rotationally-invariant and free of sampling errors.
We thank K. Krogh-Jespersen, Z.G. Soos, M. Lobanov, and R. Gaudoin for enlightening discussions. This work was funded by AFOSR and the New Jersey Commision on Science and Technology.