Faraday Discuss. 1997 106 119»133 General Discussion Dr Schoé n opened the discussion of the Introductory Lecture We have been using global optimisation methods to –nd the local minima of the potential-energy hypersurface in order to predict structure candidates for not-yet-synthesized compounds.1 Information about neither the exact composition nor the parameters or symmetry of the unit cell of the hypothetical compound is available and thus these quantities need to be determined during the global optimization»in contrast to the case where the compound has already been synthesized and e.g. powder diÜraction data or information about local coordination of the atoms/ions are available. Therefore it is crucial to employ potentials that can describe reasonably well as large a region of the energy landscape of the chemical system as possible in order to encompass a large range of potential structure candidates.However we often –nd that e.g. the transferability of potentials derived for known binary compounds to the unknown ternary compound under investigation is beset with problems. Thus I wonder (i) whether you have encountered similar problems and (ii) what route you would suggest for dealing with this problem when deriving the parameters for empirical potentials ? 1 J. C. Schoé n and M. Jansen Angew. Chem. Int. Ed. Engl. 1996 35 1286. Prof. Catlow replied You are right that transferability of potentials is a key issue. It is well known that in inorganic materials there are systematic changes in parameters on varying coordination numbers.However there are a number of cases where it has proved possible to derive reasonably transferable sets of parameters for oxide materials as in the work of Bush et al.1 I think that the best procedure to try to ensure transferability when using empirical potentials is to build into the potential functions systematic variations with key structural features e.g. coordination number as was done in the four-fold coordination. recent study of Takada et al.2 on B2O3 where boron can show both three-fold and 1 T. S. Bush J. D. Gale C. R. A. Catlow and P. D. Battle J. Mater. Chem. 1994 4 831. 2 A. Takada C. R. A. Catlow and G. D. Price J. Phys. Condens. Matter 1995 7 8659. Prof. Gillan commented I am a strong believer in the use of empirical models but it is worth commenting on the way that the role of these models is changing.With the rapid increase in the power of ab initio methods [both Hartree»Fock (HF) and density functional theory (DFT)] it is becoming possible to do with ab initio what would have been done before with empirical models. It is now possible to do accurate ab initio calculations (both static and dynamic) on systems of over 100 atoms. Good examples of this are recent simulations by Schwarz and co-workers,1 and Payne and co-workers2,3 on methanol in zeolites. 1 E. Nusterer P. E. Bloé che and K. Schwarz Chem. Phys. L ett. 1996 253 448. 2 J. D. Gale R. Shah M. C. Payne and I. Stich Abst. Pap. Am. Chem. Soc. 1997 213 20. 3 R. Shah M. C. Payne M-H. Lee and J. Gale Science 1996 271 1395.Prof. Catlow added I agree. Many problems that were only tractable a few years ago by using interatomic potential methods can now be investigated by electronic structure techniques and the range and scope of the latter methods will continue to expand. Nevertheless I am convinced for three reasons that there is a continuing role for the former techniques. First as we push simulation techniques to increasingly large and 119 120 General Discussion complex systems the interatomic-potential based approaches will continue to be needed as there will remain a domain of complexity which is beyond the reach of electronic structure techniques. Secondly there will be a need for such methods to provide approximate structures and models which can then be re–ned by electronic structure techniques.Thirdly there are a number of problems that are handled more appropriately by interatomic-potential based methods including for example the modelling of microporous crystal structures and the diÜusion of hydrocarbons within these structures. In summary there will be an on-going need to use methods employing interatomic potentials which are increasingly complementary to those based on ab initio techniques. Prof. Klein commented Both speakers (Prof. Gillan and Prof. Catlow) are correct. There is going to be a continuing need for empirical potentials to deal with large aperiodic systems. The current limitation of a few hundred atoms is inadequate to handle many interesting systems. Hence the continuing need for embedding methods which in turn will require improved empirical models.Sir John Meurig Thomas commented The point raised by Prof. Klein concerning the need for massive computational eÜort if one wishes to address the question of the electronic and structural properties of complicated zeolitic structures merits further comment. When one is interested (as I am1,2) in the role of transition metals (in the framework of open-structure molecular sieves) as selective oxidation catalysts it is quite realistic to take metal silsesquioxane (molecular) analogues3 because these have been shown4 to perform (as epoxidation catalysts) in a manner closely akin to the performance of titanosilicate catalysts. I know from Fig. 2 of his paper that Prof. Sauer has looked at polyhedral oligomeric silsesquioxanes [molecular structure (XSiO1.5)n with n\6 or 8] but more needs to be done computationally on these molecular entities with ions such as TiIV located at one of their vertices to which is attached a hydroxy group (just as exists in single-site titanosilicate solid oxidation catalysts).1 J. M. Thomas Philos. T rans. R. Soc. L ondon A 1990 333 17. 2 T. Maschmeyer F. Rey G. Sankar and J. M. Thomas Nature (L ondon) 1995 378 159. 3 H. W. Roesky R. Murugavel and V. Ahandrasekhar Accs. Chem. Res. 1996 29 183. 4 T. Maschmeyer M. Klunduk C. M. Martin J. M. Thomas B. F. G. Johnson and D. S. Shephard Chem. Commun. 1997 1847. Prof. Martin addressed Prof. Catlow In your paper you have discussed mainly isolated defects and their formation energies.However we also need calculations of the concentrations of defects. This means that we also need entropies or Gibbs energies of formation of defects to be able to use the law of mass action. In addition the interaction of defects (resulting possibly in defect clusters) and its in—uence on the energies and entropies of formation is of great interest. Could you comment on the actual possibilities in calculating such quantities and the agreement with experiment e.g. in systems where defect clusters are formed. Prof. Catlow responded It is certainly possible to calculate defect entropies using methods developed by Jacobs Gillan and Harding. Such methods were used by e.g. Tomlinson et al.1 to study defect clustering in transition metal oxides using a mass action formalism.Calculations on defect supercells which are increasingly feasible may also be used to investigate defect»defect interactions. 1 S. M. Tomlinson C. R. A. Catlow and J. H. Harding J. Phys. Chem. Solids 1990 51 477. Dr Allan communicated I have a brief comment on the calculated enthalpy of mixing for NiO»MnO presented by Prof. Catlow. Calculations on such mixtures have 121 General Discussion generally involved determining the energy of the smallest supercell consistent with the composition of interest (e.g. Ni MnO for an overall Mn content of 25%). If so care 3 may be required when comparing calculation and experiment owing to the arti–cial 4 periodicity and the possible neglect of clustering eÜects. It is worth referring to a few of our results.1 We have recently been studying the enthalpy of mixing of MnO and MgO.Using a combination of lattice statics and quasiharmonic lattice dynamics with two-body potentials,2 then the lowest energy supercell * of formula Mg6Mn2O8 indicates an enthalpy of mixing mixH at 1300 K of 4.8 eV. The larger supercell Mg12Mn4O16 [with the same Mn content (25%)] with more internal * degrees of freedom suggests a value of mixH smaller by ca. 25% i.e. 3.6 eV. A value * for mixH of 3.6 eV is also obtained in preliminary calculations using the same potentials and a new hybrid Monte Carlo»molecular dynamics method currently under development at Bristol for this type of problem. 1 J. A. Purton G. D. Barrera and N. L. Allan in preparation. 2 G. V.Lewis and C. R. A. Catlow J. Phys. C. Solid State Phys. 1985 18 1149. Prof. Cheetham opened the discussion of Prof. Sauerœs paper You describe embedded cluster calculations on the pure silica polymorph of faujasite. Have you been able to compare these calculations with periodic HF or DFT calculations on this system? The symmetry is very high and there are only –ve atoms in the asymmetric unit. Such a calculation would give us a higher level of con–dence in the embedding procedure. Prof. Sauer replied When we completed our study periodic HF and DFT results on silica faujasite were not available. At this meeting Prof. Dovesi presents a poster which shows that periodic HF calculations on silica faujasite are feasible using an improved CRYSTAL code. The use of symmetry-adapted crystal orbitals results in a blockdiagonal Fock matrix at each k point in reciprocal space.Table 1 shows the results obtained for faujasite.1 Our paper describes DFT not HF calculations. However we have published embedded cluster calculations using the HF method previously ;2 the results are also shown in Table 1. Although there are some diÜerences in the basis sets used the two sets of bond distances agree within 0.01 ”. The largest deviation of the bond angle is 5° but we should take into account that the bending potential of the SiwOwSi bond angle in silica is very —at. The sequence of observed bond angles3 is reproduced by both sets of calculations. Moreover the deviation between the two sets of calculated data is of the same magnitude as the remaining deviation from experiment owing to systematic errors connected with the HF approximation and the basis set used.In ref. 2 comparison has been made with HF calculations on the ìexpandedœ structure of sodalite using exactly the same basis set. For the SiwO bond distance and the Table 1 Observeda and calculatedb,c bond lengths and bond angles of silica faujasite SiwOwSi/degrees r(SiwO)/” QM-Potc obs.a QM-Potc obs.a QMperiodicb QMperiodicb 137 155 149 142 154 149 138 149 146 143 145 141 1.622 1.608 1.622 1.612 24.63 1.61 1.60 1.62 1.61 24.53 1.607 1.597 1.604 1.614 24.26 O(1) O(2) O(3) O(4) cell a Ref. 3. b Ref. 1. c Ref. 2. 122 General Discussion Si ” wOwSi bond angle diÜerent embedded cluster calculations yield 1.607 to 1.608 and 161.5 to 160.7° while the periodic calculation (CRYSTAL code) yields 1.612 ” and 160.6°.1 C. M. Zicovich-Wilson and R. Dovesi Chem. Phys. L ett. submitted. 2 U. Eichler C. M. Koé lmel and J. Sauer J. Comput. Chem. 1996 18 463. 3 J. A. Hriljac M. M. Eddy A. K. Cheetham J. A. Donohue and G. J. Ray J. Solid State Chem. 1993 106 66. Prof. Dovesi added Recent improvement of the CRYSTAL code has permitted a study of the silicon-only faujasite where the full symmetry of the system is maintained by using an all-electron basis set. These data together with a full optimization of the structure have recently been reported.1 1 C. M. Zicovich-Wilson and R. Dovesi Chem. Phys.L ett. submitted. Dr Hill addressed Prof. Sauer In your paper you mention the calculation of NMR chemical shifts for embedded systems. It is clear that the embedding method you use does only a mechanical and not an electronic embedding. NMR chemical shifts are clearly an electronic eÜect. How important is the inclusion of electronic eÜects in the embedding to calculate NMR chemical shifts and how large is the error you estimate for your mechanical embedding? Prof. Sauer responded Indeed in our embedding scheme the wavefunction of the cluster does not ìseeœ the electron distribution of the outer part. However the electrostatic interaction and mutual polarization between the inner and outer parts are described at the level of the interatomic potential function (shell model potential).This changes the structure of the embedded cluster and indirectly aÜects its wavefunction since a diÜerent structure will lead to a diÜerent wavefunction. As far as the calculation of NMR shifts is concerned we know that NMR shielding constants are structure-sensitive parameters and we can get them right by calculation only if the structure is right. Our calculation proceeds as follows we use the structure from the embedded cluster calculation de–ne a (generally diÜerent) cluster around the nucleus the shielding of which we are interested in and add an increasing number of coordination shells to this cluster until the calculated value of the shielding constant is converged. Speci–c information is given in ref. 1 and 2. The remaining error due to the –nite-size model is smaller than the errors connected with the CPHF method used and the basis sets employed.The latter errors are minimized in our approach by using secondary internal standards. 1 U. Eichler J. Sauer and M. Braé ndle J. Phys. Chem. B in press. 2 B. Bussemer and J. Sauer Solid State NMR in press. Prof. Sauer continued I should like to add some general comments on diÜerent embedding techniques. Embedding is an obvious idea but the diÜerent approaches are very diÜerent in their aims methods and implementations. There are theoretically very appealing ì electronic œ embedding techniques which start from the known periodic solution of the ideal crystal and look for the solution of a defect crystal. One example of a working implementation is the EMBED code of Pisani et al.1 It has been successfully used in several applications but a severe limitation of the approach is the assumption that the density matrix for the defect crystal outside the internal part is identical with the density matrix of the ideal crystal.This means that charge transfer between the internal part and the environment is forbidden. Although Pisani et al. found a method to make an a posteriori correction this is not satisfying for cases where the correction is as large as the eÜect itself. An example is the application of the EMBED code to the proton transfer from the acidic zeolite to an adsorbed ammonia molecule yielding 123 General Discussion ammonium cations.2 Comparison is made between a supercell calculation using periodic boundary conditions and the embedded cluster calculation.After applying the corrections the adsorption energies still have an error of ca. 20 kJ mol~1. Moreover the computational demands of the periodic calculation neither allow the use of an appropriate basis set (STO-3G was employed) nor can all degrees of freedom be optimized. The result is that contrary to experiments the neutral adsorption complex of ammonia is more stable than the ion-pair complex of the ammonium cation with the negatively charged framework. The error is ca. 200 kJ mol~1 but most of it is due to the poor description of the anionic framework with the STO-3G basis set (see the discussion in Section C.1 in ref. 3). Our embedding scheme is much more pragmatic.We would like to avoid periodic ab initio calculations since we are interested in very large supercells and we consider it vital to relax the structures fully. It is a lesson from gas-phase calculations that reliable reaction energies will not be obtained unless all degrees of freedom are relaxed. We only describe the cluster quantum mechanically and the outer part and the interaction between the internal part and the outer part are described by an ab initio parametrized interatomic potential function. We include mutual polarization of the inner and outer parts by use of the shell model potential. We have applied our approach to the same problem adsorption of and proton transfer to ammonia for the same zeolite.4 We used a double-zeta polarization basis set and fully relaxed all structure parameters.We –nd that the ammonium form is more stable than the neutral adsorption complex and we calculate meaningful values for the heat of adsorption. In addition we have made explicit comparison with a periodic calculation using the same basis set. We employed the CRYSTAL code in single-point calculations on the structure obtained by our combined QM-Pot approach. The error is 6»9 kJ mol~1 less than half the error of the theoretically more appealing ì electronic œ embedding scheme mentioned above. 1 C. Pisani F. Cora` R. Nada and R. Orlando EMBED93 user documentation University of Torino Italy 1993. 2 C. Pisani and U. Birkenheuer Int. J. Quant. Chem. 1995 29 221. 3 J. Sauer P. Ugliengo E. Garrone and V.R. Saunders Chem. Rev. 1994 94 2095. 4 M. Braé ndle J. Sauer R. Dovesi and N. Harrison in preparation. Prof. Hillier commented We need benchmark calculations to assess the validity of various embedding techniques. Perhaps a large plane wave calculation should be carried out to provide this information. Prof. Sauer commented There is the question of what criteria we have to assess the quality of an embedded cluster calculation. The best criteria are internal ones since they do not depend on the systematic deviation from observed data which every quantum mechanical method in combination with any choice of basis set shows. One such criterion is that the result should be stable with respect to the choice of diÜerent clusters for the same problem in particular with respect to a series of clusters of increasing size.Another criterion is agreement with fully periodic calculations using the same quantum method and basis set as used for the embedded cluster calculations. Prof. Klein said Linear scaling will enable the study of several thousand atoms. The prospects are excellent that this will be possible in the near future. However larger system often have important dynamical processes that occur in the nanosecond time domain. What are the prospects for achieving this length of simulation with ab initio methods? I suggest that at present it doesnœt look too hopeful unless some form of embedding technique is employed. Prof. Gillan responded The problems of going to larger systems and longer times are completely diÜerent.For large systems the problem is that for traditional ab initio 124 General Discussion methods the computer time is proportional to at least N2 (where N is the number of atoms) and the challenge is to develop methods in which the eÜort is proportional to N. It is now virtually certain that this will be achieved (as it has been already within tightbinding methods). However for long times the eÜort is already proportional to time and there doesnœt seem to be anything straightforward that can be done to speed up calculations. Dr Line commented Although we obviously cannot calculate the dynamics over such long timescales for large systems we surely can calculate quantities such as the activation energy for diÜusion which is an experimentally measurable quantity in quasielastic neutron scattering.Dr Schoé n asked Prof. Sauer Have you tried to go beyond the harmonic approximation in your calculation of free energies of a- and b-quartz ? Prof. Sauer replied No. This is feasible but would be computationally very demanding. For example Tsuneyuki et al.1 simulated the a»b phase transition of quartz by molecular dynamics. A 4]4]3 supercell (432 atoms) was used and 6000 time steps each for 18 diÜerent temperatures were found to be necessary. 1 S. Tsuneyuki H. Aoki M. Tsukada and Y. Matsui Phys. Rev. L ett. 1990 64 776. Prof. Hillier opened the discussion of Prof. Gavezzottiœs paper The stabilization of the 2-pyridone tautomer compared to 2-hydroxypyridine in a polar environment arises to a signi–cant extent from the preferential electronic polarization of the 2-pyridone.Would the inclusion of such polarization aÜect the conclusions of your simulation studies ? Prof. Gavezzotti responded The pyridone tautomer is stable in condensed phases and this justi–es our assumption. On the other hand polarization certainly aÜects all force –elds but there are no simple methods to account for it. Dr Line asked Have you carried out calculations of the polarization of the pyridine complexes and of the solvent eÜects of surrounding water molecules? These eÜects surely must be signi–cant for pyridine since it is polar. Prof. Gavezzotti responded The answer is unfortunately no for reasons alluded to in my reply to Prof. Hillier. Dr Price addressed Prof.Gavezzotti One very signi–cant result in your paper is that it clearly demonstrates the limitations of the assumption that organic crystal structures can be predicted by searching for the global minimum in the lattice energy. This assumption is the basis of several methods aiming at crystal structure prediction,1 including a commercial package. I have done a brief study of the possible crystal structures of 2-pyridone using a diÜerent model intermolecular potential and a diÜerent search method and obtain complementary results. I used an accurate model for the electrostatic component of the intermolecular forces using atomic multipoles up to hexadecapole to describe the atomic charge distribution derived by a distributed multipole analysis (DMA)2 of an MP2 6-31G** wavefunction for the molecule.The other intermolecular contributions were described by a diÜerent set of 6-exp atom»atom potential parameters which have been found to reproduce successfully the crystal structures of related molecules3 when used in conjunction with a DMA-based electrostatic model. Thus the two model potentials one with an accurate electrostatic model the other without any explicit representation of the electro-125 General Discussion static forces but with empirically –tted hydrogen-bonding potentials are very diÜerent in their philosophy. As shown in Table 2 the DMA-based potential also gives a minimum in the lattice energy reasonably close to the room-temperature structure with a reasonable lattice energy. The two potentials diÜer in which cell lengths are most signi–cantly in error con–rming that they are far from equivalent but yet reasonable approximations to the intermolecular potential for 2-pyridone.Using the DMA-based model potential I minimized the lattice energy starting from 100 hypothetical crystal structures for 2-pyridone. These were generated by a MOLPAK4 search for close-packed crystal structures in the 10 commonly observed packing coordinations in the –ve space groups in Table 1 of your paper. For each coordination type the 10 lowest-energy of the 25 most-dense structures were minimized using the DMA-based potential giving a mere 100 minimizations. The results in Fig. 1 (cf. Fig. 6 of your paper) con–rm that there are a large number of crystal structures within the energy range of possible polymorphs.The experimental structure was found easily with an exact match of the in–nite hydrogen-bonded chain structure as found in the structure obtained in minimizing from the experimental structure. However a P21/c structure involving hydrogen-bonded dimers was [1.3 kJ mol~1 lower in energy. Hence there are de–nitely two qualitatively very diÜerent types of crystal structure for 2-pyridone based on either dimers or catemers whose static energy diÜerence is even smaller for the more elaborate theoretically based potential than for the empirical crystal potentials. This is not an unusual case we have observed similar small energy diÜerences between minimum-energy structures with diÜerent hydrogen-bonding motifs for uracil and 6-azauracil.5 Hence the need to go to more elaborate dynamic simulation methods even to predict qualitatively which type of crystal structure might be found.My question based on the molecular dynamics simulations is what sort of simulation results might have led to a fairly con–dent prediction that a catemer would be found rather than a dimer in the experimental structure ? For example if there had been no ì—ying catemersœ in the carbon tetrachloride solution simulation would that have Fig. 1 Predicted minima in the lattice energy for 2-pyridone using a DMA-based potential showing the exact coincidence (superimposed square and diamond in the top left-hand corner) of the lowest-energy P212121 minimum with the minimum found from the P212121 experimental structure.All minima (except min from experiment) obtained using MOLPAK/DMAREL from 100 close-packed starting structures. Table 2 Results of static minimization starting from the experimental b/” a/” *a/a (%)a 5.900 5.795 5.41 5.64 » » ]1.3 ]5.3 13.645 13.564 13.82 expt 295 K expt 120 K empirical 14.37 DMA based a The percentage errors are relative to the room-temperature structure. P212121 crystal structure of 2-pyridone *c/c (%)a *b/b (%)a c/” 5.692 5.604 5.75 5.80 » » ]1.0 ]1.8 » » [8.3 [4.4 lattice energy /kJ mol~1 *V /V (%)a [86.6 (*Hsub) » » [6.1 ]2.5 » [81.4 [83.3 126 General Discussion 127 General Discussion implied that a dimer-based crystal structure would crystallize from carbon tetrachloride ? 1 H.R. Karfunkel and R. J. Gdanitz J. Comput. Chem. 1992 13 1171; A. M. Chaka R. Zaniewski W. Youngs C. Tessier and G. Klopman Acta Crystallogr. Sect. B 1996 52 165; B. P. van Eijck W. T. M. Mooij and J. Kroon Acta Crystallogr. Sect. B. 1995 51 99. 2 A. J. Stone and M. Alderton Mol. Phys. 1985 56 1047. 3 D. S. Coombes S. L. Price D. J. Willock and M. Leslie J. Phys. Chem. 1996 100 7352. 4 J. R. Holden Z. Du and H. L. Ammon J. Comput. Chem. 1993 14 422. 5 S. L. Price and K. S. Wibley J. Phys. Chem. A 1997 101 2198. Prof. Gavezzotti replied Molecular dynamics simulations for precursors to crystal structures are just tentative so far. I am afraid the question cannot be answered before (much) more experience is gained.Prof. Klein commented Given the large number of possible structures for your system the crystallization pathway may be complex. You seem to be dealing with competition between hydrogen bonding and van der Waals interactions. Perhaps you could look at an analogy to the protein folding problem to understand what is going on. There the idea of a ìrugged energy landscapeœ has evolved and the idea of funnelling down to the correct structure. In proteins the van der Waals (and entropy) terms are extremely important. That is also likely to be the case here. Prof. Cheetham said As a point of clari–cation I would like to point out that the number of hydrogen bonds per molecule is the same for both dimer and chain structures of the pyridone.Dr Schoé n commented Concerning the ruggedness of energy landscapes we have been studying the energy landscape of small polymers on a periodic lattice interacting via simple potentials1 using the lid-algorithm2 to determine the growth in the number of states and local minima as a function of the energy barriers one is able to cross. We –nd (i) many nearly degenerate deep-lying minima (ground states) and (ii) a highly rugged landscape where the number of accessible local minima grows exponentially with the height of the barrier that can be crossed starting from one of the ground states. Thus the large number of local minima found during the work presented here really should not be too surprising. Finally even for very simple ionic systems like Na`/Cl~ there exists a multitude of local minima representing very distinct structures that while not having been found experimentally in NaCl are often observed in other ionic or partly ionic AB compounds.3 1 J. C. Schoé n Habilitationsschrift Universitaé t Bonn 1997. 2 P. Sibani J. C. Schoé n P. Salamon and J-O. Andersson Europhys. L ett. 1993 22 479. 3 J. C. Schoé n and M. Jansen Comput. Mater. Sci. 1995 4 43. Prof. Jansen commented When comparing the molecular packings obtained by computer simulations with those found experimentally e.g. by crystallization from solution it is not sufficient just to discuss the energies of the –nal solids one must also consider kinetic and solvent eÜects governing nucleation and growth of the respective polymorphs. Prof.Cheetham said Dr Priceœs comment»that one typically –nds only one polymorph even though calculation predicts very large numbers»is rather depressing. It is also true for inorganic species (as Kapustinski showed many decades ago). Is it safe to ignore the entropy and base these calculations solely on enthalpy? Prof. Gavezzotti replied Apparently yes,1 but this conclusion is based on lattice vibrational entropy only. The question is one that still awaits a more complete answer. 128 General Discussion 1 A. Gavezzotti and G. Filippini J. Am. Chem. Soc. 1995 117 12 299. Dr Freeman commented Dr Price and Prof. Gavezzotti report that although the experimentally observed polymorphs are among those generated in their studies of the possible packings of molecules using potential-energy functions the experimentally observed structures are not necessarily of lowest energy.This implies of course that when powder diÜraction data can be obtained (as is the case for the majority of organic compounds) polymorph prediction (see for example ref. 1»5) can yield crystal packing information through the match of calculated and observed diÜraction pro–les. Indeed a number of studies have made use of pro–le match as an additional discriminant in simulated annealing and Monte Carlo studies.6h10 It is reasonable to conclude that as a complement to traditional methods and particularly where powder diÜraction is the only possible structural probe polymorph prediction methods are of immediate value. 1 R. J. Gdanitz Chem. Phys.L ett. 1992 190 391. 2 H. R. Karfunkel and R. J. Gdanitz J. Comput. Chem. 1992 13 1171. 3 H. R. Karfunkel and F. J. J. Leusen Speedup 1992 6 43. 4 H. R. Karfunkel B. Rohde F. J. J. Leusen R. J. Gdanitz and G. Rihs J. Comput. Chem. 1993 14 1125. 5 H. R. Karfunkel F. J. J. Leusen and R. J. Gdanitz J. Comput. Aided Mater. Des. 1993 1 177. 6 J. M. Newsam M. W. Deem and C. M. Freeman in Accuracy in Powder DiÜraction II NIST Special Publication No. 846 ed. E. Prince and J. K. Stalick National Institute of Standards and Technology Bethesda MD 1992 pp. 80»91. 7 T. M. NenoÜ W. T. A. Harrison G. D. Stucky J. M. Nicol and J. M. Newsam Zeolites 1993 13 506. 8 M. W. Deem and J. M. Newsam J. Am. Chem. Soc. 1992 114 7189. 9 C. R. A. Catlow J. M. Thomas C. M. Freeman P.A. Wright and R. G. Bell Proc. R. Soc. L ondon A 1993 442 85. 10 K. D. M. Harris M. Tremayne P. Lightfoot and P. G. Bruce J. Am. Chem. Soc. 1994 116 3543. Dr Freeman asked Could you comment on the possible role of crystal structure databases in determining packing rules ? Is it possible that energy functions could provide a convenient route to structural models and their –nal rankings could also draw on empirical relations obtained from known structures ? Prof. Gavezzotti responded I agree with your comment and did so in print.1 To answer your question The use of databases as a help in –nding the correct structure is still more of an art than a systematic exercise because ranking based on geometrical patterns is even more awkward than ranking based on energies.1 A. Gavezzotti and G. Filippini J. Am. Chem. Soc. 1996 118 7153. Dr Schoé n asked Are you familiar with the work by Hofmann and Lengauer1 on the structure prediction of molecular compounds using simple potentials ? 1 D. W. M. Hofmann and T. Lengauer Acta Crystallogr. Sect. A 1997 53 225. Prof. Gavezzotti responded I have not seen Hofmann and Lengauerœs paper. Thank you for pointing it out to me. I will consider it carefully. Dr Schoé n asked How did you calculate the contribution of the potential energy (as a function of temperature) to the total energy when calculating the heat capacity Cp ? Prof. Gavezzotti answered The C was roughly approximated by –nite diÜerences p over total energies at diÜerent temperatures. Clearly a better method is needed here.Dr Line opened the discussion of Dr Sherwoodœs paper In the section on 1H NMR shifts you report the shift of the OH stretch frequencies due to the inclusion of the point charge array relative to that with the Br‘nsted acid cluster alone. We have performed ab initio calculations1 on the water cluster found in the zeolite natrolite. We considered 129 General Discussion only one water molecule and the two neighbouring sodium atoms at 2.4 ” with the zeolite framework completely omitted. We –nd that our OH vibrational frequencies are remarkably similar to yours obtained for the point charge array (3754 and 3842 cm~1 for 6-31G** with MP2). This is therefore from the cation bonding eÜect alone since we did not consider hydrogen bonding at all. Have you tried calculating the OH shift for the point charge array alone without including the Br‘nsted acid site ? 1 C.M. B. Line and G. J. Kearley unpublished work. Dr Sherwood responded No we havenœt tried this and it isnœt clear how such a calculation would be interpreted. In considering the electrostatic environment of the water molecule it must be borne in mind that the strongest eÜects will arise from the closest zeolite framework atoms which are part of the ab initio cluster rather than being included through the point charges. To remove the cluster would leave an arbitrary rather meaningless partial contribution to the actual –eld. There is the additional problem that the structure would no longer be at equilibrium and the computation of IR vibrational frequencies would not be possible.Prof. Sauer commented It is certainly possible to calculate frequency shifts for water molecules adsorbed on metal cations in zeolites. Already very simple calculations (HF Mini-1 basis set) on free Mn`»OH complexes explain some qualitative features. For example for water adsorbed in Y-type zeolites containing diÜerently charged cations 2 cations but to lower values if trivalent cations are present. Accordingly the calculations the bending frequency of H2O is shifted to higher values for monovalent and divalent 2O»Li` and for H 2O»Mg2` yield bending frequencies higher than that of H2O while H for the H2O»Al3` complex a lower frequency than that of free H2O is calculated.1 Much better calculations are possible with todayœs computer codes.1 J. Sauer and R. Zahradnik Int. J. Quant. Chem. 1984 26 793. Dr Islam opened the discussion of Prof. Smitœs paper In your Monte Carlo simulations of hydrocarbon diÜusion in zeolites the host lattice was assumed to be rigid. How important is the incorporation of framework —exibility particularly for diÜusion of large hydrocarbon molecules? How would such a —exible model be included in your calculations ? Prof. Smit responded Since for large hydrocarbons there are no simulations reported on diÜusion in a —exible zeolite I can only speculate on the eÜects of a —exible lattice. Our simulations indicate that the diÜusion of linear alkanes is diÜerent from the diÜusion of branched alkanes. Compared with the branched alkanes the linear alkanes move freely in the channels.The —exibility of the zeolite lattice can have some in—uence on the diÜusion but I do not expect a large eÜect. The diÜusion of the branched alkanes however is an activated process in which the bulky head of the molecule has to squeeze itself through a narrow channel. In a —exible lattice this may be much easier and as a result in a —exible zeolite one would observe a lower free energy barrier. In principle the Rosenbluth free energy calculations can be done with a —exible lattice. In these calculations the hydrocarbon is a ìghost particle œ i.e. it does not in—uence the conformation of the zeolite. In the case of a rigid zeolite this is exactly what one wants to simulate. However in the case of a —exible zeolite this may give inaccurate results.If the hydrocarbon is introduced in the zeolite as a real particle the zeolite will respond to the presence of this particle. If this response is sufficiently small such that linear response theory is valid then one can compute the eÜects from —uctuations of the empty lattice and one can still use the Rosenbluth scheme. However if the perturbations 130 General Discussion of the zeolite caused by the presence of the molecule are large one has to use alternative free energy schemes in which the presence of the particle is explicitly taken into account (for example the blue moon ensemble umbrella sampling or histogram methods).1 1 D. Frenkel and B. Smit in Understanding Molecular Simulations From Algorithms to Applications Academic Press Boston 1996.Dr Yashonath asked (i) You have stated that 2-methylbutane preferentially places the tail in the entrance to the straight channel as compared to the zigzag channel. Why is this ? The dimensions of the straight and zigzag channels are about the same. (ii) You have mentioned the rotation of 2-methylbutane. Did you –nd any preference for rotation around the long molecular axis ? Prof. Smit answered (i) The reason that the tail of 2-methylbutane prefers the straight channel is that because of the bulky head group the zigzag channel is more difficult to reach than the straight channel. (ii) We have not looked at preferences for rotations. Prof. Bliek addressed both Prof. Smit and Dr Yashonath Obviously calculations on the self-diÜusivity of species in zeolitic structures are easier than those on diÜusivities involving counter-diÜusing species.How can such a situation be adequately modelled given the necessity to account for both sorbent»sorbate and sorbate»sorbate interactions. Dr Yashonath replied It is certainly easier to perform calculations on the selfdiÜusivity in zeolites as compared to the situation where a counter-diÜusing species is present. The difficulty is greater if the adsorbed —uid and the counter-diÜusing species are in the liquid phase instead of the vapour phase. The sorbate»zeolite interactions for many species are well known. In particular these are well known for hydrocarbons (see the recent work of Siepmann et al.1). Existing models for the sorbent»sorbate interactions are inadequate even for simple guest species such as rare gas atoms; see for example the recent work by Pellenq and Nicholson2 where they found that the potential of Kiselev and Du,3 used extensively in the literature yielded a signi–cantly higher value for the volume of channels.Signi–cantly more work needs to be done before one is able to obtain fairly accurate potentials for the sorbent»zeolite interactions of other molecules. More important seems to be the fact that the counter-diÜusion normally occurs signi –cantly more slowly than normal diÜusion. Experiments on counter-diÜusion indicate that the diÜusivities in these systems are signi–cantly smaller than those in previously evacuated zeolite samples.4 Molecular dynamics might prove to be inadequate in modelling counter-diÜusion since the diÜusivities are expected to be smaller and one may have to resort to other methods such as smart Monte Carlo or other more efficient sampling techniques.1 J. I. Siepmann M. C. Martin C. J. Mundy and M. L. Klein Mol. Phys. 1997 90 637. 2 J. M. Pellenq and D. Nicholson J. Phys. Chem. 1994 98 13 339. 3 A. V. Kiselev and P. Q. Du J. Chem. Soc. Faraday T rans. 2 1981 77 1. 4 C. N. Satter–eld and J. R. Katzer Adv. Chem. Ser. 1971 102 193; J. Karger and D. M. Ruthven in DiÜusion in Zeolites Wiley New York 1992 p. 433. Prof. Smit responded The simulations we have performed so far describe a single molecule at in–nite dilution. These simulations give a hopping rate of a molecule jumping from one intersection to another.Extending these calculations to higher concentrations is not straightforward. A possible way to study the eÜects of concentration and multi-component systems is to use a lattice model and use hopping rates as com-131 General Discussion puted from our MC scheme. If we assume that these hopping rates are not in—uenced by the presence of other molecules we can use a lattice model in which molecules are allowed to jump from one site to another. With such an approach it may be possible to study concentration eÜects. Prof. Catlow addressed Prof. Smit Your paper comments that the transition state theory approach exaggerates the diÜusion coefficient of the hydrocarbons by a factor of ca. 5»10 compared with corresponding MD simulations. Could you amplify this comment; are there possibilities for reducing this discrepancy ? Prof.Smit replied The hopping rate is a product of a term which includes the height of the free energy barrier and the so-called transmission coefficient.1 This product should give the same answer as a corresponding MD simulation provided someone has the patience to simulate sufficiently long to obtain reliable statistics with such an MD simulation (which is for our system of the order of microseconds!). Transition state theory assumes that all trajectories that start on top of the free energy barrier with a positive initial velocity will end up in the ìproductœ state. This implies that the transmission coefficient is 0.5 and the hopping rate can be estimated from the free energy pro–le only.Some preliminary tests have shown however that in our case far fewer trajectories end up in the product state. This implies that transition state theory for our choice of order parameter is not a good approximation and we have to perform MD simulations to compute the transmission coefficient. It may be possible that a diÜerent choice of order parameter gives a lower free energy barrier and hence a transmission coefficient that is closer to 0.5. 1 D. Frenkel and B. Smit in Understanding Molecular Simulations From Algorithms to Applications Academic Press Boston 1996. Prof. Jansen asked You have been entering the hydrocarbons into the zeolite channels by a non-physical mechanism. Does this aÜect the results ? For example the heights of the saddle points seem to be derived from a sequence of optimized con–gurations i.e.from local minima and therefore might be underestimated. How did you make sure that the hydrocarbons as a whole really can pass through the bottlenecks ? Prof. Smit answered It is correct that the growing scheme that is used in the con–gurational-bias Monte Carlo (CBMC) calculations is non-physical and hence the molecules enter the zeolite channel in an arti–cial way. It is important to note that the resulting distribution of con–gurations is the correct Boltzmann distribution only the way in which we obtain these con–gurations is not physical. If we were to use a physical scheme e.g. MD it would take many billion years of CPU time since experimentally it may take days before these systems are equilibrated.There are some zeolites however in which some of the cages are not accessible from the outside. The CBMC scheme like any other grand-canonical simulation technique would introduce molecules into these cages which would not be correct. If one knows this such artifacts can be avoided by –lling these cages with hard sphere molecules. For the systems that we have studied it is known from experiment that the molecules are adsorbed. Dr Freeman asked Could you comment on the feasibility of calculating free energy pro–les for activated events which account for sorbate»sorbate interactions ? Prof. Smit responded It is not straightforward to extend our calculations in which we assume in–nite dilution to cases in which sorbate»sorbate interactions are impor-132 General Discussion tant.An alternative approach would be the use of a lattice on which molecules can hop using the hopping rates calculated by our approach. My earlier response to Prof. Bliek is also relevant to your question. Prof. Gillan commented For simple defects in solids there are well established methods for calculating free energy pro–les and barriers. These involve doing constrained simulations with the reaction coordinate –xed at a sequence of values and hence calculating the potential of mean force. Is there any reason why such methods should not be applied to molecules in zeolites ? Prof. Smit responded No there is no reason why the methods you mentioned could not be applied to these systems. My earlier response to Prof.Islam is also relevant to your question. Prof. Sauer said Years ago Klaus Fiedler from the Central Institute of Physical Chemistry in Berlin performed MC simulations on adsorption of alkane chains and aromatic hydrocarbons in zeolites.1 He used a special importance sampling procedure. How does your con–gurational-bias MC diÜer from his approach? 1 B. Grauert K. Fiedler H. Stach and J. Jaé nchen in Zeolites Facts Figures Future ed. P. A. Jacobs and R. A. van Santen Elsevier Amsterdam 1989 pp. 701»713. Prof. Smit replied The starting point of both methods is the same molecules adsorbed in a zeolite and the interactions of these molecules with each other and with the zeolite are described with a given classical potential. The behaviour of this system is fully described by the partition function i.e.once this partition function is known all ensemble averages can be computed from this function. For most systems however this partition function cannot be computed exactly and Fiedler and co-workers have developed an MC scheme to approximate this partition function. In our CBMC scheme we cannot compute a partition function we can only calculate ensemble averages. However the calculation of these ensemble averages does not involve additional approximations. Prof. Cheetham opened the discussion of Dr Yashonathœs paper You have presented a large number of MD simulations relating to the so-called levitation eÜect. Is there any experimental validation of this eÜect ? Do you have candidate systems for the separations that you believe might be facilitated by the levitation eÜect ? Dr Yashonath answered Experimental evidence for the existence of this eÜect is something that is of considerable importance.To the best of my knowledge there have not been any experiments carried out subsequent to 1994 when the existence of an anomalous peak in the diÜusion coefficient was –rst reported based on molecular dynamics simulations.1 Among the problems in designing experiments to verify this eÜect are the in—uence of the grain boundary which can be overcome by the use of single crystals. Again there is the in—uence of the extra-framework cations and water but most of these can be circumvented. Assuming that these have been overcome it is possible to use two diÜusing species of diÜering size and then study their diÜusion in a zeolite or any other porous media of known structure.There are many porous solids whose structures are well known and this should not be a problem. However one diÜusing species should have a size with c@1 which corresponds to the linear regime while the other should have a size with cB1 (these shall be referred to as type I experiments). For simplicity we take two species which are spherical in shape. The problem is that in practice it is rarely possible to –nd two diÜerent sized sorbates without some change in other characteristics such as 133 General Discussion the mass or the interaction strength. Corrections for the change in mass could be applied but it is not clear how the changes in the interaction strength would alter the diÜusion though we have recently reported some preliminary results.2 For these reasons this design of the experiment does not seem appropriate.There is an alternative method which is probably better suited for laboratory experiment (termed type II experiments). This involves choosing one diÜusing species and two diÜerent hosts with diÜering void diameters in such a way that one void size corresponds to c@1 and the other cB1. Zeolites (preferably highly siliceous) and nanotubes may be ideal candidates. Note that the mass and even the interaction strength are unaltered if both hosts are of the same material. Yet another alternative method is to introduce one guest in a host material in which there are channels of diÜering diameters along diÜerent directions.If one could obtain from experiment the components of the diÜusion coefficient tensor then it would be possible to see if the levitation eÜect exists. The calculation on xenon3 in which you found that xenon diÜuses faster in the 10-membered channels of ferrierite as compared to the 12-membered channels of zeolite L provides a clue. If one could –nd a zeolite consisting of 10-membered channels along one direction and 12-membered cylindrical channels along another direction this would provide a good system to test the eÜect. Finally the problems with measuring the self-diÜusion coefficient have to be tackled since many of the techniques such as NMR neutron scattering etc. yield information about the diÜusion coefficient over a narrow range of D and place additional restrictions on the nature of the sorbate that can be used.It is gratifying to note that the range over which anomalous diÜusion is observed is not that narrow it extends from 5 to 7 ” in zeolite Y. However it should be noted that the intensity of the peak falls oÜ rapidly and hence one prefers to choose a sorbate»host system so as to be as close to the peak as possible. This places more stringent restrictions. In the case of zeolite Y if one were to be in the top half of the intensity then this range is narrowed to 5.3»6.3 ”. 1 S. Yashonath and P. Santikary J. Chem. Phys. 1994 100 4013. 2 B. Bhattacharjee and S. Yashonath Mol. Phys. 1997 90 889. 3 N. J. Henson A. K. Cheetham B. K. Peterson S.D. Pickett and J. M. Thomas J. Comput. Aided Mater. Des. 1993 1 41. Prof. Smit asked The levitation eÜect occurs if the size of the molecule matches the size of a channel. This occurs especially at low temperatures and with a narrow range of the diameter. Does this eÜect still exist if a —exible zeolite lattice is used? Dr Yashonath responded Yes the eÜect manifests noticeably only at low enough temperatures and in a fairly narrow diameter range. Please see the answer to Prof. Cheethamœs question (above) where we have pointed out that the range is not really that narrow. The eÜect still exists when the framework of the zeolite is modelled as a —exible one and the evidence for this comes from simulations that we have done in at least two zeolitic systems. First in the case of zeolite A it was found that the eÜect persisted when the framework of the zeolite was modelled as a —exible cage.1 There was however some decrease in the intensity of the anomalous peak. Secondly we studied a silicalite system. Here results for the rigid cage have also been reported.2 Even in this silicalite system the peak in the diÜusion coefficient persisted. These results seem to suggest that the levitation eÜect does exist when the zeolite framework is modelled as —exible. 1 P. Santikary and S. Yashonath J. Phys. Chem. 1994 98 9252. 2 S. Bandyopadhyay and S. Yashonath J. Phys. Chem. 1995 99 4286.