Molecular Dynamics Simulation of Cationic and Anionic Clays

Langmuir 2002, 18, 2933-2939

2933

Molecular Dynamics Simulation of Cationic and Anionic Clays Containing Amino Acids Steven P. Newman, Tiziana Di Cristina, and Peter V. Coveney* Centre for Computational Science, Queen Mary, University of London, Mile End, London, E1 4NS, United Kingdom

William Jones* Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom Received September 19, 2001. In Final Form: January 9, 2002 The detailed arrangements of guest species in the interlayers of a cationic and an anionic clay have been determined using molecular dynamics computer simulations. The guest species studied were the amino acids (S)-phenylalanine and (S)-tyrosine, chosen for comparison because both anionic and cationic forms are possible. A modified version of the Dreiding force field was used for all models. Our simulations provide for the first time a clear picture of the different arrangements in the interlayer of the two clay systems. Previous work using powder X-ray diffraction had allowed only possible arrangements to be proposed based on purely geometric considerations. In the anionic clay case, the interlayer anions assume a bilayerlike arrangement with their longest molecular axis oriented in an approximately perpendicular direction with respect to the hydroxide layers. In the montmorillonite system, the interlayer cations prefer a monolayer arrangement with their longest molecular axis parallel to the layers. We believe the difference to be a direct consequence of the different layer charge densities of the clays. In all cases, the water molecules are generally excluded from the hydrophobic midplane formed by the phenyl rings of the counterions. Given the poorly crystalline nature of such intercalates, computer simulations of the type reported here are likely to provide the only means of knowing such details of the interlayer arrangements.

Introduction Layered materials that are able to intercalate neutral guest molecules or to exchange inorganic and organic ions for interlayer ions have attracted considerable attention for many years.1 In principle, through the incorporation of a guest species into a layered host, novel solids may be engineered with desirable physical and chemical properties. Various layered materials such as clay minerals, graphite, transition metal dichalcogenides, and metal phosphates and phosphonates have the ability to act as host materials.2-4 A wide variety of organic and inorganic guest molecules are also known. By appropriate choice of host and guest, hybrid materials can be created that exhibit enhanced properties compared to the individual components. A key requirement of the chemistry of layered materials is detailed knowledge of the precise arrangement of guest molecules within the interlayer and the types of guestguest and guest-layer interactions. For clay systems, such information is not generally obtainable from experimental measurements. Due to the inherent structural disorder and generally low crystallinity of clays, complete structure determination using single-crystal or powder X-ray diffraction is rarely possible, for example. Structure elucida* Corresponding authors. E-mail: [email protected], [email protected]. (1) Intercalation Chemistry; Whittingham, M. S., Jacobson, A. J., Eds.; Academic Press: New York, 1982. (2) Comprehensive Supramolecular Chemistry: Solid-state Supramolecular Chemistry: Two and Three-Dimensional Inorganic Networks; Lehn, J.-M., Atwood, J. L., Davies, J. E. D., Macnicol, D. D., Vogtle, F., Alberti, G., Bein, T., Eds.; Pergamon Press: Oxford, 1996. (3) Inorganic Materials; Bruce, D. W., O’Hare, D., Eds.; John Wiley & Sons Ltd.: Chichester, 1997. (4) Clearfield, A. Chem. Rev. 1988, 88, 125.

tion is normally reduced to inferring the interlayer arrangement of the clay from the interlayer separation (which may be obtained from powder X-ray data) and geometric considerations.5 In the present work, however, molecular dynamics computer simulations are used to probe the interlayer of two complementary types of clay systems: layered double hydroxides (LDHs) and smectite clays. These are examples of naturally occurring layered minerals often referred to as anionic and cationic clays, respectively: in LDHs, the host layers are positively charged, thus requiring the presence of interlayer anions to maintain overall charge neutrality, whereas in smectite clays the host layers are negatively charged and interlayer cations are required. One important group of LDHs may be represented by the general formula

[M1-x2+Mx3+(OH)2]Ax/nn-‚mH2O, where M2+ and M3+ are divalent and trivalent cations, respectively, and A is an anion of valence n. The structure of such LDHs may be described by considering the structure of brucite, Mg(OH)2, which consists of chargeneutral layers of edge-sharing hydroxide octahedra, with Mg2+ cations occupying the octahedral sites. In a LDH, isomorphous replacement of a fraction of the Mg2+ ions with a trivalent cation, such as Al3+, occurs and generates a positive charge on the layers which necessitates the presence of interlayer, charge-balancing anions. Water of crystallization also occupies the interlayer space of LDHs. Observed M2+ and M3+ species include Mg2+, Fe2+, Co2+, (5) Costantino, U.; Coletti, N.; Nocchetti, M.; Aloisi, G. G.; Elisei, F. Langmuir 1999, 15, 4454.

10.1021/la0114528 CCC: $22.00 © 2002 American Chemical Society Published on Web 02/28/2002

2934

Langmuir, Vol. 18, No. 7, 2002

Cu2+, Ni2+, or Zn2+ and Al3+, Cr3+, Ga3+, or Fe3+, respectively.6 The combination of Mg and Al has been the most frequently studied, with a variety of Mg/Al ratios and different charge-balancing anions. In naturally occurring LDHs, the most commonly found interlayer anion is carbonate. In theory, however, there is no significant restriction to the identity of the charge-balancing anion that can occupy the interlayer region. Examples include inorganic (such as halides, oxo-anions, silicates, and polyoxometalates) as well as organic anions (such as alkyl or aryl carboxylates and sulfonates).6,7 Smectite clays consist of negatively charged, crystalline aluminosilicate sheets. Montmorillonite is one example of a naturally occurring smectite. Its structure comprises stacks of pyrophyllite-like sheets, [Al2Si4O10(OH)2], which are charge-neutral and consist of layers of SiO4 tetrahedra above and below an Al(OH)O2 octahedral layer. The SiO4 and Al(OH)O2 layers are linked via O atoms. In montmorillonite, a negative charge arises from partial replacement of Al by Mg in the octahedral (Al(OH)O2) layer, thus requiring the presence of interlayer cations to balance the layer charge. The interlayer also contains variable water content. A general formula is thus [MgxAl2-x(OH)2(Si4O10)](Cn+)x/n‚nH2O, where C is a cation of valence n. Various types of cations may be incorporated into the interlayer of smectite clays. A typical naturally occurring montmorillonite will commonly have Na+, Ca2+, or Mg2+ in the interlayer, for example. It is generally possible, however, to exchange these with other cations, such as transition-metal cations and organic cations.3 In the present work, molecular dynamics computer simulations are used to probe the interlayer properties of model MgAl-LDHs and montmorillonites containing the amino acids (S)-phenylalanine and (S)-tyrosine. The fact that amino acids are zwitterionic means that they can be incorporated into both clay systems with opposite charge. The models were chosen deliberately to echo an earlier experimental study because of the ability of simulations to provide detailed insight into the interlayer properties of the clay systems not available from the experimental measurements.8 A particularly interesting aspect of the present work is that the same modified version of the Dreiding force field can be used for the simulations of both clay systems. Although computer simulations of smectite clays and LDHs have been previously reported,9-19 this work (6) Newman, S. P.; Jones, W. In Supramolecular Organisation and Materials Design; Rao, C. N. R., Jones, W., Eds.; Cambridge University Press: Cambridge, U.K., 2001. (7) Newman, S. P.; Jones, W. New J. Chem. 1998, 22, 105. (8) Fudala, A.; Pa´linko´, I.; Kiricsi, I. Inorg. Chem. 1999, 38, 4653. (9) Aicken, M.; Bell, I. S.; Coveney, P. V.; Jones, W. Adv. Mater. 1997, 9, 496. (10) Newman, S. P.; Williams, S. J.; Coveney, P. V.; Jones, W. J. Phys. Chem. B 1998, 102, 6710. (11) Williams, S. J.; Coveney, P. V.; Jones, W. Mol. Simul. 1999, 21, 183. (12) Newman, S. P. Ph.D. Thesis, University of Cambridge, Cambridge, U.K., 1999. (13) Fogg, A. M.; Rohl, A. L.; Parkinson, G. M.; O’Hare, D. Chem. Mater. 1999, 11, 1194. (14) Boek, E. S.; Coveney, P. V.; Skipper, N. T. J. Am. Chem. Soc. 1995, 117, 12608. (15) Boek, E. S.; Coveney, P. V.; Skipper, N. T. Langmuir 1995, 11, 4629. (16) de Siqueira, A. V. C.; Skipper, N. T.; Coveney, P. V.; Boek, E. S. Mol. Phys. 1997, 92, 1. (17) Bains, A. S.; Boek, E. S.; Coveney, P. V.; Williams, S. J.; Akbar, A. V. Mol. Simul. 2001, 26, 101. (18) Keldsen, G. L.; Nicholas, J. B.; Carrado, K. A.; Winans, R. E. J. Phys. Chem. 1994, 98, 279. (19) Teppen, B. J.; Yu, C.-H.; Miller, D. M.; Scha¨fer, L. J. Comput. Chem. 1998, 19, 144.

Newman et al.

represents the first simultaneous study of such systems, thus facilitating direct comparison of their properties. Simulation Method All simulations were performed using the Cerius2 software package, generally in parallel across two nodes of a SGI Onyx2 computer.20 Model Construction. The first stage in preparing LDH and montmorillonite models for computer simulations is construction of the clay framework, that is, the positively charged hydroxide layers for LDHs and the negatively charged aluminosilicate layers for montmorillonite. The hydroxide layers of the model LDHs used in the simulations were constructed using atomic coordinates from the previously reported crystal structure of hydrotalcite, a naturally occurring LDH mineral with the composition [Mg4Al2(OH)12](CO3)‚3H2O.21 The interlayer carbonate and water molecules were removed, and a P1 superlattice was created with a two-layer repeat and unitcell parameters a ) 18.32, b ) 12.22, c ) 36.00 Å, R ) 90.0, β ) 90.0, γ ) 120.0° (equivalent to 24 hydrotalcite unit cells in the ab-plane and an interlayer spacing of 18.0 Å). The Mg/Al ratio was adjusted to 3, such that each hydroxide layer contains 18 Mg2+ and 6 Al3+ atoms, and the distribution of the two cations was ordered such that the Al atoms did not occupy adjacent hydroxide octahedra. Note that the single-crystal data suggest that in hydrotalcite the Mg2+ and Al3+ cations are randomly distributed among the octahedral sites within the close-packed hydroxide layers (2/3 and 1/3 site occupancy of the octahedral cation positions, respectively). An ordered distribution was used for the model LDHs in the present work, however, because each atom must have unit site occupancy and the periodicity of the cell makes it impossible to generate a truly random distribution of cations in the layers. Atomic charges were initially calculated using the charge equilibration method.22 This method can be applied only to neutral systems, however; thus, the required positive charge of +12 was subsequently averaged over all atoms in the cell. The composition of the model LDH framework (hereafter referred to as the model Mg3AlLDH) is thus [Mg36Al12(OH)96]12+, with the following partial charges: +0.71 for Mg, +1.35 for Al, -0.56 or -0.51 for O, and +0.23 or +0.25 for H depending on the environment of the atom. Constructing a model montmorillonite framework was not as straightforward because crystal-structure data are not available (due, at least in part, to the disordered and nonstoichiometric nature of the clay). Instead, it was necessary to adapt the model from crystal-structure data for the mineral muscovite.20 This mineral comprises aluminosilicate layers with the same basic structure as those of montmorillonite, although there is no substitution of octahedral Al for Mg in muscovite. A P1 superlattice was created with a one-layer repeat and unit-cell parameters a ) 20.78, b ) 35.98, c ) 36.00 Å, R ) 90.0, β ) 95.2.0, γ ) 90.0° (corresponding to 16 muscovite unit cells in the ab-plane). Nine Al atoms were then substituted for Mg such that the Mg atoms did not occupy adjacent hydroxide octahedra (see above comments regarding the cation distribution in the model LDH) and a charge of -9 was averaged over all atoms in the framework, following an initial calculation of atomic charges using the charge equilibration method. The composition of the model (20) Cerius2, version 4.2; Accelrys: San Diego, CA, 2001. (21) Allmann, R.; Jepsen, H. P. Neues Jahrb. Mineral., Monatsh. 1969, 12, 544. (22) Rappe´, A. K.; Goddard, W. A., III J. Phys. Chem. 1991, 95, 3358.

Simulation of Clays Containing Amino Acids

Langmuir, Vol. 18, No. 7, 2002 2935

montmorillonite framework is thus [Mg9Al55(OH)64(Si128O320)]9-, with the following partial charges: between +0.72 and +0.73 for Mg; +1.17 and +1.22 for Al; +1.26 and +1.33 for Si; +0.18 and +0.25 for H; -0.60 and -0.85 for O, depending on the environment of the atom. The second stage in preparation of the models for the simulations is construction of the interlayer molecules. The initial geometry and atomic charges of anionic and cationic forms of (S)-phenylalanine and (S)-tyrosine were calculated using the MOPAC semiempirical molecular orbital method, employing the AM1 Hamiltonian approximation (see Supporting Information).23 For water, the initial geometry and atomic charges (+0.417 for H, -0.834 for O) were taken from the TIP3P specialized water force field.24 Recent molecular dynamics simulations of bulk water have found that a combination of the Dreiding force field, as used here with the TIP3P partial charges for water, yields physical properties, such as density and radial distribution function, in good agreement with experimental values.25 The required number of the anionic or cationic forms of the amino acids were then placed within the model LDH or montmorillonite frameworks, respectively. In each case, the molecules were placed in the frameworks with their longest axis oriented in a direction approximately perpendicular to the plane of the layers in a bilayer-like arrangement. Finally, the desired number of water molecules were placed at random in the interlayer of the clays. In the work of Fudala et al., the water content of the clays is not reported.8 Simulations were performed, therefore, using models with different interlayer water contents until good agreement was obtained between the simulated and experimental interlayer spacings. Note also that Fudala et al. investigated ZnAl-LDHs, whereas our models are of MgAl-LDHs. To a first approximation, however, we do not consider the identity of the divalent cation to be an important factor in determining the interlayer structure. Simulation Details. A modified version of the Dreiding force field was used for all simulations.26 Dreiding is a generic force field, designed for application to a wide range of systems (organic, biological, and main-group inorganic molecules) containing a variety of different atoms. The philosophy behind Dreiding is to use general bonded parameters (e.g., force constants) based on the hybridization state of atoms, rather than the specific combination of atoms involved in a particular system. For the simulations described here, a number of modifications were made to the Dreiding force field to tailor it for the clay systems. First, it was necessary to include parameters for Mg, which are not present in the original Dreiding force field. This was achieved by assigning Mg the same parameters as Al (for which Dreiding does have parameters), with the exception that the equilibrium bond lengths with oxygen are different. In our previous simulations of LDHs, a Mg-O bond length of 2.100 Å, which is the mean Mg-O bond length in the crystal structure of brucite, Mg(OH)2, and a Al-O bond length of 1.697 Å, which is the value assigned by Dreiding, were used.9-12 In the present work, we retained the Mg-O bond length but used a length of 1.902 Å for Al-O, which is the mean Al-O bond length

in the crystal structure of gibbsite, Al(OH)3. This modification was made in order to counteract the contraction in the ab-plane of the hydroxide layers observed in our previous simulations of LDHs. Additionally, we used a value of 1.644 Å for the Si-O bond length, which is the average Si-O bond length in the crystal structure of muscovite, rather than a length of 1.587 Å assigned in Dreiding. In our previous simulations of LDHs, a significant distortion of the hydroxide octahedra was observed.9-12 It has been found subsequently, however, that this distortion can be reduced significantly by replacing the harmonic theta form for angle-bend energy terms, E ) (1/2)K[θ θeq]2, with a harmonic cosine theta form, E ) 1/2C[cos θ - cos θeq]2, where E is the potential energy, K is a force constant, θ is the bond angle, θeq is the equilibrium bond angle, and C is given by C ) K/(sin θeq)2. The authors of Dreiding recommend use of the harmonic cosine theta form, although the harmonic theta form is the default used in the Cerius2 implementation of the force field, presumably to reduce computation time.20 In the present work, we used the harmonic cosine theta term for all anglebend energy terms in all simulations. It is probable that the improved modeling of the hydroxide octahedra observed using this term arises because it yields a zero potential energy gradient at θ ) 180°, whereas the harmonic theta form does not. Long-range Coulombic interactions and the attractive van der Waals interactions were computed using the Ewald summation technique. Short-range repulsive van der Waals interactions were treated with a direct cutoff radius of approximately 6.6 Å for the LDH models and 7.2 Å for montmorillonite. Energy minimization was performed on all models (with no rigid bodies or cell constraints included) prior to the molecular dynamics simulations to reduce strain.27 Molecular dynamics simulations were performed in the constant composition, isothermal-isobaric (NPT) ensemble at 300 K. A time-step of 0.001 ps was used. Temperature was maintained using the Hoover thermostat implemented with a relaxation time of 0.1 ps and a cell-mass prefactor of unity. The equivalent hydrostatic pressure was set to 105 Pa (approximately 1 atm). Periodic boundary conditions were applied in three dimensions so that the simulation cell is effectively repeated infinitely in each direction; each atom in the cell therefore interacts with its periodic image as well as all the other atoms in the cell.27 The simulations were performed in two stages: first equilibration of the system and second data collection of system properties. Whether a system has equilibrated may be judged by plotting various thermodynamic quantities such as energy or temperature against time. When equilibration has been reached, these quantities fluctuate around their average values, which remain constant over time. In the present work, the total simulation time was 50 ps (equivalent to 50 000 time-steps) for each model, with the first 40 ps of the simulation considered as the equilibration period. For each simulation, a snapshot of the model was stored every 0.05 ps (every 50 time-steps) and the interlayer spacing of the model was averaged over the final 10 ps (the average interlayer spacing from the final 200 snapshots), after thermal equilibration is adjudged to have been reached.

(23) Stewart, J. J. P. MOPAC, version 6.0; QCPE No. 455, Department of Chemistry; Indiana University: Bloomington, IN, 1990. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (25) Boek, E. S.; Coveney, P. V.; Williams, S. J.; Bains, A. Mol. Simul. 1996, 18, 145. (26) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III J. Phys. Chem. 1990, 94, 8897.

Results and Discussion Mg3Al-LDHs Containing (S)-Phenylalanine and (S)-Tyrosine. A snapshot of the simulation cell of a model (27) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications: Oxford, 1994.

2936

Langmuir, Vol. 18, No. 7, 2002

Newman et al.

Figure 2. Possible hydrogen-bond interactions (broken green lines) in a model Mg3Al-LDH containing (S)-phenylalanine, following 50 ps of molecular dynamics simulation at 300 K. The criteria for hydrogen-bond assignment are that the maximum H‚‚‚A distance is 2.5 Å and the minimum D-H‚‚‚A angle is 90°, where D and A are hydrogen-bond donors or acceptors, respectively. Here, N and O are allowed to act as either donors or acceptors, whereas C is allowed to act as neither.

Figure 1. Snapshot of the simulation cell of a model Mg3AlLDH containing (S)-phenylalanine, following 50 ps of molecular dynamics simulation at 300 K: blue ) Mg and Al, red ) O, gray ) C, yellow ) N, and white ) H.

Mg3Al-LDH containing (S)-phenylalanine, following 50 ps of molecular dynamics simulation at 300 K, is shown in Figure 1. Sixty water molecules were included in the simulation cell (30 per interlayer), which corresponds to 18.4 mass % and gives an empirical formula of [Mg3Al(OH)8]((S)-phenylalanine)‚5H2O. The average simulationcell parameters in the ab-plane over the last 10 ps of the simulation are a ) 17.68, b ) 11.78 Å, and γ ) 120.0°. Comparison of these values with the initial simulationcell parameters reveals that only a small contraction of the hydroxide layers has occurred in the ab-plane during the simulation. The simulated cell parameters are within 5% of the initial values, thus suggesting that the force field is parametrized adequately to model the hydroxide layers. Furthermore, there is very little distortion of the hydroxide octahedra compared to previous simulations of LDHs, which also validates the modifications made to the force field. The average interlayer spacing over the last 10 ps of the simulation is 17.9 ( 0.4 Å, which is in good agreement with the interlayer spacing of 18.0 Å determined by Fudala et al. for a synthetic ZnAl(L-phenyalanine) LDH using powder X-ray diffraction.8 Surprisingly, the authors do not report the interlayer water content of the synthetic LDH, which hinders direct comparison between the experimental and simulation results. It has been shown previously, however, that simulations of the kind described here reproduce accurately the experimentally determined correlation between interlayer spacing and water content for LDHs containing organic anions.8

Figure 3. View along [001] of a model Mg3Al-LDH containing (S)-phenylalanine anions following 50 ps of molecular dynamics simulation at 300 K. Only the phenyl ring of the anions in one interlayer of the model is displayed (four simulation cells in the ab-plane are displayed). Centroids of the phenyl rings are indicated in black. Distances are in Å. Two short C-H‚‚‚π contacts are also indicated (the C-H‚‚‚centroid angles are 140.8 and 173.4° for the 4.06 and 2.96 Å contacts, respectively).

The snapshot shows that the (S)-phenylalanine molecules retain the interdigitated bilayer arrangement in which they were placed initially, with the longest axis of the molecule generally oriented in an approximately perpendicular direction from the plane of the hydroxide layers. The phenyl rings of the molecules thus form a hydrophobic region in the midplane of the interlayer from which nearly all of the water molecules are excluded, their preferring instead to form layers close to the surface of the hydroxide layers. For each (S)-phenylalanine molecule, the closest contact to the hydroxide layers is made through the oxygen atoms of the carboxylate group, as should be expected given the positive charge on the layers. The oxygen atoms of the carboxylate groups and water molecules thus appear to form an oxygen monolayer on the surface of the hydroxide layers. In this region, a number of close contacts occur between the water molecules, the carboxylate groups, and the hydroxide surface that are consistent with the existence of hydrogen-bond interactions (Figure 2). The hydrogen-bond, H‚‚‚O, distance is typically 2.2 Å with a variable O-H‚‚‚O angle,

Simulation of Clays Containing Amino Acids

Langmuir, Vol. 18, No. 7, 2002 2937

Figure 4. Snapshot of the simulation cell of a model montmorillonite containing (S)-phenylalanine, following 50 ps of molecular dynamics simulation at 300 K: blue ) Mg and Al, green ) Si, red ) O, gray ) C, yellow ) N, and white ) H.

where the hydrogen-bond donor is a surface OH group and the acceptor is an O atom belonging to a water molecule or carboxylate group. The amine group is also involved in what might be considered hydrogen-bond interactions with water molecules, although generally not with the hydroxide surface. The relative orientation of the interlayer anions is displayed in Figure 3. The view is along [001] with only the phenyl ring of the anions in one interlayer of the model displayed. The separation between the anions is typically 6 Å, as indicated by the selected centroid-centroid distances. Furthermore, an edge-to-face orientation is generally preferred, rather than a face-to-face orientation, possibly indicative of weak but directional C-H‚‚‚π hydrogen bonding (π-electron density is of course not accounted for in classically based simulations; any interaction is simply Coulombic between an H atom with a positive partial charge and the C atoms in the phenyl ring bearing negative partial charges).28 For Mg3Al-LDH models containing (S)-tyrosine, it was found that inclusion of 60 water molecules in the simulation cell (30 per interlayer) yielded an average interlayer spacing over the last 10 ps of the simulation of 18.0 ( 0.2 Å. This value is slightly larger than the interlayer spacing of 17.5 Å determined by Fudala et al. for a synthetic Zn3Al-LDH containing L-tyrosine using powder X-ray diffraction.8 Sixty water molecules per cell corresponds to 17.8 mass % and gives an empirical formula of [Mg3Al(OH)8]((S)-tyrosine)‚5H2O. Presumably, this small difference between simulated and observed interlayer spacing may be attributed to a small difference between the water content of the model and real systems or the presence of impurities in the latter. Again, the mean simulation-cell parameters in the ab-plane are within 5% of the initial values. The geometrical arrangement of the cell is generally comparable to that of the (S)-phenylalanine system, the only significant difference being additional hydrogen-bond-type interactions between water molecules and the OH group belonging to (S)-tyrosine. Montmorillonite Containing (S)-Phenylalanine and (S)-Tyrosine. A snapshot of the model montmorillonite containing interlayer (S)-phenylalanine cations following 50 ps of molecular dynamics simulation at 300 K is shown in Figure 4. Thirty water molecules were included in the simulation cell, corresponding to 3.99 mass %. The average simulation-cell parameters in the ab-plane (28) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond in Structural Chemistry and Biology; Oxford University Press: New York, 1999.

Figure 5. Possible hydrogen-bond interactions (broken green line) in a model montmorillonite containing (S)-phenylalanine, following 50 ps of molecular dynamics simulation at 300 K. The criteria for hydrogen-bond assignment are that the maximum H-A distance is 2.5 Å and the minimum D-H‚‚‚A angle is 90°, where D and A are hydrogen-bond donors or acceptors, respectively. Here, N and O are allowed to act as either donors or acceptors, whereas C is allowed to act as neither.

over the last 10 ps of the simulation are a ) 18.36, b ) 32.12 Å, and γ ) 88.92°. Comparison of these parameters with the initial cell parameters reveals that a significant contraction of layers occurs for this montmorillonite model: the simulated a- and b-parameters are within 12% of the initial values. This is a larger contraction than observed for the LDH model and suggests that the force field parametrization for the montmorillonite model is not entirely adequate. Despite the contraction, however, we considered that the parametrization is sufficiently accurate to enable useful information about the arrangement of the interlayer molecules to be extracted. The average interlayer spacing over the last 10 ps of the simulation is 14.4 ( 0.1 Å, which is in good agreement with the value of 14.9 Å determined by Fudala et al. for a montmorillonite clay containing L-phenylalanine cations.8 This is a significant reduction from the initial interlayer spacing of 36.00 Å. It is accompanied by a change of the interlayer arrangement of the phenylalanine cations from the initial vertical bilayer in which they were placed to a monolayer arrangement. This should be compared to the LDH models, in which the interlayer anions retain their initial vertical bilayer arrangement. The difference between the models can be explained by the different

2938

Langmuir, Vol. 18, No. 7, 2002

Newman et al.

Figure 6. View along [001] of a model montmorillonite containing (S)-phenylalanine following 50 ps of molecular dynamics simulation at 300 K. Only the interlayer cations in one simulation cell are displayed. Centroids of the phenyl rings are indicated in black. Distances are in angstroms.

charge density on the surface of the layers of the two clay systems: 6.0 × 10-3 and 15 × 10-3 e/Å2 for the model montmorillonite and LDH layers, respectively. A greater number of charge-balancing molecules are thus required per unit area of surface in the LDH system, and a more expanded interlayer arrangement (in the c-direction) should be expected due to increased intermolecular repulsion. In the montmorillonite interlayer, the phenylalanine cations adopt an orientation with their longest axis approximately parallel to the plane of the layers. In general, the phenyl ring of the cations is tilted at an angle of approximately 45° with respect to the plane of the layers. In a manner similar to that described for the LDH system, water molecules form layers close to the oxide surface of the clay and are largely excluded from the hydrophobic region formed by the phenyl rings in the midplane of the interlayer. The positively charged amine groups of the cations are also found close to the oxide surface, as should be expected given the negative charge on the clay layers. Hydrogen-bond-type interactions are observed between the amine groups, water molecules, and oxide surface (Figure 5). Furthermore, in certain cases the carboxylic acid group proton is also involved in hydrogen-bond-type interactions to the surface. The relative orientation of the cations is shown in Figure 6. The view is along the [001] direction, with only the cations displayed. Selected centroid-centroid distances are displayed to indicate the separation of the cations in the interlayer. Qualitatively, there appears to be less structure in the arrangement of the interlayer for this model, compared to the equivalent LDH example. Presumably, this is due to the greater availability of surface in the cationic clay making intermolecular interactions between the cations less important. Note that evidence is not found for the formation of dimers via the carboxylic acid groups, as is often observed in the crystal structure of organic molecules bearing this functional group. For montmorillonite models containing (S)-tyrosine, it was found that inclusion of 30 water molecules in the

simulation cell (3.95 mass %) yielded an average interlayer spacing over the last 10 ps of the simulation of 14.3 ( 0.1 Å. This value is in good agreement with the interlayer spacing of 14.6 Å determined by Fudala et al. for a montmorillonite clay containing L-tyrosine.8 Contraction of the layers in the ab-plane is similar to that of the (S)phenylalanine example. The geometrical arrangement of the cell is also generally comparable to that of the (S)phenylalanine example, the only significant difference being additional hydrogen-bond-type interactions between water molecules and the OH group belonging to (S)tyrosine. Concluding Remarks Molecular dynamics computer simulations have been used to probe the interlayer structure of a model montmorillonite and a model Mg3Al-LDH containing the amino acids (S)-phenylalanine and (S)-tyrosine. The simulations provide detailed insight into the arrangement of the counterions and water molecules in the interlayer and suggest possible guest-guest and guest-layer interactions. Such information is particularly useful since it is not generally available from experiment for these systems, due largely to their poor crystallinity. Costantino et al., for example, have recently investigated the fluorescence of the azoic dye methyl orange in a LDH.5 In such an investigation, the disposition of the interlayer molecules is an important consideration, but it could only be inferred from the measured interlayer spacing and simple geometric considerations. Photochemical reactions occurring in the interlayers of clays is another area where detailed knowledge of the interlayer arrangement is required. Shichi and co-workers, for example, demonstrated the strong dependence of photochemical products on the interlayer arrangement of LDHs containing photoactive molecules.29,30 Once again, however, the interlayer arrangements were inferred only from geometric consider(29) Shichi, T.; Takagi, K.; Sawaki, Y. J. Chem. Soc., Chem. Commun. 1996, 2027.

Simulation of Clays Containing Amino Acids

ations. Such investigations would clearly benefit from simulations of the type described in the present work. A modified version of the Dreiding force field was used for all models, demonstrating its general applicability to clay systems. Although the parametrization is not optimal, in particular for the montmorillonite model, the robustness of the force field nonetheless facilitates extraction of useful qualitative and semiquantitative structural information. The close correspondence between the simulated interlayer spacing of the models and experimental values determined previously gives confidence about the force field parametrization. Although it may be inferred from the experimentally determined interlayer separation, the simulations demonstrate clearly that the amino acids adopt significantly different arrangements in the two clay systems. In the LDH case, the interlayer anions adopt a bilayer-like arrangement with their longest molecular axis oriented in an approximately perpendicular direction with respect to the hydroxide layers. In the montmorillonite (30) Sasai, R.; Shinya, N.; Shichi, N. T.; Takagi, T.; Gekko, K. Langmuir 1999, 15, 413.

Langmuir, Vol. 18, No. 7, 2002 2939

system, the interlayer cations prefer a monolayer arrangement with their longest molecular axis parallel to the layers. In all cases, the water molecules are generally excluded from the hydrophobic midplane formed by the phenyl rings of the anions. The different packing arrangements may be explained by the different charge density of the two clay systems, the greater charge density in the LDH leading to a more expanded arrangement in the direction perpendicular to the layers. The structural difference between the two amino acids studied (an additional OH group in the para position in (S)-tyrosine) has little influence on the adopted interlayer arrangement. Acknowledgment. We acknowledge Accelrys for their generosity in supporting this project, HEFCE for funding our SGI Onyx2 computer, and W. R. Grace & Co. - Conn. for funding S.P.N. Supporting Information Available: Table of atomic charges for amino acids. This material is available free of charge via the Internet at http://pubs.acs.org. LA0114528