13876
J. Phys. Chem. 1995,99, 13876-13885
Investigation of Supramolecular Systems by a Combination of the Electronegativity Equalization Method and a Monte Carlo Simulation Technique Helge Toufar," Bart G. Baekelandt, Geert 0. A. Janssens, Wilfried J. Mortier, and Robert A. Schoonheydt Centrum voor Oppervlaktechemie en Katalyse, Katholieke Universiteit Leuven, Kardinaal Mercierlaan 92, B-3001 Heverlee, Belgium Received: April 3, 1995; In Final Form: June 26, 1995@
The electronegativity equalization method and the charge sensitivity analysis, both in the atom-in-molecule resolution, are adapted to complex supramolecular systems. The combination of these methods with the Monte Carlo simulation technique results in a tool that allows one to determine simultaneously the equilibrium (intermolecular) structure as well as the charge distribution of such systems. Furthermore, information about the reactivity of molecules in interaction with a complex environment, e.g. a solvent, can be obtained. The possibilities of the method are illustrated with examples of water and cation-water clusters as well as the interaction of water with different kinds of zeolite framework moieties.
Introduction
The description of large supramolecular systems by means of computational chemistry is a present-day objective of research. A tremendous amount of work has been done, for example, to simulate the properties of bulk water by molecular dynamics and Monte Carlo technique^.'-^ Other topics under investigation in this respect are the behavior of water and other molecules (or atoms) within the pore systems of molecular the polymerization of silicic acid in aqueous solution?.l0 and the influence of solvent molecules on the properties of large organic molecules.' In most cases, these simulations use a two-piece force field consisting of a long-range electrostatic force and a short-range dispersion force. While the basic physical laws for the calculation of these forces are straightforward, the parameters actually needed for the calculations are not. The charges required to calculate the Coulomb type interaction are often formal charges or taken from ab initio calculations on single molecules and held constant throughout the simulation. Parameters used to calculate short-range interactions, e.g. of the Lennard-Jones type, are derived from atomic polarizabilitiesand adjusted in order to reproduce intermolecular distances and/or certain thermodynamic properties, thus varying primarily at the author's discretion. While the lack of reliable and adaptable charges is a minor problem, whenever the contribution of the electrostatic forces is small (e.g. in intramolecular optimizations of covalent bonds or intermolecular optimization of apolar molecules), it becomes a considerable problem in cases of polar molecules and, especially, if ionic species become involved. Under such circumstances the effect of the mutual polarization of the interacting species cannot be ignored. Recently, a number of publications addressed these problems at different levels of accuracy. A very promising approach for larger systems is to calculate the distribution of atomic charges in a system using the principle of electronegativity equalization (EE).I23l3 Rappe et al. proposed the use of a modified EE method (QEq)I4 to calculate atomic charges for molecular mechanics force fields and molecular dynamics methods. This method was later include in the universal force field (lJFF)l5 @
Abstract published in Advance ACS Abstrucrs, September 1, 1995.
0022-3654/95/2099-13876$09.00/0
and applied to a variety of compounds.I6 Van Duin et al. included charge calculations by EE into a hydrocarbon force field.I7 An especially interesting approach was taken by Rick et al., who used a dynamical EE method in molecular dynamics simulations of liquid water.'* Rather than enforcing an equilibrium charge distribution with equalized electronegativities, these authors consider existing electronegativity differences within molecules as electronic forces acting on the atomic charges. By assigning a fictitious electronic mass to the charges, these can be handled by the MD algorithm in the same way as inert masses, although on a much faster time scale. We present here a combination of a classical Monte Carlo simulation technique (MC) with the original electronegativity equalization method (EEM). MC allows the optimization of the arrangement of the molecules within a supramolecule. EEM provides the charges required for the calculation of the electrostatic forces. A major advantage of this combination is that an EEM-based sensitivity analysis (SA)I9 can be applied easily during a simulation, allowing conclusions to be drawn about the electronic reactivity of the supramolecular system or one or more of its parts upon the perturbation by charge transfer or by changes in the extemal potential. Methods The Electronegativity Equalization Method. For a detailed description of the EEM and its a p p l i c a t i ~ n ? ~its - ~predeces~ s o r ~ and , ~ ~its rigorous derivation in the density functional t h e ~ r y ? ~as- well ~ ~ as for the SA and analogous system^,^^-^' we refer to the original literature. From the density functional theory results the following expression for the electronic energy of a molecule in the ground state:
E @ ) = Q)+ J e ( r ) 41)dr
(1)
where F(e) is the sum of the kinetic energy and the electronelectron interaction energy, and V ( I ) is the nuclear-electron potential. By integrating this expression over atomic regions (using a spherical atom approximation), one obtains the molecular electronic energy as a sum of atomic contributions: 0 1995 American Chemical Society
Investigation of Supramolecular Systems
J. Phys. Chem., VoE. 99, No. 38, 1995 13877
where 4a is the net atomic charge of an atom q = 2 - N . The unknown functional F ( e ) was replaced by a second-order expansion of the energy as a function of the atomic charge of each atom. The expansion coefficients E*, and q* have been calibrated for several atom types so as to reproduce ab initio (STO-3G)charges (obtained by a Mulliken population analysis) and energies. The partial derivative of eq 2 with respect to an atomic charge 4a defines the electronegativity Xa of the atom a in a molecule:
x*,
The second derivatives of the energy with respect to the atomic charge, i.e. the partial derivatives of the electronegativity Xa with respect to the charge 48, are the hardness kernels 7 ~ :
allowed, resulting in a total equalization of electronegativities between all atoms and all molecules. A comparison of such calculations on small systems with ab initio quantum chemical calculations shows that the charge transfer between different molecules is o v e r e ~ t i m a t e d . ~An , ~ electron ~ added to one of the molecules in the system would spread over the entire supramolecule. This would indicate, for example, bulk water to be a perfect electron conductor. Since a partially restricted charge transfer is difficult to realize within the EEM, the total avoidance of intermolecular electron transfer seems to be the best approximation. The EEM energy expression for the supramolecular system contains a summation over all atoms in all molecules of the system:
where m is the number of molecules in the system, ni is the number of atoms in molecule i, q is the atomic charge, and R is the interatomic distance. The atomic electronegativities are then given by
x*,
We can identify the expansion coefficients E*, and T,I* with the energy, the electronegativity, and the hardness of a neutral (4 = 0 ) atom in a molecule, respectively.*’ The electronegativity is the negative of the chemical potential in the density functional theory. In a molecule in an equilibrium state, the chemical potential (and the electronegativity) is equal at every point within the molecule. This means, of course, that the electronegativities of all atoms in a molecule get equalized (electronegativity equalization principle):
X a = xp =
...= x, = x -
(5)
where is the global electronegativity of the molecule. A molecule consisting of n atoms is now completely described by n equations of type (3) and the constraint of the total charge of the molecule Zq = const, with exactly n 1 unknown quantities, namely, n atomic charges and the global electronegativity. This system of equations can be solved exactly in order to determine the atomic charges and the molecular electronegativity. More information about the electronic properties of a molecule can be obtained by taking further derivatives of the energy expression in eq 2 with respect to the number of electrons N ( 4 ) and the external potential u. The resulting quantities describe the response or the sensitivity of a molecule or parts of a molecule toward perturbations in N or u, e.g. the global hardness q = (ap/aN), = (ax/aq)”, the global softness S = (a4/ the local softness Sa = (a4a/&)u, the Fukui index fa = (aqa/aq)u,and the linear response kemel p$ = ( a q a / a u ~ )Such ~. a sensitivity analysis (SA) or charge sensitivity analysis (CSA) can be a useful tool for predicting trends in molecular reactivity. It is important to note here that the EEM was already successfully applied to zeolite^,^^-^^ their interaction with adsorbed and the interaction between cations and water molecules in the gas phase?’ but without taking the supramolecular character of these systems into account explicitly and by using estimated, fixed geometries. Similar restrictions apply for the use of sensitivity analysis m e t h o d ~ . ~ * , ~ ~ Extension of the EEM to Supramolecular Systems. The simplest way to investigate a supramolecular system by means of the EEM/SA is to consider the whole system as one single molecule. A free electron flow between the molecules is then
and the hardness kernels by
So far, there is no principle difference with classical EEM as described by eq 2-4. If the expressions of eq 8 are inserted in eq 7, the electronegativity equations for the complete system can be arranged to a characteristic matrix equation: In eq 9,
+
.. . . . . . rla,n,
‘lama, 0
.
0
*
%,a,
+
.
*
.
rln,n,
X
rla,n,
rla,a,
rln,a,
.
lln,,,n,
’ ’*
1lxq=X-
x
the hardness matrix is schematically outlined for a system consisting of two molecules i and m containing ni and n, atoms, respectively. The upper left and the lower right square of the matrix 7 represent the intramolecular interactions, while the upper right and the lower left parts describe the intermolecular interactions. In this form, the equation-could be used to calculate the equilibrium charges (ie. at x = 0) for the case that all molecules are coupled to an external electron reservoir of infinite softness:
-[SI x
-
x* = -
4eq
IS1 = 1v1-l
(10)
In this case, of course, it makes no difference whether any of the molecules were charged before. This can be proven by starting the integration implicitly done in eq 5 from an arbitrary charge distribution G,:
Toufar et al.
If no electron flow between the molecules is allowed, the sum of the atomic charges of each molecule will be constant:
While the Fukui indices sum up to 1 by definition, the g values sum up to 0. Due to the constraints applied, changes on any other atom in another molecule may affect the charge distribution within a molecule but not the total charge. (iii) The molecular hardness, i.e. the change of the electronegativity of a molecule i upon changing the total charge of this molecule:
n.
C q i a = qi = const a
Furthermore, the atomic electronegativities will equalize within all molecules but not globally, i.e. between the molecules: -
-
X. =%, iB =...=X. = X i # x la . =...=Xj#Xm la
IW
-
(15)
The constraints formulated in eqs 14 and 15 can be included into the matrix eq 9 as additional rows and columns, respectively. Correspondingly, the charge vector had to be extended by the vector of these equilibrium electronegativities, and the vector of the starting electronegativities has to be extended by the vector of the molecular charges:
(iv) The intermolecular correlation hardness, Le. the change of the electronegativity of a molecule i upon changing the total change of another molecule J :
Therefore, the charge on each atom in the system is given by n.
m
and the electronegativity of a molecule in the equilibrium state is
In eq 16, lrol represents a submatrix accoLding to the squares in eq 9. Each of the additional vectors l l i , lil, O i i , and 01; contains ni elements 1 or 0, respectively. The molecule index in the first position denotes column vectors; the same index in the second position denotes row vectors. The equilibrium charges are obtained similarly to eq 10 by inverting the matrix in eq 16 and multiplying this inverse from the right with the x*/q, vector. The elements of the inverted matrix can be identified by comparison of the coefficients in eq 17 and in the definition of the inverse matrix Irl x [VI-' = 111:
A similar interpretation was recently published independently by Korchowiec et ~ 1 . 4The ~ submatrices are built from the internal softness elements, representing the change of the charge of a atom due to the change of the electronegativity on another atom if the number of electrons in all molecules remains constant. The other elements are as follows: (i) The atomic Fukui indices Le. the change of the charge on atom a; upon changing the total molecular charge of i:
(ii) The atomic correlation indices, i.e. the change of the charge on atom a,upon changing the total molecular charge of another molecule j :
Obviously, eq 9 can easily be adapted to any situation which is chemically reasonable; for example, one or more molecules could be recoupled with the external reservoir by removing the corresponding condition lines; charge transfer could be allowed within a domain of molecules by collapsing their condition lines into one, and so on. Supramolecular Systems and Internal Hardness. An alternative way to describe systems with restricted charge transfer, i.e. systems in which at least one of the molecules has no charge transfer with an external reservoir or other molecules (closed molecule), is provided by the concept of internal hardness.47 The basic idea of this concept is that by not allowing charge transfer between molecules, the degrees of freedom are reduced by 1, that is, a system which was completely determined by n variables (e.g. the atomic charges) is now determined by n - 1 variables. Transferred to a supramolecular system, this means that the number of independent variables is reduced from n (total number of atoms in the system) to n - mnct,where mnct is the number of molecules without charge transfer to the environment. The quantity used to d e t e d n e the restricted system may be chosen according to the goal of the analysis. In a closed molecule one atom may be arbitrarily selected as a reference atom. The system is then described in terms of the charge transfer from this atom to the other n - 1 atoms in the molecule.47 Thus, the internal quantities of a supramolecular system are given by
Investigation of Supramolecular Systems
where di refers to the reference atom in molecule i. Equation 25 gives the interaction of two atoms located in two closed molecules i andj. In eq 25a one of the molecules is closed (i); the other, open to the extemal reservoir. The transition between the extemal and the intemal hardness matrix is schematically presented in Figure 1. If one intends to study the bond polarization by intermolecular interactions, it might be more convenient to use the differences of the charges on bonded atoms as variables; for studies of functional groups, the systems could be formally split into subgroups and described by the charge differences between these groups. Consider, for example, the heterogeneous dissociation of a molecule. Here we are interested in the charge separation Aq, = &a(x) - &a(Y) between two fragments X and Y of a molecule. A charge separation close to 0 corresponds to a covalent linkage between the two fragments X-Y; a charge separation of Aqr = 2, to the formation of an ion pair X+Y-. The charge distribution in a closed molecule is completely determined by n - 1 independent charge separations. Deriving the energy with respect to these charge separations results in the following expressions:
where ni(0 is the number of atoms in the molecule both fragments belong to, nFa is the number of atoms in the fragment the atom a belongs to, and SFa is the sign of the electron flow with respect to this fragment, SFa = 1 corresponds to an electron donor; S F a = - 1, to an electron acceptor; and $Fa = 0, to fragments which are not affected by the electron flow. Obviously, eqs 24 and 25 are special cases of these general expressions, obtained if one assumes that the dependent atom d is the donor fragment, that the independent atom a is the acceptor fragment, and that the other atoms in the molecule are not affected. If all molecules in the system are neutral, the equilibrium charges could be calculated by using the intemal quantities with eq 10. However, if there are charged molecules in the system, an additional step is required, since the molecular charge cannot be changed within the intemal system; that is, to get the starting electronegativity vector for an equation of the type of eq 11, one has to perform an integration within the extemal system fist. Since the initial distribution of the molecular charge is not important, it is most convenient to put this charge completely on the reference atoms. For the elements of the starting vector a general expression results:
J. Phys. Chem., Vol. 99, No. 38, 1995 13879
Figure 1. Transformationof the external to the corresponding internal hardness matrix for a two-molecule system. m
with j running over all molecules m containing a dependent reference atom d, thus allowing one to calculate the equilibrium charges:
The elements of this intemal softness matrix Is'/are identical to the corresponding elements of the matrix in eq 17. The missing softness kernels for the reference atoms can be obtained by applying the closure condition, Le., by exploiting the fact that in a closed molecule the sum of all changes of atomic charges has to be 0. Again, every combination of closed and open molecules can be realized by a proper combination of intemal, extemal, and mixed submatrices. Although both approaches presented here are mathematically completely adequate, the intemal one has a major advantage from a computational point of view. While the addition of conditions to the extemal hardness matrix increases the dimensions of the matrix by mnct,the use of the intemal approach reduces the dimensions by the same amount. In the case of a water cluster, for example, this means a matrix of 2m x 2m instead of 4m x 4m, thus saving 75% of the computer memory and nearly 90% of the CPU time required to invert the matrix. On the other hand, the use of the extended form provides a lot of additional information in one step (e.g. the hardnesses and Fukui indices), which requires additional effort within the strict intemal system. Thus, the choice depends on the aim of the calculation and the special character of the system. Normal Mode Analysis in SupramolecularSystems. Additional information can be obtained by a decoupling of the softness tensor, performed by a rotation of lsdl into its principal axes repre~entation:~g
lUIT x
lsdl x
IUI = IsI, lUIT x IUI = I l l
(30)
The row vectors of IUI,which are in fact the eigenvectors of the system, represent the independent charge redistribution channels upon polarization and/or charge transfer. The corresponding eigenvalues, i.e., the diagonal elements of IsI, give the softness of the corresponding modes. Thus, the eigenvectors with the highest eigenvalues are the softest (energetically most favorable) normal modes. Analysis of these modes can be useful in the discussion of multiple site r e a c t i ~ i t i e sand ~ ~in the design of an optimal environment in order to favor reactions along certain polarization channels, thereby manipulating the selectivity of a reactive system. We notice here that the orthogonal transformation has to be performed with the fully extended intemal softness matrix, Le. containing the softness kernel elements of the matrix in eq 13
Toufar et al.
13880 J. Phys. Chem., Vol. 99, No. 38, I995 Full Integration of EEMlSA and MC
Monte Carlo simulation Energy
TABLE 1: Parameters Used in the EEM/Monte Carlo Simulation, the Electronegativity of the Neutral Atom in a Molecule, p,the Hardness of the Neutral Atom in a Molecule, TI*,and the Nonbonded Ionic Radii, rHsa
x*, eVle
element
Partial Integration of EEM/SA and MC
Energy calculation with rigid charges
4.408 77
q*, eVle2 13.773 24
8.5 32.421 05
11.082 87 90.004 88
-2.239 52 1.331 82
7.672 45 6.492 59
rHS,
A
0.6 0.7 0.4 1.3 1.3 1 .o 0.9 0.8 1.3
a EEM parameters for the alkali and earth alkali metals are not yet available. Since these atoms occur in the systems under consideration as isolated ions with constant charge only, the parameters are not required for the calculation.
intermolecular geometry
5 s ? (2)
b7
0.187
D
-0&209
0.195
W
0.202
-0.404
-0.402
0.210
t3
-0391
0 200
u . 2 (3)
Figure 2. Schematic flow diagram of the geometry optimization
LCiQ (2)
0 179
procedure.
0.199
-0 428
or the internal softness matrix of eq 27 after supplementing the elements of the reference atoms. Diagonalization of mixed matrices, as mentioned above, results in normal modes corresponding to a charge transfer to the open molecules, while for closed molecules only a polarization occurs. Force Field and Geometry Optimization. The optimization of the intermolecular geometry of the system is done by a classical Monte Carlo s i m ~ l a t i o n .There ~ ~ are two possibilities of integrating EEM into the simulation. If a full integration of EEM into the Monte Carlo formalism is applied, the total energy of the system is calculated with eq 6 after each Monte Carlo step. However, this may become very time consuming for large systems since the time required to calculate the charge increases with nearly the third power of the number of atoms. In such a case the charges were considered to be constant for a number of steps, e.g. equal to the number of particles in the system, and only the change in the intermolecular electrostaticinteraction energy of the moved molecule was considered as a criterion for the acceptance of the Monte Carlo step (partial integration). In Figure 2 a schematic flow diagram of the simulation procedure is presented. For all examples given in this paper a full integration was applied except for the Li-64H20 cluster. The internal geometries of t h e molecules were considered to be rigid during the simulation. They were taken from reliable experimental data wherever available 01' from ab initio or semiempirical calculations given in the literature. Due to the atomic resolution .of the EEM approach, the method is missing interactions caused by intra-atomic polarization effects, Le. forces caused by a deviation of the electron distribution in the atomic region from the spherical approximation. However, since these forces are considered to play a minor role in the systems under consideration, they are neglected for the moment, until a consistent way is found to estimate them within EEM. Repulsion forces are mimicked by a simple hard sphere model. The complete system was programmed in C++ and implemented on a DEC Alpha AXP 3000-600 workstation with an Open VMS operating system. Its limits are governed by the computer resources, mainly by the physical memory. A RAM
-0359
&
0229
0227
(3)
lib
0201 -0 428
m (3)
c".$"'7
0 228
0 199
553 (1)
Figure 3. Geometries of the water monomer (D), the chain trimer (A), a symmetrical energy maximum state (B), and the cyclic trimer (C) with the charges on each atom and the electronegativity (underlined) for each molecule. Values for symmetric positions are omitted. The numbering of the molecules (in parentheses) corresponds to the numbers in Tables 2-5.
of 64 Mbytes is sufficient to run systems consisting of about 1000 atoms within a reasonable time. A commercial molecular modeling package was used as a graphical interface to the program.51 The parameters actually used are given in Table 1. Examples and Discussion
Three topics are chosen as illustrative examples. The examples are kept deliberately simple in order to allow a discussion of the special features of the method without presenting excessive sets of numbers. Detailed investigations of more complex systems will be the subject of future publications. Water Clusters. First we will discuss the behavior of a water trimer (Figure 3). The Monte Carlo simulation gives a cyclic trimer (C) as the most stable configuration. Notice that the structure of case (C) is a snapshot after reaching equilibrium. To illustrate the role of the intermolecular geometry on the system's properties, we compare this structure with two nonoptimized structures, namely, a chain (A) as a local, nonsymmetrical minimum state and a symmetrical energy maximum state (e-max, B). The calculated atomic charges and molecular electronegativities are presented in Figure 3 together with the geometries in comparison with an isolated water molecule. The electronegativities and the hardnesses for these clusters are also summarized in Tables 2-4. In the symmetrical (cyclic) clusters B and C, the electronegativities are nearly equalized (exactly equalized for the
J. Phys. Chem., Vol. 99, No. 38, 1995 13881
Investigation of Supramolecular Systems TABLE 2: Global Molecular Electronegativities,xi, Global Molecular Hardnesses, q d , and Global Molecular Interaction Hardnesses, vu, for the Chain Trimer (Nonoptimized) mol. (i)
x,,eVle V I , , eVle2
q2,, eV/e2 73,eVle2 g", eV/e2
1
2
3
5.19 17.25 4.89 2.85 24.99
5.59 4.89 17.06
6.22 2.85 5.64 17.24 25.73
5.64 27.69
TABLE 3: Global Molecular Electronegativities,xi, Global Molecular Hardnesses, vu, and Global Molecular Correlation for the Energy Maximum Trimer Hardnesses, ti', (NonoDtimized) mol. (i)
x,,eVle VI,, eVle2
m,eVle2 m,,eVle2
qlCc,eV/e2
1
2
3
5.24 17.25 4.36 4.36 25.97
5.24 4.36 17.25 4.36 25.97
5.24 4.36 4.36 17.25 25.97
TABLE 4: Global Molecular Electronegativities,xi,Global Molecular Hardnesses, l a , and Global Molecular Correlation Hardnesses, vi;,for the Ring Trimer (Optimized) mol. (i)
x,,eVle V I # ,eVle2
w,,eVle2
m,, eVle2
vy, eVte2
1
2
3
5.63 17.06 5.26 5.30 27.62
5.60 5.26 17.09
5.51 5.30
5.50
16.98 27.78
27.85
5.50
idealized energy maximum cluster B). The slight differences between the water molecules in cluster C are due to the dynamic character of the Monte Carlo simulation. One would certainly get a total equalization by switching to an exact energy minimization algorithm at this point. The linear cluster A exhibits the opposite behavior, namely, a splitting of electronegativities with extreme values for the outermost molecules. In most cases the molecular electronegativities are lower than for the isolated molecule. The changes in average electronegativity are readily understood from eq 2: a supplementary negative extemal potential term, originating from the larger negative charge on the oxygens, will lower the average electronegativity. The hardness values of the different clusters show a very similar behavior. Again we find an equalization of the molecular hardnesses for symmetric clusters and a splitting for the chain. Generally, the molecules in the clusters are softer than an isolated water molecule (17 = 17.38). This is not surprising since an interaction with neighboring molecules will facilitate any change of the charge on a molecule. For the same reason, the molecules in the most stable configuration are the softest. On the other hand, the correlation hardnesses vij are the highest for the ring cluster; that is, a change of the number of electrons in one molecule in this cluster would change the electronegativity of this very molecule less than in any other configuration (vii lower), but would affect the other molecules in the cluster more (vu larger). The hardness of each of these molecules cannot be considered therefore without its environment. In the case of homogeneous systems, Le., systems where all components become identical when placed at infinite distance from each other, it is reasonable to propose that the sum of all effects caused by the change in one molecule is a measure of the reaction of the cluster to the change of one of its components:
m
rP' = zrij j
We could define this quantity as the concerted hardness of a molecule in a cluster. In this way, in the most stable cluster the concerted hardnesses of all molecules are nearly equalized and are higher than in any other configuration. This indicates a stability criterion for supramolecular systems in the sense of the maximum hardness p r i n ~ i p l e . ~ ~ - ~ ~ For the chain, the result is heterogeneous, giving a high concerted hardness for the central molecule and low values for the end members. This raises the question of definition of the hardness of the cluster as a whole. It seems reasonable that the weakest point of a system, i.e. the molecule with the lowest concerted hardness, determines the hardness of the system itself. Applying this criterion, the following ranking of hardnesses is obtained: ring (C, 27.62 eV/e2) >> e-max (B, 25.97 eV/e2) > chain (A, 24.99 eV/e2) Also other methods of representing the hardness of the clusters point toward maximum hardness principle. For instance, the hardness obtained when the whole system is considered as one molecule is ring (C, 9.25 eV/ez) >> chain (A, 8.66 eV/e2) = e-max (B, 8.65 eV/e2) and the arithmetic average of the concerted hardnesses of a system is ring (C, 27.75 eV/e2)>> chain (A, 26.14 eV/e2) > e-max (B, 25.97 eV/e2). In all these cases the most stable cluster is the hardest. Table 5 illustrates for one example (molecule 1 in the ring cluster C) what happens at the atomic level when an electron is added to one of the molecules. Within the affected molecule itself, the charge will preferably flow to both atoms, which belong to the ring (0,f = 0.336, and H1, f = 0.337), where it is most stabilized, as indicated by the Fukui indices (column 1). In the other molecules, the electrons are redistributed in order to realize this stabilization. They flow away from the atoms adjacent to the affected molecule (Le., negative g values, 0 in molecule 2 with g = -0.095 and H1 in molecule 3 with g = -0.097) toward the outer hydrogens (H2 in molecule 2 with g = 0.085 and H2 in molecule 3 with g = 0.058). The g indices of the two remaining atoms are smaller (H1 in molecule 2 with g = 0.01 and 0 in molecule 3 with g = 0.039) since these atoms are nearest neighbors to each other but in adjacent molecules, so that changing the charge on both in the same direction, as indicated by the sign of the g indices, necessarily corresponds to a destabilization of the system. It is interesting to investigate how the correlation between the molecules in a cluster affects the polarization channels, Le., the preferred ways of charge redistribution (so-called electronic vibrations) upon a nondirected perturbation of the system. As described above, the directions of these polarization channels are given by the eigenvectors of the extended intemal softness matrix (eq 21). The analysis of these eigenvectors is somewhat cumbersome for larger systems, but it tums out that both polarization modes present in the isolated water molecule (the symmetrical and the antisymmetrical 0-H polarization (H-H polarization)) survive in the clusters as well. Instead of two
Toufar et al.
13882 J. Phys. Chem., Vol. 99, No. 38, 1995 TABLE 5: Atomic Fukui Indicesf,l of the Atoms in Molecule 1 of the Ring Trimer (Optimized) and the Corresponding Correlation Indices on the Atoms of Molecule 2 ka12) and Molecule 3 ka13) of the Same Cluster
a 0
H1 H2
2
1 0.336 0.337 0.327
1.ooo
2 -0.095 0.010 0.085
3 0.039 -0.097 0.058
0.000
0.000
distinct modes we have now two groups of modes, one mainly described by 0-H polarizations in all molecules of the cluster; the other one, by H-H polarization in all molecules. The differences between the modes within a group are determined by various contributions of the individual molecules. This is illustrated in Figure 4, where the inverted eigenvalues (see eq 30) of the polarizational modes for the isolated water molecule, the three trimers discussed above, and two larger clusters (8 and 16 molecules) are given. It can be seen that the eigenvalues form a band, the width of which increases with the cluster size. The energy level of the softest mode decreases with the number of molecules for all clusters, although the effect is less pronounced for the nonoptimized trimers. This mode corresponds to a mutual enhancement of the polarizations in all adjacent molecules. Since the softest modes will contribute most to the overall behavior of the system, this means that the molecules in a cluster are more susceptible to polarization than an isolated molecule. There is some similarity with the downshift and the broadening of the IR bands, especially of those bands which correspond to bond stretch vibrations, when changing from isolated water toward water clusters and bulk water.55 This is a remarkable result in light of the recently described mapping relation between bond vibrational and electronic vibrational modes.56 Ion-Water Clusters. When a cation is added to a small water cluster, the resulting structure is of course dominated by the positive potential around this ion. The electronegativities increase dramatically and so does the polarization of the water molecules (Figure 5 ) . The rather large, singly charged ion Cs+ increases the atomic charges by more than 20%; the much smaller Be2+ causes a doubling of the atomic charges within the water molecule. The effect of Li+ is found somewhere between these extremes. The negative charge on F- has an opposite effect on the electronegativities. Due to the different geometry of the cluster, F- causes an asymmetric polarization of the water molecules. Again, the values presented in Figure 5 result from a snapshot of the MC simulation, thus causing slight deviations between atoms and molecules which would have been expected to be symmetric. The electronegativities are moderated if a next shell of water molecules is added (Figure 6, Li+-64H20). The central cation forces the molecules of the inner shell to form a slightly distorted tetrahedron. The second shell is less ordered, but still the oxygens are oriented toward the cation, thereby imposing a negative potential contribution on the first-shell molecules which partially compensates the positive potential of the cation. The final result is a decrease of electronegativities with increasing distance from the cation, leaving a gap between the first and the second coordination shell. The latter is illustrated by Figure 7, which shows the number of molecules in certain electronegativity ranges summed up over the equilibrium phase of a simulation. The gap between the electronegativities of the first and the second coordination shell is rather narrow with respect to the broad distribution of the electronegativities of the molecules
'
5.5 1 Figure 4. Eigenvalues of the polarization channels (eigenvectors of the internal charge redistribution) of water clusters of different size and structure. Both modes of the isolated molecule and the lowest and highest modes within both groups belonging th the cyclic trimer (C) are schematicallyrepresented in the inset. Solid black circles and open circles represent opposite directions of the electron flow, respectively; the areas of the circles correspond to the relative amount of the electron flows.
m
la3
23
0439
a
-0.492
-
0
-
g
m
0.245
-0 85
-0.841
074A
&b
0.428
0 B&'
@
0 431
0244
cs+
0.402 -0.855
-0.491
m
4
0 247
0 245
m -0 483 0 30*%0
1244
Li+
304
0305
-0484
0 310 -0606
0
F-
0176
0 179
g o 6 0 3
v
0.308
0 293
-0.609 0.303
12921
0.175
e22
0.483
Figure 5. Water trimers surrounding a Be2+, a Li+, a Cs', and a Fion after reaching an equilibrium state. Atomic charges on each atom are given as well as the molecular electronegativities (underlined).
within each shell. In an ideal case, one would expect a sharp distribution corresponding to the first and second coordination shell. The broadening of the distribution and the lack of a second gap between the second and third shell reflect the fact that liquid water is not an ideally ordered system and that within a shell substructures might be formed which cause a range of electronegativitiesrather than an equalization, so as, for example, the chain trimer described in the section above. Zeolite-Water Clusters. For the last example we turn to the interaction between water and silicate or aluminosilicate clusters. Interactions of this type are important for a variety of reasons. The structural chemistry of silicates and aluminosilicates is of major importance for industrial applications as well
Investigation of Supramolecular Systems
J. Phys. Chem., Vol. 99, No. 38, 1995 13883
7.53
m-
Figure 6. Complete first and parts of the second hydration shell around a Li+ ion, cut off from a LP64H20 cluster. The molecular electronegativities are noted at each water molecule (underlined). For reasons of clarity, the molecules of the first shell are drawn as bonded to the cation. skmdkhell
’
‘
’
third shell
I
a 100 3 c
c
0
50
hard sphere parameters of the model, i.e., 1.9 A for both the HO-H -0Zeo and the ZeoOH- .OH2 distance. The equilibrium structures form in all cases 6-ring structures with atoms of opposite charge adjacent to each other. In the pure silica clusters the water hydrogens are oriented toward the lattice oxygens (Figure Sc,d). In the acidic clusters one OH bond of the water molecules forms a bridge between the Bronsted hydrogen and the neighboring lattice oxygen (Figure Sa,b), while the remaining hydrogen points toward the center of the 12-ring. On the atoms in the ring a 20% increase in charge is observed with respect to the isolated water molecule. The atomic charges on the water molecule remain nearly the same when changing from the 12-ring to the small cluster, but the electronegativity differences between the water and the silicate clusters have changed significantly (Ax in Sa 0.69 eV/ e; in Sb 1.16 eV/e). This is caused by the change in the cluster composition (from a Si/Al ratio of 5 in the 12-ring to a pure aluminum cluster for the simple tetrahedron) and the different shape of the potential energy surface of both clusters. The latter corresponds to the fact that a certain minimum size of a zeolite cluster is required in order to mimic the potential gradient, Le. the electrostatic field, of a crystalline zeolite by a cluster.62For the pure silica clusters, the change of the electronegativity difference is almost negligible (Ax in 8c -0.38 eV/e; in 8d -0.48 eV/e) since the composition differs only by the OWSi ratio. A remarkable effect of the cluster size is also observed for the hardness parameters. This time, the effect is less dependent on the chemical composition than on the size; that is, we are observing similar changes when going from the Al-containing 12-ring to the small cluster (Sa 8b) as well as when changing from the pure silica 12-ring to the corresponding small cluster (8c Sd). In both cases, the hardness of the water molecule as well as the correlation hardness between water and the zeolite cluster is much smaller when interacting with the bigger cluster. Obviously, the large clusters are much more suitable to stabilize changes in an adsorbed molecule than small clusters. The Fukui indices of the water molecules and the corresponding g indices on the cluster show in more detail what is happening. While the Fukui indices of the isolated water molecule are nearly equal for all three atoms ( f ~ ~ f *o0.333), in the adsorbed species the atoms next to the zeolite surface exhibit higher values. The magnitude of this effect is strongly dependent on the size of the zeolite cluster, especially for the acid cluster (Sa and 8b). To stabilize local changes on the water molecule, the corresponding atoms of the cluster, i.e., the nearest neighbors in the “ring”, have to undergo a charge flow in the opposite direction, indicated by a negative g index. Obviously, a large cluster may compensate such charge flows more easily than a small one. This means that the expected reaction, namely, the heterolytic dissociation of the adsorbed water molecule, would be easier on the large cluster than on the small one. Therefore, the choice of sufficiently large zeolite clusters seems to be crucial for estimating some of the most interesting properties of the system, especially the molecular hardness and the local reactivities. This is especially important for zeolites, when the opposite walls of a cage or pore are close to the adsorbed molecule, thereby increasing the effect of stabilization. If the pores are large compared with the adsorbed molecule or on a flat surface, the effect will be less pronounced than in pores of a similar size as that of the molecule.
-
6.8 7 0 7 2 7 4 7 6 7 8 8.0 8 2 8.4
-
Xmol [ e W Figure 7. Frequency count of water molecules according to their electronegativity in a Li+-64H20 cluster within ranges of 0.05 eVfe.
as for fundamental research.57 Through their enormous variability, comparable only to carbon compounds, silicates are especially suited for the application of computational chemistry. We focus here on the interaction of water with zeolite clusters. This special interaction is interesting in at least two respects. On one hand, we obtain in this way a better insight into the origin of catalytic activity and selectivity in zeolite catalysis. On the other hand, the interaction between water, (ahmino-)silicate clusters, and other ingredients such as organic and inorganic cations and anions is relevant for understanding zeolite synthesis. Here, the formation of a supramolecular structure is considered to be extremely i m p ~ r t a n t . ~ * - ~ We restrict ourselves at present to the interaction of a single water molecule with rather small zeolite clusters, namely, a 12ring from the faujasite structure and a simple (H)T(OH)4 tetrahedron (T = Al, Si). The clusters were taken from a DLS61 optimized all-silica faujasite structure, terminated by 0-H groups. After substituting Si by Al, the T-0 distances were adjusted again (from 1.61 to 1.79 A). The hydrogens were placed in the T-0-T plane at the bisector of the T-0-T angle at a distance of 0.96 8, from the oxygen. In Figure 8, the optimized geometries of the water-zeolite complexes are shown. We used the FAU-12-ring as a cluster completely surrounding the water molecule at least in two dimensions. The water-zeolite distances are determined by the
-
Conclusions The combination of the Monte Carlo simulation technique with the EEM formalism allows one in a very consistent way
Toufar et al.
13884 J. Phys. Chem., Vol. 99, No. 38, 1995
0.208 -0.415 0.282 -0.437
', 1'
0.228 I
J'
0.208
0.208
0.359
0.349
0.230'1
0.202
5.87
5.11 16.74
-0.419
16.67
0.305
0.208 0.347
0.210 0.347
k4Ll :f
/-
\' 4.71 7.76
Figure 8. Optimized geometries of a single water molecular adsorbed at different zeolite clusters: (a, top left) 12-ring with 2 A1 on T-positions (light atoms) and two bridging hydroxyls; (b, bottom left) H-Al(OH)4 cluster; (c, top right) all silica 12-ring; (d, bottom right) Si(OH)4 cluster. In parts a and c the electronegativitiesof the molecules are given in roman type in the upper part, together with the molecular hardnesses (in italics) and the correlation hardness (in italics, encircled). In the inset, the charges of the atoms in the water molecule and the adjacent atoms in the zeolite cluster are given in roman type, together with the atomic Fukui indices for the water molecule (in italics) and the corresponding correlation indices for the adjacent atoms in the cluster (in parentheses). In parts b and d the electronegativities of the molecules are kiven in roman type at the left hand side, together with the molecular hardnesses (in italics) and the correlation hardness (in italics, encircled). At the molecules themselves, the charges of the atoms in the water molecule and the adjacent atoms in the zeolite cluster are given in roman type, together with the atomic Fukui indices for the water molecule (in itdlics) and the corresponding correlation indices for the adjacent atoms in the cluster (in parentheses).
to investigate some of the major 'properties of supramolecular systems which are dominated by nondispersive forces. It provides means to (i) optimize the intermolecular structure of supramolecular systems, with adjustment of the charge distribution; (ii) calculate the charge sensitivity at an atomic resolution; and (iii) investigate the behavior of the system with an external or internal perturbation at the molecular as well as the atomic level.
The calculations presented here indicate that the correlation hardnesses of supramolecular systems are maximized if the system achieves an equilibrium situation. Furthermore, within homogeneous systems an equalization of the molecular electronegativities takes place. Important parameters of the molecules in a supramolecular system, such as the atomic charges, the local reactivities, and
Investigation of Supramolecular Systems the molecular hardnesses, are strongly dependent on the intermolecular interaction, i.e., on the kind of interacting molecules as well as on their position with respect to each other. For the case of molecules adsorbed on aluminosilicate clusters, the molecular hardness, the correlation hardness, and the local reactivities are also dependent on the size of the aluminosilicate cluster. One has to keep this in mind when discussing the catalytic activity and selectivity of zeolites or other aluminosilicates on the basis of cluster calculations, since both properties are closely related to the above mentioned parameters. We see this approach as complementary with highly accurate ab initio methods. It lacks some accuracy compared to such methods, but it allows one to determine important electronic properties for large systems which are too complex to be handled by ab initio methods. This is especially important for the charge sensitivities, which are hard to calculate by other methods. Further work is in progress to improve the force field, especially to include dispersive interactions in an appropriate way. This will be done on the basis of the local and global softnesses, which give a measure of the polarizability of the atoms in the system. Acknowledgment. The authors wish to thank Prof. R. F. Nalewajski for intensive and fruitful discussions. H.T. is grateful to the European Community (Human Capital and Mobility Program) for a Postdoctoral Mandate. G.O.A.J. thanks the Institute for Scientific Research in Industry and Agriculture (I.W.O.N.L.) and the Flemish Institute for the Support of Scientific-Technologic Research in Industry (I.W.T.) for post graduate research grants. B.G.B thanks the Research Council of the Catholic University of Leuven (Onderzoeksraad) for a Postdoctoral Mandate (PDM 93/57). This research was made possible with a GOA grant of the Flemish Government. References and Notes (1) Jorgensen, W. L. J . Chem. Phys. 1982, 77 (7), 4156-4163. (2) Ahlstrom, P.; Wallquist, A.; Engstrom, S.; Jonnson, B. Mol. Phys. 1989,68, 563-581. (3) Zhu, S.-B.; Yao, S.; Zhu, J.-B.; Singh, S.; Robinson, G. W. J . Phys. Chem. 1991, 95, 621 1-6217. (4) Leherte, L.; Andre, J.-M.; Vercauteren, D. P.; Derouane, E. G. J . Mol. Catal. 1989, 54, 426-438. ( 5 ) Leherte, L.; Andre, J.-M.; Vercauteren, D. P.; Derouane, E. G. lnt. J . Quantum Chem. 1992,42, 1291-1326. (6) Jameson, C. J.; Jameson, A. K.; Baello, B. I.; Lim, H.-M. J . Chem. Phys. 1994,100 (8), 5965-5976. (7) Smit, B.; Siepmann, J. I. Science 1994, 264, 1118-1 120. (8) Lee, S. H.; Moon, G. K.; Choi, S. G.;Kim, H. S. J. Phys. Chem. 1994,98, 1561-1569. (9) Feuston, B. P.; Garofalini, S. H. Chem. Phys. Lett. 1990,170, 264270. (10) Garofalini, S. H.; Martin, G. J . Phys. Chem. 1994, 98, 1311-1316. (11) Durell, S. R.; Brooks, B. R.; Ben-Naim, A. J . Phys. Chem. 1994, 98, 2198-2202. (12) Sanderson, R. T. J . Am. Chem. Sac. 1952, 74, 272-274. (13) Politzer, P.; Weinstein, H. J . Chem. Phys. 1979, 71, 4218-4220. (14) Rappe, A. K.; Goddard, W. A., III, J. Phys. Chem. 1991,95, 33583363. (15) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.k 111; Skiff, W. M. J . Am. Chem. Sac. 1992,114, 10024-10035. (16) Casewit, C. J.; Colwell, K. S.; Rappe, A. K. J . Am. Chem. Sac. 1992, 114, 10046-10053. (17) van Duin, A. C. T.; Baas, J. M. A.; van den Graaf, B. J . Chem. Sac., Faraday Trans. 1994, 90 (19), 2881-2895. (18) Rick, S.; Stuart, S. J.; Beme, B. J. J . Chem. Phys. 1994, 101 (7), 6141-6156. (19) Baekelandt, B. G.; Mortier, W. J.; Lievens, J. L.; Schoonheydt, R. A. J . Am. Chem. Sac. 1991, 113, 6730-6734.
J. Phys. Chem., Vol. 99,No. 38, 1995 13885 (20) Mortier, W. J.; Ghosh, S. K.; Shankar, S. J . Am. Chem. Sac. 1986, 108, 4315-4320. (21) Mortier. W. J. Struct. Bondina 1987, 66, 125-143. (22) Mortier, W. J.; Van Genechtei, K.; Gasteiger, J. J . Am. Chem. Sac. 1985,107, 829-835.
(23) Mortier, W. J. In Theoretical Aspects of Heterogeneous Catalysis; Van Nostrand Reinhold Catalysis Series: New York, 1990; pp 135159. (24) Mortier, W. J. Stud. Sur$ Sci. Catal. 1988, 37, 253-268. (25) Sanderson, R. T. Science 1955, 121, 207-208. (26) Parr, R. G.; Donnelly, R. A.; Levy, M.; Palke, W. E. J . Chem. Phys. 1978, 68 (8), 3801-3807. (27) Hohenberg, P.; Kohn, W. Phys. Rev. B 1964, 136, 864-871. (28) Kohn, W.; Sham, L. J. Phys. R e v . A 1965, 140, 1133-1138. (29) Gilbert, T. L. Phys. Rev. B 1975, 12, 2111-2120. (30) Palke, W. E. J . Chem. Phys. 1980, 72, 2511-2514. (31) Parr, R. G. Ann. Rev. Phys. Chem. 1983, 34, 631-656 and references therein. (32) Parr, R. G.; Yang, W. In Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (33) Ray, N. K.; Samuels, K.; Parr, R. G. J . Chem. Phys. 1979, 70, 3680-3684. (34) Guse, M. P. J. Chem. Phys. 1981, 75, 828-833. (35) Baekelandt, B. G. Ph.D. Thesis, Leuven, 1992. (36) Baekelandt, B. G.; Mortier, W. J.; Schoonheydt, R. A. Struct. Bonding 1993, 80, 187-227. (37) Nalewajski, R. F. Struct. Bonding 1993, 80, 115-186. (38) Van Genechten, K. A.; Mortier, W. J.; Geerlings, P. J . Chem. Phys. 1987,86, 5063-5071. (39) Van Genechten, K. A,; Mortier, W. J. Zeolites 1988, 8, 273283. (40) Uytterhoeven, L.; Lievens, J.; Van Genechten, K.; Mortier, W. J.; Geerlings, P. In Adsolption in Microporous Adsorbents; preprints of the workshop in Eberswalde; 1987; Vol. 1, pp 73-80. (41) Uytterhoeven, L.; Mortier, W. J. J . Chem. Sac., Faraday Trans. 1992, 88, 2747-2751. (42) Nalewajski, R. F.; Korchowiec, I. J . Mol. Catal. 1993, 82, 383405. (43) Nalewaiski. R. F.: Koster. A. M.: Bredow. T.: Jug. K. J . Mol. Catal. 1963,82, 407-223, (441 Kollman. P.A.; Allen, L. C. J . Am. Chem. Sac. 1971.93, 49915000. (45) Kutenmacher, H.; Popkie, H.; Clementi, E. J. Chem. Phys. 1973, 58, 1689-1699. (46) Korchowiec, J.; Bredow, T.; Jug, K. Chem. Phys. Lett. 1994,222, 58-64. (47) Nalewajski, R. F.; Korchowiec, J.; Michalak, A. Proceedings Chemical Science; Gadre, S . , Ed., in press. (48) Nalewajski,R. F.; Korchowiec,J.; Zhou, Z . lnt. J. Quantum Chem.: Quantum Chem. Symp. 1988,22, 349-366. (49) Nalewajski, R. F.; Koninski, M. 2.Naturforsch. 1987,42a, 451462. (50) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 1087-1092. (51) HyperChem, Release 3 for Windows, Molecular Modeling System; Hypercube, Inc., and Autodesk, Inc.: Waterloo, Ontario, 1993. (52) Pearson, R. G. J . Chem. Educ. 1987, 64, 561-567. (53) Parr, R. G.; Chattaraj, P. K. J . Am. Chem. Sac. 1991,113, 18541855. (54) Pearson, R. G.; Palke, W. E. J . Phys. Chem. 1992,96, 3283-3285. ( 5 5 ) Hermansson, K.; Ojamae, L. lnt. J . Quantum Chem. 1992, 42, 1251- 1270. (56) (a) Janssens, G. 0. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A. lnt. J . Quantum Chem., in press. (b) Baekelandt, B. G.;Janssens, G. 0. A.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A.; Nalewajski, R. F. J . Phys. Chem. 1995, 99, 9784-9794. (57) Liebau, F. Structural Chemistry of Silicates; Springer-Verlag: Berlin, 1985. (58) Gabelica, Z.; Blom, N.; Derouane,E. G. Appl. Card. 1983,5, 227248. (59) McCormick, A. V.; Bell, A. T. Catal. Rev. Sci. Eng. 1989, 31, 97- 127. (60) k " r , G. 0. Zeolites 1992, 12, 428-430. (61) Baerlocher, C.; Hepp, A.; Meier, W. M. DLS-76:A Program for the Simulation of Crystal Structures, ETH: Zurich, Switzerland, 1978. (62) Janssens, G. 0. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A. J . Phys. Chem. 1995, 99, 3251-3258. Y
Jp950951C
