J. Phys. Chem. 1989, 93, 8263-8270 varying, radially directed electric field. The oscillations in these cases are due to dielectrophoretic forces, whereas ours are due to the electrocapillary effect. However, they also observed changes in the mode of the oscillation as the intensity of the electric field was changed. A stability analysis of the resting state of the mercury drop, leading to a selection of the mode most likely to evolve, can be undertaken" following the methods outlined by Jones et aL9
5. Conclusions The conditions under which the mercury beating heart phenomenon can be adequately described in terms of the electrocapillary effect, without reference to surface film formation, solution oxygen, or surface contaminants, have been established. This has been achieved by replacement of the chemical driving forces, ~~~~~
~
~
~~~
(1 1 ) Laidlaw, W. G.; Olson,J. To be
published.
8263
studied by Keizer et al. in their earlier extensive and elegant investigations, with an applied oscillating electrical potential. This has led to the observation that the modes of oscillation of the mercury drop can be controlled by selection of the mass of the drop and the frequency and amplitude of the applied potential and that transitions between modes of oscillation can be made to occur in a reproducible and reversible way.
Acknowledgment. The authors acknowledge the support of NSERC operating grants (W.G.L. and V.I.B.) and support from a STEP grant (C.U.). More importantly, it is a great pleasure to note the advice and assistance of our colleague, Hanna Elzanowska. Brian Kipling first brought this phenomenon to our attention, and his questions concerning its use as a demonstration prompted the investigation. Thanks are also due to Laura Glavina, for technical assistance. Registry No. Hg, 7439-97-6; 02,7782-44-7; NaOH, 1310-73-2; Na2S04,7757-82-6; PO4)-, 14265-44-2; K2S0,, 7778-80-5.
Transferability of Molecular Distributed Polarizabilities from a Simple Localized Orbital Based Method David R. Garmed and Walter J. Stevens* Center for Advanced Research in Biotechnology, National Institute of Standards and Technology, Gaithersburg, Maryland 20899 (Received: January 25, 1989; In Final Form: May 9, 1989)
A distributed model of molecular electric polarizability is presented in which the polarization of the charge density from each localized molecular orbital (LMO) is represented by a point dipole polarizability tensor located at the LMO charge centroid. Specifics of the generation of the polarizability tensors by finite-field Hartree-Fock SCF calculations are given, and the application of the model to classical calculationsof molecular interactions is briefly discussed. Large basis set calculations on a number of small molecules demonstrate that transfer of bond and lone pair polarizabilities and centroids among various molecules is possible, without explicitly invoking intramolecular interactions between induced dipoles. A catalog of these transferable LMO properties is developed for a number of different types of bonds and lone pairs. A set of molecular dipole polarizabilities is then reproduced from the catalog as a successful test of this idea.
Introduction A model is presented for distributed, transferable molecular electric dipole polarizabilities developed from finite-field S C F calculations with Boys orbital localization. The model for any system consists of a set of easily calculated independent point dipole polarizabilities, located at the localized molecular orbital (LMO) charge centroids. We use this as an approximate calculational tool to predict polarizability-based effects from nonuniform applied fields, particularly the electrostatic fields of interacting, but nonbonded, molecules. The results are roughly as successful in predicting the effects as the more rigorous method developed by Stone' which is, at present, effectively restricted to systems characterizable by linear arrays of polarizable points.2 In large systems the present method might to some extent fill in to make possible the calculations and theoretical studies performed for small systems using Stone's distributed polarizability model. The predictable effects include the electric polarization part of two-body interactions and the induced dipole moments. Also, the bulk of N-body interactions with N greater than two is, in many cases, primarily a polarization phenomenon, predicted accurately by this modeL3 The required finite-field Hartree-Fock calculations are approximately as efficient as the equivalent coupled perturbed Hartree-Fock (CPHF) calculations of the molecular polarizability, if certain simple measures are taken. Extending this model to large molecules requires some further approximation because ab initio calculations of sufficient quality NAS-NRC Postdoctoral Fellow. *Author to whom correspondence should be addressed.
to predict accurate polarizabilities are presently impossible for systems of many more than 15 atoms. The results enumerated below illustrate that it is often possible to use the simple approximation of transferring the distributed polarizabilities of functional groups from small molecules into a larger system. The combined polarizability model would then be an accurate predictor of the overall molecular polarizability of the system and of the local polarizability of importance in intermolecular interactions. This article is mostly concerned with this transferability, rather than with results on molecular interactions from the model which are discussed in greater detail el~ewhere.~ Other theories relating to distributed polarizabilities include the localized uncoupled perturbed Hartree-Fock (UCPHF) quantities used by Karl~trom.~ These are expected to be generally of somewhat poorer quality than the finite-field S C F polarizabilities used here since they do not sum to the molecular dipole polarizability. However, we have not attempted a general comparison of localized UCPHF and C P H F polarizabilities for this work. Another method, capable of handling polarization effects in nonuniform fields, attempts to saturate a single-center molecular polarizability by adding in any number of higher moment polarizabilitie~.~The problems of asymptotic convergence and (1) Stone, A. J. Mol. Phys. 1985, 56, 1065. (2) Fowler, P. W.; Stone, A. J. J . Phys. Chem. 1987, 91, 509. (3) Work in preparation. (4) Karlstrom, G. Theor. Chim. Acta 1982, 60, 535. (5) Liu, S.-Y.; Dykstra, C. E. Chem. Phys. 1986, 107, 343
This article not subject to U S . Copyright. Published 1989 by the American Chemical Society
8264 The Journal of Physical Chemistry, Vol. 93, No. 25, 1989
spherical convergence radii, which stimulated work on distributed methods, would seem to limit this to certain small molecules, although this potentially may be a computationally more efficient scheme than distributed polarizabilities, allowing also the ability to easily introduce nonlinear polarizations within one computational framework.6 The second version of the mutually consistent field developed by Otto includes a classical representation of molecular polarizability.’ This consists of a distribution of independent dipole polarizable points, designated as pseudopolarization tensors, determined from perturbations in the ab initio electrostatic potentials induced by applied fields. The details of the development and the placement of points are quite different from the method presented here, however. The total energy expression used incorrectly does not consider the change in internal energy of a molecule upon polarization, but only the change in the electrostatic interactions between molecules. Various usable semiempirical theories of molecular polarization have been proposed. One class of these relies on intramolecular induced moment interactions to modify isotropic point dipole polarizabilities, generated empirically and assigned to various atomss These can be extended by using more complex constructions of the charge density but with the same physical basks9 Other semiempirical studies have shown that polarizabilities assigned to classes of bonds and lone pairs might give approximately the same results without explicitly considering the interactions between induced dipoles.IO These polarizabilities presumably should be comparable to those defined by the ab initio results of this work. The principal advantages of a method derived from calculation are that the microscopic details can be verified as being correct by performing sample SCF calculations on whatever system is being modeled. Some insight into the physical nature of the polarization of molecules might be gained and improvements in models developed from studies using this ab inito decomposition, without the ambiguity introduced by having fitting parameters.
The Polarizability Model The method presented here is similar to that of Maestro and Moccia,” who proposed writing the C P H F dipole polarization of a closed-shell molecule as a sum of contributions from individual LMO’s using Boys localization criteria.’* An equivalent, but much simpler, finite-field method of obtaining the resulting localized polarizabilities is outlined below. The equivalence of the two algorithms has been verified by repeating the calculation on methane described in ref 11. We also add to their treatment by assigning the local polarizabilities to the L M O centroids (centers of charge) and then using these for polarization calculations involving nonuniform applied fields. The molecular polarizability can be calculated by the finite-field method from a,, = pFd Fr
where pFd is the s component of the dipole moment induced when a perturbing uniform electric field of magnitude Fr is applied along coordinate direction r. The field magnitude must be small enough so that pFd is a linear response to the perturbation to within a certain desired accuracy (see discussion below). For a singleconfiguration wave function, the one-electron density matrix and expectation values are given by sums of contributions from individual orthogonal orbitals. Making this decomposition for the perturbed and unperturbed molecular dipole moments in iind allows (6) Dykstra, C. E. J . Compur. Chem. 1988, 9, 476. (7) Otto, P.Inr. J . Quantum Chem. 1985, 28, 895. (8) Thole, B. T. Chem. Phys. 1981, 59, 341. Applequist, J.; Carl, J. R.; Fung, K.-K. J . A m . Chem. SOC.1972, 94, 2952. (9) Olson, M. L.; Sundberg, K. R. J. Chem. Phys. 1978, 69, 5400. (IO) Denbigh, K. G. Trans. Faraday SOC.1940, 36, 936. Miller, K. J.; Savchik, J. A. J. Am. Chem. SOC.1979, 101, 7206. ( 1 1 ) Maestro, M.; Moccia, R. Mol. Phys. 1975, 29, 81. (12) Boys, S. F. Rev. Mod.Phys. 1960,32,296,300. Boys, S . F. Quanrum
Ed.; Academic Theory in Atoms, Molecules and Solid Stare; Lowdin, P.-O.,
Press: New York, 1966; p 253.
Garmer and Stevens for the definition of a polarizability to be calculated from each localized orbital 4: by (ai)rs
= Mfld)s/Fr
(2)
where the induced dipole moment in the LMO charge density is given by
(3) using also the perturbed LMO (= 4: 4;). The quantity (&r&) is the s coordinate of the LMO charge centroid which is easily computed from the dipole operator matrix. Thus, the LMO induced dipole moments and polarizabilities are found from the displacements of the centroids by an applied field. The pertubation by a uniform electric field is introduced into the electronic Hamiltonian operator for a specific calculation also by using the dipole operator matrix. In order to complete the polarizability matrices, a series of four S C F runs must be performed with 0, F,, Fy, and F, fields applied in an arbitrary Cartesian coordinate system. Note that using canonical orbitals or other specified rotations of the canonical orbitals, rather than the Boys LMO’s, leads to differing sets of calculated polarizability tensors. These variously defined tensors would have to sum to the molecular polarizability by definition, as would the expectation quantities related to polarizations in a uniform field, such as the induced dipole moments. However, leaving aside considerations of possible arbitrary orbital rotations due to degenerate orbitals, the method used in this work is a unique and well-behaved decomposition for any molecular system. Several of the reasons why a localized orbital method is preferred are outlined below. The dipole polarizability of the LMO’s was chosen to be the basic element of this distributed polarizability model because a dipole is the lowest order induced moment defined for LMO charge densities. The long-range electrical interaction of a charge distribution is determined by its multipolar electrical moments, and the best first approximation usually comes from the lowest order moments. The LMO’s are naturally assigned two electron units of charge so that there is no perturbation at the monopole order by the field. Note that the induced dipoles in a model calculation give rise to higher moments for the molecule when these are evaluated at the molecular center of mass, as is usual. These polarizabilities are taken to be independent for classical interaction calculations involving a nonuniform applied field. For example, the induced molecular dipole moment is computed as
+
pYd =
C r=xy,r C Fr(e)(mi)rs
(4)
f
where F(rF) is the applied field at the centroid of LMO 4:. Thus, we do not try to model the intramolecular induced moment interactions as is done in several empirical models for the molecular p ~ l a r i z a b i l i t yand ~ ~ we ~ do not try to calculate them directly as is done in the distributed polarizability method due to Stone.’ Rather, in using the coupled Hartree-Fock localized polarizabilities independently, the interaction of any induced moment with the others is effectively approximated by the interaction taking place in a uniform applied field. This approximation has been evaluated in practical calculations, some of which are described below. The generally successful results are probably a benefit from using localized orbitals, giving a spatially separated description of the molecular charge density. Then, induced moment interactions tend not to be the most critical part of the polarization effect. There are several other reasons for choosing a localization method. The orbitals so produced are localized in space so that there are physically reasonable choices for locations in which to place a polarizable point. One of these, the centroid of the electron density contribution of the orbital, is used here. The location for such a polarizable point must be determined by some other criteria than just reproducing the induced dipole moments in a uniform applied field since these are invariant to an arbitrary translation of the molecule. Other criteria could be applied, such as reproducing the higher order induced moments in uniform external
Transferability of Molecular Distributed Polarizabilities fields obtained at little extra cost from the finite-field calculations. However, this would probably not yield unique methods nor lead to a useful transferability of properties. Another advantage of using localized orbitals and densities is that they are very often roughly in the shape of an ellipsoid sliced by one or two nodes. Zero-field and finite-field S C F calculations appear to indicate that this leads to a relatively rapid convergence of the higher static and induced moments for these densities, computed relative to the centroids. Thus, truncation of the induced moments for each LMO density at the order of dipoles introduces acceptable errors in most applications. Stone and Alderton13 explain in detail the relatively rapid convergence of their static distributed multipolar expansion as due to the compact semispherical reference charge densities. Since this appears to be an operationally equivalent situation, we will not repeat this important argument here. A change from the total molecular polarizability is that the polarizability tensors for individual LMO’s are not necessarily symmetric. This asymmetry is due to the asymmetric environment of the L M O charge densities in molecules. However, this introduces no significant extra theoretical or computational difficulty so that no effort has been made toward developing an effective symmetrization method. A driver for doing finite-field polarizability calculations using the existing S C F routines has been added to a local version of the HONDO program.14 This allows for calculations of field-induced moments of any order and of the perturbed LMO’s required for this work. Measures taken to ensure maximum efficiency with this code include the reuse of all molecular integrals in the perturbed calculations and the use of the converged unperturbed eigenvectors to start the perturbed SCF cycles. The resulting code runs approximately as rapidly as the CPHF routine in the standard HONDO package. Orbital localization is performed after each converged SCF run using the routine in HONDO which performs 2 X 2 rotations of molecular orbitals picked in a randomly chosen order, causing the orbitals to converge to a set of LMO’s. Obviously, in order to define well-behaved polarizabilities, it is required that the LMO’s from the finite-field calculations be connected to the zero-field LMO’s such that there would be no discontinuous jump in the orbitals if the SCF cycles and localization were to be carried out for increments of the field increasing from 0 to F i n magnitude. However, there generally exist several rotations leading to differing but locally optimal sets of LMOs. It is possible to guarantee that a connected set of LMO’s will be found, in most cases, by using the same order of 2 X 2 orbital rotations following each S C F run. Also, it is often necessary to use a simple algorithm identifying the nearest neighbors between the unperturbed and the perturbed centroids to make the proper assignment of perturbed with unperturbed LMO’s. This is required because different orderings of the LMO’s sometimes occur with our code, even when the localizations are properly connected. The matter of choosing field strengths and convergence tolerances for finite-field polarizability calculations has been discussed el~ewhere,’~ but the specifics are worth mentioning. The standard convergence criteria for the HONDO S C F program is a test on the maximum variation in any density matrix element from one cycle to the next. In these calculations this tolerance was set at lo-’ au, resulting in energy convergence of hartree or better. This convergence is 2 orders of magnitude tighter than is acceptable for most other applications. At this level of tolerance the induced moments, both for the individual LMO’s and for the entire molecule, are within 10” to IO-’ au of the precise values. The applied field strength of au thus leads to an uncertainty of IO-2 to au in calculated polarizabilities, from eq 2 . The possible effect of nonlinearity in the induced moments as a function of the applied field must also be counted as a source (13) Stone, A. J.; Alderton, M . Mol. Phys. 1985, 56, 1047. (14) Dupuis, M.; King, H.F.; Inf.J. Quantum Chem. 1977, 11, 613; J . Chem. Phys. 1978, 68, 3998. (15) Werner, H.-J.; Meyer, W. Mol. Phys. 1976, 31, 8 5 5 .
The Journal of Physical Chemistry, Vol. 93, No. 25, 1989 8265 TABLE I: Induced Dipole Moments (in Atomic Units) in the Optimal Formamide-Water Complex from ab Initio SCF and Classical Model Calculations“
orientation error in 102(pindJIpindl,deg
type of calculation
SCF with standard DZD/DZ (ECP) basis self-consistent LMO polarizability model noniterative LMO polarizability model self-consistent molecular u model
8.0 8.6 6.9 13.5
-1.1 7.4 -18.2
‘Details of the structure are given in ref 16. TABLE 11: Calculations of Quantities Related to the DiwleOuadruwle Polarizability Tensor As Defined in the Text
water methanol
acetic acid ethanediol propene acrolein
1.55 6.25 10.68 10.79 5.30 19.50
0.85 5.30 10.51 8.74 5.41 18.40
0.89 1.48 1.79 2.22 0.65 2.56
of error in the polarizabilities. Considering only the effect of the first nonlinear polarizability (the first hyperpolarizabiility, 6 ) in the usual expression for an induced moment from an applied uniform field Fr gives When eq 2 is used to calculate the finite-field polarizabilities, the observed values are thus related to the precise values by sobs = a r s + Y2PmsFr (6) rs Several calculations of components of B were performed on several small molecules using basis sets described in the next section, and comparisons in atomic units were made. The results seem to indicate that values of j3 have rather variable and unpredictable magnitudes when compared with a’s. However, none larger than the molecular a’s by more than an order of magnitude have been calculated. This translates to values of j3 components expected to be no greater than 100 au for molecules used in this work. Then eq 6 shows that a reasonable estimate of the possible errors from the first hyperpolarizability is less than au. The approximate balance between numerical errors and estimate bias only results from applied fields of roughly the magnitude used, and it minimizes the total expected error. Because the first hyperpolarizabilities are reduced in magnitude by approximate local symmetries, it is prudent to extend this kind of qualitative argument using a few available second hyperpolarizabilities and, thus, project that these will contribute errors of the order of au. Given these arguments and the goals of this work, it was not considered necessary to generally follow the expensive procedure of calculating induced moments at a variety of field strengths in order to make corrections to the values of aok.
Molecular Interaction Model Examples indicating that the distributed model given here is reasonably accurate and an improvement over a representation by a single molecular dipole polarizable point are given in Tables I and 11. A calculation of the electronic induced dipole moment in the globally optimal (cyclic double H bonded) formamide-water complex16was made at the SCF level and compared with various model calculations, giving the results listed in Table I. The model calculations are entirely classical and use accurate representations of the electrostatic fields of unperturbed formamide and water molecules. A calculation with a single dipole polarizable point at the center of electronic charge in each molecule was found to give a poor estimate of the induced moment. The results using the distributed LMO model were obtained in two ways. The (16) Jasien, P. G.; Stevens, W. J. J . Chem. Phys. 1986, 84, 3271.
8266
The Journal of Physical Chemistry, Vol. 93, No. 25, 198'9
calculation with complete self-consistency between polarizable points in a molecule and the total field from the other molecule (static plus polarization fields) gives a good prediction of the S C F result. The induced moment obtained when only the static field was allowed to polarize the points (noniterative model) was noticeably less accurate. This result would be expected with use of an accurate model of the molecular polarization because the supermolecular S C F procedure brings in terms corresponding to the terms in the classical self-consistent polarization scheme." The contribution of the polarization field to the induced moments is responsible for the potential singularities in the classical calculations of induced moments and adds a relatively significant extra expense to classical calculations on large systems. Therefore, proposals have been made to do away with this part of the classical polarization calculations for practical appIi~ations.'~J* This example is illustrative of the results obtained in a variety of other complexes having weak nonbonded interactions. Decompositions of the supermolecular S C F results show that short-ranged overlap effects such as charge transfer will, if neglected, generally lead to errors of the same magnitude as the observed deviations of the self-consistent LMO polarization model results. Another evaluation of the performance of the model in nonuniform fields has been made. In an analytic applied field, a distribution of dipole polarizable points behaves equivalently to an expansion in terms of certain higher order polarizabilities located at a single point. The computation of these higher order es from the model quantifies the expected deviations from ab initio results. The polarizability with the longest range of effect beyond the dipole-dipole term is the dipole-quadrupole A,,, which quantifies the quadrupole e,, induced p~larizability,'~ by an applied uniform field along coordinate direction r. The specific function of A used for comparison is a sum of the absolute values of the unique tensor components exemplified by (ASCF)
=
E
c
r=x.y,r sSf=x.y,r
14,5Fl/18
(7)
and indicated by enclosing a tensor quantity in brackets. Some SCF and model results computed by using the basis sets described in the next section are given in Table 11. These were calculated in coordinate systems diagonalizing the molecular (Y tensor. However, the relative percentages of the model contributions to the total S C F quantities are approximately independent of the coordinate systems used. These and other results indicate that an error of roughly the same absolute magnitude exists with the model calculations in molecules of differing sizes. The total magnitude of the effect of the quadrupolar response appears to increase rapidly with molecular size. For a system the size of a water molecule, the improvement in this aspect of the nondipolar polarization effect is only partially included in the model calculations; for larger molecules the relative accuracy of the representation is very much better. This seems to indicate that large molecules in nonuniform fields can be treated with roughly the same level of error as small molecules.
Transferability Tests A number of calculations of LMO polarizable points are presented here for several small molecules having "functional groups" in common. The molecules variously used are listed in Tables 111, V, and VI. Calculations on these can test the proposition that the centroids and polarizability tensors of a functional group can be transferred between otherwise different molecules. Calculations are also presented, differentiating between static and induced interactions as the source of deviations from perfect transferability. The molecular geometries are standardizations by functional groups of gas-phase structures,*Oinvolving only very small changes (17) Stone, A. J. Chem. Phys. Lert. 1989, 155, 102. (18) Niesar, U.; Corongiu, G.; Huang, M.-J.; Dupuis, M.; Clementi, E., to be published in I n t . J . Quantum Chem.. Symp. (19) Buckingham, A . D. Q.Rec. (London) 1959, 13, 183.
Garmer and Stevens for standardization in any particular molecule. This simplifies the comparisons, particularly for the centroid locations. It was observed that if this standardization was not done, the bond centroids swing along with variations in bond angles. Also, the location along the bond, given in terms of the percentage of the bond length, is well conserved. The conformation of ethanediol was taken from ab initio calculations predicting an internal hydrogen-bond-like arrangement of the hydroxyl groups.21 The calculations were performed with a standard triplecvalence basis set using the effective core potential (ECP) approximation?2 Diffuse polarization functions specified by Werner and MeyerIs were added to the non-hydrogen atoms to approximately maximize the molecular polarizability. Polarization basis functions were added to all hydrogen atoms, excepting those attached directly to the ring carbon atoms in phenol. The exponent of 0.321 1 au used for these functions matches one of the s-shell exponents and also approximately maximizes the molecular polarizability for methane. The polarizabilities obtained for these molecules with such basis sets are expected to be well within 10% of the C P H F limits, with the worst deviations expected for polarizations perpendicular to the plane in planar functional groups. The LMO centroids and polarizabilities of several of the smaller test molecules have also been checked with larger basis sets to verify that the trends observed in a series of calculations can be expected to be real. It is also worth noting that CPHF results have been found to be reliable predictors of the experimental results in large basis set calculations on closed-shell neutral systems.I5 The effect of correlating the wave function is usually a simple 5-10% increase in the magnitude of all components. The first functional group to be examined was the methyl group, consisting of a standardized tetrahedral set of connections with three C-H bond lengths of 2.06 bohrs. The bond length of the connection to the rest of the molecule was variable. In every case the molecular orbitals localized so that a centroid was positioned along each C-H connection and along the bond connecting with the rest of the molecule. For an examination of a particular C-H centroid location and polarizability, the structures and polarizabilities were transformed so that the methyl carbon was at the origin, the C-H bond in question was along the x axis, and the x y plane bisected the angle between the other two C-H bonds such that the hydrogens had negative y coordinates. The resulting C-H centroid properties for a few of the test molecules are listed in Table 111. The mean location and polarizability tensor of the centroid, weighted by the number of hydrogens of each listed type, were also computed. One observation from studying the variation among the polarizabilities is that the average polarizability, i%, is no better conserved than the individual components. By way of quantifying this, the expected standard deviation of the average polarizabilities, from the independent variability of each component, is 0.14 au. This is identical with the standard deviation actually observed for the average polarizabilities. The best centroid location is obviously along the C-H bond axis, since even when this is not mandated by symmetry, the distance of the centroid from the axis is only 0.03 bohr at most. The mean distance along the bond is 69 (1)% of the distance from the carbon to the hydrogen atom which also appears to be sufficiently well conserved for most applications. The mean polarizability tensor is almost cylindrical about the C-H axis, a condition that holds exactly only in methane. The small anisotropy between the ayy and azzcomponents was determined to be significant, however. The ayycomponent was larger in all of the 17 determinations for molecules other than methane. These results, and the off-diagonal polarizabilities, are consistent with 0-10' rotations of the approximate polarizability cylinder from the C-H bond axis in the ( 2 0 ) Harmony, M. D.; Laurie, V. W.; Kuczkowski, R. L.; Schwedeman, R. H.; Ramsay, D. A,; Lovas, F. J.; Lafferty, W. J.; Maki, A . G.J . Phys.
Chem. Ref. Data 1979, 8, 619. ( 2 1 ) Ha, T.-K.; Frei, H.; Meyer, R.; Giinthard, H. H. Theor. Chim. Acta 1974, 34, 211. ( 2 2 ) Stevens. W. J.; Basch, H.; Krauss, M. J . Chem. Phys. 1984,81. 6026.
Transferability of Molecular Distributed Polarizabilities
The Journal of Physical Chemistry, Vol. 93, No. 25, 1989 8267
TABLE 111: Some Methyl Group C-H Bond Centroid Positions and Polarizability Components in a Variety of Molecules”
centroid position molecule methane
ethane propene
acetonitrile methanol
methylamine propane
ethanol methylacetylene dimethyl ether
acetic acid meanb data std dev
X
Y 0.0
1.389 1.403 1.394 1.384 1.345 1.390 1.413 1.422 1.398 1.405 1.405 1.398 1.389 1.360 1.393 1.415 1.369 1.358
-0.007 -0.001 0.005 0.028 -0.002 -0.010 -0.019 -0.005 -0.013 -0.003 0.010 -0.003 0.017 -0.004 -0.004 0.004 0.025
1.387 0.022
0.003 0.013
polarizability tensor major components z
a
aYY
an
(YXY
aYx
0.0 0.0 0.0
-0.002
4.04 3.90 3.63 4.07 3.90 3.67 3.70 3.85 3.77 3.91 3.83 3.90 3.86 4.06 3.73 3.64 3.61 3.86
6.04 5.84 5.05 6.09 5.61 5.47 5.43 5.64 5.65 6.05 5.54 5.87 5.68 5.80 5.82 5.10 5.27 5.63
3.04 3.20 3.28 3.60 3.52 2.82 2.92 3.21 2.97 3.10 3.30 3.19 3.27 3.79 2.76 3.08 2.99 3.43
3.04 2.66 2.56 2.5 1 2.57 2.73 2.76 2.69 2.69 2.58 2.66 2.63 2.63 2.59 2.62 2.75 2.57 2.51
0.0 -0.36 -0.29 -0.72 -0.36 0.08 0.20 -0.1 1 -0.11 -0.51 -0.34 -0.38 -0.34 -0.49 -0.22 0.20 -0.29 -0.37
0.0 -0.57 -0.97 -1.09 -1.17 -0.1 1 -0.20 -0.56 -0.37 -0.68 -0.60 -0.63 -0.62 -1.45 -0.33 -0.36 -0.68 -0.85
-0.001 0.004
3.86 0.14
5.68 0.28
3.25 0.28
2.67 0.14
-0.25 0.25
-0.66 0.42
-0.002 0.0 0.0 -0.008
0.0 0.007 0.0 -0.002
0.0 0.006 0.0
0.0 -0.010
0.0
%x
“Group is assumed tetrahedral with 2.06-bohr C-H bond lengths. A standard placement of each bond was made as illustrated in Figure l a where atom A represents the carbon, atom B represents the hydrogen, and bond A-C is the connection to the rest of the molecule. In the molecules having more than one class of equivalent bonds, the bond in the plane of the heavy atoms is listed first. The unlisted components are not significantly different from zero in most cases. All quantities are in atomic units. Weighted by number of bonds. various molecules. The symmetry between r,s and s,r components, which would be observed for an isolated polarizable object, is only slightly conserved for the LMO’s. This is the largest and clearest effect of the molecular environment observed among the tabulated C-H bond LMO properties. Calculations were performed in which the polarization of the LMO’s of the rest of a molecule was decoupled from the polarization of the methyl C-H bonds. This was achieved by using the zero-field LMO’s as the initial guess for finite-field S C F calculations. The non-methyl LMO’s were frozen in the S C F cycles by zeroing their off-diagonal Fock matrix elements and in the final LMO calculations by allowing rotations only among the methyl LMO’s. The anisotropy of these decoupled methyl C-H polarizabilities was smaller and usually the reverse of that for the coupled polarizabilities, that is, usually a,, > ayy. This seems to indicate that interactions between the first-order LMO’s are responsible for the anisotropy of the coupled tensors, rather than interactions of the first-order LMO’s with the zero-order environment or the structure of the zero-order LMOs. The presence of a more polarizable group than another C-H bond is apparently effectively enhancing the polarizabilities in the direction of the group, but only by a relatively modest amount compared with what might be included from such an interaction in a calculation including explicit classical intramolecular induced moment interactions. As would be expected if this argument was correct, the overall transferability of these decoupled polarizabilities was somewhat better than for the fully coupled polarizabilities. In lieu of a lengthy discussion of particular cases, a summary of the other conclusions from these calculations follows. Although some trends can be observed and explained, the deviations from absolute transferability seem to be due to effective interactions of a complicated nature between the C-H LMO and the functional group attached to the methyl group. This is true as well for the variation between the polarizabilities when the C-H bond is in or out of the plane of the heavy atoms in the molecule, and the magnitude of the deviations is comparable. The differences between uncoupled and coupled polarizabilities are often not very large and also appear to be complicated and, thus, probably difficult to model. The final anisotropic mean polarizability tensor determined is worth testing for use in methyl groups in arbitrary molecules, and in CH,XY and CHXYZ groups as well. The data from the complete set of molecular LMO polarizability calculations are too lengthy to reproduce here. Thus, we present only a summary of the important results from consider-
ations of several other types of bonds. The centroids and tensors associated with bonds can be cataloged naturally according to the bonded atoms and their hybridization, with generally acceptable errors resulting. Given this, these should be designated as bond polarizabilities rather than as functional group polarizabilities. The first nine entries in Table IV are the mean values for those centroids and bond polarizability tensors which can be generated from calculations on the 11 molecules listed in the top section of Table VI. The average error of the centroid locations is approximately 0.02 bohr with no deviations larger than 0.05 bohr occurring. The average deviation of the polarizability components, not trivially zero because of local symmetry, is roughly 0.2 au. The calculated properties of LMO’s assigned to multiple bonds also catalog according to the bonded atoms, with deviations of roughly the same magnitude occurring in each assigned data set. The largest observed deviations were of the order of 0.6 au. These were associated with certain polarizability components of specific single bonds between atoms in the a and p positions relative to an unsaturated atom. The bulk of the components in these positions deviate from the mean values by less than this amount. It is difficult to catalog these deviations for predictive purposes because of the lack of an observable pattern in the limited data set. A consistently negative out-of-plane parallel polarizability component was found for C-0 bond L M O s in C(sp2)-GX(sp3) fragments. This implies movement of charge in a direction opposite to an applied field and can only result from an interaction between the perturbations on the LMO densities, an effect which would not be present in UCPHF calculations. Deformation densities induced by the field verify that the induced dipoles in the model reflect a real effect on the charge density in this case, not an artifact of the localization. The general molecular structure transformation scheme used to examine these centroids and polarizabilities is illustrated in Figure la. In the case of bond LMO’s, three atoms can be used to position each centroid. Although the choice of the atom C is sometimes arbitrary, various choices have little effect on the final properties. In some cases, there is local symmetry inherent in the deficition of the bond polarizability, constraining the values or relative values of some of the centroid positions and polarizabilities. An example of this is the C(sp3)-H bond which, for generality, must have equal a and a,, components, all zero off-diagonal polarizabilities a n r a centroid located on the bond axis. These averaged and symmetrized properties listed in Table IV were used
Garmer and Stevens
8268 The Journal of Physical Chemistry, Vol. 93, No. 25, 1989 TABLE IV: Symmetrized Mean LMO Polarizable Points by Type and FragmenP* centroid position polarizability tensor elements X Y 2 a ffxx ffYY ffzz ffXY ffX2 aY2 C-H bond LMO in fragment C(sp3)-H-X from 14 distinct samples' 0.0 3.81 5.64 2.90 2.90 0.0 0.0 0.0 1.399 0.0
ffY.
a,,
ffZY
0.0
0.0
0.0
1.385
-0.020
0.0
C-H bond LMO in fragment C(sp2)-H-X(sp2) from 6 distinct samples 3.31 5.40 2.65 1.88 0.11 0.0 0.0
-0.73
0.0
0.0
0.967
0.038
0.0
0 - H bond LMO in fragment 0-H-X from 6 distinct samplest 1.70 2.84 1.38 0.89 -0.14 0.0 0.0
-0.74
0.0
0.0
1.436
0.0
0.0
0.0
0.0
0.0
1.416
-0.040
0.0
0.0
0.0
1.027
0.044
0.0
C - 0 bond LMO in fragment 0-C(sp3)-X from 5 distinct samples' 2.31 4.40 1.33 1.20 0.18 0.0 0.0
-0.50
0.0
0.0
1.022
0.03 1
0.0
C - 0 bond LMO in fragment 0-C(spz)-X from 2 distinct samples 1.69 3.83 1.39 -0.13 -0.38 0.0 0.0
-0.36
0.0
0.0
1.266
0.0
f0.602
C=C bond LMO's in fragment C-C-X from 2 distinct samplest 7.19 11.16 3.61 6.80 0.0 0.0 0.0
0.0
0.0
0.0
1.463
0.0
f0.463
C=O bond LMO's in fragment C-0-X from 2 distinct samples' 2.69 4.70 1.43 1.94 0.0 70.46 0.0
0.0
0.0
0.0
0.172
-0.232
f0.516
0.19
f0.34
70.39
0.115
-0.244
0 lone pair LMO's in fragment O(sp3)-C(spz)-X(sp3) from 2 distinct samples f0.495 2.36 3.20 1.52 2.36 0.17 f0.56 70.32 -0.23
71.08
10.8 1
0.221
f0.553
0 lone pair LMO's in fragment O(spz)-C(sp2)-X from 2 distinct samples 1.67 2.07 1.69 1.24 0.74 0.0 0.0
0.0
0.0
C-C bond LMO in fragment C(sp3)-C(sp3)-X from 4 distinct samples' 3.67 6.03 2.49 2.49 0.0 0.0 0.0
C-C bond LMO in fragment C(spz)-C(sp3)-X(sp2) from 2 distinct samples 3.20 5.71 2.31 1.56 0.42 0.0 0.0 -0.32
0 lone pair LMO's in fragment O(sp3)-X(sp3)-Y(sp3) from 7 distinct samples
0.0
2.05
2.28
2.19
1.69
0.17
f0.25
70.28
0.34
"See Figure I for the definitions of the positioning; atoms A, B, and C in the fragments are listed in order. Both LMO's of each double bond or lone pair, related by reflection, are included by giving the appropriate signs for the centroids and components with f and r symbols. The averages were taken over the entire set of LMOs of each molecule, but the number of samples indicates how many sets of equivalent LMOs there were among these. *Dagger indicates a case where LMO properties had to be modified slightly as described in the text to satisfy general local symmetry requirements.
TABLE V: Some sp3 Oxygen Lone Pair Centroid Positions and Polarizability Components in a Variety of Moleculesa polarizability tensor elements centroid position molecule X Y 2 a ffrr a"" a,, water 0.165 -0.224 0.505 2.10 2.16 2.13 2.02 methanol 0.172 -0.23 1 0.515 2.04 2.71 1.71 1.70 ethanol 0.170 -0.235 0.516 2.07 2.72 1.76 1.72 dimethyl ether 0.170 -0.230 0.519 2.09 2.55 2.40 1.31 ethanediol 2.06 1.63 2.60 1.94 0.534 0.166 -0.2 17 acceptor, facing 1.64 1.92 0.500 2.00 2.43 -0.245 acceptor, opposing 0.187 2.73 1.75 0.502 2.07 1.74 0.187 -0.246 donor, facing 2.02 2.66 1.55 1.86 0.522 0.i54 -0.219 donor, opposing 2.34 3.18 1.51 2.32 -0.245 0.495 acetic acid 0.1 17 2.38 3.23 1.53 2.39 -0.244 0.494 formic acid 0.113 1.44 1.78 0.502 2.42 4.05 0.138 -0.235 phenol 'The OXY fragments were positioned as shown in Figure Ib, and a bond angle of 107' was assumed. The various C-0 bond lengths were taken as 2.703 bohrs and the 0 - H bond lengths as 1.814 bohrs. Only the lone pair with a positive z coordinate for its centroid is given for all molecules but ethanediol. I n this case, properties for both lone pair LMO's at each oxygen atom are given and labeled according to whether they are facing or opposing the other hydroxyl group, relative to the oxygen atom. The off-diagonal tensor components are smaller in magnitude than the diagonal elements.
in the comparisons with other molecules, not included in the original list, and for calculations of the molecular polarizabilities listed in Table VI. The symmetrization for C(sp3)-H has the largest effect on any particular bond polarizability component and the most deleterious effect on the estimated molecular polarizabilities, overall. The other type of LMO density encountered is very much like the classical idea of an oxygen atomic lone pair. A transformation was performed on these centroids and polarizabilities to put them into a common axis frame relative to the X-0-Y or UV(W=O) group as illustrated in Figure 1b,c. A number of the X-0-Y oxygen lone pairs are listed in Table V. Examining the centroid locations and polarizabilities showed that the lone pairs occurred
in pairs with properties related by a reflection in the XOY plane or the plane bisecting the UWV angle, respectively, for an sp3 or sp2 hybridized oxygen atom. This result was little affected by the rest of the molecular environment, even for the ethanediol structure, wherein the donor hydrogen of the internal hydrogen bond is located much closer to one lone pair centroid than the other in the acceptor group. Besides the expected dependence on the oxygen atom bonding pattern, single or double, the lone pair properties were found to be sensitive to the hybridization of the X or Y atoms for sp3 oxygen. The oxygen lone pair centroids are rotated out of the plane bisecting the XOY angle toward any unsaturated atom in t h e CY position. The components and average polarizability were
The Journal of Physical Chemistry, Vol. 93, No. 25, 1989 8269
Transferability of Molecular Distributed Polarizabilities
TABLE VI: SCF Molecular Polarizabilities from Finite-Field Calculations, Followed by Values Estimated with the Use of the Mean LMO Polarizabilities of Table IV' molecule methane ethane propane ethylene propene methanol ethanol ethanediol dimethyl ether formic acid acetic acid I-propanol 2-propanol formaldehyde acetaldehyde carbonic acid av % dev
a
axx
aYY
0111
16.1 15.3 27.3 26.6 38.5 37.8 27.6 27.6 39.5 39.0 19.4 19.6 30.9 30.9 33.6 35.2 31.0 31.6 19.8 20.1 31.2 31.5
16.1 15.3 29.3 25.3 41.8 36.6 36.9 36.9 51.2 46.9 21.7 21.3 34.5 32.0 34.9 37.5 36.3 33.0 24.2 23.8 35.4 32.5
16.1 15.3 26.4 27.2 37.7 37.8 24.3 24.8 35.9 36.8 18.7 19.6 30.0 30.8 33.9 33.4 28.5 31.4 21.1 22.9 33.8 36.1
16.1 15.3 26.3 27.2 36.1 39.1 21.5 21.1 31.3 33.2 17.8 17.8 28.1 29.8 32.1 34.6 28.1 30.5 14.2 13.7 24.5 25.8
42.0 42.1 41.7 42.1 15.9 15.3 27.6 26.7 24.0 25.0
48.0 42.8 43.3 39.1 20.7 20.8 33.0 30.7 29.2 29.6
40.3 41.9 42.2 44.2 15.4 15.1 27.5 27.1 26.1 28.0
37.8 41.7 39.5 43.1 11.7 10.1 22.4 22.2 16.6 17.3
2.1
5.2
6.5
5.1
axY
ffXI
aYz
(YY.
%X
(YZY
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
3.9
0.0
0.0
3.8
0.0
0.0
-0.1
-1.0
0.1
0.0
0.1
-0.7
-0.1
0.0
-0.9 -0.8 0.0
0.1 -0.1
0.0
-2.2
0.1
-0.1
-2.2
-0.1
0.0
-0.1
0.0
0.0
-0.1
0.0
0.0
0.0
1.2
0.0
0.0
1.o
0.0
0.0
0.3
0.0
0.0
0.5
-0.1
-0.1
0.6
-0.1
0.0
0.2
-0.1
-1.7
0.2
-0.1
-1.9
0.0
0.0
0.0
0.0
0.0
0.0
2.7
0.0
0.0
2.6
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
'The coordinate systems diagonalize the SCF tensor for each molecule. The molecules listed in the top section were used in developing the averaged LMO polarizabilities, and those in the next section were not. The statistics listed at bottom were tabulated from the entire list. ,.Y
Y
a ) Bond AB LMO's.
n
b ) A(sp3) lone p a i r s Y
c
\
\
B-A-
i
- - -->X
ZY
c) A(sp2) lone p a i r s .
Fipre 1. Definition of the standard axis systems used to develop mean LMO polarizabilities. The molecules are translated and rotated so that each three-atom fragment has its first atom in the position indicated by A, its second atom along the direction defined by vector B, and then its third atom in the plane ABC with a positive y coordinate.
consisently and significantly changed when compared with X(sp3)-0-Y(sp3) oxygen lone pairs. This naturally leads to the following three definitions of fragments containing transferable oxygen lone pairs: O(sp2), X(sp3)-O-Y(sp3), and X(sp2)-O-Y(sp3). The hydrogen atom is considered to be an sp3 hybridized species here. The mean centroid locations and polarizabilities for these types are listed in Table IV. Additional dependence on the type of X or Y was not explored in detail. Note, however, that the phenolic oxygen lone pair is somewhat different from the carboxylic lone pairs and the water lone pair is also unique. The oxygen lone pair properties in C=O bonds are transferable between carboxylic acids, aldehydes, and carbonic acid. A test of the overall effect of using the mean LMO polarizabilities is t o compare the predicted molecular polarizabilities with
the calculated SCF values. This was done by using the symmetrized mean LMO polarizabilities listed in Table IV, and the results are given in Table VI. The molecular polarizability components are sums of the values from the mean L M O tensors which have been transformed to the proper coordinate system where the three atoms defining the fragment coincide with the proper atoms in the molecule. The percentage deviation of the predicted average polarizabilities is approximately 2%on average, and that of the components is 5 6 % . The accuracy of the prediction for molecules actually used to derive the LMO polarizabilities appears to be slightly better than for the others, an expected result considering the procedure used. The predicted tensors are very close to being symmetric, more'so than would be expected considering the accuracy of the prediction of individual components of a. The factors mentioned above in connection with the small deviations from transferability of the L M O properties are also responsible for most of the errors in predicting molecular polarizabilities. The majority of the mean LMO polarizabilities have ratios of principal components of the order of 2-3: 1: I . This effect is responsible for most of the anisotropy of the molecular polarizabilities. The substantially nonplanar molecules have bonds and lone pairs oriented in all directions, thus tending to cancel out some of this anisotropy from the overall molecular polarizability. However, when the mean LMO polarizabilities are used, this local anisotropy is reproduced correctly as is important for intermolecular polarization calculations. Conclusions The type of transferability, and the fair degree of results obtained for consequent molecular properties, demonstrated here with several types of bonds and lone pairs very probably holds for most other types. The prediction of the molecular polarizabilities is similar in quality to the best semiempirical methods of
8270 The Journal of Physical Chemistry, Vol, 93, No. 25, 1989
predicting molecular polarizabilities. The level of accuracy holds fairly well even for ethanediol, having an internal hydrogen bond in the conformation considered. Extended T systems hade to be treated as independent functional units in this model. First of all, the localization for these systems is very often unique and not physically understandable. A simple example of the effect of conjugation of double bonds on polarizabilities occurs in trans-butadiene. The parallel polarizability of each double-bond L M O is increased by 50% over that for the similar LMO in ethylene, in a direction approximately along the axis of the double bond. The other components and the C(sp2)-H bond LMO properties vary by approximately 5% from the mean properties given here, while the C(sp2)-C(sp2)bond polarizability must be considered as a new entry. The centroid positions for the double-bond LMO’s vary by only 0.02 bohr from ethylene values. The LMO density type not considered at all here is the core density of atoms, associated with n - 1 principal quantum number orbitals, where n is the quantum number of the valence orbitals. The polarizability of these could be ignored in most calculations because of very small expected interaction energies. Otherwise, it could be treated as a transferable spherical polarizable point, depending only on the atom and its valence state. This is an approximation commonly employed when putting in core-valence correlation effects, via core polarization, in ECP calculations. I t has been illustrated that the LMO polarizability model can reproduce some of the features of the electrical response of molecules in nonuniform applied fields at an appropriate level of accuracy. The most obvious method of extending the model given here to assemblies of larger molecules such as proteins is to combine the polarizabilities of the functional groups in the same way as was used to get the results in Table VI. Induced dipole moment interactions could then be entirely suppressed or be allowed between centroids in separate functional groups if they are free to rotate in such a manner that the intervening distance may
Additions and Corrections change drastically. These two proposals respectively correspond to a model including only the lowest order three-body interactions in the classical electrical model and to a rational implementation of the LMO polarizabilities in a completely self-consistent electrical model. Results given here and elsewhere3seem to show that either procedure could result in a reasonably faithful representation of molecular polarization in large systems. In order to truly test and quantify these concepts, it will be necessary to carry out benchmark calculations with nonuniform applied fields, to compare with models in systems such as ethanediol. Having reliable ab initio based methods of calculating polarization effects in large molecules is the most important goal of this work. Such methods would be useful in completely classical calculations and for quantum mechanical calculations with a polarizable reaction field.3 We would also like to point out that a polarizability model could possibly be a valuable extension to the electrostatic interaction now commonly used in docking biomolecules and for indicating possible reactive sites. Francl has shown that accounting for electrical polarization from an approaching nucleophile can improve the qualitative prediction of relative reactivity in small molecules.23 The proposal is to replace the perturbation calculations used in that work with a classical distributed model calculated from a much cheaper series of finite-field S C F calculations. In large systems the transferability outlined in this work could alternatively be used. Acknowledgment. The authors thank Dr. Paul G. Jasien of Biosym Technologies for discussion and help in preparing for this work and Dr. Morris Krauss for helpful discussions and comments on the manuscript. We also thank the referees for pointing out the need to clearly show the capabilities of the proposed methods.
(23) Francl, M . M. J . Phys. Chem. 1985,89, 428.
ADDITIONS AND CORRECTIONS 1989, Volume 93 A. V. McCormick, A. T. Bell,* and C. J. Radke: Influence of Alkali-Metal Cations on Silicon Exchange and Silicon-29 Spin Relaxation in Alkaline Silicate Solutions.
Page 1740. The authors have discovered an error in the labeling of Figure 7 and in the text discussing this figure. A corrected copy of Figure 7 is shown below. The sentence beginning the second paragraph of the right-hand column should read as follows: Proton or cation exchange is expected to occur much more slowly than wL,I7 and hence relaxation via scalar coupling would be expected to contribute to 1/ T2as shown by point C and to 1/ T , as shown by point A.
I Figure 7.
I In r
