ARTICLE pubs.acs.org/JPCA
Ion Induced Dipole Clusters Hn (3 e n-odd e 13): Density Functional Theory Calculations of Structure and Energy Lulu Huang,*,† Cherif F. Matta,‡,§ and Lou Massa*,|| †
Center for Computational Materials Science, Naval Research Laboratory, Washington, D.C. 20375-5341, United States Department of Chemistry, Mount Saint Vincent University, Halifax, Nova Scotia, Canada B3M2J6 ‡ Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada, B3H 4J3 Hunter College and the Graduate School, City University of New York, New York, New York 10065, United States
)
§
ABSTRACT: We investigate anew the possible equilibrium geometries of ion induced dipole clusters of hydrogen molecular ions, of molecular formula Hn (3 e n-odd e 13). Our previous publications [Sapse, A. M.; et al. Nature 1979, 278, 332; Rayez, J. C.; et al., J. Chem. Phys. 1981, 75, 5393] indicated these molecules would have a shallow minimum and adopt symmetrical geometries that accord with the valence shell electron pair repulsion (VSEPR) rules for geometries defined by electron pairs surrounding a central point of attraction. These earlier calculations were all based upon Hartree Fock (HF) calculations with a fairly small basis of atomic functions, except for the H3 ion for which configuration interaction (CI) calculations were carried out. A related paper [Hirao, K.; et al., Chem. Phys. 1983, 80, 237] carried out similar calculations on the same clusters, finding geometries similar to our earlier calculations. However, although that paper argued that the stabilization energy of negative ion clusters Hn is small, vibration frequencies for the whole set of clusters was not reported, and so a definitive assertion of a true equilibrium was not present. In this paper we recalculate the energetics of the ion induced dipole clusters using density function theory (DFT) B3LYP method calculations in a basis of functions (6-311++G(d,p)). By calculating the vibration frequencies of the VSEPR geometries, we prove that in general they are not true minima because not all the resulting frequencies correspond to real values. By searching the energy surface of the B3LYP calculations, we find the true minimum geometries, which are surprising configurations and are perhaps counterintuitive. We calculate the total energy and binding energy of the new geometries. We also calculate the bond paths associated with the quantum theory of atoms in molecules (QTAIM). The B3LYP/6-311++G(d,p) results, for each molecule, deliver bond paths that radiate between each polarized H2 molecule and the polarizing H ion.
’ INTRODUCTION Ab-initio Hartree Fock calculations were carried out1 4 some time ago with a view toward suggesting that Hn ions (3 e n-odd e 13) held together by ion induced dipole bonds might possibly be extant in the interstellar medium. Both H2 and H are known species of the interstellar medium so a possible suggestion for the formation of their bound combinations was investigated. Such bound combinations might be formed upon the surface of carbon grains within interstellar dust clouds. The question of main interest was whether or not there existed equilibrium geometries for the presumed ion induced dipole clusters. Given the interstellar temperatures in dust clouds, less than 10 K, even a quite small well depth would prove sufficient to prevent vibrational breakup of the presumed ion induced dipole bonds. The initial calculations were indeed consistent with expectations for welldefined equilibrium structures for the suggested molecules. In fact, the optimized structures were intuitively appealing, indicating these molecules would have a shallow minimum and adopt symmetrical geometries that accord with the valence shell electron pair repulsion (VSEPR) rules for geometries defined by r 2011 American Chemical Society
electron pairs surrounding a central point of attraction. These early calculations were mainly based upon Hartree Fock (HF)5 calculations with a fairly small basis of atomic functions. A related paper6 carried out similar calculations on the same molecules finding geometries similar to our earlier calculations;7,8 however, although the paper argued that the stabilization energy of negative ion clusters Hn is small, the general set of vibration frequencies is not reported. Therefore whether or not the reported geometries are true equilibrium geometries may not be definitively inferred. In this paper we recalculate the energetics of the ion induced dipole molecules using density function theory (DFT)9 11 B3LYP12,13 method calculations in a basis of functions (6-311++G(d,p)).14 16 By calculating the vibration frequencies of the VSEPR geometries, we prove that in general they are not true Special Issue: Richard F. W. Bader Festschrift Received: April 27, 2011 Revised: June 28, 2011 Published: July 25, 2011 12445
dx.doi.org/10.1021/jp203913n | J. Phys. Chem. A 2011, 115, 12445–12450
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Optimized geometry for the molecules H n , 3 e n-odd e 13. The calculations are all DFT B3LYP/6-311++G(d,p). The molecules in the left column have the VSEPR geometries but correspond (except for H 3 and H5 ) to imaginary vibration frequencies and thus cannot be true energy minima. The geometries in the right column correspond to true energy minima having all positive vibration frequencies. In the simplest case, that of H 3 , there is only one converged geometry, which is linear and corresponds to vibration frequencies that are all positive. (a) H3 for VSEPR and new geometry, with no imaginary frequency. (b) H5 , VSEPR geometry (left) and new geometry (right), with no imaginary frequencies. (c) H7 , VSEPR geometry (left) with two imaginary frequencies and new geometry (right) with no imaginary frequency. (d) H 9 , VSEPR geometry (left) with five imaginary frequencies and new geometry (right) with no imaginary frequency. (e) H11 , VSEPR geometry (left) with five imaginary frequencies and new geometry (right) with no imaginary frequency. (f) H13 , VSEPR geometry (left) with three imaginary frequencies and new geometry (right) with no imaginary frequency.
minima because not all the resulting frequencies correspond to real values. However, searching the energy surface of the B3LYP calculations allows finding true minimum geometries, which are surprising and one may say almost counterintuitive. The total energy and binding energy of the new geometries are calculated. Also calculated are the bond paths associated with the Bader quantum topology of QTAIM17. The B3LYP/6-311++G(d,p) results, for each molecule, deliver bond paths that radiate between each polarized H2 molecule and the polarizing H ion. As with the original calculations that suggested the possible existence of these molecules in an environment such as the interstellar medium, having an extraordinarily low temperature and molecular density, here too the principle question of interest is whether or not there exist geometries based upon true energetic minima associated with the mechanism of the ion induced dipole bonds.
’ RESULTS We had previously calculated the structure of these molecules using Hartree Fock calculations and a limited basis of
atomic orbitals. To improve those calculations, we use here DFT-B3LYPcalculations that can account for correlation energy, and we employ a fairly large basis of Gaussian orbitals (6-311++G(d,p)). Gaussian09 is used for geometry optimization and vibrational frequency calculations.18 Two sets of optimized geometries occur, both of which are displayed in Figure 1. First, there are a set of VSEPR geometries listed down the left column in the figure. With the exception of H3 and H5 , these geometries are all associated with imaginary vibration frequencies and therefore cannot be true equilibrium geometries. The geometries for n = 7, 9, 11, and 13 have respectively the following number of imaginary vibration frequencies 2, 5, 5, and 3. Second, the geometries that correspond to true minima, thus having all real vibration frequencies, are listed in the right column of the figure. For both sets of geometries the calculated Mulliken charges are listed for every atom. In Figure 2 are shown the atomic charges for the ion induced dipole clusters of hydrogen calculated in a different 12446
dx.doi.org/10.1021/jp203913n |J. Phys. Chem. A 2011, 115, 12445–12450
The Journal of Physical Chemistry A
ARTICLE
Figure 2. QTAIM charges q(A) (in italic font) [au] and bond paths [au] superimposed upon equilibrium geometries of molecules Hn (3 e n-odd e 13). The small balls between two atoms are the bond critical points, and the numbers next to the points are the bond path length.
way, viz., according to the Bader methodology of QTAIM. Also shown in the same figure are the calculated bond critical points associated with the electron density of the DFT true equilibrium geometries. The charges and bond critical points are superimposed upon the equilibrium structures. The AIMALL program19 is used for QTAIM atomic properties, q(A) net charge of atom, and bond path calculations. The energetics of the VSEPR geometries and the true equilibrium geometries, have also been calculated. The energetics are displayed in Table 1 where, as in Figure 1, the quantities corresponding to VSEPR geometries are listed in the left column, and those corresponding to the true
equilibrium geometries are listed in the right column. The quantities shown for each of the clusters are the total molecular energy, the orbital energies of the HOMO and LUMO molecular orbitals, and the energy gap given by their difference. Also shown are the binding energy for each molecular structure calculated as the difference between the total energy of the full clusters, and the sum of its components considered as free H2 molecules and a free H ion. A classical electrostatic calculation was carried out seeking an intuitive reason to support the cause of stabilization energy in the new geometries of clusters Hn . Using H5 as an example, Figure 3 shows the results associated with the classical electrostatic 12447
dx.doi.org/10.1021/jp203913n |J. Phys. Chem. A 2011, 115, 12445–12450
The Journal of Physical Chemistry A
ARTICLE
(2) The angle θ between the two arms of the “V” is varied at 5° steps keeping bond lengths fixed. (3) The atomic charges on H+ and H are kept constant at their approximate QTAIM/DFT values calculated at the equilibrium geometry of q(H ) = 0.16 au = 2.5632 10 20 C and q(H+) = +0.08 au =1.2816 10 20 C. (4) As can be seen from Figure 3, there is no minimum on the “classical electrostatic potential energy surface” whereby the sum of the two attractive and the two repulsive electrostatic interactions show no hint of any potential energy well at or near the equilibrium angle of 67.6°. This zeroth-order approximation (i.e., one that considers atoms as point charges without consideration of the higher multipolar contributions) suggests that the cluster adopts the V-shape not due to these classical zero-order electrostatic interaction terms. Higher terms and, likely, quantum effects must be driving the cluster to adopt this particular geometry at equilibrium.
Table 1. Energy Quantities Corresponding to Optimized Geometries of the Ion Induced Dipole Hydrogen Clustersa VSEPR geomb H3
E(RB3LYP)
0.04249 0.13342
0.04249 0.13342
gap
0.09093
0.09093
H9
H11
1.690
8.6085
8.6085
point group E(RB3LYP)
C∞v
C∞v
HOMO
0.03966
0.03743
LUMO gap
0.12880 0.08914
0.12957 0.09214
2.89885068
3.479
2.89926879
3.741
dipole moment
0.0000
7.1075
point group
D∞h
Cs
E(RB3LYP)
4.08139642
4.08214086
HOMO
0.03657
0.03211
LUMO
0.12916
0.12741
gap binding energy
0.09259 5.346
0.09530 5.814
dipole moment
0.0000
4.2497
point group
D3h
Cs
E(RB3LYP)
5.26440858
5.26595642
HOMO
0.03125
0.02225
LUMO
0.12736
0.12488
gap
0.09611
0.10263
binding energy dipole moment
7.494 0.0000
8.475 3.4590
point group
Td
C1
E(RB3LYP)
6.44811307
6.44905801
HOMO
0.02111
0.01598
LUMO
0.12425
0.12225
gap
0.10314
0.10627
binding energy dipole moment point group H13
1.690
dipole moment
binding energy
H7
1.71642841
HOMO LUMO binding energy
H5
1.71642841
new geomb
E(RB3LYP)
10.093 0.0000 C3v 7.63203842
10.691 1.9336 Cs 7.63212087
HOMO
0.01065
0.00992
LUMO
0.12147
0.12136
gap
0.1108
0.1114
binding energy
12.826
12.880
dipole moment
0.0000
0.4463
point group
Oh
C1
a
As in Figure 1, the left column quantities correspond to VSEPR structures, and the right column quantities correspond to true minimal energy geometries. Listed for each molecule are the total energy, binding energy, the orbital energies of the HOMO and LUMO orbitals, and the energy gap corresponding to their difference. b All energies in atomic units (au), except binding energies which are listed in kcal/mol. Dipole moments are in debyes (D).
interactions between the two “Arms” of the H5 molecule. Approximations employed to implement the classical study are as follows: (1) Bond lengths are those in the DFT optimized V-geometry.
’ DISCUSSION AND CONCLUSIONS All the quantum calculations of this paper are based upon a chemical model that employs DFT Kohn Sham equations of motion implemented with 6-311++G(d,p) basis functions. Thus all the calculations are probably of comparable accuracy. It is of interest that the intuitively appealing long ago published VSEPR geometries, whose point groups are indicated in Table 1, do not correspond to true energy minimum geometries. This has been proven of course by calculating the vibration frequencies associated with that geometry. The mathematical characteristic of a minimum in the molecular energy surface is that all the associated 3N 6 vibration frequencies must be found real. This decidedly is not generally the case with the beautifully symmetric VSEPR geometries. Responding to this surprise, a search over the molecular energy surface did find in every case, a new geometry that did satisfy the criterion for equilibrium, viz., all positive vibration frequencies. Both the original supposed VSEPR geometries and the true equilibrium geometries may be found in Figure 1. The charge on each atom is given, for each of the structures in Table 2. Those charges are calculated according to the Mulliken formalism that has been so widely used. For the case of the chemical model employed here, the Mulliken charges do not always correspond to an alternation of positive and negative atomic charges, suggestive of electrostatic stability. To pick some particularly striking examples, in Table 2 notice that in the VSEPR geometry for H7 all the atoms are negative, or in the VSEPR geometry for H13 the central hydrogen ion has actually become (improbably) overall positive in charge. To get a more realistic assement of the atomic charges, they are calculated in accordance with the QTAIM formalism of Bader. In Figure 2 the QTAIM atomic charges are superimposed on the true equilibrium geometries for the molecules. They can be compared directly to the Mulliken charges. The QTAIM atomic charges in the case of each molecular structure alternate between negative and positive values in such fashion as to be consistent with electrostatic stability. Notice that the hydrogen ion remains negative in every structure, in contrast to the earlier mentioned Mulliken charge behavior. Figure 2 also indicates the occurrence of QTAIM bond paths and the position of all their bond critical points (BCP) in the molecular structure. These occur in a uniform way shared by all 12448
dx.doi.org/10.1021/jp203913n |J. Phys. Chem. A 2011, 115, 12445–12450
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Classical zeroth-order electrostatic interaction between the two “arms” of the V-shaped cluster using the QTAIM/DFT point charges obtained at the equilibrium geometry. Color scheme: blue = the two Coulombic attractive interactions between positive and negative hydrogens; green = the repulsive interaction between the two negatively charged hydrogens; red = the repulsive interaction between the two positively charged hydrogens. Black represents the sum of the four interactions, i.e., an approximation to the “classical electrostatic potential energy surface” including only atomic monopoles.
Table 2. Mulliken Atomic Charges for VSEPR and for the New Geometries of This Work atoms H3
H5
H7
H9
charges (VSEPR Geom)
charges (new geom)
atoms H11
charges (VSEPR geom)
charges (new geom)
H
0.9594
H
0.1252
Ha
0.0660
Ha1
0.1294
0.0812
Hb H
0.1067 0.7710
Hb1 Ha2
0.0637 0.1294
0.1458 0.0124
Ha1
0.0516
0.0503
Hb2
0.0637
0.1089
Hb1
0.0629
0.1097
Ha3
0.1294
0.0124
Ha2 Hb2 H
0.0516 0.0629 0.6448
0.0503 0.1097 0.7472
Hb3 Ha4 Hb4
0.0637 0.0403 0.1075
0.1089 0.0250 0.1214
Ha1
0.0513
0.0120
Ha5
0.0403
0.0250
Hb1
0.0671
0.0913
Hb5
0.1075
0.1214
Ha2
0.0513
0.0120
H
0.0314
0.3275
Hb2 Ha3
0.0671 0.0513
0.0913 0.0833
Ha1 Hb1
0.0847 0.0872
0.0069 0.1162
Hb3 H
0.0671 0.2363
0.1294 0.6576
Ha2 Hb2
0.0847 0.0872
0.0068 0.1162
Ha1
0.1443
0.0112
Ha3
0.0847
0.0328
Hb1
0.0466
0.1095
Hb3
0.0872
0.1325
Ha2
0.1443
0.0112
Ha4
0.0847
0.0329
Hb2 Ha3
0.0466 0.1443
0.1095 0.0632
Hb4 Ha5
0.0872 0.0847
0.1325 0.0354
Hb3
0.0466
0.1361
Hb5
0.0872
0.0918
Ha4
0.1443
0.0632
Ha6
0.0847
0.0354
Hb4
0.0466
0.1361
Hb6
0.0872
0.0919
0.8811
the molecular structures of Figure 2. There is a BCP between each pair of hydrogen atoms forming the hydrogen molecules and another between each of the hydrogen molecules and the dipole inducing negative hydrogen ion. Because of the nonsymmetric geometries of the newly found equilibrium structures, it is surprising perhaps that there are no bond paths between any of the hydrogen ligand molecules that would account for the “symmetry breaking” of the presumed VSEPR geometries. The bond paths that result from the calculations of this paper are analyzed further and more completely in a separate paper.20 On the other hand, DFT calculations are known to not describe accurately those forces that are very weak Van der Waals forces, as might be expected to occur here between hydrogen molecules.
H13
0.5495
The X3LYP functional21 has been put forward as representing these weak forces better than the B3LYP functional that has been used here. A search for any nascent bond paths between the hydrogen molecules was undertaken as a test of X3LYP accuracy. The results obtained were qualitatively similar to those shown here for the X3LYP functional, and thus no bond paths between hydrogen molecules occurred. However, independently of DFT, in a separate paper, using configuration interaction, which is known to be highly accurate if computationally demanding, we have found bond path connections between and among the hydrogen molecules of these clusters.22 Energetics related to the molecular structures of this paper are collected together in Table 1. The total molecular energies are listed. The HOMO 12449
dx.doi.org/10.1021/jp203913n |J. Phys. Chem. A 2011, 115, 12445–12450
The Journal of Physical Chemistry A and LUMO molecular orbital energies and the energy gap between them have been calculated. Also listed are the binding energies for each molecular structure. The VSEPR geometries have energies that are close to those of the true minima configurations, but as we have emphasized, they do not conform to that condition of minima that the vibration frequencies should all be real. The binding energies of the ion induced dipole molecules increase with the number of hydrogen molecules contained in their molecular structure. We have calculated here only as far as n-odd = 13 in the series of Hn ion induced dipole bonded molecules. It is known from mass spectroscopy experiments23 that the positive analogs of these molecules, Hn+, are known to exist for all values of n-odd from 3 up to at least 99. The corresponding negative ions have not yet been experimentally discovered. But surely they can be, as suggested now by the new DFT equilibrium geometries of this paper.
’ AUTHOR INFORMATION Corresponding Author
*L.H.: phone, (202) 404-2138; fax, (202) 404-7546; e-mail:
[email protected]. L.M.: phone, (212) 772-4016; fax, (212) 772-5332; e-mail,
[email protected].
’ ACKNOWLEDGMENT The research reported in this article was supported by the Office of Naval Research. L.M.’s studies were funded by the U.S. Naval Research Laboratory (project no. 47203-00 01) and by a PSC CUNY Award (project no. 63842-00 41). C.F.M. acknowledges Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI), and Mount Saint Vincent University for funding. We thank Douglas J. Fox of Gaussian, Inc. for invaluable discussion.
ARTICLE
(17) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: New York, 1990. (18) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, € Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; S.; Daniels, A. D.; Farkas, O.; Fox, D. J. Gaussian09; Gaussian, Inc.: Wallingford, CT, 2009. (19) Keith, A. T. AIMALL, version10.09.12; TK Gristmill Software, aim.tkgristmill.com, 2010. (20) Matta, C. F.; Huang, L.; Massa, L. J. Phys. Chem. A 2011, 10.1021/ jp203973d. (21) Xu, X.; Goddard, W. A., III. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 2673–2677. (22) Huang, L.; Matta, C. F.; Massa, L. “Hydrogen Clusters (H2n+1 ) as a Possible Chemical Source of Diffuse Interstellar Bands: Spectra, Energies, and Structures from Configuration Interaction Calculations.” (to be submitted). (23) Clampitt, R.; Gowland, L. Nature 1969, 223, 815–816.
’ REFERENCES (1) Harrison, S. W.; Massa, L. J.; Solomon, P. Nat. Phys. Sci. 1973, 245, 31–32. (2) Harrison, S. W.; Massa, L. J.; Solomon, P. Chem. Phys. Lett. 1972, 16, 57–59. (3) Harrison, S. W.; Massa, L. J.; Solomon, P. Astrophys. J. 1974, 189, 605–607. (4) Harrison, S. W.; Massa, L. J.; Solomon, P. J. Chem. Phys. 1973, 59, 263–266. (5) Roothaan, C. C. J. Rev. Mod. Phys. 1951, 23, 69–89. (6) Hirao, K.; Yamabe, S. Chem. Phys. 1983, 80, 237–243. (7) Rayez, J. C.; Rayez-Meaume, M. T.; Massa, L. J. J. Chem. Phys. 1981, 75, 5393–5397. (8) Sapse, A. M.; Rayez-Meaume, M. T.; Rayez, J. C.; Massa, L. J. Nature 1979, 278, 332–333. (9) Parr, R. G.; Yang, W. Density-functional theory of atoms and molecules; Oxford Univ.ersity Press: Oxford, U.K., 1989. (10) Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2011, 7, 669–676. (11) Dierksen, M.; Grimme, S. J. Chem. Phys. 2004, 120, 3544–3554. (12) Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100. (13) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. (14) McLean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639–5648. (15) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294–301. (16) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265–3269. 12450
dx.doi.org/10.1021/jp203913n |J. Phys. Chem. A 2011, 115, 12445–12450
