Lattice energies of small organic molecules using density functional theory
Organic crystals are of particular interest in pharmaceutical industry, because the majority of drugs are produced and distributed in crystalline form. Molecular crystals play also an important role as agrochemicals, pigment dyes, and organic photovoltaics . The applicability and also properties like the solubility depend very much on the stability of these compounds. It is therefore crucial to understand the factors that affect the stability. For comparing the stability of two or more molecular crystals that are, for example, based on the same molecular scaffold, but contain different functional groups, lattice energies can be calculated to predict the energetically most favored, and thus, thermodynamically most stable compound.
Usually, organic molecules crystallize in more than one crystal structure – a phenomenon which is referred to as polymorphism . A well-known example is paracetamol for which two different crystalline modifications are shown in Figure 1
For industrial applications, it is of great interest to understand polymorphism and to know the crystal structure of a compound, since many important chemical and physical properties depend on the structure. Computational methods can help to predict and identify the most stable polymorphic form based on the calculation of lattice energies .
In this application note, we consider a set of pharmaceutically relevant compounds, namely aspirin, two polymorphic forms of paracetamol, and one modification of ibuprofen which will be compared based on the relative stability obtained from lattice energy calculations. The molecular structures of these compounds are shown in Figure 2.
Since especially the energy differences between two polymorphic forms are typically rather small, accurate methods are needed for reliable predictions of the lattice energies. In the following, ab initio calculations are presented which have been performed using the ABINIT  module of Scienomics MAPS platform .
Geometry optimizations for relaxing the atomic positions were performed within the framework of the density functional theory (DFT) with the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) parametrization together with Fritz-Haber-Institute (FHI) pseudopotentials. Initial single-point calculations for aspirin had shown that relative energies are converged within 1.3 kJ/mol using an energy cutoff of 600 eV and within 0.6 kJ/mol using a 2x2x2 k-point Monkhorst-Pack grid for Brillouin zone integration. Subsequent to the geometry relaxation, the semi-empirical dispersion correction DFT-D2 of Grimme  has been calculated for the optimized structures.
The optimized crystal structure is shown in Figure 3 together with a plot illustrating the energy behavior during the optimization.
The lattice energies are given in Table1 and have been calculated according to Elatt = Ecryst / Z - Emol with Ecryst as the energy of the crystal, Z the number of molecules in the unit cell, and Emol as the energy of the single molecule. Emol was obtained by placing a single molecule into a large cubic unit cell with the box length a = 15 Å. In addition to the calculated values, experimental sublimation energies (taken from Ref. 6 and reference therein), which represent a measure for the lattice energies, are listed for comparison in Table 1, as well.
Van der Waals-type interactions make up an essential part of the intermolecular interactions in organic solids. This is also reflected in the calculated lattice energies. When comparing the lattice energies with and without dispersion correction, it can be seen that the dispersion energy contributes 56 - 80 % to the lattice energy. The predicted ranking and relative energies of aspirin, paracetamol I, and (±)-ibuprofen at DFT-D2 level are in excellent agreement with the corresponding experimental sublimation enthalpies, which are approximately related to the lattice energy by Hsub= -(Elatt + 2RT)
with R as the gas constant and T the temperature. In case of paracetamol, the ranking of the two polymorphs at DFT-D2 level is contrary to the experimental data, whereas without dispersion correction the ranking is the same as obtained in experiment. This discrepancy at DFT-D2 level may be attributed to the fact that the dispersion correction was not applied during optimization and that only the atomic positions have been relaxed, while the cell as not been fully optimized. The calculated results are, however, in line with previous computational studies , which also found a higher lattice energy for form II using an empirical dispersion correction and an Gaussian type orbital-based approach. The latter suffers from basis set superposition errors in contrast to plane wave-based approaches, which have been used in the present case. Another reason may be the neglect of entropic effects, whose evaluation is, however, beyond the scope of the this study.
In summary, the results demonstrate that the calculation of lattice energies is a useful approach to assess the thermodynamic stability of organic molecular crystals, in which dispersion-type interactions make up an essential part of the intermolecular interactions. Relative energies based on the semi-empirical dispersion-corrected lattice energies are in good agreement with the experimentally observed trend rendering the chosen methodology as a suitable and reliable approach for predicting the thermodynamic stability of organic crystal solid.
- G. Desiraju, J. Chem. Sci., 122 (2010) 667-675 .
-  A. Nangia, Acc. Chem. Res., 41 (2008) 595-604.
- S. Woodley et al., Nature Materials, 7 (2008) 937-946.
- X. Gonze et al., Computer Phys. Commun., 180 (2009) 2582-2615.
- MAPS, Version 3.3.2, Scienomics, Paris, France.
- T. Li et al., Pharmaceutical Research, 23 ( 2006), 2326-2332.
- S. Grimme, J. Comp. Chem., 27 (2006), 1787-1799.
- J.S. Chickos et al., J. Phys. Chem. Ref. Data, 31 (2002 ) 537-698.