- Email: [email protected]
Charge distributions in water and ion–water clusters
23 February 2001
Chemical Physics Letters 335 (2001) 327±333
www.elsevier.nl/locate/cplett
Charge distributions in water and ion±water clusters A.A. Rashin a,*, I.A. Topol b, G.J. Tawa b, S.K. Burt b b
a BioChemComp Inc., 543 Sagamore Ave, Teaneck, NJ 07666, USA Advanced Biomedical Computing Center, SAIC Frederick, NCI-Frederick Cancer Research and Development Center, Frederick, MD 21702, USA
Received 14 June 2000; in ®nal form 14 September 2000
Abstract We applied sequential atomic multipole ®tting and the Gauss electric ®eld ¯ux theorem to calculate total molecular charge, charge distribution and charge transfer in a number of water, ion, and ion±water clusters. We have found that: (a) charge redistribution upon formation of ionic and of hydrogen bonds are similar and they are drastically dierent in covalent bonding; (b) polarization accounts for charge re-distribution upon formation of ionic or hydrogen bonds without invoking a charge transfer; (c) the `charge leak' from molecular cavities, encountered in the hydration theories, is usually compensated by the electron density of the surrounding water; (d) the ¯ux theorem provides alternative means of de®ning atomic radii in molecules in solution or clusters. Ó 2001 Elsevier Science B.V. All rights reserved.
1. Introduction A central role of electrostatic interactions in chemistry follows from the famous `electrostatic theorem' of Feynman [1]. It states that `the force on nucleus a due to the other nuclei and the electrons, in some particular Born±Oppenheimer nuclear con®guration, is just what would be computed from classical electrostatics from the location of the other nuclei and the electron charge density' [2]. The corresponding statement for the interaction energy is given by the Hellman±Feynman theorem [2,3]. Thus, descriptions of molecules, bonding, intermolecular interactions, and interactions of molecules with their environment on the basis of the electron density distributions attracted signi®cant attention (e.g., [3±8]). Many
*
Corresponding author. Fax: +1-201-692-0215. E-mail address: [email protected] (A.A. Rashin).
authors developed classical representations of molecular charge distributions to allow classical calculations of forces and energies in polyatomic molecules and their complexes (e.g., [9±14]). Recently it has been claimed on the basis of experimental electron density studies that hydrogen bonds in water have signi®cantly covalent character [15] (for reviews, see [16,17]) due to the charge transfer. This would make attempts at classical calculations more dubious, and add to the long dispute on the nature of the hydrogen bond (for a brief discussion see [18]). Relatedly the ionic character of a bond can often be presented through the magnitudes of charges ascribed to the atoms forming such a bond. Fully ionic bond would be represented by charges 1. However, because atomic electron clouds in molecules interpenetrate, often only a partial charge is ascribed to each atom participating in the ionic bond [3,19] suggesting a partially covalent character of ionic bonds. The diuse character of the atomic electron
0009-2614/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 0 ) 0 1 3 3 0 - 0
328
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
clouds creates problems in the dielectric theories of ionic or polar hydration [18,20]. In such theories hydration eects are due to the polarization of the solvent represented as a dielectric outside the solute cavity. The solvent polarization is caused by the charge distribution of the solute contained in this cavity. Diuse electron densities of many solutes extend far beyond the boundary of any reasonably sized cavity, leading to internal contradictions in the theory [18,20,21]. Here we use the multipole ®tting of molecular electrostatic potentials [13] and the charge integration based on the Gauss electric ®eld ¯ux theorem [20] to study the quantum mechanical charge distributions in ionic and uncharged molecules, in their complexes, and in clusters. Our study leads to new answers to the questions of `covalency' of the hydrogen bond, of `charge transfer' upon formation of ionic and hydrogen bonds, of the `charge leak' from solute cavities in the theories of hydration, and possibly to a new de®nition of atomic radii in crystals, clusters and condensed media.
2. Computational methods 2.1. Calculation of the total charge in a spherical volume According to the Gauss theorem the ¯ux of the electric ®eld through a closed surface equals 4p times the total charge inside the surface [20,22]. Thus one can place a closed surface around any individual molecule or such a molecule in the cluster, calculate quantum mechanically the ¯ux of the electric ®eld of the entire cluster (or of an individual molecule) through this surface, and obtain the total charge inside the surface from the value of the ¯ux [20]. One can similarly construct a surface of a small spherical probe anywhere in the electron density, and obtain the total charge inside it. We constructed the surfaces of molecular cavities and of the probe surfaces from analytical representations of curved fragments of the molecular surface [23]. In most calculations of the ¯ux through the surface of the molecular cavity we used the probe with radius
0.8 A for the construction of molecular surfaces [20]. However, in some cases when we wanted to exclude from our results electron densities from volumes formed by toroidal surfaces (a usual part of the molecular surfaces), we used very small probe radii (0.1 A). 2.2. Multipole ®tting The electrostatic potential of a complex polyatomic molecule and its changes can be accurately represented by the potential of its atomic multipoles representing its charge density. This allows for the classical implementations of the Hellman± Feynman theorem. The multipole ®tting to the electrostatic potential around individual molecules and clusters has been done with a modi®ed algorithm of Kong and Yan [13]. In this algorithm the electrostatic potential, EP, on a large set of points around a molecule or a cluster is obtained from a quantum mechanical calculation. In the ®rst ®tting step point charges on the nuclei are RMS ®tted to reproduce EP. Then the ®rst dierential EP is calculated as a dierence between EP from the quantum calculation and the one produced on the same set of points by the ®tted charges. In the next ®tting step point dipoles on the nuclei are RMS ®tted to the ®rst dierential EP. Then the second dierential EP is calculated and point quadrupoles are RMS ®tted to reproduce it in the third ®tting step. At each step DRMS is calculated to evaluate the agreement between EP's from the quantum mechanical electron density and from the ®tted multipoles. In our modi®ed algorithm any of the multipoles can be ®xed, and only other multipoles ®tted. One can calculate point charges for individual molecules and ®x their values in a cluster of such molecules. Then one can ®t in the complex only dipoles and quadrupoles. Dierences in their values in individual molecules and in their cluster provide the measure of polarization of individual molecules upon the complex formation. DRMS's obtained when all multipoles are ®tted can be compared to DRMS's obtained by ®tting of only higher multipoles when the point charges are ®xed to the values obtained for individual molecules. If
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
329
these compared DRMS are the same or very similar, it indicates that the electron density redistribution can be interpreted as a polarization with no necessity to invoke a charge transfer in the complex formation. 2.3. Quantum mechanical calculations Quantum mechanical calculations were performed with Gaussian94 DFT programs [24]. Cluster structures (described in detail in [25,26]) were optimized with 6-311+G(d,p) and 6-311++G(d,p) basis sets in B3LYP approximation. The structure of methane and argon dimers were optimized in the MP2 approximation to accurately account for the electron correlation which dominate the dispersion interactions. DFT calculations of electron density distributions were carried out with DZVPD basis set [27]. We used B3P86 exchange-correlation functional in Gaussian94 program [24] for our analysis of charge distributions, because it has been shown to produce signi®cantly more accurate dipole moments [26,28] than MP2 or B3LYP calculations. 3. Results and discussion Charges found with our calculations inside atomic or molecular cavities of dierent sizes are shown in Fig. 1. The data in this ®gure stresses an often forgotten fact that for any atom or molecule the electron density extends beyond the van der Waals or ionic radius. As a reasonable measure of how far the electron density extends, we chose the radius of the cavity with only 1% of the electron density being still outside it. Li has the most compact electron density of all atoms or molecules considered by us. For Li it extends beyond the for Ar2 it extends by 0.5 standard radius by 0.4 A, for Na by 0.7 A, for O by 0.8 A, for F by 0.9 A, for H, CH3 ; NH4 and NO it extends by 1 A, A, 2 To describe the and for Cl it extends by 1.2 A. charge attributed to each atom or molecule, forming van der Waals, ionic, or hydrogen bonds, the electron density in some contacting volumes of space is usually assigned to the bonded atoms or
Fig. 1. Diuse character of the electron density of ions, neutral atoms and molecules. For NO2 the same radii are used for N and O; for NH 4 the radius of H is ®xed at 1.16 A [20]; for CH3 but it does not make dierence at the radius of H is ®xed at 1A, larger C radii. Charges inside cavities with radii closest to the ionic or van der Waals radii are marked with bold short arrows.
molecules. It is quite obvious that the volume assigned to one member of the complex would contain a signi®cant part of the charge density from another member(s) of the complex. Thus charges in volumes assigned to dierent members of a complex might not indicate any real charge transfer but rather represent an artifact of the interpenetration of electron clouds. Instead of dividing the total electron density of a complex into parts contained in some volumes, we can analyze how the charge distribution expresses itself in the electrostatic potential (EP) of the complex outside its electron density distribution. EP can be approximated by the potential from a set of multipoles on atoms of the complex. If there has been a charge transfer upon the complex formation, it should lead to the change in the atomic charges (monopoles) of the complex compared to the atomic charges of the isolated constituent atoms or molecules prior to the complex formation. Results of our calculations of changes in the multipoles and EP upon formation of a number of complexes are shown in Table 1. We found that in all cases the change in EP upon the complex formation from its components can be very accurately represented by changes in only dipoles or quadrupoles with the monopoles kept unchanged
330
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
Table 1 Electrostatic potentials of complexes from DFT calculation and from ®tted multipoles Complexa
Mono-®t RMSb
Mono+dipole ®t RMSb
All ®t RMSb
Monomeric charges RMSc
Dipole Fit RMSc
Dipole+quadrupole ®t RMSc
Li F Li F 2 H4 O2 H6 O3
0.000901 0.000150 0.000344 0.006993
0.000043 0.000022 0.000039 0.005777
0.000030 0.000009 0.000008 0.004843
0.004743 0.004352 0.001997 0.007725
0.000079 0.000026 0.000043 0.005835
0.000014 0.000007 0.000009 0.004903
Four ionic/molecular complexes are represented here; Li F dimer is linear; water trimer forms a symmetric triangular structure. The sequential ®t ®rst to monopoles, then to dipoles, and then to quadrupoles was performed as described in methods. c Atomic monopoles are ®xed at e for LiF complexes, and at the monomeric atomic charges for water complexes (qO 0:680e; qH 0:340e), and then dipoles and quadrupoles are ®t as described in Section 2. a
b
from their values in the isolated components. This result strongly suggests that changes in the electron density upon the formation of most (if not all) ionic or hydrogen bonds are due to the polarization (expressed in the changed dipoles and quadrupoles) but not to the charge transfer. It has been noted [3,7] that the electron densities of the ionic or hydrogen bonded complexes is
usually within a few percent of the sum of the isolated atoms or molecules forming the complex. Usually the electron densities are represented either as two-dimensional density maps [3,29] which are often dicult to interpret in terms of actual numbers or percentages, or given in one particular point along the bond [3,7]. We have chosen to present the numeric values of the dierential elec-
Fig. 2. Comparison of the dierential electron density for the ionic, covalent and hydrogen bonds. The electron densities were calculated at points equally distant along the bond. For Li F , the nuclei are at points 7 and 21. For O2 , the nuclei are at points 6 and 16. For the water dimer, H2 O2 , the nucleus of H-atom in one water molecule is at point 1, and the nucleus of O-atom of another water molecule is at point 16. % change was calculated as the ratio of the the electron density of the complex minus the sum of densities of isolated components divided by this sum.
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
tric densities along the entire bonds in three representative ionic, covalent, and hydrogen bonded dimers (Fig. 2). We ®nd that indeed the electron density in the most typical hydrogen bond of the water dimer behaves very similarly to the electron density in the most ionic bond [3,19] of Li F . In both complexes the electron density along the bond equals the sum of the electron densities of individual water molecules or of individual ions within a few percent. Both complexes also exhibit a wave-like variation of the magnitude and sign of the deviation of the electron density in the complex from the sum of the
331
electron densities of individual water molecules or ions. These main properties of the electron densities in the hydrogen and ionic bonds are in a sharp contrast to the behavior of the electron density upon formation of a typical covalent bond of O2 (Fig. 2). The electron density along the covalent bond is more then twice higher than the sum of the electron densities of the isolated atoms. Upon formation of covalent bonds, the electron density of individual atoms or molecules is pulled from their periphery onto the bond [3]. Thus two types of direct analysis of electron densities support the notion that the hydrogen
Table 2 Total charge inside individual water cavities with varying atomic cavity radii in dierent water clusters Clustera
Water#b
Oxygen radiusc
Hydrogen radiusc
Total charge inside cavityd
Oxygen radiusc
Hydrogen radiusc
Total charge inside cavityd
H2 O1
1
1.4 1.5 1.5 1.5 1.7 1.9
1.0 1.0 1.1 1.16 1.1 1.1
0.333 0.264 0.235 0.219 0.143 0.080
1.9 1.9 2.2 2.2 2.3 2.3
1.5 1.9 1.7 2.2 1.7 2.2
0.051 0.027 0.023 0.006 0.012 0.004
H4 O2
1
1.5 1.5 1.5 1.5 1.6 1.7
1.1 1.16 1.26 1.16 1.16 1.1
0.117 0.068 )0.040 0.193 0.121 0.063
1.6 1.6 1.6 1.7 1.8
1.1 1.16 1.26 1.16 1.1
0.067 0.021 )0.08 0.054 )0.047
1.4 1.5 1.4 1.5 1.4 1.5
1.1 1.1 1.1 1.1 1.1 1.1
0.156 0.063 0.143 0.063 0.142 0.059
1.5
1.16
0.011
1.5
1.16
0.007
1.5
1.16
0.002
2
H6 O3
1 2 3
H8 O4 ttt
1
1.4 1.4
1.0 1.1
0.157 0.043
1.4 1.45
1.16 1.1
)0.044 0.001
H10 O5
1
1.4 1.4
1.0 1.1
0.126 0.003
1.45
1.1
)0.046
H12 O6 B
1 2 5
1.45 1.4 1.4 1.4 1.5 1.5
1.1 1.0 1.0 1.1 1.1 1.14
1.45 1.35 1.325 1.45
1.1 1.1 1.1 1.16
0.006 )0.029 )0.005 0.020
6 a
0.009 0.164 0.075 )0.076 0.044 0.004
Water clusters and their nomenclature are taken from [25]. In most cases we studies all or a few molecules in each cluster. c The radii of O and H were varied to ®nd the values for which the total charge in the molecular cavity is zero within 1%. d Total charge in the cavity is calculated with the Gauss ¯ux theorem as described in Section 2. b
332
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
bond is similar to the strongly ionic bond, and it does not involve any signi®cant charge transfer. Both bonds are apparently electronically very dierent from the covalent bond. This suggests that recent results of measurements of the properties of the electron densities in hydrogen bonded systems [15±17] might be better interpreted not in the terms of the `covalency' of hydrogen bonds, but in terms of the overlap of electron densities [17] which to a dierent degree is characteristic for all kinds of bonds (not only covalent ones). Our results also support the idea that accurate representation of electrostatic interactions can be achieved using empirical potentials with charges obtained with proper DFT EP ®tting augmented by an adequate empirical polarization scheme. As mentioned in Section 1, the diuse character of the atomic electron clouds creates problems in the dielectric theories of ionic or polar hydration (see [18,20]). Diuse electron densities of many solutes extend (`leak') far beyond the cavity of any reasonable size which boundary they are supposed to polarize. To mend this contradiction in the theory [18,20] it was suggested to represent the volume polarization of the solvent, caused by the solute's charge outside the cavity, by additional surface polarization charges [8,21]. However, eects of the signi®cant penetration of the
solvent electron density into the cavity occupied by the solute remain unaccounted for. To analyze the `leak' problem, we used our method of Gauss ¯ux integration to calculate total charge inside water cavities of dierent size in a variety of water clusters [25,26]. The results of our calculations are summarized in Table 2. These results show that in all cyclic or three-dimensional water clusters water cavities with radii within 10% for of those used in continuum theories (1.5 A oxygen, and 1.16 A for hydrogen, [20]) contain zero charge (within 1%) for a single water molecule. For an isolated water molecule the total charge in such a cavity is about 0:8. The charge outside the cavity is also large for the water dimer (Table 2). In larger water clusters the loss of the electron density from an individual water cavity is compensated by the penetrating diuse electron density of surrounding waters. Our results indicate that the published concerns and solutions to the problem of the `charge leak' might have been misplaced or incomplete, and that dierent approaches might be required. The problem is not that there is not enough charge inside the solute cavity, but how the solute's and solvent's charges polarize each other. We will address this question elsewhere. We applied the same Gauss ¯ux integration to ionic cavities in ion±water clusters studied earlier [26] (Table 3). We varied the radii of ionic cavities
Table 3 Total charge inside ionic cavities with varying radii in ion±water clusters.
a
Clustera
Ionb
Cavity radiusc
Charge inside cavityd
Clustera
Ionb
Cavity radiusc
Charge inside cavityd
F 6WD
F
1.423* 1.500 1.505 1.510 1.600 (1.330)
)0.8061 )0.9843 )0.9959 )1.0139 )1.2299
Li 6WA
Li
1.316* 0.800 0.750 0.700 (0.680)
0.0388 0.9737 1.0076 1.0433
Cl 6WD
Cl
1.700 1.937* 2.030 2.100 (1.810)
)0.3630 )0.8030 )1.0045 )1.1941
Na 6WA
Na
1.200 1.130 1.100 1.000 (0.980)
0.9350 1.0024 1.0316 1.1400
Ion±water clusters are taken from [26]. We considered two anions and two cations in this study. c The ionic cavity radii were varied to ®nd the values for which the total charge in the ionic cavity is e within 1%; cavity radii from [30] are marked by *; standard ionic radii are shown in parenthesis. d Total charge in the cavity is calculated with the Gauss ¯ux theorem as described in Section 2. b
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
in the clusters until they contain the total charge of 1e or 1e (within 1%). We ®nd that cavities containing charges of 1e or 1e have radii only 10±15% o the classical ionic radii. This magnitude of the deviation can be expected because of the small size of the clusters. However, the procedure can provide a new way to determine ionic radii in crystals and clusters. 4. Conclusions Our investigations have shown that: · The changes in the charge distribution upon formation of simple ionic bonds and of water±water hydrogen bonds are similar and drastically dierent from the corresponding changes in covalent bonding. · Partial charges on atoms are preserved upon formation of ion±ion or water clusters, and the concomitant charge redistribution can be fully accounted for by the polarization represented by induced atomic dipoles and sometimes quadrupoles, thus eliminating the need to invoke a `charge transfer'. · In water clusters the loss of the electron density from an individual water cavity is compensated by the penetrating diuse electron density of surrounding waters. · In ion±water clusters spherical volumes containing full ionic charge of e have radii similar to the classical ionic radii. Acknowledgements We thank the sta and administration of the Advanced Biomedical Computing Center for their support of this project. This project has been funded in whole or in part with federal funds from the National Cancer Institute, National Institutes of Health, under contract No. N01-C0-56000 and from grant No. 1R43 GM50663-01 A1. The content of this publication does not necessarily re¯ect the views or policies of the Department of Health and Human services, nor does mention of trade
333
names, commercial products, or organizations imply endorsement by the US Government.
References [1] R.P. Feynman, Phys. Rev. 56 (1939) 340. [2] R.G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989. [3] R.F.W. Bader, Atoms in Molecules. A Quantum Theory, Clarendon, Oxford, 1994. [4] A.D. Buckingham, P.W. Fowler, J.M. Hudson, Chem. Rev. 88 (1988) 963. [5] A.D. Buckingham, P.W. Fowler, J. Chem. Phys. 79 (1983) 6426. [6] A.D. Buckingham, P.W. Fowler, Can. J. Chem. 63 (1985) 2018. [7] J. Bentley, J. Phys. Chem. A 102 (1998) 6043. [8] D.M. Chipman, J. Chem. Phys. 106 (1997) 10194. [9] U.C. Singh, P.A. Kollman, J. Comput. Chem. 5 (1984) 129. [10] W.A. Sokalski, A. Sawaryn, J. Chem. Phys. 87 (1987) 526. [11] D.E. Williams, J. Comput. Chem. 9 (1988) 745. [12] K.M. Merz, J. Comput. Chem. 13 (1992) 749. [13] J. Kong, J.M. Yan, Int. J. Quant. Chem. 46 (1993) 239. [14] J. Li, T. Zhu, C.J. Cramer, D.G. Truhlar, J. Phys. Chem. A 102 (1998) 1820. [15] E.D. Isaacs, A. Shukla, P.M. Platzman, D.R. Hamann, B. Barbiellini, C.A. Talk, Phys. Rev. Lett. 82 (1999) 600. [16] A. Hellemans, Science 283 (1999) 614. [17] S. Borman, Chem. Eng. News 5 (10) (1999) 36. [18] A.A. Rashin, Prog. Biophys. Mol. Biol. 60 (1993) 73. [19] B.K. Vainshtein, V.M. Fridkin, V.L. Indenbom, Modern Crystallography, vol. 2, Nauka, Moscow, 1979. [20] A.A. Rashin, K. Namboodiri, J. Phys. Chem. 91 (1987) 6003. [21] B. Mennucci, J. Tomasi, J. Chem. Phys. 106 (1997) 5151. [22] J.D. Jackson, Classical Electrostatics, Wiley, New York, 1975. [23] A.A. Rashin, M.A. Bukatin, Biophys. Chem. 51 (1994) 167. [24] A.M. Frisch et al., GA U S S I A N 94, Gaussian Inc., Pittsburg, 1995. [25] G.J. Tawa, I.A. Topol, S.K. Burt, R.A. Caldwell, A.A. Rashin, J. Chem. Phys. 109 (1998) 4852. [26] I.A. Topol, G.J. Tawa, S.K. Burt, A.A. Rashin, J. Chem. Phys. 111 (1999) 10998. [27] A.A. Rashin, M.A. Bukatin, J. Andzelm, A.T. Hagler, Biophys. Chem. 51 (1994) 375. [28] A.A. Rashin, L. Young, I.A. Topol, S.K. Burt, Chem. Phys. Lett. (1994) 182. [29] J.J. Alkorta, H.O. Perez, Villar, J. Mol. Graphics 12 (1994) 3. [30] A.A. Rashin, B. Honig, J. Phys. Chem. 89 (1985) 5588.
Chemical Physics Letters 335 (2001) 327±333
www.elsevier.nl/locate/cplett
Charge distributions in water and ion±water clusters A.A. Rashin a,*, I.A. Topol b, G.J. Tawa b, S.K. Burt b b
a BioChemComp Inc., 543 Sagamore Ave, Teaneck, NJ 07666, USA Advanced Biomedical Computing Center, SAIC Frederick, NCI-Frederick Cancer Research and Development Center, Frederick, MD 21702, USA
Received 14 June 2000; in ®nal form 14 September 2000
Abstract We applied sequential atomic multipole ®tting and the Gauss electric ®eld ¯ux theorem to calculate total molecular charge, charge distribution and charge transfer in a number of water, ion, and ion±water clusters. We have found that: (a) charge redistribution upon formation of ionic and of hydrogen bonds are similar and they are drastically dierent in covalent bonding; (b) polarization accounts for charge re-distribution upon formation of ionic or hydrogen bonds without invoking a charge transfer; (c) the `charge leak' from molecular cavities, encountered in the hydration theories, is usually compensated by the electron density of the surrounding water; (d) the ¯ux theorem provides alternative means of de®ning atomic radii in molecules in solution or clusters. Ó 2001 Elsevier Science B.V. All rights reserved.
1. Introduction A central role of electrostatic interactions in chemistry follows from the famous `electrostatic theorem' of Feynman [1]. It states that `the force on nucleus a due to the other nuclei and the electrons, in some particular Born±Oppenheimer nuclear con®guration, is just what would be computed from classical electrostatics from the location of the other nuclei and the electron charge density' [2]. The corresponding statement for the interaction energy is given by the Hellman±Feynman theorem [2,3]. Thus, descriptions of molecules, bonding, intermolecular interactions, and interactions of molecules with their environment on the basis of the electron density distributions attracted signi®cant attention (e.g., [3±8]). Many
*
Corresponding author. Fax: +1-201-692-0215. E-mail address: [email protected] (A.A. Rashin).
authors developed classical representations of molecular charge distributions to allow classical calculations of forces and energies in polyatomic molecules and their complexes (e.g., [9±14]). Recently it has been claimed on the basis of experimental electron density studies that hydrogen bonds in water have signi®cantly covalent character [15] (for reviews, see [16,17]) due to the charge transfer. This would make attempts at classical calculations more dubious, and add to the long dispute on the nature of the hydrogen bond (for a brief discussion see [18]). Relatedly the ionic character of a bond can often be presented through the magnitudes of charges ascribed to the atoms forming such a bond. Fully ionic bond would be represented by charges 1. However, because atomic electron clouds in molecules interpenetrate, often only a partial charge is ascribed to each atom participating in the ionic bond [3,19] suggesting a partially covalent character of ionic bonds. The diuse character of the atomic electron
0009-2614/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 0 ) 0 1 3 3 0 - 0
328
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
clouds creates problems in the dielectric theories of ionic or polar hydration [18,20]. In such theories hydration eects are due to the polarization of the solvent represented as a dielectric outside the solute cavity. The solvent polarization is caused by the charge distribution of the solute contained in this cavity. Diuse electron densities of many solutes extend far beyond the boundary of any reasonably sized cavity, leading to internal contradictions in the theory [18,20,21]. Here we use the multipole ®tting of molecular electrostatic potentials [13] and the charge integration based on the Gauss electric ®eld ¯ux theorem [20] to study the quantum mechanical charge distributions in ionic and uncharged molecules, in their complexes, and in clusters. Our study leads to new answers to the questions of `covalency' of the hydrogen bond, of `charge transfer' upon formation of ionic and hydrogen bonds, of the `charge leak' from solute cavities in the theories of hydration, and possibly to a new de®nition of atomic radii in crystals, clusters and condensed media.
2. Computational methods 2.1. Calculation of the total charge in a spherical volume According to the Gauss theorem the ¯ux of the electric ®eld through a closed surface equals 4p times the total charge inside the surface [20,22]. Thus one can place a closed surface around any individual molecule or such a molecule in the cluster, calculate quantum mechanically the ¯ux of the electric ®eld of the entire cluster (or of an individual molecule) through this surface, and obtain the total charge inside the surface from the value of the ¯ux [20]. One can similarly construct a surface of a small spherical probe anywhere in the electron density, and obtain the total charge inside it. We constructed the surfaces of molecular cavities and of the probe surfaces from analytical representations of curved fragments of the molecular surface [23]. In most calculations of the ¯ux through the surface of the molecular cavity we used the probe with radius
0.8 A for the construction of molecular surfaces [20]. However, in some cases when we wanted to exclude from our results electron densities from volumes formed by toroidal surfaces (a usual part of the molecular surfaces), we used very small probe radii (0.1 A). 2.2. Multipole ®tting The electrostatic potential of a complex polyatomic molecule and its changes can be accurately represented by the potential of its atomic multipoles representing its charge density. This allows for the classical implementations of the Hellman± Feynman theorem. The multipole ®tting to the electrostatic potential around individual molecules and clusters has been done with a modi®ed algorithm of Kong and Yan [13]. In this algorithm the electrostatic potential, EP, on a large set of points around a molecule or a cluster is obtained from a quantum mechanical calculation. In the ®rst ®tting step point charges on the nuclei are RMS ®tted to reproduce EP. Then the ®rst dierential EP is calculated as a dierence between EP from the quantum calculation and the one produced on the same set of points by the ®tted charges. In the next ®tting step point dipoles on the nuclei are RMS ®tted to the ®rst dierential EP. Then the second dierential EP is calculated and point quadrupoles are RMS ®tted to reproduce it in the third ®tting step. At each step DRMS is calculated to evaluate the agreement between EP's from the quantum mechanical electron density and from the ®tted multipoles. In our modi®ed algorithm any of the multipoles can be ®xed, and only other multipoles ®tted. One can calculate point charges for individual molecules and ®x their values in a cluster of such molecules. Then one can ®t in the complex only dipoles and quadrupoles. Dierences in their values in individual molecules and in their cluster provide the measure of polarization of individual molecules upon the complex formation. DRMS's obtained when all multipoles are ®tted can be compared to DRMS's obtained by ®tting of only higher multipoles when the point charges are ®xed to the values obtained for individual molecules. If
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
329
these compared DRMS are the same or very similar, it indicates that the electron density redistribution can be interpreted as a polarization with no necessity to invoke a charge transfer in the complex formation. 2.3. Quantum mechanical calculations Quantum mechanical calculations were performed with Gaussian94 DFT programs [24]. Cluster structures (described in detail in [25,26]) were optimized with 6-311+G(d,p) and 6-311++G(d,p) basis sets in B3LYP approximation. The structure of methane and argon dimers were optimized in the MP2 approximation to accurately account for the electron correlation which dominate the dispersion interactions. DFT calculations of electron density distributions were carried out with DZVPD basis set [27]. We used B3P86 exchange-correlation functional in Gaussian94 program [24] for our analysis of charge distributions, because it has been shown to produce signi®cantly more accurate dipole moments [26,28] than MP2 or B3LYP calculations. 3. Results and discussion Charges found with our calculations inside atomic or molecular cavities of dierent sizes are shown in Fig. 1. The data in this ®gure stresses an often forgotten fact that for any atom or molecule the electron density extends beyond the van der Waals or ionic radius. As a reasonable measure of how far the electron density extends, we chose the radius of the cavity with only 1% of the electron density being still outside it. Li has the most compact electron density of all atoms or molecules considered by us. For Li it extends beyond the for Ar2 it extends by 0.5 standard radius by 0.4 A, for Na by 0.7 A, for O by 0.8 A, for F by 0.9 A, for H, CH3 ; NH4 and NO it extends by 1 A, A, 2 To describe the and for Cl it extends by 1.2 A. charge attributed to each atom or molecule, forming van der Waals, ionic, or hydrogen bonds, the electron density in some contacting volumes of space is usually assigned to the bonded atoms or
Fig. 1. Diuse character of the electron density of ions, neutral atoms and molecules. For NO2 the same radii are used for N and O; for NH 4 the radius of H is ®xed at 1.16 A [20]; for CH3 but it does not make dierence at the radius of H is ®xed at 1A, larger C radii. Charges inside cavities with radii closest to the ionic or van der Waals radii are marked with bold short arrows.
molecules. It is quite obvious that the volume assigned to one member of the complex would contain a signi®cant part of the charge density from another member(s) of the complex. Thus charges in volumes assigned to dierent members of a complex might not indicate any real charge transfer but rather represent an artifact of the interpenetration of electron clouds. Instead of dividing the total electron density of a complex into parts contained in some volumes, we can analyze how the charge distribution expresses itself in the electrostatic potential (EP) of the complex outside its electron density distribution. EP can be approximated by the potential from a set of multipoles on atoms of the complex. If there has been a charge transfer upon the complex formation, it should lead to the change in the atomic charges (monopoles) of the complex compared to the atomic charges of the isolated constituent atoms or molecules prior to the complex formation. Results of our calculations of changes in the multipoles and EP upon formation of a number of complexes are shown in Table 1. We found that in all cases the change in EP upon the complex formation from its components can be very accurately represented by changes in only dipoles or quadrupoles with the monopoles kept unchanged
330
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
Table 1 Electrostatic potentials of complexes from DFT calculation and from ®tted multipoles Complexa
Mono-®t RMSb
Mono+dipole ®t RMSb
All ®t RMSb
Monomeric charges RMSc
Dipole Fit RMSc
Dipole+quadrupole ®t RMSc
Li F Li F 2 H4 O2 H6 O3
0.000901 0.000150 0.000344 0.006993
0.000043 0.000022 0.000039 0.005777
0.000030 0.000009 0.000008 0.004843
0.004743 0.004352 0.001997 0.007725
0.000079 0.000026 0.000043 0.005835
0.000014 0.000007 0.000009 0.004903
Four ionic/molecular complexes are represented here; Li F dimer is linear; water trimer forms a symmetric triangular structure. The sequential ®t ®rst to monopoles, then to dipoles, and then to quadrupoles was performed as described in methods. c Atomic monopoles are ®xed at e for LiF complexes, and at the monomeric atomic charges for water complexes (qO 0:680e; qH 0:340e), and then dipoles and quadrupoles are ®t as described in Section 2. a
b
from their values in the isolated components. This result strongly suggests that changes in the electron density upon the formation of most (if not all) ionic or hydrogen bonds are due to the polarization (expressed in the changed dipoles and quadrupoles) but not to the charge transfer. It has been noted [3,7] that the electron densities of the ionic or hydrogen bonded complexes is
usually within a few percent of the sum of the isolated atoms or molecules forming the complex. Usually the electron densities are represented either as two-dimensional density maps [3,29] which are often dicult to interpret in terms of actual numbers or percentages, or given in one particular point along the bond [3,7]. We have chosen to present the numeric values of the dierential elec-
Fig. 2. Comparison of the dierential electron density for the ionic, covalent and hydrogen bonds. The electron densities were calculated at points equally distant along the bond. For Li F , the nuclei are at points 7 and 21. For O2 , the nuclei are at points 6 and 16. For the water dimer, H2 O2 , the nucleus of H-atom in one water molecule is at point 1, and the nucleus of O-atom of another water molecule is at point 16. % change was calculated as the ratio of the the electron density of the complex minus the sum of densities of isolated components divided by this sum.
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
tric densities along the entire bonds in three representative ionic, covalent, and hydrogen bonded dimers (Fig. 2). We ®nd that indeed the electron density in the most typical hydrogen bond of the water dimer behaves very similarly to the electron density in the most ionic bond [3,19] of Li F . In both complexes the electron density along the bond equals the sum of the electron densities of individual water molecules or of individual ions within a few percent. Both complexes also exhibit a wave-like variation of the magnitude and sign of the deviation of the electron density in the complex from the sum of the
331
electron densities of individual water molecules or ions. These main properties of the electron densities in the hydrogen and ionic bonds are in a sharp contrast to the behavior of the electron density upon formation of a typical covalent bond of O2 (Fig. 2). The electron density along the covalent bond is more then twice higher than the sum of the electron densities of the isolated atoms. Upon formation of covalent bonds, the electron density of individual atoms or molecules is pulled from their periphery onto the bond [3]. Thus two types of direct analysis of electron densities support the notion that the hydrogen
Table 2 Total charge inside individual water cavities with varying atomic cavity radii in dierent water clusters Clustera
Water#b
Oxygen radiusc
Hydrogen radiusc
Total charge inside cavityd
Oxygen radiusc
Hydrogen radiusc
Total charge inside cavityd
H2 O1
1
1.4 1.5 1.5 1.5 1.7 1.9
1.0 1.0 1.1 1.16 1.1 1.1
0.333 0.264 0.235 0.219 0.143 0.080
1.9 1.9 2.2 2.2 2.3 2.3
1.5 1.9 1.7 2.2 1.7 2.2
0.051 0.027 0.023 0.006 0.012 0.004
H4 O2
1
1.5 1.5 1.5 1.5 1.6 1.7
1.1 1.16 1.26 1.16 1.16 1.1
0.117 0.068 )0.040 0.193 0.121 0.063
1.6 1.6 1.6 1.7 1.8
1.1 1.16 1.26 1.16 1.1
0.067 0.021 )0.08 0.054 )0.047
1.4 1.5 1.4 1.5 1.4 1.5
1.1 1.1 1.1 1.1 1.1 1.1
0.156 0.063 0.143 0.063 0.142 0.059
1.5
1.16
0.011
1.5
1.16
0.007
1.5
1.16
0.002
2
H6 O3
1 2 3
H8 O4 ttt
1
1.4 1.4
1.0 1.1
0.157 0.043
1.4 1.45
1.16 1.1
)0.044 0.001
H10 O5
1
1.4 1.4
1.0 1.1
0.126 0.003
1.45
1.1
)0.046
H12 O6 B
1 2 5
1.45 1.4 1.4 1.4 1.5 1.5
1.1 1.0 1.0 1.1 1.1 1.14
1.45 1.35 1.325 1.45
1.1 1.1 1.1 1.16
0.006 )0.029 )0.005 0.020
6 a
0.009 0.164 0.075 )0.076 0.044 0.004
Water clusters and their nomenclature are taken from [25]. In most cases we studies all or a few molecules in each cluster. c The radii of O and H were varied to ®nd the values for which the total charge in the molecular cavity is zero within 1%. d Total charge in the cavity is calculated with the Gauss ¯ux theorem as described in Section 2. b
332
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
bond is similar to the strongly ionic bond, and it does not involve any signi®cant charge transfer. Both bonds are apparently electronically very dierent from the covalent bond. This suggests that recent results of measurements of the properties of the electron densities in hydrogen bonded systems [15±17] might be better interpreted not in the terms of the `covalency' of hydrogen bonds, but in terms of the overlap of electron densities [17] which to a dierent degree is characteristic for all kinds of bonds (not only covalent ones). Our results also support the idea that accurate representation of electrostatic interactions can be achieved using empirical potentials with charges obtained with proper DFT EP ®tting augmented by an adequate empirical polarization scheme. As mentioned in Section 1, the diuse character of the atomic electron clouds creates problems in the dielectric theories of ionic or polar hydration (see [18,20]). Diuse electron densities of many solutes extend (`leak') far beyond the cavity of any reasonable size which boundary they are supposed to polarize. To mend this contradiction in the theory [18,20] it was suggested to represent the volume polarization of the solvent, caused by the solute's charge outside the cavity, by additional surface polarization charges [8,21]. However, eects of the signi®cant penetration of the
solvent electron density into the cavity occupied by the solute remain unaccounted for. To analyze the `leak' problem, we used our method of Gauss ¯ux integration to calculate total charge inside water cavities of dierent size in a variety of water clusters [25,26]. The results of our calculations are summarized in Table 2. These results show that in all cyclic or three-dimensional water clusters water cavities with radii within 10% for of those used in continuum theories (1.5 A oxygen, and 1.16 A for hydrogen, [20]) contain zero charge (within 1%) for a single water molecule. For an isolated water molecule the total charge in such a cavity is about 0:8. The charge outside the cavity is also large for the water dimer (Table 2). In larger water clusters the loss of the electron density from an individual water cavity is compensated by the penetrating diuse electron density of surrounding waters. Our results indicate that the published concerns and solutions to the problem of the `charge leak' might have been misplaced or incomplete, and that dierent approaches might be required. The problem is not that there is not enough charge inside the solute cavity, but how the solute's and solvent's charges polarize each other. We will address this question elsewhere. We applied the same Gauss ¯ux integration to ionic cavities in ion±water clusters studied earlier [26] (Table 3). We varied the radii of ionic cavities
Table 3 Total charge inside ionic cavities with varying radii in ion±water clusters.
a
Clustera
Ionb
Cavity radiusc
Charge inside cavityd
Clustera
Ionb
Cavity radiusc
Charge inside cavityd
F 6WD
F
1.423* 1.500 1.505 1.510 1.600 (1.330)
)0.8061 )0.9843 )0.9959 )1.0139 )1.2299
Li 6WA
Li
1.316* 0.800 0.750 0.700 (0.680)
0.0388 0.9737 1.0076 1.0433
Cl 6WD
Cl
1.700 1.937* 2.030 2.100 (1.810)
)0.3630 )0.8030 )1.0045 )1.1941
Na 6WA
Na
1.200 1.130 1.100 1.000 (0.980)
0.9350 1.0024 1.0316 1.1400
Ion±water clusters are taken from [26]. We considered two anions and two cations in this study. c The ionic cavity radii were varied to ®nd the values for which the total charge in the ionic cavity is e within 1%; cavity radii from [30] are marked by *; standard ionic radii are shown in parenthesis. d Total charge in the cavity is calculated with the Gauss ¯ux theorem as described in Section 2. b
A.A. Rashin et al. / Chemical Physics Letters 335 (2001) 327±333
in the clusters until they contain the total charge of 1e or 1e (within 1%). We ®nd that cavities containing charges of 1e or 1e have radii only 10±15% o the classical ionic radii. This magnitude of the deviation can be expected because of the small size of the clusters. However, the procedure can provide a new way to determine ionic radii in crystals and clusters. 4. Conclusions Our investigations have shown that: · The changes in the charge distribution upon formation of simple ionic bonds and of water±water hydrogen bonds are similar and drastically dierent from the corresponding changes in covalent bonding. · Partial charges on atoms are preserved upon formation of ion±ion or water clusters, and the concomitant charge redistribution can be fully accounted for by the polarization represented by induced atomic dipoles and sometimes quadrupoles, thus eliminating the need to invoke a `charge transfer'. · In water clusters the loss of the electron density from an individual water cavity is compensated by the penetrating diuse electron density of surrounding waters. · In ion±water clusters spherical volumes containing full ionic charge of e have radii similar to the classical ionic radii. Acknowledgements We thank the sta and administration of the Advanced Biomedical Computing Center for their support of this project. This project has been funded in whole or in part with federal funds from the National Cancer Institute, National Institutes of Health, under contract No. N01-C0-56000 and from grant No. 1R43 GM50663-01 A1. The content of this publication does not necessarily re¯ect the views or policies of the Department of Health and Human services, nor does mention of trade
333
names, commercial products, or organizations imply endorsement by the US Government.
References [1] R.P. Feynman, Phys. Rev. 56 (1939) 340. [2] R.G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989. [3] R.F.W. Bader, Atoms in Molecules. A Quantum Theory, Clarendon, Oxford, 1994. [4] A.D. Buckingham, P.W. Fowler, J.M. Hudson, Chem. Rev. 88 (1988) 963. [5] A.D. Buckingham, P.W. Fowler, J. Chem. Phys. 79 (1983) 6426. [6] A.D. Buckingham, P.W. Fowler, Can. J. Chem. 63 (1985) 2018. [7] J. Bentley, J. Phys. Chem. A 102 (1998) 6043. [8] D.M. Chipman, J. Chem. Phys. 106 (1997) 10194. [9] U.C. Singh, P.A. Kollman, J. Comput. Chem. 5 (1984) 129. [10] W.A. Sokalski, A. Sawaryn, J. Chem. Phys. 87 (1987) 526. [11] D.E. Williams, J. Comput. Chem. 9 (1988) 745. [12] K.M. Merz, J. Comput. Chem. 13 (1992) 749. [13] J. Kong, J.M. Yan, Int. J. Quant. Chem. 46 (1993) 239. [14] J. Li, T. Zhu, C.J. Cramer, D.G. Truhlar, J. Phys. Chem. A 102 (1998) 1820. [15] E.D. Isaacs, A. Shukla, P.M. Platzman, D.R. Hamann, B. Barbiellini, C.A. Talk, Phys. Rev. Lett. 82 (1999) 600. [16] A. Hellemans, Science 283 (1999) 614. [17] S. Borman, Chem. Eng. News 5 (10) (1999) 36. [18] A.A. Rashin, Prog. Biophys. Mol. Biol. 60 (1993) 73. [19] B.K. Vainshtein, V.M. Fridkin, V.L. Indenbom, Modern Crystallography, vol. 2, Nauka, Moscow, 1979. [20] A.A. Rashin, K. Namboodiri, J. Phys. Chem. 91 (1987) 6003. [21] B. Mennucci, J. Tomasi, J. Chem. Phys. 106 (1997) 5151. [22] J.D. Jackson, Classical Electrostatics, Wiley, New York, 1975. [23] A.A. Rashin, M.A. Bukatin, Biophys. Chem. 51 (1994) 167. [24] A.M. Frisch et al., GA U S S I A N 94, Gaussian Inc., Pittsburg, 1995. [25] G.J. Tawa, I.A. Topol, S.K. Burt, R.A. Caldwell, A.A. Rashin, J. Chem. Phys. 109 (1998) 4852. [26] I.A. Topol, G.J. Tawa, S.K. Burt, A.A. Rashin, J. Chem. Phys. 111 (1999) 10998. [27] A.A. Rashin, M.A. Bukatin, J. Andzelm, A.T. Hagler, Biophys. Chem. 51 (1994) 375. [28] A.A. Rashin, L. Young, I.A. Topol, S.K. Burt, Chem. Phys. Lett. (1994) 182. [29] J.J. Alkorta, H.O. Perez, Villar, J. Mol. Graphics 12 (1994) 3. [30] A.A. Rashin, B. Honig, J. Phys. Chem. 89 (1985) 5588.