I. INTRODUCTION
Quantum-chemical calculations have been successfully used to complement the experimental powder X-ray diffraction (PXRD) data at several stages of the crystal structure solution process: (1) Optimization of the molecular geometry in the gas phase to obtain more accurate starting molecular structure models suitable for structure solution by real-space methods. At this stage, different levels of theory can be applied: molecular mechanics (MM), semi-empirical (SE) methods, Hartree–Fock methods, density functional theory (DFT). (2) Minimization of the energy of crystal structure and estimate of reasonable bond distances and angles to be used as chemical restraints in the Rietveld refinement. (3) Energy minimization as the most reasonable approach to compute optimal positions for the H atoms (the accurate location of the H position is a difficult task from PXRD data) to be refined, with unit cell and non-H atomic positions, by the Rietveld method. (4) Selecting the correct structure from a set of potential structures that differ in the torsion angle rotation of functional group (Naelapää et al., Reference Naelapää, van de Streek, Rantanen and Bond2012) or in the space group, and could not be distinguished on the basis of the PXRD data alone. (5) Geometry optimization with dispersion-corrected DFT (DFT-D) approach in the solid state for refining crystal structure and providing results whose accuracy is comparable with that of single crystal, especially when, the Rietveld refinement yields inaccurate molecular geometry (Smrcok et al., Reference Smrcok, Mach and Le Bail2013). (7) Assessment of the correctness of experimental crystal structures (van de Streek and Neumann, Reference van de Streek and Neumann2010).
Many packages for quantum-chemistry calculations support the application to periodic systems, expanding the electron orbitals in crystalline solid either in periodic combination of atomic orbital basis sets modulated with plane waves [e.g., GAUSSIAN (Frisch et al., Reference Frisch, Trucks, Schlegel, Scuseria, Robb, Cheeseman, Scalmani, Barone, Mennucci, Petersson, Nakatsuji, Caricato, Li, Hratchian, Izmaylov, Bloino, Zheng, Sonnenberg, Hada, Ehara, Toyota, Fukuda, Hasegawa, Ishida, Nakajima, Honda, Kitao, Nakai, Vreven, Montgomery, Peralta, Ogliaro, Bearpark, Heyd, Brothers, Kudin, Staroverov, Keith, Kobayashi, Normand, Raghavachari, Rendell, Burant, Iyengar, Tomasi, Cossi, Rega, Millam, Klene, Knox, Cross, Bakken, Adamo, Jaramillo, Gomperts, Stratmann, Yazyev, Austin, Cammi, Pomelli, Ochterski, Martin, Morokuma, Zakrzewski, Voth, Salvador, Dannenberg, Dapprich, Daniels, Farkas, Foresman, Ortiz, Cioslowski and Fox2013), CRYSTAL (Dovesi et al., Reference Dovesi, Orlando, Erba, Zicovich-Wilson, Civalleri, Casassa, Maschio, Ferrabone, De La Pierre, D'Arco, Noel, Causa, Rerat and Kirtman2014)] or in plane wave basis sets alone [e.g., VASP (Kresse and Hafner, Reference Kresse and Hafner1993), ABINIT (Gonze et al., Reference Gonze, Amadon, Anglade, Beuken, Bottin, Boulanger, Bruneval, Caliste, Caracas, Cote, Deutsch, Genovese, Ghosez, Giantomassi, Goedecker, Hamann, Hermet, Jollet, Jomard, Leroux, Mancini, Mazevet, Oliveira, Onida, Pouillon, Rangel, Rignanese, Sangalli, Shaltaf, Torrent, Verstraete, Zérah and Zwanziger 2009 ), NWCHEM (Valiev et al., Reference Valiev, Bylaska, Govind, Kowalski, Straatsma, van Dam, Wang, Nieplocha, Apra, Windus and de Jong2010), QUANTUM ESPRESSO (Giannozzi et al., Reference Giannozzi, Baroni, Bonini, Calandra, Car, Cavazzoni, Ceresoli, Chiarotti, Cococcioni, Dabo, Dal Corso, Fabris, Fratesi, de Gironcoli, Gebauer, Gerstmann, Gougoussis, Kokalj, Lazzeri, Martin-Samos, Marzari, Mauri, Mazzarello, Paolini, Pasquarello, Paulatto, Sbraccia, Scandolo, Sclauzero, Seitsonen, Smogunov, Umari and Wentzcovitch2009), CPMD (CPMD, 1990)]. Several packages are freely available or can be acquired for academic use at relatively low cost. The use of these programs is in general quite easy and requires the preparation of input data files where the user supplies the atomic positions (Cartesian or fractional coordinates), crystallographic data (e.g., unit-cell dimensions, space group, number of atoms per unit cell or in the asymmetric unit), aspects related to the computational method (e.g., density functional, basis set), and type of calculation (e.g., geometry optimization, vibrational analysis, single-point energy). These programs generate output information, which can be analysed, and the final molecular structure, which can be displayed by many programs.
EXPO2014, the evolution of EXPO2013 (Altomare et al., Reference Altomare, Cuocci, Giacovazzo, Moliterni, Rizzi, Corriero and Falcicchio2013), is a computer program which is able to carry out almost automatically all the steps of the structure solution by PXRD data. It includes new additional features to automate the process of creation of input files and elaboration of the output files created by the quantum-chemistry packages. Quantum-chemistry programs have parallel computing capabilities and can be run on platforms ranging from modern laptops to supercomputers. A typical not expensive hardware could be a multi-core Linux workstation and the time required for the calculation is in the order of magnitude of 100 h for a DFT-optimization on small molecules on a single CPU.
PXRD is a powerful and rapid analytical method used for phase identification of crystalline materials, for providing information about unit cell, space groups (Bergmann et al., Reference Bergmann, Le Bail, Shirley and Zlokazov2004), and structure model. Powder diffraction also allows for rapid, non-destructive analysis of multi-component mixtures without needing extensive sample preparations. Unfortunately, some problems (diffraction peak overlap and broadening, preferred orientation, and background) make the structure solution process by PXRD a not straightforward task. Ab initio structure solution from powder data is still now a challenge in many cases.
In EXPO2014, either ab initio procedures by Direct Methods or real-space techniques can be used to deal with PXRD structure solution (David et al., Reference David, Shankland, McCusker and Baerlocher2002; Dinnebier and Billinge, Reference Dinnebier and Billinge2008). Direct Methods are very effective and fast when high-quality data are available, especially synchrotron radiation data, and heavy atoms are present. Their complexity rises with the number of atoms. On the other hand, real-space methods can be applied on large molecular system with a complexity bound to the number of internal degrees of freedom. One of the drawbacks of the real-space methods is the necessity to input an accurate 3D (three-directional) structural model of the asymmetric unit with bond distances, bond angles and ring conformations that should match as closely as possible that of the correct and unknown crystal structure. A starting model for real-space methods can be prepared following essentially two different approaches: (1) searching for similar crystal structure in database [e.g., Cambridge Structure Database (CSD) (Groom et al., Reference Groom, Bruno, Lightfoot and Ward2016), Inorganic Crystal Structure Database (ICSD) (Belsky et al., Reference Belsky, Hellenbrandt, Karen and Luksch2002), Crystallography Open Database (COD) (Grazulis et al., Reference Grazulis, Chateigner, Downs, Yokochi, Quiros, Lutterotti, Manakova, Butkus, Moeck and Le Bail2009), Pauling File (Villars et al., Reference Villars, Berndt, Brandenburg, Cenzual, Daams, Hulliger, Massalski, Okamoto, Osaki, Prince, Putz and Iwata2004), ICDD databases with atomic coordinates: PDF-4 + , PDF-4/Minerals, PDF-4/Organics (Faber and Fawcett, Reference Faber and Fawcett2002)], (2) performing a geometry optimization by using quantum-chemistry packages on the basis of three different levels of quantum-chemistry theories described on the next section, for each of them EXPO2014 provides a specific support. But the use of theoretical calculations is not limited to this preliminary step and can be applied also in the following steps of the process of crystal structure determination from powders.
Experimental powder diffraction data, especially standard laboratory data, could have low/medium-quality information because of preferred orientation, severe reflection overlap, low-resolution, absence of heavy atoms, and rapid decay of scattering power (organic compounds). The quality of data can be sufficient to locate heavy atoms or significant part of molecule by Direct Methods, or to position all molecular fragments in the cell by real-space techniques, but not enough to refine the structure. The refinement of structure can require the use of restraints and constraints to stabilize the refinement especially when the number of refined parameters is large. The hydrogen atoms are normally placed in calculated positions during structure refinement calculations rather than explicitly refined, leaving a great uncertainty in the position of terminal H-atoms (e.g., hydrogens in methyl group, hydroxyl group, and aminic group), frequently important in the definition of hydrogen bond network. In this scenario, the DFT-optimization can be used as an independent source of information for obtaining accurate bond distances and angles that should be used as restraints in the Rietveld analysis, and in an accurate location of hydrogen atoms. At the end of the structure refinement, the DFT-D method can be used to assess the correctness of the experimental crystal structure.
II. EXPERIMENTAL
A. Geometry optimization through MM force field
In MM, the forces acting on atoms are not described on the basis of a solution of the electron wavefunction, but of a simplified geometrical model responding to chemical knowledge of fragments of the structure, with parameters derived either from experiment or from quantum-chemistry calculations. Such models and their parameters are named “force fields” and are usually specific to compound classes, though a few of them attempt to cope with all the periodic table.
The total potential energy is typically taken to be the sum of the bond stretching energy E str, the bending energy E bend, the twisting (or torsion) energy E tor, and the energy of interaction between non-bonded atoms E non-bond. The last contribution includes van der Waals, steric, and electrostatic interactions between atoms not chemically bound.
$$E_{{\rm tot}} = E_{{\rm str}} + E_{{\rm bend}} + E_{{\rm tor}} + E_{{\rm non - bond}} . $$
To optimize the molecular geometry and find the most stable conformation, each of the forces f i acting on a nucleus must vanish:
$$f_i = -\partial E_{{\rm tot}} /\partial q_i = 0,$$
where q i are the nuclear coordinates of the ith atom.
In EXPO2014, two MM force fields can be used: (1) the Merck molecular force field (MMFF) for organic molecules and biomolecules (Halgren, Reference Halgren1996); (2) the universal force field (UFF), which accounts for all the periodic table. This is widely used for system containing inorganic elements (Rappe et al., Reference Rappe, Casewit, Colwell, Goddard and Skiff1992). In EXPO2014, the support for MM is provided using the class OBForceField from the library Open Babel (O'Boyle et al., Reference O'Boyle, Banck, James, Morley, Vandermeersch and Hutchison2011). The geometry optimization by MM can be carried out by the EXPO2014 menu “Optimize Geometry” (Figure 1). Clicking on “Tools” > “Optimize Geometry” > “Optimize”, the MMFF94 force field will be applied and, in case of failure (e.g., inorganic elements), an UFF calculation will be executed. The initial molecular geometry can always be restored clicking on “Modify” > “Undo”. The geometry optimization can also be applied only on an atom selection (“Tools” > “Optimize Geometry” > “Optimized Selected Atoms”). For example, if one is interested in optimizing only the position of H atoms and fixing the position of all other atoms you can choose “Select” > “Select Type” > “H atoms” and then “Tools” > “Optimize Geometry” > “Optimize Selected Atoms”.
Menu for geometry optimization in EXPO2014.
The two salient features of MM calculations on small to medium-sized molecules are that they are fast and they can be very accurate, especially on organic compounds, for which the force fields are very detailed. MM calculations might still take only seconds. MM, even when it does not provide very accurate geometries, can supply reasonably good structures for other kinds of calculations: SE, ab initio or DFT; and this is undoubtedly one of its main application.
B. Geometry optimization with Molecular Orbital PACkage (MOPAC)
MM force fields do not include an explicit treatment of the electrons and therefore do not well account for a variety of effects. In particular, the molecular potential energy cannot always be well described as a sum of bond and angles terms, and the parametrization of forces for elements with a wide variance of nearest-neighbour geometries or bond lengths may be unfeasible.
The determination of the electron structure is typically a quantum problem, the solution of which is numerical and subject to different possible levels of approximation. As usual, in quantum mechanics, it starts with the Schrödinger equation:
$$H\Psi = E\Psi, $$
where H, the Hamiltonian operator, summarizing the kinetic energy, all the particle–particle interactions inside a molecule and interaction with the external field (usually because of other molecules, like in the condensed state). The lowest energy state can be estimated by the minimization of the energy associated with the wave function. The wavefunction is a complicated many particle function, usually described in terms of orbitals, single-particle functions composed to yield a correct form for the wavefunction. Nowadays, the most profitable descriptions of the molecular energy, in the context of geometry resolution are two: (1) a model simplified Hamiltonian with parameters derived from experimental datasets and quantum calculations of high accuracy, denoted as SE, fast, and especially accurate for compounds belonging to the same classes as those found in the experimental datasets. (2) The density functional energy, describing all the energy in terms of the orbitals and their associated electron density, much slower than SE methods but covering all molecules.
MOPAC (MOPAC2016, 2016) is one of the most widely used semi-empirical software packages and is designed for a wide range of functionality. We tested version 2012 and 2016. MOPAC is not distributed with EXPO2014, so a copy of the software needs to be separately downloaded. For more information or to obtain a copy, please visit http://www.openmopac.net/MOPAC2016.html, the software is currently free for academic use.
An interface to MOPAC is provided in EXPO2014 and allows users to perform a gas-phase geometry optimization. This approach is useful to generate molecular models with accurate angles and bond lengths, suitable as a starting model for structure solution by real-space methods. The advantage of SE calculations is that they are much faster than self-consistent field method (SCF) or DFT methods and can be applied routinely to large system. A MOPAC optimization is performed selecting “Tools” > “Optimize Geometry” > “MOPAC” from the principal menu. The calculation is performed by choosing some options in the MOPAC dialogue window (Figure 2). This dialogue window allows the user to: (1) define the title of the calculations, e.g., the name of molecule. It is then reported in the second line of the MOPAC input file. (2) Specify the SE Hamiltonian used in the electronic part of the calculation to obtain molecular orbitals and the associated energy. Not all chemical elements are supported by every selected Hamiltonian. (3) Write the charge of the system under study. (4) Allow the user to specify the number of singly occupied orbitals. Singlet is the default, and specifies a closed shell (most organic and low atomic weight molecules); doublet specifies one singly occupied orbital (e.g., organic radicals); triplet specifies two singly occupied orbitals; and so on. If the multiplicity is unknown, trying each sensible case is advisable. (5) Specify if the molecular geometry is described using the Cartesian coordinate definition or the Gaussian Z-matrix format. The position of molecule in the cell unit will be lost if the geometry is defined in Z-matrix format. (6) Specify the folder where the MOPAC executable file is located. EXPO2014 automatically locates the position where MOPAC exists if MOPAC is installed in a standard directory (C:/Program Files/MOPAC/on Windows, /opt/mopac/on Linux) or the environment variable PATH includes the installation directory; otherwise the MOPAC executable must be specified directly by the user. By selecting Submit, MOPAC will be started and at the end of the geometry optimization, the new geometry will be displayed and an output file will be produced in the working directory. The output file can be viewed inside EXPO: (“Tools” > “Optimize Geometry” > “View MOPAC Output”). Previously generated MOPAC output files can be read and the molecular geometry can be displayed through the menu “File” > “Import Structure” and selecting “MOPAC Output File (*.out)” in the file extension menu.
(Colour online) Dialogue for generating input file for MOPAC.
It is important to note that these calculations require the inclusion of hydrogen atoms that are normally excluded in structure solution calculations from PXRD. Selecting “Add Hydrogens” in the menu “Tools”, it is possible to assign sensible positions of H-atoms. Wrong starting geometry, missing hydrogens or wrong multiplicity are possible causes of incorrect result. In this situation, the use of “Modify” > “Undo Optimize Geometry” to restore the initial geometry is suggested. Then the output file can be read to find a possible cause of error. Additional calculations can be performed by clicking the Preview button and changing the input file. For example, the keyword NOOPT can be used to perform only an energy calculation; the calculated energy will be written in the output file. The use of the keyword NOOPT OPT-H is devoted to optimize the positions of hydrogen atoms while freezing the positions of all other atoms.
C. Geometry optimization by DFT
An approach that has gained increasing prominence in recent years is DFT, in particular in the Kohn–Sham formulation (Kohn, Reference Kohn1999; Sholl and Steckel, Reference Sholl and Steckel2009). This is the only formulation for which there is ample coverage in the mentioned quantum-chemistry programs. The electron density is expressed in terms of a set of orbitals, and the energy is obtained from functionals of the density matrix (kinetic energy) or of the electron density alone. Some part of this, usually denoted E XC[ρ], is fruit of approximation and several prescriptions for it can be found in the DFT literature. The correct choice of the density functional is important to the calculation of several molecular properties, but not so much in the geometry optimization, except when two geometries must be distinguished which differ in energy by a small value, typically <20 kJ mol−1.
DFT is rather efficient; the computational suites optimize the code to scale almost linearly with the number of electrons, so as to make affordable computations with a large number of electrons on modern hardware. It can be used for crystalline solids with hundred atoms per unit cell, which is a sufficient number for the application in the field of powder diffraction. EXPO2014 is able to create input files for well-known quantum-chemistry programs, selecting “File” > “Import Structure” from its main menu. Input files for GAMESS-US (Schmidt et al., Reference Schmidt, Baldridge, Boatz, Elbert, Gordon, Jensen, Koseki, Matsunaga, Nguyen, Su, Windus, Dupuis and Montgomery1993) and NWChem contain keywords for geometry optimization by DFT calculations to perform using the B3LYP functional (Stephens et al., Reference Stephens, Devlin, Chabalowski and Frisch1994) and a standard 6–31G* basis set (Hehre et al., Reference Hehre, Ditchfield and Pople1972).
One of the most widely used applications of the DFT is the validation of experimental molecular crystal structures. This type of calculation is generally performed by using plane-wave DFT with GGA (generalized gradient approximation) functionals such as PBE (Perdew et al., Reference Perdew, Burke and Ernzerhof1996) and PW91 (Perdew and Wang, Reference Perdew and Wang1992). The long-range dispersive interactions (London or attractive van der Walls forces), which play an incredibly important role in determining the structure and the properties of molecular crystal and are not included in the DFT Hamiltonian, are described applying specific correction schemes such as DFT-D3 (Grimme et al., Reference Grimme, Antony, Ehrlich and Krieg2010) and TS (Tkatchenko and Scheffler, Reference Tkatchenko and Scheffler2009). The procedure is simple: a geometry optimization is performed using, as starting molecular geometry, the experimentally determined coordinates of the crystal structure. The final calculated crystal structure is compared with the experimental one. By using “Tools” > “Optimize Crystal Structure by DFT-D”, EXPO2014 can create input files for the optimization of crystal structure by DFT-D3 with functional PBE. Currently, the supported software are NWChem, ABINIT, CRYSTAL, and QUANTUM ESPRESSO. The input files generated for these software can be directly used to run the crystal structure validation or can require only minor modifications. At the end of the DFT-optimization, the output file generated by the quantum-chemistry software can be imported in EXPO2014 and the calculated crystal structure can be visually compared with the experimental one by using “Tools” > “Overlay Structures” to generate an overlaid (Figure 3). The quality of the agreement between the experimental and the theoretical structure can be evaluated in terms of root-mean-square displacement (RMSD) calculated for all non-H atoms. According to a recent paper (van de Streek and Neumann, Reference van de Streek and Neumann2010), the value of RMSD never exceeds 0.25 Å for correct crystal structure. RMSD information is also available in EXPO2014 by “Tools” > “Overlay info”.
Example of overlaid between experimental (black) and optimized (grey) crystal structures of ibuprofen. The optimization was performed by using DFT-D3 in NWChem.
III. CONCLUSIONS
EXPO2014 has recently added features to support the increasing application of quantum-chemistry calculations as complementary tools to the structure determination from powder diffraction data. EXPO2014 can assist the user in preparing input files for quantum codes such as MOPAC, GAMESS-US, NWChem, GAUSSIAN, ABINIT, CRYSTAL, and QUANTUM ESPRESSO. Input files can be directly run for geometry optimization of molecules or to automate the procedure of crystal structure validation. Only for MOPAC a graphic dialogue has been developed and present the basic features required to run basic quantum calculations. In addition, a quick optimization through the force fields provided by the Open Babel can be performed.
ACKNOWLEDGEMENTS
The authors thank Dr. Caterina Chiarella of the Institute of Crystallography (Bari, Italy) for technical support.
