Lattice energies and thermal correction

Within this section the determination of the sublimation free energy from periodic DFT in combination with the mixed “exp/theory” method is discussed. The resulting ΔGsubl are subsequently used for predicting solubilities via the sublimation cycle. As lattice energies are calculated by DFT for the static crystal lattice and the static ideal gas, zero-point energy as well as thermal enthalpy and entropy corrections are required and will be discussed first. Assuming that the vibrational lattice and molecular modes are decoupled the zero-point energy difference between the crystal and the gas, ∆E_ZPE, can be derived from the lattice frequencies as molecular vibrations cancel out. Experimental vibrational lattice frequencies, νi, from literature^{147, 181} were used to calculate the zero-point energy contribution by ∆E_ZPE=E_ZPE^g −E_ZPE^s =-h⁄ ∑ ν2 _i. The obtained ∆E_ZPE (see Table 4.3) is in excellent agreement with published calculated values for naphthalene¹⁴⁷ and benzoic acid¹⁸².

Figure 4.1: Temperature-dependent solid state (experiment¹⁸³ - ○; calculations via the

“exp/theory”-method - straight line) and calculated ideal gas (dashed line) (a) heat capacities CP

and (b) CP/T as well as the corresponding thermodynamic (c) enthalpy and (d) entropy functions of naphthalene.

Thermal enthalpies, entropies and Gibbs energies are related to the integrals of the heat capacities in the specific state of matter. Integrals of CP from zero Kelvin give the enthalpy while thermal entropies are related to integrals of CP/T. Within the mixed

“exp/theory” method the experimental solid state heat capacities from literature^{168, 183} were interpolated using while ideal gas heat capacities and thermodynamic functions have been calculated via the RRHO (Table 2.2) as described in section 2.7.1 using the def2-TZVP-BP86-D3 level of theory for calculating the harmonic vibrational frequencies. This is exemplarily shown in Figure 4.1 for the crystalline solid and ideal gas of naphthalene in the whole temperature range from 0 K to 298 K.

Integrals of the experimental solid state^{168, 174} and calculated ideal gas heat capacities between zero Kelvin and 298 K were combined with the ∆E_ZPE to give the final thermal enthalpy correction, Hcorr, and the sublimation entropy, TΔSsubl, via eq. (44) which are summarized in Table 4.3 for both molecules. There are only minor differences of the calculated Hcorr to the 4.9 kJ·mol^-1 of the frequently used 2RT-approximation (see eq.

(55)). The calculated sublimation entropies via the mixed “exp/theory” method are in good agreement with the primary experimental data from measured vapor pressures (Table 4.1) with absolute deviations of 3.3 and 0.4 kJ·mol^-1 for naphthalene and benzoic acid, respectively. However, when used to calculate solubilities these relatively small errors can affect the predictions accuracy due to the exponential relation within eq. (10).

Table 4.3: Summary of all zero-point, thermal enthalpy and entropy corrections which are required to calculate the sublimation free energy, ΔGsubl, via eq. (44).

Naphthalene Benzoic acid

1 ∆EZPE 2.30^a 2.75^a

2 H^ig ± stdev 20.67 ± 0.30^b 21.41 ± 0.29^b

3 ∫ CP,exp dT 24.79^c 24.03^c

4 Hcorr -6.42 ± 0.30^d -5.37 ± 0.29^d 5 TS^ig ± stdev 102.91 ± 0.57^b 106.09 ± 0.51^b

6 ∫ (CP,exp/T) dT 49.44^c 49.34^c

7 TΔSsubl 53.47 ± 0.27^e 56.74 ± 0.23^e

a Calculated from experimental lattice vibrational frequencies^{147, 181} via ∆E_ZPE=−h⁄ ∑ ν2 _i. b Average of BP86-, B3LYP- and MO6-def2-TZVP method together with the standard deviation.c Calculated from integrals of the solid-state heat capacities^{168, 174}. d Sum of 1, 2 and 3 according to eq. (44). e Sum of 5 and 6 according to eq. (44).

Lattice energies have been determined via eq. (42). Geometry optimizations of the ideal gas molecules as well as of crystalline solid have been performed within the TURBOMOLE software package (V7.1)¹³⁹. Further details are given in section 2.7.2.

Calculations have been performed using three GGA type density functionals, B-P86^134,

135, BLYP^{135, 184} and PBE¹³⁶ as well as the semi-empirical GGA-type B97-D¹⁸⁵. They are

combined with various Gaussian def2-type basis sets by Weigand and Ahlrichs¹³³ and the D3 dispersion correction by Grimme et al.¹³⁷. A k-point sampling was performed from k

= 1 to k = 21 using uniform k-points in all three dimensions and no significant change for k > 3 was observed in test runs for benzoic acid and naphthalene (see Table 7.22 in the appendix). Hence, lattice energies have been calculated using a uniform k-points mesh of k = 3x3x3 for all calculations.

Lattice energies have been calculated for the experimental crystal structures. They contain the necessary information of the dimensions of the unit cell as well as starting atomic coordinates for the periodic DFT calculations. Two representatives of crystal structures of benzoic acid and naphthalene, NAPHTA04¹⁸⁶ and BENZAC02¹⁸⁷, have been used (see Figure 4.2). Details on the unit cell dimensions are given in Appendix Table 7.16.

Figure 4.2: Crystal structures of (left) naphthalene – NAPHTA04¹⁸⁶ and (right) benzoic acid – BENZAC02¹⁸⁷.

The molecular crystals differ in the type of molecular interactions that dominate. In case of benzoic acid cyclic hydrogen-bonds between the carboxylic acid groups form strong dimer structures.⁹⁶ In contrary, naphthalene forms a layered structure which is dominated by van der Waals interactions.¹⁸⁸

Sublimation free energies can be determined via eq. (43) from calculated lattice energies when combined with thermal and entropy corrections from Table 4.3. Figure 4.3 (a) and (b) show the absolute deviations of the calculations from the experiment summarized in Table 4.1. All calculations overestimate the sublimation Gibbs energies in comparison to experiment. For both substances (NAPHTA04¹⁸⁶ and BENZAC02¹⁸⁷) there is a drastic increase in accuracy when using the split valence triple-zeta (TZV) basis set (def2-TZVP) over the smaller split valence double-zeta (SV) one (def2-SVP). An inclusion of a second polarization term for all hydrogen atoms does not change the overall quality of the calculation. This trend is consistent for all three functionals and slightly larger for benzoic acid. It is related to the strong anisotropic electrostatic interactions of the hydrogen bonds between the valence electrons of the carboxyl groups which are more likely to be affected by basis set superposition errors (BSSE).¹⁸⁹ This issue will be discussed in more detail later on in section 4.2 for a larger set of molecular crystals.

For both molecules, the PBE and the semi-empirical B97-D functional perform best whereas the BLYP functional leads to the largest deviation from experiment. In case of the weakly bound naphthalene, the differences between the functionals are significantly larger than for the hydrogen-bonded benzoic acid. Even though for naphthalene the PBE functional performs best, the semi-empirical B97-D functional gives the best overall performance whereas the BLYP functional leads to the largest deviations to experiment.

For def2-TZVP/B97-D3 deviations are 2.9 kJ mol^-1 for naphthalene and 1.6 kJ mol^-1 for benzoic acid. Based on the above, the B97-D method in combination with a def2-TZVP basis set will be used in the following for solubility predictions.

Figure 4.3: Differences in the calculated sublimation Gibbs energies and experiment from Table 4.1 of (a) naphthalene (NAPHTA04) and (b) benzoic acid (BENZAC02). ΔGsubl,calc combines calculated lattice energies with thermal corrections from Table 4.3. For Elatt BLYP (■), B-P86 (▲), PBE (◆) and B97 (○) were used in combination with the “D3” dispersion correction def2-type basis sets with increasing size from left to right.

As the unit cell geometry is not optimized, the quality of the result is depending on the experimental unit cell parameters. Lattice energies of two structurally related experimental unit cell geometries for naphthalene (NAPHTA04¹⁸⁶ and NAPHTA23¹⁹⁰) as well as benzoic acid (BENZAC01¹⁹¹ and BENZAC02¹⁸⁷) give energy difference that are smaller than 1.5 kJ mol^-1. Further details are given in Appendix Figure 7.2. These rather minor differences can be minimized when optimizing the unit cell parameters in addition to the molecular geometries. A combined optimization of the unit cell and molecular geometries, however, significantly increases computational times.