- Email: [email protected]
The electrostatic potential in the semiempirical molecular orbital approximation
29 July 1994
ELSEVIER
CHEMICAL PHYSICS
Chemical Physics Letters 225 (1994) 11-17
The electrostatic potential in the semiempirical molecular orbital approximation Peter L. Cummins *, Jill E. Gready Department ofBiochemistry, University of Sydney, Sydney, NSW 2006, Australia Received 15 March 1994; in final form 5 May 1994
AbStIYICt A method is described for calculating the molecular electrostatic potential (MEP) within the MNDQ, AM1 or PM3 semiempirical approximations. The MEPs are used to obtain atomic charges from which dipole moments are calculated for comparison with the exact semiempirical values. The results for a wide range of neutral molecules and ions show excellent agreement between the dipole moments calculated from MEPderived atomic charges and the exact semiempirical dipole moments. The neutral molecules and cations showed a mean deviation in the dipole of less than 0.0 1 D. The maximum errors ( z 0.04 D) were obtained for some molecules containing third row or heavier elements. By contrast, all previously reported methods based on tbe semiempirical approximations for generating MEP-charges in general yield significantly larger deviations from tbe exact dipoles. The maximum error for the anions tested was 0.07 D, with a mean deviation only slightly higher than for neutral species and cations. The correct limiting behaviour is obtained when the MEP surfaces are calculated at large distances, i.e. tbe exact dipole is approached asymptotically. Quadrupole moments of the MEPderived charge distributions are reported for some nonpolar and polar molecules.
1. Introduction The molecular electrostatic potential (MEP) is defined in terms of the interaction energy between a molecule and a unit positive charge. If this point charge is exterior to the molecule, i.e. in a region where the molecular electron density is relatively small, the MEP may be used to rationalise chemical reactivities [ 1 ] or to obtain the atomic charges used in molecular dynamics or Monte Carlo simulations of biomolecular systems [ 2-41. Quantum mechanics provides a rigorous definition of the electrostatic potential as the expectation value of the Coulomb operator. In the density functional and ab initio formalisms, the calculation of MEPs is straightforward from
the molecular orbitals expressed in terms of a suitable basis set. However, in the semiempirical methods (MNDO, AM 1 and PM3) [ 5-91, properties are not computed as integrals over basis functions, but simplified functional forms are used to approximate the kinetic energy, one-electron and two-electron integrals required to evaluate the total energy within a minimal-basis HF formalism [ 5, lo]. As the basis is not strictly defined for the calculation of properties, it is not surprising that a number of diverse approaches have been proposed for the evaluation of the MEP and MEPderived atomic charges using semiempirical methods [ 11-l 61. In addition, several approaches for combining MNDO, AMI and PM3 with molecular mechanics (QM/MM) or continuum self-consistent reaction field (SCRF) meth-
* Correspondingauthor. E-mail: [email protected] 0009-2614/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIOOO9-2614(94)00617-Y
P.L. Cummins, J.E. Gready/Chemical PhysicsLetters225 (1994) I l-l 7
12
ods have also been proposed in the effort to model condensed phases [ 17-241. In this Letter, we have attempted to define the molecular electrostatic potential in a manner consistent with all approximations involved in the semiempirical method. To this end, we derive simplified functional forms for the electrostatic potential based on the core-electron attraction terms used in MNDO, AM1 and PM3. The dipole moments calculated from the MEP-derived atomic charge distributions are compared with the exact dipole moments as a test of the method, and also with other methods [ 1l-141 for obtaining MEP-derived charges.
2. Method The neglect of diatomic differential overlap (NDDO) approximation [ 25-271 on which AMl, MNDO, and PM3 are based requires that all oneelectron terms arising from diatomic overlap be neglected. Consequently, there are no three-centre contributions from atomic orbital products Q& where ,u and Yare centred on different atoms, as there would if the electrostatic potential were computed from a suitable finite-basis set expansion [ 11,15 1. Thus in the NDDO approximation, the MEP can be expressed solely in terms of atom-centred contributions (all equations given in atomic units), vi = C PpvINv + Zil& w
3
(1)
where Zi is the core charge of atom i, R, the distance from atom i, Ppy are elements of the first order density matrix, and Ipy are the one-electron integrals of the Coulomb operator involving only those orbitals fiy and @,,centred on the ith atom, i.e. Zpv=
s
I&-R’1 -‘$v(R’)@p(R’) CN’>
(2)
where Riis the position vector. All the semiempirical methods (MNDO, AM 1 and PM3 ) make use of twoelectron repulsion integrals of the form (ss, pv) [ 5, lo] in order to approximate one-electron integrals such as those given by Eq. (2 ). Since this approximation underestimates the core-electron attractions, the core-core repulsions must likewise be reduced from their Coulomb values in order to main-
tain the balance between the attractions and repulsions in the molecule. Purely empirical functions of the interatomic distances are then added to the corecore repulsion terms in order to tine-tune this balance [ 561: these have no explicit role in defining the electrostatic potential. In terms of the two-electron integrals, Eq. ( 1) may be cast in the form
where the subscript m indicates the transformation of (ss, ,uv) from the molecular Cartesian to local (diatomic) coordinate system. The total MEP at a point is simply given by the sum of V, over all atoms. However, since (ss, VP) underestimates core-electron attractions they are unsuitable for calculation of the MEP. In our interaction energy method [ 141, the Coulomb interaction terms in the fitting procedure used to obtain atomic charges were replaced by the semiempirical equivalent, i.e. (ss, ss), in order to compensate approximately for this error. While this method lead to atomic charges and molecular dipole moments no worse than those obtained by Ferency et al. [ 131, its deficiency was that an electrostatic potential could not be rigorously defined. However, simple analytic forms for the molecular electrostatic potential may be obtained by considering the approximations involved in the evaluation of the two-centre electron repulsion integrals (ss, ,uv) in MNDO, AM1 and PM3. Classically, the two-centre electron repulsion integrals (pv, na) are equivalent to the interaction energy between the continuous charge distributions @,&,and e&p In semiempirical treatments, these charge densities are replaced by sets of discrete charges [ 5,10 1. At large internuclear separations, i.e. large compared with the dimensions of the point-charge arrangement used to represent the atomic orbital products, the interactions between multipoles are accurately given by the sum of Coulomb interactions between charges. Since these interaction terms may diverge as internuclear separations are reduced, a modified Coulomb interaction is introduced which gives the correct one-centre limit, i.e. the semiempirical one-centre two-electron integrals. In MNDO, AM 1 and PM3, the Dewar-Sabelli, Ohno, Klopman (DSOK) equation [ 28-3 1 ] is used to calculate the interaction energies between charges that represent the multipoles of @,+, and &&,. If R is the
P.L. Cummins, J.E. Gready /Chemical Physics Letters 225 (1994) 1l-l 7
distance between any pair of these charges (one associated with nucleus X, the other with Y ) , the DSOK equation is given by a screened Coulomb potential of the form f(R)=
(Rz+Cuz) -1’2,
WI
a=px+py .
(4b)
Thus p are parameters chosen to reproduce the onecentre repulsion integrals (for X =Y) at zero internuclear separation. Consequently, a simplified analytic form for the Coulomb potential may be recovered by setting a! =px =pr = 0 in the two-electron integrals (ss, ,uv). Luque et al. [ 23 ] have tested both cw=O and cr=px in a self-consistent reaction field (SCRF) formalism for modelling solvent effects and found that the former choice gave heats of formation of various solutes in better agreement with experiment, despite their assertion that the latter is the more rigorous. Chudinov et al. [ 2 1 ] point out correctly that cr=O yields a’true’ electrostatic potential in their SCRF treatment. This choice was also used by Wang and Ford [ 201 in their solvation studies, but strangely not in their recent MEP calculations where a more complicated treatment is given [ 16 1. Setting a=0 in (ss, pv), we obtain the following expressions for the one-electron integrals for ( p = px or p,., z = pZorbitals ) : (SS,SS)=l/Ri,
(W
(SS,SZ)=f[(Ri-DI)-'-(Ri+DI)-'1
,
(5b)
(SSyZZ)=f{f[(Ri+2D2)-’ + (Ri_2D2)-‘]
- (SS, SS)} 3
(SS,pp)=f[(Rf+4D~)-“2-(ss,ss)]
(SC)
(5d)
where the parameters D, and D2 are determined by the condition that the dipole and quadrupole moments of each atomic orbital product q!~& are equal to those for the corresponding point-charge distributions. Values for D, and D2 are readily calculated using the formulae [ 10 ] D1 = (2”+ 1) (4C”&)“+1’2 ’ 3’9 &#,+ &)zn+2
(ha)
D~=[&,(2n+1)(n+1)]“2(~,)-‘,
(6b)
where c,,, and & are respectively the exponents of
13
the Slater s and p atomic orbitals for principal quantum number n. It is emphasized that Eqs. (3) and ( 5) would yield the correct Coulomb potential only in regions where the molecule’s electron density is small enough to neglect, i.e. at a point outside the molecule, since then RiB D. For the first and second row atoms, for example, the D parameter values for AM1 range from 0.2 to 0.5 A, compared with typical van der Waals (vdW) radii ranging from 1.2 to 1.7 A. The extent of the electron density will depend on the magnitude of the density matrix elements P#,,. Note also that in the limit of Ri+a, a in Eq. (4) may be neglected and the interaction energy method [ 141 would correspond to the present electrostatic potential method (EPM). In the present work, we have used Eq. ( 5 ) to compute the molecular electrostatic potentials for a wide range of neutral and charged molecules, and derive sets of atom-centred charges by a least-squares fitting procedure. This fitting procedure has been widely used to obtain point-charge representations of molecules and is described in detail elsewhere [ 2,12,14]. Briefly, the electrostatic potential is computed at a large number of points on a number of surfaces surrounding the molecule. A set of atomic charges are then derived which reproduce these potentials in a least-squares sense. Using a number of surfaces ensures a good fit both close to and far away from the molecule. In the present calculations, the atomic radii ( Rk) used in the calculation of the kth surface are givenby (inA) & = c&w +0.3+
1) )
(7)
where r~is a scaling factor, RvdW is the vdW radius of an atom and k= 1,2, 3 and 4, i.e. four surfaces were generated and used in the fitting procedure for obtaining atomic charges. For elements of unknown vdW radii, we have used twice the covalent or atomic radii in the present calculations. Most calculations were performed using a value of Q= 1.5 for scaling the vdW radii in order to be reasonably consistent the other MEP-derived charge calculations [ 11- 14 1. In some cases, we varied G in order to ascertain how well the method performs for surfaces very close to the molecule. A surface density of at least two points per A2 was chosen to give dipole moments to within a precision of f 0.01 D. All calculations were per-
P.L. Cummins, J.E. Gready/Chemical PhysicsLetters225 (1994) II-I 7
14
formed at the AM 1 optimised geometries of the neutral molecules and ions.
3. Results and discussion Here we present AM1 results only, as we reached similar conclusions with test calculations for the other NDDO-based semiempirical methods (MNDO and PM3). In Table 1, the dipole moments calculated from the MEPderived point-charge representation using the electrostatic potential method EPM (i.e. Eqs. ( 3 ) and ( 5 ) ) are compared with the exact AM 1 and experimental values for neutral molecules. Also included in Table 1 for comparison are the dipoles calculated from semiempirical MEP-derived atomic charges reported for the various other methods in the
literature [ 1l-l 41. We have recalculated MEPs for the interaction energy method (IEM) [ 141, since a minor error, which we have now corrected, was detected in the surface generation code. This error lead to dipole moments which differ from the ones given in Table 1 by Y/o- 10%. The MEPs obtained using the EPM and IEM methods were calculated with a value of cr= 1.5 in Eq. (7). The results in Table 1 indicate the IEM consistently predicts dipole moments smaller than the exact values. However, the agreement between the EPM and exact molecular dipole moments is excellent. The mean absolute difference between EPM and exact dipoles for neutral molecules is less than 0.0 1 D, with a maximum deviation not exceeding 0.04 D. On average, the deviations tend to be larger for molecules containing third-row or heavier elements. In contrast, the other methods often yield
Table 1 AM 1 ( MEPderived and exact SCF) and experimental dipole moments (D ) for neutral molecules Molecule
EPM .
IEM ”
HF Hz0 NH, CHsO CH,OH NH&HO CHOOH CH3CN CH,CO HCl FCl
1.73 1.86 1.84 2.33 1.63 3.70 1.48 2.89 1.35 1.40 0.87 1.50 0.08 0.64 4.37 2.36 3.60 1.95 2.05 0.49 1.75 1.88 2.30 5.64 0.61 2.16 3.17 5.30
1.64 1.76 1.76 2.24 1.58 3.70 1.42 2.78 1.39 1.35 0.84 1.45 0.08 0.61 4.24 2.27 3.48 1.88 1.99 0.48 1.69 1.82 2.18 5.54 0.60 2.14 3.13 5.21
CH3Wz
NO NNO ClNO HCN imidaxole pyrrole pyridine furan HsSiC& HsS PH3
ZnO BrI CH3GeN CH,HgBr HSO
Ref. [ 121 1.28 1.36 1.93 1.01 2.75 1.39 2.61
’ Electrostatic potential method (Eqs. (3) and (5) in text). cSeeRef. [ 141 for sources of experimental dipole values.
Ref. [ 131 1.75 1.89 1.75 2.32 1.62 3.68 1.46 2.84 1.38 1.45
Ref. [ 161 1.82 1.74 2.29 1.62 3.60 1.42 2.82
b Interaction energy method, Ref. [ 141.
Exact
Exp. c
1.74 1.86 1.85 2.32 1.62 3.70 1.48 2.89 1.35 1.38 0.91 1.50 0.09 0.64 4.36 2.36 3.60 1.95 2.04 0.49 1.76 1.86 2.29 5.66 0.63 2.19 3.15 5.32
1.82 1.85 1.47 2.33 1.70 3.73 1.41 3.92 1.42 1.11 0.89 1.31 0.16 0.17 2.98 3.8 1.74 2.19 0.66 1.18 0.97 0.58 0.72 3.99 3.10
P.L.. Cummins, J.E. Gready / Chemical Physics Letters 225 (1994) 1l-1 7
significantly larger deviations from the exact dipoles. The deorthogonalization (finite basis set) approach [ 11,15 1, also adopted by Besler et al. [ 12 ] for the calculation of atomic charges, gives by far the largest deviation of 0.95 D for NH&HO, while the IEM [ 141 and Ferenczy et al. [ 13 ] methods give maximum deviations of 0.17 and 0.10 D, respectively. The recently published results of Wang and Ford [ 16 ] also yield a maximum deviation of 0.10 D. It is worth noting that the numbers of molecules tested in these other studies [ 12,13,16] are quite small and the maximum errors may well be larger, particularly for molecules containing heavier elements. The dipole moments calculated using the EPM and IEM at the centre of mass for cations are given in Table 2 along with the exact SCF values. A deviation of 0.02 D is obtained for CH3Zn+ and OPCl$, while the largest deviation of 0.04 D is obtained for ClNO+ and HgCl+. Note that the IEM overestimates the dipole for HCNH+ and FCl*+ whereas all other cations give IEM dipoles which are less than the exact values. The average deviation for cations is much the same as for neutral molecules, i.e. not more than 0.01 D. The EPM, IEM and exact dipoles calculated at the centre of mass for anions are given in Table 3. The largest deviations of 0.06 and 0.07 Dare obtained for Table 2 AM1 (MePderived and exact SCF) dipole moments (D) for cations calculated at the center of mass Cation
EPM ’
IEM b
Exact
H,O+ CH,NHf CHsOH: HCNH+ NCO+ HCl+ imidazole*H+ fluorobenzene.H+ FCl+ FCl’+ ClNO+ H,S+ 0PCl: CH,Zn+ ZnOH+ H&l+
2.33 2.33 2.16 0.50 1.00 2.13 1.63 4.23 1.72 2.22 2.83 1.85 0.56 0.13 3.35 4.14
2.19 2.16 2.01 0.55 0.96 2.06 1.56 4.14 1.72 2.26 2.67 1.78 0.55 0.09 3.26 4.04
2.34 2.33 2.16 0.50 0.99 2.11 1.63 4.23 1.76 2.23 2.79 1.85 0.58 0.15 3.35 4.18
‘Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 14 1.
15
Table 3 AM1 (MEPderived and exact SCF) dipole moments (D) for anions calculated at the center of mass Anion
EPM ’
IEM b
Exact
OHNH,CH,Ophenolate CHOOH,NCONHHCO? CNNONO1 OCNPH,HP@SCNHSCH,Ss20:ClOAlO-
0.97 1.39 1.39 3.70 1.70 3.12 2.62 0.81 0.87 0.79 3.05 2.64 0.29 0.06 1.63 3.32 4.78 3.25 1.94
0.94 1.36 1.28 3.55 1.60 3.00 2.52 0.84 0.88 0.75 2.96 2.54 0.27 0.03 1.58 3.22 4.60 3.07 1.85
0.97 1.39 1.38 3.70 1.69 3.14 2.63 0.84 0.89 0.73 3.05 2.62 0.29 0.13 1.61 3.28 4.76 3.29 2.00
’ Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 141.
AlO- and SCN-, respectively. However, considering all anions studied, the mean deviation remains quite small at 0.02 D, which is marginally larger than for neutral molecules and cations. Interestingly, the IEM gives a larger dipole than the EPM for CN- and NO-, whereas all other IEM values are underestimated. For anions the IEM dipole is occasionally the closer to the exact value, although the differences between IEM and EPM are reasonably small in such instances. All results for ions (Tables 2 and 3) were obtained with 0=1.5inEq. (7). In Fig. 1 we have plotted the ratio of MEP-derived to exact dipole as a log function of the vdW radii scaling factor CJ(Eq. (7)) for H20, FCl and CN-. All four MEP surfaces were used in the fitting procedure for each r~value. As the potential surfaces are moved further from the molecule (increasing a) the ratio of MEP to exact dipoles approaches unity, i.e. the exact dipole is approached asymptotically. As expected, the exact dipole is obtained for both methods in the limit as g tends to infinity, since formally the EPM is the limiting case of the IEM. However, the stability of the EPM over a wider range of scaling factors for neutral molecules is well demonstrated. In fact the
P.L. Cummins, J.E. Gready/ ChemicalPhysicsLetters 225 (1994) 11-17
16
Table 5 AM1 (MEPderived) and experimental quadrupole moments ( 10mz6esu cm*) calculated at the center of mass Molecule C2H2 (332 C2H4
benzene HCN NNO NH3 H,O xx d YY zz pyrrole XT YY zz pyridine xx YY zz
0:o
015
Table 4 Mean titted (i.e. cakulated from atomic charges) and exact (Eqs. (3) and (5) ) electrostatic potentials (eV) and corresponding root-mean-square (rms) deviations of the fitted potentials from the exact potentials as a function of the van der Waals radii scalingfactoru (Eq. (7)) Electrostatic potential (eV)
(fitted)
(exact)
llllS
difference
Hz0
CN-
2.5 1.5 1.0 0.5 2.5 1.5 1.0 0.5
0.0392 0.0878 0.1616 0.3461 1.4958 2.2179 2.8970 4.1593
0.0392 0.0876 0.1613 0.3263 1.4958 2.2183 2.8987 4.1584
Exp. =
4.04 5.70 1.22 -6.09 0.70 -5.29 - 1.98 - 1.32 -0.03 1.35 -7.10 5.10 2.00 -1.38 -4.42 5.80
3.67 -5.36 1.11 -5.59 0.61 - 5.06 -1.90 -1.25 -0.03 1.28 -6.51 4.67 1.84 -1.11 -4.28 5.39
3.0 -4.3 1.5 -5.6rt2.8 2.1 -3.0 -2.8 -2.4 0.3 2.7 -12.4k2.3 6.6k 1.2 5.8f 1.6 -6.2f 1.5 -3.5kO.9 9.7fl.l
-
1:o
Fig. 1. The dipole ratio, exact dipole/MEP-derived dipole, as a function of the log of the vdW radii scaling factor rr defined in Eq. (7) for (a) H20, (b) FCland (c) CN-. EMP=electrostatic potential method (-); IEM=interaction energy method (---) (see text).
0
IEM b
a Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 141. cSeeRef. [ 141 for source of experimental values. d X= axis perpendicular to molecular plane; z= dipole axis.
]og(o)
Species
EPM *
0.0006 0.0027 0.0083 0.0474 0.0067 0.0241 0.0624 0.2797
EPM reproduces the Hz0 dipole moment quite accurately even at radii close to half the vdW radius, i.e. the order of atomic or covalent radii. The effeo tive vdW radii of anions are larger than for cations or neutral molecules, hence the deviation tends to be larger for some anions at smaller 0 values. The mean electrostatic potentials calculated from the MEPderived (i.e. fitted) atomic charges and the exact values (Eqs. (3) and ( 5) ) and the root-mean-square (rms) deviations between fitted and exact potentials are given for several u values in Table 4. The rms deviations are small compared with the mean potentials for all values of a, and thus the atomic charges are able to reproduce the exact MEP quite close to the molecule with good precision. Although one of the strengths of the present method is that a set of atomic charges with a dipole very close to the exact MNDO, AM1 or PM3 value is virtually guaranteed, in some cases the exact values deviate significantly from experiment [ 91. This is particularly true for molecules containing the third row (Si, P, S and Cl) or heavier atoms, the cyanide (CN) functional group and atoms in unusual valence states (e.g. NNO). The agreement for molecules containing only first (i.e. hydrogen) and second row atoms
P.L. Cummins, J.E. Gready /Chemical Physics Letters 225 (1994) 11-I 7
bound according to their usual valence rules is, however, generally quite good. While dipole moments were used in the paramaterisation of MNDO, AM1 and PM3 [ 5-7 1, quadrupoles were not due to lack of accurate experimental data for large sets of molecules. Consequently, the agreement between AM 1 and observed quadrupoles is expected to be more variable, as is shown by the MEP-derived (cr= 1.5 ) and experimental quadrupole moments given in Table 5. As for dipole moments, the EPM quadrupole moments are generally larger than the IEM values.
17
[ 31 S.J. Weiner, P.A. Kollman, D.T. Nguyen and D.A. Case, J. Comput. Chem. 7 ( 1986) 230. [4] C.I. Bayly, P. Cieplak, W.D. Cornell and P.A. Kollman, J. Phys. Chem. 97 (1993) 10269. [S] M.J.S. Dewar and W. Thiel, J. Am. Chem. Sot. 99 (1977) 4899. [6] M.J.S. Dewar, E.G. Zoebisch, E.F. Healy and J.J.P. Stewart, J. Am. Chem. Sot. 107 (1985) 3902. [7] J.J.P. Stewart, J. Comput. Chem. 10 (1989) 209. [8] W. Thiel, Tetrahedron 44 (1988) 7393. [ 91 J.J.P. Stewart, J. Comput. Aided Mol. Design 4 ( 1990) 1. [lO]M.J.S.DcwarandW.Thiel,Theoret.Chim.Acta46 (1977)
[ 111 g. Luque and M. Grozco, Chem. Phys. Letters 168 ( 1990) 269.
4. Conclusion
[ 121 B.H. Besler, K.M. Men Jr. and P.A. Kollman, J. Comput. Chem. 11 (1990) 431.
[ 131 G.G. Ferenczy, C.A. Reynolds and W.G. Richards, J.
In this Letter we have described an efficient method for calculating the molecular electrostatic potential (MEP) and MEP-derived atomic charges within the semiempirical MNDO, AM1 and PM3 approximations. We tested the method for neutral molecules and ions, and found excellent agreement between the dipole moments calculated from the MEP-derived atomic charges and the corresponding exact SCF dipole moments. As the semiempirical dipole moments are generally quite accurate at least for molecules containing first and second row elements, the method should be suited for the calculation of atomic charges in simulations of large bioorganic molecules or as an alternative when the cost of density functional or ab initio calculations is prohibitive.
Acknowledgement This work was financed by the National Health and Medical Research Council.
References [ 1] P. Politzer and J.S. Murray, Rev. Comput. Chem. 2 ( 1991) 273. [2] U.C. Singh and P.A. Kollman, J. Comput. Chem. 5 (1984) 129.
Comput. Chem. 11 ( 1990) 159.
[ 141 P.L. Cummins and J.E. Gready, Chem. Phys. Letters 174 (1990) 355.
[ 15 ] C. Aleman, F.J. Luque and M. Orozco, J. Comput. Aided Mol. Design 7 (1993) 721.
[ 161 B. Wang and G.P. Ford, J. Comput. Chem. 15 ( 1994) 200. [ 171 M.J. Field, P.A. Bash and M. Karplus, J. Comput. Chem. 11 (1990) 700. [ 181 V.V. Vasilyev, A.A. Bliznyuk and A.A. Voityuk, Intern. J. Quantum Chem. 44 (1992) 897. [ 191 C.J. Cramer and D.G. Truhlar, J. Comput. Aided Mol. Design6 (1992) 629. [20] B. Wang and G.P. Ford, J. Chem. Phys. 97 (1992) 4162. [ 211 G.E. Chudinov, D.V. Napolov and M.U. Basilevsky, Chem. Phys. 160 (1992) 41. [ 221 G. Rauhut, T. Clark and T. Steinke, J. Am. Chem. Sot. 115 (1993) 9174. [ 23 ] F.J. Luque, M.J. Negre and M. Orozco, J. Phys. Chem. 97 (1993) 4386. [ 241 J. Gao, F.J. Luque and M. Orozco, J. Chem. Phys. 98 (1993) 2975. [25] J.A. Pople, D.P. Santry and G.A. Segal, J. Chem. Phys. 43 (1965) S129. [ 261 J.A. Pople, D.L. Beveridge and PA. Dobosh, J. Chem. Phys. 47 (1967) 2026. [27] J.A. Pople and D.L. Beveridge, Approximate molecular orbital theory (McGraw-Hill, New York, 1970). [28]M.J.S. Dewar and G. Klopman, J. Am. Chem. Sot. 89 ( 1967) 3089. [29] M.J.S. Dewar and N.L. Sabelli, J. Phys. Chem. 66 (1962) 2310. [ 301 K. Ohno, Theoret. Chim. Acta 2 ( 1964) 2 19. [ 3 1] G. Klopman, J. Am. Chem. Sot. 87 ( 1965) 3300.
ELSEVIER
CHEMICAL PHYSICS
Chemical Physics Letters 225 (1994) 11-17
The electrostatic potential in the semiempirical molecular orbital approximation Peter L. Cummins *, Jill E. Gready Department ofBiochemistry, University of Sydney, Sydney, NSW 2006, Australia Received 15 March 1994; in final form 5 May 1994
AbStIYICt A method is described for calculating the molecular electrostatic potential (MEP) within the MNDQ, AM1 or PM3 semiempirical approximations. The MEPs are used to obtain atomic charges from which dipole moments are calculated for comparison with the exact semiempirical values. The results for a wide range of neutral molecules and ions show excellent agreement between the dipole moments calculated from MEPderived atomic charges and the exact semiempirical dipole moments. The neutral molecules and cations showed a mean deviation in the dipole of less than 0.0 1 D. The maximum errors ( z 0.04 D) were obtained for some molecules containing third row or heavier elements. By contrast, all previously reported methods based on tbe semiempirical approximations for generating MEP-charges in general yield significantly larger deviations from tbe exact dipoles. The maximum error for the anions tested was 0.07 D, with a mean deviation only slightly higher than for neutral species and cations. The correct limiting behaviour is obtained when the MEP surfaces are calculated at large distances, i.e. tbe exact dipole is approached asymptotically. Quadrupole moments of the MEPderived charge distributions are reported for some nonpolar and polar molecules.
1. Introduction The molecular electrostatic potential (MEP) is defined in terms of the interaction energy between a molecule and a unit positive charge. If this point charge is exterior to the molecule, i.e. in a region where the molecular electron density is relatively small, the MEP may be used to rationalise chemical reactivities [ 1 ] or to obtain the atomic charges used in molecular dynamics or Monte Carlo simulations of biomolecular systems [ 2-41. Quantum mechanics provides a rigorous definition of the electrostatic potential as the expectation value of the Coulomb operator. In the density functional and ab initio formalisms, the calculation of MEPs is straightforward from
the molecular orbitals expressed in terms of a suitable basis set. However, in the semiempirical methods (MNDO, AM 1 and PM3) [ 5-91, properties are not computed as integrals over basis functions, but simplified functional forms are used to approximate the kinetic energy, one-electron and two-electron integrals required to evaluate the total energy within a minimal-basis HF formalism [ 5, lo]. As the basis is not strictly defined for the calculation of properties, it is not surprising that a number of diverse approaches have been proposed for the evaluation of the MEP and MEPderived atomic charges using semiempirical methods [ 11-l 61. In addition, several approaches for combining MNDO, AMI and PM3 with molecular mechanics (QM/MM) or continuum self-consistent reaction field (SCRF) meth-
* Correspondingauthor. E-mail: [email protected] 0009-2614/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIOOO9-2614(94)00617-Y
P.L. Cummins, J.E. Gready/Chemical PhysicsLetters225 (1994) I l-l 7
12
ods have also been proposed in the effort to model condensed phases [ 17-241. In this Letter, we have attempted to define the molecular electrostatic potential in a manner consistent with all approximations involved in the semiempirical method. To this end, we derive simplified functional forms for the electrostatic potential based on the core-electron attraction terms used in MNDO, AM1 and PM3. The dipole moments calculated from the MEP-derived atomic charge distributions are compared with the exact dipole moments as a test of the method, and also with other methods [ 1l-141 for obtaining MEP-derived charges.
2. Method The neglect of diatomic differential overlap (NDDO) approximation [ 25-271 on which AMl, MNDO, and PM3 are based requires that all oneelectron terms arising from diatomic overlap be neglected. Consequently, there are no three-centre contributions from atomic orbital products Q& where ,u and Yare centred on different atoms, as there would if the electrostatic potential were computed from a suitable finite-basis set expansion [ 11,15 1. Thus in the NDDO approximation, the MEP can be expressed solely in terms of atom-centred contributions (all equations given in atomic units), vi = C PpvINv + Zil& w
3
(1)
where Zi is the core charge of atom i, R, the distance from atom i, Ppy are elements of the first order density matrix, and Ipy are the one-electron integrals of the Coulomb operator involving only those orbitals fiy and @,,centred on the ith atom, i.e. Zpv=
s
I&-R’1 -‘$v(R’)@p(R’) CN’>
(2)
where Riis the position vector. All the semiempirical methods (MNDO, AM 1 and PM3 ) make use of twoelectron repulsion integrals of the form (ss, pv) [ 5, lo] in order to approximate one-electron integrals such as those given by Eq. (2 ). Since this approximation underestimates the core-electron attractions, the core-core repulsions must likewise be reduced from their Coulomb values in order to main-
tain the balance between the attractions and repulsions in the molecule. Purely empirical functions of the interatomic distances are then added to the corecore repulsion terms in order to tine-tune this balance [ 561: these have no explicit role in defining the electrostatic potential. In terms of the two-electron integrals, Eq. ( 1) may be cast in the form
where the subscript m indicates the transformation of (ss, ,uv) from the molecular Cartesian to local (diatomic) coordinate system. The total MEP at a point is simply given by the sum of V, over all atoms. However, since (ss, VP) underestimates core-electron attractions they are unsuitable for calculation of the MEP. In our interaction energy method [ 141, the Coulomb interaction terms in the fitting procedure used to obtain atomic charges were replaced by the semiempirical equivalent, i.e. (ss, ss), in order to compensate approximately for this error. While this method lead to atomic charges and molecular dipole moments no worse than those obtained by Ferency et al. [ 131, its deficiency was that an electrostatic potential could not be rigorously defined. However, simple analytic forms for the molecular electrostatic potential may be obtained by considering the approximations involved in the evaluation of the two-centre electron repulsion integrals (ss, ,uv) in MNDO, AM1 and PM3. Classically, the two-centre electron repulsion integrals (pv, na) are equivalent to the interaction energy between the continuous charge distributions @,&,and e&p In semiempirical treatments, these charge densities are replaced by sets of discrete charges [ 5,10 1. At large internuclear separations, i.e. large compared with the dimensions of the point-charge arrangement used to represent the atomic orbital products, the interactions between multipoles are accurately given by the sum of Coulomb interactions between charges. Since these interaction terms may diverge as internuclear separations are reduced, a modified Coulomb interaction is introduced which gives the correct one-centre limit, i.e. the semiempirical one-centre two-electron integrals. In MNDO, AM 1 and PM3, the Dewar-Sabelli, Ohno, Klopman (DSOK) equation [ 28-3 1 ] is used to calculate the interaction energies between charges that represent the multipoles of @,+, and &&,. If R is the
P.L. Cummins, J.E. Gready /Chemical Physics Letters 225 (1994) 1l-l 7
distance between any pair of these charges (one associated with nucleus X, the other with Y ) , the DSOK equation is given by a screened Coulomb potential of the form f(R)=
(Rz+Cuz) -1’2,
WI
a=px+py .
(4b)
Thus p are parameters chosen to reproduce the onecentre repulsion integrals (for X =Y) at zero internuclear separation. Consequently, a simplified analytic form for the Coulomb potential may be recovered by setting a! =px =pr = 0 in the two-electron integrals (ss, ,uv). Luque et al. [ 23 ] have tested both cw=O and cr=px in a self-consistent reaction field (SCRF) formalism for modelling solvent effects and found that the former choice gave heats of formation of various solutes in better agreement with experiment, despite their assertion that the latter is the more rigorous. Chudinov et al. [ 2 1 ] point out correctly that cr=O yields a’true’ electrostatic potential in their SCRF treatment. This choice was also used by Wang and Ford [ 201 in their solvation studies, but strangely not in their recent MEP calculations where a more complicated treatment is given [ 16 1. Setting a=0 in (ss, pv), we obtain the following expressions for the one-electron integrals for ( p = px or p,., z = pZorbitals ) : (SS,SS)=l/Ri,
(W
(SS,SZ)=f[(Ri-DI)-'-(Ri+DI)-'1
,
(5b)
(SSyZZ)=f{f[(Ri+2D2)-’ + (Ri_2D2)-‘]
- (SS, SS)} 3
(SS,pp)=f[(Rf+4D~)-“2-(ss,ss)]
(SC)
(5d)
where the parameters D, and D2 are determined by the condition that the dipole and quadrupole moments of each atomic orbital product q!~& are equal to those for the corresponding point-charge distributions. Values for D, and D2 are readily calculated using the formulae [ 10 ] D1 = (2”+ 1) (4C”&)“+1’2 ’ 3’9 &#,+ &)zn+2
(ha)
D~=[&,(2n+1)(n+1)]“2(~,)-‘,
(6b)
where c,,, and & are respectively the exponents of
13
the Slater s and p atomic orbitals for principal quantum number n. It is emphasized that Eqs. (3) and ( 5) would yield the correct Coulomb potential only in regions where the molecule’s electron density is small enough to neglect, i.e. at a point outside the molecule, since then RiB D. For the first and second row atoms, for example, the D parameter values for AM1 range from 0.2 to 0.5 A, compared with typical van der Waals (vdW) radii ranging from 1.2 to 1.7 A. The extent of the electron density will depend on the magnitude of the density matrix elements P#,,. Note also that in the limit of Ri+a, a in Eq. (4) may be neglected and the interaction energy method [ 141 would correspond to the present electrostatic potential method (EPM). In the present work, we have used Eq. ( 5 ) to compute the molecular electrostatic potentials for a wide range of neutral and charged molecules, and derive sets of atom-centred charges by a least-squares fitting procedure. This fitting procedure has been widely used to obtain point-charge representations of molecules and is described in detail elsewhere [ 2,12,14]. Briefly, the electrostatic potential is computed at a large number of points on a number of surfaces surrounding the molecule. A set of atomic charges are then derived which reproduce these potentials in a least-squares sense. Using a number of surfaces ensures a good fit both close to and far away from the molecule. In the present calculations, the atomic radii ( Rk) used in the calculation of the kth surface are givenby (inA) & = c&w +0.3+
1) )
(7)
where r~is a scaling factor, RvdW is the vdW radius of an atom and k= 1,2, 3 and 4, i.e. four surfaces were generated and used in the fitting procedure for obtaining atomic charges. For elements of unknown vdW radii, we have used twice the covalent or atomic radii in the present calculations. Most calculations were performed using a value of Q= 1.5 for scaling the vdW radii in order to be reasonably consistent the other MEP-derived charge calculations [ 11- 14 1. In some cases, we varied G in order to ascertain how well the method performs for surfaces very close to the molecule. A surface density of at least two points per A2 was chosen to give dipole moments to within a precision of f 0.01 D. All calculations were per-
P.L. Cummins, J.E. Gready/Chemical PhysicsLetters225 (1994) II-I 7
14
formed at the AM 1 optimised geometries of the neutral molecules and ions.
3. Results and discussion Here we present AM1 results only, as we reached similar conclusions with test calculations for the other NDDO-based semiempirical methods (MNDO and PM3). In Table 1, the dipole moments calculated from the MEPderived point-charge representation using the electrostatic potential method EPM (i.e. Eqs. ( 3 ) and ( 5 ) ) are compared with the exact AM 1 and experimental values for neutral molecules. Also included in Table 1 for comparison are the dipoles calculated from semiempirical MEP-derived atomic charges reported for the various other methods in the
literature [ 1l-l 41. We have recalculated MEPs for the interaction energy method (IEM) [ 141, since a minor error, which we have now corrected, was detected in the surface generation code. This error lead to dipole moments which differ from the ones given in Table 1 by Y/o- 10%. The MEPs obtained using the EPM and IEM methods were calculated with a value of cr= 1.5 in Eq. (7). The results in Table 1 indicate the IEM consistently predicts dipole moments smaller than the exact values. However, the agreement between the EPM and exact molecular dipole moments is excellent. The mean absolute difference between EPM and exact dipoles for neutral molecules is less than 0.0 1 D, with a maximum deviation not exceeding 0.04 D. On average, the deviations tend to be larger for molecules containing third-row or heavier elements. In contrast, the other methods often yield
Table 1 AM 1 ( MEPderived and exact SCF) and experimental dipole moments (D ) for neutral molecules Molecule
EPM .
IEM ”
HF Hz0 NH, CHsO CH,OH NH&HO CHOOH CH3CN CH,CO HCl FCl
1.73 1.86 1.84 2.33 1.63 3.70 1.48 2.89 1.35 1.40 0.87 1.50 0.08 0.64 4.37 2.36 3.60 1.95 2.05 0.49 1.75 1.88 2.30 5.64 0.61 2.16 3.17 5.30
1.64 1.76 1.76 2.24 1.58 3.70 1.42 2.78 1.39 1.35 0.84 1.45 0.08 0.61 4.24 2.27 3.48 1.88 1.99 0.48 1.69 1.82 2.18 5.54 0.60 2.14 3.13 5.21
CH3Wz
NO NNO ClNO HCN imidaxole pyrrole pyridine furan HsSiC& HsS PH3
ZnO BrI CH3GeN CH,HgBr HSO
Ref. [ 121 1.28 1.36 1.93 1.01 2.75 1.39 2.61
’ Electrostatic potential method (Eqs. (3) and (5) in text). cSeeRef. [ 141 for sources of experimental dipole values.
Ref. [ 131 1.75 1.89 1.75 2.32 1.62 3.68 1.46 2.84 1.38 1.45
Ref. [ 161 1.82 1.74 2.29 1.62 3.60 1.42 2.82
b Interaction energy method, Ref. [ 141.
Exact
Exp. c
1.74 1.86 1.85 2.32 1.62 3.70 1.48 2.89 1.35 1.38 0.91 1.50 0.09 0.64 4.36 2.36 3.60 1.95 2.04 0.49 1.76 1.86 2.29 5.66 0.63 2.19 3.15 5.32
1.82 1.85 1.47 2.33 1.70 3.73 1.41 3.92 1.42 1.11 0.89 1.31 0.16 0.17 2.98 3.8 1.74 2.19 0.66 1.18 0.97 0.58 0.72 3.99 3.10
P.L.. Cummins, J.E. Gready / Chemical Physics Letters 225 (1994) 1l-1 7
significantly larger deviations from the exact dipoles. The deorthogonalization (finite basis set) approach [ 11,15 1, also adopted by Besler et al. [ 12 ] for the calculation of atomic charges, gives by far the largest deviation of 0.95 D for NH&HO, while the IEM [ 141 and Ferenczy et al. [ 13 ] methods give maximum deviations of 0.17 and 0.10 D, respectively. The recently published results of Wang and Ford [ 16 ] also yield a maximum deviation of 0.10 D. It is worth noting that the numbers of molecules tested in these other studies [ 12,13,16] are quite small and the maximum errors may well be larger, particularly for molecules containing heavier elements. The dipole moments calculated using the EPM and IEM at the centre of mass for cations are given in Table 2 along with the exact SCF values. A deviation of 0.02 D is obtained for CH3Zn+ and OPCl$, while the largest deviation of 0.04 D is obtained for ClNO+ and HgCl+. Note that the IEM overestimates the dipole for HCNH+ and FCl*+ whereas all other cations give IEM dipoles which are less than the exact values. The average deviation for cations is much the same as for neutral molecules, i.e. not more than 0.01 D. The EPM, IEM and exact dipoles calculated at the centre of mass for anions are given in Table 3. The largest deviations of 0.06 and 0.07 Dare obtained for Table 2 AM1 (MePderived and exact SCF) dipole moments (D) for cations calculated at the center of mass Cation
EPM ’
IEM b
Exact
H,O+ CH,NHf CHsOH: HCNH+ NCO+ HCl+ imidazole*H+ fluorobenzene.H+ FCl+ FCl’+ ClNO+ H,S+ 0PCl: CH,Zn+ ZnOH+ H&l+
2.33 2.33 2.16 0.50 1.00 2.13 1.63 4.23 1.72 2.22 2.83 1.85 0.56 0.13 3.35 4.14
2.19 2.16 2.01 0.55 0.96 2.06 1.56 4.14 1.72 2.26 2.67 1.78 0.55 0.09 3.26 4.04
2.34 2.33 2.16 0.50 0.99 2.11 1.63 4.23 1.76 2.23 2.79 1.85 0.58 0.15 3.35 4.18
‘Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 14 1.
15
Table 3 AM1 (MEPderived and exact SCF) dipole moments (D) for anions calculated at the center of mass Anion
EPM ’
IEM b
Exact
OHNH,CH,Ophenolate CHOOH,NCONHHCO? CNNONO1 OCNPH,HP@SCNHSCH,Ss20:ClOAlO-
0.97 1.39 1.39 3.70 1.70 3.12 2.62 0.81 0.87 0.79 3.05 2.64 0.29 0.06 1.63 3.32 4.78 3.25 1.94
0.94 1.36 1.28 3.55 1.60 3.00 2.52 0.84 0.88 0.75 2.96 2.54 0.27 0.03 1.58 3.22 4.60 3.07 1.85
0.97 1.39 1.38 3.70 1.69 3.14 2.63 0.84 0.89 0.73 3.05 2.62 0.29 0.13 1.61 3.28 4.76 3.29 2.00
’ Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 141.
AlO- and SCN-, respectively. However, considering all anions studied, the mean deviation remains quite small at 0.02 D, which is marginally larger than for neutral molecules and cations. Interestingly, the IEM gives a larger dipole than the EPM for CN- and NO-, whereas all other IEM values are underestimated. For anions the IEM dipole is occasionally the closer to the exact value, although the differences between IEM and EPM are reasonably small in such instances. All results for ions (Tables 2 and 3) were obtained with 0=1.5inEq. (7). In Fig. 1 we have plotted the ratio of MEP-derived to exact dipole as a log function of the vdW radii scaling factor CJ(Eq. (7)) for H20, FCl and CN-. All four MEP surfaces were used in the fitting procedure for each r~value. As the potential surfaces are moved further from the molecule (increasing a) the ratio of MEP to exact dipoles approaches unity, i.e. the exact dipole is approached asymptotically. As expected, the exact dipole is obtained for both methods in the limit as g tends to infinity, since formally the EPM is the limiting case of the IEM. However, the stability of the EPM over a wider range of scaling factors for neutral molecules is well demonstrated. In fact the
P.L. Cummins, J.E. Gready/ ChemicalPhysicsLetters 225 (1994) 11-17
16
Table 5 AM1 (MEPderived) and experimental quadrupole moments ( 10mz6esu cm*) calculated at the center of mass Molecule C2H2 (332 C2H4
benzene HCN NNO NH3 H,O xx d YY zz pyrrole XT YY zz pyridine xx YY zz
0:o
015
Table 4 Mean titted (i.e. cakulated from atomic charges) and exact (Eqs. (3) and (5) ) electrostatic potentials (eV) and corresponding root-mean-square (rms) deviations of the fitted potentials from the exact potentials as a function of the van der Waals radii scalingfactoru (Eq. (7)) Electrostatic potential (eV)
(fitted)
(exact)
llllS
difference
Hz0
CN-
2.5 1.5 1.0 0.5 2.5 1.5 1.0 0.5
0.0392 0.0878 0.1616 0.3461 1.4958 2.2179 2.8970 4.1593
0.0392 0.0876 0.1613 0.3263 1.4958 2.2183 2.8987 4.1584
Exp. =
4.04 5.70 1.22 -6.09 0.70 -5.29 - 1.98 - 1.32 -0.03 1.35 -7.10 5.10 2.00 -1.38 -4.42 5.80
3.67 -5.36 1.11 -5.59 0.61 - 5.06 -1.90 -1.25 -0.03 1.28 -6.51 4.67 1.84 -1.11 -4.28 5.39
3.0 -4.3 1.5 -5.6rt2.8 2.1 -3.0 -2.8 -2.4 0.3 2.7 -12.4k2.3 6.6k 1.2 5.8f 1.6 -6.2f 1.5 -3.5kO.9 9.7fl.l
-
1:o
Fig. 1. The dipole ratio, exact dipole/MEP-derived dipole, as a function of the log of the vdW radii scaling factor rr defined in Eq. (7) for (a) H20, (b) FCland (c) CN-. EMP=electrostatic potential method (-); IEM=interaction energy method (---) (see text).
0
IEM b
a Electrostatic potential method (Eqs. (3) and (5) in text). b Interaction energy method, Ref. [ 141. cSeeRef. [ 141 for source of experimental values. d X= axis perpendicular to molecular plane; z= dipole axis.
]og(o)
Species
EPM *
0.0006 0.0027 0.0083 0.0474 0.0067 0.0241 0.0624 0.2797
EPM reproduces the Hz0 dipole moment quite accurately even at radii close to half the vdW radius, i.e. the order of atomic or covalent radii. The effeo tive vdW radii of anions are larger than for cations or neutral molecules, hence the deviation tends to be larger for some anions at smaller 0 values. The mean electrostatic potentials calculated from the MEPderived (i.e. fitted) atomic charges and the exact values (Eqs. (3) and ( 5) ) and the root-mean-square (rms) deviations between fitted and exact potentials are given for several u values in Table 4. The rms deviations are small compared with the mean potentials for all values of a, and thus the atomic charges are able to reproduce the exact MEP quite close to the molecule with good precision. Although one of the strengths of the present method is that a set of atomic charges with a dipole very close to the exact MNDO, AM1 or PM3 value is virtually guaranteed, in some cases the exact values deviate significantly from experiment [ 91. This is particularly true for molecules containing the third row (Si, P, S and Cl) or heavier atoms, the cyanide (CN) functional group and atoms in unusual valence states (e.g. NNO). The agreement for molecules containing only first (i.e. hydrogen) and second row atoms
P.L. Cummins, J.E. Gready /Chemical Physics Letters 225 (1994) 11-I 7
bound according to their usual valence rules is, however, generally quite good. While dipole moments were used in the paramaterisation of MNDO, AM1 and PM3 [ 5-7 1, quadrupoles were not due to lack of accurate experimental data for large sets of molecules. Consequently, the agreement between AM 1 and observed quadrupoles is expected to be more variable, as is shown by the MEP-derived (cr= 1.5 ) and experimental quadrupole moments given in Table 5. As for dipole moments, the EPM quadrupole moments are generally larger than the IEM values.
17
[ 31 S.J. Weiner, P.A. Kollman, D.T. Nguyen and D.A. Case, J. Comput. Chem. 7 ( 1986) 230. [4] C.I. Bayly, P. Cieplak, W.D. Cornell and P.A. Kollman, J. Phys. Chem. 97 (1993) 10269. [S] M.J.S. Dewar and W. Thiel, J. Am. Chem. Sot. 99 (1977) 4899. [6] M.J.S. Dewar, E.G. Zoebisch, E.F. Healy and J.J.P. Stewart, J. Am. Chem. Sot. 107 (1985) 3902. [7] J.J.P. Stewart, J. Comput. Chem. 10 (1989) 209. [8] W. Thiel, Tetrahedron 44 (1988) 7393. [ 91 J.J.P. Stewart, J. Comput. Aided Mol. Design 4 ( 1990) 1. [lO]M.J.S.DcwarandW.Thiel,Theoret.Chim.Acta46 (1977)
[ 111 g. Luque and M. Grozco, Chem. Phys. Letters 168 ( 1990) 269.
4. Conclusion
[ 121 B.H. Besler, K.M. Men Jr. and P.A. Kollman, J. Comput. Chem. 11 (1990) 431.
[ 131 G.G. Ferenczy, C.A. Reynolds and W.G. Richards, J.
In this Letter we have described an efficient method for calculating the molecular electrostatic potential (MEP) and MEP-derived atomic charges within the semiempirical MNDO, AM1 and PM3 approximations. We tested the method for neutral molecules and ions, and found excellent agreement between the dipole moments calculated from the MEP-derived atomic charges and the corresponding exact SCF dipole moments. As the semiempirical dipole moments are generally quite accurate at least for molecules containing first and second row elements, the method should be suited for the calculation of atomic charges in simulations of large bioorganic molecules or as an alternative when the cost of density functional or ab initio calculations is prohibitive.
Acknowledgement This work was financed by the National Health and Medical Research Council.
References [ 1] P. Politzer and J.S. Murray, Rev. Comput. Chem. 2 ( 1991) 273. [2] U.C. Singh and P.A. Kollman, J. Comput. Chem. 5 (1984) 129.
Comput. Chem. 11 ( 1990) 159.
[ 141 P.L. Cummins and J.E. Gready, Chem. Phys. Letters 174 (1990) 355.
[ 15 ] C. Aleman, F.J. Luque and M. Orozco, J. Comput. Aided Mol. Design 7 (1993) 721.
[ 161 B. Wang and G.P. Ford, J. Comput. Chem. 15 ( 1994) 200. [ 171 M.J. Field, P.A. Bash and M. Karplus, J. Comput. Chem. 11 (1990) 700. [ 181 V.V. Vasilyev, A.A. Bliznyuk and A.A. Voityuk, Intern. J. Quantum Chem. 44 (1992) 897. [ 191 C.J. Cramer and D.G. Truhlar, J. Comput. Aided Mol. Design6 (1992) 629. [20] B. Wang and G.P. Ford, J. Chem. Phys. 97 (1992) 4162. [ 211 G.E. Chudinov, D.V. Napolov and M.U. Basilevsky, Chem. Phys. 160 (1992) 41. [ 221 G. Rauhut, T. Clark and T. Steinke, J. Am. Chem. Sot. 115 (1993) 9174. [ 23 ] F.J. Luque, M.J. Negre and M. Orozco, J. Phys. Chem. 97 (1993) 4386. [ 241 J. Gao, F.J. Luque and M. Orozco, J. Chem. Phys. 98 (1993) 2975. [25] J.A. Pople, D.P. Santry and G.A. Segal, J. Chem. Phys. 43 (1965) S129. [ 261 J.A. Pople, D.L. Beveridge and PA. Dobosh, J. Chem. Phys. 47 (1967) 2026. [27] J.A. Pople and D.L. Beveridge, Approximate molecular orbital theory (McGraw-Hill, New York, 1970). [28]M.J.S. Dewar and G. Klopman, J. Am. Chem. Sot. 89 ( 1967) 3089. [29] M.J.S. Dewar and N.L. Sabelli, J. Phys. Chem. 66 (1962) 2310. [ 301 K. Ohno, Theoret. Chim. Acta 2 ( 1964) 2 19. [ 3 1] G. Klopman, J. Am. Chem. Sot. 87 ( 1965) 3300.