Chemical Physics Letters 394 (2004) 37–44 www.elsevier.com/locate/cplett
Atom–bond pairwise additive representation for intermolecular potential energy surfaces F. Pirani
a,b
, M. Albertı´
a,1
, A. Castro a, M. Moix Teixidor a, D. Cappelletti
b,c,*
a
c
Dipartimento di Chimica, Universita` di Perugia, 06123 Perugia, Italy b INFM, Sezione di Perugia, Italy Dipartimento di Ingegneria Civile ed Ambientale, Universita` di Perugia, via G. Duranti 93, 06125 Perugia, Italy Received 14 June 2004 Available online 17 July 2004
Abstract A method has been developed to describe the force field of atomic species interacting with hydrocarbon molecules, either aliphatic or aromatic, of use for molecular dynamics simulations. The potential energy surfaces are represented by a simple analytical form written as a sum of atom–bond interaction contributions, for which a new potential model, [n(x),m], is proposed. The prototypical systems, methane and benzene, interacting with rare gases He, Ne, Ar, Kr and Xe, are analyzed as test cases. The method appears suitable for extensions to more complex systems, including modifications for treating ion–ion and ion–molecule interactions. 2004 Elsevier B.V. All rights reserved.
1. Introduction Model simulations of molecular aggregates, both in the gas phase and in the condensed matter conditions of interest in several applications, including biochemistry, are carried out by molecular dynamics methods, leading to the interpretation of macroscopic properties using the tools of statistical mechanics. Often the limiting factor in the ability of these computations to model real systems is the lack of a realistic description of the interaction potential energy (force field). Specifically, a challenging problem concerns the intermolecular part of the potential energy, arising from Ônon-chemicalÕ (often referred to as Ônon-bondedÕ) interaction components, which includes exchange-repulsion at short-range and electrostatic, induction and dispersion contributions at long-range. Current methods are based on the assumption of the pairwise atom–atom-like additivity of the in*
Corresponding author. Fax: +39-075-585-3864. E-mail address:
[email protected] (D. Cappelletti). 1 On leave from the Dept. de Quimica Fisica, Universitat Barcelona, Barcelona, Spain. 0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.06.100
teraction potential energy. These include Ôunited atomÕ descriptions, where the interaction sites coincide with the various functional groups of the molecule, the main advantage being the reduction of the number of centers considered, and Ôall atomÕ descriptions where each constituent atom is explicitly taken into account. In this context, three basic issues should be addressed: (i) the need of including many-body effects in the interaction; (ii) the parameterization of the potential energy function, which must be simple and at the same time reliable in a wide range of configurations (both radial and angular); (iii) the development of a variety of tools to anticipate the relevant features of the interaction, which relate to the basic parameters appearing in the expression for the potential energy surface. This last general issue is crucial, since experimental investigations and ab initio calculations are reliable only for a limited number of systems [1,2], while the use of combination (mixing) rules appears appropriate only in particular cases [3]. Some of us, in 1985 [4], proposed a new method to scale intermolecular interactions, employing a unique fundamental physical property of the interacting
38
F. Pirani et al. / Chemical Physics Letters 394 (2004) 37–44
particles, namely the polarizability, to account both for the strength of induced dipoles (i.e. the attraction: an idea already well known although difficult to quantify) and also for the average atomic and molecular size (exchange-repulsion). Subsequent studies [5], demonstrating the existence of a linear correlation between average polarizability and atomic or molecular volumes, confirmed the validity of these basic ideas. The availability of additional data from new experiments permitted to put this method into firmer grounds [6] and explicit relationships between average polarizability and important potential features (like the position of the potential minimum, the well depth and the long range dispersion coefficient) were proposed to describe the interaction energy in systems involving atoms and simple molecules. The effects of induction and electrostatic contributions were also included [7] permitting an extension to systems involving ions. In a successive contribution, Halgren [8] suggested a similar approach to determining atom–atom potential parameters, where the role of the atomic polarizability was further exploited, and where the direct implications for molecular mechanics force fields were also assessed. These methods [6–8] are more generally applicable than the older combination rules, whose validity is limited to specific families of systems. The extension of their practical usefulness, has stimulated the present work, addressed to a generalization of the use of polarizability as scaling factor for the intermolecular interactions. Recently, we have demonstrated [9] that for atom– molecule systems the main features of the intermolecular interaction can be obtained, at a given geometry of the molecular aggregate, by a combination of atom–molecular bond pairwise contributions, without explicit reference to the behavior of the atoms in the molecules (this has been accounted for by others [10,11]). In particular the concept of bond polarizability additivity, first introduced by Denbigh [12] and subsequently discussed by others (see, for instance, Smith and Mortensen [13]), has been exploited to represent both the molecular repulsion, in terms of a size which is mainly ascribed to the molecular bonds nearest to the approaching probe atom, and the molecular attraction, as due to many dispersion centers on the molecular frame. The distribution of the dispersion on the molecular frame allows one to go beyond the simple atom–atom additive picture of the interaction and to include effectively higher order contributions. This approach directly accounts for the specific role of the ÔchemicalÕ interaction between the two atoms in the bond and is founded on the assumption that the properties of any bond in a molecule (binding energy, bond length and bond polarizability) are not significantly influenced by the presence either of other bonds or of external particles generating a perturbative polarization [12,13].
In the present work we further develop these ideas making an explicit use of the bond polarizability tensor components (each one tied to its reference system), to account more directly for the anisotropy in the electronic charge distribution in the molecular frame. These considerations lead to a new formulation of the potential energy surface (PES), expressed as additive atom–bond pair potentials. For an accurate description of the pair potential radial and angular dependencies we introduce, in Section 2, a novel anisotropic two-body function. In the case of closed shell atomic species (either neutral or ionic) interacting with hydrocarbon molecules (both aliphatic and aromatic) the full PES can be represented in a compact analytical form (Section 3). We report the potential parameters for the two-body (atom–bond) interaction terms concerning rare-gas–CH and raregas–CC aromatic bond, and we discuss in Section 4 the properties of rare-gas–methane and rare-gas–benzene systems whose PESs can be built as a sum of atom-bond contributions. Methane and benzene have been chosen because prototypical examples of aliphatic and aromatic hydrocarbons for which some PESs exist in the literature for comparison. The present work should be intended as the presentation and test of a method which will be applied in the near future to more challenging and interesting cases, such as those involving ions interacting with aromatic molecules.
2. Potential model A large number of flexible and accurate two body potential functions, suggested by semiempirical methods [14,15] or built up empirically combining different functions at different distance ranges [16] have been extensively used to fit scattering and spectroscopic data, especially when analyzed simultaneously or in combination with other microscopic or macroscopic properties. Such functions provide accurate values for the potential energy but pay the price of an elevate number of parameters and often of relatively long times needed for their computation. By contrast, molecular mechanics methods exploit a few, very simple and computationally ÔlightÕ, potential functions, essentially the Buckingham (exp,6) and the Lennard–Jones (LJ) in the (12,6) version (sometimes (9,6) or buffered (14,7) forms have been proposed [8]). These functions, although satisfactorily reproduce the features of the potential well, exhibit known limits in describing both the short-range repulsion and the long-range attraction. In particular the LJ (12,6) is too repulsive at short-range and asymptotically overestimates the attraction as much as a factor two. Several attempts have been carried out in the past in order to mitigate the LJ model limitations. A noticeable
F. Pirani et al. / Chemical Physics Letters 394 (2004) 37–44
correction has been suggested by Maitland and Smith with their [n(x),6] function [17]. Such a modification involves the use of an exponent n(x) of the repulsive-term which varies with the reduced distance x = r/rm, where r is the two body distance, rm the position of the potential well and its depth: " nðxÞ 6 # V ðrÞ 6 1 nðxÞ 1 ¼ f ðxÞ ¼ : ð1Þ e nðxÞ 6 x nðxÞ 6 x The same authors propose n(x) = 13 + c(x1), where c is an additional adjustable form parameter, depending on the specific case at hand. Investigations have been carried out with 2 6 c 6 10 [17,18]. Fig. 1 compares the behaviors of [n(x),6] and LJ (12,6) potential models. As
800 600 400
can be seen, a significant improvement in the long range attraction is achieved with the [n(x),6] function using a value of c = 10, condition which however generates an excessive fall off of the repulsive wall. Lower c values provide a more reliable short range behavior of the [n(x) 6] function, but determine a too large attraction. To be more quantitative it can be reminded that the analysis of total cross-section measurements, performed on several neutral–neutral systems, indicates [6] that f(x = 1.7) 0.057 ± 0.006, while values calculated with LJ (12,6) and [n(x),6] with c = 10 and 4 are, respectively, 0.083, 0.059 and 0.066. It is worthy to stress at this point that the incorrect representation of the long-range attraction can strongly affect the description of large systems, such as those of chemical interest, since most contributions to the intermolecular potential arise from long-range interacting pairs. We propose here a modification of the Maitland– Smith extension of the Lennard–Jones potential, to be denoted [n(x),m], which represents a further improvement in the directions discussed above and whose applicability is extended also to include ion–neutral and ion–ion systems. The [n(x),m] function is defined as V ðrÞ ¼ f ðxÞ e "
200 0 0.6
0.8
1.0
-0.2
[n(x),6] (γ=10) [n(x),6] (γ=4)
-0.4
nðxÞ m # m 1 nðxÞ 1 ¼ ; nðxÞ m x nðxÞ m x
LJ(12,6)
-0.6
nðxÞ ¼ b þ 4x2 :
-0.8
1.0
1.2
1.4
0.00
-0.01
-0.02
2.0
2.5
ð2Þ
where as before, m = 6 for neutral–neutral systems, while m = 4 for ion-induced dipole, m = 2 for ion-permanent dipole and m = 1 for ion–ion cases. To meet most of the requirements discussed above we introduce, in the n(x) term, a more pronounced dependence on x:
[n(x),m] (m=6, β=9)
f(x)
39
3.0
x Fig. 1. The reduced form of the potential f(x) is plotted against x = r/ rm. Comparison between [n(x),6] (with c = 10 and 4), LJ (12,6) and [n(x),m] (m = 6 and b = 9) functions.
ð3Þ
Here, b is a constant parameter within wide classes of systems, and depends on the nature and hardness of the interacting particles. Specifically, we recommend b = 10 for systems involving rare gas atoms, alkali ions and alkaline earth dications; b = 9 for systems involving less hard neutral atoms and halogen ions; b = 8 is suggested for soft and strongly polarizable atoms (alkali) or for those cases where the long-range attraction is dominated by electrostatic forces like the ion-permanent dipole and ion–ion components. This model is simple computationally, also regarding its derivatives, and is defined in practice by only two parameters. It is capable to cover an ample phenomenology and exhibits further useful features, discussed below. The f(x) calculated with m = 6 and b = 9 (Eqs. (2) and (3)) is plotted in Fig. 1 for a comparison in a wide x range with the LJ (12,6) and [n(x),6] functions. Fig. 2a shows the reduced form of the [n(x),m] potential (with b = 9) for m = 6 (induced dipole–induced
40
F. Pirani et al. / Chemical Physics Letters 394 (2004) 37–44
f(x)
(a)
72, coincident with those of a LJ (12,6). When m is given the values 4, 2 and 1, the r/rm and k(x = 1) values read 0.87 and 48, 0.83 and 25, 0.79 and 12, respectively. The use of b = 8 in a [n(x),m] function with m = 1 provides 0.78 and 11. In this case the reference (alkali halides) values are 0.76 and 8 [19]. As a final step we now focus on the angular dependence in anisotropic systems. The atom (or ion)-homonuclear diatom case is easily accounted for by defining angular dependent e(a) and rm(a) parameters, in Eq. (2). Specifically,
0.0
-0.5
-1.0 1.0 (b)
2.0
3.0
50
n(x)
40
20 10
(c)
2.0
3.0
90 m=6 m=4 m=2 m=1
k(x)
60
ð4Þ
rm ðaÞ ¼ rm? sin2 ðaÞ þ rmk cos2 ðaÞ;
ð5Þ
where e^, ei, rm^ and rmi are, respectively, the well depth and location for the parallel and perpendicular approaches of the probe atom (or ion) to the diatom and a is the Jacobi angle, defining the geometry of the system. We note that for future extensions to the case of heteronuclear diatoms, Eqs. (4) and (5) must be modified in order to include additional angular terms (for example depending on cos(a)) in order to account for the different behavior at a = 0 and a = p. In the next section we will exploit this anisotropic potential model to account for the specific role of any atom-bond pair within the overall interaction in an atom–molecule system.
30
1.0
eðaÞ ¼ e? sin2 ðaÞ þ ek cos2 ðaÞ;
30
3. The atom–bond pairwise additive representation of the PES
0 1.0
1.5
2.0
x Fig. 2. Features of the [n(x),m] potential function (b = 9): (a) Dependence of f(x) on m. m = 6 (continuous line, induced dipole-induced dipole), m = 4 (dashed line, ion-induced dipole), m = 2 (dotted line, ion– permanent dipole) and m = 1 (dash-dotted line, ion–ion). (b) behavior of the n(x) = b + 4x2 function. (c) behavior of the reduced force constant k(x).
dipole), m = 4 (ion-induced dipole), m = 2 (ion-permanent dipole) and m = 1 (ion–ion) cases. The behaviors of n(x) (Eq. (3)) and of the reduced force constant (k = d2f(x)/dx2) are plotted in panels (b) and (c) of the same figure. These results suggest that, as expected, the potential well area widens as m decreases, that is as the strength of the attraction increases. The location of the inner zero of the potential, r, and the force constant k at the minimum of the potential energy (i.e. the curvature of the potential function) k(x = 1), are important ÔshapeÕ parameters of the interaction. For a [n(x),m] with b = 9 and m = 6, the values of the r/rm ratio and k(x = 1), are, respectively, 0.89 and
Pairwise additive functional forms derive from the many body expansion formulation of the potential energy U X ð1Þ X ð2Þ X ð3Þ U¼ Ui þ Uj þ Uk þ ; i
j
(1)
k
where the U is the potential energy of the separated atoms while U(2), U(3) and higher order terms define both the intramolecular and intermolecular potential energy. In the standard atom–atom pairwise potentials the expansion is truncated to the second order terms. This approximation is quite crude especially for the intermolecular part. The desirable inclusion of three-body effects typically involves a large number of anisotropic U 3k terms. The use of atom–bond interaction components allows us to maintain the additive formulation of the intermolecular potential V and naturally includes high-order effects [9]. A more realistic description of the interaction anisotropy is then achieved since such a formulation indirectly accounts for the electronic charge distribution along the molecular frame.
F. Pirani et al. / Chemical Physics Letters 394 (2004) 37–44
Accordingly, V is defined as combination of atom– bond components Vab, represented with the anisotropic [n(x),m] potential model. Specifically X V ¼ V ab ðrab ; aab Þ; ð6Þ
Table 2 ˚ for three Well depth De in meV and equilibrium distance Re in A relevant cuts of the CH4–rare gas PES. System
4. Results and discussion The potential parameters for the rare-gases–CH and rare-gases–CCar pairs, obtained with the procedure illustrated in [9] and using bond polarizability data from [12], are given in Table 1. Tables 2 and 3 report the values of binding energy De and equilibrium distance Re (from the molecular center of mass) for some limiting geometries of rare gas–methane and rare gas–benzene systems, calculated using the parameters of Table 1 and Eqs. (2)–(6). The b parameter of the [n(x),m] function has been fixed to 10 in all cases, for the reasons discussed above. The tables report also data selected from literature for a comparison. The Table 1 Potential parameters for the CC aromatic (CCar), and CH–rare-gas pairs System
Rm^
Rmi
e^
CCar–He CH–He CCar–Ne CH–Ne CCar–Ar CH–Ar CCar–Kr CH–Kr CCar–Xe CH–Xe
3.583 3.234 3.629 3.316 3.879 3.641 3.997 3.782 4.163 3.976
4.005 3.480 4.015 3.549 4.189 3.851 4.284 3.987 4.425 4.175
0.860 1.364 1.706 2.544 3.895 4.814 4.824 5.667 5.654 6.233
ei
0.881 1.016 1.862 1.960 4.910 3.981 6.365 4.781 7.849 5.384 ˚ and in The perpendicular and parallel component of Rm and are in A meV, respectively.
Vertex Re
b
where rab is the distance of the incoming probe atom a from the dispersion center localized on the b bond and aab is the angle between rab and the bond axis. From another point of view the present method represents a polyatomic molecule as an ensemble of diatoms, coinciding with the bonds, each one playing the role of an anisotropic interaction center. The CCar (carbon carbon aromatic) and CH bonds, in the hydrocarbon molecules of interest here, are assumed to be independent diatomic sub-units having electronic charge distributions of near cylindrical symmetry. Moreover, in each atom–bond pair the reference point is set to coincide with the geometric bond center, since the dispersion center and the bond center coincide for CCar and are almost coincident for CH. In other words, the electronic charge distribution around the CH bond shows an ellipsoidal shape whose center is close to the bond center [9] (see also [20]).
41
Face De
Re
Edge De
Re
De
CH4–He
3.84 4.01 4.12
2.51 1.92 1.79
3.36 3.41 3.36
4.78 4.32 3.36
3.45 3.53 3.59
4.26 3.83 3.11
Present [9] [21]
CH4–Ne
3.91 4.02 4.08
4.89 4.16 4.11
3.44 3.48 3.44
8.97 8.50 6.51
3.53 3.58 3.61
8.03 7.68 6.61
Present [9] [22]
CH4–Ar
4.18 4.18 4.18
10.5 10.6 10.5
3.77 3.75 3.72
17.2 17.7 17.9
3.85 3.83 3.92
15.7 16.5 14.4
Present [9] [23]
CH4–Kr
4.31 4.28
12.89 13.6
3.90 3.88
20.36 21.4
3.99 3.95
18.71 20.2
Present [9]
CH4–Xe
4.47 4.48
15.01 14.91
4.05 4.10
24.11 22.56
4.14 4.17
22.14 20.89
Present [9]
agreement appears to be quite good considering the different origin and uncertainty of the information. We estimate the uncertainty associated with the present De and Re to be, respectively, 10% and 3%. For the CH4–He and CH4–Ne systems the present method seems to overestimate De and underestimate Re although the anisotropy appears to be consistent with the experiments. In these cases the experimental information comes from the analysis of diffraction oscillations [21,22] in the differential cross-sections, an observable sensitive to the anisotropy and to the onset of the repulsion but scarcely sensitive to the strength of the interaction. For CH4–Ne, a successive simultaneous analysis of total and differential cross-section results [24] indicates that De must be increased of 15% and Re decreased of 2%. For the CH4–Ar system a thorough comparison between the present PES and the analytical fit to the ab initio calculated data of [23] is shown in Fig. 3. Specifically, in Fig. 3a the two spherically averaged PESs are compared with the experimentally determined isotropic potential of [24]; Fig. 3b shows the angular dependence of the interaction energy while Fig. 3c compares the radial dependence of the potential V for three limiting configurations of the complex. The overall agreement is very good. The average interaction energy for the CH4–Kr and CH4–Xe systems, obtained from the present anisotropic PESs appears to be in full agreement with the experimentally determined isotropic potentials ˚ ) and De (meV) of [24]. In particular the present Re (A are 4.10 and 16.2, 4.28 and 18.4, while the experimental values read 4.02 and 17.0, 4.24 and 19.6 for CH4–Kr and –Xe, respectively. Moreover, also the average interactions of C6H6–He, Ne and Ar compare well with the
42
F. Pirani et al. / Chemical Physics Letters 394 (2004) 37–44
Table 3 ˚ for three relevant cuts of the C6H6–rare gas PES Well depth De in meV and equilibrium distance Re in A System
Out-of-plane (perpendicular) Re
In-plane vertex De
Re
In-plane edge De
Re
Ref. De
C6H6–He
3.25 3.20 3.23 3.17 4.99 3.31
9.8 10.2 10.2 – 5.5 8.4
4.79 4.88 4.93 – 7.29 –
5.3 5.1 4.4 – 0.7 –
5.23 5.4 5.49 – 7.29 –
3.2 2.6 2.2 – 0.7 –
Present [32] [9] [25] [26] [27]
C6H6–Ne
3.30 3.30 3.31 3.30 3.32
19.5 20.2 20.2 18.7 19.8
4.86 4.87 4.96 – –
10.3 10.6 9.1 – –
5.28 5.34 5.47 – –
6.3 5.6 4.9 – –
Present [32] [9] [28] [29]
C6H6–Ar
3.57 3.58 3.59 3.58 3.50 3.55 3.41 3.55
44.1 43.3 43.5