Marie-dose Huron and Pierre Claverie
Calculation of the Interaction Energy of One Molecule with Its Whole Surr Applic~ti~n to) Pure Polar Compounds Marle-Jos6 Huron* lnstifut Francais do Petrole, 92502 Rueil-Malmalson,France
and Pierre Claverie Institut de Biologie Physico-Chimique, Laboratoire de Chimie Quantique, 75005 Paris, France (ReceivedFebruary 20, 1973; Revised Manuscript Received May 7, 1974)
A numerical method is proposed for calculating the total interaction energy of a polar molecule with its entire surrounding. The total dispersion and repulsion energies are expressed as products of energy integrals (calculated numerically outside the effective volume 121 of the molecule) by calibration constants Kzj characterizing packings of spheres. The electrostatic interaction energy is calculated numerically, considering the solvent as a homogeneous continuum medium with a macroscopic dielectric constant D, and developing the reaction potential on basis of harmonic functions. To take into consideration the nonperfect representat ion by a continuum model of the behavior of the molecules in the immediate vicinity of the central molecule, an empirical anomaly factor f was introduced. The model when applied to pure polar compounds gives satisfactory results in a comparison with the vaporization energies, using a different factor f for associated and for nonassociated molecules.
I. Introduction In the first article of this series,l we proposed a method for calculating the dispersion and repulsion energies, which enabled us to treat the case of nonpolar molecules, and we announced the treatment of the case of polar molecules. This case requires an evaluation of the electrostatic energy and, in principle, of the polarization energy. The difficulty in creating a discrete model of the solvent for calculating the electrostatic energy is much greater than for dispersion and repulsion, because the respective positions of the molecules play an essential part. For dispersion and repulsion, we were able to use a model essentially based on the packing of spherical atonis present in the solvent molecule. The atoms were treated independently and their proportions in the molecule were fixed only by the chemical formula of the solvent. This was possible because the dispersion and repulsion energies are sums of atom-atom contributions that are respectively all attractive and all repulsive. Such a decomposition has no meaning for the electrostatic energy since the molecule acts as a whole. This is why we have developed, in the second article2 of this series, a method where the solvent is considered as an infinite continuum dielectric. The important role of the positioning of the molecules has another elonsequence; Le., a calculation of the work required to go from a hypothetical completely discharged state to a charged state gives the variation of the free energy and not the variation of the internal energy. This is in agreement with our analysis in Appendix C of our first article,l because the passage from an uncharged state to a charged state is accompanied by a change in the spatial distribution of the molecules, including the position of the molecules. This third article deals with the applying of the first two articles to polar molecules. In the first part, we analyze a way of calculating thie energy terms. In the second part we give the resultb we have obtained for calculating the energies between identical polar molecules (alcohols, ethers, The Journal of Physical Chemistry. Vol. 78,
No. 18, 1974
and others), and we make a comparison with vaporization energies.
11. Calculating the Different Parts S% the Energy
Dispersion and Repulsion. The energies of dispersion and repulsion are expressed as the sums of products of energy integrals by calibration constants K;,1 E, + E, =
Ec i
J
''
e i j da v
{Kij Kolv
+
E
where i is a summation on each atom bf the solute molecule and j a summation on each distinct species of atoms of the solvent. The potential of atom-atom interaction eij has the folllowing expression
e i j = eijd
+
eijr = k,kj{-0.214(4~i~j)3/Yij5 +
1
47,000~i c j exp[- 1 ~ . 3 5 ~ ~ ~ / '1( , ~ R (2)~ ~ ~ ) * ~
R; and Rj are the van der Waals radii of atoms i and j . The factor ki OF kj, a characteristic of the species of atom, i s introduced to increase the depth of the binary potential eij for atoms other than C and I(. Since Kitaigorodskii's formula is expressed with the universal constants as a function of z = r/Ro, the minimum of the curve representing eij(z) is independent of the nature of the atoms. However there is no reason for the physical value of this minimum to be the same for every pair of atoms. Indeed, Mitaigorodskii, et ~ 1 .had , ~ to introduce a formula of the type (2) in order to reproduce the sublimation energy of CQ2, with factor hi greater than I for the oxygen (ko = 1.36). In the same way the potentials used for a nitrogen atom by Ahmed and Kitaigorodskii4 correspond to k~ m 1.18.
Total Interaction Energy of a Polar Molecule
1863
The results for hydrocarbons show that kc = 1 and kH = 1 are suitable. For oxygen and nitrogen, we adapted k o and k N by calcu tatrng the vaporization energy of molecules 0 2 and N2 with ey 1, assuming that for these molecules there are only dispersion m d repulsion energies. Thus we obtained k o = 1.242 and kN = 1.145. For sulfur, the adaptation was made for the CS2 molecule and gave the value K S = 2.40. In the repulsion term, the multiplicative factor c, is introduced for each atom
ci = 1 - q i / N v ( i ) = P i / N v ( i )
(3 1
q Lis the net charge of atom i, NJi) is the number of electrons in the valence shell of the neutral atom I , and P, is the electronic popLllation of the valence shell of atom i. The variation of cI, which is proportional to PL, is a way of taking into iiccount the reduction in repulsion for short distances actcompanying the decrease in the electronic density. This effecf is important for atoms H with high positive charges (group!; OH, NH). Without being sufficient to completely represent the interaction of the hydrogen bond, this change improves the description of liquids with such bonds. In Table IA, in Appendix A, we give the values of the quantities characterizing the packings of spheres.* They enable us to calculate the calibration constants Kyd and KLJrfound in eo 1. Electrosta tic Energy. The expression of the electrostatic free energy for 3 molecule is2
type of molecule if we want the theory to predict values of energy and not just to reproduce already known experimental values. We proceed by fitting the coefficient f using various experimentally known cases. The success of this fitting may be checked in the following way: after choosing several types of molecules (for example, molecules associated by hydrogen bonding and nonassociated molecules), we must find a single value or a narrow range of values o f f for each type, so that we can reproduce the results concerning all the molecules of a given type. Factor f certainly varies with the temperature,; and this makes the situation still more difficult, since some values will also have to be assumed for its derivative. Consider what becomes the energy in the presence of factor f. We call A F the physical free energy and AF the free energy calculated according to the continuum model without any correction. Similarly we write AE and AE as the corresponding internal energies given by the Gibbs-Helmholtz relation
AE = A F - T a (AF) aT AS A F = fAF eq 5 becomes AE = f A F - T
a (feF)
A E = f A E - TAF--
dT
Using
where W,' is a harmofiic basis function in the volume QI of the molecule, and a coefficient expressed as a function of the dielectric constant and of surface integrals containing the inside and outside basis functions We1 and W,E.2 N A is the number of atoms in the central molecule, N I the the reaction number of basis functions used in 01,and potential. In eq 4 the coefficient */2 cannot be considered as a rigorous one as eplained in Appendix B. We have introduced, in the electrostatic energy, a corrective factor f, because in a continuum model the macroscopic dielectric constant does not give a perfect representation of the statis1,ical behavior of the solvent molecules in the immediate vicinity (of the solute molecule. Bjerrum5 proposed to attributing this fact to the observation of several anomalous effects concerning electrolyte solutions, in particular the heat of dilution. He then introduced an anomaly coefficient f by which the free energy, calculated by using the macroscopic dielectric constant, should be multiplied. Bjerrum pointed out that the effect of this coefficient could be simulated by using a suitably modified value of the ionic radius, and Fowlerfi retained this suggestion. However this procedure is noidvery feasible when dealing with molecules instead of spherical ions, because it would be necessary to find a systematic waysfdefining a modified van der Waals radius for each atom i n the solute molecule. We did not attempt to find such a systematic procedure, and we simply used an overall anomaly coefficient f for each solvent molecule. This faletor milst be fixed a priori according to the
eq 7 becomes
In the application of section 111, we use eq 10 rather than eq 8, assuming for a given type of molecules a constant value for d log fld log T. There is no rigorous theoretical argument to justify this choice. However some qualitative arguments make a constant or almost constant value more plausible for d log fld log T than for df/dT. d log f/d log T = k gives flfo = (T/To)k,which appears more suitable than a linear dependence o f f on T. For example, i f f must decrease with T,5a linear dependnce would imply that f vanishes and becomes negative for some temperatures T; while the power dependence with a negative k gives a vanishing value only for T infinite. Bjerrums considered that d log f/d log T (and not dfldT) could be a constant. Polarization Energy. A part of the polarization energy, namely, the energy of the solvent molecules polarized by the solute molecule, is already included in the "electrostatic energy" calculated according to the continuum model.7 This is due to the fact that the dielectric constant takes into account the polarizability of the solvent molecule (temperature independent part of the dielectric constant). The other part of the polarization energy, namely, that of the solute polarized by the solvent, was not calculated in the present work. It is well known that the polarization energy is significantly smaller than the other terms of energy. The Journalof Physical Chemistry. Voi. 78. No, 18, 1974
1864
The Journal of Physical Chemistry, Vol. 78, No. 18, 1974
Marie-Jose Huron and Pierre Ciaverie
Total Bnteractior'i Energy of a Polar Molecule Since the electrostatic energy is not evaluated with sufficient accuracy to make relevant an accurate evaluation of the polarization energy, we do not see any use in completing the evaluation of the polarization energy. 111. Results The model described above has been applied to the study of some pure polar compounds. The calculated energy Ecalcdi s compared with the vaporization energy Eexp
Energy E,:l is not multiplied by the factor y2 in (11))because is idrr?ady iincluded in A F (see eq 4) and consequently in AE and Eel. For the heteroatoms, we used the following van der Waals radii*: Ro = 1.50 8, RN = 1.55 8, in NO2, RN, 1.60 in N2, Rs = 1.80 A. A correction of 7% was made on the total dispersion energy to take into account the nonadditivity of this energy.l The net charges required for calculating the electrostatic energy are calculated using the Del Re method for alcohols, esters, aniline, arid phenol; the CNDO method for acetone, nitromethane, and nitrobenzene; the dipole moments for H20, S02, and SHzpand the quadrupole moments for COZ and 6S2.9 To develolp the potential as indicated in ref 2 for calculating Eelt we tised A', = 16 basis functions in volume QI; and N E = ~ N basis A functions in volume QE (i.e., 4 basis functions per atom). The vaporization enthalpies, the volumes, the dielectric constants, arid their derivative were taken from ref 10-14. The values we used are indicated in Table I. The molecules we have studied can be separated into two groups: (1) the group of associated molecules (containing OH or NH2) for which the values at 25' f = 1.63 and d log f/d log T = -0.50 seem suitable; and (2) the group of nonassociated molecules for which at 25' we took f = 0.33 and d log fld' log I' = -0.50. The energies calculated with these values are correct except for water (see Table I). For molecules bCH3)2CO, CH3N02, and C G H ~ N Owe ~, took into account the contribution of the atomic dipoles in calculating the electrostatic energy (see eq 48 in 2). The electrostatic energy calculated with the CNDO charges alone would be inuch Power (1.632 kcal/mol instead of 4.662 kcal/mol for C 3 & J 0 2 at 25") 1.682 kcallmol instead of 4.289 kcal/mol for C#&,N02 a t 25'). Hence the atomic dipoles play an important part in the CNDO representation of the electrostbtic distribution, and we have to take them into account. For nitrobenzene tine dipole moment with the CNDO charges and the atomic dipoles is too high. This could be an explanat,ion for the excessive values found for the calculated energies. The only marked disagreement with experiment occurs for water. It could be argued that the association between water molecules is rather special with respect to other associated liquids, and that water would deserve a particular fitting of the anomaly coefficient f. Before resorting to such an explanation, however, it would be necessary to improve the electrostatic representation of the molecule; we used indeed a three-charge model fitted to reproduce the dipole moment, but the higher multipole moments are not accu-
1865
rately represented in this way,15 and this may produce a nonnegligible error for a molecule which, simultaneously, is highly polar and has a small size. Therefore, it cannot be excluded that the disagreement may be decreased without an ad hoc modification of the anomaly coefficient f.
IV. Conclusion The method developed in this paper seems suitable for the treatment of pure polar compounds. It takes into account the shape of the molecules and gives, (qualitatively, the variation of the energy with temperature. The complete calculation on a CDC 6600 lasts about 27 sec for a triatomic molecule, 89 sec for CH3N02, and 300 sec for CsH5N02. T o treat very large polar molecules, the simplifications suggested in paper 112 should be tried in order to obtain short computation times. The application to dilute solutions could be tried, using the table given in Appendix A, and following the lines indicated in the conclusion of paper I; thus, the solution energy of a gas (the solute) would appear as the difference between the solute-solvent interaction and the cavitation energy (necessary energy to create a cavity in the solvent to accomodating the solute molecule). Appendix A. Packing of Spheres To calculate the calibration constants KIJdand KLJr,l we need the characteristics of the packing of hard spheres with radius R, around a sphere with radius R,. Table IA gives these characteristics for many values of ratio R,/R,. The packing of spheres is generated by a computer program which gives the coordinates of each sphere center rel'ative to the center o f the central sphere plus the value of r k o and [(4RLRJ)1'2/rko]6, where rho is the distance from the center of the central sphere to the center of a sphere K belonging to the hth shell. If we represent & = I m (R,olr~o)6as a function of ratio R,/R,, we find several discontinuities in the curve near the following values: 0.90, 1.20, P.46,1.70. These discontinuities correspond to a sudden change in the value of nl which is the number of spheres belonging to the first shell, i e . , spheres for which rl0 = R, 4- R,. This phenomenon is related to the existence of more or less big holes left between the tangent spheres in the packing. We need nl, n2, rl0, and r2'3 to calculate the repulsion energy.l n2 includes only the spheres of the second shell that contribute to the repulsion energy. rzo/Rv0 is a mean value of the real distances from these n2 molecules to the center of the central sphere. r~oIRLlo is such that the repulsion of these n2 molecules, all a t the same r2O distance, is equal to the repulsion of these n2 molecules at their real distances. All the characteristics given here will be useful when dealing with the general case of a central molecule surrounded by molecules of another type (mixtures). Appendix B. Relation between Reaction Potential and Total Electrostatic Free Energy Following the principle of the Debyedetermining free energy,6J6we may express AF with eq (212 of Appendix C in paper I1
where ( is a charging parameter, i.e., we consider the hypothetical state where each charge has its value multiplied by The Journal of Physical Chemistry, Vol. 78, No. I S , 1974
Marie-Jose Hiiron and Pierre Claverie
1.oo
0.565 0.647 0.667 0.678 0.684 0.687 0.697 0 .TO6 0,741 8.774
1.10 1.20 1.20 1.30 1.10 1.185 1.20 1.185 1.20 I .20 1.50 1.50 1.55 I .50
0.800
0.833 0.847 0.876 0.882 0.889 0.904 0.912 0.937 0 941 0.960
1.60
1.60 1.55 1.50 1.60 1.70 I .50 1.77
I
0.968
0.983 1.000 1.ooo I ,017
1.80
1,033 1.041 S.062 1.067 1.097 1.306 I .133 I . 142 1 180 1. ,200 1.250 1.292 1,350 1.41’7 1.435 1.454 1.461 1.475 1.500 1.545 I 600 1.650 1,700
I .55 1.77
I .770
1”77
1.70 1.60 1.70 1.77 1.70 1.77 I .77 1.80 1.50
1.55 1.60 1.70 I .70
1.60 1.90 1.77 I .80 I .70 I .60 I .65 1.70
I
1.77 1.70 1.80 1.77 1.90 1.60 1.70 1.70 1.60 1.55 1.50
1.so
1.77 1.77 1.70 1.80 1.77 1.70 1.60 1.70 1.77 1.55 1.80 1.77 1.50 1.70 1.60 1.50 1.55 1.60 1.50 1.55 1.50 1.50 1.20 1.20 1.185 1.20 1.185 1.10
1.30 1.20 1.20 1.10
1.oo 1.oo 1.00 1.oo
5.873 6.934 7.175 6.643 6.715 6.754 6 370 6.971 7.460 8.045 8.515 9.245 9.614 10,492 10.733 10.963 13.491 13.537 13.707 13.728 13 .850 13.898 13.993 14.1OOa 14.454a 14.217 14,323 14.377 14.540 14.577 14.872 14 ,968 15,149 15 .252 15.860 16.684 16.985 17.182 17.412 18.216 18.584 19.179 20.604 20.629 20,752 20.927 21.414 21.653 23.091 23.515
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
1.041 1.024 1.021 1,019 1.018 1.018 1.016 1.015 1.011 1.008 1.006 1.004 1.003 1.002 1.002 1.001
6 12
12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 14 14 14 14 14 14 14 18
3 3 3
3 3
3 6 6 6 6 0
1.ooo 1.ooo
0
1.000
0 0 0
1.000 1.ooo 1.ooo
8 2 4
1.469 1.453
1.ooo 1.ooo 1.000
4
3
1.436 1.381 1.364 1.317 1.308 1.240 1.218 1.154 1.134 1.839 1.367 1.310 1.317 1.275 1.110
3
1.lo6
2 2 2
1.029 1,197 1.187 1.167 1.122 1.064 1.087
1.021 1.024 1.028 1.031 1.035 1.041
22 22
3
1.001
1.018 1.os9
18 18 18 18 18
1 P 3
1.001
1.000 1 .ooo 1.001 1.001 1.002 1.002 1.003 1.004 1.006 1.008 1.011 1.015 1.016 1,018
12
1.137 1.074 1.067 1.347 1.331 1.331 I ,316 1.301 1.245 1.201 1,154 4.155 1.133 1.090 1.080 1.070
1
2 2 4 1.
I 2 1 1 I
3 3 4
4
2 2 2 3
1 I
1.211
1.063
14.100 is th(s value given by our calculations; 14.454 is the value given by Hirschfelder for a close packed crystal (see ref 11 in ref I);in lhe calculations we used 14.454.
E. U((,S2)is
the; total interaction energy of the assembly of
N molecdesi. In U we separate the electrostatic part Ue,
u u v y . .. ‘ $ q i ’ l . . . x q j v . . ) = I
[Xuuvel(. . . q i u . . .q;.
which can be written ~k~,,~’(.
. . q i u . .. q y . . . )
(B2)
Therefore
0 *lJ
where p and v are summations on the molecules of assembly; qL@(with i .= 1 . . np) and qJu(with j = 1 . . . n,) are respectively the net charges of atoms i and j belonging to the molecules p and if (continuous distribution or multipole charge distriibution could be assumed as well). The essential point is that the electrostatic interaction between molecules p and Y varies linearly with the set of charges of evory molecule, Le. The Journal of Physical Chemistry, Vol. 78, N o . 18, 7974
5
c Cu,,By. u v u#v
. 1 5
(B3)
Total lnteractior Energy of a Polar Molecule The potentid LTNuCluel can be expressed as a function of the electrostatic potential cpu(r,Q’)created by all the molecules except v at point i (instantaneous potential depending on all position and angular coordinates Q’ of all molecules except u )
1867
the intermediary of the dielectric constant used in the continuum model, the exact value of the integral in (B9) is not %(&j))
1
To get
%, it would be necessary to have
(cp(Fi))E.
E
= 1. However the reaction potential \kI(D(())is a rather insensitive function of D when t varies from 0 to 1. This can be verified for a dipole in a sphere, assuming I)([) = Do (D Do)t2. (cp(ij))+
+ -
61, and the weight function fw((,n)(see ref
1, CS), we obtain
The overestimate due to the use of 1/2 is at a maximum for small values of D (Le., for negligible electrostatic energies) and significantly decreases when D increases. For these reasons, we shall keep the factor Yz. At any rate, the anomaly factor f we have introduced to correct the continuum model can compensate for the overestimate resulting from this approximation.
References and Notes
If we deal with n pure liquid, all molecules are identical, and the mean value OS the potential is the same for all molecules. Thus we maiy drop index v and write the mean value (p(Fj))~. All the terms inside the brackets are equal; therefore Z V [ E j ]== NZj, where N is the total number of mole-
Compare eq 139 with eq 4 in this paper. The reaction potential @‘(?,)I is nothing but ( ( p ( F J ) ) ~ , and our continuum model provides a way of approximately obtaining this poowewer, since ( p ( i J ) ) its a function of t by
(1) M. J. Huron and P. Ciaverie, J. Phys. Chem.. 76, 2123 (1972). in eq C-6 to C-12 in Appendix C, 5 after a parentheses is an index which is incorrectly placed: it shouid read ( ) t instead of ()(; Table Ii, for benzene should read A = 1.306 and 1.352 instead of 1.332 and 1.376, respectively. (2) M. J. Huron and P. Ciaverie, J. Phys. Chem., 78, 1853 (1974). (3) A. 1. Kitaigorodskii. K. V. Mirskaya, and V. V. Nauchtei, Sow. Phys. Crystallogr.., 14, 769 (1970). (4)) N. A. Ahmed and A. I. Kitaigorodskii, Acfa Crystalbgr. Sect. 13, 28, 739 (1972). N. Bjerrum, Trans. Faraday Soc., 23, 445 (1927). R. H. Fowler, “Statistical Mechanics,” Cambridge Universiry Press, New York, N. Y., 1966, p 546. B. Linder, Advan. Chem. Phys., 12, 225 (1967); see p 267. A. Bondi, J. Phys. Chem., 68, 441 (1964). (a) 6.Del Re, J. Chem. SOC., 4031 (1958); (b) J. A. Pople and M. Gordon, J. Amer. Chem. SOC.,89,4253 (1967). (a) “international Critical Tables,” McGraw-Hill, New York, N. Y., 1926: (b) API, 44 Tables, “Selected Values of Properties of Hydrocarbons and Related Compounds,” American Petroleum InStltute Research Project No. 44, Thermodynamics Research Center, Texas A & M University. J. Polak and G. C. Benson, J. Chem. Thermodyfl., 3,235 (1971). J. G. Helpinstii and M. van Winkle, Ind. Eng. Chem., Process Des. DeweC op., 7, 213 (1968). P. Pascal, “Nouveau traite de chimie minerale,” tomes V et Xil, Mason et Cie, Paris, 1956. A. A. Maryott and E. R. Smith, Naf. Bur. Stand. 0. S., Circ., 589 (1958). D. Eisenberg and W. Kauzmann, “The Structure and Properties of Water,” Clarendon Press, New York, N. Y., 1969, Chapter 1, sections 1-2
H. Eyring, D. Henderson, B. J. Stover, and E. M. Eyring, ”Statistical Mechanics and Dynamics,” Wiley, New York, N.Y., 1964.
The Journal of Physical Chemistry, Vol. 78, No. 18, 1974
