The structure and energetics of glycine polymorphs based on first principles simulation using density functional theory

摘要

The structure and energetics of glycine polymorphs based on first principles simulation using density functional theory C. M. Freeman,a J. W. Andzelm,a C. S. Ewig,a J.-R. Hilla and B. Delleyb a Molecular Simulations Inc., 9685 Scranton Road, San Diego, CA 92121, USA. E-mail: clive@msi.com b Paul Scherrer Institut, Zurich, Badenerstr. 569, CH-8048, Zurich, Switzerland Received (in Cambridge, UK) 3rd August 1998, Accepted 2nd October 1998 Density functional theory (DFT) calculations for the crystal structures of polymorphs of the glycine zwitterion are reported; with unit cell parameters constrained at experimentally determined values, energy minimized configurations for three known glycine polymorphs are in good agreement with crystallographically determined structures, and the calculated energies are in qualitative agreement with observed lattice stabilities.Molecular orbital calculations, employing Hartreendash;Fock,1 density functional theory (DFT),2 or even an appropriate combination of both methods3 are standard techniques for the structural and energetic description of molecular systems.Significant developments have also recently taken place in the use of first principles methods to represent ionic and covalently bound periodic systems, such as silicon and zeolites, typically employing plane-wave and plane-wave/pseudo-potential approaches. 4ndash;7 However, the simulation of molecular solids using first-principles methods has received considerably less attention8 and the use of such methods to address the relative energies of crystal polymorphs less still.The lack of suitably refined protocols underlies this comparatively low level of application. Nevertheless, the results presented here show that such calculations are now within the compass of contemporary quantum mechanical methods and computational resources. To explore the ability of periodic DFT method calculations to describe molecular crystals, calculations were conducted targeting the known polymorphic structures of the glycine zwitterion and using DMol3.9ndash;11 The level of theory employed has been shown to be effective in the description of isolated molecules. 12dagger; Initial investigations addressed the selection of a suitable quantum mechanical Hamiltonian and basis set. These calculations employed the b-glycine polymorph, which possesses a unit cell containing two symmetry related molecules and a total of 20 atoms.13 Here results for the BLYP14ndash;17 functional are summarized.A minimal basis resulted in a 0.34 Aring; RMS deviation between minimized and observed structures, a double numeric basis (0.33 Aring;), a double numeric basis with polarization functions on heavy atoms (0.22 Aring;) and double numeric basis with polarization functions on all atoms (0.21 Aring;). Of the basis sets sampled then, optimal structural agreement is achieved on full geometry optimization with a so-called double numeric basis set with polarization functions on all atoms (DNP).The DNP basis set has been demonstrated to provide an efficient route to molecular polarizabilities and charge distributions as a result of its reasonable representation of the tail of the wavefunction.10 However, the ability of this computational approach to describe the relative energetics of polymorphic crystal structures is perhaps of more interest. To address the calculation of such relative polymorph energies, crystallographically determined structural starting models for a- and g-glycine were employed.18,19 As an additional test of the simulation protocol two incorrectly packed arrangements of glycine molecules were used as starting points; S1, containing four molecules in the unit cell dimensions of the a-polymorph, and S2, containing two molecules with cell dimensions of the b-polymorph. Energetic results for the systems considered are collected in Table 1. The similarity of the relative energies obtained in the DFT calculations for the a-, b- and g-polymorphs is mirrored in the experimental observation that all polymorphs are formed at temperatures close to room temperature, with the a-polymorph having the greatest observed stability.18 The energy differences obtained in the simulations are small and entropic contributions to the lattice stability have not been considered in the present calculations.The inclusion of vibrational entropy would, of course, be possible with increased computation times through calculation of the systemrsquo;s dynamical matrix and appropriate integration of the resulting phonon spectrum.However, the fact that the relative energies of incorrectly packed molecular arrangements S1 and S2 are higher than the observed polymorph energies is of significant interest. For the polar b- and gpolymorph crystals such a procedure assumes that the macroscopic dipolar energy of the structures is not significant.20,21 This appears to be a reasonable approximation.22 Fig. 1 shows starting structures and optimized structures for each of the polymorphs investigated. Good agreement is Table 1 The relative BLYP energies, per glycine molecule, of the a-, b- and g-glycine polymorphs. Also tabulated are the relative energies of hypothetical packing arrangements S1 and S2 and the root mean squared (RMS) displacements between the experimental structure and energy minimized structure Polymorph a b g S1 S2 Energy/kcal mol21 0.0 2.3 2.0 5.8 15.9 RMS/Aring; 0.19 0.21 0.30 0.33a 1.20a a RMS between starting and energy minimized structure.Fig. 1 Structural agreement between energy minimized (dashed line) and experimental structures (solid line) for a-, b- and g-glycine polymorphs. Hypothetical packing arrangements S1 and S2 (dashed line) show comparison with starting points (solid line). Chem. Commun., 1998, 2455ndash;2456 2455obtained between simulation and experiment for the known crystal models.The lowest RMS displacement between optimized and experimental structure (0.19 Aring;) is obtained for the apolymorph, a structure determined using neutron diffraction data with the lowest crystallographic R value of the polymorphs (R = 0.032).19 The starting geometries of incorrectly packed arrangements, S1 and S2, are altered by the optimization procedure to a greater extent than the observed polymorphs. In neither case do the resulting structures resemble the experimental polymorph structures.For S1 the structural deviation has an RMS value of 0.33 Aring; implying that, in this case, energy minimization is able to locate a stationary point in the vicinity of the starting structure. S1 was generated using the unit cell and asymmetric unit of the apolymorph and applying the symmetry operations of the P21/c spacegroup. This procedure yields a packing configuration similar for half the molecules to that of the a-polymorph (spacegroup P21/n). For S2 larger structural changes are evident in both the RMS displacement (1.2 Aring;) and in Fig. 1. This polymorph was generated by randomly placing two molecules within the cell of the b-polymorph. The present calculations employ the experimental unit cell dimensions as a constraint. This reduces computation times through a reduction in the number of variables used to describe the system. However it is important to note that cell parameter data, through indexed powder diffraction patterns, are among the most readily obtained structural information for solid state materials.The current calculations, which focus on the optimization of the geometry within the unit cell, provide valuable information that augments experimental observation and leads to the full description of the crystal. The magnitudes of partial charges are key to the accuracy of models employing classical mechanics. Table 2 collects the Mulliken23 charges for each of the atoms of the glycine molecule in the calculations conducted.The calculations are performed within a P�1 triclinic simulation cell. However, for the a-, b- and g-polymorphs symmetry is maintained during the optimization, hence charges are reported for a single molecule only. S1 was produced with the P21/c spacegroup and also has four equivalent molecules, one of which is listed. For S2, symmetry is not present and charges for bh molecules are listed. Interestingly for the experimental polymorphs there is little difference between the atomic partial charges obtained, despite the differences in local environment that each crystal form necessarily imposes on its constituent molecules.For the hypothetical packing arrangements S1 and S2, differing charge distributions are exhibited on several atoms, most notably O1, O2 and N. Also tabulated in Table 2 are gas phase derived atomic charges using the same basis set and Hamiltonian for the structure of the unique molecule of the a-polymorph.The magnitude of charges on oxygen and nitrogen atoms of the isolated molecule are lower than in the crystal structure calculations. The unambiguous definition of atomic partial charges from supplied molecular orbitals is not straightforward. 24,25 In particular Mulliken analysis is sensitive to the choice of basis set. However, these findings, using a particular charge definition method23 and uniform basis set and calculation type, indicate that the practice of transferring partial charges from gas phase calculations to the condensed phase may lead to an inaccurate description of the charge distribution exhibited by molecules in the crystalline state.Of course, no such ambiguity affects the first-principles calculations where charge density and its variation with environment are both consequences of the molecular orbitals obtained in the solution of the Schrouml;dinger equation. The calculations demonstrate that energy minimization leading to a determination of the relative energetics of molecular crystals is practical using first principles methods.Dagger; Notes and references dagger; The Brillouin zones of the polymorphs were sampled at the G point only.A detailed report on the effect of the choice of Hamiltonian and basis set on the structures and energies for molecular solids is in preparation. Geometry optimization (varying all atomic coordinates within the fixed cell) required calculation times of a few days on a single processor of an Origin 200 SGI (180 MHz) workstation.Dagger; The computed structures are available from the authors by email. 1 A. Hinchliffe, Computational Quantum Chemistry, Wiley, Chichester, 1988. 2 R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, OUP, 1989. 3 A. D. Becke, J. Chem. Phys., 1993, 98, 5698. 4 E. Wimmer, J. Comput.-Aided Mater. Des., 1993, 1, 215. 5 F. Haase, J. Sauer and J. Hutter, Chem. Phys. Lett., 1997, 266, 397. 6 R. Shah, J. D. Gale and M.C. Payne, J. Phys. Chem., 1997, B101, 4787. 7 I. Stich, J. D. Gale, K. Terakura and M. C. Payne, Chem. Phys. Lett., 1998, 283, 402. 8 M. Marchi, J. Hutter and M. Parrinello, J. Am. Chem. Soc., 1996, 118, 7847. 9 B. Delley, J. Chem. Phys., 1990, 92, 508. 10 B. Delley, J. Phys. Chem., 1996, 100, 6107. 11 DMol3 (Version Cerius2, 3.5), Molecular Simulations Inc., 9685 Scranton Road, San Diego, CA 92121, USA, 1997. 12 B. B. Laird, R. B. Ross and T. E. Ziegler, Chemical Applications of Density Functional Theory, ACS Symp.Ser. No 629, ACS, Washington, 1996. 13 Y. Iitaka, Acta Crystallogr., 1960, 13, 35. 14 C. Lee and C. Sosa, J. Chem. Phys., 1994, 100, 9018. 15 C. Lee, W. Yang and R. G. Parr, Phys. Rev., 1988, B37, 785. 16 R. Colle and D. Salvetti, Theor. Chim. Acta, 1975, 37, 329. 17 A. D. Becke, Phys. Rev. 1988, A38, 3098. 18 Y. Iitaka, Acta Crystallogr., 1961, 14, 1. 19 P.-G. Jonsson and A. Kvick, Acta Crystallogr., 1972, B28, 1827. 20 M. W.Deem, J. M. Newsam and S. K. Sinha, J. Phys. Chem., 1990, 94, 8356. 21 J. L. Derissen, P. H. Smit and J. Voogd, J. Phys. Chem., 1977, 91, 1474. 22 B. P. van Eijck and J. Kroon, J. Phys. Chem. B, 1997, B101, 1096. 23 R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833. 24 K. B. Wiberg and P. R. Rablen, J. Comput. Chem., 1993, 14, 1504. 25 J. Meister and W. H. E. Schwartz, J. Phys. Chem., 1994, 98, 8245. Communication 8/06102E Table 2 Mulliken (ref. 23) derived partial electronic charges obtained with the BLYP Hamiltonian and DNP basis sets for a-, b- and g-glycine polymorphs and hypothetical packing arrangements S1 and S2 (molecules 1 and 2). Also tabulated are partial electronic charges for an isolated molecule in the configuration of the unique molecule of the a-glycine polymorph Polymorph Isolated molecule a b g S1 S2: 1 S2: 2 O1 20.534 20.640 20.618 20.654 20.630 20.687 20.597 C1 0.401 0.553 0.555 0.546 0.530 0.550 0.514 N 20.222 20.410 20.456 20.425 20.494 20.410 20.447 C2 20.167 20.168 20.185 20.146 20.129 20.185 20.183 H1 0.234 0.352 0.356 0.348 0.357 0.269 0.331 H3 0.258 0.286 0.334 0.293 0.320 0.402 0.329 H2 0.265 0.344 0.344 0.353 0.351 0.260 0.283 H5 0.137 0.133 0.149 0.165 0.118 0.149 0.184 O2 20.509 20.650 20.654 20.629 20.622 20.538 20.606 H4 0.137 0.200 0.174 0.147 0.198 0.196 0.184 2456 Chem Commun., 1998, 2455ndash;2456