Minimal Atomic Multipole Expansion for Molecular Electrostatic Potentials: MAME, Lab Reports of Chemistry

Download Minimal Atomic Multipole Expansion for Molecular Electrostatic Potentials: MAME and more Lab Reports Chemistry in PDF only on Docsity! Rules for minimal atomic multipole expansion of molecular fields E. V. Tsipera) Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854; Center for Computational Material Science, Naval Research Laboratory, Washington, DC 20375; and School of Computational Sciences, George Mason University, Fairfax, Virginia 22030 K. Burke Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, New Jersey 08854 ~Received 24 October 2003; accepted 21 November 2003! A nonempirical 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. © 2004 American Institute of Physics. @DOI: 10.1063/1.1640995# The problem of representing the electrostatic potential outside a molecule using atomic charges or higher atomic multipoles is very important for understanding intermolecu- lar 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 Mulliken1 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 potential-derived ~PD! charges, which avoid representation of the density by producing the ‘‘best’’ fit to the molecular potential directly.6,7 Atomic dipoles and quadrupoles7 are often used to increase accuracy in solvation problems8 and force field calculations.9 Induced atomic di- poles appear naturally in electronic polarization of molecular solids10 to account for the small part of molecular polariza- tion 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 CHelpG13 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 counter-intuitive, such as nega- tive charges on hydrogens in alkanes.15 Higher multipoles only increase the redundancies inherent in distributed multi- pole 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,14,16 or by introducing restraints.17 Our approach does not use sampling and eliminates redundancies before they appear. We approximate the true molecular potential, f~r!, as a sum of multipoles of strength qk centered at nuclear posi- tions ri , f~r!'fapprox~r!5( i ( k qkfk~r2ri!, ~1! where fk(r) is the potential due to the kth multipole of unit strength: fk(r)51/r for charge, (n"r)/r 3 for a dipole in the direction n, and so on. Since ¹2fapprox50 everywhere ex- cept at ri , but ¹ 2f54pr(r), the atomic multipole expan- sion can only be accurate in regions where r~r!'0. Further- more, f on any closed surface S on which r50, determines f~r! everywhere outside S. We therefore choose S to be an isodensity surface, r(r)5 f , where f is sufficiently small to ensure negligible charge beyond S, but with sufficient poten- tial on S for a determining fit ~Fig. 1!. We chose fapprox to minimize s25S21 R S dS@fapprox~r!2f~r!# 2 ~2! over S which leads to a system of linear equations (kCmkqk5bm , where Cmk5S 21 R S dSfm~r!fk~r!, ~3a! bm5S 21 R S dSfm~r!f~r!. ~3b! Atomic multipoles defined in this way are fully rotationally invariant, which is an issue with some PD schemes.13 The error s can be compared to f̄ , f̄25S21 R S dSf2~r!. ~4! 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 ev- erywhere 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 ap- plication in biochemistry; and water, to see how general- MAME rules apply to a small polar molecule. All den- sities and potentials are produced on a cubic mesh by the GAUSSIAN 98 program,18 at the B3LYP/aug-cc-pVTZ level JOURNAL OF CHEMICAL PHYSICS VOLUME 120, NUMBER 3 15 JANUARY 2004 11530021-9606/2004/120(3)/1153/4/$22.00 © 2004 American Institute of Physics Downloaded 31 Mar 2004 to 128.6.71.63. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp ~6-31111G** for pentane!. Surface integrals ~3! are com- puted by triangulation of S. The program runs within a few seconds, and is available on request. Figure 1 shows f~r! on S for n-pentane. Red spots ~f.0! show an excess of positive charge near each hydro- gen, 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 in- wards. 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 pro- ton. This unique property of hydrogens is well-known in x-ray structure analysis, which systematically underestimates the C–H bond lengths for this reason. Mulliken charges are intuitively meaningful but produce large errors in the potential ~Table I!. PD charges are nega- tive on some hydrogens and still give significant errors. Add- ing 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 a.u.! 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 mol- ecule 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 ~52%! in the glycine zwitterion, (NH3) 1 – CH2 – COO 2, down from 4% with charges alone and 4%–6% with stan- dard PD schemes. The glycine zwitterion is highly polar with dipole m510.3 D, which MAME recovers within 0.1% ac- curacy. Table II lists MAME results for glycine in its standard form, NH2 – CH2 – COOH, and illustrates the need for special treatment of lone pairs. In the zwitterion, the NH3 group is well-described by a charge on N and three dipoles on hydro- gens, similar to methyls in pentane. The NH2 group in gly- cine 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 sp3 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 sym- metry 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. Com- puting a multipole expansion of three charges we describe FIG. 1. ~Color! Electrostatic potential over isodensity surface S of pentane. Here f 5531024 a.u. produces S at ;1.4 Å from the hydrogens and leaves 20.2e charge outside. Missing charge is negligible for f 51024, with S at 1.8 Å. TABLE I. Partial atomic charges in n-pentane. f 5531024 a.u., error s as in ~2!, (%)5(s/f̄), ef̄53.6 mHartree ~1 mH527 meV;kT at 300 K!. ‘‘m’’ indicates atomic dipoles, ‘‘mr’’—dipoles restricted along H–C bonds. Method q~H!, range q~C!, range es , mH ~%! CD charges Mulliken 10.11...10.14 20.59...20.11 9.3 ~260! ZINDO 10.03...10.04 20.15...20.04 3.5 ~99! PD charges CHelp 20.04...10.04 20.11...10.15 3.5 ~97! CHelpG 20.04...10.04 20.16...10.16 3.2 ~87! MK 20.03...10.06 20.22...10.13 3.1 ~86! PD charges plus dipoles CHelp1m 20.76...10.10 20.66...12.05 2.8 ~78! CHelpG1m 20.32...20.30 10.65...10.86 1.8 ~51! MK1m 20.27...20.20 10.50...10.64 1.8 ~49! MAME Charges 20.01...10.09 20.34...10.13 2.6 ~72! m~H! ~m50.07...0.09! 20.01...10.03 0.5 ~15! mr(H) ~m50.06...0.09! 20.02...10.01 1.6 ~45! TABLE II. MAME for glycine without and with lone pair multipoles. f 51024 a.u., ef̄515 mH. NH2 CH2 C5 5O –OH es , mH ~%! mr(H) 20.05 20.03 10.79 20.54 20.17 4.1 ~27! mr(H) 10.06 10.11 10.37 20.67 10.13 1mr(N,O) 20.65 10.19 20.71 1.6 ~11! 1ur(O) 20.53 20.87 1154 J. Chem. Phys., Vol. 120, No. 3, 15 January 2004 E. V. Tsiper and K. Burke Downloaded 31 Mar 2004 to 128.6.71.63. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp