4740
J . Phys. Chem. 1990, 94, 4140-4746
Determlnatlon of Net Atomic Charges Using a Modified Partial Equallzation of Orbltal Electronegativlty Method. 2. Application to Ionic and Aromatic Molecules as Models for Polypeptides Kyoung Tai No,+'* J. Andrew Grant,: Mu Shik Jhon,s and Harold A. Scheraga*J Baker Laboratory of Chemistry, Cornel1 University, Ithaca, New York 14853-1301, and Department of Chemistry, Korea Advanced Institute of Science and Technology, P.O. Box 150, Cheong Ryang Ri, Seoul, Korea (Received: October 2, 1989; In Final Form: January 4 , 1990)
To finalize a scheme intended to calculate directly approximate charge distributions for proteins, parameters for a modified PEOE method are obtained for both ionic and aromatic functional groups. For ionic groups, the electrostatic potentials of several ionic molecules (computed with the 6-3 1G** basis set) are used as constraints;for aromatic groups, the experimental gas-phase dipole and quadrupole moments of several aromatic molecules are used as constraints. The electrostatic potentials calculated with the modified PEOE point charges of ionic molecules are quite close to the ab initio electrostatic potential surfaces, and the dipole and quadrupole moments calculated with the modified PEOE charges of aromatic molecules agree well with experimental data. The method enables the charge distribution, described as a set of atom-centered monopoles, to be calculated directly. This avoids assumptions concerning transferability of charges from small model compounds and describes in an approximate fashion the redistributionof charge upon assembly of amino acid residues into polypeptide molecules. The computational requirements for the procedure are negligible.
Introduction In our preceding paper,l we developed methods for the direct calculation of net atomic charges of neutral molecules on the basis of the concept of partial equalization of orbital electronegativity (PEOE). The original formulation of the PEOE methodz4 describes the orbital electronegativity in terms of purely atomic properties, viz., ionization potentials and electron affinities. This is probably the reason that this prescription fails to reproduce experimentally measured dipole and quadrupole moments accurately. The inability of this method to describe such a fundamental feature of the molecular charge distribution suggests that it would be unsuitable for modeling the electrostatic component of the total inter/intramolecular potential energy function for biomacromolecules. To improve the description of the permanent multipole moments, the original method was modified to determine the parameters that describe the electronegativity function by utilizing information obtained from experimental dipole and quadrupole moments. This is similar in spirit to a suggestion first made by Abraham and Smith.$ The resulting calculated dipole and quadrupole moments of several types of neutral molecules, chosen to resemble the functional groups found in polypeptides, were found to be in good agreement with observed values.] In addition, however, polypeptides may contain several kinds of either ionic or aromatic side chains as well as the ionic terminal amino +H3N and carboxyl COT groups. It is the purpose of the present work to extend our modified PEOE methods to such groups to produce a comprehensive consistent scheme to determine approximate molecular charge distributions for proteins. It is perhaps prudent to reiterate that the use of atom-centered monopoles, although conceptually and computationally simple, is only a very approximate representation of molecular charge distributions. There is an increasing body of literature highlighting the importance of molecular anisotropy, although it remains unclear as to when the modeling of such effects will be practical for biomacromolecular systems. Central to the parametrization of orbital electronegativity in the preceding paper1 was the availability of experimentally determined dipole and quadrupole moments, as well as accurately measured molecular geometries. For ionic species, such reliable
gas-phase data are not available. Approximate gas-phase dipole moments estimated from properties of ionic solution^^*^ are not considered to be of sufficient accuracy. Furthermore, it remains difficult to determine precisely the geometry of the solution species, and even then the problem remains of how this is related to the gas-phase geometry. In the present work, therefore, parameters describing orbital electronegativity for ionic species are obtained from ab initio electrostatic potentials. For the aromatic molecules, parametrization of orbital electronegativity is carried out in terms of experimental dipole and quadrupole moments, as in the previous paper on neutral mo1ecules.l Combined with the results for neutral molecules, this finalizes a scheme for the rapid and direct computation of molecular charge distributions of peptides/proteins. The directness of the approach circumvents assumptions about the transferability of charges from a small model compound to a large biomolecule. The method accounts, at least in an approximate way, for the redistribution of charge between molecular fragments, particularIy important in the vicinity of ionic groups. The necessity to employ ab initio calculations in parametrizing orbital electronegativities for ionic groups lacks the elegance and clarity characterizing the procedure adopted for neutral molecules. It is also the source of some uncertainty, because of the imponderable nature of the accuracy achieved in the ab initio calculations. A full discussion of the expected reliability and problems associated with the calculation of net atomic charges for ionic species is presented in the Results and Discussion section.
Methods Determination of the PEOE Parameters of the Atoms in Ionic Groups. In the modified PEOE methods described previously,' the net atomic charge, Qi, of the ith atom is a function of the following independent parameters: the orbital electronegativity constants (a,b),a set of damping factors v), and a single scaling factor xH+. For reasons already outlined, it is necessary to de( I ) No, K. T.; Grant, J. A.; Scheraga, H. A. J . Phys. Chem., preceding paper in this issue. (2) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (3) Morticr, W. J.; Genechten, K. V.;Gasteiger, J. J. Am. Chem. SOC. 1985, 107, 829.
(4) Hammarstrom, L. G.; Liljefors, T.; Gasteiger, J. J. Compuf.Chem. 1988, 9, 424.
*Towhom requests for reprints should be addressed. 'On leave from the Department of Chemistry, Soong Si1 University, Sang Do 1-1, Dong Jak Gu. Seoul, Korea. 'Cornell University. 'Korean Advanced Institute of Science and Technology. 0022-3654/90/2094-4140$02.50/0
( 5 ) Abraham, R. J.; Smith, P. E. J . Compur. Chem. 1987, 9, 288. (6) (a) Debye, P. Polar Molecules; The Chemical Catalog Co.: New York, 1929. (b) Guggenheim, E. A. Trans. Faraday SOC.1949, 45, 714. (7) (a) BBttcher, C. J. F. Theory ojPolarization; Elsevier: Amsterdam, 1973; Vol. 6. (b) Myers, A. B.; Birge, R. R. J . Chem. Phys. 1981, 74, 3514.
(c) Makosz, J . J. J. Chem. Phys. 1987, 87, 6053.
0 1990 American Chemical Society
Determination of Net Atomic Charges. 2 TABLE I: Dependence of PW Charges and Dipole Moments on Basis Set (Charges a d Dipole Moments in e and D, Respectively) STO-3G 6-31G 6-31G* 6-31G** PD PD PD PD 6-31G**' N -0.678 -0.818 -d -0,845 0.462 -d 0.427 0.471 0.000 0.000 0.000 0.000 0.000 0.120 0.274 0.308 0.249 0.040 0.053 0.036 0.025 -0.3 13 -0.068 -0.347 -0.313 0.264 0.324 0.311 0.314 2.185 2.164 2.205 2.529 2.190 -0.244 -0.232 -0.204 -0.205 -0.713 -0.912 -0.824 -0.825 0.855 H 0.870 0.686 1.095 1.460 0.577 1.673 1.343 1.316 DMb -0.444 -0.334 0.032 -0.119 CH$OzC H 0.057 0.038 -0.061 -0.020 0.734 1.144 0.867 0.909 C(0) -0.862 0 -0.729 -0.960 -0.858 3.425 DM' 2.466 3.877 3.429 3.823 glycine C" 0.085 0.053 0.043 0.065 0.066 0.044 0.051 zwitterion H(Ca) 0.031 0.819 C' 0.635 1.020 0.813 -0.763 -0.660 -0.865 -0.763 0 -0.400 -0.503 -0,379 -0.404 N 0.314 0.325 0.346 0.306 H(N) DMb 12.196 11.859 12.813 12.223 12.225 0.283 0.149 0.236 0.245 alanine C' 0.07 1 0.021 0.113 0.073 zwitterion H(Cu) C5 -0.359 -0.317 -0.349 -0.359 0.118 0.106 0.115 0.115 H(Ce) 0.883 0.672 1.158 C' 0.883 -0.772 0 -0.669 -0.890 -0.772 -1.199 N -1.033 -1.374 -1.171 0.483 0.572 0.509 0.516 H(N) 11.561 DMb 1 I .202 12.008 1 1.470 1 1.468 PD: potential-derived. Dipole moment calculated from point charges of each basis set. eDipole moment calculated directly from ab initio wave function for each molecule. dSince NH4+ has high symmetry, the PD charges are sometimes quite unrealistic. The net atomic charge of the H atom is 1.74.
termine this parameter set for the ionic molecules by using information obtained from a b initio calculations. Closed shell restricted HartreeFock (RHF) calculations were therefore carried out for the molecules NH4+, CH3NH3+, HCO;, CH3C02-, +H3NCH2CO2-(glycine zwitterion), and +H3NCH(CH3)C02(alanine zwitterion). The basis sets used for the calculations are given in Table I. Having computed the R H F wave function, it is a fairly straightforward matter to calculate the electrostatic potential at a point in space by
where r, ZA,RA, R,,, and 4 represent a point around a molecule, the nuclear charge of atom A, the position of nucleus A, the density matrix element, and the atomic orbital, respectively. Similarly, the electrostatic field (the derivative of the potential) can be obtained, with little further computational effort. All the calcuiations were performed using the GAUSSIAN 86 computer program.I0 The number of points (N,) used for the calculation of the potential or the field around the molecule was 700-1000, depending upon the molecular size. The points were chosen to lie in a spherical shell between the van der Waals surface of the molecule and a distance up to 3 A from this surface. To avoid choosing symmetrically identical points, the origin of the coordinate system was taken at (x, + d/4, y , + 2d/4, z, + 3d/4), where (x,, (8) Cox, S. R.; Williams, D. E. J . Comput. Chem. 1981, 2, 304. (9) (a) Kollman, P. A. J . Am. Chem. SOC.1978, 100, 2974. (b) Singh, U. C.; Kollman, P. A. J . Comput. Chem. 1984, 5, 129. (10) Frisch, M. J.; Binkley, J. S.; Schlegel, H. B.; Raghavachari, K.; Melius, C. F.; Martin, R. L.;Stewart, J. J. P.; Bobrowicz, F. W.; Rohlfing, C. M.; Kahn, L. R.; Defrees, D. J.; Seeger, R.; Whiteside, R. A,; Fox, D. J.; Fleuder, E. M.; Pople, J. A. GAUSSIAN 86; Carnegie-Mellon Quantum Chemistry Publishing Unit: Pittsburgh, PA, 1984.
The Journal of Physical Chemistry, Vol. 94, No. 11, 1990 4741 TABLE 11: Potential-Derived (PD) and Field-Derived (FD) Net Atomic Charges for Formamide Obtained at Seven1 Shell Thickness by Using the 6-31G** Basis Set (Charges in e and Dipole Moments in D) ( d , - 4):A 0- 1 0-2 0-3 1-2 2-3 Potential-Derived Net Atomic Charge (PD Charge) C 0.732 0.728 0.728 0.717 0.732 0 -0.598 -0.588 -0.588 -0.586 -0.590 N -1.010 -1.003 -1.002 -0.983 -1.005 -0.007 -0.007 -0.007 -0.005 -0.010 H(C) H(N)b 0.442 0.440 0.439 0.433 0.440 DMc 4.268 4.265 4.264 4.260 4.260 Field-Derived C 0.764 0 -0.599 N -1,027 H(C) -0.008 0.444 H(N)b DMC 4.270
Net Atomic Charge (FD Charge) 0.758 0.758 0.728 0.322 -0.589 -0.484 -0.597 -0.597 -0.993 -0.495 -1.021 -1.024 -0.002 0.093 -0.007 -0.008 0.443 0.444 0.435 0.283 4.267 4.269 4.265 4.251
In the calculation of the point charges, the electrostatic potentials or fields were obtained for points located between d , and d2. d , and d2 represent the shortest distances from the van der Waals surface of a molecule. bCis and trans H atoms were not distinguished. CDipole moment calculated from the PD and FD point charges. The dipole moment calculated directly from the wave function is 4.251 D.
yc, zc) are the coordinates of the center of mass, and d is the grid interval. A set of potential-derived (PD) net atomic charges was subsequently obtained in the conventional manner,8s9namely, by minimizing the function
where I/lfab(ri)and I/lfpD(ri)are the electrostatic potentials at point ri of the kth molecule calculated from the ab initio wave function and the PD point charges, respectively, where
(3) NAk,Qj,k,e, and Rj,kare the number of atoms of the kth molecule,
the net atomic charge of the jth atom of the krh molecule, the electronic charge, and the position of the jth atom in the kth molecule, respectively. A set of field-derived charges (FD) can be obtained in an identical manner, replacing values and expressions appropriate for the components of the electric field for the electrostatic potential in eqs 2 and 3 . The dependence of these charges on the choice of the outermost boundary defining the space sampled is summarized in Table 11. Likewise, variation of the PD charges with respect to the basis sets employed is given in Table I. The net atomic charge distribution determined in this manner for the chosen set of ionic molecules is a function of the set of modified PEOE parameters as described previously. The classification of atomic species and the damping factors used in the construction of this parameter set are detailed in Table 111. The optimum values of the elements of this parameter set are obtained simply by finding the minimum of the function F = CIFkl k
(4)
where the function Fk is described in eq 2, the superscript k indexing the molecules in the chosen set. Clearly the individual elements ai of the parameter space emerge when the following condition is satisfied: (8F/8ai),,;,+, = 0 for all i (5) As in the case of neutral molecules,' the electronegativity of the ith atom, xi("),in an ionic group, in the nth iteration, is described by a linear function of the net atomic charge, Qi("): xi(") = ai + biQi(") (6)
4742
No et al.
The Journal of Physical Chemistry, Vol. 94, No. I I , 1990
TABLE 111: Atomic Species and Damping Factors Used for Ionic CrouDs in the Modified PEOE Method Atomic Species C atom connected to Nt C(Nt) c(0-1 C atom in COC atom connected to both N t and C(0-), C? C" in single amino acid H atom connected to N t H(Nt) H atom connected to C(N+,C') H(C3 H(C02-)b H atom in H C 0 2 00 atom in C 0 2 N+ N t sp3 atom damDing. factor' (A-B)~ f6 Nt-H(N) h N+-Ci, Nt-C(Nt) fs 0--c(0-) h C( O-)-Ci, C ( O-)-C" fit C (O-)-H (C02-) ~
~~
'Appears only in zwitterions of single amino acids. The subscript m indicates that the C" atom is in a single (monomeric) amino acid. Does not appear in peptide molecules. 'Damping factors between C(N+), C ( 0 - ) , C;, and the atoms not included in an ionic group are the same as those in neutral molecules.' dDamping factors for the atomic pair A and B.
where ai and bi are the electronegativity parameters of the ith and ( n ) represents the nth iteration. The iterative nature of the PEOE procedure dictates that a set of initial charges be assigned to each atom in the molecule under consideration, to facilitate the computation of the initial electronegativity function, viz. xi(") = ai + biQ(0) (7) In the case of neutral molecules these charges {Q'O))were chosen to be zero. However, for the ionic species, it was necessary to assign finite initial charges to avoid large, sometimes oscillatory, changes in the electronegativity function, often preventing convergence of the procedure. The choice of these charges was found to be somewhat arbitrary; however, numerical values, given in Table IV for all of the modified PEOE parameters for ionic species, gave accurate final charge distributions on the basis of experience acquired in reproducing the a b initio dipole moments and electrostatic potential. The final net atomic charge Q(")was obtained by convergence within IO steps and was insensitive to small changes in Q(O). Larger changes in Q'O),e.g., 0.1, did yield significant changes in Q"). The results of these comparisons are presented in Table V. In carrying out these calculations, the PEOE parameters associated with atoms not comprising ionic functional groups (the methyl group in the alanine zwitterion, for example) were taken from the DP5 and DQP5 sets given in the preceding paper.' The optimized ab initio geometries were taken from ref 11 for the ionic molecules. The distances N-H, C-N, and C-H and the bond angles HNH, NHC, and H C N in the cationic groups are 1.01, 1.51, and 1 .OO A and 109.47, 111.6, and 109.47', respectively. The distances C-O, C-C, C-H (in HCO,), and C-H and the bond angles HCO (in HC02-), HCC, and OCC in the anionic groups are 1.231, 1.522, 1.127, and 1.00 A and 114.5, 109.47, and 114.5', respectively. For formamide, the distances N-H, C-0, N-C, and C-H and the bond angles HCO, OCN, CNH, and CNH,are 1.001, 1.219, 1.352, and 1.098 Aand 122.5, 124.7, 1 18.5, and 1 20°, respectively. The methylammonium ion is in the staggered form, and the dihedral angle HC-CO of CH3-COOis 60'. The bond angles at C" in zwitterionic glycine and alanine were assumed to be tetrahedral, and the distances C"-C@ and C@-Hwere taken as 1.526 and 1.09 A, respectively. The dihedral angles NC"-C'O and HN-C'C' in glycine were taken as -18.9 and 60.0°, respectively, and the dihedral angles HC@-C"C', C@CpC'O,HC"-C@H,HN-C"C8, and NC"-CBH of alanine were taken as 180, -104, -55.2, -174, and 62.3', respectively. Aromatic Molecules. In obtaining modified PEOE parameters for use in describing aromatic molecules, the technique established ( I I ) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2 I74
TABLE I V ODtimized PEOE Panmeters for Ionic CrouDs" DPb DQP atom a; b, ai b; ood C(W 8.846 6.994 8.660 6.893 0.35 (or 0.30, 0.25, 0.20) c(o-) 5.165 2.894 5.159 3.005 0.20 c: 7.796 1.577 7.772 2.008 0.3Y (or 0.30, 0.25) 7.443 8.459 7.067 8.445 0.35 H(N+) H(C(Nt)l. 7.774 19.097 7.963 19.067 0.00 (or 0.051 L
8.895 9.937 9.024 9.962 0.10 14.510 9.339 14.664 9.324 -0.60 14.397 14.658 15.722 14.277 -0.40 0.488 0.467 0.679 0.703 0.485 0.466 0,704 0.683 0.786 0.805 "The parameters ai and bi are those of eq 6. bDP5 parameters were used for the atoms not included in the ionic group. 'DQP5 parameters were used for the atoms not included in the ionic group. dInitial charges in the PEOE calculation. The initial net atomic charge of the H atoms bonded to C(N+) or C can be either 0.00 or 0.05; 0.05 gave a slightly better result than did 0.00. 'The sum of the initial net atomic charges of C i (or C(N+)) and the H atoms bonded to it, H(C(N+)) and H(CZ)),must be 0.35. The values in parentheses are other possible initial net atomic charges of C: (or C(Nf)).
for the neutral molecules] was employed, with only a minor modification. That is, instead of using only the total molecular dipole moment as part of the fitting function, sometimes the components of the dipole moment were used, these being experimentally available for some in the set of molecules chosen for the parametrization. This leads to a slightly modified function to be minimized in order to obtain the parameter set (a), viz. a,b,c
s=E E I
wpj[P,,/""(a)- Pj,Pb"]2+
/
c k
aa,bb,cc
E/ I
we.k[ek.//a'c{a)- ek,/?bs12 (8)
The notation in this equation is explained in the preceding paper.' The dipole contribution to this function now includes an additional summation over its Cartesian components. The principal rotation axes a, b, and c were obtained by rotating the molecule by the two independent angles 4, and 42 in the original (arbitrarily chosen) Cartesian frame until the difference between the calculated and observedI2 moments of inertia, I/"'c, I r k , satisfied
(~(E.[I,C''~ - I P ~ ~ ] ~ ) / ~ ~ , , J=, ,O, /
m = 1, 2; n = I , 2 ( 9 )
The accuracy in the experimental geometry of a molecule determined by the value (El(I/"Ic- I r k ) * ) was used as a criterion for determining the weighting factors wp and we. The ratio of these factors, w,/we, was always chosen to be 10. A description of the nature of the modified PEOE parameter set is provided in Table VI, and the optimized numerical values are given in Table VII. Results and Discussion Prior to making use of the electrostatic potential- or field-derived net atomic charges, the importance of the choice of the location of the grid points was checked. A set of calculations for formamide was therefore carried out, calculating the potential or field at points (12) Toluene: Rudolph, H. D.; Dreizler, H.; Jaeschke, A.; Wendling, P. 2.Naturforsch. 1967, 220, 940. o-Xylene: Rudolph, H. D.; Walzer, K.; Krutzik, I. J. Mol. Spectrom. 1973, 47, 314. Phenol: Pedersen, T.; Larsen, N. W.; Nygaard, L. J. Mol. S r r c t . 1969.4, 59. Aniline: Lister, D. G.; Tyler, J. K.; Hog, J. N.; Larsen, N. W. J. Mol. Struct. 1974, 23, 253. Pyridine: Sorensen, G . 0.;Mahler, L.; RastrupAndersen, N. J. Mol. Struct. 1974, 20, 119. 2-Methylpyridine: Dreizler, H.; Rudolph, H. D.; Miider, H.2.NaturJorsch. 1970, 25a, 25. 4-Methylpyridine: Rudolph, H. D.; Dreizler, H.; Seiler, H. Z . Naturforsch. 1967, 22a, 1738. Pyrrole: Nygaard, L.; Nielsen, J. T.; Kirchheiner, J.; Maltesen, G.; RastrupAndersen, J.; Ssrensen, G. 0. J. Mol. Slruct. 1%9, 3, 491. N-methylpyrrole: Arnold, W.; Dreizler, H.; Rudolph, H. D. 2. Naturforsch. 1968, 23a, 301. Imidazole: Griffiths, J. H.; Wardley, A,; Williams, V. E.; Owen, N. L.; Sheridan, J. Nature 1967, 216, I301
Determination of Net Atomic Charges. 2
The Journal of Physical Chemistry, Vol. 94, No. 1I , 1990 4743
TABLE V Net Atomic Charges and Dipole Moments of Ionic Molecules Calculated from Modified PEOE, PD, and Mulliken Analysis (Charges in e, Dipole Moments in D, and Quadrupole Moments in D A)
PEOE NH,' CH3N H3'
HCOi
CH,COi
glycine zwitterion
alanine zwitterion
DP -0.289 0.322 2.68 0.117 0.078 -0.3 15 0.321 2.184 (3.89, -1.94) 2.49 0.765 -0.794 -0.177 1.310 (5.60, 3.15) 3.63 -0.139 -0.001 0.717 -0.788 3.445 (6.93, 0.39) 4.43 0.158 0.019 0.876 -7.83 -0.41 1 0.302 12.~78 (6.73, 2.10) 7.67 0.149 0.020 -0.0 14 0.012 0.877 -0.783 -0.4 12 0.303 11.561 (10.39, -1.94) 14.13
DQP -0.340 0.350 2.19 0.147 0.075 -0.4 19 0.349 2. I87 (3.98, -1.99) 2.66 0.766 -0.795 -0.177 1.317 (5.61, 3.16) 3.59 -0.116 -0.006 0.713 -0.789 3.476 (6.92, 0.42) 4.59 0.192 0.014 0.875 -0.784 -0.5 15 0.329 12.077 (6.77, 2.08) 7.87 0.175 0.015 0.003 0.010 0.874 -0.784 -0.5 17 0.330 1 1.578 (10.35, -1.94) 13.34
PD" 6-31G** -0.845 0.47 I 0.45 (0.04)j 0.249 0.040 -0.313 0.3 14 2.185 (3.77, -1.88) 2.32 (0.16) 0.855 -0.825 -0.248 1.316 (5.91, 3.16) 2.64 (0.51) -0.119 -0.020 0.909 -0.862 3.425 (7.77, 0.05) 2.21 (0.17) 0.065 0.051 0.819 -0.763 -0.404 0.314 12.225 (7.39, 1.59) 3.41 (1.15) 0.245 0.07 1 -0.359 0.1 18 0.883 -0.772 -1.199 0.516 1 1.468 (10.36, -1.97) 10.97 (2.64)
net at. charge ( N O N M S ) b
0.398 -0.005 0.797 -0.768 -1.164 0.507 12.219 (7.57, 1.44) 5.24 0.205 0.023 -0.290 0.091 0.828 -0.763 -0.412 0.299 1 1.474 (10.07, -1.93) 12.55
Mulliken' 6-31G** -0.653 0.413 1.14 -0.207 0.217 -0.641 0.399 1.263 (2.205)' (4.61, -2.31) 6.55 0.583 -0.736 -0.112 1.379 (1 ,460)' (4.61, -2.31) 7.03 -0.398 0.079 0.695 -0.768 3.556 (3.823)c (6.53, 0.81) 6.77 -0.187 0.183 0.753 -0.710 -0.624 0.375 12.016 (1 2.196)c (6.82, 1.12) 12.26 -0.101 0.179 -0.351 0.146 0.783 -0.721 -0.610 0.368 11.613 (11.561)c (9.61, -2.31) 20.63
"The point charges of H atoms bonded to the same atom are assumed to be equal. b N e t atomic charge set obtained at an NONMS. 'Dipole moments in parentheses are calculated from wave functions. dRoot mean square value defined as NQ
rms = [{x(vab(ri)i- 1
Wi))*1/W"*
and the values in the table are given as rms X IO4. 'Two values of diagonal components of quadrupole moment. /Percent error defined as NQ
NQ
PE = 100{xlVab(ri) -v(ri)~~/~~~vab(ri)~~ i- 1 i-I encapsulated by spherical shells of differing thickness with respect to the van der Waals surface. The results are presented in Table 11. In general, there is only slight variation between charge sets irrespective of the specific method or shell size chosen. An exception to this trend is the result obtained for the set of points chosen to lie between 2 and 3 A from the van der Waals surface. These charges differ considerably from all the other sets, although they closely resemble a set of charges found to be at a local minimum for a set of similarly derived charges in the previous paper.' Notwithstanding the peculiarities of the charge set just discussed, the uniformity of the charges presented in Table I1 suggests that preference need not be given to a particular choice of technique or sampling. The choice made, however, for the remainder of the study to determine charge distributions for ionic species was to employ the PD method, sampling the region of space 0-3 A from the van der Waals surface. It is clear, however, that the choice is arbitrary, and the situation remains that sampling a specified region with respect to the van der Waals surface may
lead to an erroneous description of the potential in regions not included as part of the fitting procedure. The implications of misrepresenting the long-range part of the potential in describing the interactions between distant fragments of a protein have been succinctly outlined by Price et al." The potential-derived charges for the chosen set of ionic molecules, calculated by using small (STO-3G, 6-3 1G) and moderate (6-31G*, 6-31G**) size basis sets, are presented in Table I. There are considerable variations in the charge sets between the small and the moderate size basis sets. It is doubtful whether the charges calculated by using the small basis sets are of much utility as a result of considerable inaccuracies in these wave functions and presumably, therefore, deficiencies in the calculated electrostatic potential. This is not surprising given the small variational freedom afforded by these basis sets and their lack (13) Price, S.L.;Harrison, R.J.; Guest, M. F.J . Comp. Chem. 1%9, 10,
552.
4744
The Journal of Physical Chemistry, Vol. 94, No. 1 1 , I990
TABLE VI: Atomic Species and Ihmping Factors Used for Aromatic Molecules in the Modified PEOE Method Atomic Species car,a C in aromatic molecules N atom in aromatic molecules as in pyridine \ N//
Damping Factor
f,, between atoms in aromatic molecules between H atom and the atoms in aromatic molecules fi2
f,, between atoms in aromatic molecules and atoms not in aromatic system except H
'The superscript a r represents the carbon atom in an aromatic molecule. TABLE VII: Optimized PEOE Parameters for the Aromatic Molecules D P DQPb atom card
-N= -N-
I
11, /I 2 /-I3
a: 12.784 16.977 18.226
b: 16.551 5.423 6.642 0.48 0.74 0.54
b: 16.887 7.879 8.984
a:
11.721 18.920 17.948 0.45 0.70 0.55
"DP5 parameters are used for the atoms not included in the aromatic system. DQP5 parameters are used for the atoms not included in the aromatic system. ?The parameters ai and bi are those of eq 6. dThe superscript ar represents the carbon atom in an aromatic molecule.
of polarization functions. The point charges obtained by using the 6-31G* and 6-31G** basis sets differ little from each other, the addition of the single polarization function to the hydrogen atoms in set 6-31G** not appearing to play a significant role. It is the 6-31G** basis set PD calculations from which the modified PEOE parameters were obtained in the manner already described in the Methods section. Before commenting further on the results of the PD charges obtained using this basis set, we outline reasons for its choice, reemphasizing that despite the fact that it is widely used, it is perhaps best described as being of moderate quality. Further, some ab initio calculations have been undertaken on ionic species known to be problematic, largely because of the requirement for additional diffuse functions necessary to obtain very accurate descriptions of the extra electron in anions. Even addressing these difficulties, the influence of electron correlation remains imponderable at the Hartree-Fock level of calculation. Despite these apparent flaws, it is felt that, given the inherent approximations and limitations both in the PD procedure and in the monopole descriptions of charge distributions that we are seeking, the increased computational expense of more accurate ab initio calculations is not warranted. Should deficiencies become apparent in our point charge distributions involving ionic groups, the above points may be investigated. Returning to a discussion of the nature of the 6-31G** PD charge distribution, a key feature observed for the anionic species is that the dipole moment associated with the charge set is always less than that obtained directly from the wave function. This reflects a further problem in using monopole distributions since, in the case of nonbonding valence electrons, the region of greatest electron density is displaced from the center of the atom. A more realistic description might be to locate the center of charge at a point equidistant from the nonbonding electron pairs. The PD charge distribution of the glycine zwitterion is surprisingly different from that of the similar alanine zwitterion. For example, the charge on the N atom in alanine is almost 3 times greater in magnitude than in glycine, a similar discrepancy being present when comparing the charges on the C" atoms. The presence of a near-optimal (nonminimum) solution (NONMS) was identified
No et al.
in the fitting procedure which, for the alanine molecule, produces a charge distribution closely related to that of the globally optimized glycine charge distribution. These results suggest problems both in seeking a unique determination of the PD charges and in transferring them between similar sites in different molecular environments. It is to be recalled that an advantage of PEOE procedures in general is that they involve the direct calculation of approximate charge distributions. It is possible, therefore, to account for the distortion in charge distributions when molecular fragments, such as amino acid residues, are assembled, modeling the redistribution of charge, a feature not accounted for by transferable charge schemes. The optimal modified PEOE parameters, obtained by using inforniation from the PD charge distributions as described previously, are given in Table IV. The table includes the initial atomic charges required to start the iterative procedure. As for neutral molecules, atoms [C(O-), C"] have low values of the electronegativity coefficient b. The damping factors, characterizing certain properties of the bond type, lie in a narrow range of values, 0.46-0.70. Table V illustrates the point charge distributions calculated by using the modified PEOE parameter set, the PD method, and the often used Mulliken population analysis. In the presentation of the results from the PD method for glycine and alanine zwitterions, the values of charge sets found at an NONMS are also given. As already noted, it is striking how the globally optimized charge set of one molecule resembles more closely the charge set obtained at an NONMS in the fitting function of the other. Table V also reports the quality with which each charge set reproduces the ab initio electrostatic potential in terms of the function F (see eq 2). As expected, the PD charges reproduce the electrostatic potential most closely (smallest value of E ) . It is fairly obvious that deviations from the ab initio electrostatic potential for these charges arise mainly because of their inability to describe complex anisotropic features of the potential, particularly in regions close to the van der Waals surface. The typical deviations in the electrostatic potential, calculated by using PD charges, from the ab initio potential are between 0.45 X and 10.79 X hartree. These deviations correspond to percentage errors in the range 0.04-2.5%. Analogous errors computed for the modified PEOE charges are in the range 0.25-3.4%, whereas those for the Mulliken charges are 0.1-5.0%. Calculations of the electrostatic potential surface are illustrated graphically in Figure 1 for the alanine zwitterion ion in the plane accommodating the N, C", and C(0) atoms. The surfaces obtained from the PEOE charge sets, DP (Figure lb) and DQP (Figure IC), are almost coincident with the PD surface (Figure la). This encourages us to believe that the modified PEOE procedure has sufficient flexibility to reproduce any electrostatic potential used as part of the constraints in the parametrization of the method. The surface calculated from the Mulliken charge distribution (Figure Id) is in poorer agreement with that obtained from the PD charges, most notably in the region of the methyl group. This, together with the previously mentioned poor numerical result describing the agreement of this potential surface to the ab initio surface, seems to be further indication of the folly in misusing Mulliken population analysis in modeling inter/intramolecular interactions. In a more quantitative study describing the computation of electrostatic potentials, Price et aI.l3 assert that application of Mulliken analysis in computing, say, relative energies of large molecules is based on unfounded optimism. The modified PEOE parameters for the aromatic molecules are given in Table VII, and multipole moments calculated by using the resultant point charge distribution given in Table VIII. Agreement with experiment for the calculated dipole moments of toluene and o-xylene and the calculated quadrupole moment of benzene is very good, particularly for the DQP parameter set. Likewise, the DQP parameter set leads to very accurate predictions for the dipole moments of both pyrrole and N-methylpyrrole. This is satisfying since it indicates that the charge distribution between the methyl group and the aromatic ring is modeled accurately.
Determination of Net Atomic Charges. 2
-2
-4
.
.
-4
.
-2
0 ANGSTROM
'
The Journal of Physical Chemistry, Vol. 94, No. 1 1 , I990 4145
2
-4
4
-2
a
2
'
4
-4
-2
0 ANGSTROM
2
0 ANGSTROM
2
4 C
4
d Figure 1. Electrostatic potential surfaces of the alanine zwitterion in the plane containing the N, C, and C(0)atoms, calculated for several net atomic charge sets (Table V): (a) potential-derived(PD) charge set obtained with the 6-31G**basis set; (b) modified PEOE charge set calculated with the DP parameter set; (c) modified PEOE charge set calculated with the DQP parameter set; (d) 6-31G** Mulliken point charge set. The unit of electrostatic ANGSTROM
b
potential is kJ/mol. TABLE VIII: Calculated Dipole and Quadrupole Moments Compared with Experimental Data (Dipole Moments in D and Quadrupole Moments in D A)
benzene toluene o-xylene phenol aniline pyridine 2-methylpyridine 4-methylpyridine pyrrole N-methyl yrrole imidazole
r
modified PEOE DP DOP 0 (3.21, 3.21, -6.24)' 0 (2.81, 2.81, -5.62) 0.40 0.38 0.64 0.63 1.12 (pa = 0.06, pb = 1.12) 1.08 (pa = 0.05, @b = 1.08) = 0.96, pb = 1.29) 1.61 (Ia 1.32 (pa = 0.78, pb = 1.06) 2.15 (10.56, -4.04, -6.52) 2.16 (10.43, -4.38, -6.05) 1.99 (pa = 0.63, pb = 1.90) 1.95 (/La = 0.47, pb = 1.86) 2.5 1
1.50 2.12
4.50 (pa = 4.50, pb = 0.12)
2.51 1.73
exotl data 0 (2.80, 2.80, -5.60)
0.38 0.64 1.23 (pa 0.13, fib = 1.22) 1.87 (pa = 1.07, pb = 1.53) 2.23 (9.7, -3.5, -6.2) 1.86 (pa = 0.72, pb = 1.72) 2.70 1.74 2.12
3.80 (pa = 3.7, pb = 0.9)
OThe values in parentheses represent the diagonal quadrupole moments. *Geometrywas taken from the crystal structure. The moments of inertia from the crystal data are Ia = 54, Ib = 54, and 1, = 106, and the experimental values are In = 52, = 58, and I, = 110 (I in units of amu/A2).
For those molecules containing nitrogen or oxygen atoms that possess a nonbonding electron pair, viz., phenol, aniline, pyridine, and 4-methylpyridine, Table VI11 shows that calculated dipole moments are always smaller than the experimental values. This is a limitation in atom-centered monopole distributions, already alluded to in the discussion pertaining to ionic species, namely, the unsuitability in choosing the atomic site as the expansion center for charge density associated with nonbonding electrons. A possible remedy is the introduction of a new pa-
rameter describing the displacement of the point charge, along some appropriate line connected to the atom center. In contrast to pyridine, the predicted dipole moment of 2-methylpyridine is slightly overestimated. This can be attributed to the withdrawing of electron density from the ring by the methyl group and hence to a concomitant decrease in the electron density of the nitrogen nonbonding electron pair. Unlike the case for pyridine, the nitrogen atom is now a more suitable site to utilize as an expansion center for the local charge distribution. It is a satisfying feature
4146
The Journal of Physical Chemistry, Vol. 94, No. I I I990 I
-0.783 0.303
P-
0.020
0.303I H o.Oi2
I
-0.014
, I
0.012
H
0.012
Figure 2. Partial charges of the alanine zwitterion computed by the modified PEOE method, with the DQP parameter set.
of our PEOE method that it models, qualitatively at least, features of the hyperconjugation effect associated with the methyl group bonded to an aromatic ring. The predicted dipole moment of imidazole is very poor, most likely because of the necessity, in the absence of the gas-phase geometry, to use the geometry of the crystal structure. In the crystal structure of imidazole, there are intermolecular hydrogen bonds that modify the geometry with respect to the gas phase, probably lengthening the N-H bond, and hence are the likely source of the error in the (overestimated) predicted dipole moment. For the aniline molecule the N atom possess an atypical valence angle,I2 being well described by neither sp2 nor sp3 hybridization states. Consequently, according to our classification of electronegativity parameters (Table VI), none are suitable for describing this atom. As a result, the modified PEOE charge distributions fail badly to reproduce the experimental dipole moments for this molecule. Ultimately this molecule was excluded from the optimization of the revised PEOE parameters given in Table VII. Briefly therefore, with the exception of aniline and imidazole, the modified PEOE parameters are found to lead to charge distributions capable of reproducing experimental dipole and, where available, quadrupole moments.
Summary The work presented here completes a scheme designed to facilitate the direct computation of approximate molecular charge
No et al. distributions in proteins, by including functional groups that are ionic or aromatic in nature. No assumptions are made concerning the transferability of net atomic charges between molecular fragments; hence the method accounts in an approximate way at least for redistribution of charge density upon their assembly. It should be clarified, however, that, when the PEOE formalism is used, the charge density is assumed to depend only on the connectivity and not on conformation. While this is computationally advantageous, the method is unable to model changes in the charge density occurring as a result of distortions in protein geometry significantly different from that of the standard geometries of the molecules used in obtaining the modified PEOE parameters. Recognizing the undoubted superiority of approximating molecular charge density as a set of distributed m~ltipoles,’~ we nonetheless believe that the modified PEOE method will be of utility in describing charge distributions in proteins. The modified PEOE method leads to the direct calculation of atomic charges, requiring almost negligible computational resources, with the fundamental parameters calculated largely from experimental data. Figure 2 provides an example of a set of conformation-independent charges for the alanine zwitterion; a similar diagram was presented in paper I for the corresponding terminally blocked neutral molecule. Acknowledgment. We thank Drs. R. L. Williams and D. R. Ripoll for many helpful discussions. This work was supported by the US.-Korea Cooperative Science Program between the National Science Foundation and the Korea Science and Engineering Foundation (NSF Grant INT-87-05307), by the National Institute of General Medical Sciences, National Institutes of Health (GM-143 12), and by the National Science Foundation (DMB84-01811). Support was also received from the National Foundation for Cancer Research. K.T.N. thanks the Korea Science and Engineering Foundation for support. The research was conducted using the Cornell National Supercomputer Facility, a resource of the Cornell Theory Center, which receives major funding from the National Science Foundation and the IBM Corp., with additional support from New York State and members of the Corporate Research Institute. Registry No. NH4, 14798-03-9; CH3NH3+,17000-00-9; HC02, 7 I 47-6; CH3C02-,7 1-50- 1 ; glycine, 56-40-6; alanine, 56-41-7; formamide, 75- 12-7; benzene, 7 1-43-2; toluene, 108-88-3; o-xylene, 95-47-6; phenol, 108-95-2; aniline, 62-53-3; pyridine, 110-86-1; 2-mettylpyridine, 10906-8; 4-methylpyridine, 108-89-4; pyrrole, 109-97-7; n-methylpyrrole, 36-54-8; imidazole, 288-32-4. (14) Stone, A. J.; Price, S. L. J . Pbys. Cbem. 1988, 92, 3325.
