J. Phys. Chem. 1988,92, 867-871 pancreatic ribonuclease A. The resulting L values were inconclusive, but the larger T values suggested collective motions of regions larger than single residues, in accord with molecular dynamics calculations. A similar refinement of an oligonucleotidelM(with 12 nucleoside units) considerably reduced the noise level of the difference Fourier map, allowing 15 additional low-occupancy water positions to be identified. Advances in measurement techniques and procedures for restrained refinement'@ will doubtless soon produce a ratio of observations to adjustable parameters large enough to permit a determination of anisotropy for at least some of the atoms in other macromolecules. One strategy might be fmt to refine all the atoms isotropically and then to select judiciously perhaps 50 or 100 atoms with large Ui, for anisotropic refinement. When Ui, is small, the pdf cannot be very anisotropic (in the absolute sense) since Ui,is an average of three positive numbers. On the other hand, for pdf s with very large second moments the Gaussian formalism is likely to be quite inadequate. More likely, the pdf would have (109) Konnert, J. H.; Hendrickson, W. A. Acta Crysrallogr., Secf.A 1980, 36, 344-350.
867
to be modeled as a superposition of two or more pdPs, requiring location of their centroids and weights. Unfortunately, the larger the rms displacements, the more rapidly the scattering from the center in question falls off with scattering angle. It is just when the anharmonic terms become important that they may become impossible to determine. There seems no way out of the dilemma. In any event, the next decade will surely see more use of crystallographic information in the field of large-amplitude molecular motions, a field that has traditionally been reserved for molecular spectroscopyl10and gas-phase electron diffraction.
Acknowledgment. This work has been supported by the Swiss National Science Foundation, and owes much to discussions held while J.D.D. and V.S. were enjoying guest privileges at the Chemistry Laboratories of the California Institute of Technology, Pasadena, CA. J.D.D. is grateful for the award of a Sherman Fairchild Scholarship w%ich made his stay at Caltech possible. We have benefited from many discussions with Professors Hans-Beat Biirgi and Emily Maverick. (110) Ito, M. J. Phys. Chem. 1987.91, 517-526.
ARTICLES Molecular Dynamics Studies on Zeolites. 3. Dehydrated Zeolite A P. Demontis, G. B. Suffritti, Istituto di Chimica Fisica. Uniuersith di Sassari, I-071 00 Sassari, Italy
S. Quartieri, Istituto di Mineralogia e Petrologia, Universith di Modena. I-41 100 Modena, Italy
E. S. Foiq and A. Gamba* Dipartimento di Chimica Fisica ed Elettrochimica, Universith di Milano, I-201 33 Milano, Italy (Received: November 13. 1986)
A model potential in two forms (harmonic and anharmonic) is proposed to be used in molecular dynamics simulations of silicate frameworks. This model is applied to the calculation of structural and vibrational properties of the anhydrous phase of Linde Zeolite 4A. Our system is formed by a cubic box corresponding to one crystallographic cell containing 662 atoms without constraints. The results are compared with experimental data. The proposed models satisfactorily reproduce the main features of the aluminosilicate framework structure and dynamics.
Introduction Recently, a number of theoretical studies on zeolites have been carried out in order to explain the physical and chemical properties of several In particular, attention was devoted to (1) Dempey, E. J . Phys. Chem. 1969, 73, 3660. (2) Mortier, W. J. J. Phys. Chem. 1975, 79, 1447. (3) Cohen de Lara, E.; Tan, T. N. J. Phys. Chem. 1976,80, 1917. (4) Ogawa, K.; Nitta, M.; Aomura, K. J. Phys. Chem. 1978, 82, 1665. (5) No, K. T.; Jhon, M. S . J. Korean Chem. Soc. 1979.23, 374. (6) No, K. T.; Chon, H.; Ree, T.; Jhon, M. S. J. Phys. Chem. 1981,85, 2065. (7) Texter, J. Phys. Reu. 1981, 23, 4407. (8) Cheetham, A. K.; Eddy, M. M.; Jefferson, D. A.; Thomas, J. M. Nature (London) 1982, 229, 24.
statistical methods such as Monte Carlo and molecular dynamics (MD) techniques for the determination of the microscopic properties of zeolite^.^^-^^ (9) Gellens, L. R.; Mortier, W. J.; Lissillour, R.; Le Beuze, A. J. J. Phys. Chem. 1982,86,2509. (10) Klinowski, J.; Ramdas, S.; Thomas, J. M.; Fyfe, C. A.; Hartman, J. S . J. Chem. Soc., Faraday Trans. 2 1982, 78, 1025. (11) No, T. K.; Bae, D. H.; Jhon, M. S. J. Phys. Chem. 1986,90, 1772. (12) Sauer, J.; Hobza, P.; Zahradni'k, R. J. Phys. Chem. 1980.81.3318. (13) Beran, S. Chem. Phys. Lerr. 1981,84, 111. (14) Beran, S . Chem. Phys. Left. 1982, 92, 86. (15) Soto, J. L.; Myers, A. L. Mol. Phys. 1981, 42, 971. (16) Soukoulis, C. M. J. Phys. Chem. 1984,88,4898. (17) Demontis, P.; Suffritti, G. B.; Alberti, A.; Quartieri, S.; Fois, E. S.; Gamba, A. Gazz. Chim. Iral. 1986, 116, 459.
0022-3654/88/2092-0867$01 .50/0 0 1988 American Chemical Society
868 The Journal of Physical Chemistry, Vol. 92, No. 4, 1988
TABLE I: Potential Constants for Silicate Frameworks (Harmonic Model)
bond Si-0 AI-0
0-(Si)-0 0-(AI)-0
Na-0
k, kcal mol-’A-2 500.0 500.0 103.0 103.0 30.0
rn. A 1.605 1.760 2.61786 2.87070 2.385 if rsxp< 2.5 2.614 if relp > 2.5
In a previous workla we proposed a model potential used in two different approximations, one being a simple force constant potential and the other consisting of a short-range atom-atom potential, in line with a model used by Clarke et al.20 for MD simulations of ionic fluids. The common feature of both approaches is the finite range of the atom-atom potentials, so that a small and possibly predetermined number of interactions per particle is considered in the simulation, and no truncation problems arise. Consequently, the computer time is greatly reduced and accuracy is gained. In spite of the difficulty, stated in ref 1 1 , in obtaining potentials suitable for the representation of the zeolite framework, our model satisfactorily reproduced the main features of aluminosilicate framework structure and dynamics when applied to anhydrous and hydrated The purpose of this work was to perform MD simulations of dehydrated Linde Zeolite 4A in order to verify the ability of our model to match the structure and the spectroscopic properties of zeolites with various topologies. Method and Models In the atom-atom approximation, the potential energy of a system of atoms can be expanded in a power series around equilibrium distances of the single-potential functions between particles k! where Uoij = U(roij)is the value of the atom-atom potential function a t equilibrium distance roijand
This expansion is valid for any analytically regular potential function and can include both short- and long-range effects. If only the second-order terms are retained and Uoljare neglected, because they have no influence on the dynamics of the system, the usual force constant model is obtained. This model has been widely adopted for molecules and crystals when, for instance, harmonic vibration frequencies had to be computed. For silicates, the following assumptions were made: (a) The potentials for A1-0, Si-0, 0-0, M-O (M = Na+, Ca2+, or any cation present in the crystal) interactions were represented by
U(r) = y2k(r - ro)2
(3)
The values of k and r, adopted for each interaction are reported in Table I. In particular, the values for k were derived from spectroscopic data of Natrolitez1-22 taking into account also some common features of vibrational spectra of zeolites,23while the
Demontis et al. equilibrium bond lengths ro were deduced from structural data.24 The Na-0 distances in anhydrous Zeolite 4A can be split into two groups centered at about 2.35 and 2.60 A, r e s p e c t i ~ e l yso ,~~ two different equilibrium distances were used. This ad hoc assignment was made necessary because of the peculiar coordination of Na(2) atoms which are linked, by very different bond lengths, to only three coplanar oxygen atoms (see Figure 3 , ref 25). (b) N o other possible contacts were considered. (c) The initial topology of the framework bonds was retained during the M D simulation. (d) Only first neighbors were considered as interacting atoms. In spite of the simplicity of this model, the simulations reported below show that it can satisfactorily describe the main features of the framework structure and dynamics. Evidently, bond making and breaking, as well as diffusion of the framework atoms, cannot be described, so that phase transitions and defect dynamic simulations are beyond the limits of the model, but at ordinary temperatures these phenomena are neglibible in silicate frameworks. On the other hand, this potential form allows a great reduction in computing time and simulations of systems with a very large number of particles. It is to be remembered, in fact, that just one crystallographic cell of a zeolite may contain about a thousand atoms. If a more realistic description of the potential is required, terms of order higher than two in eq 1 have to be considered. For instance, by including terms up to the third order, eq 2 becomes
U(r) = Uo+ ‘/’k(r - ro)’ + A / 6 ( r - ro)3
(4)
and, assuming Uoand A to be negative, they can be adjusted so that Uorepresents the bond energy and the relative maximum of the function at r = (ro- 2 k / A ) > ro is zero. This happens for A = -(-2k3/3Uo)’/’. If the evaluation of this function is truncated at r,, = ro - 2 k / A where the first derivative is zero, eq 4 approximates a realistic interaction allowing anharmonic vibration, bond breaking, and diffusion. However, in this approximation, rmaxis determined by the values of k and Uo.In order to make rmaxan adjustable parameter, the expansion in eq 1 has to be extended at least up to the fourth order. Nevertheless, the lack of the first-order term makes the resulting function too symmetrical with respect to the minimum, so that a fifth-order expansion is preferable. For the complete determination of the involved coefficients, five constraints must be imposed; in this work they were the following: (ii) the value (i) the value of the function in the minimum (Uo); = 0; (iiii) of r-, where the function must be zero; (iii) U’(rmax) k assigned from spectroscopic data; (v) an arbitrary value rl < ro where the function is zero. We set rl = 2ro - rmaxin order to get simple forms for the expansion coefficients. In particular, the fifth-order expansion was represented by
U(r) = Uo+ y2k(r -
+ A / 6 ( r - r0)3+ B/24(r - ro)4+ C/12O(r - ro)5 ( 5 )
with
A = --
+ 0.25k
1
(18) Demontis, P.; Suffritti, G. E.;Quartieri, S.; Fois, E. S.; Gamba, A. In Dynamics of Molecular Crystals; Laswmbe, J., Ed.; Elsevier: Amsterdam, 1987; p 699. (19) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. Zeolites, in press. (20) Clarke, J. H. R.; Smith, W.; Woodcock, L. W. J. Chem. Phys. 1986, 84, 2290. (21) Pechar, F.; Rykl,D. Infrared Spectra of Natural Zeolites; Academia
Nakladatelstoi Czlkoslovenske’ Akademie Ved: Praha, 1985. (22) Pechar, F.; Gregora, I.; Rykl, D.Collect. Czech. Chem. Commun. 1981,46, 3043.
(23) Flanigen, M. Zeolite Chemistry and Catalysis; ACS Monograph 171; Rabo, J. A., Ed.; American Chemical Society: Washington, DC, 1976; Chapter 21. (24) Alberti, A., University of Sassari, private communication. (25) Pluth, J. J.; Smith, J. V. J . Am. Chem. Soc. 1980, 102, 4704.
The Journal of Physical Chemistry, Vol. 92, No. 4, 1988 869
Molecular Dynamics Studies on Zeolites
TABLE U: Potential Constants for Silicate Frameworks (Anharmonic Model)" bond uo k A B Si4 -95.0 500.0 -936.586 AI-0 0-(Si)-0 0-(A1)-0 Na-0
-85.0 -5.844 -7.027 -25.0
465.0 100.0 100.0
100.0
"Energies are in kcal mol-' and distances in
-888.023 -337.758 -308.024 630.86
Lmax
VO > 0.25k - ro)21
rm
-29863.2 1
2.65 2.8073 3.21 3.52 3.10
< rl, (6)
The models using eq 3 or eq 4 and 5 will be referred to below as harmonic and anharmonic models, respectively. By use of an anharmonic model for silicates, the above assumptions (a)-(d) were modified in the following way: (a') Quation 4 was adopted for Ala,Si-0, and 0-0 contacts, while for M - O (M = Na+, Ca2+,etc.) eq ( 5 ) was used. The values of the selected parameters are reported in Table 11. In this case the flatness of the potential near the minimum made the introduction of two different values of ro for Na-0 contacts unnecessary. (b') The same as (b). ( d )The framework bond topology may be checked every n step a t will. In order to get a direct comparison with harmonic model, the initial topology was left unchanged in this work. (d') Atoms were considered interacting for r Irmx. In Tables I and I1 some constants are slightly different from those used in ref 18 in order to get better IR frequencies, but the new values are considered valid for all zeolites; indeed, they were adopted also for hydrated Natrolite in ref 19. In the present M D calculations all the atoms were left free to move without geometry and symmetry constraints under the action of the potential described by the harmonic or anharmonic model. The periodicity of the crystals was simulated by the usual minimum image convention, and the equations of motion were integrated by means of a modified Ver1et"salgorithm.% At each time step the atomic coordinates were referred to the asymmetric unit by symmetry group operations and stored in order to obtain their frequency distributions. The mean coordinates were then calculated from the frequency distributions by 1
J-xD(x)dx
Jk) dx
(7)
where D(x) is the unnormalized frequency distribution and x is a generic coordinate. The computer program was also coded to obtain the temperature factors of the atoms as follows: temperature, or DebyeWaller, factors were calculated according to ref 27, assuming symmetrical motion of the atoms about the equilibrium positions
T(s) = exp{-'/zsT(u.uT)s)
(8)
where s is a vector of the reciprocal lattice of the crystal (in Cartesian coordinates) and u = r - (r) is the displacement of an atom from the equilibrium position, denoting with r and (r) the current and mean position, respectively. The elements of the symmetric 3 X 3 matrix B = ( u,uT),commonly called anisotropic temperature factors matrix, can be computed from MD trajectories through the following formula Bi, = ( U P , ) = (rirj) -
TO
1.605 1.760 2.6 1186 2.87070 2.45
A.
Finally, in order to get correct repulsive behavior for r the following inequality must be verified:
(x) =
520.99
C
(ri)(rj)
Figure 1. Plot of atomic positions in the large cage of dehydrated Zeolite 4A.
where the brackets denote time averages. If the s vector is expressed in reciprocal crystallographic coordinates, the so-called /3 matrix is obtained. I R spectra were calculated on the basis of linear response from the Fourier transform of the total dipole moment autocorrelation function. As no charge is explicitly contained in the potentials, in order to compute the dipole moment of the system, charges have to be given to the atoms. Standard charges, i.e., +4e for Si, +3e for Al, -2e for 0, and +e for Na, were assigned. This choice was made in several papers related to electrostatic computations; it may affect intensities, but not frequencies, of calculated spectra. An investigation on the effects of charges on spectral properties is in progress. In order to compare theoretical and experimental spectra, the absolute calculated intensities were transformed into relative transmittances following the Lambert-Beer equation. To estimate the contribution of the motion of each atomic species to the spectrum, Fourier transform of the total velocity autocorrelation function (power spectrum) and of the autocorrelation functions of the velocities of atoms of a particular species were evaluated, as usual. By comparing total and particular power spectra, it is possible to estimate the frequency bands corresponding to the motion of single atomic species.
MD Simulations The structure of Zeolite 4A (see Figure 1) was thought to be the most reliably established of all synthetic zeolites characterized to date, but there has been some discussion on the space group to be attributed to its dehydrated form.25q33-35-40The results derived from the structural refinement of the dehydrated form of Zeolite 4A in the Fm3c space group ( a = 24.555 A) by Pluth (28) Gordon, R. G. Adu. Magn. Reson. 1968, 3, 1 . (29) Berne, B. J. In Physical Chemistry-An Advanced Treatise; Eyring, H . , Henderson, D., Jost, W., Eds.; Academic: New York, 1971; Vol. VI11
(9) 1977, 67, 404.
(26) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76,631. (27) Willis, B. T. M.; Pryor, A. W. Thermal Vibrations in Crystallography; Cambridge University Press: Cambridge, 1975.
(32) Berens, P. H.; Wilson, K. R. J . Chem. Phys. 1981, 74, 4872. (33) Reed, T. B.; Breck, D. W. J . Am. Chem. SOC.1956, 78, 5912. (34) Subramanian, B.; Seff, K. J . Phys. Chem. 1977.81, 2249. (35) Bursill, L. A.; Lodge, E. A,: Thomas, J. M.; Cheetham,A. K. J. Phys. Chem. 1981, 85, 2409.
870 The Journal of Physical Chemistry, Vol. 92, No. 4, 1988
Demontis et al.
TABLE III: Experimentpl and Calculated Structural Parameters of Zeolite A (SpaceGroup Fm3c )" 822
833
9.7 3.14 2.49
7.9 6.16 5.36
5.9 8.31 7.99
0 1.57 1.52
0 3.04 3.03
1.2 5.09 4.85
0.090 42 0.090 91 0.093 49
9.4 3.17 2.54
6.5 8.32 8.04
8.7 6.10 5.32
0 3.03 3.06
0 1.53 1.48
2.0 5.08 4.84
0.1 13 67 0.1 17 96 0.12687
0.246 63 0.246 3 1 0.246 29
15.5 23.8 7.83
18.4 14.3 7.50
5.5 9.84 9.75
0 0.48 2.05
0 4.02 3.97
0 5.91 6.0
o.o0ooo 0.000 06 000 12
0.14459 0.13899 0.132 13
0.145 91 0.141 92 0.135 13
23.9 4.04 2.02
11.1
7.68 6.57
10.8 7.70 6.64
2.23 2.13
2.28 2.16
5.79 5.14
0.053 79 0.052 05 0.053 05
0.058 65 0.057 43 0.058 40
0.171 52 0.179 22 0.18532
11.6 19.7 9.72
11.9 19.8 9.98
13.6 21.3 10.7
3.6 17.5 8.02
1.2 16.2 7.48
0.1 16.3 7.59
0.099 60 0.094 38 0.092 13
0.099 60 0.094 43 0.092 14
0.099 60 0.094 61 0.092 16
15.7 11.9 7.07
15.7 12.0 7.1 1
15.7 11.9 7.04
5.4 8.07 4.96
5.4 8.03 4.92
5.4 8.12 4.97
0.00000 0.002 77 01424
0.21650 0.203 33 0.19896
0.211 10 0.20056 0.19551
32 27.3 15.3
19 27.9 16.4
0.25000 0.250 19 0.250 17
0.10600 0.10791 0.10668
0.10600 0.108 16 0.10695
15 12.4 7.51
15 11.6 7.71
xla
Ylb
ZIC
anhar
0.000 00 000 53 00021
0.093 16 0.092 90 0.095 53
0.184 99 0.18762 0.18928
exptl har anhar
0.000 00 000 55 000 21
0.187 15 0.18871 0.19045
0.Ooo 00 0.002 86 0.000 46
atom
811
812
813
823
Si exptl har AI
O(1)
exptl har anhar
00)
exptl har anhar
o(3)
exptl har anhar Na( 1) exptl har anhar Na(2) exptl har anhar Na(3) exptl har anhar
28 372. 534. 24 23.7 17.1
0 81.3 -146.
-2 5.50 5.31
0 79.5 -143. 4 6.44 6.30
-I 1 16.6 11.4 -2 1.67 3.01
"Anisotropic displacement factors multiplied by lo4 and given by 3
T = exp(-&&jhihj) I
and Smith25were chosen as reference data for the M D studies, but a M D simulation of Zeolite 4A in the Pmjm space group ( a = 12.2775 A)34was performed, too. The AI, Si ordering of the framework tetrahedral sites is respected in the former space group and lost in the latter. Both structures are characterized by fractional occupancies in some extraframework sodium sites. The M D system used for the simulations in the Fm% space group was a cubic box of 24.555-A side, corresponding to one crystallographic cell, containing 662 atoms (96 Si, 96 AI, 96 Na, 384 0). The M D system used for the simulation in the Pm3m space group was a cubic box of 24.584-A side, corresponding to eight primitive crystallographic cells, containing the same number of atoms (192 Si, Al;96 Na; 384 0)as for the FmJc space group. The crystal chemistry of the tetrahedral site was simulated alternating Si and AI atoms in the framework tetrahedra, in line with the Si, A1 ordering found in the Fmfc structure (see Figure 1). In both structures the cation sites with partial occupancies were located a t the maximum possible distance among them, in order to simulate a minimum repulsive energy distribution. In this work, the rhombohedral phase recently described by some author^^^*^ for Zeolite 4A was not considered, due to the present lack of an accurate structure refinement in this space group. On the other hand, this phase differs from the cubic one only as regards the ordering in the Na(2) site distribution, which induces a very slight distortion in the aluminosilicate framework. (36) Filippini, G.; Gramaccioli, C. M.; Simonetta, M.; Suffritti, G. B. Acta Crystallogr., Sect. A 1976, A32, 259. (37) (a) Della Valle, R. G.; Pawley, G. S. Acta CrystaUogr.,Sect. A 1984, A40, 297. (b) Heiser, G. A.; Shukla, R. C.; Cowley, E. R. Phys. Rev. B Condens. Matter 1986, 33, 2158. (38) Demontis, P.; Suffritti, G. B.; Fois, E. S.; Gamba, A. Chem. Phys. Lett. 1986, 127, 456. (39) Miecznikowski, A.; Hanuza, J. Zeolites 1985, 5, 188. (40) Takaishi, T.; Ohgushi, T.; Nonaka, K. In Zeolites: Synthesis, structure, technology and Application; DrTaj, B., Hocevar, S., Pejornik, S., Eds.; Elsevier: Amsterdam, 1985; p 467.
Moreover, the rhombohedral-cubic transition takes place at 335 K, a temperature close to that of the simulated system. ps, and The time step used in M D simulations was 2 X the fluctuations of total energy were less than 0.1%. After equilibration, for the Fmjc space group a run 21.6 ps long was performed at room temperature for both harmonic and anharmonic models. Each trajectory required about 91 and 103 min, respectively, on CRAY X-MP12 vector computer. For the Pm3m space group a run 11.3 ps long was performed at room temperature in the anharmonic model.
Results and Discussion The results of the simulations of Zeolite 4A structure in the Fm3c space group are collected in Table I11 where the mean values of the calculated crystallographic coordinates and anisotropic temperature factors are reported and compared with experimental data. The agreement between calculated and experimental coordinates is generally quite good, the standard error being 0.1 1 and 0.20 A for harmonic and anharmonic models, respectively. The distribution functions for both models are all Gaussian-like except for the x coordinate of Na(2). This finding can be explained by the peculiar situation of this atom, lying in the plane of the 8-ring and linked to the ring oxygen atoms by three coplanar bonds, so that in our model Na(2) moves almost freely up and down the 8R plane. In this respect, the dynamics of this atom (and of the other corresponding to it by crystal symmetry) is not well reproduced by the present simulation and probably long-range potential functions would be necessary in order to correct this shortcoming. Anyway, the mean x coordinate of Na(2) is reasonably good. As the simulations were made without any symmetry constraint, the unimodal and regular shape of the distribution functions auggests that the calculated structures keep the experimental symmetry of the crystal. Moreover, the coordinates of some special position atoms like Na(1) and Na(3) deserve a comment. The
The Journal of Physical Chemistry, Vol. 92, No. 4, 1988 871
Molecular Dynamics Studies on Zeolites
‘--r
I
I
cm-‘
Figure 2. IR spectrum of Zeolite 4A calculated with the harmonic model (continuous line). The experimental IR spectrum (thin) taken from ref 23 is reported for comparison, where dashed bands are attributed to linkages external to tetrahedra and the others to internal ones.
three N a ( 1) coordinates are equal by crystal symmetry; the calculated ones are close within 0.008 and 0.0007 8, for the harmonic and anharmonic model, respectively. Likewise, the y and z coordinates of Na(3) are close within 0.007 8, for both models. These results confirm that the crystal symmetry is well reproduced by the .MD simulations. From the comparison of experimental and calculated anisotropic temperature factors, reported in Table 111, it appears that the calculated ones are of the correct order of magnitude, but in general they are smaller than the experimental values. This underestimation can be explained by recalling that, in the present M D simulations, long-wave (acoustic) phonons cannot be reproduced because of the relatively small dimensions of the system. Indeed, these phonons give a large contribution to the temperature factors. This is the M D analogue of the Brillouin zone sampling problem in lattice dynamics36 and could be solved by performing M D simulations with systems containing a very large number of crystallographic cells (see ref 37a,b). The lack of acoustic phonon contribution may explain the relatively large values of off-diagonal temperature factors, because in cubic crystals acoustic vibrations are likely to be mainly polarized along crystallographic axes. In Table I11 the calculated PI’, bl2,and P I 3 of Na(2) are very large in correspondence with the exaggerated amplitude of motion of this atom along x coordinate (see above). The symmetry constraints on the elements of the @ matrix are generally well reproduced except for some off-diagonal elements, but this effect may be due to the lack of acoustic phonon contribution. In order to obtain a more direct description of the atomic thermal motions, an analysis of thermal vibration ellipsoids was performed by ORTEP computer program.41 It uses as input data the anisotropic temperature factor coefficients derived from experiment and M D simulations and gives as results the lengths and orientations, in the Cartesian space, of the principal thermal ellipsoid axes. These data, available as supplementary material, confirm the above-mentioned underestimation of the calculated thermal parameters and emphasize the direction dependency of the atomic thermal motions. In particular, we stress the excellent agreement between the orientation of calculated and experimental thermal ellipsoids of N a ( l ) , Na(3), and O(3) atoms. The MD simulation in the PmSm space group started from the structural data refined by Subramanian and Seff34and, even in the equilibration phase, showed that the structure described in (41) Johnson, C. K. “ORTEP A FORTRAN
Thermal Ellipsoid Plot
Program”, ORNL-3794, Oak Ridge National Laboratory, Oak Ridge, TN,
1965.
this model by a primitive cell was not stable. Indeed, assuming this crystal symmetry, the distribution functions of the atom coordinates were not unimodal; moreover, the single atomic coordinates were shifted toward the corresponding ones of the F m k cell. The trajectory performed at room temperature confirmed this trend. This result strengthens the opinion that M D is a helpful technique for space group assignment in experimentally difficult and/or doubtful cases. The IR spectrum simulated in the F m k cell in the framework of the harmonic model is reported in Figure 2, together with the experimental The agreement of the frequencies and of the overall trend of intensities is satisfactory. The results for the anharmonic model show correct bandwidths and frequencies, unlike the intensities which are not well reproduced. However, it is known that the calculated intensities may depend on the length of the M D run38 and on the charges assigned to the atoms. From the power spectra of Na+ (not shown) it emerges that they are involved in many vibrational modes of the crystal, in a wide range of frequencies, which is larger for anharmonic model, but in any case lower than 700 cm-I. This finding agrees with the assignment of bands, with frequency higher than 750 cm-’ to internal tetrahedra s t r e t c h i n g ~ . ~ ~ Miecznikowski and H a n ~ z maintain a ~ ~ that the frequency of the vibrations of Na+ cations should fall in the range 160-225 cm-’. Indeed, the simulated spectra show in this range a few high-intensity components, but they are blurred by the low-frequency motion of Na(2) atoms (see above). We are aware that a classical method is not the most suitable for describing accurately properties having a quantum origin. Moreover, in the case of a crystal as complex as Zeolite 4A, it is not easy to derive a detailed analysis of normal modes, like in the “simple” case described in ref 11. What we intend to say is that molecular dynamics not only is founded on sound physical principles but also is able to “take a look” at the overall I R spectrum, giving both a reliable representation of the frequency pattern and some characterization of the single atomic species motion. Thus, in this respect, lattice (or pseudolattice) dynamics and M D appear to be complementary techniques.
Concluding Remarks The results obtained in this work confirm that M D can satisfactorily reproduce the structural and dynamic features of zeolites with as different topoiogies as Natrolite (ref 18 and 19) and Zeolite 4A (this work). However, some work must be done in order to refine the potential models and to extend the calculations to exchanged and/or hydrated zeolites. Moreover, some theoretical problems still need to be solved as regards the simulation of IR spectra, a field where some recent developments were a c h i e ~ e d . ~ ~ , ~ ~ Anyway, it is doubtless that MD could become an useful predictive tool in the study of crystallographic and spectroscopic problems for which the experimental data show some deficiencies or are difficult to interpret. Acknowledgment. We are grateful to Prof. A. Alberti for useful discussion and critical reading. The Minister0 della Pubblica Istruzione Italian0 (M.P.I. 40% and 60%) is acknowledged for financial support. We thank the “Centro Interdipartimentale di Calcolo ed Informatica Applicata” (CICAIA) of the University of Modena for the computer facilities. Supplementary Material Available: A listing (Table IV) of the observed and calculated lengths and orientations of the principal thermal ellipsoid axes (2 pages). Ordering information is given on any current masthead page. (42) Koszykowsky, M. L.; Noid, D. W.; Marcus, R. A. J . Phys. Chem. 1982, 86, 2113.
(43) Noid, D. W.; Koszykowsky, M. L.; Marcus, R. A. Annu. Reu. Phys. Chem. 1981, 32, 261.
