Davide
Presti
a,
Alfonso
Pedone
*a,
Maria Cristina
Menziani
a,
Bartolomeo
Civalleri
b and
Lorenzo
Maschio
*b
aDipartimento di Scienze Chimiche e Geologiche and INSTM research unit, Università di Modena e Reggio Emilia, Via G. Campi 183, I-41125 Modena, Italy. E-mail: alfonso.pedone@unimore.it; Fax: +39 059 373543; Tel: +39 059 2055043
bDipartimento di Chimica, Università di Torino and NIS Centre of Excellence, Via P. Giuria 7, I-10129 Torino, Italy. E-mail: lorenzo.maschio@unito.it; Fax: +39 011 670 7855; Tel: +39 011 6707564
First published on 18th October 2013
The molecular crystal of oxalyl dihydrazide differentiates into five polymorphs that are governed by inter- and intramolecular hydrogen bonds. The complex mixture of such interactions with long range dispersive forces makes its computational characterization very challenging; thus it represents an ideal benchmark for ab initio methods when striving for a description of polymorphism in molecular crystals. Indeed, a complete experimental energetic profile of this system is still lacking, and it is here investigated by means of periodic dispersion-corrected DFT and Local second order Møller–Plesset Perturbation theory (LMP2) calculations. In this work, the empirical dispersion correction schemes proposed by Tkatchenko and Scheffler (TS) [Tkatchenko et al., Phys. Rev. Lett., 2009, 102, 073005] and Grimme (D2) [Grimme, J. Comput. Chem., 2006, 27, 1787] have been used in combination with the PBE semilocal functional for geometry optimizations. We observed that PBE-TS provides a remarkable improvement in predicting the crystal structure of oxalyl dihydrazide polymorphs with respect to commonly used DFT-D functionals. The relative stabilities of the five forms have then been computed at the PBE-TS/D2, PBE0-D2, B3LYP-D2 and B3LYP-D3(BJ)+gCP level on the PBE-TS hydrogen-optimized geometries and benchmarked against high level periodic LMP2 calculations. PBE-TS, B3LYP-D2 and B3LYP-D3(BJ)+E(3) (that is including three-body corrections) achieve good predictions of the stability ordering, though the broadness of the energy range is slightly larger than in the case of LMP2.
Ab initio modelling of molecular crystals is challenging since the crystal packing is governed by intermolecular interactions, such as long range dispersive forces (vdW) or hydrogen bonding. Dispersive interactions are quantum mechanical in nature and correspond to the multipole moments induced in response to instantaneous fluctuations in the electron density. These long range correlation effects are naturally accounted for by correlated wave-function4–6 and Quantum Monte Carlo (QMC)7 methods. In recent years, post-Hartree–Fock methods have become available also for crystalline systems,5,8–11 and in this context molecular crystals have been specially used as a cross-benchmark of different proposals.12–18 On the other hand, Density Functional Theory (DFT)19,20 has become the most widely used method for evaluating energies and electronic properties of molecules and condensed phases since it generally provides accurate results in a reasonable computational time. However, most of the exchange-correlation functionals (including hybrid functionals) do not include the long range electron correlation effects responsible for van der Waals forces, and thus usually fail to predict structures and energetics of molecular crystals. To overcome this problem, several dispersion-correction schemes have been developed from London-type pairwise corrections21–27 based on either empirical parameters (like the well known DFT-D2 scheme21 and its enhanced D3 version)28 or non-empirical parameters (like the exchange-dipole moment model,29 the Tkatchenko–Scheffler (TS) scheme26 or the TS scheme with the inclusion of many-body dispersion (MBD) effects30) and non-local DFT functionals (e.g. vdw-DF31 and VV10).32
Among these, the Tkatchenko–Scheffler model26 is a simple approach that partly recovers a non-empirical character. It has been already successfully applied to investigate various systems by different research groups,33–40 and will be adopted in the present work in the computational study of the oxalyl dihydrazide (ODH) molecular crystals. ODH exploits a peculiar kind of conformational polymorphism with five known polymorphs (i.e. α, β, γ, δ and ε) whose stabilities are determined by the competition between different types of hydrogen bond (HB) and vdW forces between π-stacked molecules. The α phase is experimentally the denser one (1.76 g cm−3) and involves only the intermolecular hydrogen bonds. The other phases (β, γ, δ, ε), instead, present mixtures of inter- and intramolecular hydrogen bonds and are less dense than the α one. According to density measurements, the α phase was proposed to be the most stable even if there are still no certain proofs, due to the incompleteness of the experimental data.41 The authors, in fact, only point out that the stability ordering is α, δ, ε > γ > β (without any specification about the relative stability between α, δ and ε). Such ordering is justified by the endothermic phase transition of α, δ, ε into γ that occurs at different temperatures under 250 °C, and the great difficulty to obtain the receding metastable β phase.
A few theoretical works have been carried out to shed light on the ODH polymorphs. Karamertzanis et al.42 used a GGA dispersion-corrected D2-PW91 functional to predict the relative stability of the ODH, with reasonably good results. The energies span in a range of about 15 kJ mol−1, that is slightly wider than the usually accepted one (10 kJ mol−1) for experimentally observed polymorphs.43 Recently, Wen et al.44 showed that the energy ranking of the five polymorphs is remarkably sensitive to diverse functionals and basis sets,44 and compared B3LYP-D* results with their fragment-based Hybrid Many-Body Interaction method (HMBI)11,16,44,45 – that combines MP2 level short-range one- and two-body terms with MM two- and many-body long-range terms. They obtained structures in good agreement with the experimental ones, but underlined the inefficacy of DFT-D in predicting the energetics, specially if compared to HMBI. In fact, HMBI single-point energies (SPEs), performed on the B3LYP-D* optimized structures, achieved the correct stability ordering of the five ODH polymorphs within an energy span of 5 kJ mol−1.
Here, we then combine the well-known PBE46 functional with the TS dispersion correction and compare it with the PBE, PBE047 and B3LYP48 functionals, augmented with the DFT-D2 approach21 and the B3LYP functional augmented with the D3(BJ) correction.28 We investigate also the combination of B3LYP-D3(BJ) with the geometrical-CounterPoise (gCP) correction.49 The latter theoretical approach has recently been successfully applied to study the cohesive energy of molecular crystals.50 As a further assessment of the TS scheme and to compare with HMBI results, full periodic Local second order Møller–Plesset perturbation theory (LMP2) calculations of ODH polymorphs are also presented for the first time.
In the Tkatchenko–Scheffler correction (TS),26 at difference with the DFT-D2 scheme,21 the electron density is used to calculate the Hirshfeld volumes of an atom in isolation and in the bonding environment. It takes advantage of the relationship between polarizability and volume to compute the effective dispersion coefficient for an atom in a molecule, thus accounting for the relative variation in dispersion coefficients of differently bonded atoms. Computed coefficients are then used in a London-type pairwise correction term. This approach is more rigorous from a physical point of view, and the determination of the dispersion coefficients could be considered nearly nonempirical.
An energy cut-off of 800 eV and a 7 × 4 × 3 k-points mesh (distance between k-points: 0.04 Å−1) have been employed in the calculations. They have been carefully individuated as a good compromise between computational cost and accuracy (see ESI†). In geometry optimizations, tolerances have been set on 5 × 10−6 eV per atom, 5 × 10−4 Å, 1 × 10−2 eV Å−1 and 2 × 10−2 GPa for the total energy convergence, the maximum ionic displacement, the maximum force and the maximum stress, respectively.
The starting experimental structures of ODH were taken from the original work of Ahn et al.41 (available in the Cambridge Crystal Structure Database with reference codes VIPKIO01–VIPKIO05). The polymorphs, all belonging to the P21/c symmetry, contain two molecules (28 atoms) per unit cell in the α, β, γ, δ, ε phases and four molecules (56 atoms) in the γ phase, respectively. Since three-dimensional structures from databases were obtained from X-ray diffraction (XRD) techniques, we carried out geometry optimizations of the hydrogen atoms with the PBE-TS method to obtain refined experimental structures. Hereafter, we will refer to them as revised experimental geometries.
Full relaxations of cell parameters and atomic positions have then also been performed at both PBE-TS and PBE-D2 levels of theory. The fully optimized structures have been compared with the experimental ones in order to assess the accuracy of these dispersion corrected functionals in the prediction of the crystal structure of the five oxalyl dihydrazide polymorphs. The last SCF ground-state energies have also been collected to calculate the relative stability of the five polymorphs at the fully relaxed geometry.
Single-point energy calculations on the revised experimental structures have also been carried out with the PBE0-D2, B3LYP-D2 and B3LYP-D3(BJ) functionals by using a development version of the CRYSTAL code.53,54 The B3LYP-D3(BJ) calculations were performed with and without the gCP correction and including both two- and three-body corrections. All electron basis sets have been employed. The large QZVP55 basis set has been considered for PBE0-D2, B3LYP-D2 and B3LYP-D3(BJ) calculations, while the 6-31G(d) basis set has been used when the gCP correction is included.
Periodic Local-MP2 energies have been computed by means of the CRYSCOR program6,8,56 on the revised experimental structures. A p-aug-6-31G(d,p) Gaussian basis-set has been chosen. This has already proven to be a good compromise between computational cost and accuracy of the results in the study of molecular crystals14 and graphene-based nanostructures.57 The local correlation approach58 allows one to exploit the fundamentally local character of dynamic electron correlation, leading to quantum chemical methods that scale linearly with the size of the system. CRYSCOR implements this kind of approach together with fast integral evaluation techniques59 for correlated calculations (MP2) in periodic systems. It employs the periodic Hartree–Fock solution in the basis of Gaussian-type orbitals (GTO) provided by the CRYSTAL code.53,54
All calculations (DFT and LMP2) have been performed at 0 K temperature. Zero-point energy (ZPE) effects have been evaluated through vibrational frequencies calculations with CASTEP – employing the linear-response phonon DOS (Density Of States)60 approach on the fully optimized PBE-TS structures that are, however, very close to the revised experimental ones. As is shown below, ZPE effects have in this case a negligible impact on relative stabilities.
Table 1 reports the relative deviation % of the optimized cell parameters and volume computed at 0 K, with respect to the experimental ones. B3LYP-D* results from the recent paper of Wen et al.44 and D2-PW91 plane-waves calculation42 are also reported for comparison.
Method/basis-set | a (Å) | b (Å) | c (Å) | γ (°) | Volume (Å3) |
---|---|---|---|---|---|
α polymorph | |||||
Experimenta | 3.622 | 6.832 | 9.129 | 99.30 | 223.0 |
B3LYP-D*/6-31G(d,p)b | −3.66% | −1.24% | −1.44% | −1.81% | −5.80% |
B3LYP-D*/TZPb | −1.33% | −0.26% | −1.73% | −1.58% | −2.90% |
D2-PW91/PWc | −0.80% | −2.28% | −1.14% | −0.61% | −4.04% |
PBE-D2/NCPP | −4.14% | −1.86% | −0.54% | −1.44% | −6.11% |
PBE-TS/NCPP | 0.67% | −1.24% | 0.53% | −0.32% | 0.01% |
ε polymorph | |||||
---|---|---|---|---|---|
a Experiment from ref. 41. b Ref. 44. c Ref. 42. | |||||
Experimenta | 5.364 | 3.841 | 12.319 | 109.00 | 240.0 |
B3LYP-D*/6-31G(d,p)b | 2.52% | −4.91% | 0.23% | 3.19% | −4.52% |
B3LYP-D*/TZPb | 3.86% | −3.56% | 1.36% | 3.76% | −1.25% |
D2-PW91/PWc | −0.39% | 0.57% | −1.44% | 3.01% | −3.43% |
PBE-D2/NCPP | 2.78% | −2.02% | −0.68% | 3.83% | −2.76% |
PBE-TS/NCPP | −0.62% | 1.96% | −0.82% | 1.59% | −0.60% |
The overall accuracy given by PBE-TS is remarkable: the error in predicting cell volumes is notably smaller with respect to the other methods and less than 0.6% for all phases but β whose deviation is around 4%. Moreover, in most cases, PBE-TS functional yields very small deviations on cell parameters.
On the other hand, the Grimme dispersion (D2) scheme tends to overestimate the long range interaction and consequently to give rise to more compact cells.
For the β polymorph – which is the less stable phase, and was experimentally obtained and characterized only once41 – the agreement of computational methods with experiments is rather poor. Although the a parameter is considerably overestimated by the PBE-TS (~7%), b and c compensate with an opposite trend, giving the smaller volume deviation (about 4.5%). Conversely, the PBE-D2 leads to a better agreement with experiment of single cell parameters but a larger discrepancy in the cell volume.
The PBE-TS fully optimized structures of the five polymorphs have been reported in Fig. 1 which also shows the complex intermolecular hydrogen bond network that constitutes the three dimensional structure.
Hydrogen bond | Type | PBE-TS (Å) | Deviation (Å) | PBE-D2 (Å) | Deviation (Å) | Exp.a (Å) |
---|---|---|---|---|---|---|
α polymorph | ||||||
In-plane | N–H⋯N | 2.882 | −0.037 | 2.833 | −0.086 | 2.919 |
N–H⋯O | 2.953 | −0.059 | 2.937 | −0.075 | 3.012 | |
Out-of-plane | N–H⋯O | 3.094 | −0.082 | 3.037 | −0.139 | 3.176 |
γ polymorph | ||||||
---|---|---|---|---|---|---|
In-plane | N–H⋯O | 2.851 | −0.002 | 2.798 | −0.055 | 2.853 |
N–H⋯O | 2.834 | −0.036 | 2.807 | −0.063 | 2.870 | |
N–H⋯N | 3.141 | −0.078 | 3.121 | −0.098 | 3.219 | |
N–H⋯N | 3.105 | −0.043 | 3.092 | −0.063 | 3.184 | |
Out-of-plane | N–H⋯O | 3.045 | −0.074 | 3.019 | −0.100 | 3.119 |
N–H⋯O | 3.029 | −0.079 | 2.986 | −0.122 | 3.108 | |
δ polymorph | ||||||
---|---|---|---|---|---|---|
In-plane | N–H⋯O | 2.821 | −0.027 | 2.801 | −0.047 | 2.848 |
N–H⋯N | 3.086 | −0.076 | 3.060 | −0.102 | 3.162 | |
Out-of-plane | N–H⋯O | 3.072 | −0.075 | 3.046 | −0.101 | 3.147 |
For each ribbon there are two different intermolecular N–H⋯O HBs – one in-plane and one out-of-plane (the latter between parallel ribbons). The N–H⋯N hydrogen bonds from the –NH groups and from the nitrogen of –NH2 groups link adjacent ribbons along c. As always, the PBE-TS gives shorter HB distances that are in agreement with the cell compacting along a and c, thus indicating an overestimation of the interactions between adjacent perpendicular ribbons.
In conclusion, the overall accuracy of PBE-TS in predicting cell parameters and hydrogen bond distances (mean abs. dev. = 0.06 Å and max. abs. dev. = 0.08 Å) of oxalyl dihydrazide appears higher than that obtained with the PBE-D2 functional (mean abs. dev. = 0.11 Å and max abs. dev. = 0.14 Å).
The PBE-TS and PBE-D2 energy rankings, reported in Fig. 2, are consistent with those obtained by D2-PW91 calculations.42 The difference between the most and the less stable phases corresponds to 12–15 kJ mol−1 per molecule, in accordance with the commonly accepted value (10–15 kJ mol−1 per molecule).43,44 PBE-TS leads to the expected experimental stability:41 α, δ, ε > γ > β, which is in agreement with the density ordering of the five polymorphs (α is the most dense – 1.76 g cm−3, β is the less dense – 1.59/1.66 g cm−3) and the experimentally observed endothermic transformations of α, δ, ε into γ while heated below 250 °C – temperature at which ODH decomposes to N,N′-dioxalylhydrazidylhydrazine (DOHH) and gaseous hydrazine through a solid state chemical reaction.
The PBE-D2 method, instead, reverses the stability of the two less stable phases, giving β > γ.
The relative stability of the ODH polymorphs as computed for the revised experimental structure (i.e. optimized hydrogen atoms) and the fully relaxed ones is only slightly different, thus indicating that apart from the H atoms positions the potential energy surface of this system is very flat.
We have finally evaluated the zero-point energy (ZPE) contribution as described in the “computational details” section. Results are shown in Fig. 3. It is possible to observe that the inclusion of the ZPE correction does not change the energy ranking, and the energy range is only slightly broadened (different from that observed by Wen et al.).44
Fig. 3 PBE-TS and ZPE corrected PBE-TS relative stabilities of oxalyl dihydrazide at the fully optimized geometry. |
A summary of the relative stabilities given by the aforementioned approaches is reported in Table 3 and shown in Fig. 4, together with LMP2 results.
Relative stability, ΔE (kJ mol−1 per molecule) | ||||
---|---|---|---|---|
Method | β | γ | δ | ε |
a Basis set: QZVP. b Basis set: 6-31G(d). c Basis set: p-aug-6-31G(d.p). | ||||
PBE-TS | 14.23 | 13.26 | 10.72 | 6.85 |
PBE-D2 | 11.38 | 12.08 | 10.65 | 6.46 |
PBE0a-D2 | 26.21 | 8.86 | 9.35 | 6.21 |
B3LYPa-D2 | 19.45 | 9.85 | 10.36 | 7.93 |
B3LYPa-D3(BJ) | 20.35 | 7.42 | 7.86 | 4.25 |
B3LYPa-D3(BJ)+E(3) | 19.22 | 6.63 | 7.20 | 3.37 |
B3LYPb-D3(BJ)+gCP | 7.30 | 2.18 | 3.95 | −4.07 |
LMP2c | 9.84 | 6.60 | 7.80 | 5.28 |
Fig. 4 Relative stabilities of oxalyl dihydrazide at the PBE-TS revised experimental (hydrogen optimized) geometry. |
The predicted energy ranking obtained by LMP2 calculations is α > ε > γ > δ > β and an energy range below 10 kJ mol−1. These results are slightly different with respect to those reported by Wen et al.44 In fact, their HMBI calculations without including the three body correction provide the δ phase to be more stable than the γ one whereas the inclusion of the three-body dispersion terms reverses the stability order of the α and ε polymorphs. In both cases, the energy range among the five polymorphs is below 5 kJ mol−1. The difference between this value and our LMP2 result can be ascribed to the different approximations inherent to the HBMI and LMP2 approaches, different basis set adopted, and treatment of BSSE effects.
The energy differences obtained at the DFT level of theory underline a marked dependency from the exchange-correlation functional form. In particular, the PBE-TS functional provides a stability ordering more similar to the one obtained by Wen44 since the γ and δ stability is reversed with respect to the LMP2 calculations but the difference between the two phases is very small. Moreover, the energy range among the phases is wider and more similar to that computed at the LMP2 level. On the other hand, PBE0-D2 and B3LYP-D2 provide a stability ordering similar to LMP2 but with a much wider stability range of about 26 and 20 kJ mol−1, respectively.
B3LYP-D3(BJ)/QZVP gives relative stabilities similar to B3LYP-D2 whereas the inclusion of the three-body contribution (E(3)) stabilizes all the other phases with respect to α without changing the stability order. In fact, the E(3) contribution is very small, being about 1.0 kJ mol−1 in agreement with the results reported by Wen et al. in ref. 44.
However, it has to be noted that the Axilrod–Teller–Muto term61,62 is not the most important many-body contribution to dispersion interactions in solids. The polarizability of molecules in solids is significantly influenced by the environment through collective electrodynamic response.63 Very recently, Tkatchenko et al.63–65 have demonstrated that these many-body effects cannot be neglected if one aims at a more accurate prediction of solid-state properties and thus further investigations including such terms should be carried out in the future.
Since Wen et al.44 computed the relative stability of the five polymorphs by employing BSSE corrected total energies we decided to investigate the recently proposed gCP correction,49 which mimics the Boys–Bernardi66 counterpoise correction, in the prediction of the energy ranking of ODH. The gCP correction, as fitted for the 6-31G(d) basis set, is added to the B3LYP-D3(BJ) results. Oddly, the ε phase becomes the most stable one while all other polymorphs are stabilized with respect to α. When compared with other methods, this suggests that the gCP is probably not well balanced.
In spite of the larger stability range, the PBE-TS and B3LYP-D2 functionals seem to be, amongst the considered functionals the most compatible with LMP2 and HMBI results. Interestingly, the B3LYP-D3(BJ)+E(3)/QZVP method gives results very similar to LMP2 for all polymorphs but the β phase.
The predicted energy ranking that can be deduced from the different methods adopted in the present work appears to be: α > ε > γ ~ δ > β. Overall, computed results motivate present and future efforts to assess the applicability of dispersion correction schemes to molecular crystal polymorphism.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c3ce41758a |
This journal is © The Royal Society of Chemistry 2014 |