Chemical Physics 59 (1981) 341-350 North-Holland Publishing Company IMOLECULAR POLARCIZABILITLES CALCULATED WITH A MOIXEIEID DIPOLE INTERACTiON B.T. THOLE TheoretfcnlChemists Group. The University of Groningen. 9747 AG Groningen, The Netlrerlands Received 16 February 1981 Tire point dipole interaction model for moleculu poltibihty recentiy proposed by Applequist, Carl, and Fung is modified by replacing the point dipole interaction by an interaction between smeared out dipoles. Rules are developed to indiczte plausible forms fer this modified interaction. The polarizabihties of a wide range of chemically different molecules can be calculated, using for each atom one polarizability independent of its chemical enviromnent. The errors are comparable to experimental uncertainty. Special care is taken to produce a model that tends to avoid infinite polarizabilities without use of cutoffs at short distances. 1. Mroduction The accurate calculation of the effects of polarizability on electrostatic interaction is an interesting problem in various fields. Traditionally the polarizability of molecules was calculated to predict optical rotation [l] and London dispersion forces [2]. More recently, inclusion of polarizability in molecular dynamics calculations has been studied [3]. Finally in many cases the stabilization of e.g. a charged or a dipolar moIecule surrounded by a polarizable solvent is important. It has been shown that this effect cau be included into quantum mechanical calculations on such a molecule, where the solvent is apprcximated by a set of interacting polarizabilities [4]. In all these cases one tries to avoid simulating the polarizability by introduction of a dielectric, because this macroscopic concept cannot be used but with the greatest care on the molecular level. The interactive dipole model of Applequist et al. [S], based on the early work of Silberstein [6] seems to be a good starting point to treat the polarizability of arbitrary configurations of polarizable atoms. This model is surprisingly successful in predicting mean poiarixabilities and 0301-0104/S1/0000-0000/%02.50 @ North-Holland Publishing Company is easily extended to the other fields mentioned, either by a matrix inversion method or by iterative solution of the equations for the interaction of the induced dipoles. Applequist et al. [S] used isotropic atom polarizabilities. They observed that their model predicts mean polarizabilities very well but that the predicted anisotropy is too large. Birge [7] showed that the large auisotropy can be removed by inclusion of the efIect of electron repulsion, which essentially decreases the polarizability component in the direction of the bond and gives errors in the anisotropy comparable to the experimental uncertainty. Birge needs 14 parameters to describe the atomic polarizabilities in different chemical environments and their tendencies to resist excessive polarizability along the bonds. The transferability of these parameters to other classes of molecules is rather bad [S]. This property is shared by all previous parameterizations. It is therefore interesting to discover trends in these parameters. This will give more insight into the physical process of polarization and may give new parameters that are more readily transferable. This paper shows that a suitably chosen modification of the interaction 312 B.T. T/de / Molecuiar polan;a&i&ies decreases the errors in predicted polarizabilities, while the number of parameters is strongly reduced. In the next section the theory of the point dipole interaction is repeated, and extended for our purposes. The third section then gives a procedure for modifying the interaction. 2. The point dipole model The modified dipole interaction model developed in this paper is based on the early work of Silberstein 163 and the more recent application of Silberstein’s model by Applequist et al. [5]. The nomenclature of rhe latter authors will be used and the reader is referred to refs. [S, 8] for a detailed discussion of the point dipole interaction model. The difference between the procedures presented here and those outlined by Applequist et al. is the modification of the dipole field tensor. The atomic polarizability tensor remains isotropic. This is in contrast with Birge’s model E7), which modies the atomic polarizability tensor and keeps the field tensor unchanged_ The molecule is considered as an arrangement of N atoms each of which has a polarizability. The induced dipole moment at atom p i.e. pp. can be calculated as a function of the applied electric field EP at atom p: where ocPis the atomic poiarizability tensor of atom p and T, is the dipoIe field tensor x2 Tpq= (t~)a-3(r~) yx i . xy xi? Y2 1 YZI7 (2) m zy z’_I where i is the unit tensor, r, is the distance between atoms p and q, and x, y and z are the cartesian components of the vector connecting atoms p and g_Eq. (1) can be rearranged into the single matrix equation A&=& (3) where a is a 3N x3N matrix containing the inverse of the atom polarizability tensors along the 3 x 3 “diagonals” which are ccupled by the dipole field tensors. Inversion of A yields B, a symmetric matrix called relay matrix, which plays the role of a super polarizabity; p,=BE (4) t?j=~-‘=(&-‘+~)-‘. (5) The molecular polarizability is the response to a uniform applied field @ and is thus obtained by contraction of B to a 3 X 3 tensor OL,,~: L,=[++=~,&- (6) Diagonalization of or,01 yields the three components of the polarizabiiity. In more general cases the whole matrix e is needed [4], and in moIecular dynamics eq. (3) may not be solved by matrix inversion but by an iterative pro- cedure. A well known property of the method presented so far, is that it may lead to infinite moIecular polarizabilities. For a diatomic moIecule the values of the polarizabilities parallel and perpendicular to the bond axis can be derived from the above: CYI/= (QI*+cLui-4cK&3/r3)/(l -4cYNYBIr6), (7) QJ_= (cY*+cYr3--2a*~n/r’)/(l -~,aQ/P). (8) When r approaches (~u~Gc&“~; “11goes to infinity. This is caused by the cooperative interaction (head to tail) between two induced dipoles in the direction of the bond. The generd criterion for a physically meaningful solution is that the matrix A, and therefore also its inverse 6, be positive definite. A matrix a is positive definite if all its eigenvalues are positive or, equivalently, if JAp>O for all vectors CL. That A and 6 must be positive delirite can be seen in two ways. Fit, when b has? negative eigenvalue -b then, when we t&e B equal to the corresponding eigenvector we obtain $.=~~=-_b& (9) which means that & and 2 are in opposite directions. Here we will only corsider frequen- E3.T. T?de / hfdecular p&&abilities 343 ties of the electric iie!d that are so low that 6 can follow the variations in &_This means that 6 and 8 may not be in opposite directions. That such a system is metastabIe is seen from the reaction of the total energy to a small perturbation SC in the dipoles. The total energy U is the energy of interaction between the dipoles, plus the energy needed to polarize the atoms. The latter energy is equal to $XP&cYQLP’~~and therefore: = _fi: .& +&‘& (10) After a perturbation S@ we have *u = _sg’ . ~+~fi*j$+$~fiLf~~fi+. (11) In the equilibrium state, where (3) is satisfied, the first order contribution to SU vanishes. For the equiIibrium to be stable, the second order change in U must be positive for ali 6fi, so S’U = S;‘liSE > 0 for all Sfi, 02) which means that A must be positive definite. When A is not positive definite then d2U < 0 for some SE and the system is metastable. A small perturbation will make all dipoles grow to i&&y. Further insight is gained by the introduction of the very useful tensor t, defined by: (T,), = s,r+ - 3XiXiP = (Cry*)-“‘(~;jLi -3 - 3 ftiiljll -5) = (apc&pZtij(U), (13) where u = .~/(a,cy,)~‘~ and 6, is the Kronecker de1ta. T,(x) = (~~s!-1’zt{x/(~~p,)“6~, (14) t is a shape fitnction which does nor depend on p and CJ.In this paper only isotropic polarizabiities will be used. Then 6 is diagonal and we have f = &-Uz~&-I/z, (151 W=~-l+f,~-‘lz(~tI)~-1/2~ (16) Now $&, = (~-1/2Cr)‘(~+~)(~-l/2~)_ (17) When it-i is positive definite then the rhs of eq. (17) is positive for all G and therefore A is positive definite. Thus, 4 is positive definite when ail eigenvalues oft are larger than -1. 3. The modified dipole intenarfioa When in the point dipole interaction atoms are coming closer than a certain limit, the polarizability becomes i&mite. In reality this is never observed. There must be some process which damps the interaction. This process is probably already at work between bonded atoms because Applequist needs rather small atomic polarizabilities to describe molecules and he observes that his predicted anisotropies are too large. The idea developed in this section is that the damping effect may be simulated by changing the interaction :ensor T such that it does not behave as re3 for small r. A set of more or less tentative principles will be given which were used as guidelines to investigate which changes to T might be applied for obtaining a good fit to experimental data, while avoiding unphysical forms. The principles were chosen so as to supply possible starting points for future theoretical analysis. The following list is given in the order of roughly incre&ng tenta- tiveness: (a) T must remain a tensor to ensure that rotation of the coordinate system does not change the result of the calculation, (b) The most essential idea of the model is that the interaction between any two polarizabilities is of the same form. In a classical macroscopic system, consisting of any set of conducting bodies or bodies of some dielectric constant, when all the dimensions of the system, i.e. the dimensions of the bodies and the distances between them are expanded by some factor, then alI the polarizabiities increase by the third power of that factor. In classical electrostatics this transformation is equivalent to a change of the unit of length. Simple extrapolation of the scaliig property to the atomic world is strictly speaking not allowed because atoms have internal structure and are not 344 B.T.TkdeJ MolecularpoIarizabil!ties homogeneous as macroscopic bodies are. Nevertheless, we may hope that the macroscopic scaling principle is still approximately valid. When we write T in the form (14) with arbitrary t and if in eqs. (5) and (14) we multiply all x by a factor of k and al! polarizabilities by k3 then & changes to k3& and ? to Ic-% and so 8, which denotes the total polarizability, changes to 61 ={(,p&‘+-~)}--1= pg. (18) This shows that in the point dipole case the scaling property exists and that it is preserved when (14) is satisfied. We may consider t as the interaction between two atoms with unit polarizability. In this model it will be assumed that t(u) remains independent of the atoms involved so that we need only study the interaction between unit polarlz2bilities to know all the interactions. Note that (cx~,)“~ has the dimension of length and that it may be associated with the average radius of atoms p and q. so that we may say that ‘; is only a function of an effective distance of interaction u. The quantum mechanical detail of atoms p and q apart from their ‘-radius“, is thus neglected. The use of (c+Q)*‘~ as the scaling distance is arbitrary. Forms lie (a;” tar)/2 and {((up+ %.)/P are also possible, but they do not allow the elegant forms of eqs. (15) through (17). (c) r&) is constrained to the tensor form q(u) = -a’+9(2l)/arsazl, (1% where c_”is some spherically symmetric posential. This assumption is made to retain the interpretation of t ^osa dipole-dipole interaction where a dipole is thought to be built up from hvo infinitesimally shifted monopoles, which in our case will not be point charges. Then 4, is the corresponding monopole-monopole interaction. 190rk-g th2t d/du&u) = L(.$(u)/u and &+/au; = 8, it is easily seen that (20) (21) (di In the pomt dipole case cp= u-l SO v'= --u I q" = 2~~’ and tij = sijll -3 - 3 u~rijlc -5. (22) These relations describe the limiting behaviour for large U, when the atoms do not touch. (e) For a homonuclear diatomic molecule, solution of (l)-(6) for two atoms of unit polarizability shows 2 2 ali=- and cyI=p l-cp’fli’ We assume that ail= oL when II = 0, making the compound atom isotropic. Then cp”= cp’/u for small U. Solution of this differential equation shows that in this region we have to first order cy’= CU,with c an arbitrary constant. This behaviour can also be obtained from the assumption that p is the potential of some “well behaved” charge distribution p. We start with the well known relation si(u) = 3 I 4iipu2 drt, 0 saying that the field at u for a spherically symmetric charge density is equal to that of the charge contained in the sphere with radius u, concentrated at the centre. Because of (20) u 40’= -1 4q.m’ drc. (25) II When p is finite when u approaches zero we have to first order in t&s region: U 9’= -1. U?-I 47i_p(0)rc2 du =-$srp(O)rr (26) 0 and (2”= -$rp (0) (27) and we see that c =--gap(O), a negative quantity when p is everywhere positive. Furthermore for large u, when p decreases fast enough U al I 4?j-p2 dn - I 47rpu’dic = Q, Gw 0 0 the total charge. So in this region +Y’=-Q/u’ which shows that Q = 1. It would be more elegant to consider ‘p as the interaction between two charge distributions, instead of that between a point charge and a charge distribution; a(u) = II du: du;p(ul) p(u~+c)/u~=- (29) This may be necessary for further theoretical investigation. For the present, however, we are satisfied with a procedure which gives us smooth curves which roughly follow physical intuition.We conclude that t is derivable from some density p, which must be finite at the origin, fall off rapidly for large r, and must contain a total charge of unity. Thus the unit point charge giving the point dipole equations can be considered to be “smeared out”. (f) In the previous section it was shown that it is the descaled tensor $ which governs the stability of a polarizable system. The requirement that ali its eigenvalues be greater than -1 is rather complicated. However, from eq. (23) we see that the polarizability of a diatomic molecule remains finite when y” and cp’/rr are less than unity, independent of the values of the CX’Sand u’s In polyatomic molecules favourable interactions between pairs of atoms may push the polarisability to infinity even though each pair is separately stable. Therefore unity is an upper limit. However the interaction decreases rapidly with distance and geometrical restrictions do not allow favourable interaction between every possible number of atoms. Therefore the limit of the instability region will probably not be impractically low. Numerical experience indicates that, with the forms studied in this pa2er, infinite polarizabilities become improbable when (pe and cp’/u remain below 0.6. Concluding we assume that infinite polarizabilities will become highly improbable if not impossible when c9’ and cp’/u are below 0.5 for every value of U, independent of the values of the atomic polarizabilities. 4. Calculations .. Some forty more or less different functions to fit the polarizabiities of the sixteen molecules considered by l3irge [7] were examined. In each case the atomic polarizabilities and the parameters in the interaction tensor were varied simultaneously to minimize the root mean square (T of the relative errors in the 48 components. Directly in the beginning, it became clear that per atom only one isotropic polarizability independent of its valence state is needed. Introduction of different LY’Sfor one atom did not significantly improve the fits. Thus it seems that the apparent additive polar&ability of an atom is strongly correlated to the geometry. For example a carbonyl 0 is certainly different from a hydroxyl 0, but this difference is sufficiently displayed in the geometry of its surroundings. The scaling principle describes this effect very weil. Further it appeared that one must be careful not to use a too flexible fit function because this flexibility often leads only to wildly fluctuating functions giving artificially good results which are physica!ly meaningless. This is why the optimum form of T could not be obtained by simply using a general function with many param eters. Only smooth shapes of cp’/zl and q” could be accepted. After the restriction to a tensor t derivable from some well behaved density p, functions of the types presented in table 1 were tried. This table gives an impression of the variation in the results when p is varied. Of the one parameter forms pr ... pa. p4 is the simplest and, surprisingly, the best. This function is a simple conical charge density with 1.662 .& as the radius of its base (for an atom of unit polarizability). ps and & have been tested to see whether there was a preference for a more curved shape of p. This preference is rather small as is seen from the small decrease of c going from p4 to ~6, and from the small coefficient with which ps is mixed in. Case 7 was tried in order to study the behaviour when principle (c) was relaxed i.e. when e was not constrained to be the second derivative of a potential. For this purpose the factors cp’/rc- r$’ and q’/u in (21) were chosen to be of the form p-: but they were allowed to have different radii a and b. Not much use was made of this freedom to decrease o or to make a and b different. Therefore, it may be concluded that the factors 346 B.T. Thole /Molecular polarizabilities Table 1 Charge densities and the corresponding values of the fitting error and the parameters pl= a3/8?; exp (-au) &u =-4u-3{2-(a’u”t2au+2) exp (-au)} pz= 3aj4s exp (-I&) pi/u = -ue3{1 -exp (-aus)) p,=3u/4n{l-(l+l/auS) exp (-llnu3)) (2;/u=-a{l-exp(-l/au3)} PJ=3/3i+z-uu)/aJ uca =o uza L/u=-_(4a-3u)/a4 u 1.662 ++’is twice as large as PI/U, the mean polarizability is approximately conserved. When the atoms penetrate, the interaction is effectively damped with respect to the point dipole case. The anisotropy is also damped because p” and p’/u approach each other. The interpretation may be that when two atoms touch, a new body is formed with an elongated shape giving an anisotropic polarizability, as in classical electrostatics. When the atoms penetrate further the shape becomes less elongated until at u = 0 a composite spherical atom is formed. Of course the model does not pretend to predict the polarizabilities of atoms as built up from smaller atoms. For 14~1 the individuality of the constituent atoms disappears and probably only quantum mechanics can describe the _-dual transition to the polarizability of the compound atom. This region is unimportant for most applications but it is gratifying that the polarizability is damped in this region, and therefore no problems are to be expected when, accidentally, atoms come into close contact in e.g. MD calculations. 5. Discussion The present model, using only five parameters gives an improved fit with respect to Birge using 14, and Applequist using eight parameters. Part of the improvement is due to the complete optimization used, where Applequist and Birge used partial optimizations. The error in the componenB is 6% and in the mean polarizabilities 3.5% which are about the experiaental errors. Tne parameter values of the model are transferable to different classes of molecuIes in contrast to those of Applequist and Birge. This fact, together with the small number of parameters indicates that the physical process behind the interaction of polarizabilities is reasonably described. The ideal situation in this kind of semiciassical model would be to predict experimental molecular polarizabilities from experimental atomic polarizabilities. The point dipole interaction is strongly amplifying and therefore atomic polarizabities are needed which are much smaller than the additive values, especially for H with its short bonding distances. Conversely the smeared-out dipole interaction is on the average mildly damping and therefore i? is accompanied with atomic polarizabilities which are about 20% higher than the additive ones. Hartree-Fock poiarizabilities are considered to be 20 to 60% higher than experimental values. But especially the hydrogen value, which is exact, shows that the smeared-out dipole model approaches the ideal more closely than the point dipole model (table 3). Table 3 Polarizabilities of atoms in polyatomic molecules at 5893 8, Atom H(alkane) H(alcoho1) HMdehyde) H(amide) C(alkane) C(carbonyI) C(nitrile) N(amide) N(nitrilej O(&oholj O(ether) O(carbonyl) ” Ref. [SJ. Polarizability (A3) additive Applequist Hartreemodel”’ et al. ‘) Fock” 0.407 0.135 0.667” 0.405 0.135 0.167 0.161 1.027 0.878 1.74 1.027 0.616 0.928 0.36 0.530 1.00 1.236 C.520 0.604 0.465 0.73 0.651 0.455 0.841 0.434 this work 0.514 1.405 1.105 0.862 ‘) Ref. [14] from variation-perturbation theory. <) Exact. It is interesting to note that Applequist’s atomic values for C, N and 0 decrease in chemical environments with shorter bond distances. In the present model this trend is produced automatically by damping the interaction more strongly with shorter distances. The present model is connected to some notions recently suggested by other authors. Oxtoby and Gelbart [9] proposed to describe the interaction of the induced dipoles of penetrating atoms as the interaction between penetrating dipole densities. This is equivalent to the principle that the dipole field tensor must be the second derivative of the interaction between smeared out char’ge distributions. This notion has been criticized [lo, 111 but its success suggests that it has some value. Maybe any approximation which uses the point dipole interaction for large distances until the atoms “touch” and then damps the interaction in a smooth way, making it zero for r = 0, will do reasonably. The freedom to choose the density function and the generally decent behaviour of the polarizabity are further responsible for the good results. Thus the value of this kind of approximation may be largely empirical. Bounds and Hinchcliffe [12] and Winicur [13] noted that the anisotropies of diatomic molecules are very similar as functions of the internuclear distance and that they seem to be determined by the “shape” of the molecule, defined by a suitably chosen measure of its length-to-width ratio. This is strongly connected to the scaling principle of the present model which uses rW(crsq)-“6 as a shape parameter. Inclusion of the effects of smearing out and the scaling principle into the interaction tensor is thus a generalization of existing notions of the apparent ‘cclassical” behaviour of the molecular polarizability. The small number of parameters, their transferability to different systems, and the consistent physical picture give confidence to the calculations in diEcult cases such as formation of hydrogen bonds, strong deformation of bonds and close contacts between different molecules or different parts of luge molecules. Because each atom has only one polarizability, ambiguity in the assi,onment of the valence state to the atom has no consequences. The model is simple enough to be applied to large systems. No exclusion of bonded interactions or arbitrary cutoffs at short ccntacts occurring otherwise are needed. Hence the present method is very well s*Gted to be incorporated into empirical ener,y expressions, as they are used, e.g. in molecular dynamics calculations, energy optimization of crystal and protein structures, etc. These computational schemes might be improved considerably in this way. Acknowledgement The author is indebted to H.J.C. Berendsen and P.Th. van Duijnen for their valuable suggestions and the help in preparing the manu- script. Appendix: Details of the computation Probably any computer program using Applequist’s point dipole interaction model in spectroscopy, molecular dynamics, Monte Carlo or energy minimization, can be easily adapted to the smeared-out dipoie interaction. When r is the distance vector between two atoms with polarizability ccr and a~ and s = 1.662(ara~)“” then when r>s we have, exactly as in Applequist’s -model & = &Jr3 - 3*,ri/r’. (A-I) When r < s the interaction has to be changed into r, = (Ic3-3v~)6i;/r3-33t;4(rir;/r5), where u = r/s. (AZ) A FORTRAN program, running on a CDC Cyber 170/670 computer, under NOS, and which calculates and fits molecular polariiabilities, can be obtained from the author. In MD and MC calculations, polarizability is, for practical reasons, often treated by changing the magnitude or the position of the point charges in the molecule as a function of the electric field. Artificial devices are then usually needed to prevent excessive polarization. It may be possible to mimic the smeared-out dipole model also in these calculations. An obvious guess would be to use a mcdified electric field, viz. the field of a smeared-out charge distri- bution: E:=rJr' [r>s]. (A.3) Ei=(4u3-3u4)ri/r3 [rcs]. For the calculation of the total energy the corresponding potential is given by p = l/r, [r>s]. (A.3 dp =(u'-2u3+2u)/r [r