J. Phys. Chem. B 2008, 112, 5153-5162
5153
Molecular Dynamics Simulations of Dendrimer-Encapsulated r-Keggin Ions in Trichloromethane Solution Ralf Brodbeck, Thorsten To1 nsing, and Dirk Andrae* Theoretische Chemie, Fakulta¨t fu¨r Chemie, UniVersita¨t Bielefeld, Postfach 10 01 31, D-33501 Bielefeld, Germany
Dirk Volkmer Institut fu¨r Anorganische Chemie II, UniVersita¨t Ulm, Albert-Einstein-Allee 11, D-89081 Ulm, Germany ReceiVed: October 22, 2007; In Final Form: January 28, 2008
We report on molecular dynamics simulations of dendrimer-encapsulated R-Keggin ions in trichloromethane solution. The simulations were done within the NVE ensemble at temperatures around T ) 300 K. The eight examined systems are model compounds for dendrizymes, a hybrid material where a polyoxometalate ion (the core) is surrounded by amphiphilic cationic dendrimers (the shell) such that the complete system may exhibit enzyme-like regioselectivity and substrate selectivity, e.g., in olefin oxidation. The influence of dendrimer type, dendrimer generation, and number of dendritic cations bound by electrostatic interaction to the polyoxometalate core on the structure and dynamics of the shell has been studied. It is shown that the resulting distribution of trichloromethane molecules within the shell may serve as an indicator for the shell’s permeability for small molecules. The dendritic shell causes a size exclusion effect that influences the access of small molecules to the central polyoxometalate ion, i.e., to that part where the enzyme-like reaction of a dendrizyme is supposed to take place.
1. Introduction During the last 30 years, several polyoxometalates (POMs) of molybdenum and tungsten were found to be potent regioselective oxidation catalysts for a variety of homogeneous catalytic reactions (see reviews1-5 and references therein). Since catalytic transformations employing POM catalysts often yield product mixtures, substrate selectivity became a highly desired additional property for POMs. However, the known catalytically active POMs have no sterically demanding side groups or other structural elements that could help to preferentially recognize a desired target substrate. Because of the chemical inertness of the M(VI)-O bonds (M ) Mo or W), the covalent attachment of organic ligands to the metal-oxygen framework most often leads to product mixtures that display completely rearranged cluster cores with concomitant loss of their original catalytic activity.6 To circumvent the difficulties of attaching substrate recognizing groups directly to the metal-oxygen cluster via chemical modification, Volkmer and co-workers7 combined a catalytically active POM with amphiphilic cationic dendrimers. The resulting ion clusters, members of a class of compounds known as “dendrizymes”,7,8 are composed of a central POM, the “dendrizyme core”, and surrounding cationic dendrimers, the “dendrizyme shell”. Such an ion cluster indeed shows regioselectivity and substrate selectivity in olefin epoxidation in trichloromethane solution, with hydrogen peroxide as the oxidizing agent, but slow degradation occurs under reaction conditions.9 Both experimental and theoretical methods were applied to characterize the first dendrizyme, to elucidate its structure and to understand its substrate selectivity on the atomic scale.7 * Corresponding author. E-mail:
[email protected].
Molecular dynamics (MD) simulations using the Amber94 force field were performed for a dendrizyme model compound in the gas phase, i.e., without solvent molecules. The structure of the POM was kept frozen, and Mulliken partial charges were used to parametrize the nonbonding interactions between POM and the dendrimers. The present work reports on results obtained from MD simulations under periodic boundary conditions for eight dendrizyme model systems in trichloromethane solution. The internal degrees of freedom of all molecules in the simulation cell were active. A force field for the POM, to be used within the generalized Amber force field (GAFF) approach, was newly parametrized from quantum chemical reference data. Timeaveraged structures and radial and spatial distribution functions, obtained from a detailed analysis of the MD trajectories are presented and discussed below. The dendrizyme core of the model compounds studied here was chosen from the R-Keggin ions [(XO4)W12O36]n- (1: X ) P, n ) 3; 2: X ) Zn, n ) 6; see Figure 1). These heteropolyoxotungstates can be considered as the smallest possible POMs to represent the core of a dendrizyme model compound. Both Keggin ions, being of nearly spherical shape, are also quite comparable in their radial size and volume. On the other hand, their size still permits application of sophisticated quantum chemical methods, such as first-principles or ab initio techniques, to provide data for the derivation of force field parameters. Because of their nearly spherical shape, the central heteroatom X (P or Zn) of the Keggin ions provides a natural origin for the radial and spatial distribution functions obtained from the MD simulations. For the dendrizyme shell, two types of dendrimers of the first and second generation were chosen. One type is derived from
10.1021/jp710215u CCC: $40.75 © 2008 American Chemical Society Published on Web 04/01/2008
5154 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Brodbeck et al. van der Waals and electrostatic contributions, accounting for interactions between nonbonded atoms:
UFF )
∑
Kr(r - req)2 +
bonds
{1 + cos(nθ - θ0)} +
Figure 1. Ball-and-stick representation of the R-Keggin ions [(XO4)W12O36]n- (point group Td, 1: X ) P, n ) 3; 2: X ) Zn, n ) 6). Labels on symmetry unique atoms indicate the atom types defined in our force field.
an N-methyl pyridinium ion (“pyridinium dendrimers”, 3a and 3b in Figure 2), the other from a quaternary dimethyl ammonium ion (“ammonium dendrimers”, 4a and 4b in Figure 2). Both types differ mainly in the conformational flexibility of their dendritic molecular chains. In the former dendrimers, only hindered rotations of aryl rings are possible, whereas the latter dendrimers have a comparably larger flexibility due to the methylenoxy links between their aryl rings. From the Keggin ions and dendrimers, a total of eight different non-charged ion clusters can be built up (see Table 1 for their sum formulas, and corresponding shorthand labels ngkx are used in the following, with association number n ) 3 or 6, dendrimer generation k ) 1 or 2, and h or w for dendrimer type x). For illustration, structures of these compounds (time-averaged over the trajectory obtained from the MD simulation) are shown in Figures 3 and 4. The influence of various parameters on the structure and the dynamics of the dendrizyme shell can be studied within this set of eight compounds: (i) the influence of the association number, induced by a change of the central Keggin ion (sequences 3gkx vs 6gkx); (ii) the influence of conformational flexibility, due to variation on the dendrimer type (sequences ngkh vs ngkw); and finally (iii) the influence of the dendrimer generation (sequences ng1x vs ng2x). This paper is structured as follows. Computational details are briefly discussed in the next section. Then, in Section 3, the results from the MD simulations of dendrizyme model compounds in trichloromethane solution are presented, with emphasis on the structure of the dendrizyme shell and its permeability for small molecules. A summary with conclusive remarks follows as Section 4. Supporting Information, related to the generation of the POM force field and to the MD simulations, is provided in electronic form. 2. Computational Details 2.1. The Force Field. A series of classical MD simulations was performed for the eight dendrizyme model systems in trichloromethane solution, using the Amber suite of programs.10 A prerequisite for this approach is the availability of an empirical force field determining the potential energy UFF of the full system under study (POM, dendrimers, and solvent molecules). Within the GAFF,11 this potential energy is represented as a sum over harmonic contributions, to represent deformation of bonds, angles, and dihedral coordinates for bonded atoms, and
∑
Kφ(φ - φeq)2 +
2 qiqj
{ [( ) ( ) ] }
angles iej
∑
∑
Vn
ij
σij
κ
Rij
12
-2
dihedrals σij 6
+
Rij
κ′Rij (1)
Here, req and φeq are equilibrium structural parameters, Kr, Kφ, and Vn are force constants, n is a multiplicity, and θ0 is a phase angle. van der Waals interactions are represented by a 12-6 Lennard-Jones potential, with parameters for different atom types (ij, σij, i * j) being constructed from parameters ii and σii according to the Lorentz-Berthelot rules. For 1-4 interactions, scaling factors were applied as recommended by Kollman and co-workers:12 κ ) 2.0 for van der Waals interactions, and κ′ ) 1.2 for electrostatic interactions (κ ) κ′ ) 1.0 otherwise). Nonbonded interactions were calculated with the particle mesh Ewald summation technique,13 using an atom-based cutoff radius Rc ) 1.5 nm. Orthorhombic simulation cells and periodic boundary conditions were used in all simulations (see Table 1 for cell dimensions). In order to properly describe the thermal energy transfer between core and shell of the dendrizyme model system, the Keggin ion’s internal degrees of freedom should be active during the MD simulation. Therefore, a force field parametrization of the two different Keggin ions was necessary.14 We added four new atom types to the GAFF: “w” (tungsten), “o1” (Oµ1 oxygen, ∠(W-O-W) ≈ 125°), “o2” (Oµ2 oxygen, ∠(W-OW) ≈ 150°), and “ot” (terminal oxygen). Force field parameters for all other centers of the Keggin ions (phosphorus, zinc, phosphate, and zincate oxygens) were already available within the GAFF (Lennard-Jones parameters for phosphate oxygens were also used for the zincate oxygens, and bonding parameters for the zincate oxygens were explicitly given). The LennardJones parameters ii and σii for the four new atom types are given in Table 2. Parameters for the oxygen centers were taken from those for chemically analogous organic oxygen types (ether, ketone), and parameters for tungsten were obtained in the following way: the parameter ii was calculated with semiempirical rules from Halgren,15 and the parameter σii was taken as the average of density isovalue surface radii (F ) 0.002 a.u.) obtained from atomic MCSCF calculations for low-lying states with electron configurations 5d46s2, 5d56s1, and 5d6 (see Supporting Information for further details). The parameters for bonding interactions in the Keggin ion are given in Table 2. All parameters required for bonding interactions and LennardJones terms in the dendrimers were available in the GAFF. As can be seen from eq 1, complete definition of the potential energy also requires a set of partial charges. These partial charges were obtained from the molecular electrostatic potential (MEP) with the RESP method16,17 (which has several advantages over other methods used for this purpose). The MEP of the Keggin and dendrimer ions was calculated from the electron densities at the corresponding equilibrium structures (Tdsymmetric for the Keggin ions, Cs-symmetric for the pyridinium ions, and C2V-symmetric for the ammonium ions). The equilibrium structures were obtained from Kohn-Sham density functional calculations18 using the B3LYP19,20 gradient-corrected hybrid density functional together with effective core potentials (ECPs) and associated basis sets of Stevens and co-workers21,22 (no ECP for hydrogen). The basis sets for the Keggin ions were
MD Study of Dendrimer-Encapsulated R-Keggin Ions
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5155
Figure 2. Structural formulas of first- and second-generation dendrimers derived from the N-methyl-3,5-bis(tert-butyl)pyridinium ion (3a, 3b) and from the bis(3,5-dimethoxybenzyl)dimethylammonium ion (4a, 4b).
TABLE 1: Parameters of MD Simulations (NVE Ensemble) for the Eight Dendrizyme Model Systems in Trichloromethane Solutiona dendrizyme model system and shorthand label
Nsol
cell dimensions nm
F g/cm3
δt fs
∆ttraj ns
ttot ns
Tav K
σT K
(C34H48N)3[PW12O40] (C74H96N)3[PW12O40] (C52H60NO12)3[PW12O40] (C116H124NO28)3[PW12O40]
3g1h 3g2h 3g1w 3g2w
443 601 1107 1357
3.909, 3.914, 4.129 4.207, 4.510, 4.578 5.363, 5.368, 5.365 5.867, 5.752, 5.744
1.51 1.48 1.48 1.46
1.8 2.0 1.8 2.0
0.80 0.80 0.80 0.80
2.83 2.53 2.92 2.56
301.79 300.07 301.05 300.38
3.36 2.81 2.21 1.93
(C34H48N)6[ZnW12O40] (C74H96N)6[ZW12O40] (C52H60NO12)6[ZnW12O40] (C116H124NO28)6[ZnW12O40]
6g1h 6g2h 6g1w 6g2w
536 912 1120 1217
4.226, 4.422, 4.200 5.069, 5.184, 5.122 5.554, 5.232, 5.507 5.411, 5.867, 5.746
1.47 1.45 1.47 1.46
1.8 1.8 2.0 2.0
0.72 0.72 0.80 0.80
2.10 2.62 2.70 2.15
301.05 301.80 300.41 300.21
2.97 2.84 2.04 1.94
a
Parameters: number of solvent molecules (Nsol) in the orthorhombic simulation cell, cell dimensions, equilibrated solution density (F), used time step (δt), trajectory length (∆ttraj), total simulation time (ttot), equilibrated average temperature (Tav), and its standard deviation (σT).
augmented by polarization and diffuse functions on all centers (this basis set is termed SBKJC+(d,f)). In the case of the dendrimer ions, polarization functions on all non-hydrogen centers were added (this basis set is termed SBKJC(d)). Equilibrium structures, corresponding total electronic energies, and the RESP partial charges for the Keggin and dendrimer ions are given as Supporting Information, as are technical details of the Kohn-Sham density functional calculations. 2.2. Molecular Dynamics Simulations. An important first decision, required for any MD simulation, is the choice of the start configuration, i.e., the initial positions of all atoms (nuclei). For the dendrizyme model systems studied here, isomeric forms, differing only by the mode of electrostatic association of the dendrimers to the central POM, are separated by high-energy barriers (insurmountable even at simulation temperatures around T ) 1000 K). The thermodynamically most important isomers of the ion clusters should be those lowest in energy, where
attractive electrostatic interactions are strong. Therefore, the dendritic cations have been placed around the POM anions in such a way that regions with highest magnitude of the electrostatic potential are close together, at a distance roughly corresponding to the sum of their van der Waals radii. The dendrizyme model compounds with these start conformations were then embedded in an orthorhombic cell, filled with trichloromethane molecules according to a preoptimized liquid structure (available within Amber). Trichloromethane molecules in the space occupied by the dendrizyme were removed. In the next step, a simultaneous relaxation of total energy and density was performed through simulation within the NPT ensemble. Temperature and density equilibration took place under tight coupling to a Berendsen thermostat and barostat23 (with initial coupling constant τ ) 0.1 ps, and later τ ) 0.5 ps), using an isothermal compressibility of κ ) 1 × 10-5 bar-1. Subsequent temperature and density equilibration under weak
5156 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Brodbeck et al.
Figure 4. Time-averaged structures of the dendrizyme model compounds with ammonium dendrimers (hydrogen atoms omitted for clarity). Figure 3. Time-averaged structures of the dendrizyme model compounds with pyridinium dendrimers (hydrogen atoms omitted for clarity).
coupling (τ ) 2.0 ps) led to the target average pressure of P ) 1 bar. The coordinates and velocities thus obtained were used
to start an equilibration within the NVT ensemble with weak thermostat coupling (τ ) 2.0 ps), keeping the cell dimensions fixed. This was run until the average temperature of T ) 300 K was reached. The resulting coordinates and velocities were used to start the production simulations within the NVE ensemble. The integration of the equations of motion was done
MD Study of Dendrimer-Encapsulated R-Keggin Ions TABLE 2: Force Field Parameters for the r-Keggin Ions [(XO4)W12O36]n- (1: X ) P, n ) 3; 2: X ) Zn, n ) 6) to Be Used with the GAFF11 (See Eq 1) involved atom types
bonding parameters
i-j
Kr kcal mol-1 (pm)-2
req pm
zn-o w-o w-ot w-o1 w-o2
200.0 200.0 620.0 300.0 300.0
190.4 247.0 171.0 193.0 192.0
i-j-k
Kφ kcal mol-1 degree-2
φeq degree
o-zn-o zn-o-w o-w-o2 o-w-ot o-w-o1 p-o-w w-o-w w-o1-w w-o2-w ot-w-o1 ot-w-o2 o1-w-o2 o1-w-o1 o2-w-o2
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
109.47 119.4 83.4 170.8 71.4 125.8 89.3 127.9 152.6 102.2 103.3 125.0 115.0 130.0
i-j-k-l
V2 kcal mol-1
θ0 degree
w-o1-w-o2
0.61
191.3
atom type
Lennard-Jones parameters
i
ii kcal mol-1
σii pm
w o1 o2 ot
0.54 0.17 0.17 0.21
222 168 168 166
with the velocity Verlet algorithm, as implemented in the Sander program of the Amber suite. No constraints (i.e., the Shake algorithm) or restraints were applied. Coordinates were saved every 20 time steps to form a trajectory (in the NVE ensemble) of 0.72 ns or 0.80 ns length. On average, this is time enough for a CH3Cl molecule to diffuse a distance of roughly 2.3 nm, the same order of magnitude as the box dimensions (see Table 4 in the Supporting Information for translation diffusion coefficients of trichloromethane). The total simulation time (in all three ensembles) ranges from 2.10 ns to 2.92 ns. The values of the time steps and other important technical parameters of the MD simulations are given in Table 1. Averaged temperatures (averaged over the complete NVE trajectory) and their standard deviations are also included. We checked for all NVE simulations the constancy of total energy by monitoring it as a function of time. Figure 1 in the Supporting Information shows a plot of the total energy versus the simulation time of the saved trajectories. No drifts or discontinuities are visible, i.e., on the examined time scale, the thermal equilibrium was kept in all NVE simulations. 2.3. Trajectory Analysis. The trajectories obtained from the MD simulations were analyzed with respect to the structure of the dendrizyme shell and the distribution of trichloromethane within the shell. For this purpose, time-averaged dendrizyme structures and radial and spatial distribution functions were used. The time-averaged structures were obtained from averaging the position vectors of every dendrizyme atom over the trajectory.
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5157 Radial distribution functions gXA(r) give the time-averaged frequency to find a center of type A in a spherical shell of thickness ∆r at distance r from a center of type X. However, the full spatially resolved information can be obtained only from spatial distribution functions gAXB-D(r). These provide the frequency to find a center D in a (cubic) volume element ∆r at position r relative to a center of type X. The reference coordinate system for r is defined by the plane through centers A, X, and B (with X as the apex), together with a plane’s normal. The radial distribution functions gXA(r) yield, primarily, information about the angularly averaged distance between the centers X and A. They are also convenient to analyze the influence of dendrimer type, dendrimer generation, and number of coordinated dendrimers on the radial extension of the dendrizyme shell. With X ) P or X ) Zn, and A taken as an atom of trichloromethane, these functions are particularly well suited to study the average distribution of trichloromethane within the dendrizyme shell, as developed during the simulation (trichloromethane was not initially present in the space occupied by the dendrizyme model compounds). Even though spatial distribution functions provide much more detailed information, we will use mainly the radial distribution functions to discuss the distribution of trichloromethane within the dendrizyme shell, and its permeability for substrate molecules. The calculation of radial distribution functions was done with the Ptraj program, which is part of the Amber suite, using a shell thickness of ∆r ) 20 pm. Therefrom follows an inherent uncertainty of all radii r extracted from our analysis, such that all radii are to be understood as r ( (1/2)∆r ) r ( 10 pm. In order to make the comparison between radial distribution functions gXA(r) for different dendrizyme model compounds an max easier task, rescaled functions g˜XA(r) ) gXA(r)/gmax XA , with gXA ) max gXA(r) and X ) P or X ) Zn, are frequently presented in the following (except for cases where center A belongs to trichloromethane). Spatial distribution functions were calculated with a modified version of the Sidan program.24,25 For the purpose of visualization, images of spatial distribution functions were generated with the Molekel program26,27 from data converted to the Gaussian cube format. 3. Results and Discussion All data for the eight dendrizyme model compounds in trichloromethane solution, being presented and discussed in the following, were obtained from MD simulations within the NVE ensemble (at temperatures around T ) 300 K). 3.1. Structure of the Dendrizyme Shell. In this subsection, the structure of the dendrizyme shell is discussed on the basis of time-averaged structures and radial distribution functions. As already mentioned in Section 2.2, different coordination isomers are separated by high-energy barriers. Therefore, energetically most favorable start configurations were chosen, such that parts of the ions with opposite high magnitude values of the MEP are in close contact to each other, at a distance roughly corresponding to the sum of their van der Waals radii. Figures 3 and 4 show time-averaged structures for the eight compounds studied here, as obtained from the trajectories. All compounds are stable under the simulation conditions (temperature, pressure, solvent polarity, ionic strength), i.e., the ion clusters neither collapse nor do they dissociate into separate parts. The radial size of these compounds reflects the different steric demand of the various shell-forming dendrimers, since the R-Keggin ions 1 and 2 are almost of identical spherical size. As a measure of radial size, we define the radius R as R ) maxr{gXA(r) ) 0.1} (X ) P, Zn), where A may be any type of
5158 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Brodbeck et al.
TABLE 3: Radial Size R of the Dendrizyme Model Systemsa system R/nm
3g1h 1.43
3g2h 2.05
3g1w 1.97
3g2w 2.29
6g1h 1.73
6g2h 2.01
6g1w 2.01
6g2w 2.45
a
See text for details of definition; system labels as given in Table 1.
center present in the dendrimers. The resulting radii, ranging from R ≈ 1.4 nm to R ≈ 2.5 nm, are given in Table 3. As expected, the outermost centers of the dendrizyme shells are the hydrogen atoms of the terminal groups (methoxy or tertbutyl groups, respectively). Interestingly, the radial sizes are quite similar for the four compounds ng2h and ng1w (R ≈ 2 nm), even though their shell structures turned out to be quite different (see below). The compound 3g1h has the smallest radial size, whereas the largest value was found for compound 6g2w. The rescaled radial distribution functions g˜ XN(r) (X ) P or X ) Zn, and N from the dendrimers), shown in Figure 5, indicate an influence of dendrimer type, dendrimer generation, and number of coordinated dendrimers on the structure of the innermost part of the dendrizyme shell. The nitrogen atoms are found, on average, in a radial shell around X with radii from r ) 560 pm to r ) 860 pm. For the dendrizyme model compounds with the pyridinium dendrimers, ngkh, the first maximum appears at smaller radii (600 pm to 660 pm) than for the compounds with ammonium dendrimers, ngkw (680 pm to 720 pm). The two methyl groups on the nitrogen center seem to hinder a closer contact between the ammonium dendrimers and the Keggin ions. Also, in the case of the pyridinium dendrimers, some pyridinium rings are always arranged tangentially to the nearly spherical central Keggin ion, thus leading to smaller X-N distances for these compounds. In four cases (6g2h, 6g1w, 3g2w, 6g2w), the radial distribution functions are bimodal (show two maxima), whereas a unimodal distribution is found for 3g1h and 6g1w, and even a trimodal one is found for 6g1h and 3g2h. The spatial arrangement of the nitrogen centers around X, as derived from the time-averaged structures, can help to partly understand these findings. All compounds with three dendrimers show a nearly trigonal-planar arrangement, whereas distorted trigonal prisms are found for 6g1h, 6g1w, and 6g2w, and a distorted pentagonal pyramid is found for 6g2h. No constraints were imposed to force the X-N distances in a specific dendrizyme to be identical, and the spatial arrangement found for the time-averaged structures still resemble those in a compound’s start configuration. Hence, the N-methylpyridinium and the ammonium groups behave like anchors for the dendritic molecular chains, which move freely during the MD simulations. Turning now to a discussion of other parts of the dendrizyme shell, we focus on the aryl carbon centers. Two issues are of primary interest: (i) the orientation of the innermost aryl (phenyl and pyridyl) rings relative to the nearly spherically shaped central Keggin ion, and (ii) the radial distribution of the aryl carbon centers. Information on the ring orientation can be obtained again from the time-averaged structures (see Figures 3 and 4). In the following, “tangential aryl ring orientation” means that the normal vector of the ring plane is oriented in parallel to the radius vector from the Keggin ion center X to the center of the aryl ring. In a “radial aryl ring orientation”, the normal vector of the ring and the radius vector from X to the center of the ring are orthogonal. In the compounds 3g1h and 3g2h, with three pyridinium dendrimers, a tangential orientation of aryl rings is frequently found, or even predominates. As the association number is increased (compounds 6g1h and 6g2h),
Figure 5. Rescaled radial distribution functions g˜ XOPOM(r) and g˜ XN(r) (X ) P, Zn) of the Keggin ion’s oxygen centers (XO4 oxygens not shown) and the dendrimer’s nitrogen centers for all examined dendrizyme model systems.
a larger fraction of aryl rings adopts a radial orientation. A similar situation is found for the compounds with the ammonium dendrimers. In the compounds 3g1w and 3g2w, still some rings are in tangential orientation relative to the central Keggin ion, but several rings tend toward a radial orientation, due to the greater conformational flexibility of the molecular chains in these dendrimers. The overcrowding of the dendrizyme shell, accompanied by an increase of the association number (compounds 6g1w and 6g2w), necessarily increases the fraction of aryl rings with radial orientation even further. In compounds with six dendrimers (X ) Zn), close pairs of aryl rings can be found in the dendritic shell, indicating attractive van der Waals interactions between these rings. In all dendrizyme model compounds considered in this work, at least half of the terminal substituted phenyl groups are directed into the bulk solvent phase (the other terminal groups are folded back and interact with the central Keggin ion or inner parts of neighboring dendrimers). Hence, an attractive interaction between dendrizyme shell and solvent seems to be important, leading to solvation of the terminal parts of the dendritic molecular chains. This relates favorably with the experimentally observed solubility of prototype dendrizyme compounds in trichloromethane and other organic solvents of comparable polarity. On the other hand, the importance of an atomistic solvent representation in liquid-phase MD simulations is emphasized. MD simulations in the gas phase (vacuum) gave significantly different results:7 nearly all terminal groups of the dendrimers folded back toward the POM surface. In the vacuum, such an arrangement is energetically favorable because the electrostatic energy is decreased, and competitive interactions between dendrimer and solvent molecules are absent. Consequently, a porous structure, with voids in the dendrizyme shell, was observed. The radial distribution functions of carbon centers in phenyl rings, gXCph(r), start off at distances between 600 pm and 700 pm, as shown by the rescaled radial distribution functions in Figure 6. The global maxima are found between 750 pm and 1030 pm. The higher the amount of tangentially oriented phenyl groups in the innermost part of the dendrizyme shell, the smaller is the radius for the first maximum of gXCph(r). For the compounds 3g1h, 3g1w, and 3g2w, a unimodal distribution is found, whereas a bimodal distribution is found for all compounds with X ) Zn, i.e., the sequence 6gkx (the compound 3g2h is a borderline case). Similar observations can be made
MD Study of Dendrimer-Encapsulated R-Keggin Ions
Figure 6. Rescaled radial distribution functions g˜ XCph(r) of the phenyl carbon centers for all eight examined dendrizyme model systems.
Figure 7. Rescaled radial distribution functions g˜ XCpy(r) of the pyridinium carbon centers for the four dendrizyme model systems ngkh.
in the case of the carbon centers in pyridine rings, and their rescaled radial distribution functions are presented in Figure 7. The bimodality of the functions gXCpy(r) is caused by different orientation of the pyridine rings (tangential or radial), the only exception being 3g1h, which shows a unimodal distribution due to the exclusively tangential orientation of all pyridine rings in this case (this leads also to a unimodal distribution for the phenyl carbon centers, see Figure 6). The phenyl groups directed into the bulk solvent are responsible for the slow decay in the long-range part of the functions gXCph(r), observed for all dendrizyme model compounds except 3g1h (see Figure 6). As could have been expected, the greatest extent of this function is found for compounds with dendrimers of second generation. For a given dendrimer generation, the longer and more flexible dendritic molecular chains of the ammonium dendrimers extend farther into the bulk solvent than those of the pyridinium dendrimers. 3.2. Permeability of the Dendrizyme Shell. A theoretical means to quantify the substrate selectivity of a dendrizyme would be very valuable as a complementary source of information with respect to experimental data. The ease of access of substrate molecules to the dendrizyme core is expected to be greater for smaller molecules, whereas larger molecules are expected to diffuse at a smaller rate (or even not at all) through the dendrizyme shell. This kinetically controlled size exclusion effect may serve to partly explain substrate selectivity of a dendrizyme. Other effects, such as (differential) thermodynamical effects, e.g., free enthalpies of absorption of different
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5159
Figure 8. Radial distribution functions gXHCHCl3(r) of the trichloromethane hydrogen centers for all eight examined dendrizyme model systems.
substrates to the dendrizyme shell, or active transport of substrate molecules by the dendrimer branches through the shell, are being neglected. A quantity or function measuring the permeability of the dendrizyme shell may indicate the lesser or greater ease of access of smaller molecules to the dendrizyme core. We suggest to use radial distribution functions gXA(r) for this purpose, where, as before, X is the central heteroatom of the Keggin ion (X ) P or X ) Zn), but A is now an atom of the solvent molecules (we chose the trichloromethane hydrogen atom, A ) HCHCl3). The nearly spherical shape of the Keggin ions supports the use of radial distribution functions even further, since it permits angular integration of spatial distribution functions without too much loss of relevant information. Therefore, a comparative study of the radial distribution functions gXHCHCl3(r) (one for each of the eight systems under study) over the radial range of the dendrizyme shell is expected to provide some useful information for the problem of estimating the shell’s permeability. The radial distribution functions gXHCHCl3(r) for the eight systems under study are shown in Figure 8. We remind the reader that values of gXA(r) larger (smaller) than the long-range limit indicate a higher (lower) frequency to find A at that radial distance r from X than in the bulk solvent. All radial distribution functions show a more or less pronounced structure over a radial range from 450 pm to about 850 pm: a single maximum, a single maximum with a shoulder, or two maxima can be observed. The width of this range fits well to the expected size of a sphere containing a trichloromethane molecule (internuclear distance dC-Cl ≈ 180 pm). Molecules of trichloromethane are thus in all cases capable of approaching the central Keggin ion, despite the presence of the dendritic shell. In almost all cases, the maximum value of gXHCHCl3(r) in this short-range region is larger than 1, indicating solvation, partially at least, of the Keggin ion by solvent molecules, the only exceptions being the two systems 6gkw where the maximum values in this radial range are smaller than 1. Beyond the minimum around 850 pm, i.e., outside the sphere containing the closest lying solvent molecules, all radial distribution functions increase again and eventually approach the bulk value, the limiting long-range value. Depending on the system, this limiting value is practically approached at radii ranging from 1400 pm (for 3g1h) to 2200 pm (for 6g2w), reflecting the average radial system size R discussed above (see Table 3). Within this radial range, between the minimum around 850 pm and the radius where the bulk value is approached, the volume exclusion effect of the
5160 J. Phys. Chem. B, Vol. 112, No. 16, 2008
Brodbeck et al.
TABLE 4: Position Fluctuationsa si (in Units of 100 pm) for Selected Atom Types i in Each of the Eight Dendrizyme Model Systems (System Labels as Given in Table 1) atom type system
N
CNMe
HNMe
Cph
Hph
Cpy
Hpy
CtBu
HtBu
3g1h 3g2h 6g1h 6g2h
2.3-6.4 2.0-3.2 1.9-4.2 2.1-2.8
2.4-6.1 2.2-3.4 1.9-4.2 2.1-2.7
2.5-6.2 2.3-3.6 2.1-4.3 2.1-3.0
2.9-9.8 1.9-5.2 1.9-5.4 1.7-4.2
2.9-10.5 1.9-5.5 1.9-5.4 1.7-4.4
2.4-7.5 1.9-3.4 1.8-4.5 1.9-3.2
2.4-8.2 2.0-3.5 1.8-4.7 1.9-3.3
3.5-10.8 2.2-6.0 2.1-6.3 1.9-4.7
4.2-11.2 2.7-6.3 2.5-6.7 2.0-5.4
N
CMe
HMe
Cph
Hph
CCH2
HCH2
OOCH2
OOCH3
2.9-3.3 2.7-5.0 1.8-4.8 1.6-2.2
2.8-6.6 2.6-10.4 1.8-9.5 1.5-7.9
2.8-7.0 2.7-10.7 1.9-9.9 1.7-8.2
2.6-5.5 2.4-9.4 2.0-8.2 1.5-6.7
2.7-6.0 2.4-9.9 2.2-8.6 1.5-7.4
2.9-4.3 2.6-7.5 1.9-7.6 1.6-4.9
2.9-4.7 2.6-7.8 1.9-6.8 1.5-5.2
2.9-4.0 2.6-6.9 2.6-6.3 1.6-4.6
2.8-6.1 2.5-9.9 2.6-8.9 1.6-7.3
atom type 3g1w 3g2w 6g1w 6g2w
a By definition, si2 ) 〈xi2〉 - 〈xi〉2 + 〈yi2〉 - 〈yi〉2 + 〈zi2〉 - 〈zi〉2, where 〈‚〉 denotes the time average over the complete trajectory (ranges are obtained because several atoms belong to a chosen atom type).
dendrizyme shell on the frequency to find a trichloromethane molecule is clearly visible. Taking the maximum value of gXHCHCl3(r) in the radial range between 450 pm and 850 pm as an indicator of solvent permeability of the dendrizyme shell, the following detailed results are obtained for the eight systems under study: (i) For constant dendrimer generation k and constant number of dendrimers n in the dendrizyme shell, the permeability of the dendrizyme shell is always clearly higher for the system ngkh (with pyridinium dendrimers 3a or 3b) than for the system ngkw (with ammonium dendrimers 4a or 4b). This result supports the naı¨ve expectation that the dendrizyme shell formed by the conformationally more rigid pyridinium dendrimers is of lower density and has a smaller volume exclusion effect than the one obtained from the conformationally more flexible ammonium dendrimers. (ii) An increase of the number of dendrimers n (from n ) 3 to n ) 6, through a change from X ) P to X ) Zn), keeping dendrimer generation k and dendrimer type x constant, leads to a reduction of the permeability of the dendrizyme shell in most cases. An exception is the pair of systems 3g1h/6g1h: the latter system shows a slightly larger maximum value in gXHCHCl3(r) than the former, thus indicating a slightly larger permeability. This may be partly understood from an inspection of the timeaveraged structures of these two systems (see Figure 3): the fraction of shielded Keggin ion surface does not differ much in both systems, because the dendrimers are mainly tangentially oriented in the former system, but predominantly radially oriented in the latter system. Thus, the naı¨ve expectation that an increase in the number of dendrimers necessarily leads to a decrease in permeability is not true in general. (iii) An increase of the dendrimer generation k for constant number of dendrimers n and constant dendrimer type x leads almost always to a reduction of the permeability of the shell, as expected. The effect is most clearly observed for the systems with n ) 6 dendrimers in their shell, it is less pronounced for the pair of systems 3g1h/3g2h, whereas the pair of systems 3g1w/3g2w is an exception. An understanding of this result can be gained from data on position fluctuations for methyl group carbon and hydrogen atoms, given in Table 4. The upper values of the intervals given there, belonging to carbon and hydrogen atoms of the terminal methyl groups, indicate a higher mobility of these groups in 3g2w than in 3g1w. Thus, even though the Keggin ion in 3g2w is, on average, better shielded than that in 3g1w, the high conformational flexibility of the secondgeneration dendrimers 4b in 3g2w does not efficiently prohibit access of trichloromethane molecules to the dendrizyme core.
On the basis of these results, the eight systems studied here can be arranged into a sequence of decreasing the solvent permeability of the dendrizyme shell as follows:28 6g1h > 3g1h > 3g2h > 3g2w > 3g1w > 6g2h > 6g1w > 6g2w. This result differs considerably from the sequence 3g1h > 3g2h > 3g1w > 3g2w > 6g1h > 6g2h > 6g1w > 6g2w, naı¨vely expected from the estimated sterical demand of the isolated dendrimers. The radial distribution functions discussed so far necessarily suffer from a loss of information due to angular integration of the fully three-dimensional spatial distribution functions. In order to try to find additional details of the spatial arrangement of the trichloromethane molecules within the dendrizyme shell, we examined the pair of spatial distribution functions gOXO-A(r), with A ) Caryl and A ) HCHCl3, respectively. The former probes the dendrizyme shell, and the latter probes the solvent. Their combined analysis is expected to provide insight into the distribution of solvent molecules within the dendrizyme shell. A combined presentation of isovalue surfaces for these two spatial distribution functions is given in Figure 9 for the four systems 3gkx, and in Figure 10 for the four systems 6gkx. The grainy appearance of these isovalue surfaces is due to technical details of the trajectory analysis (discrete nature of raw data, finite size of volume elements for data analysis). The isovalues for the carbon distributions (in red) are in the range between 50 and 120, and those for the hydrogen distributions (in blue) lie between 8 and 12. These values were carefully chosen to obtain good contrast for both distributions. In particular, isovalues for hydrogen distributions corresponding to the almost constant value in the bulk solvent had to be avoided. The central O-X-O unit (in black), on which the definition of the coordinate system is based, is partly visible in all figures. Inspection of sequences of isovalue surfaces, including those shown in Figures 9 and 10, reveals that spatial regions with a high probability to find trichloromethane molecules are always separated from spatial regions with a high probability to find dendritic carbon centers. Cavities formed by the dendrimers and filled with trichloromethane molecules were not observed in any of the eight dendrizyme model systems included in this study. Hence, previous assumptions as to the presence of cavities in the dendrizyme shell7 cannot be confirmed. The trichloromethane molecules (and appropriately sized substrate molecules) can diffuse only to uncovered or unshielded parts of the surface of the dendrizyme core. Therefore, for a molecule of given size, the permeability of the dendrizyme shell is expected to be proportional to the fraction of the surface of the dendrizyme core left unshielded by the dendrizyme shell. (The situation partly resembles the one for water permeability of soil.
MD Study of Dendrimer-Encapsulated R-Keggin Ions
J. Phys. Chem. B, Vol. 112, No. 16, 2008 5161
Figure 10. Isovalue representation of the spatial distribution functions gOZnO-Caryl(r) (in red) and gOZnO-HCHCl3(r) (in blue) of aryl carbon centers and trichloromethane hydrogen centers, respectively, for dendrizyme model systems with Keggin ion 2. The central O-Zn-O unit (in black) is partly visible.
Figure 9. Isovalue representation of the spatial distribution functions gOPO-Caryl(r) (in red) and gOPO-HCHCl3(r) (in blue) of aryl carbon centers and trichloromethane hydrogen centers, respectively, for dendrizyme model systems with Keggin ion 1. The central O-P-O unit (in black) is partly visible.
4. Conclusions
It is known that the porosity of a soil determines the amount of water the soil can hold, and also affects the permeability of the soil, i.e., the ability of water to flow through the soil.) We conclude that the experimentally observed substrate selectivity of a prototype dendrizyme system in olefin epoxidation9 may be partly due to the size exclusion effect of the dendrizyme shell. Certainly, the accessible part of the POM surface depends on the type and number of dendrimers in the surrounding shell, and is subject to fluctuations due to dynamical processes of dendritic molecular chains within the dendrizyme shell.
This work presents results from MD simulations within the NVE ensemble at temperatures around T ) 300 K for eight dendrizyme model compounds in trichloromethane solution. These dendrizyme model compounds are non-charged ion clusters composed of a central POM (the core), an R-Keggin ion [(XO4)W12O36]n- (1: X ) P, n ) 3; 2: X ) Zn, n ) 6, see Figure 1), and surrounding dendritic ammonium or pyridinium ions (forming the shell, see Figure 2). Similar compounds, with a redox active POM in the center, can exhibit enzyme-like regioselectivity and substrate selectivity in olefin epoxidation,9 thus mimicking in part the structure-activity profile of a
5162 J. Phys. Chem. B, Vol. 112, No. 16, 2008 metalloenzyme: the POM controls the regioselectivity of the reaction, whereas the dendritic shell influences (controls) access of potential substrates. A true dendrizyme satisfying all necessary conditions (regioselectivity, substrate selectivity, stability under reaction conditions, etc.) has yet to be made. In this work, we sought theoretical criteria, derived from MD simulations on dendrizyme model compounds, allowing us to discriminate different dendrizyme shells with respect to their influence on the accessibility of the dendrizyme core for small molecules. Simulations involving many different potential substrates were outside the scope of the present work. Such studies, involving an additional minority compound in the simulation, remain to be done in the future. The MD simulations used an improved, newly parametrized Keggin ion force field, and the solvent was explicitly included. Hence, the results are expected to be more reliable than those obtained from gas-phase simulations.7 Within the chosen set of eight compounds, the influence of dendrimer type, dendrimer generation, and number of coordinated dendrimers on the structure and dynamics of the dendrizyme shell could be studied. The radial distribution function gXA(r), where X is the central heteroatom of the Keggin ion and A is the hydrogen atom of trichloromethane, was found to be a source of useful data: the height of the nearest maximum of this function (see Figure 8) can serve as a measure of the permeability of the dendrizyme shell. The analysis of the MD simulation data also included timeaveraged structures and spatial distribution functions for dendrimer carbon centers and hydrogen centers of trichloromethane, respectively (their isovalue surfaces are shown in Figures 9 and 10). This lead us to conclude that substrate selectivity is influenced by the size exclusion effect of the dendrizyme shell, depending primarily on the permeability of the dendrizyme shell and on the unshielded part of the surface of the dendrizyme core. No indications were found for the formation or presence of solvent-containing cavities in the dendrizyme shell. We expect our conclusions to remain valid as long as thermodynamical effects are of lesser importance. However, thermodynamical effects immediately gain importance as soon as the dendritic molecular chains take active part in the transport of substrate molecules toward the dendrizyme core, or in the catalytic reaction itself. Maybe only then can the term “dendrizyme” be appropriately used for members of the class of compounds studied here. Acknowledgment. The authors would like to thank Paul Ko¨gerler (formerly at Iowa State University, Ames, Iowa, now at RWTH Aachen, Germany) for providing access to the Amber suite of programs and for generously granting computer time at the Ameslab. Financial support from the German Research Foundation (DFG, AN 335/2) is gratefully acknowledged. We also thank the Leibniz Computing Center Munich (LRZ) and the Rechenzentrum Garching (RZG) for granting access to their computing facilities. Supporting Information Available: Technical details of the Kohn-Sham density functional calculations and tables providing the following data: exponents of polarization and diffuse Gausstype basis functions used in calculation on the Keggin ions; electronic energies at equilibrium structures, Cartesian coordinates of equilibrium structures, atom types and RESP partial charges for the Keggin ions 1, 2, and the dendrimers 3a-4b; Amber force field parameters for trichloromethane; translation diffusion coefficients of trichloromethane; a plot demonstrating the energy conservation in the NVE ensemble simulations. This
Brodbeck et al. information is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Hill, C. L.; Prosser-McCartha, C. M. Coord. Chem. ReV. 1995, 143, 407-455. (2) Hill, C. L., Ed. J. Mol. Catal. A 1996, 114, 1-371 (special volume on polyoxometalates in catalysis). (3) Neumann, R. Prog. Inorg. Chem. 1998, 47, 317-370. (4) Neumann, R. Liquid phase oxidation reactions catalyzed by polyoxometalates. In Modern Oxidation Methods; Ba¨ckvall, J.-E., Ed.; Wiley-VCH: Weinheim, Germany, 2004. (5) Hill, C. L., Ed. J. Mol. Catal. A 2007, 262, 1-242 (special volume on polyoxometalates in catalysis). (6) Gouzerh, P.; Proust, A. Chem. ReV. 1998, 98, 77-112. (7) Volkmer, D.; Bredenko¨tter, B.; Tellenbro¨ker, J.; Ko¨gerler, P.; Kurth, D. G.; Lehmann, P.; Schnablegger, H.; Schwahn, D.; Piepenbrink, M.; Krebs, B. J. Am. Chem. Soc. 2002, 124, 10489-10496. (8) Brunner, H. J. Organomet. Chem. 1995, 500, 39-46. (9) Volkmer, D. To be submitted for publication. (10) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. Amber 7; University of California: San Francisco, CA, 2002. (11) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157-1174. (12) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (13) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577-8593. (14) Attempts to generate a force field for a larger class of tungsten POMs, ranging from the Lindqvist ion, [(O)W6O18]2-, over the R-Keggin and R-Wells-Dawson ions to the Preyssler ion, [(PO4)5W30O90]15-, were only partly successful.29 (15) Halgren, T. A. J. Am. Chem. Soc. 1992, 114, 7827-7843. (16) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129145. (17) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269-10280. (18) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (19) Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652. (20) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B: Condens. Matter 1988, 37, 785-789. (21) Stevens, W. J.; Basch, H.; Krauss, M. J. Chem. Phys. 1984, 81, 6026-6033. (22) Stevens, W. J.; Krauss, M.; Basch, H.; Jasien, P. G. Can. J. Chem. 1992, 70, 612-630. (23) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (24) To¨nsing, T. Ph.D. Thesis, Department of Chemistry, University of Bielefeld, 2004 (in German; electronic version available at http://bieson.ub.unibielefeld.de/volltexte/2004/562/). (25) To¨nsing, T.; Oldiges, C. Phys. Chem. Chem. Phys. 2001, 3, 55425549. (26) Portmann, S.; Lu¨thi, H. P. Chimia 2000, 54, 766-770. (27) Flu¨kiger, P.; Lu¨thi, H. P.; Portmann, S.; Weber, J. Molekel 4.3; Swiss Center for Scientific Computing: Manno, Switzerland, 2000-2002. (28) We do not expect a change of the criterion for substrate permeability to induce much rearrangement in this sequence. Instead of the nearest maximum value of gXHCHCl3(r) we could have used an integral of 4πr2gXHCHCl3(r) over some suitably chosen radial range. (29) Brodbeck, R. Ph.D. Thesis, Department of Chemistry, University of Bielefeld, 2006 (in German; electronic version available at http:// bieson.ub.uni-bielefeld.de/volltexte/2007/1086/).