5754
Ind. Eng. Chem. Res. 2007, 46, 5754-5765
Transferable Force Field for Water Adsorption in Cation-Exchanged Titanosilicates James P. Larentzos,*,† William F. Schneider,‡ and Edward J. Maginn‡ Geochemistry Department, Sandia National Laboratories, P. O. Box 5800 MS 0754, Albuquerque, New Mexico 87185-1454, and Department of Chemical and Biomolecular Engineering, UniVersity of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556-5637
A new, transferable classical force field potential for titanosilicate ion-exchange materials is reported. Planewave pseudopotential density functional theory (DFT) simulations are used to compute partial charges, where the atomic boundaries are partitioned through the Wigner-Seitz approach. The repulsion-dispersion interactions are modeled with a 12-6 Lennard-Jones potential. The force field is parametrized with the acidexchanged titanosilicate and transferred to the sodium-, cesium-, and strontium-exchanged titanosilicates and their niobium-substituted counterparts. Grand canonical Monte Carlo simulations are conducted to compute the sorptive properties of water in the acid-exchanged titanosilicate. Excellent agreement of the saturated water loading, positions, and occupancies with experimental neutron diffraction results is observed, while the adsorption energies at low water loading match DFT calculations. In transferring the force field to the sodium-, cesium-, and strontium-exchanged titanosilicates, canonical replica exchange Monte Carlo simulations are used to enhance cation sampling. The Wigner-Seitz force field parametrization effectively predicts the adsorption properties and provides a molecular-scale picture of the origins of titanosilicate selectivity in these materials. Overall, the methodology for translating electronic structural information derived from DFT calculations to classical force field based simulations can easily be extended to alternative crystalline materials such as zeolites and clays. 1. Introduction Processing of high-level nuclear waste requires efficient separation of cesium, strontium, and the actinides in order to minimize disposal volume. Crystalline inorganic ion exchangers are promising materials to perform these separations, because of their resistance to radiation degradation, their high thermal and chemical stabilities, and their high selectivity for targeted ions. Unfortunately, trace concentrations of cesium, strontium, and actinides in the presence of high sodium and hydroxide ion concentrations make the separation process through ion exchange especially challenging. Thus, there is a pressing need to develop and synthesize novel ion exchangers that exhibit a high selectivity toward targeted ions even in the presence of large concentrations of competitive ions and over a broad range of pH, while maintaining the material’s structural integrity under conditions of high radiation fields and extreme acidic or alkaline media. Titanosilicates are a potentially attractive class of ionexchange materials, especially with regard to nuclear waste remediation. In particular, synthetic titanosilicates with the mineral sitinakite topology1,2 of ideal formula M2Ti2O3SiO4‚ xH2O (M ) alkali metal cation) have been shown to be ideally suited for cesium and strontium uptake.3-5 There are numerous experimental studies characterizing the structure and performance of these materials.1,2,4,6-21 Previous experimental X-ray diffraction studies1 as well as molecular-simulation studies22 have attributed the high cesium selectivity to the fact that cesium fits nearly perfectly at the tunnel centers, where it is able to attain a highly stable, * Corresponding author. Phone: (505) 845-7486. Fax: (505) 8447354. E-mail:
[email protected]. † Sandia National Laboratories. ‡ University of Notre Dame.
8-coordinate complex with eight bonds to the silicate framework oxygen atoms. Further experimental studies have shown cesium selectivity can be enhanced with partial heteroatom substitutions in the framework, such as a 25% substitution of niobium for titanium, thus reducing the number of counterbalancing cations required for charge neutrality and creating more room for water to increase the cesium coordination from a maximum of 8 to a maximum of 12.9,15,17,19,22 Despite the enhancement in cesium selectivity with partial substitution of niobium into the framework, there is a decrease in selectivity for strontium that has been attributed to a reduction in the strontium coordination number from 10 in the non-niobium phase to 7 in the niobium phase.17,21 To better understand the origins of selectivity in the various cation-exchanged phases of the niobium and nonniobium titanosilicates, it is important to characterize the structural features of the extraframework cations and water molecules. Experimental studies have largely focused on X-ray diffraction to characterize these materials, which provide an averaged snapshot of cation and water structure in the channels, whereas molecular simulations have the advantage of providing instantaneous snapshots of cation and water structure. To this end, it is our goal to turn to molecular simulations to provide additional theoretical insight into cesium and strontium selectivity within the M2Ti2O3SiO4‚xH2O class of ion exchangers. Our long-term objective is to understand the nature of adsorption, diffusion, and ion-exchange selectivity in this class of materials. Performing these calculations in a fully quantum mechanical manner is infeasible, hence the need for accurate classical potentials. Classical molecular simulations have proven to be an extremely useful tool for examining zeolites and related crystalline porous materials, where much success has been achieved in accurately predicting properties such as cation and adsorbate structure, adsorption isotherms, isosteric heats of adsorption, adsorption energies, and adsorbate diffusion
10.1021/ie070276g CCC: $37.00 © 2007 American Chemical Society Published on Web 07/21/2007
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5755
coefficients.23-27 Accurate prediction of these properties through classical simulations requires the development and validation of quality force field potentials that describe the nature of the molecular interactions. A significant number of studies have been dedicated to the development of force fields, with the intention that they will be transferable over a wide range of structurally related materials. For titanosilicates, there are only a handful of classical simulation studies in which force fields have been developed, which have mainly focused on titanium silicate-1 (TS-1)28-34 and the Engelhard titanosilicates (ETS-4 and ETS-10).35-39 In a previous work, we developed a computationally efficient, fixed charge force field for the sitinakite based titanosilicates. The repulsion-dispersion interactions between water and the framework oxygen atoms were modeled with a Buckingham potential with parameters taken from the work of de Leeuw and Parker.40 Partial charges were empirically adjusted to reproduce experimentally determined water positions in the Nb-substituted titanosilicate, H0.5NaNb0.5Ti1.5O3SiO4‚2H2O. This force field was successful in accurately predicting preferred cation and water locations in the sodium- and cesium-exchanged phases of the non-niobium and niobium-substituted ion exchangers.22 To our knowledge, this is the only modeling study of the sitinakite topology titanosilicates. While this force field was successful for this limited class of materials, empirical adjustment of partial charges to match experimental data suffers from transferability and accuracy issues. Furthermore, the Buckingham potential is susceptible to instabilities at small atomic separations (the wellknown “Buckingham castastrophe”). For these reasons, we have developed a new force field based on the 12-6 Lennard-Jones potential with partial charges derived from density functional theory (DFT). Given the importance of Coulombic interactions in these systems, the determination of partial charges in this manner is essential for obtaining accurate results without recourse to previous experimental data. It also is a generalized methodology that can easily be applied to other crystalline materials, such as zeolites and clays. We begin by developing a potential model of the acid-exchanged titanosilicate, H2Ti2O3SiO4‚xH2O. We then use this model to compute isotherms and preferred adsorption sites of water and compare the results against neutron diffraction experiments.13 Finally, we extend the model to the sodium-, cesium-, and strontium-exchanged systems to test the force field, demonstrate its transferability, and provide a molecular-scale picture of the origins of titanosilicate selectivity for cesium and strontium. 2. Structural and Computational Details 2.1. Titanosilicate Structure. In each of the cation-exchanged materials, the titanosilicate framework consists of TiO4 cubanelike clusters that are bridged to each other in the ab plane by silicate groups and in the c direction by framework oxygen atoms. The negatively charged framework encloses onedimensional channels running along the c direction, with an accessible pore diameter of ∼3.5 Å, and is neutralized by countercations residing within the channels. The crystal structure of the acid-exchanged titanosilicate was determined through neutron diffraction studies,13 where it has been found to have a structural formula of H2Ti2O2SiO4‚1.5H2O. The unit cell is tetragonal, with lattice constants a ) b ) 11.0343 Å and c ) 11.8797 Å, space group P42/mbc, Z ) 8. Within each pore, an extensive hydrogen-bonding network is formed as water molecules are found experimentally to occupy two channel sites, denoted as O5 and O6 crystallographically. The hydrogen-bonding network acts to stabilize the framework
Figure 1. Crystallographic representation of the H-TS structure. The framework comprised Si (yellow), Ti (purple), and O (red). The O5 and O6 water sites are located within the channels and form a hydrogen-bonded network stabilizing the framework.
and maintain its structural integrity. An illustration of the ab plane for the acid-exchanged structure is given in Figure 1. When exchanging sodium into the acid-exchanged material, the crystal structure undergoes a phase change to P42/mcm, Z ) 4 symmetry with lattice constants of a ) b ) 7.8082 Å and c ) 11.9735 Å.1 The structural formula is reported as Na1.64H0.36Ti2O2SiO4‚1.8H2O. Powder diffraction data indicate that sodium ions occupy two sites, denoted as Na1 and Na2 crystallographically. The nonexchangeable Na1 site is located within the framework and between the silicate groups, whereas the exchangeable Na2 site resides within the channel, coordinating with four framework oxygen atoms and two water molecules. In addition, a significant amount of water is present in the unit cell and distributed over two sites, denoted as O(W1) and O(W2) crystallographically. In exchanging sodium for cesium, the space group symmetry remains P42/mcm, but the lattice constants expand slightly to a ) b ) 7.8258 Å and c ) 11.9815 Å.1 The structural formula is reported as Na1.49Cs0.2H0.31Ti2O3SiO4‚H2O, Z ) 4. Because of its size, cesium cannot fit into the framework positions occupied by sodium, making it impossible for complete ion exchange and leaving residual sodium ions in the Na1 framework sites as well as in the channels. In addition, the presence of sodium ions within the channels severely limits the titanosilicate uptake capacity for cesium. Powder diffraction experiments have shown cesium to partially occupy two sites, denoted as Cs1 and Cs2. Cs1 is located in the center of the tunnel and coordinates with eight framework silicate oxygen atoms. Cs2 is less prevalent but is stabilized by coordinating with four framework oxygen atoms and two water molecules. In partially substituting 25% of the titanium with niobium, the number of counterbalancing cations residing in the titanosilicate channels is reduced, resulting in an enhancement of selectivity for cesium. This creates more room for Cs1 cesium cations to coordinate with water, increasing the cesium coordination number up to 12. Cesium also occupies a less stable, 6-coordinated Cs2 site. In addition, water is observed to occupy two channel sites, the O(W1) and O(W2) sites. While the water sites have the same z-coordinate, the O(W1) adsorbs closer to
5756
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
the framework, whereas the O(W2) water site is located near the center of the channel and coordinates with Cs1 cesium. The stuctural composition is reported as Cs0.2H0.3NaNb0.5Ti1.5O3SiO4, with lattice constants of a ) b ) 7.8397(4) Å and c ) 12.0321(6) Å, space group P42/mcm, Z ) 4.17 Ion exchange of strontium into the acid-exchanged titanosilicate causes the crystal structure to undergo a phase change to space group symmetry P42/mmc, Z ) 8. The structural composition becomes Sr0.55H0.90Ti2O3SiO4‚2.2H2O.21 Crystallographically, strontium partially occupies two sites, denoted as Sr1 and Sr2. The Sr1 strontium are located in a highly stable, 10-coordinated environment, while the Sr2 strontium form either a 7- or 9-coordinated complex. The Sr1 and Sr2 sites are also occupied by water, which are denoted as the O(W1) and O(W2) sites. Overall, water populates six crystallographic sites, O(W1)O(W6), and coordinates heavily with the strontium cations. Exchange of strontium into the non-niobium and niobiumsubstituted framework causes a reversal in selectivity character as compared to the cesium-exchanged structures, where the nonniobium titanosilicate is more selective for removing strontium than the niobium-substituted structure. The coordination number of strontium is reduced to nine in the niobium phase where 16% of the titanium atoms are replaced by niobium and is further reduced to seven in the 25% niobium-substituted material. In the 16% substituted material, strontium occupies only one site, the Sr1 site, whereas water occupies two sites, the O(W1) and O(W2) sites, which act to stabilize the strontium sites. The structural composition is Sr0.49H0.70Ti1.68Nb0.32O3SiO4‚1.5H2O, with lattice constants of a ) b ) 7.8461(2) Å and c ) 11.9463(5) Å, space group P42/mcm, Z ) 4.21 2.2. Force Field Development. The main focus of this work is to develop a transferable force field that will accurately predict cation and water sorptive behavior in titanosilicate materials. We begin with the acid-exchanged titanosilicate (hereafter referred to as H-TS), since high-quality neutron diffraction experiments have pinpointed the locations of all protons and water molecules. After testing the model with H-TS, the force field is extended to the sodium- (Na-TS), cesium- (Cs-TS), and strontium-exchanged (Sr-TS) titanosilicates as well as the niobium-substituted cesium (Cs-NbTS) and strontium (Sr-NbTS) materials to demonstrate its transferability. For these systems, neutron diffraction experiments were not conducted. Thus, powder diffraction experiments were used, which makes structural characterization of extraframework species especially difficult. To develop a simple and robust analytical model, the interatomic energetic interactions are modeled with a pairwise additive short-range 12-6 Lennard-Jones potential and a longrange Coulomb potential, as given by the expression
Vij ) 4ij
[( ) ( ) ] σij rij
12
-
σij rij
6
+
qiqj 4π0rij
(1)
where rij is the separation distance between atoms i and j, σij represents the collision diameter, ij is the well depth, 0 is the permittivity of a vacuum, and qi is the charge of atom i. 2.2.1. Repulsion-Dispersion Parameters. For the shortranged van der Waals parameters, we have adopted the Kiselev convention by assuming that the framework oxygen atoms shield the effects of the framework silicon and titanium.41 The sodium, cesium, and strontium Lennard-Jones parameters are based on previous studies that were able to accurately compute each respective cation’s hydration free energy in aqueous solution.42,43 Water is treated with a rigid simple point charge (SPC) model.44
Table 1. 12-6 Lennard-Jones O-Owater Interaction Parameters Used in Previous Studies material
water model
σOw-Oframework (Å)
Ow-Oframework (K)
Zeolite-A62 Faujasite63 Zeolite-4A64,65 Zeolite-4A65 Montmorillonite46 Montmorillonite47 Silicalite54 Silicalite54 Silicalite54
TIP4P TIP4P SPC/E TIP3P TIP4P SPC/E TIP4P TIP5P MSPC/E
3.022 3.310 2.495 2.487 3.1536 3.166 3.077 3.06 3.06
78.02 70.05 282.10 279.29 77.08 78.18 85.41 86.78 83.46
Table 2. 12-6 Lennard-Jones Cross Parameters and Partial Charge for Sodium, Cesium, Strontium, and Water; O Represents Either the Framework or Water Oxygen Atom cation
σi-O (Å)
i-O (K)
qi
Oframe Owater Hwater sodium42 cesium42 strontium43
3.169 3.169 0.000 2.876 3.526 3.245
78.201 78.201 0.000 62.74 62.74 62.74
Wigner-Seitz -0.82 0.41 1.00 1.00 2.00
The SPC model is one of the simplest models that has been used to reproduce the thermodynamic phase behavior of water.45 Although the water model neglects polarization effects, it has been successfully used to match aqueous-phase hydration energies42 as well as the adsorption properties of zeolites and clays. We have chosen to use the SPC model given its simplicity and computational efficiency. The Lennard-Jones (LJ) parameters for the water-water interaction are taken from the SPC water model, thus leaving only the collision diameter, σij, and the well depth, ij, of the interaction between the water oxygen and framework oxygen atoms as the remaining adjustable parameters. There have been numerous classical simulation studies that have examined the sorptive behavior of simple molecules in porous materials.27 Despite the success of these studies, the van der Waals parameters tend to be limited in their transferability. For example, a comparison of selected water-framework LJ parameters used in previous studies of zeolites and clay materials is presented in Table 1, where it is clearly illustrated that Lennard-Jones parameters vary drastically depending on the material and water model used. Since it is not evident which of these LJ parameters to incorporate into our model, we have chosen to apply the same parameters as the pure-component oxygen in the SPC model to the water-framework interaction, justifiable by previous successful studies involving the clay montmorillonite.46,47 Thus, for the Owater-Owater and OframeworkOwater LJ parameters, we apply the SPC oxygen pure-component parameters. The Lennard-Jones force field parameters for all titanosilicate materials are summarized in Table 2. 2.2.2. Electrostatic Parameters. The H-TS partial charges that are used in the electrostatic contribution to the potential energy are derived through density functional theory (DFT) and transferred to Na-TS, Cs-TS, Cs-NbTS, Sr-TS, and SrNbTS. The electronic structure for one dehydrated unit cell of H-TS containing 96 atoms (8 H, 16 Si, 16 Ti, and 56 O) is determined by constraining the nuclear positions of all framework atoms (Si, Ti, O1, O2, O3, and O4) and all protons (H1) to the crystallographic coordinates. Periodic plane-wave DFT simulations are performed with the Vienna ab initio simulation package (VASP).48 Frozen-core electronic states are described using the accurate projector-augmented wave (PAW)49 approach, electron exchange and correlation are treated within the generalized gradient approximation (GGA),50 and plane waves are
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5757
included to a 400 eV energy cutoff.49 A 2 × 2 × 2 MonkhorstPack mesh is used to sample the first Brillouin zone. The unit cell parameters are computed by scaling the cell basis matrix parameter while constraining all nuclei to their fractional locations. The minimum energy cell parameters are all within 0.05 Å of the experiment, in agreement with the level of accuracy expected from the GGA. Additional DFT calculations are conducted on the unrelaxed experimental structure. Upon converging the electronic structure, the electron density is analyzed and partial charges are derived. To model the electrostatic potential, partial charges are placed at the atomic centers and assigned based on a partitioning of the electron density. However, given the arbitrary nature of defining atomic boundaries, there are many choices for the partition of the electron density. In this work, we have adopted a simple approach in which a spherical “Wigner-Seitz” radius is defined around each atomic nucleus. The number of electrons can be computed through integration of the electron density over the Wigner-Seitz sphere volume. The atomic charge can then be deduced by accounting for the nucleus valency. Hereafter, we refer to this as the Wigner-Seitz parametrization. For silicon and titanium, the Wigner-Seitz cutoff radii are taken to be equivalent to their covalent radii,51 which are 1.11 Å and 1.36 Å, respectively. To easily account for ion exchange without adjustment of the framework charges, all protons are constrained to have a +1 charge. Thus, the only undetermined Wigner-Seitz radii are those of the oxygen atoms in the system. To determine the oxygen Wigner-Seitz radius, a charge neutrality constraint is applied. From crystallography, there are four distinct framework oxygen sites in the asymmetric unit cell, which are labeled O1, O2, O3, and O4. Each of these oxygen sites have different bonding environments, and thus, we expect that they will have differing partial charges. The O1 and O3 oxygen sites are involved in bridging a Ti atom in a cubane cluster to a Si atom in the silicate group. The O2 oxygen site connects three Ti atoms together to form the cubane clusters. Finally, the O4 site bridges two Ti atoms of two different cubane clusters together along the c-axis of the unit cell.13 We assume that the Wigner-Seitz radii on each of these four oxygen sites must be equivalent, but because the oxygen sites have different bonding environments and, thus, different electron densities, they each have a different charge. Thus, the final step is to vary the oxygen Wigner-Seitz radius until charge neutrality is satisfied, where the Wigner-Seitz radius for oxygen was found to be 1.35 Å. The Wigner-Seitz radius is reasonable when considering that the optimized Wigner-Seitz radius is slightly greater than the ionic radius of oxygen, which is reported to be ∼1.28 Å.51 One potential criticism of the Wigner-Seitz methodology is the fact that the spheres do not completely fill space and can overlap. Thus, the total number of electrons contained within the spheres does not necessarily correspond with the total number of electrons in the system. Bader charge analysis52,53 is a popular alternative methodology that accounts for every electron in the system by identifying atomic boundaries with the minimum (electron) charge density surfaces. We find that the Bader approach coupled with the PAW has difficulty partitioning charge between oxygen and hydrogen because of the low electron density near the latter. The very reasonable agreement between the predicted and experimentally measured properties gives us confidence in the Wigner-Seitz derived charges. The charges of each distinct crystallographic site in H-TS as determined from the Wigner-Seitz methodology are shown
Table 3. Partial Charges for Each Titanosilicate Framework Atom. The Charges on the O2 (and O5 for Sr-TS) Do Not Include the Proton-Smearing Model Corrections atom
H-TS
Si Ti Ti/Nb O1 O2 O3 O4 O5 O6
2.02 1.07
Na-TS/Cs-TS 2.02 1.07
2.02 1.07
-0.89 -1.00 -0.89 -0.60
-0.89 -1.00
-0.89 -1.00 -0.60 -0.89 -1.00 -0.60
-0.60
Sr-TS
Cs-NbTS/Sr-NbTS 2.02 1.32/1.23 -0.89 -1.00 -0.60
in Table 3. Desbiens et al. have summarized the Si partial charges derived from ab initio computations for zeolite systems, where the charges range from +0.52 to +3.3 elementary charge units.54 The Wigner-Seitz parametrization yields a Si charge of +2.02 elementary charge units, which is within the range of previously accepted Si charges. In our previously derived empirical force field,22 our assumption that the Si charge is equivalent to the Ti charge was made based on the fact that their formal charges are equivalent. Clearly, the DFT parametrization shows that the Ti charge is significantly less than the silicon charge (by ∼1 elementary charge unit), suggesting that the electron density localizes over the Ti atoms to a greater extent than previously assumed. The Wigner-Seitz parametrization also shows that the most negative charge is on the O2 oxygen atoms contributing to the cubane-like clusters, followed by the O1/O3 silicate oxygen atoms and, finally, by the O4 cubane cluster bridging oxygen atoms. We note that the derived O1 and O3 partial charges are identical, which is reasonable given that they are structurally similar, where they each bridge the cubane-like clusters to the silicate groups to form an eight-membered ring. Furthermore, in Na-TS, Cs-TS, Cs-NbTS, and Sr-NbTS, the space group symmetry changes to P42/mcm, and the O1 and O3 sites become equivalent (denoted as the O1 site in these structures).1,17 The partial charges that are used in each cation-exchanged system are summarized in Table 3. We emphasize that, although the crystallographic labeling of the atom types has changed from one system to the next, the bonding environments are identical. Thus, the partial charges are transferred from the H-TS system to the Na-TS, Cs-TS, Cs-NbTS, Sr-TS, and Sr-NbTS systems without modification. In partially substituting niobium into the titanosilicate framework, powder diffraction techniques are unable to resolve their precise locations. Thus, the experiments suggest that niobium is randomly distributed over all octahedral sites. Since the niobium and titanium arrangement is uncertain, a “composite” atom is placed in all octahedrally coordinated sites that could be populated by either titanium or niobium. The partial charges for the composite Ti/Nb atom are determined by smearing the additional charge imposed by niobium over all the octahedral titanium sites in the cubane-like clusters. This approach of using a composite or average “T-atom” has been used previously with great success in modeling aluminum substitution for silicon in zeolites.55 2.3. Simulation Models. For the classical simulations, the unit cells of each cation-exchanged titanosilicate are expanded into nearly cubic-shaped supercells. A 2 × 2 × 2 supercell is used for H-TS and Sr-TS, while a 3 × 3 × 2 supercell is used for Na-TS, Cs-TS, Cs-NbTS, and Sr-NbTS. In general, cations and water molecules fractionally occupy crystallographic sites, and a single unit cell thus contains noninteger numbers of each. Given a sufficiently large supercell, the experimental
5758
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
compositions could be captured exactly, but at considerable computational expense. Thus, the compositions are approximated in the simulations as Na1.75H0.25Ti2O3SiO4‚xH2O (NaTS), Na1.50Cs0.25H0.25Ti2O3SiO4‚xH2O (Cs-TS), Cs0.25H0.25NaNb0.5Ti1.5O3SiO4‚xH2O (Cs-NbTS), Sr0.55H0.90Ti2O3SiO4‚ xH2O (Sr-TS), and Sr0.49H0.70Ti1.68Nb0.32O3SiO4‚xH2O (SrNbTS). An unknown water content, x, is given since this quantity will be computed through grand canonical Monte Carlo simulations. Finally, all framework atoms (Si, Ti, O1, O2, O3, and O4, Na1) are frozen at crystallographically determined positions, while the water molecules and exchangeable cations are permitted to move freely, without symmetry constraint, throughout the titanosilicate channels. For H-TS, the position of the acid-exchangeable proton position, H1, has been determined through neutron diffraction experiments. There are 16 H1 protons per unit cell (128 per supercell) to counterbalance the net negative charge on the framework. The H1 protons are associated with the framework O2 oxygen atoms and are involved in hydrogen bonding to the O5 water oxygen site. Thus, for the H-TS DFT calculations, we explicitly treat the exchangeable H1 protons and fix them to the experimentally determined positions. However, in the classical simulations, the protons must be treated in a classical sense without sacrificing the important physics of the system. There have been a few studies in which protons have been modeled in classical simulations. Auerbach and co-workers examined the acid phase of Zeolite Y, where they treated the proton explicity by reducing the charge to +0.2 elementary charge units and derived repulsion-dispersion Lennard-Jones parameters to mimic the short-range behavior.56 In our previous study with the sodium- and cesium-exchanged sitinakite-based titanosilicates,22 we chose an alternative approach. Since powder diffraction is unable to detect the proton positions in these materials, we assumed the protons delocalize over all framework oxygen atoms (O1, O2, O3, and O4 sites) in the system. Thus, in treating the protons in classical force field based simulations, they were incorporated into a united atom with the framework oxygen atoms, and the overall sum of the proton charge was essentially smeared over all framework oxygen atoms. We extend this same approach in our current study. However, in the H-TS material, neutron diffraction experiments indicate that the protons associate with the cubane-like O2 oxygen atoms. Thus, in each cation-exchanged material, a united atom site is used, where the proton charge is smeared over all O2 (O2 and O5 in Sr-TS) framework oxygen atoms. A charge of +1 elementary charge units is assumed for each proton. 2.4. Simulation Techniques. Classical grand canonical Monte Carlo (GCMC) simulations and canonical replica exchange Monte Carlo (REMC) simulations are conducted to simulate titanosilicate materials. All classical simulations are performed with a locally developed software package. Full periodic boundary conditions are employed with an interaction cutoff of half of the smallest simulation box length. The long-range Coulombic interactions are evaluated with the Ewald summation method. 2.4.1. Replica Exchange Monte Carlo Simulations. For the sodium-, cesium-, and strontium-exchanged materials, there exist rough energy landscapes exhibiting deep energy wells that can cause cations and water molecules to become trapped in local minima. It is difficult to properly sample configurations in these cases with conventional simulations. By heating the system to higher temperatures, the cations and water molecules are given enough thermal energy to overcome the deep energy wells and efficiently sample configuration space. In previous studies,22,57
we used a simulated annealing molecular dynamics protocol to overcome the energy barriers, where the system was heated to higher temperatures and systematically quenched to low temperatures. In this work, we utilize canonical ensemble (NVT) replica exchange Monte Carlo (REMC) simulations58,59 to efficiently probe configuration space. In this method, NVT Monte Carlo simulations are performed independently and simultaneously on multiple processors at various temperatures. Periodically, attempts are made to swap the configuration of the high-temperature processor with that of the low-temperature processor. A total of 16 independent realizations of the system are simulated in the NVT ensemble. Each realization is conducted at a different temperature, where Ti represents the simulated temperature on processor i. In replicating the system over 16 processors, the temperature ranges from 300 to 2440 K. Intermediate temperatures are chosen to ensure swap acceptance rates of at least 3% between adjacent processors. A Ti+1/Ti temperature ratio of 1.15 satisfied the acceptance rate criteria at low water loadings. At higher water loadings, the temperature ratio was reduced to 1.1, giving a temperature range between 300 and 1250 K. Finally, we emphasize that the hightemperature simulations are not necessarily physically reasonable in an experimental sense but are used to improve sampling of the configurational space. The predicted properties of the system are computed from the 300 K realization. To investigate the preferred adsorption sites and the influence of water on cation distribution, water is systematically added to the system followed by a sequence of replica exchange Monte Carlo simulations. To begin, the exchangeable cations are randomly inserted into the dehydrated structure. Canonical Monte Carlo simulations with periodic replica exchange attempts every 100 Monte Carlo steps are performed for a total of 30 million Monte Carlo steps to equilibrate the dehydrated system. Statistics are collected every 5000 Monte Carlo steps and analyzed to determine the preferred cation distribution. Next, a predetermined amount of water is randomly introduced into the system, and another 30 million Monte Carlo steps with replica exchanges are performed to determine the preferred cation and water distribution. The process is repeated until Nmax water molecules per titanosilicate unit cell are attained, where Nmax is the approximate number of water molecules observed experimentally at saturation. Thus, the actual water loadings that are simulated range from 0 to Nmax, with M intermediate water loadings. For each of the M intermediate steps, a replica exchange Monte Carlo simulation is conducted and the cation and water distributions are examined. 2.4.2. Grand Canonical Monte Carlo Simulations. Grand canonical Monte Carlo (GCMC) simulations are used to examine the adsorption behavior and the extent of water loading in the titanosilicate materials. The GCMC algorithm is useful for adsorption studies since the chemical potential, volume, and temperature are fixed, while their conjugate variables (number of particles, pressure, and energy) fluctuate about their mean value at equilibrium. The chemical potential can be related to the fugacity, where the equilibrium condition is satisfied when the fugacities of the gas phase and the adsorbed phase are equivalent. Four Monte Carlo moves are considered: translations, rotations, insertions, and deletions. Translation and rotation moves are attempted with equal probability 90% of the time, while the insertion and deletion moves are each attempted 5% of the time. To improve insertion efficiency and speed the convergence, a configurational biasing algorithm was imple-
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5759 Table 4. Calculated Water Positions in H-TS, with Experimental Sites Given in Parantheses; The Shifted Coordinates are Reported beneath the Raw Averaged Coordinates atom
X
Y
Z
dev. (Å)
std. dev. (Å)
pop. frac.
exptl. occ.
water (O5)
0.307 0.307 0.502 0.500
0.052 0.052 0.002 0.000
0.499 0.500 0.099 0.099
0.37 0.37 0.23 0.23
0.23
95.5%
100.0%
0.64
46.7%
50.0%
water (O6)
Figure 2. H-TS adsorption isotherm. Error bars are determined by block averaging over simulation data sets for each state point.
entire adsorption isotherm is computed to ensure sufficient sampling. Error bars are then obtained by block averaging over the M simulation data sets for a given state point. All GCMC simulations are conducted for 10 million equilibration steps, followed by 10 million production steps where averages are computed. Finally, to compute the adsorption energy, an additional 10 million NVT Monte Carlo steps are conducted on the final equilibrated configuration from the GCMC simulations. 3. Results and Discussion
Figure 3. Water adsorption energy as a function of loading for H-TS.
Figure 4. Population fraction of the experimental Na-TS sites as a function of water loading.
mented,60 where the water molecule is inserted through a combination of energy and orientational biasing. One problem with using GCMC for systems of this type is that cations are often trapped in deep energy minima; this makes thermal equilibration moves difficult to carry out for the cations, leading to poor sampling. Thus, initial configurations for the GCMC simulations are obtained from the final configurations of the REMC simulations conducted at M different water loadings. For each of the M different initial configurations, the
3.1. Force Field Development: The Acid-Exchanged Titanosilicate. We test the Wigner-Seitz parametrized force field by comparing water adsorption properties of H-TS computed through GCMC simulations to experimental neutron diffraction data. Seven independent simulations are conducted over a range of state points to obtain an adsorption isotherm for the acid-exchanged system. The saturation pressure of water at 300 K is ∼3.56 kPa. The adsorption isotherm is computed by varying the ratio of the gas-phase pressure to the saturation pressure from 10-6 to 1 at 300 K, increasing the pressure by an order of magnitude at each intermediate state point. At saturation (P/Psat ) 1.0), the parametrization predicts 12.00 ( 0.78 water molecules per unit cell. There is outstanding agreement with the neutron diffraction measurement of 12 water molecules per unit cell, which happens to be the only experimental data point on the isotherm to compare against. The adsorption isotherm is shown in the Figure 2. In characterizing the positions and relative occupancies of the water molecules, an averaged coordinate for the water sites at saturation is computed and compared with the asymmetric unit cell water coordinates determined through neutron diffraction in Table 4. In addition, an idealized coordinate for each site is reported, where the raw averaged coordinate is shifted onto a nearby site with higher point group symmetry. A complete description of the methodology used to compute the averaged position and occupancy of each site is given in a previous work.57 Overall, excellent agreement is obtained with the WignerSeitz parametrization, where the oxygen atoms of the water molecules deviate from the experimental results by 0.37 Å and 0.23 Å for the O5 and O6 sites, respectively. The standard deviation is also computed by measuring the fluctuations about the averaged positions. The magnitude of this fluctuation gives an indication of the relative mobility of water and the strength of adsorption. On the basis of the standard deviation, the GCMC simulations indicate that the O5 water is more tightly bound and, thus, less mobile than the O6 water molecules. In addition, a population fraction is computed for each site, which is essentially the fraction of time the sites are occupied with a water molecule over the entire simulation trajectory. The population fractions for the O5 and O6 sites are 95% and 47%, respectively, which compare well with the reported experimental occupancies of 100% and 50%, respectively.
5760
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Figure 5. Population fraction of the experimental (a) Cs-TS and (b) Cs-NbTS sites as a function of water loading.
The high-loading region of the isotherm is accurately predicted with the force field parametrization. This is not surprising, given that, at high pressures, the adsorption of water is insensitive to the partial charges used and depends strongly on the collision diameter parameter, σ, between the water molecule and the framework oxygen atoms. To rigorously test the potential parameters, the low-loading region of the isotherm must be accurately predicted before the force field can be used with confidence. Unfortunately, to our knowledge, there have been no experimental studies of the Henry’s law region of the adsorption isotherm, nor have there been any studies of the isosteric heat of adsorption to give some indication of the binding strength as water is adsorbed. Thus, we direct our attention to DFT calculations to provide a comparison with the low-loading region of the adsorption isotherm. We expect DFT within the GGA to provide an accurate description of the water binding energy. The GGA captures the major contributions to chemical bonding, including contributions such as polarization not explicitly captured in the classical model. Purely nonbonded dispersion interactions are not correctly predicted but are also not a major component of water-oxide bonding. For example, the GGA is known to predict the binding energy of water to oxide materials to within a few kcal/mol of the experiment.61 Thus, we expect that the binding energy of water to the titanosilicate framework will be a useful metric to test the force field at low water loadings. By demonstrating that classical simulations provide internally consistent binding energies as DFT within the GGA, the “effective” potential produced through the Lennard-Jones and Coulombic interactions can be shown to be suitable in modeling all the chemically significant interactions (including electrostatic, polarization, and covalency effects) captured in the GGA. Adsorption energies for a given configuration ri at a water loading of 1 molecule/unit cell are calculated through the following relation corresponding to n water molecules
∆Eads(ri) ) EH2O+TS(ri) - (ETS + nEHIG2O)
(2)
where EH2O+TS(ri) is the energy of the hydrated titanosilicate with nuclear coordinates ri, ETS is the energy of the dehydrated titanosilicate, and EHIG2O is the energy of a gas-phase water molecule behaving ideally. In this equation, the dehydrated unit cell and a single water molecule in the
ideal gas state are the reference systems. The average adsorption energy is then computed through the expression
〈∆Eads(r)〉 )
1
N
∑ ∆Eads(ri) N i)1
(3)
where N is the number of water configurations. For the hydrated system, 10 different water configurations are investigated, where each water molecule is located near the O5 adsorption site. These water configurations are obtained by running a classical canonical Monte Carlo simulation at 300 K with a water loading of 1 water molecule/unit cell. Through DFT simulations, the average adsorption energy is found to be -56.9 ( 6.3 kJ/mol. To test the force field, the classical static energy is calculated on a reference system that comprised a 2 × 2 × 2 supercell of the dehydrated acid-exchanged structure. For the hydrated system, the same configurations of O5 water molecules as for the above DFT calculations are used and replicated into the supercell. The classical energy calculations indicate an average adsorption energy of -57.9 ( 2.4 kJ/mol. The Wigner-Seitz parametrization agrees well with the DFT simulations and indicates that the low-loading behavior is correctly captured with the proposed force field. The adsorption energy obtained from classical molecular dynamics simulations as a function of water loading is shown in Figure 3. The water-framework and water-water contributions to the total adsorption energy are also shown. At low water loading, the adsorption energy remains essentially constant at approximately -50 kJ/mol until a loading of nearly 8 water molecules/unit cell, at which point the O5 adsorption sites become fully populated. Note that the adsorption energy is ∼5 kJ/mol higher (less negative and weaker) than the adsorption energy calculated from DFT, which is likely due to the fact that the adsorption energy is determined at 300 K from Monte Carlo simulations, whereas the DFT adsorption energy is determined at a “frozen” state with energy minimizations. Through classical simulations, the adsorption energy at low loading is found to be entirely due to the water-framework interactions, where the O5 water molecules accept hydrogen bonds from the neighboring O2-H1 united atom site. As the pressure is increased, the O6 water site begins to populate, where they accept hydrogen bonds from the neighboring O5 water molecules. As a result, the water-water adsorption energy decreases at the expense of a slight increase in waterframework energy. Overall, the total adsorption energy decreases with the formation of the hydrogen-bonding network, and an
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5761 Table 5. Simulated Sodium and Water Sites in Na-TS; The Coordinates have been Shifted onto Nearby, High-Symmetry Sites atom sodium (Na2) sodium (OW2) water (OW1) water (Na2) water (OW2)
X 0.415 0.407 0.270 0.427 0.473
Y
Z
Na-TS 0.415 0.074 0.407 0.301 0.270 0.500 0.427 0.054 0.473 0.293
dev. (Å)
std. dev. (Å)
0.28 0.29 0.06 0.11 0.50
0.43 0.45 0.30 0.69 0.53
overall interaction energy of approximately -68 kJ/mol is observed at saturation. Since the Wigner-Seitz parametrization accurately predicts the water loading at saturation as well as the adsorption energy at low loading, we believe this is a well-suited parametrization to compute properties of the titanosilicate materials. Thus, it will be used for all subsequent simulations on these materials. We anticipate that future experiments will provide additional insight into the adsorption behavior in the Henry’s law region to rigorously test the force field. 3.2. Transferability of the Force Field. 3.2.1. Na-TS. Replica exchange simulations are conducted at various hydration levels of Na-TS, ranging from 0 to 7 water molecules/unit cell. To investigate the effect of water adsorption on cation distribution, a plot of the population fraction of the Na-TS sites as a function of water loading is shown in Figure 4. In the dehydrated structure, the exchangeable sodium cations are observed to occupy the experimental O(W1) and O(W2) sites, where each site is ∼25% occupied. We refer to these sodium cations as the O(W1) and O(W2) sodium cations. In addition, since the framework sodium cations are fixed, the Na1 site is fully occupied with sodium. As water is slowly introduced into the system, a redistribution of the sodium cations is observed. Initially, water begins to fill the O(W1) site, the preferred water adsorption site. Meanwhile, the O(W1) sodium cations are forced to completely depopulate from the O(W1) site and redistribute to the Na2 site. The O(W2) sodium cations remain in place up to a water loading of 3 water molecules/unit cell. At loadings above 3 water molecules/unit cell, water fully occupies the O(W1) site. At this stage, the O(W2) site becomes the preferred water adsorption site and begins to populate, causing the O(W2) sodium cations to redistribute to the Na2 site. Experimentally, the water loading at saturation is reported as 7.4 water molecules/unit cell. Near saturation (7 water molecules/unit cell), water is predicted to fully occupy the O(W1) site and partially occupy the disordered O(W2) and Na2 sites. Similarly, the exchangeable sodium cations partially occupy the Na2 and O(W2) sites. This behavior where sodium and water adsorb to both the Na2 and O(W2) sites was observed in previous simulation studies22 and justified experimentally because of similar X-ray scattering lengths for sodium and water. Finally, the adsorption isotherm of water in Na-TS is computed through grand canonical Monte Carlo simulations. At saturation, a water loading of 6.6 ( 0.8 water molecules is computed, which agrees (within statistical error) with the experimentally measured loading of 7.4 water molecules/unit cell. The predicted sodium and water positions are given in Table 5. Overall, the Wigner-Seitz parametrization is in better agreement with the experiment compared to a previously used empirical force field for Na-TS.22 3.2.2. Cs-TS and Cs-NbTS. The site occupancies as a function of water loading for Cs-TS and Cs-NbTS are illustrated in Figure 5. In the dehydrated Cs-TS material,
Figure 6. Side view along the a-axis of Cs-NbTS. When the cesium cation is located at the tunnel center with fractional coordinates of (1/2, 1/2, 1/4), there is a vacant region of the titanosilicate at fractional coordinates (1/2, 1/2, 3/4). It is this region where the simulations predict additional water adsorption.
cesium cations are predicted to predominantly locate in the Cs1 site (population fraction of 45%), indicating the high stability of this site due to the high coordination environment. In the Cs1 site, cesium coordinates with eight framework oxygen atoms. As water adsorbs, the O(W1) site rapidly populates with water and approaches full occupancy at a water loading of 4 water molecules/unit cell. In addition, there is a non-negligible amount of water locating in the Cs1 and Cs2 sites at a water loading of 4 per unit cell. As the water populates these sites, cesium steadily depopulates from the Cs1 site and begins to prefer the 6-coordinated Cs2 sites. Interestingly, in computing the adsorption isotherm of CsTS, there is an indication that there is, in fact, a higher water loading than previously thought through powder diffraction experiments. A water loading of 5.8 ( 0.3 water molecules/ unit cell is observed through grand canonical Monte Carlo simulations, which is significantly higher than the experimental estimate of 4 water molecules/unit cell. At 4 water molecules/ unit cell, the cesium population fractions are approximately 35% and 8% for the Cs1 and Cs2 sites, respectively. Experimentally, the occupancies of these sites are 20% and 10%, respectively. At a water loading of 6 water molecules/unit cell, the simulated population fractions are observed to be 29% and 10%, respectively, which is in better agreement with the experiment. Although we cannot compare these structures directly since there is slightly more cesium in the simulated structure as compared to the experimental structure, we do approach the experimental occupancies as water loading increases to 6 water molecules/ unit cell. In addition, we note that Cs-TS is essentially the same as Na-TS, where one sodium cation is ion-exchanged for cesium. Thus, we expect to have a water loading similar to that for the Na-TS structure, but with slightly less water because of the larger ionic radius of cesium. Given the difficulties of X-ray diffraction experiments to detect sodium in the presence of cesium in the Cs-TS structure, we believe that the experimentally determined water loading is too low. Furthermore, we have previously demonstrated that thermogravimetric analysis (TGA) in conjunction with powder diffraction experiments is necessary to provide accurate water loadings.57 Thus, future TGA experiments may be helpful in reconciling the differences between simulation and experiment. For the saturated structure determined through grand canonical Monte Carlo simulations, we report the simulated cesium cation and water coordinates in Table 6 for the Cs-TS and Cs-NbTS structures to give an indication of the accuracy of the force field.
5762
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Table 6. Simulated Cesium and Water Sites in Cs-TS and Cs-NbTS; The Coordinates Have Been Shifted onto Nearby, High-Symmetry Sites atom
X
Y
Z
cesium (Cs1) cesium (Cs2) water (OW1)
0.500 0.500 0.281
Cs-TS 0.500 0.255 0.500 0.164 0.281 0.500
cesium (Cs1) water (OW1) water (OW2)
0.500 0.293 0.584
Cs-NbTS 0.500 0.254 0.293 0.500 0.584 0.468
dev. (Å)
std. dev. (Å)
0.06 0.40 0.08
0.41 0.26 0.40
0.05 0.12 0.40
0.27 0.46 0.69
In partially substituting 25% of the titanium crystallographic positions with niobium, the number of counterbalancing cations residing in the tunnels is reduced. This allows Cs1 cesium cations to coordinate with additional water molecules besides the eight framework oxygen atoms. Thus, the increased selectivity for cesium is explained by the higher coordination environment, which is reported experimentally to be between 8 and 12. As shown in Figure 5b, essentially all of the cesium prefers the Cs1 site in dehydrated Cs-NbTS, where it binds to eight framework oxygen atoms in the energetically favorable (1/2, 1/2, 1/4) position. As water begins to adsorb, the O(W1) adsorption site rapidly fills with water. At a loading of 4 water molecules/unit cell, the O(W1) site is near full occupation and water begins to fill near the experimental Cs2 site, where a favorable hydrogen-bonding network is created with the existing water molecules in the O(W1) site. Cesium remains in the Cs1 site and is unaffected by the increased water loading, especially at loadings as high as 6 water molecules/unit cell. Experimentally, 4 water molecules/unit cell adsorb to the framework structure. However, it is critical to note that the experimental water loading for Cs-NbTS was determined by comparing the unit cell volume increase with Cs-TS. Since Cs-TS and Cs-NbTS have similar compositions, the volume increase was attributed solely to the 25% niobium substitution, indicating that the hydration level remains the same. However, as mentioned above for Cs-TS, we have reason to believe that the experimentally measured water loading is too low, especially since TGA experiments have not been conducted. Interestingly, grand canonical Monte Carlo simulations predict a water loading of 5.6 ( 0.2 water molecules/unit cell at saturation, suggesting that the hydration level remains constant with substitution of niobium into the titanosilicate framework. Furthermore, this discrepancy with experiment can be explained by re-examining the acid-exchanged titanosilicate structure, presented above. In the acid-exchanged structure, there is an absence of large countercations in the tunnels. Thus, water adsorbs to two sites. The first water site hydrogen bonds to the framework oxygen atoms of the cubane-like structure, while the second water site creates a highly stable hydrogen-bonding network with the framework water molecules. Now considering the Cs-NbTS structure, there is only one cesium cation and one proton present in each unit cell, based on the experimentally determined composition. The protons are known to bind to the cubane framework O2 oxygen atoms. Since there are 2 equivalent Cs1 sites located at the center of the tunnels and at c ) 1/4 and 3/4, on average, the Cs1 site is half occupied with cesium. Thus, half of the tunnel will contain a vacancy in the Cs1 site, while the other half is occupied with cesium, as illustrated in Figure 6. The vacancy section of the tunnel has an environment that mimics the acid-exchanged phase, where water can easily form a hydrogen-bonding network. In fact, in the absence of cesium,
electron density is observed near this site,17 which suggests that the simulation results may be reasonable. From an experimental point of view, the increase in cesium selectivity with niobium substitution has been attributed to an increase in coordination number up to 12, where water stabilizes cesium. However, it is extremely important to consider that X-ray diffraction techniques are limited to providing the averaged atomic positions. Since the O(W1) and O(W2) sites are each partially occupied, the instantaneous coordination number is reported to vary between 8 and 12. To obtain a coordination number of 12 (8 framework, 4 water), all of the water must migrate out of the O(W1) site and into the nearby O(W2) site. Alternatively, a coordination number of 8 (8 framework, 0 water) is obtained when all of the water is located in the O(W1) site. While it is certainly possible to have instantaneous coordination numbers of 12 (or 8) depending on the positioning of the niobium in the octahedral sites, the statistically averaged coordination number based on partial occupancies of the O(W1) and O(W2) sites is 10. Through simulations, we observe a cesium coordination number of 8 from the framework oxygen atoms but only observe small contributions (
