Atomic Electron Populations by Molecular Orbital Theory W. J. HEHRE R. F. STEWART AND J. A. POPLE Dept. of Chemistry Carnegie-Mellon University Pittsburgh Pennsylvania 15213 U.S.A. Received 20th September 1968 The detailed distribution of electrons in molecules can bz broken down into atoiiiic orbital populations if an atomic orbital basis is used for a molecular orbital wave function (LCAO). A total assignment of all electrons in a molecule to atomic orbitals (and hence to particular atoms) can be made by finding gross populations as suggested by Mulliken. These are of some interest in discussing various molecular properties associated with charge density particularly if a minimal basis is used consisting of only atomic orbitals in the valence and inner shells. It has now become possible to evaluate atomic orbital populations for a range of molecules both by semi-empirical and ab initio methods.The semi-empirical methods (based on zero differential overlap) utilize some experimental atomic data but can be applied to large systems with little difficulty. Ab initio methods generally require the evaluation of a large number of difficult integrals but they are becoming available for many small systems. The main aim of this paper is to compare atomic electron populations obtained by the semi-empirical INDO molecular orbital theory with some values from full SCF calculations using a minimal basis set of exponential-type functions. Agreement is moderately good but only if the exponents in the full calculations are chosen carefully to minimize the total energy. INTRODUCTION Molecular orbital wave functions and associated electronic charge distributions in linear combination of atomic orbital (LCAO) theory may be obtained either by approximate semi-empirical methods or by full ab initio computations.Semi- empirical methods frequently make use of experimental atomic data in determining parameters and attempt to determine properties of molecules in terms of those of the constituent atoms. Such theories by virtue of their approximate nature have the advantage of being relatively easy to apply to quite large molecules of chemical interest. The full quantum-mechanical calculations however involve much more computation and are only readily available for small molecules but they do lead to well-defined upper bounds for the total energy. Some of the approximate methods make use of the ab initio methods for calibration purposes.Clearly a comparative study of the prediction of these various kinds of wave functions for some simple molecules will be of value in illuminating their deficiencies and pointing the way to more satisfactory theories. This paper is primarily concerned with a comparison of the electron distributions given by the semi-empirical INDO scheme and by full LCAOSCF calculations using a minimal basis set of Slater-type atomic orbitals. 2.-C 0 M P UT A T I 0 NA L D ETA1 L S In all the wave functions discussed here molecular orbitals $i are written as linear combinations of atomic orbitals @/, 15 16 ATOMIC ELECTRON POPULATIONS where 4!t are Slater-type exponential functions. If we deal with a minimal basis set and no atoms heavier than fluorine these have the form 41s(Cl,O = (c;/4* exp (-Cd 42S(C2,0 = (C337+ r exp (-C2r) 42pK2,r) = cts; In>* y exp ( - 123) cos 0 (2.2) and similar expressions for the other 2p functions.For closed-shell ground-states of diamagnetic molecules the electron density (corresponding to a single-determinant wave function using the molecular orbitals (2.1)) is where Ppv is the density matrix occ The overlap integral is the total magnitude of the overlap charge distribution $P+v and if this is " assigned " equally to the orbitals p and v the Mulliken gross atomic orbital population is then a measure of the total number of electrons associated with the orbital 4 in the molecule.2 A gross atomic population QA may be obtained by summing the q,L for all orbitals on a particular atom A.These populations then give a theoretical distribution of electron density throughout the molecule. and full details will not be given here. Only valence electrons are considered explicitly so that two electrons are assigned to inner shell 1s orbitals and these are included in an effective core. The scheme is formally based on a set of Slater-type atomic orbitals with values for the C-exponents chosen by Slater's rules for atoms (except hydrogen for which 5 = 1-2). However the atomic electro-negativities are obtained from experimental atomic energy levels rather than by computation so that some of the deficiencies of a calculation using a limited basis set are avoided. The ab initio calculations obtain the LCAO coefficients in (2.1) by minimizing the calculated total energy.This leads to the Roothaan self-consistent equations involving a set of one- and two-electron integrals. To simplify the evaluation of these integrals the exponential functions (2.2) are replaced by linear combinations of K gaussian functions a procedure first suggested by Foster and Boys.' Thus The semi-empirical INDO method has been specified elsewhere W. J. HEHRE R. F. STEWART AND J . A . POPLE 17 Here yls and g2p are the gaussian-type functions sls(a,r) = (2a/.Y exp (-ar2) gZp(a,r) = (128a5/n3)4r exp (- ar2) cos 0. (2.9) (The 2s exponential function is written as a linear combination of 1s gaussians). The constants d and a are chosen to minimize the integral (2.10) subject to it normalization constraint and are independent of the c-exponents. This procedure was suggested by 0-ohata Taketa and Huzinaga6 Calculations have been carried out for K = 3 4 5 and 6.Values of d and a for K = 3 4 and 5 were obtained by methods described elsewhere and are given in table 1. Values for K = 6 were taken from ref. (6). (Some preliminary calculations have already been published for K = 4 using the exponents and coefficients proposed by 0-ohata Taketa and Huzinaga. However the values given in table 1 lead to lower &-values. Also the E value for the 2s function with K = 5 is less than the value with K = 6 given in ref. (6)). Using these expansions the various integrals can be reduced to integrals involving gaussian functions and these can be evaluated by methods originally introduced by The self-consistent equations were solved in an iterative manner.Sub- routines from programmes QCPE 47 and QCPE 92 were used in the computation.1°9 '' Calculations have been carried out both for a standard set of [-exponents (using Slater's rules but with = 1.2 for hydrogen) and by optimizing the valence 5's to give a minimum calculated energy. For the latter inner shell cl were chosen equal to optimum atomic values l 2 (rounded to two decimals) and the valence i2 (c for hydrogen) were varied in steps of 0.01 until the lowest energy was found. 3.-RESULTS AND DISCUSSION The first point to study is the efficiency with which a gaussian expansion of the type used in this paper reproduces the results of a calculation based directly on Slater exponential orbitals. Table 2 gives results for H F and NH with K = 3 4 5 and 6 using the same nuclear geometry and orbital c-values as the original exponential-type calculations reported in the literature.' 3 9 l4 Convergence to the exponential result with increasing K is fairly rapid for energies dipole moments and populations.The differences between the gaussian and exponential results are much less for the binding energy (difference between molecular and atomic values using the same exponents) than for the total energy. This suggests that much of the difference arises from the atomic inner shell description which is largely unaltered in the mole- cule. Atomic populations with the 3-gaussian functions differ from the exponential values by only about 0.02. This suggests that this simplest level of calculation is sufficient to study the relation between the ab initio and semi-empirical populations.Results for a series of simple molecules and ions with K = 3 are listed in tables 3-9. For the neutral molecules the nuclear geometry is chosen according to the standard model A defined previo~s1y.l~ For the ions model B is used l6 except for CH; where the HCH angle is taken as tetrahedral. The following points may be noted about these results. I . The optimum values of the c-exponents vary considerably from one molecule to another. This is most marked for hydrogen where the range is 1.14-1.48. For carbon in all neutral molecules and positive ions the optimum exponent is substantially K als 3 0.109818 0.405771 2,227660 4 0.088019 0.265204 0.954620 5.21 6850 5 0,074453 0.197573 0.578651 2.07 1740 11.30570 TABLE 1 .-COEFFICIENTS AND EXPONENTS FOR GAUSSIAN FIT OF SLATER ORBITALS dlS &Is 012s d2s E2s U2P d2P 0.444635 0.0601 83 0,458 179 0*080098 0.422307 0.535328 3.31 x 0.156762 0.596039 6.87 x lo-’ 0.235919 0.566171 0.1 54329 2.581 580 - 0.059945 0.919238 0.1 62395 0.291 626 0.061257 0,477008 0.065439 0.26323 1 0.532846 4.38 x lo-’ 0.1 60728 0.580559 2 .7 0 ~ lo-’ 0.164372 0.551787 0.2601 41 2*000240 - 0.054721 0.466262 0,285746 0.056752 11.615300 -0.01 1984 1.798260 0.0571 32 0.193573 0.04249 1 0.1 47266 0.056063 0.1 65345 0.482571 0.088064 0.54805 1 0.127331 0.483493 0.331815 6.88 x 0-194473 0.369824 2.70 x 0.307982 0.366774 0.113540 1,673710 - 0.056859 0.864327 0.1 23547 0.0221 40 8,984940 - 0.01 5964 3.320390 0.020790 TABLE 4.-cOMPARISON OF ATOMIC POPULATIONS Q FOR NEUTRAL MOLECULES H,X,Y QH Qx QY standard optimized INDO standard optimized INDO standard optimized 5.909 5.903 5.971 8 *09 1 8.097 0.866 0.755 0.725 9.1 34 9.245 9.275 0.850 0.781 0.829 8.301 8.437 8.342 0.849 0.832 0.912 7.452 7-503 7.264 0-87 1 0.978 1 -009 6.51 8 6.088 5.963 0.877 0.979 1.018 6.370 6.063 5.946 5.997 0-866 0.917 1.001 6-267 0.816 0.813 0.947 6.184 6.1 87 6.053 0.886 1.001 1.031 6.299 5.838 5.676 9.042 0.890 0.973 1 a045 6-155 5.823 5-672 8.065 0.790 0.792 0.950 6.145 6.073 5.91 1 7.064 0.848 0.829 0.902 7.302 7.277 7.027 9.001 8.022 0.858 0.824 0.948 7.120 7-093 6-928 6.1 66 HOF 0.840 0.750 0.789 8.172 8.227 8.094 8.987 2P 2 .6 9 ~ 10-4 2.90 x 10-5 3 . 7 2 ~ INDO 8.029 9-1 59 9.23 1 8.230 8-238 7.1 34 7.1 39 9.065 9.169 8.082 8.124 9-023 9.1 17 W . J. HEHRE R. F . STEWART AND J. A . POPLE 19 greater than the standard Slater value for the free atom (1.625).The same applies to nitrogen except for the negative ion NH;. 2. The atomic populations after optimization of exponents differ considerably from those using standard exponents. For the neutral molecules the optimized calculations lead to a wider variation of hydrogen populations corresponding to a greater variation in polarity. For both neutral molecules and ions the values of Q H for the optimized set agree better with the INDO charges than do the standard TABLE 2 total energy binding energy dipole moment atomic populati2,ns molecule basis set (hartrees) (hartrees) (D) H 3 gaussians - 98.5072 0.0484 0.815 0.866 4 gaussians - 99.1922 0.0476 0.870 0.851 5 gaussians - 99.3893 0.0468 0.818 0.847 6 gaussians - 99.4656 0.0468 0.88 1 0.846 Slater b - 99.4785 0.0468 0.878 0.846 3 gaussians - 55.4453 0.3055 1-729 0.854 4 gaussians - 55.8372 0.3013 1 -763 0.846 5 gaussians - 55.9522 0.3006 1-767 0-845 6 gaussians - 55.9974 0.3001 1 a765 0.845 Slater - 56.0052 0.3001 1.764 0.845 a ref.(6); b ref. (13); c ref. (14). A 9.1 34 9.149 9.153 9.1 54 9.154 7.439 7.462 7.465 7.466 7.467 TABLE 3.-oPTIMIZED EXPONENTS FOR NEUTRAL MOLECULES HmXnY 5iW) 1.19 1.33 1 a28 1 *25 1.18 1.18 1 -23 1.31 1-19 1.21 1.37 1.28 1-31 1.33 12(C) 1 -65 1.75 1 -74 1.69 1.67 1 -78 1.75 1 -69 1 -97 2.26 2.28 2.22 2.55 2.53 1.94 2.5 5 2-23 1.95 1.96 2.56 1 -96 2-25 2.23 2.55 set. For the entries listed in tables 4 and 8 the root mean square values of QH- QH(IND0) are 0.109 for the standard set and 0.084 for the optimized set. 3. Dipole moments calculated with the optimized exponents are in better agreement with experimental values than are those with standard exponents.However calculated values still show insufficient polarity in most cases and the INDO results are superior. 4. The optimum exponents for hydrogen are linearly correlated with the electron population of the atom. As the atom becomes more positive the orbital becomes more contracted. This is illustrated in fig. 1. 20 ATOMIC ELECTRON POPULATIONS 5. Many detailed features of INDO charge densities are reproduced by the optimized calculations but less satisfactorily by the standard set. In the paraffins CH4 and C2H6 the optimized results and INDO give a nearly neutral hydrogen whereas the standard exponents lead to considerable positive character. Both the TABLE 5.-cOMPARISON OF DIPOLE MOMENTS FOR NEUTRAL MOLECULES HmXnY HmXnY co HF H20 H3N HjCF H2CO HCN H2NF HNO HOF standard 0.8 15 0-985 1.520 1,602 0-651 0-922 2.065 1353 1 -343 1 -275 optimized 0.730 1 -465 1.841 1 -677 0.981 1 ~626 2.457 1-783 1.601 1-714 INDO 0.941 1.966 2.1 34 1-876 1.692 1 -999 2.455 2.433 1.756 2.244 expt.0-13 a 1.8195 1.846 C 1.468 d 1-855 2.339 f 2.986 g QA. L. McClellan TabZes of Experimental Dipole Moments (W. H. Freeman and Co. San Fransisco b R. Weiss Physic Rev. 1963 131 659. C G. Birnbaum and S. K. Chatterjie J. Appl. Physics 1952 23 220. d D. K. Coles W. E. Good J. K. Bragg and A. H. Sharbaugh Physic. Rev. 1951 82,877. e M. Larkin and W. Gordy J. Chem. Physics 1963 38,2329. f J. M. Shollery and A. H. Sharbaugh Physic. Reu. 1951 82,95. 9 B. N. Bhahacharya and W. Gordy Physic.Rev. 1960,119 144. Calif. 1963) p. 48. TABLE 6.-TOTAL ENERGIES FOR NEUTRAL MOLECULES H,X,,Y total cnergy standard optimized - 1.1191 - 1.1 192 - 107.4868 - 107.4884 - 147.6193 - 147.6197 - 195.9269 - 195.941 6 - 111.2186 - 111.2176 - 98.5432 - 98.5588 - 74.9437 - 74.95 18 - 55443 1 - 55.4441 - 39.7102 - 39.721 6 - 78.2692 - 78.2947 - 77.0444 - 77.0584 - 75.831 6 - 75.8459 - 137.1213 - 137.1514 - 1 12.3205 - 112-3379 - 91.6521 - 91.6620 - 152.841 1 - 152.8480 - 128.0387 - 128,0426 - 172.3290 - 172.3430 optimized set and INDO lead to a marked increase of positive charge on hydrogen along the series CH4 NH3 H20 and HF. The same is true for the hydrocarbon series C2H6 C2H4 C2H2. In methyl fluoride the hydrogen is more negative than in the isoelectronic hydrocarbon ethane for all three calculations.The same applies to the formaldehyde-ethylene pair. This provides some further backing W. J. HEHRE R. F. STEWART AND J . A . POPLE 21 for the alternating inductive charge displacements in a-systems suggested previously on the basis of CND0/2 calc~lations.~~ Among the ions the large negative charge on N in NH; predicted by the standard calculation is much reduced when optimized exponents are used. INDO predicts the nitrogen to be approximately neutral. TABLE 7.-OPTIMIZED EXPONENTS FOR IONS (XH,) * XHni i 1 (HI M C ) 4XN) MO) 52(F) H 1 -40 CHf 1 *41 1 072 CHZ 1-33 1.77 CH 1.15 1.53 NH+ 1 -44 2-04 NH; 1.36 2.05 NHZ 1.33 2-06 NH 1.16 1.80 OH+ 1 -45 2.35 OH 1.41 2.33 OH 1-37 2-31 OH- 1-14 2.09 FH+ 1 *48 2.66 TABLE 8.-cOMPARISON OF ATOMIC POPULATIONS Q FOR IONS (XH,)' ion CH+ CH CH; NHf NH; NH NH OH+ OH; OH,+ OH- FH+ standard 0.695 0-674 1 so87 0.594 0.621 0.637 1-124 0-536 0.572 0.602 1-239 0.488 QH optimized 0-705 0.756 1.012 0.637 0.689 0-71 5 0-967 0.583 0.604 0.621 1.017 0.529 INDO 0.848 0.849 1.131 0.697 0-734 0.752 1.108 0.562 0.601 0.638 1.113 0-454 standard 5.305 5.978 6-739 6.406 7,136 7.454 7-75 1 7.464 7.856 8.195 8.761 8.512 ex optlm1zed 5.295 5.731 6.965 6.363 6.934 7.140 8.066 7.417 7.793 8.138 8-983 8.471 TABLE TO TOTAL ENERGIES FOR IONS (XH,)* ion XHm Hi- CH+ CH CH; NHf NHZ NH NH OH+ OH; OH$ OH- FH+ electronic state "$ IZ+ 1A 1 ' A 1 2rI 2A1 1A 1 l A 1 3c- 3B1 total energy standard optimized - 0.5307 -0.5510 - 37.41 75 - 374476 - 38.7270 - 38.7821 - 38.8748 - 38.8957 - 53.8128 - 53.8452 - 55.1 626 - 55.2008 - 55.8451 - 55.8825 - 54.5625 - 54.621 1 - 74.0397 - 74.071 1 - 74.6549 - 74.6823 - 75.3162 - 75.3369 - 74.01 89 - 74.1 276 - 98.1 934 - 98.221 8 INDO 5.1 53 5.454 6.606 6.303 6.799 6.994 7.784 7.43 8 7.797 8.085 8.887 8.546 22 ATOMIC ELECTRON POPULATIONS Overall we may conclude that the ab iizitio calculations reported here do provide some support for the charge distributions predicted by the simpler semi-empirical method.However it is most important to optimize exponents if such a comparison is to be made. The correlation between optimum exponents and electron populations suggests that exponent variation is necessary to give an adequate account of the + 1.30 0 * 1.20 11h3 i ~ l . l ~ l ~ l ~ l ~ - 4 0 *SO * 6 0 - 7 0 -80 .90 1.00 1-10 Atomic population FIG. 1.-Atomic Population of Hydrogen in Various Molecular Environments as a function of Cls(H).basic electronegativities of the atomic orbitals. No exponent variation is used in the semi-empirical approach but this kind of effect is partly allowed for by the use of empirical electronegativity parameters. Valuable discussions with Dr. M. D. Newton are acknowledged. This research was supported by U.S. Army Research Office-Durham Grant DA-ARO-D-3 1-124- G722 and National Science Foundation Grant GP-8472. J. A. Pople D. L. Beveridge and P. A. Dobosh J. Chem. Physics 1967 47 2026. R. S. Mulliken J. Chem. Physics 1955,23 1833 1841 2338,2343. J. C. Slater Physic Reu. 1930 36 57. C. C. J. Roothaan Rev. Mod. Physics 1951 23 69. J. M. Foster and S. F. Boys Rev. Mod. Physics 1960 32 303. K. 0-ohata H. Taketa and S. Huzinaga J.Physics SOC. (Japan) 1966 21 2306. R. F. Stewart Small Gaussian Expansions of Atomic Orbitals J. Chem. Physics to be published. W. J. Hehre and J. A. Pople Chem. Physics Letters. 1968 2 379. S. F. Boys Proc. Roy. SOC. A 1950,200 542. lo I. G. Csizmadia M. C. Harrison J. W. Moskowitz S. Seung B. T. Sutcliffe and M. P. Barnett Quantum Chemistry Program Exchange (Indiana University Bloomington Indiana U.S.A.) l 1 D. R. Davis and E. Clementi Quantum Chemistry Program Exchange (Indiana University Bloomington Indiana U.S.A.) l2 E. Clementi and D. L. Raimondi J. Chem. Physics 1963 38 2686. l 3 B. J. Ransil Rev. Mod. Physics 1960 32 239 245. l4 W. E. Palke and W. N. Lipscomb J. Amer. Chem. Soc. 1966 88 2384. Correction given in l5 J. A. Pople and M. Gordon J. Amer. Chem. Soc. 1967,89,4253. M. D. Newton and W. N. Lipscomb J. Amer. Chem. Soc. 1967,89,4261. J. A. Pople D. L. Beveridge and P. A. Dobosh J. Amer. Chem. SOC. 1968 90,4201.