129

An Approach to Computing Electrostatic Charges for Molecules

U. Chandra Singh and Peter A. Kollman*

School of Phurmag, Department of Pharmaceutical Chemise, University of California,San Francisco, California
94143
Received 10 February 1983; accepted 17 May 1983

We present an approach for deriving net atomic charges from ab initio quantum mechanical calculations
using a least squares fit of the quantum mechanically calculated electrostatic potential to that of the partial
charge model. Our computational approach is similar to those presented by Momany [ J . Phys. Chem., 82,
592 (1978)], Smit, Derissen, and van Duijneveldt [Mol. Phys., 37, 521 (1979)], and Cox and Williams [ J .
Comput. Chem., 2, 304 (1981)], but differs in the approach to choosing the positions for evaluating the
potential. In this article, we present applications to the molecules H,O, CH,OH, (CH&O, H,CO, NH,,
(CH,O),PO;, deoxyribose, ribose, adenine, 9-CH3 adenine, thymine, 1-CH, thymine, guanine, 9-CH3
guanine, cytosine, 1-CH, cytosine, uracil, and l-CH, uracil. We also address the question of inclusion of
“lone pairs,’’ their location and charge.

INTRODUCTION cule than those calculated from the actual wave-

One of the fundamental problems in simula- function.* Such Mulliken populations are also
tions on complex systems using simple analytical often very basis-set dependent. A second ap-
potential functions is how to accurately repre- proach is to use the molecular electrostatic
sent the molecular-chargedistribution. The abil- potential evaluated at points in space around the
molecule as a guide and to fit this to point-charge
ity of molecular mechanics without the use of
models. Important work in this area has been
atomic charges to reproduce the structure and
done by Scrocco and T ~ m a s i Smit,
,~ Derrison,
energies of complex hydrocarbons is impressive.’
and van Duijneveldt; and Cox and Williams?
However, the pathway for similarly accurate
Our work follows closely the approach by Cox
simulations of polar and charged molecules is not
and Williams, with some improvements in meth-
so clear. Some approache use an atom-centered
odology. We also present a number of applica-
multipole expansion,2 and given a sufficiently
tions that help answer the questions of “how,
large expansion, these can describe the molecu-
when, and whether” one should include lone-pair
lar-charge distribution nearly as accurately as
electron sites in the electrostatic model. Finally,
the wavefunction describing the electronic struc-
we present partial charge models for fragments
ture. However, for carrying out molecular me-
of nucleic acids and study Watson-Crick hy-
chanics, Monte Carlo, or molecular-dynamics
drogen bonding of the bases using these par-
simulations of large and complex molecules, a
tial charges. When combined with appropriate
simpler description of the charge distribution is
Lennard-Jones dispersion and repulsion parame-
required. Most studies3 use only the monopoles
ters, the partial charges calculated here are shown
of the atom-centered multipole expansion (net
to yield H bond energies in very good agreement
atomic charges), and this seems a reasonable
with gas phase enthalpies of association.
compromise for complex systems. The question
remains: How can one arrive at an appropriate
set of site charges for such an expansion? The
simplest procedure is to use Mulliken popula-
tions from quantum mechanical calculations, but
these are based on a simplified model of describ- METHODS
ing the electron distribution and often yield One of us (U.C.S.) has combined a number of
rather different multipole moments for the mole- the best features of several ab initio quantum
mechanical programs to create the program
*Towhom correspondence should be addressed. GAUSSIAN 80 UCSF? This program has as its

Journal of Computational Chemistry, Vol. 5, No. 2, 129-145 (1984)

0 1984 by John Wiley & Sons, Inc. CCC 0192-8651/84/020129-17$04~OO
130 Singh and Kollman

"basis" Gaussian 80.' To this is added: (1) the then calculate the electrostatic potential at these
efficient self-consistent-field (sCF) clamping and points. A ,%-matrixprocedure is used to input an
extrapolation features from the program GAMESS initial partial-charge model in which charges or
(from the National Resource of Computational their location can be varied in any combination,
Chemiatrylo, (2) the molecular-properties pack- subject to the absolute constraint of units of
age from GAUSSIAN 79," (3) the electrostatic charge and any selected restraints on the calcu-
potential program from the Pisa (Italy) group of lated partial-charge dipole moment or quad-
Tomasi et al.,12 (4) the program by Connolly to rupole moment components of the molecule.
calculate molecular surfaces,13 (5) a Levenberg- A Levenberg-Marquardt nonlinear optimization
Marquardt nonlinear least squares program to fit procedure14 then fits the potential calculated
the potential to an analytical model,14 and (6) from the partial charges to that calculated
the Morokuma energy-decomposition analysis.l5 quantum mechanically. The algorithm is robust
Thus this program contains many of the essen- and seems to have a rather wide radius of con-
tial features for studying molecular-chargedistri- vergence to the optimal-charge model.
butions and interactions. We sought to establish an appropriate crite-
Our goal is to find simple, analytical, partial- rion for the radii employed in determining the
charge models that beat reproduce the charge surface points to be used for the electrostatic
distribution of complex molecules. To find such potential calculations. For H,O, we examined
charges, we first carry out quantum mechanical the use of successive shells of various multiples
calculations on the molecules of interest. We of the van der Waals radii,13 ranging from 1.2 to
then use the Connolly surface alg~rithm'~ to 2.0 times the radius (smaller radii lead to posi-
calculate a shell or a number of shells of points tive potentials in the lone-pair region). The
at a specified distance from the molecule and calculated point-charge models are rather insen-

Table I. Atom-centered point-chargemodels'for H,O.

Basis Set
Property STO-3G' Clemb 4-31Gc 6-31C* 6-31G** DZPPf

qOg -0.616(-0.614)h -0.792 -0.938(-0.940)h -0.812 -0.794(-0.786)h -0.752

qiii 0.308 0.396 0.469 0.404 0.397 0.376

,j 1.732 2.229 2.639 2.283 2.242 2.118

9,' 1.300 1.674 1.981 1.714 1.678 1.590

a, -0.059 -0.076 -0.090 -0.078 -0.076 -0.073

QC -1.241 -1.598 -1.891 -1.636 -1.602 -1.517

01 0.51 0.57 1.08 1.81 1.54 2.15

s 4.1 3.9 5.6 10.8 10.5 15.3

'Reference 16.
bReference 17.
'Reference 18.
dReference19.
eReference 19.
'Reference 20.
gCharge on oxygen.
hCharge on oxygen found by Cox and Williams (ref. 7).
'Charge on hydrogen.
jDipole moment from charge model in debyes (the experimental value is shown in Table IV).
kQuadrupole-momentcomponents in buckinghams relative to center of mass and in principal axis system
(experimental values shown in Table IV).
'Root mean square deviation between quantum mechanically calculated and point-charge electrostatic
potential in kcal/mol.
mPercentdenation between quantum mechanically calculated and point-charge electrostatic potentials.
Electrostatic Charges for Molecules 131

sitive to which shell(s) is used within the range Table IV. Constraint values and weights.
of 1.2-2.0 times the van der Waals radius, but Pvparty.
rnol*cule H20 Q130H (CH,)O ,3320
the percent deviation between the potential
calculated with the analytical model and that wb 1.659(250)d

calculated quantum mechanically increases as the llYb

1.85(440)' -1.15kd I .300( 250)=

factor decreases. To be sure that no artifacts "Zb 2.300( 250)f
q.s 2.639(20)h - 3.30tzo)~ 0.0401 1 0 9
from being too close to the atoms entered into
Q@ -0.050 (20) -2.20 (20) -0.051 (10)
the calculations, we used shells of 1.4, 1.6, 1.8, QcS -2.589 (20) -1.20 ( 2 0 ) 0.011 (10)
and 2.0 times the van der Waals radius in all of
the remaining calculations, using symmetry "The principal axis of the moleculesis the 2 &. H20,
where possible to reduce the number of sur- (CH,),O and, CH,O are in the XZ plane with each
molecule at its experimental geometry. Methanol has its
face points. The density of points in these C - 0 bond along the X axis. The geometries used were
surfaces depended on the size of the molecule, taken from the Chem. Soc. Special Publications, Nos. 11
and 18,1958and 1965,exce t that ideal tetraedral CH,
and varied 1-5 points/A to yield typically gryg were assumed with f;(C--H) = 1.10 A.
200-300 points/molecule. pole-moment com nent in debyes, with constraint
weight in parentheses. r m e n t a l values of the dipole
moment was used for constraint.
Table 11. Effect of shells on the point charges." 'Experimental value from ref. 21.
calculated components found with 6-31G** ba-
S h e l l 8 range (van der WMlS meale sis set.
'Propertyb 1.2-1.8 1.4-2.0 1.6-2.2 1.8-2.4
"Experimental value from ref. 22.
fExperimental value from ref. 22.
-0.794 -0.794 -0.788 -0.786
gQuadrupolemoment in buckinghams with constraint
weight in parentheses. The experimental value of the
0.397 0.397 0.394 0.393 quadrupole moments were used for constraint.
2.235 2.235 2.217 2.212 hExperimental value from ref. 23.
1.678 1.678 1.664 1.660
'Experimental value from ref. 22.
jExperimental value from ref. 22.
-0.077 -0.076 -0.076 -0.076

-1.601 -1.602 -1.588 -1.584

2.22 1.54 1-16 0.86

11.9 10.5 9.8 8.9

RESULTS
"Usin 6-31G** basis (ref. 6) on H,O.
bSee'fable I for notation.
"Create shells using van der W e radius of each H2O
a t o q times scale factors which range 111 lncrements of
0.2 A between the tabulated limits. For H,O we addressed a number of questions:
What is the basis-set dependence of the derived
partial charges? How does this depend on the
Table 111. 5-point charge model for H,O. inclusion of lone pairs and/or constraints? How
propertya STC-3G 4-31C
B.8i.
6-316
s$
6-3lC** DZPP
do a 4- and 5-point model compare for H,O?
In Table I we report the values found for
-0.130 -0.700 -0.538 -0.544 -0.372
the point-charge fit (atom-centered point-charge
0.250 0.419 0.491 0.489 0.517
model) for H,O with various basis sets. The
-0.185 -0.068 -0.223 -0.217 -0.331
calculated charges are compared with those de-
1.732 2.629 2.231 2.189 2.048
termined by Cox and Williams7 and the agree-
1.290 2.059 2.275 2.256 2.440

0.016 -0.012 0.008 -0.021 -0.024

ment is very good. Our root mean square fits to
-1.274 -2.047 -2.283 -2.235 -2.416
the potential are not as low as those of Cox and
0.468 0.908 0.466 0.465 0.443
Williams because we have chosen points closer to
103.73 106.12 70.85 69.28 68.04 the molecule to use in our fitting procedure (they
0.24 0.72 0.63 0.58 0.46 use the van der Waals radius + 1.2 A). Thus we
1.8 3.2 3.7 3.5 3.3 reinforce their conclusion that, whereas the
point-charge models are rather insensitive to
~

"See Table I for notation.

bSee Table I for notation: scale factor 1.4-2.0 times the location of the points at which to evaluate
van der Waals radius. the electrostatic potential, the quality of the fit
'Charge on lone pair.
d a t a m eof lone pair from oxygen (in A). worsens if the points are chosen closer to the
"Angle in degrees that lone pairs make with HOH molecule and improves if the points are chosen
bisector (an angle of 109.5' corresponds to approxi-
mately "tetrahedral" lone pair with respect to the hy- farther away. Our range of points of 1.4-2.0
drogens as well as the other lone pair). times the van der Waals radius gives electro-
 Singh and Kollman

Table V. &point charge H 2 0 model with constraints? Table VIII. Charge models for dimethyl ether with
B l * i S set lone pairs.
Pipper t y STO-3G 4-31G 6-3IG* 6-31C** DZPP
-. Basis setb
Propertya STO-3G 4-316 6-31G*
'(0
-0.433 -0.628 -0.453 -0.453 -0.341

qH 0.289 0.506 0.508 0.506 0.522 0.8464 ( 1.101)'

40 -0.0546 ( 2.999)' 0.1259 ( 2.239)'
-0.3529 (-0.307) -0.1573 (-0.151) -0.1874 (-0.203)
-0.072 -0.192 -0.282 -0.279 -0.352
qLP
U 1.777 1.986 1.911 1.903 1.901
2,
qcr
I
0.1200
0.1051
( 0.125)
( 0.091)
0.1178
0.0890
( 0.088)
( 0.033)
0.1167
0.0828
( 0.103)
( 0.055)
%pd -0.4005
1.3593
(-0.550)
( 1.3179)
-0.1112
2.2126
(-1.5030
( 1.7111)
-0.1580
1.7856
(-1.129)
+1.403 2.430 2.450 2.433 2.516 ( 1.5205)
Qa 2.4238 ( 2.998) 4.9243 ( 3.750) 4.2915 ( 3.544)
.Q
0.048 0.235 -0.105 -0.113 -0.066 Qb -1.7776 (-2.093) -3.3991 (-2.5407) -2.7776 (-2.313)
Qb -a.6463 (-0.904) -1.5252 (-1.203) -1.5139
Qc -1.451 -2.195 -2.344 -2.320 -2.460 3op 0.428
109.85
( 0.326)
(108.75)
0.906
116.0
( 0.329)*
(102.60)
0.739
111.82
(-1.231)
( 0.329)
(102.97)
d 0.710 0.615 0.526 0.517 0.462 0 0.304 ( 0.59) 0.569 ( 2.49) 5.65 ( 1.46)
4.5 ( 9.0) 5.1 (22.5) 6.1 (15.9)
8 100.97 62.41 65.66 65.38 66.55

2.35 2.12 1.04

0 0.56 4.45
"See Table VII for notation.
s 4.5 23.6 14.1 12.9 7.4 bSee Table I for notation.
"Models derived with constraints (Table IV) in
"See Table IV for values chosen for constraint and parentheses.
Tables I and I1 for other notation. dCharge on lone pairs.
eDistance from lone pairs from oxygen.
fDistance fixed as in 6-31G* model; otherwise lone
pair will collapse very near oxygen, with very large and
Table VI. 4-point charge H,O model. opposite charges on lone pair and oxygen.
88816 setb
gAnglethat lone pairs make with COC bisector.
Propertya STO-SG 4-31G 6-3lG* 6-31Gt* DZPP h~ deviation of quantum mechanical and model
potential in kcal/mol.
0.401 0.548 0.539 0.537 0.543 'Relative percent-rms deviation.
-0.401 -0.548 -0.539 -0.537 -0.543

1.777 2.308 2.073 2.048 1.901

1.782 2.358 2.420 2.406 2.508

0.015 -0.082 -0.037 -0.051 -0.049

-1.797 -2.276 -2.383 -2.355 -2.459 Table IX Charge models for methanol.
0.248 0.188 0.283 0.279 0.332 m&eP
Property No Lone Pair. With lone Pairs
e 72.14 61.35 66.42 65.49 65.92
iCb 0.149 -0.176
2.45 1.35 1.24
a

s
1.44

11.4 12.7 8.1 7.5 7.4

1.05
$*
";.9
-0.656
-0.001
0.042
0.054
0.101

"See Tables I and I11 for notation.

p 0.424

1.657
0.386
-0.233
1.639
bSeeTable I for notation. 2 -1.144 -1.149

2 0.744
0.484
-1.228
1.200
-0.019
-1.219
3
8*
1.40
13.4
0.638
6.1
d1
Table VII. Partial-charge model for dimethyl ether. elm
0.619
119.4
02" 87.3
Basis S e t b ~

Propert ya STC-3G 6-316'

~~

4-316
~
"6-31G** basis set used (ref. 19).
(-0.2243) -0.2908 (-0.5399) -0.4187 (-0.4038) -0.3453 :Charge on carbon.
(-0.1025) 0.1107. ( 0.1423) 0.1192 ( 0.0095) 0.0253
Charge on oxygen.
dCharge on hydrogen in XZ plane.
( 0.6679) 0.0528 ( 0.0567) 0.0392 ( 0.0799) 0.0600 'Charge on hydrogens in and out of XZ plane.
( 0.0733) 0.0408 ( 0.0355) 0.0254 ( 0.0562) 0.0437 Charge on hydroxyl hydrogen.
( 1.3735) 1.3276 ( 2.2230) 1.7281 ( 1.7968) 1.5337
gCharge on lone pair.
hDipole-moment components in debyes.
( 2.5524) 3.0545 ( 5.0853) 3.8433 ( 4.4203) 3.6163 Quadrupole-momentcomponent in buckinghams.
(-1.4362) -1.7583 (-2.9991) -2.2787 (-2.5267) -2.0843 jrms deviation between quantum mechanical and model
(-1.1162) -1.2962 (-2.0862) -1.5646 (-1.8936) -1.5319
electrostatic potentials.
kRelative percent-rms deviatioq.
( 1.21) 1.32 ( 1.45) 2.89 1.83
( 1.17)
'Lone-pair-oxygen distance in A.
8 (18.2) 19.8 (13.1) 26.1 (12.8) 20.1 "Angle lone pairs makes with C-0 bond; one lone
pair above and one below COH plane.
"Molecular property; values in parentheses without "Angle lone pairs make with 0-H bond; one lone pair
dipole and quadrupole constants described in Table IV. above and one below COH plane.
bBasis set; see Table I for notation.
'Charge on oxygen.
dCharge of carbon.
eCharge on two hydrogens in XY plane.
Char e on the four out-of-planehydrogens.
gDipofemoment of molecule in debyes.
hQuadrule-moment components of molecule in
bykinghams in principal axis system.
rms deviation of quadrupole moment and point-charge
model.
Relative percentage of rms deviation.
Electrostatic Charges for Molecules 133

Figure 1. Quantum mechanical potential for l-CH, cytosine. Potential points > 5
kcal/mol are green; potential points between - 5 and + 5 kcal/mol are blue, and potential
values < - 5 kcal/mole are red.

Table X. Formaldehyde point-charge models.'

Model
atom centered 6 pointb 5 point'
property no constraint constraint' no constraint constraint no constraint constraint

-0.463 -0.435 0.152 -0.143

0.421 0.421 0.250 0.283 0.262 0.262

0.021 0.007 0.057 0.035 0.055 0.042

-0.258 -0.105 -0.186 -0.173

2.79 2.56 2.80 2.50 2.80 2.56

-0.190 -0.068 -0.023 0.101 -0.025 0.067

0.064 0.020 -0.030 -0.099 -0.032 -0.075

0.126 0.048 0.053 -0.002 0.057 0.008

1.00 1.62 0.50 1.67 0.50 1.38

6.5 10.6 3.3 10.9 3.3 9.0

- - 0.459 0.696 0.541 0.547

en - - 109.2 113.2 111.2 110.8

"6-31G** basis set.

bCharge on all atoms and two lone pairs.
'No charge on oxygen and two lone pairs.
dSee Table IV for constraint values and weights.
"Charge on oxygen.
Charge on carbon.
Charge on hydrogen.
hCharge on lone pair.
'Dipole moment in debyes.
JQuadrupole moment in buckinghams.
rms deviation of quantum mechanically calculated and model electrostatic potentials.
'Percent-rms deviation between potentials.
mDistance of lone pair from oxygen.
C-0-lone-pair angle.
134 Singh and Kollman

Table XI. Energies of water dimer and H,CO . . . HOH (kcal/mol).

Water D i m e r H2CO ...
HOH

Geometry" Atom Centeredb 5 point' Geometryd ab i n i t i o e Atom Centeredf Extendedg

Linear, 8-0 -4.7 -4.0 l i n e a r , 8-0 -3.5 -3.3 -2.8

e-3o -4.9 (-5.2) -4.4 8-65 - -3.6 -4.1
8-40 -4.5 -4.8 (-5.2) 8-60 -4.8 -3.7 (-4.1) -t.4 (-4.7)
8-75 -4.1 -4.8 8-90 -3.7 -2.2 -2.1
+30 -3.9 -3.7 +90,6-90 -1.4 -2.6 -0.4
+30,B-90 -4.3 -3.8
bifurcated (-3.81 1-3-71
cyclic 1-3-61 [-3.51

"See ref. 25 for description of geometry and results of ab znitio calculations.

atom-centered charges from DZPP basis set. Value in parenthesis is complepd optimized geometry.
Oxygen van der Waals parameter with well depth of 0.12 kcal/mol and radius of 1.8 A was used.
5-point charge model; DZPP basis set. Value in parenthesis is after complete optimization of geometry.
Oxygen van der Waals well depth of 0.12 kcal/mol and radius of 2.0 A was used.
0-H . . . 0 linear, 0 = angle between C-0 and H-0 vectors, with external water hydrogen trans to
-
C-0 . 0. + = W.means H+ bond approaching in a direction toward oxygen.
"Results of ab inztro calculabons m g 6-31G** basis set.
'Using atom-centered model 6-31G** for H,CQ DZPP for H 0;oxygen van der Waals parameters as in
atom-centered water dimer; values in parentheses are from compfete geometry optimization.
gUsing6-point 6-31G** charge model for H,CO; 5-point DZPP for H,O; van der Waals parameters as used in
extended model for (H,O),.

static potentials of sufficient range ( - 20 to + 20

kcal/mol, where 1 kcal/mol = 4.18 KJ/mol) to
Table XII. Charges for 1-CH, cytosine. describe rather well typical interaction distances
Property STO-3G-Ma
B m i e Set/llodel
STO-3C-UAb CLKK-Mc CLKH-UAd
and energies of strong and weak interactions,
with the exception of proton aflinitie.
Charge8 NI
c2
-0.187
0.859
-0.572
0.939
-0.053
0.910
-0.723
1.061 For this reason we believe that our use of
N3 -0.860 -0.791 -0.948 -0.848
c4 0.935 0.629 1.126 0.640 Connolly molecular surface to determine the
C5 -0.576 -0.230 -0.727 -0.220
C6 0.185 0.377 0.043
-0.529
0.460
-0.543
electrostatic-potential-point choices is a useful
02 -0.508 -0.518
N4
cin
-0.834
-0.289
-0.743
0.235
-1.079
-0.514
-0.923
0.286
and general one. It leads to partial charges in
EN4A
w14B
0.351
0.329
0.338
0.335
0.428
0.396
0.405
0.406
essential agreement with that from cubic grid
85
86
0.153
0.098
0.237
0.195 algorithm, but more easily ensures that the col-
EClA
EClB
0.116
0.115
0.167
0.174 lection of points chosen for potential evaluation
wcic _....
n.115 0.174
"Xf 3.517 3.509 3.293 3.272 contains approximately equal contributions from
-4.181 -4.056 -4.872 -4.667
13.336
-8.518
13.003
-8.559
13.504
-7.161
13.663
-7.942
all regions of the molecule.
-4.788
0.56
-4.444
1.64 (4.01-1
-6.343
0.90
-5.721
2.74
We further explore the dependence of charge
3.7 11.0 (27.0)j 5.7 17.2 model on choice of points for potential evalua-
Bs~o-3Gbasis set; all-atom model, i.e. all atoms in- tion in Table 11. As one can see, the point-charge
cluded in the fit of the model potential to the quantum model derived is not very dependent on the
mechanicall calculated one? "range" of the shells of points chosen.
bsT0-3G gasis set, united-atom model; hydrogen on
carbon not used in the fit of the model potential to the We next examine the ability of 5-point charge
quantum mechanical one. models for H,O to fit the electrostatic potential.
'Clementi basis set (ref. 17), all atom model.
dClementi basis set ref. 17), united-atom model. These models contain a partial charge on each
6
"Atom charges foun from fit to potential, see ref. 25
for numbering convention; geometries taken from ref. 30.
hydrogen, one on the oxygen and one in each
'Dipolemoment components in debyes; molecule in lone-pair direction with variable distance d from
XY plane N 1 4 2 bond along X axis; geometry from ref. the oxygen and angle (relative to the HOH bisec-
30. tor line). Table 111 contains the results of fitting
EQuadrupole-momentcomponents in buckinghams in
principal axis system. the quantum mechanical electrostatic potentials
hRoot mean square deviation between model and to such models. Not surprisingly, the fit of the
q y t u m mechanical in kcal/mol.
Percent-m deviation. partial-charge model to the quantum mechanical
ju and s values found by taking AA charges and potential is substantially improved (note (I and s
summing C-H hydrogen and carbon charges and plac-
ing on carbon to calculate potentid. values) over the &point model. Most intriguing
are the differences in the lone-pair angle 6 found
in the various models. The cruder basis set STO-
Electrostatic Charges for Molecules 135

Figure 2. (a) Difference between united-atom model and quantum mechanical potential for 1-CH, cytosine.
Absolute differences < 1 kcal/mol are blue; absolute differences between 1 and 2 kcal/mol are red, and
differences > 2 kcal/mol are green. (b) Difference between all-atom model and quantum mechanical potential.
Color code same as for (a).

3G and 4-31G calculations find 8 near the tetra- The other three basis sets yield quadrupole mo-
hedral value of 109.5O , whereas the better basis ments in reasonable agreement with experiment
sets 6-31G*, 6-31G**, and double zeta plus and dipole moments that are somewhat high but
polarization (DZPP) all find 8 between 50" and closer to the experimental value than 4-31G. To
70 O . The reason for this difference can be traced derive a 5-point charge model that simulta-
to the molecular electric moments calculated with neously describes well the water molecules' di-
the various basis sets. The s ~ o - 3 Gbasis calcu- pole and quadrupole moment requires lone pairs
lated a quadrupole moment far less than experi- that are inverted and on the same side of oxygen
ment, whereas the 4-31G basis gives the poorest as are the hydrogens.
agreement with the experimental dipole moment. Another option in the fitting program we use
136 Singh and Kollman

Table XIII. Charges for cytosine.“ basis gives results that do not differ greatly from
property STO-36-AA
Basis S e t l l b d e l
STO-3G-UA CLEM-AA CLEM-UA
those with the most extensive 6-31G * basis.
Charges The addition of “classical” lone pairs (Table
N1 -0.642 -0.778 -0.769 -1.060
c2 0.982 0.969 1.121 1.126 VIII) leads to a significant improvement in the
N3 -0.892 -0.782 -1.002 -0.839
c4 1.014 0.672 1.252 0.720 fit of the charge model to the electrostatic poten-
CS -0.623 -0.241 -0.796 -0.237
C6
02
0.351
-0.538
0.412
-0.536
0.327
-0.575
0.525
-0.573
tial. The lone pair remains in the classical direc-
N4
HI
-0.880
0.310
-0.775
0.370
-1.155
0.367
-0.983
0.481
tion, independent of basis set. However, the
HN4A
wi40
0.361
0.342
0.343
0.345
0.446
0.419
0.416
~ .~
0.424
charge on the oxygen varies widely depending on
H5
H6
0.152
0.063
0.234
0.131 basis set and whether a dipolar constraint has
3.873 3.841 3.869 3.827
-4.138 -4.048 -4.770 -4.615 been applied.
4.203 5.043 4.512 5.555
-0.033 -1.138 1.492 -0.158 Methanol has lower symmetry than H20 and
-4.170 -3.905 -6.005 -5.397
0
s
0.56
3.4
1.42
8.8
0.80
4.6
2.33
13.3
(CH,),O, so it is of interest to see what charge
models can be derived for this molecule. Since
“See Table XI1 for rotation. Cox and Williams have derived atom-centered
charge models for this molecule, Table IX con-
tains only the results of the 6-31G** models
is the ability to constrain the partial-charge with and without lone pairs.
model to yield specified values of the dipole or When lone pairs are added to methanol, not
quadrupole moment components; the result of surprisingly, the charge models have a much
putting on such constraints in the 5-point H,O better fit to the quantum mechanical potential.
model is described in Table IVYwith constraint The lone pairs are in the classical direction and
values and weights used in Table V. are in a similar location for all basis sets. How-
Interestingly enough, all models except s ~ o - 3 G ever, such models only converge to reasonable
now have inverted lone pairs, and because the values of lone-pair charge and distance when
constant has significantly altered the lone-pair there is no charge on oxygen. With a charge on
position in the 4-31G model, the percent-rms oxygen, the charge on oxygen gets large and
deviation for this model is the largest of all. positive and the lone pair very negative, with
However, this 4-31G model is in much better lone-pair “collapse” onto the atom center.
agreement with the models based on the more
extended basis sets than it was without the con-
H ,CO
straints (Table III), so that this constrained op-
timization method may be a useful approach in We next sought to extend our studies to for-
cases where there is some experimental data on maldehyde, H,CO, in order to determine lone-
the molecular multipole moments, but only a pair-containing charge models for this molecule,
modest (4-31G or double-zeta basis) level of the- as well as to see if the lone-pair inversion found
ory can be applied to them. for H,O with the extended basis sets occurred in
A final choice of models considered here is a
4-point model, with only hydrogen and lone-pair
partial charges. Such models are described in
Table VI. Comparing the tables, one sees a sig- Table XIV. Charga for 1-CH, uracil.”
nificant worsening of the quality of fit of the 4- Property Basis S e t h d e l
STC-3G-AA SPC-3C-UA CLEK-AA CLEM-UA
compared with the 5-point model. Charges
NI -0.159
~ .. -0.570 -0.012 -0.702
c2 0.775 0.896 0.845 1.065
N3 -0.768 -0.758 -0.957 -0.970
(CH,),O AND CH,OH c4
c5
0.834
-0.529
0.572
-0.205
1.018
-0.705
0.628
-0.237
C6 0.160 0.367 0.065 0.471
02 -0.472 -0.494 -0.504 -0.540
We next extended our studies to the molecules 04 -0.474 -0.394 -0.533 -0.408
ClM -0.329 0.239 -0.599 0.264
(CH,),O and CH,OH to see how the point-charge 83 0.334 0.347 0.403 0.429
85 0.146 0.224
models for these molecules differed from those 86 0.098 0.180
RClA 0.133 0.195
for H,O. In the case of (CH,),O, we developed HClB
RClC
0.126
0.126
0.190
0.190
point-charge models using STO-~G,4-31G, and -1.628
-3.264
-1.667
-3.130
-1.671
-3.724
-1.738
-3.510
6-31G* basis sets. The point-charge models for -8.807
5.636
-8.460
5.146
-9.968
7.593
-9.2364
6.4667
(CH3),0 without lone pairs are presented in Ta- 3.171
0.56
3.314
1.70
2.375
0.90
2.7697
2.70
7 .a
ble VII, with and without dipole constraints. As 8 5.4 16.3 23.2

one can see, with dipole constraints the s ~ o - 3 G “See Table XI1 for notations.
Electrostatic Charges for Molecules 137

Figure 3. Quantum mechanical potential for 9-CH3 guanine. See Figure 1 for color code.

H,CO as well. The partial atom-centered charges Table XV. Charges for uracil."
we derive for this molecule differ only slightly Property B a s h Set/Model
STO-3G-M STO-3G-UA CLM-M CLEkUA
from those found in the calculations of Cox and
Williams? For example, with a 6-31G** basis N1
c2
-0.573
0.836
-0.762
0.883
-0.702
0.973
-1.040
1.075
set, we found partial charges of -0.463, 0.421, N3
c4
-0.713
0.831
-0.675
0.540
0.882
1.019
-0.851
0.602
and 0.021 on oxygen, carbon, and the hydrogens, c5
C6
-0.518
0.276
-0.181
0.373
-0.680
0.288
-0.205
0.498

compared with values of - 0.459,0.423, and 0.018 02

04
-0.489
-0.474
-0.495
-0.385
-0.531
-0.534
-0.547
-0.402
HI 0.304 0.378 0.354 0.478
found previously. When we carry out these calcu- H5 0.137 0.202
H6 0.069 0.123
lations adding a constraint to force the charge "r -1.288
-3.112
-1.347
-2.987
-1.214
-3.520
-1.303
-3.323
UY
model to have the experimental dipole and Qa -12.142
9.654
-11.846
9.095
-12.939
11.583
-12.507
10.675
Qb
quadrupole moments, there are relatively small QC
0
2.488
0.52
2.752
1.36
1.356
0.85
1.832
2.14
4.7 12.3 7.1
( - 10%) changes change in the charges on the 17.8

atoms. With the s ~ o - 3 G

basis, the oxygen charges "See Table XI1 for notation.
change from -0.28 to -0.33, there is relatively
little charge in the carbon charges, and the hy-
drogen charges change from -0.02 to 0.00. In Table XVI. Charges for l-CH, thvmine.'
addition, as was found for H,O, putting on the property Basis Serlnodel
STO-3G-AA STO-3G-UA CLEM-M CLEM-UA
dipole and quadrupole constraints enables one to
Charges
force the more limited basis-set calculations to N1
Cl
-0.233
0.849
-0.739
1.113
-0.066
0.941
-0.942
1.416
have charges more closely approximating that N3
C4
-0.851
0.809
-1.012
0.980
-1.092
0.918
-1.399
1.248
found with more extended basis sets. c5
C6
-0.176
0.034
-0.595
0.551
-0.072
-0.181
-0.826
0.772
02 -0.488 0.529 -0.529 -0.600
We first considered a model for the charge 04 -0.464 -0.472 -0.503 -0.523
C5H -0.382 0.097 -0.627 0.124
distribution in H,CO in which lone pairs were HC5A
HC5B
0.119
0.111
0.184
0.170
placed in an sp2 direction around the oxygen HC5C
CIH
0.111
-0.251 0.236
0.170
-0.467 0.250
(Table X). Interestingly, in the case of s~o-3G HClA
HClB
0.118
0.103
0.172
0.145
HCLC 0.103 0.145
and 4-31G, the lone pairs nearly collapsed onto H3 0.355
0
0.370 0.448
n.rLs
0.479
"6 ....l.U
the oxygens, with large and compensatingcharges -0.981
-3.238
-0.975
-3.175
-0.947
-3.480
-0.944
-3.581
in the oxygen and lone-pair charges. However, -10.251
7.344
-10.368
7.652
-11.728
9.253
-1 1.681
9.382
2.907 2.716 2.475 2.299
the model in which the lone pairs were held fixed 0.58 1.23 1.04 2.08
6.1 12.9 9.7 19.4
a t the distance found in the 6-31G** calculation
had a fit to the quantum mechanical electro- "See Table XI11 for notation.
138 Singh and Kollman

Table XVII. Charges for thymine.’ examined the hydrogen bond directionality in
Property
SO-3C-M
Basis setlnodel
STC-3C-UA CLEWM CLEM-UA (H20), and H2C0 HOH with both atom-
~~~ ~~~

centered and lone-pair models and the results are

Charges
NI
c2
-0.624
0.901
-0.858
1.034
-0.797
1.070
-1.208
1.309
presented in Table XI. For (H20)2both models
A3
C4
-0.788
0.787
-0.862
0.858
-0.989
0.898
-1.131
1.036
correctly reproduce the fact that the trans linear
CS
C6
-0.179
0.154
-0.474
0.485
-0.095
0.104
-0.595
0.677 dimer is the lowest energy structure26and the
02 -0.504 -0.523 -0.555 -0.589
04 -0.458 -0.451 -0.497 -0.491 bifurcated structure is not a local minimum and
C5U
HCSA
KCSB
-0.370
0.117
0.107
0.084 -0.615
0.181
0.165
0.096
is -1.5 kcal/mol higher in energy than the
RCSC
HI
0.107
0.317 0.367
0.165
0.383 0.472 linear str~cture.2~ The lone-pair model has more
0.341 0.432
K3
116
0.335
0.100
0.410
0.173 “lone-pair directionality”, i.e. the lowest energy
-0.720 -0.726 -0.554 -0.561
“x
-3.122 -3.071 -3.495 -3.411 occurs at a larger 8, but both models suggest
:
VY

F
-11.337
8.739
2.598
-11.489
8.994
2.494
-12.339
10.606
1.733
-12.399
10.753
1.646
that the correct orientation of the external hy-
8
0.54
5.3
0.89
8.9
0.83
7.6
1.43
13.2
drogen on the proton-donor molecule is at least
as important as a lone-pair approach in that the
“See Table XI1 for notation. B = -30 (indeed all 8 < 0) and 8 = 30, = 90 +
geometries are higher in energy than any with
8 < 0. The importance of the external hydrogen
on H bond directionality has also been found in
static potential less than 1%worse than that of ab initio calculations as
the collapsed lone pairs. In the case of the 6-31G In the case of H2C0 - - - HOH, we examined
and 6-31G** basis sets, we found no such col- the angular dependence of the interaction energy
lapse, but there still was significant sensitivity of for HOH approach in the H2C0 plane, as well as
the charges to lone-pair placement. In contrast the “ T ” approach of the OH bond toward the
to the charges derived without lone pairs, in carbonyl oxygen. As one can see, both models
which the 6-31G* and 6-31G** basis sets results reproduce the lone-pair directionality and prefer-
were nearly identical, there are significant dif- ence of H-0 approach in the plane, although
ferences between the charges derived from these the model including lone pairs has a slightly
two basis sets when lone pairs were included. larger minimum energy 8. The biggest difference
These differences are reduced when dipole and between the models comes in the calculated rela-
quadrupole constraints are applied to the calcu- tive stabilities of T and a H bonds, where the
lation. When one turns to 5-point charge models atom-centered model finds a much smaller en-
(no charge on oxygen) (Table X), we find physi- ergy difference between H bonds than the ex-
cally reasonable charges, which, for the 6-31G* tended model. Ab initio calculations using the
and 6-31G** basis sets, give essentially identical 6-31G** basis set2’ on a number of complexes
quality fits to the 6-point charge model. How-
ever, with sm-3G and 4-31G basis sets,the extra
flexibility in the charge model does improve the
fit. Only in the case of s ~ o - 3 Gdo the charges on Table XVIII. C h a r m for 9-CH3 manine.’
the carbon and hydrogens differ greatly in the property STO-3G-M
Baslm SetIModel
STO-U: UA CLEU-M CLELI-UA
two models. In none of these models did the lone
pairs “invert”; in all models these lone pair ----
Charecs
NI -0.729 -0.746 -0.881 -0.952
C2 0.871 0.842 0.988 0.972
charges remained in intuitively reasonable posi- N3 -0.709 -0.702 -0.726 -0.712
c4 0.391 0.490 0.333 0.504
tions. In Table X, we present the results of the c5
C6
-0.060
0.690
-0.088
0.714
0.003
0.771
-0.088
0.865
above model calculations using the 6-31G** ba- N7
C8
-0.543
0.266
-0.575
0.428
-0.615
0.157
-0.703
0.570
sis set. N9
N2
-0.022
-0.778
-0.360
-0.758
0.128
-0.998
-0.466
-0.979
Given that one has these charge models, can 06
H1
-0.458
0.336
-0.459
0.340
-0.492
0.420
-0.500
0.432
H8
they be used to study intermolecular interac- c9n
0.046
-0.459 0.216
0.132
-0.747 0.263
W12A 0.339 0.333 0.414 0.405
tions? mN2B 0.325 0.324 0.392 0.394
HC
.9A. ~ ~ n..
. i.
f.
i3
. 0.239
. ...
To carry out such studies, one includes repul-
~

mc9n 0.164 0.241

mcgc 0.164 0.241
sion and dispersion parameters in addition to the “x -3.416
5.357
-3.471
5.294
-3.187
5.784
-3.293
5.699
electrostatic models. One can evaluate the en- 3
Qb
25.931
-19.728
26.062
-19.984
29.991
-2t.831
29.566
-21 -784
ergy at specific geometries, as well as carrying -6.203
0.89
-6.078
1.22
-8.160
1.22
-7.782
1.90
5.6 7.7 7.1 11.2
out energy refinement. In these studies, we used
the AMBER molecule mechanics p r ~ g r m . We ~* ‘See Table XI1 for notation.
Electrostatic Charges for Molecules 139

Figure 4. (a) Difference between united-atom model and quantum mechanical potential for 9-CH3.guanine. See
Figure 2(a) for color code. @) Difference between all-atom model and quantum mechanical potential for 9-CH3
guanine. See Figure 2(a) for color code.

including fonnaledhyde (Table XI) suggest that dimethyl phosphate (Tables XII-XXIV). We ex-
the extended model more accurately represents amined both the N-CH, bases of most direct
the relative energies of u and R H bonds. relevance to DNA and RNA, as well as the bases
themselves. Because of the size of these mole-
cules, we used only the more limited s~o-3Gand
NUCLEIC ACID BASES, SUGARS,
the minimal atomic basis of Clementi in these
AND PHOSPHATES
studies. The s~o-3Gand Clementi basis sets give
We next determined charge models for the rather different results for H,O, with the latter
nucleic acid bases, deoxyribose and ribose, and yielding multiple moments more characteristic
140 Singh and Kollman

Table XIX. Charges for manine.' XXV. With some of the charge models, we added
Property STO-3G-AA
Basis SetlModel
STO-3G-UA CLEM-AA CLEM-UA
10-12 H-bond potentials with well depths of 0.5
Charges kcal/mol and varying R* (minimum-energy dis-
N1 -0.691 -0.716 -0.819 -0.893
c2 0.783 0.797 0.898 0.940 tance), in order to ensure reasonable H-bond
length (H * - 0 and H * - N) of 1.8-1.9 A. As
N3 -0.723 -0.721 -0.768 -0.761
c4 0.540 0.557 0.567 0.616
c5 -0.010 -0.039 0.110 0.026
C6
N7
0.655
-0.562
0.685
-0.586
0.683
-0.653
0.774
-0.725
one can see, the agreement of our models with
0.640
C8
N9
0.357
-0.517
0.463
-0.585
0.325
-0.646 -0.848 the best available experimental AH values is
N2 -0.720 -0.725 -0.953 -0.969
06 -0.456 -0.458 -0.477
0.416
-0.487
0.426
very good and better than all previous point-
H1 0.335 0.339
H8
H5
0.037
0.328 0.345
0.112
0.408 0.459
charge potentials. Also interesting is the ratio of
HNZA
HN20
0.331
0.312
0.330
0.314
0.415
0.382
0.412
0.390 Watson-Crick H-bonding energies between AT
-3.419 -3.459
Px
pz
-3.569
-4.995
-3.583
L.992 5.184 5.173 and GC; the calculations using the electrostatic
20.636 19.787 22.791 21.468
2
Qc
-14.430
-6.206
- 13.8 12
-5.975
-14.505
-8.286
-13.699
-7.769
potential derived charges and our earlier Mulli-
0 0.92 0.95 1.26 1.46 ken charges32closely reproduce the experimental
5.4 5.7 7.0 8. I
ratio of -
6076, whereas the H-bond energies
'See Table XI1 for notation. calculated using the charges reported by Orn-
stein and Rein,34 Renugopalakrichnan,
Sakshminarayanan, and Sasi~ekharan,~~ and
Pullman and Pullman33 are not as accurate in
of more extended basis sets (Table I), so it is
encouraging that the multipole moments of the
bases are calculated to be similar for the STo-3G Table XX. Charges for 9-CH3 adenine.a
and Clementi basis sets. property
STO-3G-M
Basis SetlMadel
STO-30UA CLEM-M CLEM UA
The partial charges for the purines with N9-H
compared with N9-methyl are quite similar ex- Charges
81 -0.774 -0.760 -0.817 -0.859
cept for the charges at these groups and the C2
N3
0.661
-0.728
0.371
-0.717
0.607
-0.787
0.700
-0.776
neighboring C8 and C4 atoms. Similarly, for the c4
c5
0.346
-0.097
0.695
-0.151
0.587
-0.019
0.580
-0.104
C6 0.769 0.813 0.855 0.933
pyrimidines, the N1-H and N1-CH, and their N7 -0.343 -0.599 -0.609 -0.713
C8 0.263 0.488 0.180 0.575
neighbors C2 and C6 differ most between the N9 -0.063 -0.451 -0.191 -0.340
N6 -0.768 -0.793 -0.971 -1.004
bases and their N1-CH, analogs. H7
...
H8
-0.032
0.062
0.043
0.143
Because our current molecular mechanical C9M
HN6A
-0.431
0.333
0.230
0.335
0.298
0.410
0.185
0.411
simulations on proteins and nucleic acids use a HN68
HC9A
0.324
0.158
0.340 0.386
-0.043
0.411

HC9B 0.159 -0.037

united-atom (UA) approach for aliphalic and HC9C 0.159
~ ~ . . -0.037
2.274 2.227 2.357 2.358
aromatic CH, groups, we have also determined a 0.264 0.204 0.302 0.321
15.185 15.294 15.229 13.643
united-atom model for the charges of the bases -10.406
-4.779
-10.512
-4.782
-9.521
-5.708
-10.304
-5.339
and their N-CH, analogs as well as the all-atom 8
0.79
8.8
1.23
11.05
1.20
12.1
1.50
14.8
(AA) models. This is an alternative to using the
all-atom model, summing up the charges on the 'See Table XI1 for notation.
C-H hydrogens and placing the net charges on
the carbon. We compared the fit to the potential
using these two approaches for 1-CH, cytosine Table XXI. Charges for adenine?
(Table XII), and it is clear that the united-atom Bash setlnodel
Property STO-3G-M STO-3GUA CLEM-M CLEM-UA
model is superior to the all-atom model in which Charges
the C-H charges are summed up (in parenthe- Nl
C2
-0.769
0.648
-0.768
0.592
-0.824
0.607
-0.860
0.720
ses in Table XII). Figure 1 shows the quantum N3
c4
-0.757
0.680
-0.750
0.746
-0.828
0.744
-0.830
0.763
mechanically calculated potential for 1-CH, cy- c5
C6
-0.022
0.728
-0.081
0.793
0.078
0.829
-0.002
0.857
N7 -0.574 -0.616 -0.656 -0.740
tosine, and Figures 2(a) and 2(b) show the dif- C8 0.373 0.528 0.329 0.687
N9 -0.568 -0.685 -0.693 -0.897
ference between the quantum mechanically N6 -0.777 -0.819 -1.008 -0.994
82 -0.023 0.051
calculated potential and the united atom [Fig. H8
H9
0.049
0.338 0.363
0.128
0.415 0.467
2(a)] and all-atom [Fig. 2(b)] models. Figures 3 HN6A
HN6B
0.345
0.328
0.349
0.346
0.428
0.400
0.420
0.407
and 4 show corresponding potentials for 9-CH3 u, 2.166
-0.147
2.140
-0.142
2.298
0.207
2.294
0.183
guanine, and Figures 5 and 6 show corresponding 2 9.772
-5.200
11.461
-6.936
12.732
-5.809
13.862
-7.556
-4.572 -4.525 -6.922 -6.306
potentials for 1-NH, deoxyribose. Qc
0.98 1.04 1.31 1.60
B 9.9 10.5 11.8 14.5
The results of Watson-Crick H-bond calcula-
tions using these charges are presented in Table 'See Table XI1 for notation.
Electrostatic Charges for Molecules 141

Figure 5. Quantum mechanical potential for 1-NH, deoxyribose. See Figure 1 for color code.

this regard. The charges reported by Ornstein heavily closer atom-atom interactions. The re-
and Rein,34and Renugopalakrichnan, Sakshmin- sults in Table XXV show that the H-bond en-
arayanan, and S a s i ~ e k h a r a ndo
~ ~not reproduce ergies we calculate are in satisfactory agreement
the dipole moments28 as well for the isolated with experiment with either model.
monomers,38but the charges reported in ref. 35
give qualitatively reasonable H-bond energies by
comparison with experiment, although the GC H
bond energy is significantly too strong in com-
parison with the AT structure. Table =I. Charges of 1-NH, ribose."
Langlet, Claverie, and C a r ~ nhave ~ ~ carried Basis SetlModel

out more precise calculations of the Watson- property STD-3G AA STO-36 UA CLEM-A4 CLEK-UA
Charges
Crick H bonding in nucleic acid bases, and the 01' -0.343 -0.413 -0.452 -0.503
CI' 0.373 0.537 0.543 0.726
agreement between our simple model and theirs C4'
C2'
0.100
0.101
0.197
0.083
0.098
0.181
0.276
0.094
is encouraging; both give quite good agreement C3'
N
0.303
-0.828
0.271
-0.894
0.446
-1.074
0.315
-1.145
with experiment. Given their finding that the 03'
C5'
-0.520
0.180
-0.542
0.175
-0.683
0.271
-0.683
0.174
Hoogsteen base pair of AT is slightly more stable HI'
HNA
0.054
0.312 0.325
0.057
0.372 0.386
HNB 0.308 0.317 0.376 0.387
than Watson-Crick, we examined this dimer as H2' 0.008 -0.027
02' -0.546 -0.512 -0.682 -0.669
well and found results qualitatively similar to H3' 0.007 -0.037
H4' 0.061 0.043
those of ref. 35 (see parentheses in Table XXV). H03'
05'
0.303
-0.463
0.316
-0.488
0.391
-0.596
0.410
-0.600
In Table XXV we have carried out the cal- H5 '
H5'
0.001
0.016
-0.036
-0.006
culations with E = l and a distance-dependent 805'
H02 '
0.289
0.324
0.300
0.329
0.379
0.436
0.390
0.442
dielectric E = R . Elsewhere we discuss37why we 1.723
-1.740
1.749
-1.754
2.452
-1.917
2.497
-1.930
0.836 0.848 0.863 0.868
prefer the latter approach for solution-phase 5.375 4.949 7.536 7.254
-3.658 -3.253 -5.344 -5.128
simulations, but it is clear that the former model -1.717 -1.696 -2.192
1.02
-2.126
1.10
0.69 0.75
is more appropriate for comparison with the 6.9 7.6 8.2 8.9

gas-phase data of Yanson, Teplitsky, and 'Notation as in Table XII; AMBER (ref. 24) optimized
Sukh0dub.3~On the other hand, the E = R model geometry was used, with C(3') endo sugar pucker; 01-C1'
bond is along the 2 axis and Ol'-Cl-N are in the X Z
does indirectly include polarization effects for +
plane, with N in the X direction. We used a model where
attractive interactions in that it weighs more the base is replaced by -NH,.
142 Singh and Kollman

Table XXIII. Charges of 1-NH, deoxy ribose." We have also compared, in the case of l-NH,
Property STO-3G-M
Basis SetlXodel
STO-3GUA CLEM-M CLEM-UA
deoxyribose, the difference in charges when one
uses a C3' endo rather than a C2' endo geometry
Charges
01 -0.368 (-0.371)b -0.386 -0.494 -0.475 (Table XXIII, all-atom approximation). The H-
CI' 0.558 ( 0.470) 0.502 0.861 0.730
C4' 0.036 ( 0.029) 0.185 0.122 0.253 bonding-site charges change rather little, with
C2' -0.307 (-0.351) -0.047 -0.373 -0.103
C3' 0.233 ( 0.417) 0.172 0.276 0.190 the exception of 05', which changes from - 0.404
N -0.869 (-0.811) -0.876 -1.172 -1.151
03'
C5'
-0.508
0.118
(-0.541)
( 0.171)
-0.514
0.153
-0.665
0.206
-0.666
0.144
to -0.477. However, the atoms to which it is
HI' 0.009 ( 0.017) - -0.015 bonded, C5' and H05', become 0.070 more posi-
HH(\ 0.305 ( 0.298) 0.313 0.374 0.376
HHB
H2'
0.329
0.081
( 0.309)
( 0.071)
0.327
-
0.415
0.071
0.409
tive, almost exactly compensating for this charge
H2' 0.081 ( 0.102) - 0.080
H3' 0.025 (-0.001) -- 0.016 shift. There are larger polarizations in the C2',
84' 0.056 ( 0.047) 0.043
H03' 0.306 ( 0.306) 0.313 0.412 0.423 C3', C4' charges, but since these are not likely to
05' -0.404 (-0.477) -0.425 -0.488 -0.489
H5'
H5'
0.003
0.039
(-0.004)
( 0.023)
-
-
-0.039
0.017
form strong electrostatic interactions, it is not
H05' 0.279
1.783
( 0.296) 0.285
1.799
0.354
2.389
0.360
2.426
clear how significant these shifts are. More stud-
Ux
-1.619 -1.619 -1.899 -1.901 ies are required to more precisely delineate the
2 0.215
8.843
0.234
9.081
0.121
9.957
0.131
10.271
validity of the approximation of using a set of
9 -1.594
-7.249
-1.657
-7.423
-1.372
-8.585
-1.508
-8.763
charges derived from one conformation to repre-
$ 0.68 0.85 1.08 1.25
7.4 9.4 9.8 11.2
#
sent others.
'Notation as in Table XI1 and XXII, except sugar The question of whether one should include
pucker is C(2') endo. lone pairs in partial-charge representations of
bCharges derived using a C(3') endo conformation.
molecules has been addressed here, but no single
definitive answer to this question has been given.
Obviously, inclusion of lone pairs improves the
DISCUSSIONS AND CONCLUSIONS fit of the model to the electrostatic potential.
We have presented applications on a variety of In the case of (H,O),, the directionality of
molecules of a method for generating point- the hydrogen bond, the angle the OH bond of
charge models for molecules, which can then be the proton donor makes with the bisector of the
used in simulations. It is clear from the results of HOH angle of the acceptor, goes from 8 = 25 to
Watson-Crick hydrogen-bonding studies that 8 = 65 upon addition of lone pairs, and the latter
this method, when used with a suitable ab initio is closer to the experimental angle. The lone
basis set, can give very good agreement with pairs in the optimum model, interestingly, are
experimental enthalpies of H-bond formation for inverted from the tetrahedral direction and point
nucleic acid base pairs. It is superior to Mulliken in the direction of the hydrogens. This seemingly
populations in this regard and is, we believe, a strange result is consistent with point-charge
preferred method of approach for generating such models for water based on ab initio dimer
charges. However, even given a set of atomic surfaces39 and with such models designed to fit
charges that reproduce well the molecular elec- water liquid properties.4o Such models, unlike
trostatic potentials, one still faces the task of those in which the lone pairs are in the "tradi-
determining a set of nonbonded parameters that tional" direction, can given excellent agreement
can be used with the charges to give good geome-
tries and energies for complex formation. It
should also be kept in mind that in most typical Table XXIV. Charees for dimethvl DhosDhate.'
simulations, the intermolecular potential con- Basis SetlUodel
tainselectrostatic, dispersion, and repulsion com- STO-3G-M SIO-3GUA STO-3G*-Mb SX-3G*-UA

ponents, but no polarization or charge transfer. r

h.-.-.
".."_~_s
P 1.560 1.596 0.867 0.913
-0.509 -0.534 -0.377 -0.409
Thus the parameters used must compensate for OCA
OCB -0.508
0.032
-0.535
0.086
-0.375
0.034
-0.410
0.108
CA
the absence of some components. For example, CB
OPA
0.037
-0.846
0.088
-0.849
0.037
-0.650
0.110
-0.655
one often u6es electrostatic charges that lead to OPE
HCAI
-0.848
0.003
-0.851 -0.652
0.005
-0.657

larger dipole moments than the gas-phase value HCA2

.._.._
"C.7
0.025
-.- - .
n.mL
0.032
_..-.
0.023
HCBl 0.002 0.004
for the molecule, in order to simulate polariza- HCB2 0.023 0.030
HCB3 0.015 0.024
tion effects. In any case,the charges and non- 0
*
1.12
1.2
1.17
1.45
1.08
1.2
1.17
1.3
bonded parameters must be tested in some model
cases to ensure that they lead to results in good 'Notation as in Table XII; geometry from ref. 31.
bsTo-3G* basis set was used;see ref. 32.
agreement with experiment or accurate calcula-
tions, and we have tried to do that here and
elsewhere.
Electrostatic Charges for Molecules 143

Figure 6. (a) Difference between united-atom model and quantum mechanical potential for 1-NH, deoxyribose. Color
code as in Figure 2(a). (b) Difference between all-atom model and quantum mechanical potential for 1-NH, deoxyribose.
Color code as in Figure 2(a).

with both the dipole and quadrupole moments lone-pair-containingpartial-chargemodels. These

for water. charges may collapse into the atom,extend too
In the case of CH,OH, (CH,),O, and CH,O, far from the atom, or become unrealistically large
we have found numerical problems in deriving in magnitude, unless appropriate restraints are
144 Singh and Kollman

Table XXV. Hydrogen-bond energies calculated with different charge models.

Model AE(AT)~ R(H.. .B) AT^ AE(GC)~ R(R.. .B) ,GCg

Chargesa -
cb -
R*‘

STO-UA~ 1 2.0 11.3 (11.8)g 1.92.1.86 21.2 1.87,1.90.1.87

STO-NUAI 1 2.0 11.9(13.5)g 1.91,1.91 21.1 1.87.1.90,1.87
STO-IJA~ R 2.0 13.4 1.86,1.82 21.7 1.84.1.86.1.83
STO-WAi R 2.0 13.7 1.87,1.81 22.1 1.84.1.84.1.83
Clelnl 1 2.1 13.1 1.91.1.98 22.6 1.93,1.95,1.93
Clem R 2.1 15.2 1.87,1.93 23.2 1.90.1.90,1.90
Pullman‘ 1 2.0 8.0 1.93,1.94 18.6 1.90,1.93,1.89
Fulllpan R 2.0 8.5 1.90.1.91 19.9 1.88,1.90,1.86
Rein 1 - 5.3 1.89,1.92 11.1 1.85.1.90.1.81
Rein R - 5.6 1.82,1.88 12.0 1.81.1.84.1.72
Sasim 1 2.0 9.1 1.91,1.91 23.6 1.87,1.91,1.86
Sasi R 2.0 10.0 1.88,1.89 25.8 1.85.1.89.1.84
pAKn 1 - 6.8 1.86,1.89 13.1 1.82,1.88,1.81
PAK R - 6.6 1.78,1.84 11.7 1.77,1.81,1.73
Claver i e ” 12.9 (13.6)g - 23.7 -
Yanson (exp)P 13.0 - 21.0 -
~

aSet of charges used.

bDielectric model.
9&12 potential minimum energy R; all well depths 0.5 kcal/mol.
dEnergy of AT Watson-Crick H bond in kcal/mgl.
“H6 . . . 0 4 Hl-H3 H-bond distances in AT in A.
fEnergy of GC hydrogen bond in kcal/mol,
g06 ... H4, H1 . .. H3, and H2 ... 02 H-bond distances in GC (A).
hThis study; united atom charges sTo-3G basis set.
’This study; all atom charges; sT0-3G basis set.
jThis study; Clementi basis set (ref. 17).
kCharges from Pullman and Pullman (ref. 33).
’Charges from Omstein and Rein (ref. 34).
mCharges from Renugopalakrichnan, Sakshminarayanan, and Sasisekhraran (ref. 35).
”Charges from Singh and Kollman (ref. 29) based on Mulliken population.
“Semiempiricalanalysis of Langlet, Clavene, and Caron (ref. 36).
PExperimental data of Yanson, Teplitsky, and Sukhodub (ref. 37).
qcalculated energy for Hoogsteen base pair.

placed on the fitting procedure. We note that our generating such monopoles. This approach can
software allows such restraints to be imposed. be employed for either atom-centeredor lone-pair
Often, the “physically realistic” set of charges models and for either all-atom or united-atom
differs from those that are not so only slightly in models. It is a very simple generalization to
quality of fit to the quantum mechanical poten- determine charge models with higher multipo1.es
tial. as well. Elsewhere we present a complete set of
For H2C0 we have also compared the lone-pair charges appropriate for simulations of proteins
models with those without lone pairs for H2C0 and nucleic acids.41
. - HOH. Again, the qualitative features do not
+

change upon inclusion of lone pairs, but the P.A.K. is pleased to acknowledge research support by
the NSF (CHE-80-26560) and NIH (CA-25644). We are
angle the C-0 bond makes with the HO bond also very grateful to F. van Duijnveldt for his program to
increases slightly from 50 to 60 O .
O carry out the least square fitting to quantum mechanical
In summary, there is no compelling qualitative potentials and to Kitty Ghio for the PISA electrostatic
potential program.
reason to include )one pairs in simulations in-
cluding N and 0 lone pairs (? stronger case can
be made concerning the inclusion of lone pairs
for second-row atoms, such as S).41 In the long References
term,as more accurate models develop, lone pairs 1. N. Allinger, J . Am. Chem. Soc., 99, 8127 (1977); 0.
or higher multipoles than monopoles in the elec- Ermer and S. Lifson, zbid., 95, 4121 (1973).
trostatic expansion can be employed, but for the 2. A. Pullman, C. Zakrewska, and D. Perahia, Int. J.
Quantum Chem., 16,393 (1979).
current status of simulations, atom-centered 3. L. Dunfield, A. Burgess, and H. Scheraga, J . Phys.
point monopoles seem reasonable and capable of Chem., 82, 2609 (1978).
reproducing well many molecular properties and, 4. D. M. Hayes and P. Kollman, J . Am. Chem. Soc., 98,
3335, 7811 (1976).
as such, should be useful. In this article we have 5. E. Scrocco and J. Tomasi, Ado. Quantum Chem., 11,
presented a simple and systematic approach to 115 (1978).
Electrostatic Charges for Molecules 145

6. P. H. Smit, J. L. Derissen, and F, B. van Duijneveldt, 25. P. Kollman, in Modern Theoretical Chemistry, H. F.
Mol. Phys., 37, 521 (1979). Schaefer, Ed., Plenum, New York, 1977, Vol. 4.
7. S. R. Cox and D. E. Williams, J. Comput. Chem., 2, 26. T. R. Dyke and J. Muenter, J. Chem. Phys., 60,2929
304 (1981). (1974).
8. U. C. Singh and P. Kollman, “GAUSSIAN 80 UCSF,” 27. H. Popkie, H. Kistenmacher, and E. Clementi, J.
Quantum Chemistry Exchange Program Bull., 2, 17 Chem. Phys., 69, 1325 (1973).
(1982). 28. D. Hankins, J. Moskowitz, and F. Stillinger,J. Chem.
9. J. S. Binkley, R. A. Whiteside, R. K i r b a n , R. Phys., 63, 4544 (1970).
Seeger, D. J. Defrees, H. B. Schlegel, S. Topiol, L. R. 29. P. Kollman, P. Weiner, and A. Dearing, Bwpolymers,
Kahn, and J. A. Pople, “GAUSSIAN 80,” Quantum 20, 2583 (1981).
Chemistry Program Exchange, (1980). 30. V. Bloomfield, D. Crothers, and I. Tinoco, Physical
10. M. Dupuis, D. Spangler, and J. J, Wendoloski, Chemistry of Nucleic Acids, Harper and Row, New
“GAMESS,” NRCC Software Catalog, 1980, Program York, 1974.
No. QGO1, Vol. I. 31. M. Newton, J. Am. Chem. SOC.,96,256 (1973).
11. P. Marsh and D. E. Williams, “GAUSSIAN 79,” 32. J. B. Collins, P. Schleyer, J. Binkley, and J. Pople, J.
Quantum Chemistry Program Exchange, (1980). Chem. Phys., 64, 5142 (1976).
12. Subroutine courtesy of R. Bonaccosi from the Uni- 33. B. Pullman and A. Pullman, Prog. Nuckic Acid Res.
versity of Pisa. Mol. Biol., 9, 327 (1969).
13. M. Connolly, Quantum Chemistry Program Ex- 34. R. Ornstein and R. Rein, Bipolymers, 18, 1277, 2821
change, (1982). (1979).
14. R. Fletcher, Practical Methals of Optimization, 35. V. bnugopalakrichnan, A. Sakshminarayanan, and
Wiley, New York, 1980, Vol. 1, Chap. 5. V. Sasisekhraran, Biopolymers, 10,1159 (1971).
15. K. Kitaura and K. Morokuma, Int. J. Quantum 36. J. Langlet, P. Claverie, and F. Caron, in Zntennokcu-
Chem., 10,325 (1976). Zur Forces, €’roc. 14th Jerusalem Symp., B. Pullman,
16. W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Ed., Reidel, Dordrecht 1981.
C h . Phys., 61,2657 (1969). 37. I. Yanson, A. Teplitsky, and L. Sukhodub, Biopoly-
17. E. Clementi, J. M. Andre, M. C. Andre, D. Klint, and mers, 18,1149 (1979).
D. Hahn, Acta Phys. Cherrz., 27,493 (1969). 38. The calculated dipole moments in debyes for the
18. W. J. Hehre, R. F. Stewart, and J. A. Pople, J. N-CH, bases are for A,T, C, G: 2.2,3.3,5.4,6.3 [this
Chem. Phys., 61,2657 (1969). study, UA model s ~ o - 3 Gbasis (ref. 16)]; 2.4,
19. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta., 3.7,5.7,6.6 [this study UA model Clementi basis (ref.
17)]; 1.8,2.2,4.4,4.7 (ref. 32); 1.1,3.4,4.4,5.4 (ref. 33);
2.6,3.8,8.1,7.8 (ref. 34); 3.0,3.9,6.5,7.2 (ref. 35). Indi-
Chem. Phy;, 63;2823 (1970). rect experimental data in ref. 30 suggests that the
For H (6s2p/2slp), E. Clementi and H. Popkie, J. dipole moments for 9-CH, adenine and 1-CH,
Chem. Phys., 67, 1077 (1972). -
thymine are 3.0 and 3.9 D, respectively.
39. H. Popkie, H. Kistenmacher, and E. Clementi, J.
21. W. S. Benedict and L. D. Kaplan, J. Chem. Phys.,
30,388 (1959). Chem. Phys., 69, 1325 (1973).
22. W. H. Flygare, J.Am. Chem. SOC.,94, 330 (1972). 40. W. Jorgensen, J . Chem. Phys. Lett.,77, 4156 (1982).
23. W. Flygare and R. C. Benson, Mol. Phys., 20, 225 41. S. J. Weiner, U. C. Singh, K. Ghio, G. Alagona, D.
(1971). Case, P. Weiner and P. Kollman, J. Amer. Chem.
24. P. Weher and P. Kollman, J. Comput. Chem., 2,287 SOC.(In prw).
(1981).