Effects of Empirical Dispersion Energy on the Geometrical Parameters and Relative Energy of a Salicylideneaniline Molecular Switch in the Solid State

The geometries of the enol (E) and keto (K) forms of a crystalline salicylideneaniline molecular switch, (E)-2-methoxy-6-(pyridine-3-yliminomethyl)phenol (PYV3), have been determined using periodic density functional theory (DFT) calculations with a variety of exchange-correlation functionals (XCFs). They are compared to X-ray diffraction (XRD) data as well as to geometries obtained using empirical dispersion energy in the form of the second iteration of Grimme’s scheme, either with its original parameters (DFT-D2) or with parameters revised for the solid state (DFT-D*). Using DFT, a good agreement with experiment on the unit cell parameters is obtained with the PBEsol, PBEsol0, and ωB97X XCFs. DFT-D2 contracts the unit cell too much with all considered XCFs, whereas DFT-D* lessens this effect thus allowing B3LYP, PBE, and PBE0 to achieve reasonable agreement with respect to XRD data. When considering molecular geometries, both DFT and DFT-D* have a similar effect on the bond lengths, both systematically underestimating (overestimating) the length of the single (double) bonds (within 0.003 A), as well as on valence angles attaining differences of 2° with respect to XRD data. The errors on the torsion angles are less spread out with DFT-D* (averaging 1°) than DFT for which only PBEsol, PBEsol0, and ωB97X perform well. Finally, the relative keto–enol energies, ΔEKE, have been calculated, showing that the inclusion of dispersion energy stabilizes the keto form more than it does the enol form. This results in the PBE and PBEsol XCFs wrongly predicting the keto form as the most stable form.


Introduction
Thermochromic and photochromic compounds have been extensively studied and still attract a lot of attention as they have many potential applications [1][2][3][4][5][6][7][8]. With the focus moving from the liquid phase to the more practical solid one, computational methods well set for quantum calculations of inorganic solids are now being challenged by the molecular crystalline state. In particular, accurate description of the intramolecular parameters (defined by the fractional coordinates of the asymmetric unit) and the intermolecular ones (defined by the unit cell parameters) is required as a starting point for the prediction and study of their properties. Density functional theory (DFT) was recently shown to be an efficient tool granted that the appropriate exchange-correlation functional (XCF) is used [9,10]. In [9], Ruggiero and co-workers assessed the reliability of a range of XCFs for the optimization of three pyridine carboxylic acid crystals and they highlighted the performances of ωB97X [11]. In [10], we showed the effectiveness of three XCFs (HSEsol [12], PBEsol0, and ωB97X [11]) in optimizing the molecular and crystal structures of three salicylideneanilines. Still, with respect to the XCFs used so far, a more precise description of the solid state can be obtained by modifications to the DFT energy by adding London dispersion interactions (Equation (1)) in the form of empirical terms as proposed by Grimme [13].
where E disp is the empirical dispersion energy. We consider the second iteration of this scheme, containing less empirical parameters than the first one, and known as DFT-D2. The general expression for E disp reads with the first summation running over all lattice vectors, g, and the second one running over all atom pairs (excluding the self-interaction case, when i = j for g = 0). s 6 is a scaling factor depending on the functional, C ij 6 = C ii 6 C jj 6 is the dispersion coefficient for the ij pair, R ij,g is the distance between atom i in the reference cell (g = 0) and atom j in cell g, and f dmp R ij,g is a dampening function specified as In the latter expression, R vdw is the sum of the van der Waals radii of atoms i and j and d defines the steepness of the function. Since this scheme was parameterized for clusters, Ugliengo and coworkers proposed to scale the coefficients published by Grimme in order to better describe the dispersion energy in molecular crystals, leading to DFT-D* (these modifications include a decrease of the scaling factor by a factor 0.95 and a scale up of the atomic van der Waals radii by 1.05 and 1.30 for heavy and hydrogen atoms, respectively) [14]. This modification results in smaller dispersion energies.
In this paper, we evaluate the effects of these additional empirical dispersion energy terms (DFT-D2 and DFT-D*) on the optimized geometrical parameters of a molecular crystal from the salicylideneaniline family, (E)-2-methoxy-6-(pyridine-3-yliminomethyl)phenol (PYV3, Figure 1). This compound is in fact a molecular switch that can commute between an enol (E) and a keto (K) form, allowing in parallel to evaluate the effect on the relative energy of these two forms. The theoretical structures obtained with DFT-D2 and DFT-D* XCFs are compared to X-ray diffraction (XRD) data as well as to those results obtained with more traditional DFT functionals [10]. This paper is organized as follows: Section 2 summarizes the key computational aspects (additional details can be found in [10]), whereas Section 3 presents and discusses the results before conclusions are drawn in Section 4. to be an efficient tool granted that the appropriate exchange-correlation functional (XCF) is used [9,10]. In [9], Ruggiero and co-workers assessed the reliability of a range of XCFs for the optimization of three pyridine carboxylic acid crystals and they highlighted the performances of ωB97X [11]. In [10], we showed the effectiveness of three XCFs (HSEsol [12], PBEsol0, and ωB97X [11]) in optimizing the molecular and crystal structures of three salicylideneanilines. Still, with respect to the XCFs used so far, a more precise description of the solid state can be obtained by modifications to the DFT energy by adding London dispersion interactions (Equation (1)) in the form of empirical terms as proposed by Grimme [13].
where Edisp is the empirical dispersion energy. We consider the second iteration of this scheme, containing less empirical parameters than the first one, and known as DFT-D2. The general expression for Edisp reads with the first summation running over all lattice vectors, g, and the second one running over all atom pairs (excluding the self-interaction case, when i = j for g = 0). is a scaling factor depending on the functional, = is the dispersion coefficient for the ij pair, , is the distance between atom i in the reference cell (g = 0) and atom j in cell g, and , is a dampening function specified as In the latter expression, is the sum of the van der Waals radii of atoms i and j and defines the steepness of the function. Since this scheme was parameterized for clusters, Ugliengo and coworkers proposed to scale the coefficients published by Grimme in order to better describe the dispersion energy in molecular crystals, leading to DFT-D* (these modifications include a decrease of the scaling factor by a factor 0.95 and a scale up of the atomic van der Waals radii by 1.05 and 1.30 for heavy and hydrogen atoms, respectively) [14]. This modification results in smaller dispersion energies.
In this paper, we evaluate the effects of these additional empirical dispersion energy terms (DFT-D2 and DFT-D*) on the optimized geometrical parameters of a molecular crystal from the salicylideneaniline family, (E)-2-methoxy-6-(pyridine-3-yliminomethyl)phenol (PYV3, Figure 1). This compound is in fact a molecular switch that can commute between an enol (E) and a keto (K) form, allowing in parallel to evaluate the effect on the relative energy of these two forms. The theoretical structures obtained with DFT-D2 and DFT-D* XCFs are compared to X-ray diffraction (XRD) data as well as to those results obtained with more traditional DFT functionals [10]. This paper is organized as follows: Section 2 summarizes the key computational aspects (additional details can be found in [10]), whereas Section 3 presents and discusses the results before conclusions are drawn in Section 4.

Computational Aspects
Calculations were performed with the Crystal14 package [15] using various XCFs: B3LYP [16], PBE [17], PBE0 [18], PBEsol [19], PBEsol0, and ωB97X [11] with and without including London dispersion interactions. This is achieved by using Grimme's empirical dispersion corrections [13]. First, the original parameters were employed (DFT-D2), followed by those revised for the crystalline state (DFT-D*) [14]. The steepness d of the damping function was set to 20 and a cutoff distance of 25 Å was set for the lattice summations. The default integration grid was used; the truncation criteria for the bielectronic integrals (TOLINTEG keyword) were set to 8 8 8 8 16 for all XCFs, except for ωB97X where 7 7 7 7 16 was used; the maximum order of the shell multipole is kept to the default, as well as all convergence criteria for both the SCF cycles and the optimization. A Monkhorst-Pack shrinking factor of 6 was used, yielding 64 integration points in the irreducible Brillouin zone. Pople's 6-31G(d,p) basis set was used as taken from Basis Set Exchange [20,21].

Crystal Structures and Molecular Geometries
First, we start by comparing the optimized structures of PYV3 obtained using the three 'levels' of DFT, i.e., DFT, DFT-D2, and DFT-D*, and a selection of XCFs to X-ray diffraction (XRD) data obtained at 108 K. According to the electron density map and to key bond lengths, at this temperature molecules are under the enol form [22,23]. The focus is first set on the unit cell volumes (Table 1) and the unit cell parameters ( Figure 2), then on key bond lengths and angles (Figures 3 and 4). Although we compare the calculated pure enol form at 0 K to the XRD structure, the experimental variations between the low and room temperature structures (a~3% increase of volume, up to 0.01 Å differences in bond lengths, and 0.6 • in angles) are well within the accuracy that can be obtained by computations.
Considering the unit cell volume, the use of the reference DFT XCFs yield either large errors with respect to XRD data (B3LYP 17%, PBE 12%, and PBE0 11%) or errors within a satisfying range (≤2%): PBEsol < −1%, PBEsol0 < −1%, ωB97X < −2%. Addition of empirical dispersion energy overall leads to smaller unit cell volumes. In the case of DFT-D2, the unit cells are over-contracted, meaning relative errors ranging from −15 to −8%, while the DFT-D* results tone down this effect to yield volumes either in closer agreement with the experimental values for those XCFs performing poorly with DFT (B3LYP −5%, PBE −3%, and PBE0 −6%), or still over-contracting the unit cell volume for PBEsol, PBEsol0, and ωB97X (errors between −11.6 and −11.2%). Table 1. Differences of unit cell volume [∆V = V(E) − V(XRD), in Å 3 , and ∆V/V = ∆V/V(XRD), in %] as calculated with different XCFs and models for PYV3. The 6-31G(d,p) basis set was used for all calculations.
1 TOLINTEG = 7 7 7 7 16. Figure 2 shows the unit cell parameter (a, b, and c) variations with respect to XRD data. When using conventional DFT XCFs (B3LYP, PBE, and PBE0), variations on the unit cell parameters, all are positive (except for PBE) and larger than 3% (Figure 2a). On the other hand, when using XCFs optimized for the solid state (PBEsol and PBEsol0) as well as ωB97X the differences with respect to XRD data are smaller than 3% and not systematically positive (leading to error cancellation on the estimated unit cell volume, e.g., for PBEsol: the error on a, −2.7%, is partially cancelled by the one on b, 1.9%, while c is fairly accurately estimated, 0.5%).
In the case of DFT-D* ( Figure 2b) the a parameter is strongly underestimated (errors between −6 and −10%) for all XCFs while the errors on b are positive (except for ωB97X) allowing for error This explains why the volume is closer to experiment with these three functionals (−5.6 to −3.3%) than with PBEsol, PBEsol0, and ωB97X (−11.6 to −11.2%). Indeed, for the latter XCFs, the errors on the a and c parameters are both negative (−2.3 to −3.1%) and are not compensated by the errors on b. Finally, the DFT-D2 model, Figure 2c, results in strong and systematic underestimations of all parameters, ranging from −2.5 to −7.7%. A distinction can be made for the a parameter as determined with B3LYP, PBE, and PBE0 for which the error is smaller,~3%, than with the other XCFs,~7%, while the errors for b and c are similar among the XCFs. This results in B3LYP, PBE, and PBE0 volumes that are slightly better than the PBEsol, PBEsol0, and ωB97X ones.
Crystals 2017, 7, x FOR PEER REVIEW 4 of 8 estimated unit cell volume, e.g., for PBEsol: the error on a, −2.7%, is partially cancelled by the one on b, 1.9%, while c is fairly accurately estimated, 0.5%).
In the case of DFT-D* ( Figure 2b) the a parameter is strongly underestimated (errors between −6 and −10%) for all XCFs while the errors on b are positive (except for ωB97X) allowing for error compensation on the volume when c is accurate (between −0.8 and 0.4%): B3LYP, PBE, and PBE0. This explains why the volume is closer to experiment with these three functionals (−5.6 to −3.3%) than with PBEsol, PBEsol0, and ωB97X (−11.6 to −11.2%). Indeed, for the latter XCFs, the errors on the a and c parameters are both negative (−2.3 to −3.1%) and are not compensated by the errors on b. Finally, the DFT-D2 model, Figure 2c, results in strong and systematic underestimations of all parameters, ranging from −2.5 to −7.7%. A distinction can be made for the a parameter as determined with B3LYP, PBE, and PBE0 for which the error is smaller, ~3%, than with the other XCFs, ~7%, while the errors for b and c are similar among the XCFs. This results in B3LYP, PBE, and PBE0 volumes that are slightly better than the PBEsol, PBEsol0, and ωB97X ones. Accurate unit cell parameters constitute a good starting point for further investigation of thermodynamic or optical properties, at least as long as the bond lengths and angles are also accurately modelled. Figure 3 shows key bond length variations with respect to XRD data s along the O-C2-C1-C7-N-C8 path (Figure 1). Moreover, two bond length ratios are considered, the N-and Cratios, which describe the π-conjugation and are defined as Accurate unit cell parameters constitute a good starting point for further investigation of thermodynamic or optical properties, at least as long as the bond lengths and angles are also accurately modelled. Figure 3 shows key bond length variations with respect to XRD data s along the O-C 2 -C 1 -C 7 -N-C 8 path (Figure 1). Moreover, two bond length ratios are considered, the N-and C-ratios, which describe the π-conjugation and are defined as The comparison to experiment is not presented for DFT-D2 since we just showed its poor reliability for unit cell parameters. Figure 3a, for DFT, and Figure 3b, for DFT-D*, highlight the systematic underestimation (overestimation) of the single (double) bond lengths, with the exception of C 1 -C 7 and C 7 =N with ωB97X. Still, in the worst case scenarios, these errors are smaller than 0.025 Å and generally of the order of 0.01 Å. This results in N-/C-ratios that are slightly underestimated, especially in the case of the two GGA XCFs, PBE and PBEsol, whereas the ωB97X's ratios are slightly overestimated. The DFT-D* results are extremely similar to the DFT ones since the differences do not exceed 0.003 Å. This affects the ratios by about the same amount. The comparison to experiment is not presented for DFT-D2 since we just showed its poor reliability for unit cell parameters. Figure 3a, for DFT, and Figure 3b, for DFT-D*, highlight the systematic underestimation (overestimation) of the single (double) bond lengths, with the exception of C1-C7 and C7=N with ωB97X. Still, in the worst case scenarios, these errors are smaller than 0.025 Å and generally of the order of 0.01 Å. This results in N-/C-ratios that are slightly underestimated, especially in the case of the two GGA XCFs, PBE and PBEsol, whereas the ωB97X's ratios are slightly overestimated. The DFT-D* results are extremely similar to the DFT ones since the differences do not exceed 0.003 Å. This affects the ratios by about the same amount.  As for the valence and torsion angles, variations with respect to XRD data are presented in Figure 4 for two valence angles (C1-C7-N and C7-N-C8), the three torsion angles involving the C7=N bond (C2-C1-C7-N, C1-C7-N-C8, and C7-N-C8-C9), one torsion angle associated with the H-bond (C2-O-N-C7), and one torsion angle describing the global torsion of the molecule (C2-C1-C8-C9). As before, only DFT and DFT-D* results are discussed (Figure 4a,b, respectively). The absolute errors on the valence angles are within a [−2°; 2°] range for all XCFs using both types of DFT (~2% of relative error). For the torsion angles, in the case of DFT without dispersion, they are systematically underestimated, except for ωB97X. Furthermore, the error on the molecular torsion angle C2-C1-C8-C9 and, to a smaller extent on the angle C7-N-C8-C9, is quite large for those XCFs that poorly perform for the unit cell parameters (with errors ranging between −5.7° and −10.5° for B3LYP, PBE, and PBE0). For the other XCFs (PBEsol, PBEsol0, and ωB97X), the deviations are much smaller: −3.1°, −3.1°, and −1.0°, respectively. This highlights the key role of the torsion angles on the unit cell parameters or, in other words, their interdependence. It is however difficult to say whether the error on the torsion angles drives the error on the unit cell or vice versa.
When using the DFT-D* model, the amplitude of the errors is greatly reduced compared to standard DFT, in particular for B3LYP, PBE, and PBE0. Note that the scale of Figure 4a for DFT goes from −12 to 4° while that of Figure 4b for DFT-D*, from −2 to 4°. For DFT-D*, the errors on the torsion angles are mainly positive, except for the C1-C7-N-C8 angle. For the important molecular torsion (C2-C1-C8-C9), the error is around 1° for all XCFs, except for PBE0 (<1°) and ωB97X (2°), which is on par with the errors of PBEsol, PBEsol0, and ωB97X with DFT for the same angle. Accurate molecular torsions are obtained in parallel to accurate unit cell parameters. As for the valence and torsion angles, variations with respect to XRD data are presented in Figure 4 for two valence angles (C 1 -C 7 -N and C 7 -N-C 8 ), the three torsion angles involving the C 7 =N bond (C 2 -C 1 -C 7 -N, C 1 -C 7 -N-C 8 , and C 7 -N-C 8 -C 9 ), one torsion angle associated with the H-bond (C 2 -O-N-C 7 ), and one torsion angle describing the global torsion of the molecule (C 2 -C 1 -C 8 -C 9 ). As before, only DFT and DFT-D* results are discussed (Figure 4a,b, respectively). The absolute errors on the valence angles are within a [−2 • ; 2 • ] range for all XCFs using both types of DFT (~2% of relative error). For the torsion angles, in the case of DFT without dispersion, they are systematically underestimated, except for ωB97X. Furthermore, the error on the molecular torsion angle C 2 -C 1 -C 8 -C 9 and, to a smaller extent on the angle C 7 -N-C 8 -C 9 , is quite large for those XCFs that poorly perform for the unit cell parameters (with errors ranging between −5.7 • and −10.5 • for B3LYP, PBE, and PBE0). For the other XCFs (PBEsol, PBEsol0, and ωB97X), the deviations are much smaller: −3.1 • , −3.1 • , and −1.0 • , respectively. This highlights the key role of the torsion angles on the unit cell parameters or, in other words, their interdependence. It is however difficult to say whether the error on the torsion angles drives the error on the unit cell or vice versa.
When using the DFT-D* model, the amplitude of the errors is greatly reduced compared to standard DFT, in particular for B3LYP, PBE, and PBE0. Note that the scale of Figure 4a for DFT goes from −12 to 4 • while that of Figure 4b for DFT-D*, from −2 to 4 • . For DFT-D*, the errors on the torsion angles are mainly positive, except for the C 1 -C 7 -N-C 8 angle. For the important molecular torsion (C 2 -C 1 -C 8 -C 9 ), the error is around 1 • for all XCFs, except for PBE0 (<1 • ) and ωB97X (2 • ), which is on par with the errors of PBEsol, PBEsol0, and ωB97X with DFT for the same angle. Accurate molecular torsions are obtained in parallel to accurate unit cell parameters.

Keto-Enol Energies
After investigating the geometrical structures, we now turn to studying the effect of the methods used on the relative energy of the keto and enol forms (Table 2). Experimentally, XRD measurements show only the enol form of PYV3 to be present both at low and room temperatures. This highlights the fact that the enol form (E) is more stable than the keto one (K). The relative energy difference, ΔEKE = E(K) − E(E), is thus expected to be positive. Table 2 shows the relative energies computed with the selected XCFs. Although ΔEKE varies significantly depending on the XCF, all DFT values are positive, ranging from 1 to 13 kJ/(mol of asymmetric unit), in agreement with experiment. Addition of dispersion energy has a more important stabilizing effect on the keto form compared to the enol one, leading overall to less positive ΔEKE values, if not negative ones like in the case of PBE and PBEsol [−2.5 and −2.7 kJ/(mol asym. unit) for DFT-D2 and −1.3 and −1.9 kJ/(mol asym. unit) for DFT-D*, respectively]. Since we know from experiment that ΔEKE is positive, these two XCFs are not considered any further. DFT-D2 values, ranging from 2.7 to 10.2 kJ/(mol asym. unit), vary little from the DFT-D* ones (ranging from 3.4 to 9.9 kJ/(mol asym. unit)). Still, as expected, the dispersion energies obtained with DFT-D2 are larger than with DFT-D*, which favors the keto form hence decreasing the relative energy difference ΔEKE (with the exception of ωB97X). The energy difference can be used to analyze the K/E Boltzmann distributions at 298.15 K to further quantify the effects of the dispersion energy. In the case of ωB97X, barely no change is observed, with the keto population increasing from 1% with DFT to 2% with DFT-D2/D*. The effects are much larger for the other functionals, and they correspond to at least a 10% increase of the keto form upon adding dispersion, i.e., from 1 to 14% for PBE0, from 3 to 25% for B3LYP, and from 5 to 20% for PBEsol0.

Keto-Enol Energies
After investigating the geometrical structures, we now turn to studying the effect of the methods used on the relative energy of the keto and enol forms (Table 2). Experimentally, XRD measurements show only the enol form of PYV3 to be present both at low and room temperatures. This highlights the fact that the enol form (E) is more stable than the keto one (K). The relative energy difference, , is thus expected to be positive. Table 2 shows the relative energies computed with the selected XCFs. Although ∆E KE varies significantly depending on the XCF, all DFT values are positive, ranging from 1 to 13 kJ/(mol of asymmetric unit), in agreement with experiment. Addition of dispersion energy has a more important stabilizing effect on the keto form compared to the enol one, leading overall to less positive ∆E KE values, if not negative ones like in the case of PBE and PBEsol [−2.5 and −2.7 kJ/(mol asym. unit) for DFT-D2 and −1.3 and −1.9 kJ/(mol asym. unit) for DFT-D*, respectively]. Since we know from experiment that ∆E KE is positive, these two XCFs are not considered any further. DFT-D2 values, ranging from 2.7 to 10.2 kJ/(mol asym. unit), vary little from the DFT-D* ones (ranging from 3.4 to 9.9 kJ/(mol asym. unit)). Still, as expected, the dispersion energies obtained with DFT-D2 are larger than with DFT-D*, which favors the keto form hence decreasing the relative energy difference ∆E KE (with the exception of ωB97X). The energy difference can be used to analyze the K/E Boltzmann distributions at 298.15 K to further quantify the effects of the dispersion energy. In the case of ωB97X, barely no change is observed, with the keto population increasing from 1% with DFT to 2% with DFT-D2/D*. The effects are much larger for the other functionals, and they correspond to at least a 10% increase of the keto form upon adding dispersion, i.e., from 1 to 14% for PBE0, from 3 to 25% for B3LYP, and from 5 to 20% for PBEsol0.

Conclusions
Density functional theory has been challenged for the geometry optimization of molecular switches in their solid crystalline state. By performing comparisons with X-ray diffraction (XRD) data, a recent contribution [10] has demonstrated that HSEsol [12], PBEsol0, and ωB97X [11] can already be effective but, in this work, we investigated whether the addition of empirical dispersion energy, as proposed by Grimme [13], could further improve these results. First, we have shown that the use of the original dispersion parameters (DFT-D2) over-contracts the unit cell for the selected XCFs (B3LYP [16], PBE [17], PBE0 [18], PBEsol [19], PBEsol0, and ωB97X). On the other hand, the down-scaled parameters proposed for the solid state (DFT-D*) [14] decrease this effect so that the B3LYP, PBE, and PBE0 XCFs achieve a rather good agreement with XRD data when considering the unit cell volume, though, mostly due to error compensations. Looking at the molecular geometries, the main conclusions are: (i) inclusion of dispersion energy has almost no effect on the bond lengths, though systematically underestimating (overestimating) the length of the single (double) bonds, with the maximum difference between DFT and DFT-D* attaining 0.003 Å; (ii) the valence angles are also barely affected when using DFT-D* compared to DFT with relative errors with respect to XRD data of 2% or less in both cases; and iii) in the case of the torsion angles, the use of DFT-D* XCFs improves the results since the variations with respect to XRD data are less spread out. The average errors with DFT-D* are of the order of 1 • whereas with DFT, only the PBEsol, PBEsol0, and ωB97X XCFs perform well. This means that an accurate description of the unit cell parameters leads to accurate molecular torsions but that accurate molecular torsions do not constitute a sufficient condition to fully describe the intermolecular interactions, and to reproduce the XRD unit cell. Finally, for all XCFs, the relative keto-enol energy differences ∆E KE have been calculated, showing that the inclusion of dispersion stabilizes the keto form more than the enol form. As a consequence, the PBE and PBEsol XCFs incorrectly predict the keto form to be the most stable. The other functionals, with exception of ωB97X, predict a decrease of ∆E KE but the overall value remains positive. Overall, these results show PBEsol0 and ωB97X XCFs to be reliable in predicting molecular crystal structures and that there is no clear advantage of using empirical energy dispersion corrections as originally proposed [13] or later reparameterized [14] for solids.