3358
J . Phys. Chem. 1991,95, 3358-3363
Charge Equilibration for Molecular Dynamics Simulations Anthony K. Rap++ and William A. Goddard III*,l BioDesign, Inc., Pasadena, California 91 101, Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523, and Materials Simulation Center, Beckman Institute (139-24),1 California Institute of Technology, Pasadena, California 91 I25 (Received: October 4,1989)
We report here an approach for predicting charge distributions in molecules for use in molecular dynamics simulations. The input data are experimental atomic ionization potentials, electron affinities, and atomic radii. An atomic chemical potential is constructed by using these quantities plus shielded electrostatic interactions between all charges. Requiring equal chemical potentials leads to equilibrium charges that depend upon geometry. This charge equilibration (QEq) approach leads to charges in excellent agreement with experimental dipole moments and with the atomic charges obtained from the electrostatic potentials of accurate ab initio calculations. QEq can be used to predict charges for any polymer, ceramic, semiconductor,or biological system, allowing extension of molecular dynamics studies to broad classes of new systems. The charges depend upon environment and change during molecular dynamics calculations. We indicate how this approach can also be used to predict infrared intensities, dielectric constants, and other charge-related properties.
1. Introduction Knowledge of the charge distribution within molecules is essential for determining the electrostatic energies (including hydrogen bonding) in molecular mechanics and molecular dynamics calculations.I4 Unfortunately, reliable charge distributions are known only for a few organic molecule^.^^^ Thus, currently there is no effective approach to estimate the charges for inorganic systems (ceramics, zeolites, high- T, superconductors), and current estimates of charges for polymers and large organic systems are quite uncertain. For biological molecules, the 20 standard amino acids and four standard bases have been assigned charges2+ that are expected to be reasonably accurate; however, charges are not available for nonstandard amino acids, unusual bases, and various cofactors and substrates. An additional serious problem is that current approachesI4 to molecular mechanics and molecular dynamics use fixed charges that cannot readjust to match the electrostatic environment. Since the charges are not allowed to respond to the environment, the tradition is to incorporate a dielectric constant in the interaction potential, leading to additional uncertainties in the calculations. We propose here a general scheme for predicting charges of large molecules based only on geometry and experimental atomic properties. The charge equilibration (QEq) approach allows the charges to respond to changes in the environment, including those in applied fields, and can be applied to any material (polymer, ceramic, semiconductor, biological, metallic). In section 11, we derive the basic equations for the charge equilibration approach. The scaling parameter X relating atom size to crystal atomic radii is determined in section 111 by comparing theory and experiment for the alkali-metal halide diatomic molecules. In section IV, we discuss hydrogen atoms, which require an extension of the simple scheme of section 11. Finally, in section V we apply the QEq method to a number of molecules and compare our results with experiment or ab initio theory. The concepts involved in the QEq approach rest upon earlier ideas of Pauling, Mulliken, Margrave, Parr, Pearson, Mortier, and others. Section VI summarizes the relationship between QEq and some of these earlier ideas. In section VII, we mention some possible extensions utilizing the ability of QEq to allow polarization of the charge distribution. 11. Charge Equilibration
A. Charge Dependence of Atomic Energy. In order to estimate the equilibrium charges in a molecule, we first consider how the energy of an isolated atom changes as a function of charge. Using a neutral reference point, we can write the energy of atom A as7
’* BioDcsign, BioDesign, Inc., and Colorado State University. Inc., and California Institute of Technology. (Contribution No. 8340.
Including only terms through second order in (1) leads to
so that
(
$)Ao
= y2(1P + EA) = xi (3)
where IP and EA denote the ionization potential and electron affinity and xA is referred to as the electronegatiuiry. To understand the physical significance of the second-derivative quantity a*E/dQZ,consider the simple case of a neutral atom with a singly occupied orbital, 4A, that is empty for the positive ion and doubly occupied for the negative ion. The difference between the IP and EA for this system is
IP - EA = J o u
(4)
where s),is the Coulomb repulsion between two electrons in the $A orbital (the self-Coulomb integral). We refer to this atomic repulsion quantity pAA as the idempotential (self-Coulomb) for less awkward reference to it in later discussions. Of course, the ( I ) Williams, D. E.; Cox, S.R. Acta Crysrallogr., Secr. B 1984,40,404. Williams, D. E.;Houpt, D. J. Ibid. 1986,42,286. Williams, D. E.;Hsu, L. Y. Acra Crysrallogr.,Sect. A 1985.41,296. Cox, S.R.; Hsu, L. Y.;Williams, D.E. Ibid. 1981, 37, 293. (2)Weiner, S.J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.;Weiner, P. J. Am. Chem. SOC.1984,106,765-784. Weiner, S.J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Compur. Chem. 1986,7,230-252. (3)Brooks, R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.;Karplus, M. J . Compur. Chem. 1983, 4, 187. (4)Jorgensen, W. J.; TiradeRives, J. J . Am. Chem. Soc. 1988,110, 1657. (5) Cox, S.T.; Williams, D. R. J . Compuf. Chem. 1981,2, 304-323. (6)Chirlian, L. E.;Francl, M. M. J . Compuf.Chem. 1987,8, 894-905. (7) Iczkowsky,R. P.; Margrave, J. L. J. Am. Chem. Soc. 1%1,83, 3547. (8) Parr, R. G.; Pearson, R. G.J. Am. Chem. Soc. 1983,I05,1503-1509.
0022-3654/91/2095-3358.$02.50/00 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 8, 1991 3359
Charge Equilibration for Molecular Dynamics Simulations TABLE I: Atomic Parameters" element x . eV J . eV Li 3.006 4.772 C N 0
5.343 6.899 8.741 10.874 2.843 4.168 5.463 6.928 8.564 2.421 7.790 2.331 6.822 2.183 4.5280b
F Na Si
P S CI K Br
Rb I
cs H
10.126 1 1.760 13.364 14.948 4.592 6.974 8.000 8.972 9.892 3.84 8.850 3.692 7.524 3.422 13.8904b
R. A
L au
1.557 0.759 0.715 0.669 0.706
0.4174 0.8563 0.9089 0.9745 0.9206 0.4364 0.7737 0.8257 0.8690 0.9154 0.4524 1.0253 0.5162 1.0726 0.5663 1.0698
2.085 1.176 1.102 1.047 0.994 2.586 1.141 2.770 1.333 2.984 0.371
-
or XA(QI.-QN)
+'
/ ~ ~ A Q A ~
= 14.4/RoA or
RoA
(6)
UEAO + xOAQA) A
+
=
XN
(8)
=
CQi (=I
(9)
CD = -D
(10)
4 = -Qmt Di = xp - xp for i 1 2
(1 1)
and Cii =
Qi
Cij = Jij- J i j for i 1 2
(12)
The inequalities in (5) are implemented in our programs as follows. We first solve (lo)-( 12) for the charges and check the inequalities in ( 5 ) . If any atom is outside its range, we fix its charge at the boundary. Defining D for the nonfixed atoms as Di = xYF - xyF for i # 1
where
we solve the reduced set of equations. We find that this procedure works reliably for all cases considered. C. Shielding Corrections. In order to solve the QEq equations (lo), we must first specify the form for the Coulomb potential J A B between unit charges on centers A and B separated by a distance R . For large separations = 14.4/R (14) (where 14.4 converts units so that R is in angstroms and J is in electronvolts). However, for distances where the charge distributions on centers A and B overlap, the simple Coulomb law (14) 0, (14) leads to is no longer valid. Indeed, as R JAB(R)
JABW -
whereas it should lead to a finite value related to J A A and J B B , as illustrated in Figure 1. This overlap or shielding correction to (14) will be quite large for bonded atoms. There are a number of ways of evaluating the shielding of the two charge distributions. We have chosen to express the shielding as the Coulomb integral between atomic densities. We could obtain the atomic densities from accurate (spherically averaged) Hartree-Fock (HF) or local-density calculations on atoms. However, in the current implementation of QEq, we describe the atomic density in terms of a single Slater orbital. For an atom whose outer valence orbital is ns, np, or nd, we construct a normalized ns Slater orbital of the form
AFBQ~Q~J~~
which we rewrite as EQ(QI-.QN)
..e
where
(1')
etc. and take E A ( Q ) = m outside these ranges. B. Electrostatic Balance. In order to calculate the optimum charge distribution, we need to evaluate the interatomic electrostatic energy, C A < B Q A Q B J A B , where J A B is the Coulomb interaction between unit charges on centers A and B ( J A B depends on R A B , the distance between A and B). This leads to a total electrostatic energy of A
=
leads to a total of N simultaneous equations for the equilibrium self-consistent charges that are solved once for a given structure. These QEq equations can be written as
= 14.4/&
= UEAO + XIQA+ Y*QA~JO,A)+
B F A J ~ ~ Q(7')~
N Qtot
where the conversion factor 14.4 allows R i to be in angstroms and PAto be in electronvolts. This equation leads to PH = 0.84 A,$=1.42A,&= 1.22A,R$=1.08A,R~i=2.06A,@ = 1.60 A, and PLi= 3.01 A. Comparing with bond distances of = diatomics RLH = 0.74 A, R& = 1.23 A, RONN = 1.10 A, 1.21 A, = 2.20 A, R'& = 1.63 A, and PLiLi = 3.08 A, we see that this characteristic atomic distance corresponds roughly with the homopolar bond distance. Use of a quadratic relation such as (1') is expected to be valid only in a restricted region. In particular, the x and J are clearly invalid outside the range corresponding to emptying or filling the valence shell of electrons. Thus we restrict the ranges to -7 < Q L i < + I -4 < Qc < +4 -2 < Qo < +6 (5)
E(Qi-Qiv)
A Q A+
leads to an
Adding the condition on total charge
where the x i and &A can be derived directly from atomic data. However, the atomic IP and EA must be corrected for exchange interactions present in atoms but absent in molecule^.^ (The atomic states contain unpaired spins, whereas the molecules for which we will use x A and J A generally have all spins paired.) This leads9 to the generalized Mulliken-Pauling electronegativities and idempotentials in Table I. The idempotential is roughly proportional to the inverse size of the atom, and indeed, one can define a characteristic atomic size PAby &A
+~
xi = x2
optimum shape of the orbital changes upon adding an additional electron, and an accurate description of the electron affinity requires configuration interaction so that the &A derived from (4) may differ somewhat from the &A calculated with a Hartree-Fock wave function. Using (2) and (4) leads to ~ Q A
xOA
=
QA
where x A is a function of the charges on all the atoms. For equilibrium, we require that the atomic chemical potentials be equal, leading to N - 1 conditions
"Reference 9. bValues for QH = 0; see eqs 20 and 21.
EA(Q) = EAO +X
-
(suggesting that J A A ( R ) fiAas R 0). Taking the derivative of E with respect to atomic-scale chemical potential of the form
Y~CQAQBJAB (6') A,B
(9) Rap@, A. K.;Goddard, W. A., 111. Generalized Mulliken-Pauling Electronegativities. 1. Main Group Elements (Groups 1, 13-17). J . Phys. Chem., submitted for publication.
@$
= N,Fle-fr
(15)
3360 The Journal of Physical Chemistry, Vol. 95, No. 8, 1991
TABLE 11:
Rappi and Goddard
Charge Eauilibration Results
metal halide NaCl NaBr NaI KCI KBr KI RbCl RbBr Rbl CSCl CsBr CSI LIF LiCl LiBr Lil Na F KF RbF CsF
Q,Xp"
0.792 0.757 0.708 0.800 0.783 0.740 0.784 0.768 0.753 0.743 0.134 0.735 0.837 0.731 0.694 0.647 0.879 0.821 0.781 0.697
OReference 25. bFrom eq 17 with
QQM~
QA-0.5
0.766 0.745 0.709 0.775 0.768 0.754 0.763 0.757 0.747 0.769 0.767 0.763 0.791 0.939 0.902 0.841 0.665 0.662 0.653 0.655
0.776 0.756 0.720 0.784 0.777 0.764 0.771 0.766 0.757 0.777 0.776 0.773 0.803 0.958 0.921 0.860 0.671 0.667 0.657 0.660
€ = I 2.504 2.561 2.663 2.095 2.165 2.299 1.918 1.977 2.088 1.874 1.935 2.054 6.033
c=2 0.682 0.644 0.585 0.737 0.708 0.663 0.741 0.713 0.672 0.763 0.739 0.703 0.748 0.737 0.685 0.608 0.666 0.695 0.695 0.711
13.513 13.563 14.936 1.751 1.530 1.435 1.427
€=o 0.395 0.368 0.328 0.447 0.423 0.387 0.459 0.435 0.400 0.479 0.457 0.424 0.399 0.379 0.351 0.310 0.411 0.450 0.458 0.473
&,,, = 0.4913.
9.0 *a*,
14 0.420 0.391 0.350 0.473 0.448 0.411 0.484 0.460 0.424 0.505 0.482 0.449 0.427 0.406 0.377 0.333 0.434 0.473 0.481 0.496
c =
----14.4 I R
Shielded Potentials
We require that the calculated momentt0
1
~ M = X
QM
lead to the experimental dipole
( 1 /4.80324)QuR~x
(19)
where RMx is the experimental bond distance (the constant 4.803 24 allows Q to be in electron units, R in angstroms, and p in debyes). The only variable here is the scaling parameter A. The best value of X is 0.4913, which leads to an average error of 0.0018 e (see Table 11). Rounding off to X = 1/2 leads also to an average error of 0.0018 e, and hence (17) becomes {A
3.0
1
2 . 0 ~ ' " ' ' " ' 0.0 1 .o 2.0
3.0
Distance
4.0
5.0
(A)
Figure 1. Shielded potentials for Is-7s Slater orbitals. Here was (carbon). Also included is the taken from eq 17 with RA = 0.759~~ unshielded Coulomb potential, 14.4/R.
where N , is the normalization constant. From (15), the average size of the atom is
RA
(r) = (2n
+ 1)/(2{A)
(16)
Consequently, we choose the valence orbital exponent {A for atom A by the relation
= X(2n + I)/(%) (17) where RA is the covalent radius in atomic units (ao = 0.529 17 A) for atom A, which we select from experimental crystal structure data (seeTable I). An adjustable parameter X is included in (17) to account for the difference between an average atom size as given by ( I 6) and the crystal covalent radius RA. We require that the same X be used for all atom of the periodic table and in section I11 determine I by comparing the predicted and experimental dipole moments of the alkali-metal halide diatomics. The diatomic Coulomb integral JAB involving these Slater functions is evaluated exactly for {A and tBat the various distances. {A
111. Alkali-Metal Halides
In order to determine the scaling factor A that adjusts atomic radii to Coulomb shielding distance, we considered the 12 alkali-metal halide molecules MX, where M = Na, K, Rb, or Cs and X = CI, Br, or I. For these systems, (8) and (9) reduce to
Qx = -QM
= (2n
+ 1)/(4RA)
( 17')
(with RA in units of u,,). We did not use M = Li and X = F in the fits because the errors were larger for these first-row elements. However, the results for these eight cases are also listed in Table 11. Including these cases, the average error increases to 0.15 e.
IV. Hydrogen The Mulliken-like definition' * for electronegativity leads for hydrogen to x! = '/2(IP + EA) = 7.17 eV, which is not consistent with the Paulint2 or otherI3 empirical values for electronegativities. With xH,the hydrogen is more electronegative than C (x$J"P= 5.34) or N (x&*P = 6.90), whereas the Pauling scale (based on chemical experience) has h drogen much more el=tropositive than C (x; = 2.1, while xc = 2.5) and slightly more electronegative than boron (x', = 2.0). As discussed in ref 9, the problem with x! is that the effective EA for H is much smaller than the atomic value because the H orbital involved in a bond cannot expand to the value achieved in a free H- ion. Consequently, we redefine :x and JOHH for hydrogen, allowing EAH to be a variable. From an examination of the charges on H in the molecules LiH, CH4, NH3, H 2 0 , and HF, we find that an accurate description of Q H is obtained if the effective charge parameter tH is allowed to be charge-dependent:
2
(20) = ZOH + QH Here s", = 1.0698 is based on (17) where RH = 0.371 A. The idempotential JHH becomes charge-dependent: {H(QH)
JHH(QH)
= (1 + QH/&)J&H
(21)
(IO) Huber, K.; Herzberg, G. K. Cohtranrs of Diatomic Molecules; Van Nostrand-Reinhold Co.: New York, 1979. (11)Mulliken, R. S.J. Chem. Phys. 1935, 3, 573. (12) Pauling, L. Nature of rhe Chemical Bond, 3rd ed.;Cornell University Press: Ithaca, NY, 1960. ( I 3) Sanderson, S.T. Chemical Bonds and Bond Energy; Academic Press: New York, 1976.
Charge Equilibration for Molecular Dynamics Simulations
The Journal of Physical Chemistry, Vol. 95, No. 8, 1991 3361 0.19 (0.16)
-0.05(0.04)
0240.28
4.78
0.18
029
Figure 2. (a) Predicted charges for Ala-His-Ala. The N and 0 termini are charged as appropriate for a peptide. Comparisons with charges from AMBER^ are given in parentheses. (b) Same as (a) except that His is protonated.
TABLE 111: Charges on Hydrogen' compd HF HZO NH, CHI LiH
exptl 0.419 0.3258 0.2678 0.150" -0.76v
QEqb
0.462 0.353 0.243 0.149 -0.767
HFC 0.462 0.398 0.338 0.124 -0.682c
-0.44 o,1
