Modeling amino acid side chains. 2. Determination of point charges

Christophe Chipot, Janos G. Angyan, Bernard Maigret, and Harold A. Scheraga. J. Phys. Chem. , 1993, 97 (38), pp 9788–9796. DOI: 10.1021/j100140a042...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. 1993,97, 9788-9796

9788

Modeling Amino Acid Side Chains. 2. Determination of Point Charges from Electrostatic Properties: Toward Transferable Point Charge Models Christophe Chipot, JBnos G. hgyhn,' and Bernard Maigret' Laboratoire de Chimie Thtorique, Unite de Recherche Associte au CNRS No. 510, Universitt de Nancy I, BP. 239, 54506 Vandoeuure- 12s- Nancy Cedex, France

Harold A. Scheraga' Baker Laboratory of Chemistry, Cornell University, Zthaca, New York 14853-1 301 Received: March 2, 1993'

An analysis of the charge transferability from one naturally occurring amino acid side chain to another structurally similar one is presented. Point charges fitted to the self-consistent-field ab initio electrostatic potential and to the potential created by distributed multipole series are compared, using split-valence type 6-3 1G** wave functions. S C F electrostatic potential-derived net atomic charges are transferable only for polar fragments, where the contribution of the high-order terms of the multipole expansion is hidden by large low-order moments. For nonpolar fragments, where the high-order terms of the multipole expansion dominate, a limited number of off-atom charges must be included in order to reproduce the higher-order moments and, hence, attain a reasonable transferability of the charge models. Conversely, a recently proposed method, consisting of an analytical fit of the multipolar part of the exact molecular electrostatic potential, allows calculations to continue to be made within the framework of atom-centered point charges. This analytical procedure, nevertheless, provides a consistent transferability from one molecule to another structurally related one, together with poorer, but still acceptable, reproduction of the electrostatic potential.

1. Introduction

For an accurate description of intermolecular interactions in classical molecular dynamics (MD) or Monte Carlo (MC) simulations, it is necessary to have a proper parametrization of the nonbonded van der Waals contribution, which involves repulsion and dispersion terms, together with the electrostatic part representing the Coulomb and induction forces. The electrostatic interactions, generally described in terms of atomcentered monopole expressions, are essential for a reliable representation of strong, anisotropic intermolecular forces, e.g., hydrogen bonding. Unfortunately, the notion of atomic charge is not defined rigorously in quantum chemistry. Therefore, any assignment of charges on a given set of atomic sites involves conceptual difficulties. Several alternative procedures for computing point chargeshave been proposed,1-25 but none of them overcomes these difficulties completely satisfactorily. Traditional approaches involving the partitioning of the electron density between the constituent atoms of a particular molecule1~z~26 have been shown to reflect the polarity of both atoms and functional groups reasonably well but are inadequate to reproduce molecular electrostatic properties11233(e.g.multipole moments, electrostatic potential, and its successive derivatives). In contrast, several works based on the use of the molecular electrostatic potential have recently been reported in the literature6,10.11,1~18,24 and are more promising. The electrostatic potential corresponds to the expectation value of a one-electron operator and is available directly by quantum chemical computation. In the first part of this study,24we considered various naturally occurring amino acid side chains. A complete model molecule was proposed by simply adding an extra hydrogen atom to each side chain -R group of the amino acid residue-NH-CHR-CO(see Figure 1 of ref 24). A full geometry optimization of these species was carried out, using the extended split-valence type 6-31G** basis set.27s28In the specific case of the negatively *Abstract published in Advance ACS Abstracts, August 15, 1993.

charged side chains (viz. aspartate and glutamate), diffuse functions were in~luded28.2~ (i.e. 6-31++G**). The net atomic charges were fitted to the ab initio self-consistent-field (SCF) electrostatic potentials and fields evaluated on a grid of regularly spaced points surrounding the molecule.6J6 At that stage, we observed that (i) there is no discrepancy between potential- and field-derived point charges, (ii) the sensitivity of the charges to the grid density is negligible, as soon as the effects of electron cloud penetration are avoided and the space surrounding the molecule is properly amp led,^^.^^ and (iii) the lack of charge transferability, especially in aliphatic compounds (viz. alanine, valine, and isoleucine side chains), suggests that the method is inadequate to describe the molecular charge distribution of nonpolar compounds. In a recent we have shown that a simple atom-centered point charge model is not suitable for reproducing the electrostatic potential of saturated hydrocarbons, where the first nonvanishing terms of the multipole expansion correspond to high-order moments (viz. typically octopoles, hexadecapoles, etc.). These observations have, therefore, led us to reinvestigate the problem of net atomic charges from a different point of view. Since the conventional fit of atom-centered charges to the molecular electrostatic potential has proved to be inappropriate to ensure an acceptable transferability, which is a major requirement for the development of force fields, other possibilities have been explored. On the one hand, the inclusion of off-atom point charges increases the accuracy of the fit and leads to an improved transferability. On theother hand, the fit to the potential created by distributed multipole moments (DMM), using the procedure introduced by Ferenczy,21 seems to offer a reasonable compromise between transferability and accuracy. A brief description of the methodology together with the technical details of the computations are recapitulated in section 2. The net point charges of the naturally occurring amino acid side chains derived from ab initio SCF electrostatic potentials and DMM potentials,using 6-3 1G** wave functions, are analyzed

0022-365419312097-9788%04.00/0 0 1993 American Chemical Society

Modeling Amino Acid Side Chains in section 3. The inclusion of a limited number of additional off-atom sites, as an attempt to improve our model, is discussed in section 4.

2. Overview of the Methods: Computational Details Molecular Electrostatic Potential-Derived Charges. The procedure employed in this work consists of fitting a series of point charges to the quantum mechanically calculated electrostatic potential, using a constrained numerical least-squares technique. A complete description of the method is presented in the first part of this study." In this section, we consider the influence of the grid definition on the values of the charges. Our sampling algorithm is very similar to the one originally employed by Cox and Williams6 and, more recently, in the CHELPG program.16 The selection of points begins with the definition of a parallelepiped-shaped grid of regularly spaced points, containing the molecule. The dimensions of the grid are defined by means of an arbitrarily chosen exclusion radius rmsx on all sides. Points found to lie further than rmax from any nucleus are systematically excluded. In addition, points inside the van der Waals envelope are rejected. Because of the short-range effects of electron cloud penetration, which generally prevents a correct representation of the SCF potential by point charge models, the size of this envelope has been d0ubled.2~ This fact is particularly critical since regions where the effects of penetration are nonnegligible often overlap with the range of action of high-order multipole moments (uiz. typically hexadecapoles). It seems therefore that such moments may not be taken into account correctly in the course of the fitting procedure. As has been shown r e c e n t l ~ ?this ~ may lead to quite significant errors in the description of the molecular charge distribution of nonpolar systems (e.g. saturated hydrocarbons), where higher moments are the first nonvanishing terms of the multipole expansion. The use of this a l g ~ r i t h m ~may J ~ .be~ criticized ~ for two major reasons: (i) since the exclusion radius r,, is a constant, whatsoever the atom, the grid is expected to be thicker in the vicinity of atoms having a small van der Waals radius. In other words, are the net point charges affected by an uneven thickness of the grid? (ii) Moreover, is there any noticeableinfluenceof the point density on the charges? The largest van der Waals radius used in the course of this study is 1.8 A, which corresponds approximately to that of the sulfur atom. Because van der Waals radii have been doubled in order to avoid the effects of electron cloud penetration and in order to take into account the multipolar contributions of higherorder than dipoles, it is necessary to choose an optimum exclusion radius. It appeared that rmax = 4.0 A is a fairly good compromise and yields satisfactory re~ults.~4 However, it may be observed that, in the vicinity of hydrogen atoms, the grid thickness is about 1.6 A whereas, in the vicinity of sulfur atoms, the thickness drops to 0.4 A. A series of tests involving different exclusion radii (rmX = 5.0, 6.0, and 10.0 A) have been carried out but did not show any particular dependenceof the net charges on r- in this range. In order to confirm the integrity of our results, we have carried out additional calculationsof potential derived net atomic charges, employing a different sampling algorithm. The GEPOL program's32 was used to generate van der Waals surfaces around the molecule. Since, in this procedure, the number of points per van der Waals sphere is limited to 60,30 the use of a simple envelope proved to be insufficient to yield acceptable charges. For this reason, a series of 17 concentric layers, separated by a grid step of Ar = 0.1 A, have been defined, so that the average number of points for these grids and those used in the course of the derivation of the SCF electrostatic potential are roughly equivalent. It is worth noting that, for consistency, the van der Waals radii have been doubled, and the furthermost envelope is located at approximately 4.0 A from any nucleus.

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9789 Our criteria for the quality of the fit are the root-mean-square deviation (rmsd) between the ab initio quantum mechanically calculated and the point charge-derived electrostatic potentials r

.

Nmt

1112

together with the corresponding mean absolute deviation

I

\-I

where Npt denotes the number of points in the grid. P c F (r,,) is the reference SCF molecular electrostatic potential evaluated at a point p, located at r,, whereas VQ (r,) is the potential created by the net atomic charges. The values of the rmsd are given in atomic units. However, in order to correlate our results with those available in the literature, which are generally expressed in kcal mol-', these values should be multiplied by a factor of 627.5095 (hartrees to kcal mol-'). The mean absolute deviation A is given in percent. The constrained least-squares fit of the point charges to the SCF molecular electrostatic potential was carried out with the GRID program.33 For each side chain, the exclusion radius and the grid step were set at 4.0 and 0.5 A, respectively. The electrostatic potential was computed on a grid of points, using the ab initio package GAUSSIAN 92.34 Table I reveals a marked homology of the two sets of potential derived net atomic charges (qv), even for molecules involving atoms with large van der Waals radii (e.g. cystine side chain). Values of the rmsd, as well as mean absolute deviations, are also fairly similar, hence indicating that the two samplings are equivalent. This leads to the conclusion that there is essentially no real influence of both the definition and thickness of the grid on the point charges, as long as the space surrounding the molecule is properly sampled.23.24 These remarks raise another question: what criterion should be used to decide whether the sampling is optimum? Recent have shown that the influence of the grid density on the charges becomes insignificant after a certain threshold, which may be characterized by a stabilization of both the rmsd and mean absolute deviation. From this threshold, increasing the density of points (e.g. by reducing the grid step Ar) does not modify either the values of the charges or the error in the fitting procedure. Derivation of Point Charges from Distributed Multipole Moments. An alternative approach for the evaluation of the molecular electrostatic potential is the use of the distributed multipole moment (DMM) appro xi ma ti or^,^^^^ which is known to provide a good representationof this quantity. In this approach, the electrostatic potential can be split into a short-range contribution (Le. electron cloud interpenetration part) and alongrange contribution (Le. multipolar part). In the present a p proximation, the effects of penetration are systematically neglected, so that the multipoleexpansionconverges formally toward the correct potential only for points located outside the actual charge distrib~tion.3S.3~As with SCF electrostatic potentials, DMM potentials can be used in a numerical fit procedure in order to obtain point charge models.23 However, in the present study, we have applied an analytical procedure, proposed recently by Ferenczy,zl avoiding the explicit evaluation of the potential on a predefined grid of points. This procedure consists of a least-squares fit of an array of net atomiccharges to the potential created by a distributed multipole series. The starting point for this type of computation is a distributed multipole analysis (DMA),35-38which has been shown to be a powerful tool to provide a compact and practical

Chipot et al.

9790 The Journal of Physical Chemistry, Vol. 97, No. 38, 1993

TABLE I: Influence of the Grid Defiition on the Potential-Derived Net Atomic Charges. of Several Amino Acid Side Chains,

Using the Split-Valence 6-316" Basis Set GRIDc

GEPOLd

Ag

Ph (D)

0.140 0.140 0.140 0.140 0.049 37.34 0.000

Arginine; flpyD= 4434,

CI HI H2 H3

c2

H4 H5

c3

H6

H7 NI H8 c 4

N2 N3 Hg

-0.425 0.138 0.116 0.116 0.203 0.007 0.007 -0,005 0.080 0.080 -0.488 0.322 0.909 -1.085 -0.902 0.507

0.145 0.145 0.145 0.145 0.048 23.52 0.000

Hio Hi1 Hi2

rmsd ( l P 3 au) A (%)

P

(D)

$FD = 4434,

GEPOLd = 4595 0.517 0.465 0.438 0.226 0.148 5.780

0.515 0.466 0.438 0.219 0.151 5.782

atom Cystine; H2 H3 H4 H5 H6

GRIP

HI H2 H3

c2

01 0 2

rmsd (10-3 au)

A (W P (D)

CI s1 s2

c2 HI

-0.303 0.020 0.036 0.020 0.935 -0.854 -0.854 0.325 0.221 3.799

Cystine; flp?ID = 2824, -0.488 -0.085 -0.085 -0.488 0.214

= 2545 -0.278 0.015 0.032 0.015 0.925 -0.854 -0.854 0.314 0.205 3.846

rmsd ( l P 3 au)

qnyL = 3028 -0.550 -0.079 -0.079 -0.550 0.229

P

(Dl

Histidine D; CI Hi H2 H3 c 2

NI C3 N2 c 4

H4 H5 H.5

rmsd (10-3 au) A (%) P

(D)

GEPOLd

flZD= 2824, NavyL = 3028

A (%)

Aspartate, 6-31GS*; flp?lD = 2447, c1

4595 -0.466 0.146 0.127 0.127 0.215 0.012 0.012 -0.021 0.087 0.087 -0.526 0.331 0.950 -1.103 -0.907 0.512

GRIDC

atom Arginine;

HI H2 H3 H4 rmsdf(10-3 au)

4v

4v

4v

atomb

0.182 0.177 0.214 0.182 0.177 0.803 56.70 2.362

= 3295, -0.693 0.205 0.205 0.188 0.229 -0.478 0.247 -0.573 0.046 0.380 0.121 0.123 0.391 12.52 4.037

0.200 0.200 0.229 0.200 0.200 0.880 70.92 2.397

flpyL = 3148 -0.745 0.217 0.217 0.200 0.256 -0.474 0.229 -0.568 0.043 0.378 0.126 0.120 0.385 8.545 4.023

= 4.0 A and Ar = 0.5 A. Outermost van der Waals envelope at 4.0 a In electronic charge units ( e a ) . Indices given in Figure 1 of ref 24. r,, A; envelope separation of 0.1 A. e Number of points requested for the fitting procedure. /Root-mean-square deviation between the ab initio computed and the point charge derived electrostatic potential. 8 Mean absolute deviation. * Dipole moment computed from potential derived net atomic charges.

representation of the charge distribution. A detailed description of our approach is given in previous papers.21s2s In this method, the explicit evaluation of the electrostatic potential, or field, on a discrete set of points is replaced by an analytical integration over different regions of the space surrounding the molecule. The main difference between this procedure and the conventional fit of net point charges to the molecular electrostatic potential is that partial charges are fitted to the potentials generated by individual multipole series rather than total charges to the sum of potentials. This approximation implies that the nondiagonal site-site couplings are neglected,Zs so that the global fitting procedure can be split into several independent fits to the site multipole series. Inorder toassess thequality ofthe fit, theelectrostatic potential created by the net atomic charges derived from distributed multipoleseries has been computed on a grid of points.24 Because this regular grid is identical to the one used in the conventional SCF molecular electrostatic potential derivation procedure (the exclusion radius and the grid step have been set at rmx = 4.0 A and Ar = O S A, respectively),a quantitative comparisonbetween the two sets of charges is possible. For each amino acid side chain, the rmsd and the mean absolute deviation A, given by eqs 1 and 2, were evaluated. Since electrostatic potential-derived point charges are, by definition, the most adequate to describe the molecular electrostatic potential, in a least-squares sense, we obviously do not expect to obtain a better rmsd with the charges derived from DMM potentials. However, a small difference in the rmsd is evidence for the ability for both sets to reproduce the electrostatic potential with comparable accuracy. The net atomic charges derived from DMMs were obtained from the MULFIT program.39 The distributed multipole analysis was carried out with a program written in the OSIPE environment.4 using the weight functions introduced by Vignt-Maeder and Cla~erie.~' The one-electron density matrix over a basis of Gaussian type orbitals (GTO) was extracted from the checkpoint files generated by GAUSSIAN 92.34 The distributed multipole

expansion has been limited to the hexadecapole (/ = 4). The boundaries of integration were fixed at p1 = rvdW + 2.0 A and p2 = rvdw + 5.0A, so that a sufficient number of neighboring sites could be taken into account during each individual fit.21

3. Results and Discussion The molecular charge distribution of the naturally occurring amino acid side chains has been reproduced from ab initio electrostatic potentials and from the potentials created by distributed multipole series,using the extended split-valencetype 6-31G** basis set. The geometry of each molecule presented in this work has been fully optimized at the same level of approximation. The Cartesian coordinates were reported in the first part of this s t ~ d y . 2 ~ The whole set of amino acid side chains may be divided roughly into threecategories: nonpolar, polar, and intermediatemolecules. Table I1 demonstrates that the best fits correspond to polar compounds, as indicated by an average small rmsd and mean absolute deviation (e.g.arginine, aspartate, or lysine sidechains). Conversely, the worst fits are obtained with typical nonpolar molecules (e.g. isoleucine, leucine, or valine side chains). The alanine side chain (i.e.CHI, in our model) is an exception because of its tetrahedral symmetry. It may be observed that the potentialderived net atomic charges of these aliphatic species do not reflect any consistent transferability of the carbon-hydrogen bond polarity: the charges of the alanine side chain reveal a clear Ch-H*+ polarity,42 which is unexpectedly inverted in the case of the -CH2- groups of the valine and isoleucine side chains. Concerning the ionic compounds,it appears that only the charges of the polar part are transferred from one molecule to another structurally similar one (e.g.the-COO-fragment of the aspartate and glutamate side chains). In this particular case, themultipole expansion is obviously dominated by the total charge. The results for nonionic, but strongly polar, speciesmay be partially explained by the magnitude of the molecular dipole moments that are reported in Table 111. In brief, for such polar compounds, the

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9791

Modeling Amino Acid Side Chains

TABLE II: Comparison between Potential-Derived (qv) and Distributed Multipole-Derived (QDMM) Net Point Charges’ of Amino Acid Side Chains, Using the Split-Valence 631G** Basis Set; Influence of the Inclusion of Additional Sites Qv atom*

atomic centers only

atomic and bond centersc

Alanine; Nptd = 2263 C -0.561 HI 0.140 H2 0.140 Ha 0.140 H4 0.140 rmsdc ( lkaau) 0.049 Af(%I 37.34 Arginine; NWt= 4434 -0.425 -0.585 0.138 0.156 0.116 0.133 0.116 0.133 0.203 -0.651 0.007 0.126 0.007 0.126 -0.005 -0.763 0.080 0.183 0.080 0.183 H7 -0.488 -0.722 NI 0.322 0.433 HE 0.909 3.234 C4 -1.085 -1.021 N2 -0.902 -0.880 0.507 0.525 0.515 0.528 0.466 0.505 0.438 0.501 0.460 0.460 0.623 -1.038 -0.755 XC4NI -0.898 0.219 0.108 rmsd (lk’au) A (%I’ 0.15 0.08 Asparagine; Npt = 2926 -0,598 -0.592 Hi 0.161 0.188 H2 0.163 0.206 Ha 0.161 0.188 c2 0.965 2.546 0 -0.639 -0.447 N -1.112 -1.158 H4 0.454 0.490 H5 0.445 0.489 XClCl -0.590 -0.840 XCP -0,481 XClN rmsd ( l k a au) 0.175 0.105 A (%) 4.55 2.04 c1

4v

0.172 -0.043 -0.043 -0.043 -0.043 0.558 133.16 0.126 -0.003 -0.033 -0.033 0.127 -0.037 -0.037 0.362 -0.041 -0.041 -0.779 0.365 1.322 -1.117 -1.114 0.477 0.484 0.485 0.487

atomic and bond centersc

4v

Aspartate, 6-31++G**; Npnt= 2447 0.037 0.138 -0.033 c 2 1.075 2.687 1.172 0 1 -0.927 -0.701 -0.989 0 2 -0.927 -0.701 -0.989 XClCl -0.293 XClOI -0.852 XClO, -0.852 rmsd ( l k a au) 0.416 0.170 0.748 A 0.28 0.10 0.48 Ha

Cysteine; Npnt = 2234 -0.361 -0.935 0.132 0.252 0.206 0.321 0.132 0.252 S -0.323 -0.403 H4 0.215 0.218 XCS 0.295 rmsd ( l k 3 au) 1.056 0.938 A (W 46.43 45.51

C HI Hz Ha

0.764 0.55

0.096 0.015 0.049 0.015 -0.436 0.260 2.243 124.50

Cysteine (+ Lone Pairs); Nmt= 2234 -0.617 -0.621 0.126 0.167 0.173 -0.005 0.211 0.217 0.027 0.167 0.173 -0.005 1.935 2.113 0.204 H4 0.043 0.038 0.121 XCS -0.085 -0.953 -1.003 -0.234 XS’ -0.953 -1.003 -0.234 rmsd ( l k a au) 0.313 0.314 0.677 A 19.46 19.55 33.01 C HI H2 Ha S

-0.104 0.015 0.038 0.015 0.961 -0.715 -1.034 0.433 0.390

0.753 21.64

Aspartate, 6-31G**; Npt = 2447 -0.303 -0.628 -0,124 HI 0.020 0.099 -0.035 Hz 0.036 0.112 -0.022 Ha 0.020 0.099 -0.035 c2 0.935 1.822 1.053 01 -0.854 -0.718 -0.918 0 2 -0.854 -0,718 -0.918 Xcica -0.083 XCPI -0.492 X C A -0.492 rmsd ( l k 3 au) 0.325 0.130 0.681 A 0.22 0.09 0.46 Aspartate, 6-31++G**; Nmt= 2447 -0.340 -0.706 -0.100 0.037 0.138 -0.033 0.044 0.142 -0.029

Cystine; Npnt = 2824 -0.488 -0.939 -0.085 -0.325 s2 -0.085 -0.325 cz -0.488 -0.939 Hi 0.214 0.287 H2 0.182 0.248 Ha 0.177 0.253 H4 0.214 0.287 Hs 0.182 0.248 Hs 0.177 0.253 XClSl 0.403 XSA 0.148 XSICI 0.403 rmsd (10-3 au) 0.803 0.493 A (56) 56.70 29.31 CI SI

0.148 -0.189 -0.189 0.148 0.022 0.004 0.016 0.022 0.004 0.016

1.097 56.59

Cystine (+ Lone Pairs); Npnt= 2824 -0.576 -0.632 0.205 1.668 1.532 0.811 sz 1.668 1.532 0.811 CZ -0.576 -0.632 0.205 HI 0.185 0.194 -0.008 Hz 0.148 0.158 -0.044 Ha 0.160 0.174 -0.018 H4 0.185 0.194 -0.008 Hs 0.148 0.158 -0.044 H6 0.160 0.174 -0.018 XClSl 0.072 Xsisa -0.002 0.072 -0.806 -0.760 -0,480 -0.779 -0.738 -0.466 -0.806 -0.760 -0.480 -0.779 -0.738 -0.466 0.875 0.154 rmsd ( l e a au) 0.170 54.89 7.75 7.68 A (96) CI SI

atom*

qDMM

xs

c1

CI HI Hz

atom*

4DMM

atomic centers only

atomic centers only

atomic and bond centers0

Glutamine; Npnt= 3428 -0.244 -0.586 0.057 0.119 HI 0.067 0.134 H2 0.067 0.135 Hs 0.071 -0.712 CZ 0.004 0.178 0.000 0.171 0.824 2.014 -0.623 -0.548 -1.117 -1.084 0.429 0.470 0.466 0.473 0.504 -0.368 -0.469 -0,431 XC;N rmsd (lkaau) 0.234 0.1 17 A (W)’ 6.49 4.63 c1

qDMM

0.192 -0.058 -0,039 -0.039 -0.059 -0.01 1 -0.013 0.967 -0.7 17 -1.059 0.398 0.437

0.716 35.15

Glutamate, 6-31G**; Npnt= 2973 -0,109 -0,613 0.206 -0.043 0.054 -0.120 0.006 0.105 -0,059 0.006 0.105 -0.059 cz 0.175 -0.781 -0.127 H4 -0.085 0.106 -0.030 H5 -0.085 0.106 -0.030 ca 0.812 1.859 1.053 01 -0.852 -0.681 -0.926 0 2 -0.824 -0.714 -0.909 XCICl 0.596 XCICI -0.090 XC30l -0.573 -0.478 XCIOl rmsd ( l e 3 au) 0.322 0.068 0.652 A (46) 0.23 0.04 0.46 CI HI H2 Ha

Glutamate, 6-31G++G**; Npnt= 2978 -0.046 -0.671 0.395 -0.063 0.060 -0.172 -0.012 0.116 -0.097 -0.012 0.116 -0.097 c2 0.248 -0.903 -0.225 H4 -0.104 0.139 -0.002 H5 -0,104 0.139 -0,001 c3 0.899 2.713 1.159 01 -0.919 -0.663 -0.993 0 2 -0.887 -0.683 -0.967 XClCl 0.723 XCIC3 -0,298 XC30, -0.928 XC30, -0.858 rmsd ( l k 3 au) 0.416 0.112 0.784 A 0.29 0.07 0.60 CI HI Hz H3

c 1

Hi H2 H3 c2

Ni

ca N2 c 4

H4 H5 Hs XCIC, XCZG XCaNi

Histidine D; NPnt= 3295 -0.693 -0.681 0.205 0.145 0.205 0.145 0.188 0.175 0.229 -0,301 -0.478 -0.931 0.247 -0.733 -0.573 -1.430 0.046 -0.715 0.380 0.434 0.121 0.224 0.123 0.225 0.530 -0,019 0.302

0.217 -0.035 -0.035 -0,029 -0.002 -0.523 0.495 -0.639 0.062 0.362 0.024 0.101

Chipot et al.

9192 The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 TABLE 11 (Continued) av

qv atomic atomic centers and bond atomb only centersC q D M M Histidine D; Npt = 3295 0.560 XN~C, 1.028 XC~N~ 1.042 XN~C~ rmsd (lt3 au) 0.391 0.104 0.901 A 12.52 3.63 44.66 Histidine E; Npnt= 3346 -0.902 -0,703 HI 0.254 0.153 H2 0.254 0.153 H3 0.210 0.155 c2 0.549 -0.727 NI -0.601 -1.465 c3 0.250 -0.611 N2 -0.471 -0.939 c4 -0,242 -0.474 H4 0.123 0.215 H5 0.185 0.233 H6 0.390 0.445 XClCI 0.593 XCIC, 0.045 XC~N~ 1.044 XN~C~ 0.987 XC~N~ 0.462 XN~C, 0.437 rmsd (lt3 au) 0.385 0.11 1 A 12.84 4.55

c1

atomic atomic centers and bond atomb only centersc Leucine; Npt = 3596 c1 -0,581 -0.635 HI 0.132 0.134 H2 0.132 0.134 H3 0.144 0.135 c2 0.551 -0.418 H4 -0.035 0.078 c3 -0.581 -0.635 0.178 HS 0.132 0.134 -0.018 H6 0.132 0.134 -0.018 H7 0.144 0.135 -0.581 -0.635 -0.057 c4 0.213 Hs 0.132 0.134 -0.687 H9 0.144 0.135 0.499 Hio 0.132 0.134 -0.464 XCICI 0.347 -0.171 XCIC, 0.347 0.027 XCIC, 0.347 0.146 rmsd (lt3 au) 0.108 0.070 0.353 A (%I 158.10 41.93

0.922 45.97

Histidine (Cation); Npl = 3541 -0.564 -0.623 0.130 Hi 0.199 0.202 0.020 0.202 0.020 H2 0.199 0.226 0.026 H3 0.213 -0.433 0.114 c2 0.259 NI -0,199 -0.939 -0.459 -0.099 0.498 c3 0.009 -0,957 -0.399 N2 -0.188 c4 -0,228 -0.989 -0.088 0.417 H4 0.379 0.486 0.265 0.101 Hs 0.263 0.207 H6 0.259 0.334 0.412 0.494 H7 0.398 0.200 XClCl 0.571 XCIC, 0.625 XC~N~ 0.317 XN~Ci 0.307 c1 XC~N~ 0.812 HI XNE, rmsd (lt3 au) 0.204 0.904 Hz 0.096 0.07 0.69 A( I ) 0.13 H3

c1

Isoleucine; Npl = 3656 Cl 0.178 -0.744 HI -0.056 0.114 H2 -0,043 0.108 c2 -0.199 -0.622 H3 0.038 0.119 H4 0.042 0.123 H5 0.039 0.115 c3 0.178 -0.744 H6 -0.043 0.108 H7 -0.056 0.114 c4 -0.199 -0.622 Hs 0.043 0.123 H9 0.039 0.115 HIO 0.038 0.119 XCICI 0.519 XCIC, 0.534 xc,c. 0.519 rmsd (lt3 au) 0.331 0.076 A (%) 230.89 100.08

c2

Lysine; Npt = 4157 -0.263 -0.666 0.100 0.174 0.156 0.067 0.067 0.156 0.178 -0.691 -0,032 0.148 -0,032 0.148 -0.890 0.067 0.169 -0.016 -0.016 0.169 -0.844 0.350 -0.010 0.217 -0.010 0.217 -0,433 -1.061 0.320 0.453 0.343 0.459 0.320 0.453 0.459 0.516 0.568 0.691 0.088 0.436 0.06 0.30 Methionine; Npl = 3368 -0.338 -0.709 0.104 0.141 0.124 0.192 0.124 0.192 -0,204 -0.801 0.243 0.158 0.243 0.158 -0.143 -0.309 -0,643 -1.043 0.199 0.278 0.260 0.298 0.199 0.278 0.567 -0.006 0.436 0.587 0.385 75.18 28.83

av

qDMM

0.128 -0.051 -0,051 -0.056 0.178 -0,089 0.128 -0.05 1 -0,050 -0.056 0.128 -0.05 1 -0,056 -0,050

0.493 866.2 0.115 -0.004 -0.03 1 -0.031 0.155 -0.048 -0.048 0.076 -0,036 -0,036 0.251 -0.002 -0,002 -0.118 0.248 0.261 0.248

atomic atomic centers and bond atom* only centersc qDMM Methionine (+ Lone Pairs); Npt = 3368 S 2.050 2.461 1.082 c3 -0.697 -0.618 0.172 H6 0.166 0.154 -0.042 H7 0.224 0.184 -0.004 H8 0.166 0.154 -0.042 XCICI 0.592 XClS -0.129 0.018 k, -0,981 -1.089 -0.628 xs XS’ -0,981 -1.089 -0.628 rmsd (lt3 au) 0.260 0.106 0.618 A (%) 16.69 18.17 50.13 Phenylalanine; Npt = 3630 -0.559 -0.568 0.153 0.147 0.148 0.143 0.148 0.143 0.365 -0.090 -0.298 -0.382 0.163 0.173 -0.112 -0.148 0.145 0.158 -0.201 -0.266 0.152 0.171 -0.112 -0.148 0.145 0.158 -0.298 -0.382 0.163 0.173 0.171 0.091 XCIC, 0.183 XCIC, 0.000 XCf, 0.000 XC,G 0.183 XGC, 0.091 rmsd (lt3au) 0.105 0.085 A 56.47 45.79 ~~

N CI

0.656 0.50

c2

c3 c 4 c5

0.126 01 -0.040 0 2 -0.013 HI -0.013 0.181 -0.014 -0.014 -0.365 0.150 -0.010 0.022 -0.010

H4 0.149 H5 -0.063 S -0.069 0.133 -0.051 -0.045 -0.054 0.149 -0.069 -0.063 1.478 0.133 236.00 -0.045 -0.054 Methionine (+ Lone Pairs); Npt = 3368 -0.051 CI -0.281 -0.700 0.189 HI 0.075 0.139 -0.053 Hz 0.081 0.158 -0.043 H3 0.081 0.158 -0.043 0.437 c2 0.119 -0.685 0.123 513.20 H4 -0.011 0.146 -0.041 Hs -0.011 0.146 -0.041

0.188 -0.027 -0,035 -0.035 0.003 -0.155 0.102 -0.083 0.098 -0.119 0.101 -0.083 0.098 -0,155 0.102

C HI H2 H3

Proline; Npl = 3924 -0.667 -1.392 0.145 -0.651 -0.070 -0.721 -0.049 -0.622 0.057 1.353 0.755 2.999 -0,640 -0.561 -0,604 -0.301 0.466 0.334 0.028 0.136 0.017 0.082 0.025 0.118 0.052 0.116 0.168 0.050 0.147 0.039 0.073 -0.022 0.454 0.509 0.768 0.507 0.568 -0.396 -1.233 -0.840 -1.198 0.122 0.582 15.78 123.88 Serine; NWt = 2464 0.278 -0.337 -0.027 0.083 0.036 0.144 -0.027 0.083

0.923 63.69 -0.852 0.398 0.099 0.095 0.241 0.886 -0.716 -0.682 0.376 -0.049 -0.094 -0.05 1 -0.039 -0.034 -0.027 -0,018 0.465

0.742 174.40 0.628 -0.120 -0.069 -0,120

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9793

Modeling Amino Acid Side Chains

TABLE II (Continued) 4v

atomb 0

H4 xco

rmsd (lt3 au) A

c1

-0.685 -0.799 0.426 0.433 0.394 0.542 0.246 23.71 16.31

HI Hz 0

H3 c2

H4 H5 H6

XClCZ

xc20

rmsd (lt3 au) A (W)

c1 Hi Hz Ho c2

c3

N

c4

c5

0.445 21.84

-0.641 0.165 0.103 -0.833 0.425 -0.657 0.122 0.140 0.135 0.534 0.508 0.209 11.76

-0.723 0.404

-0.655 0.148 0.145 0.148 -1.036 -1.309 -1.597 2.452 3.661 0.248 0.490

Cg Hg

-0.263 0.177 -0.207 0.156 -0,142 0.150 -0.281 0.169

c7

0.726 40.33

H7 Cs Hs c9

0.504 -0.050 -0.095 -0.703 0.397 -0.012 -0.022 -0.017 -0.002

H9 Xclc2 Xcg,

XC~C, XC,N XNC~ Xc4c2

XC,~, XC,G XGC7

XCC,

0.575 29.75

Tryptophan; NPt = 4203 -0.584 0.165 0.168 0.165 0.063 -0.154 -0.536 0.198 0.153 0.194 0.410

atomic atomic centers and bond atomb only centersC Tryptophan; Npt = 4203

qDMM

Threonine; Nnp(= 2964 0.497 -0.012 -0.082 -0.736 0.423 -0.226 0.044 0.031 0.060

4v

4v

atomic atomic centers and bond only centersc Serine; Npt = 2464

0.225 -0.035 -0.042 -0.035 -0.182 0.026 -0.475 0.241 -0.040 0.110 0.347

XC,C, rmsd(1t'au) A (%I

0.121 12.68

-0.370 0.191 -1.131 0.210 -1.251 0.216 -1.097 0.238 0.582 -0.766 1.122 1.417 0.034 -3.566 -1.208 0.704 1.079 1.123 -0.223 0.061 9.24

qDMM

-0.107 0 0.107 Hg -0.148 c.5 0.102 H7 -0.095 C7 0.102 HE -0.220 XcIc, 0.120 Xclc7 XCZCl XCIC4

xcs, %SO

XC,G XGC7

rmsd (lt3 au) A (%)

Hi H2 H3

CZ Cj H4

C4

-0.576 0.158 0.148 0.158 0.288 -0.198 0.166 -0.407 0.188 0.434

-0.565 0.146 0.148 0.146 1.034 -0.191 0.171 -0.531 0.185 0.191

-0.638 -0.670 0.449 0.477 -0,264 -0.763 0.192 0.233 -0.278 0.037 0.179 0.163 -0,042 -0,635 -0.497 0.239 0.046 -0.008 0.597 0.089 0.115 0.066 23.21 18.38

qDMM

-0.619 0.428 -0.241 0.142 -0,105 0.107

0.999 173.00

Valine; Npnt= 3252 0.155 22.09

Tyrosine; NPt = 3800

Ci

atomic atomic centers and bond atomb only centersc Tyrosine; NPt = 3800

CI HI H2 CZ H3

0.200 -0.036 -0.036 -0.036 -0.057 -0.094 0.103 -0.299 0.12 0.421

H4

Hs C3

Hg H7 H8 Xclcl XcIc,

0.385 -0.081 -0.081 -0.249 0.048 0.041 0.048 -0.249 0.048 0.048 0.041

-0.725 0.112 0.112 -0.619 0.116 0.123 0.116 -0.619 0.116 0.116 0.123 0.514 0.514 0.079 44.77

0.166 -0,067 -0.067 0.129 -0.050 -0,044 -0.050 0.129 -0.050 -0.050 -0.044

H5 rmsd(1t3au) 0.307 0.479 3 Hs C5 A (%) 213.84 307.60 * In electronic charge units (ecu). b Indices given in Figure 1 of ref 24. Bonds involving heavy atoms only. Number of points requested for the H4

fitting procedure. 8 Root-mean-sqauredeviation between the ab initio computed and the point charge derived electrostatic potential. f Mean absolute deviations. influence of higher-order moments is hidden by the long-range effects of the low-order moments, which roughly correspond to the predominant term of the multipole expansion. Therefore, the fit of the net atomic charges to the electrostatic potential does not involve significant errors. However, it has been pointed out that saturated hydrocarbons are inadequately described by net atomic charges,25 because such a simplistic model is unable to reproduce correctly the high-order multipole moments, uiz. hexadecapoles,which are particularly large in these compounds. As a result, a noticeable lack of charge transferability is observed for nonpolar groups; e.g.we expect the charge borne by the carbon atom in the -CH2- fragment of the glutamate side chain to be relatively close to the one in the terminal -CH3 fragment of the aspartate side chain. Unfortunately, this is not true, and the situation becomes even worse for nonpolar molecules,as indicated in Table 11. On the contrary, the net atomic charges, fitted to the potential created by the distributed multipole series, are consistently transferable from one molecule to another structurally similar one,25with no particular variations in magnitude [uiz.qc = 0.172 ecu (electronic charge unit) for the alanine side chain and 0.129 and 0.133 ecu for the terminal carbons of the valine and isoleucine side chain, respectively]. Our results show an invariant Cd+-Hb polarity, unless the carbon is connected directly to a polar group (uiz. the charge borne by the carbon of the terminal -CH3 fragment in the aspartate side chain and of the - C H r fragment in the glutamate side chain are qc = -0.124 and -0.127 ecu, respectively). This interesting fact is in apparent contradiction with chemical intuition and would suggest that hydrogen is more electronegative than carbon." It should be understood that only

methods based on a partitioning of the molecular electron density can provide charges correctly reflecting the electronegativity differences of bonded atoms. In contrast, net atomic charges fitted to ab initio SCFor DMM potentials depend on the influence of all surrounding multipoles and, therefore, are not necessarily correlated with electronegativities. Compared with electrostatic potential derived charges, the present sets of charges generally yield larger errors in the fitting procedure, for both aliphatic and ionic compounds (see Table 11). In addition, a less accurate reproduction of the dipole moments computed from DMM potential-derived charges is observed throughout Table 111. At this stage, one may wonder why the net atomic charges fitted to the electrostatic potential failed to be transferable from one molecule to another structurally similar one and why this problem is apparently solved with DMM potential-derived charges. (i) We have concluded that the fit of net atomic charges to the electrostatic potential is particularly reliable in the case of ionic species because the molecular dipole moment is large enough to conceal the effects of higher-order multipole moments. The description of such moments, uiz. hexadecapoles, which are conspicuously predominant in saturated hydrocarbons (e.g. alanine, isoleucine, and valine side chains), appears to be most difficult if only simple atomic charge models are considered. As a result of the constraint imposed on point charges to reproduce high-order moments, very large rmsd and mean absolutedeviations are usually found when treating aliphatic species. Relaxing the fitting procedure significantly would imply the inclusion of a large number of extra off-atom point charges. This is not exactly the case, as will be seen in section 4.

Chipot et al.

9794 The Journal of Physical Chemistry, Vol. 97, No. 38, 1993

TABLE III: Dipole Moment. of Amino Acid Side Chains from Potential-Derived and Distributed Multipole-Derived Net Atomic Charga, Using the Split-Valence 6-31C** Basis Set; Comparison with SCF Vaues and Influence of the Inclusion of Additional Sites YV

atomic centers side chain alanine arginine asparagine aspartate,C6-31G** aspartate,' 6-31++G**

{ :$%(+ lone pairs) { (+ lone pairs)

only 0.000 5.782 4.064 3.799 4.348 1.812 1.769 2.362 2.387 3.888 5.509 6.149 4.037 3.538 3.025 0.057 0.093 9.067 1.780 1.787 0.290 1.455 1.837 1.792 1.893 1.435 0.055

atomic and bond centersb

PDMM

5.785 4.061 3.794 4.342 1.816 1.77 1 2.377 2.383 3.894 5.495 6.133 4.039 3.545 3.026 0.072 0.093 9.069 1.792 1.793 0.291 1.463 1.833 1.794 1.896 1.439 0.06 1

O.Oo0 5.838 4.107 3.901 4.476 1.800 1.775 2.440 2.391 3.929 5.604 6.265 4.041 3.571 3.088 0.073 0.094 9.158 1.759 1.792 0.282 1.492 1.782 1.769 1.919 1.416 0.06 1

PEP O.OO0 5.783 4.061 3.793 4.335 1.783 2.391

3.891 glutamine 5.494 glutamate,c6-31G** 6.133 glutamate,' 6-31++G** 4.040 histidine D 3.537 histidine E 3.025 histidine (cation)c 0.069 isoleucine 0.090 leucine 9.071 lysine methionine 1.800 methionine (+ lone pairs) 0.292 phenylalanine 1.466 proline 1.835 serine 1.793 threonine 1.893 tryptophan 1.439 tyrosine 0.060 valine a Indebyeunits. b Bondsinvolvingheavyatomsonly. Chargedspecias; dipole moment calculated in the stondord orientotiod3 framework of Gaussian 92.34

(ii) It is known that the electrostatic potential may be expressed as an exact sum of the potential generated by a finite set of multiples (Le. the long-range contribution) and a penetration part decreasing exponentially from the nucleus (i.e. the shortrange contribution), which is neglected in the DMM approximation. Since the multipole expansion offers a correct representation of the potential outside the actual molecular charge distribution,35J7we expect the net atomic charges obtained from the distributed multipole series to be practically similar to those fitted to the exact electrostatic potential.23 This is not what is observed in Table 11, where most sets of net atomic charges obtained from both methods are substantially different. In fact, the DMM potential-derivedcharges have been fitted analytically, neglecting nondiagonal site-site couplings.25 In Table IV, we compare the charges obtained from the distributed multipole series by means of an analytical and a numerical fit. In both cases, the DMMs were computed up to the hexadecapole. For the numerical procedure, use was made of a grid of regularly spaced pointsz4(identical to that employed for potential derived net atomic charges), on which the DMM potentials were evaluated. It may be seen that the charges obtained from the numerical fit are in excellent agreement with those derived from the ab initio SCF electrostatic potential, hence confirming that the distributed multiple series provide a good representation of the exact potential outside the molecular charge distribution. From this, we conclude that the transferability of the net atomic charges obtained analytically from the DMM potentials is artificial, because global fits have been approximated by a series of partial fits, neglecting sit-ite couplings.25 So far, we have learned that electrostatic potential-derived charges are satisfactory to mimic the ab initio SCF potential with reasonableaccuracy,dependingon thenatureof the molecule. However, they do not reflect any particular transferability,

especially in the case of nonpolar fragments. On the contrary, net atomic charges fitted to DMM potentials by means of an analytical procedure are consistently transferable but less successful in reproducing the molecular electrostatic potential. This suggests that a set of point charges suitable for the reproduction of electrostatic quantities is apparently inappropriate to describe bond polarity. That this is not completely true may be seen in the following section.

4. Refmement of the Model The criteria defining an acceptable point charge model are its ability to describe the molecular electrostatic potential (or its derivatives), its charge transferability from one molecule to another structurally similar one, and finally, its tractability for a possible application in MD or MC simulations. As noted in the previous paragraphs, the problem of transferability is directly related to the notion of polarity. For nonpolar species, usually characterized by large high-order multipole moments (uiz. hexadecapoles), a simple net atomic charge model is obviously inadequate to reproduce small electrostatic potentials. Unfortunately, the inclusion of a sufficient number of extra sites to relax the fitting procedure significantly is likely to yield a model of very low tractability. For this reason, we have adopted a compromise, which consists of introducing one additional site per chemical bond involving heavy atoms (i.e. other than hydrogen). In doing so, we do not suggest that the location of the dummy charges is independent of the electron density in the bond. In principle,the exact position of the extra off-atom point charges should be optimized along each chemical bond. However, in a recent paper,25such a simple model was applied to a series of saturated hydrocarbons and led to very satisfactory results. This corresponds to the partitioning of the electron charge distribution to the various contributions due to localized orbitals, as has been proposed by C l a ~ e r i eIn .~~ the present work, our goal is to obtain consistently transferable charges from one fragment to another structurally similar one, hoping that the molecular electrostatic potential will not be misrepresented. From Table 11, it may be observed that our new procedure provides far more reliable point charge models, reproducing the electrostatic potential with greater accuracy (uiz. reduction of the rmsd by a factor of 4.7 for the glutamate side chain) and providing a complete transferability from side chain to side chain. This important fact may be illustrated by an example: if the alanine side chain (i.e. methane) is assumed to be the prototype of all saturated hydrocarbons, we expect the charge borne by the carbon (qc = -0.561 ecu) to be perfectly transferred on each side chain involving aliphatic parts. Indeed, the terminal methyl groups of such side chains fulfill this requirement (e.g. qc = -0.585, -0.6 13,-O.622,a n d 4 5 6 8 ecu for the arginine, glutamate, isoleucine, and phenylalanine side chain, respectively). In addition, it is worth noting that the molecular dipole moments calculated from these point charges are practically identical to those computed at the SCF level (see Table 111). The electrostatic potential of molecules involving lone pairs is often wildly misreproduced by simple net atomic charge models. Stone has argued that the best possible description of lone pairs consists of including a dipole on the atom carrying the lone pair.35 Since the electrostatic part of most available force fields is described in terms of monopole expression^,^^ the convenient representation of lone pairs by point charges is still widely used. Such an approach has been proposed in the first part of this study," where a pair of adequately located dummy charges (see Figure 2 of ref 24: dsx = 0.65 A and LXSX = 135.5') was included in the conventional point charge model of each sulfurcontaining amino acid side chain (uiz. cysteine, cystine, and methionine side chain). In Table 11,we have assessed the influence of the addition of extra off-atom charges on the reliability of the

The Journal of Physical Chemistry, Vol. 97, No. 38, 1993 9795

Modeling Amino Acid Side Chains

TABLE rV: Influence of Site-Site Cou lings on the Multipole Derived Net Atomic Charges' of Several Amino Acid Side Chains, Using the Split-Valence 6-316** Basis k t qDMM

atomb C HI H2 H3 H4

rms~V((10-~ au) A8 (96)

Ph (D)

0.172 -0.043 -0,043 -0.043 -0.043 0.558 133.16 0.000

-0,560 0.140 0.140 0.140 0.140 0.050 36.80 0.000

atomb

analyticale numericald Arginine; Npt = 4434

Hio HII Hi2 rmsd (lt3 au) A (96) c (D)

-0.426 0.138 0.116 0.116 0.198 0.008 0.008 -0.001 0.079 0.079 -0.490 0.323 0.907 -1.083 -0.900 0.507

0.484 0.485 0.487 0.764 0.55 5.838

0.514 0.466 0.437 0.218 0.15 5.782

Aspartate, 6-31G**; Npt = 2447 c1

Arginine; Npt = 4434 0.126 -0.003 -0.033 -0.033 0.127 -0.037 -0,037 0.362 -0.041 -0.041 -0.779 0.365 1.322 -1.117 -1.114 0.477

qDMM

~DMM

analyticale numericald Alanine; Npl = 2263c

Hi H2 H3

c2

01 0 2

rmsd (lP3au) A (96) c (D)

-0.124 -0.035 -0.022 -0.035 1.053 -0.918 -0.918 0.681 0.41 3.799

-0.303 0.020 0.036 0.020 0.935 -0.854 -0,854 0.324 0.22 3.799

atomb

analyticale numericald Cystine; Nm1= 2824

H2 H3 H4 H5 Hs rmsd ( l t 3au) A (96) P (D)

CI SI s 2

c2 Hi

0.148 4.189 -0.189 0.148 0.022

c 2

N1 C3 N2

H4 -0.488 -0.085 -0.085 -0.488 0.214

0.182 0.177 0.214 0.182 0.177 0.801 58.72 2.363

Histidine D, Npl = 3295 C1 HI H2 H3

c 4

Cystine; Npt = 2824

0.004 0.016 0.022 0.004 0.016 1.097 56.59 2.440

Hs Hs rmsd (10-3 au) A (96) P (D)

0.217 -0.035 -0.035 -0.029 -0,002 -0.523 0.495 -0.639 0.062 0.362 0.024 0.101 0.901 44.66 4.041

-0.689 0.204 0.204 0.187 0.228 -0.477 0.246 -0,573 0.045 0.380 0.121 0.123 0.391 13.51 4.037

a In electronic charge units (ecu). b Indices given in Figure 1 of ref 24. Complete neglect of the site-site couplings K&,,,,,.25 Least-squares fit on a grid of regularly spaced points. * Number of points requested for the fitting procedure. f Root-mean-square deviation between the ab initio-computed and the point charge-derived electrostatic potential. Mean absolute deviation. Dipole moment computed from distributed multipole derived net atomic charges.

model. It appears that the rmsd for sulfur-containing molecules is reduced significantly when a modeled lone pair is included and, to a lower extent, when dummy charges are placed at the center of selected chemical bonds.

5. Conclusions The lack of transferability from one molecule to another structurally similar one, which constitutes a major drawback of potential-derived net atomic charges, has led us to reinvestigate the problem of fitted charge models from another perspective. In the present work, distributed multipole expansionswereconsidered to generate point charge models, using an analytical procedure. With this method, theexplicit evaluation of electrostaticquantities, uiz. potential or its successive derivatives, on a grid of points is no longer necessary. Net atomic charges derived analytically from DMM potentials are completely transferable and, therefore, in much better agreement with chemical intuition than those obtained from a global fit to the molecular electrostatic potentials. It has been observed that atom-centered point charges fitted to the ab initio SCF potential are satisfactory to reproduce such a quantity, in a least-squares sense. However, the accuracy of the fit together with the charge transferability depends on the polarity of the molecule. Generally very low errors in the fitting procedure and acceptable transferability of the polar parts of polar molecules (e.g. -COO- fragment of the aspartate and glutamate side chain) are obtained. This may be explained by the fact that the multipole expansion is dominated by the first few moments (Le. charges and dipoles), which mask the contribution of higher-order moments. In contrast, large rmsd and mean absolute deviationsareoften encountered with nonpolar species, viz. saturated hydrocarbons, the electrostatic potential of which is generally small and difficult to describe by simple atom-centered point charge models.25 Consequently, these charges do not reflect any real consistency from one aliphatic side chain to another similar one (e.g.C*+-H6 polarity in alanine, valine, and isoleucine side chains).

We have shown that the net atomic charges fitted analytically to the DMM potentials are totally transferable because the couplings between the DMM sites have been neglected. However, the reproduction of the electrostatic potential is less reliable with these charges than with the classical potential-derived charges. This suggests that, if a given set of net atomic charges is able to reproduce some electrostatic properties, uiz. potentials, it does not imply, afortiori, that this set is suitable to describe the polarity of the chemical bonds accurately. In order to overcome the difficulties raised by the reduced accuracy of the method developed by Ferenczy,21 a different model has been proposed, consisting of an array of atom-centered point charges, supplemented by a limited number of adequately located off-atom dummy charges. This enhanced charge model has been fitted to the full ab initio SCF potential. The results obtained with this alternative scheme are very satisfactory, providing a completetransferabilityfrom one side chain to another structurally related one, together with a significant decrease of both rmsd and mean absolute deviation. These results may certainly be improved by including some additional dummy charges that will contribute to a better reproduction of the high-order multipole moments but, unfortunately, at the expense of lowering the tractability of the model. One of the great virtues of Ferenczy's method is that it provides a reasonable way to understand the "mechanism" of the fitting procedure. As has already been discussed in a previous paper,25 the series of distributed multipole moments on a given site is fitted by the ensemble of atomic charges. This means that, if high-order moments at this particular location contribute significantly to the potential (e.g. in aliphatic compounds), the charges on relatively distant sites should be modified in order to attain the best possible fit in a least-squares sense. Nevertheless, we may consider that the given atom supplemented by its first neighbors constitutes a quasi-independent subunit, sufficient to mimic the multipole moments up to order 1 = 2 or 3, therefore explaining the transferability of Ferenczy's charges.

9196 The Journal of Physical Chemistry, Vol. 97, No. 38, I993 The transferability can be maintained as long as partial point charges are fitted to individual distributed multipole series. As has been shown global fits involve couplings of individual fits: vir. an optimal fit to the DMMs on site A may be counterbalanced by an unfavorable set of charges obtained from theseriescenteredonsiteB. Asaresult,global-fit procedures generally yield smaller errors in the reproduction of the electrostatic potential but much poorer transferability of the charges. Furthermore, the models that include off-atom point charges, which offer more flexibility to describe the anisotropy of the charge distribution, have proven to be more transferable and, obviously, ensure greater numerical accuracy. Thus, depending on the purpose, different options may be considered: (i) If atom-centered point charge models, suitable for the reproduction of electrostatic quantities, e.g. potential, field, are required,conventional molecular electrostatic potentialderived net atomic charges appear to be particularly suitable. (ii) If, however, the description of these electrostatic properties is of lesser concern, but if reasonablechargetransferability is required, analytical fits to the DMM potentials provide acceptable results. (iii) Finally, if the two preceding conditions have to be fulfilled simultaneously, it becomes necessary to include some additional off-atom charges in the original model. This may be problematical for the parametrization of force fields, where it is desirable to avoid an extensive use of dummy charges.

Acknowledgment. We thank Drs. F. Colonna and G. G. Ferenczy for many stimulating discussions and the continuous interest and encouragmentof Prof. J.-L. Rivail. C.C. is indebted to Roussel Uclaf Company (France) for his doctoral fellowship. We thank the Groupement Scientifique “Mod6lisation Mol&ulaire” IBM-CNRS for their generoussupport. Thecomputations presented in this paper were carried out on the IBM-3090/600 and on the Serial IBM-RS/6000 cluster of the Cornell National Supercomputer Facility, Cornel1 University, Ithaca, NY. This work was also supportedby the US.National Science Foundation (Grants DMB-90-15815 and INT-91-15638) and the action incitative CNRS-NSF LET92 No. 21 MDRI. References and Notes (1) Uwdin, P. 0. J. Chem. Phys. 1950, 18, 365. (2) Mulliken, R. S.J . Chem. Phys. 1955, 23, 1833. (3) Hmze, J.; Jaff€, H. H. J. Am. Chem. Soc. 1962, 84, 540. (4) Bader, R. F. W.; Beddall, P. M.; Cade, P. E. J. Am. Chem.Soc. 1971, 93, 3095. ( 5 ) Gasteiger, J.; Marsili, M. Terrahedron 1980, 36, 3219. (6) Cox, S. R.; Williams, D. E. J. Compur. Chem. 1981, 2, 304. (7) Guillen, M. D.; Gasteiger, J. Terrahedron 1983, 39, 1331. (8) Nalewajski, R. F. J . Am. Chem. Soc. 1984, 106, 944. (9) Mortier, W. J.; Genechten, K.V.; Gasteiger, J. J. Am. Chem. Soc. 1985, 107, 829. (10) Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987, 8, 894. (11) Williams, D. E.; Yan, J. M. Adv. A?. Mol. Phys. 1987, 23, 87. (12) (a) Kim, S.;Jhon, M. S.;Scheraga, H. A. J. Phys. Chem. 1988,92, 7216. (b) Wee, S. S.; Mm, S.; Jhon, M. S.;Scheraga, H. A. J . Phys. Chem. 1990,94, 1656. (13) Dinur, U.; Hagler, A. T.J . Chem. Phys. 1989, 91, 2949. (14) Beslcr, B. H.; Merz, Jr., K. M.; Kollman, P. A. J. Compur. Chem. 1990,11,431. (15) Woods,R. J.; Khalil, M.;Pell, W.; Moffat, S. H.; Smith, Jr., V. H. J. Compur. Chem. 1990,11, 297.

Chipot et al. (16) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990,11,361. (17) Ferenczy, G. G.; Reynolds, C. A.; Richards, W. G. J. Comput. Chem. 1990, 11, 159. (1 8) Luque, F. J.; Illas, F.;Orozco, M. J. Compur. Chem. 1990,11,4 16. (19) (a) No, K. T.; Grant, J. A.; Scheraga, H. A. J. Phys. Chem. 1990, 94,4732. (b) No, K. T.; Grant, J. A.; Jhon, M. S.;Scheraga, H. A. J. Phys. Chem. 1990, 94,4740. (20) Colonna, F.;hgyBn, J. G.; Tapia, 0. Chem. Phys. Lett. 1990,172,

55. (21) Ferenczy, G. G. J. Comput. Chem. 1991,12, 913. (22) Aida, M.; Corongiu, G.; Clementi, E. Int. J. Quantum Chem. 1992, 42, 1353. (23) Colonna, F.; Evleth, E.; h g y i n , J. G. J. Compur. Chem. 1992, 13, 1234. (24) Chipot,C.; Maigrct, B.;Rivail, J.-L.;Scheraga,H.A. J. Phys. Chem. 1992,96, 10276. (25) Chipot, C.; h g y i n , J. G.; Ferenczy, G. G.; Scheraga, H. A. J . Phys. Chem. 1993,97,6628. (26) Hall, G. G. Adv. At. Mol. Phys. 1985, 20, 41. (27) (a) Hariharan, P. C.; Pople, J. A. Chem. Phys. Lett. 1972,16,217.

(b) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.;Gordon, M. S.;DeFreea, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654.

(28) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Inirio Molecular Orbital Theory; Wiley: New York, 1986. (29) Chandrasekhar, J.; Andrade, J. G.; Schleyer,P. v. R. J. Am. Chem. Soc. 1981,103, 5609. (30) Pascual-Ahuir, J. L.; Silla, E. J. Comput. Chem. 1990, 11, 1047. (31) Silla, E.; Villar, F.; Nilsson, 0.;Pascual-Ahuir, J. L.; Tapia, 0. J . Mol. Graphics 1990,8, 168. (32) Silla, E.; Tuilon, I.; Pascual-Ahuir,J. L. J . Compur. Chem. 1991.12, 1077. (33) Chipot, C. GRID (Fortran program performing charge fitting to molecular electrostaticpotentials for fields; availablefrom theauthors), 1992. (34) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.;Schlegel, H. B.; Robb, M. B.; Repbgle, E. S.;Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley,

J.S.;Gonzalez,C.;Martin,R.L.;Fox,D.J.;Defreea,D.J.;Baker,J.;Stewart, J. J. P.; Pople, J. A. GAUSSIAN 92; Gausian Inc.: Pittsburgh, PA, 1992. (35) Stone, A. J. Chem. Phys. Let?.1981,83, 233. (36) (a) Price, S.L.; Stone, A. J. Mol. Phys. 1982,47, 1457. (b) Price, S.L.; Stone, A. J. Chem. Phys. Lett. 1983,98,419. (c) Price, S.L.; Stone, A. J. Mol. Phys. 1984,51,569. (d) Price, S.L. Chem. Phys. Lert. 1985,114, 359. (e) Price, S. L.; Stone, A. J.; Alderton, M. Mol. Phys. 1984, 52, 987. (37) Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047. (38) Stone, A. J.; Price, S.L. J . Phys. Chem. 1988, 92, 3325. (39) Ferenczy,G.G.MULFIT(Fortranprogramperformingchargefitting to multipole series; available from the author: Chemical Works of Gedeon Richter Ltd., H-1475 Budapest, P.O. Box 27, Hungary), 1992. (40) (a) h g y i n , J. G.;Colonna, F.; Jansen, G. To be submittedto Chem. Phys. Len. (b) Colonna, F. OSIPE (Open Structured Interfaceable ProgTanuning Environment),Dynamiquedu InteractionsMolkulaires,UniversitC Rerre et Marie Curie, 4, Place Jussieu, - Tour 22, 75005 Paris, France. (41) Vignt-Maeder, F.; Claverie, P. J. Chem. Phys. 1988, 88, 4934. (42) Righini, R.; Mali, K.; Klein, L. L. Chem. Phys. Lett. 1981,80,301. (43) Frisch, M.; Foresman, J.; Frisch,A. E. Gaussian92 User’sGuide and Gaussian 92 Programmer’sGuide; Gaussian Inc.: Pittsburgh, PA, 1992; pp 47, 341. (44) FliszBr, S. Charge Distributions and Chemical Effects; Springer-Verlag: New York, 1983; p 184. (45) Claverie, P. In Infermolecular Interacrions: From Diaromics Io Biomlvmers: Wilev: New York. 1978: D 69. (46j Momany, F. A.; McGuire, R. F; Burgess, A. W.; Scheraga. H. A. J. Phys. Chem. 1975, 79,2361. (47) (a) Weiner, S. J.; Kollman, P. A.; Case, D. A,; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, Jr., S.;Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (b) Weiner, S.J.;Kollman, P.A.;Nguyen,D.T.;Case, D.A. J. Compur. Chem. 1986, 7, 230. (48) Brooks, B. B.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swammathan, S.;Karplus, M. J . Compur. Chem. 1983,4, 187. (49) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Pcstma, J. P. M. Biopolymers 1984, 23, 1513.