Volatility, thermodynamic properties and dispersion interactions of sulfur-containing tricyclic molecular materials

Abstract

Organic heterocyclic molecules play the role of precursors for various sophisticated materials from semiconductors to pharmaceuticals. Knowledge of their volatility is always an important factor to minimize hazards associated with their use and potential toxicity. Nevertheless, thermodynamic properties of polycyclic heterocycles have been very scarcely studied experimentally. In silico approaches can help provide the required data in many cases, but established benchmarks assessing the performance of first-principles models of the sublimation equilibrium for molecular crystals do not cover sulfur-based heterocyclic materials. This work aims at filling these obvious knowledge gaps at both experimental and computational sides. In this work, reference experimental sublimation data are established for four nitrogen- or sulfur-based heterocyclic compounds containing a structural motif of three fused rigid rings. Vapor pressure measurements and calorimetric experiments across broad temperature ranges are carried out to provide reliable reference data for stringent benchmarking of first-principles models of the crystal cohesion. Since the selected materials possess a very limited to no potential for hydrogen bonding, other non-covalent interactions such as dispersion or π–π stacking govern their cohesion. An accurate description of the dispersion interactions in heterocyclic polyaromatic molecules is challenging within density functional theory (DFT). Accuracy of popular post hoc dispersion corrections to DFT is benchmarked for the target heterocyclic materials, revealing that adopting the latest D4 dispersion model along with a lower-tier DFT functional does not necessarily lead to improvements of the computational accuracy over older dispersion models for this class of materials.

---

Introduction

Applications of organic molecular crystals containing aromatic heterocyclic cores, in particular those containing nitrogen or sulfur, can be found in fields ranging from electronics to pharmaceuticals.^1,2 When compared to ionic or covalent solids, design and customization of organic molecular crystals can fully exploit the vast possibilities of contemporary synthetic chemistry.³ Yet, their flexibility also comes at a cost. Insufficient understanding of the non-covalent interactions in molecular crystals that govern their cohesion, and optionally polymorphism, can induce unexpected behavior for an otherwise well-characterized material. Costly consequences may then follow, including patent litigation and drug recalls in the pharmaceutical sector,^1,3 or insufficient performance of optoelectronic devices based on organic-semiconductor (OSC) materials.^4,5 Polymorphism is influenced by an interplay of non-covalent interactions in the bulk material, and thus, accurate models of these interactions are crucial in addressing issues linked to crystal cohesion, structural stability, and polymorphism of molecular crystals.⁶

Non-covalent interactions encompass a broad variety of effects, including dispersion, hydrogen bonding, Coulomb interactions of permanent molecular multipoles, and induction interactions. Most of these interactions are commonly intensified by the presence of delocalized π-electron systems and heteroatoms with lone electron pairs or low-lying vacant orbitals.^2,7,8 Notably, the presence of rigid planar moieties causes dimers of aromatic molecules to be arranged in sandwich, T-shaped, or parallel-displaced configurations,⁹ giving birth mainly to dispersion interactions and improper hydrogen bonding.¹⁰ For instance, π–π and X–π (X = H, C, N, S) interactions have been reported to contribute substantially to binding energies in small benzene clusters,¹¹ or in crystals of monosubstituted benzenes¹² and heterocyclic molecules.^13,14 In contrast, the effect of induction on pair interaction energies or aromatic systems is rather small due to the induction-exchange coupling unless the dimers organize into a T-shaped configuration where induction becomes more pronounced,^10,15 or unless suitable substituents that polarize the delocalized electron clouds are present.¹⁶

Although dispersion-bound molecular clusters exhibit rather small interaction energies, dispersion interactions play a crucial role in stabilization of bulk molecular materials.¹⁷ Due to the ubiquity of dispersion interactions, it is necessary to accurately capture such effects when properties of crystals are modelled in silico. Techniques for describing the dispersion interactions range from DFT with semi-classical dispersion corrections^18–21 to more elaborate and costly ab initio wavefunction methods. Popular representatives of the latter include the second-order perturbative MP2 theory (namely its dispersion-corrected MP2C^22,23 and MP2D²⁴ flavors), and the coupled-clusters theory, including the CCSD(T) method representing the gold-standard²⁵ for modelling non-covalent interactions, or its more affordable approximations such as the domain-based pair-natural-orbital (DLPNO) formulation,²⁶ as well as symmetry-adapted perturbation theory (SAPT),^27,28 or diffusion Monte Carlo simulations.^6,9

While some ab initio methods have been shown to minimize the computational uncertainties related to interaction energies down to the sub-kJ mol⁻¹ scale, the cost of most of the high-performing methods scales too steeply with the system size to be applied for periodic systems.^29,30 On the other hand, periodic DFT methods are favored for their good performance-to-cost ratio, and as such, dispersion-corrected DFT remains a popular option for describing dispersion interactions in organic molecular crystals.^31–39

Within DFT, dispersion interactions can be treated self-consistently within the density functional itself, or by using post hoc dispersion corrections. Examples of the former approach include non-local functionals, such as vdW-DF^40,41 or LC-VV10,⁴² which can be tuned to high accuracy for particular systems (documented by the 0.6 kJ mol⁻¹ errors in vdW-DF2 binding energies in the S66 benchmark);⁴³ however, they may also suffer from a limited transferability, especially in the diverse realm of organic molecular crystals.⁴⁴

Post hoc corrections include the XDM⁴⁵ and MBD models⁴⁶ and numerous methods from the DFT-D family^18–20 that belong to the most popular dispersion-correction strategies. A widespread DFT-D3 model²⁰ takes into account coordination numbers of individual atoms to predict their dispersion interactions. Adopting suitable damping schemes, such as the Becke–Johnson model D3(BJ), proved to be important for correcting the otherwise non-physical behavior of too close molecular contacts.⁴⁷ Its ultimate successor, DFT-D4,¹⁸ builds upon these improvements by introducing atomic charge-dependent functions and describing three-body dispersion effects using the Axilrod–Teller–Muto term.^48,49 DFT-D4 was reported to have consistently outperformed its predecessor on larger organic, metal-containing, or periodic systems.^18,50

Notably, the favorable accuracy of various DFT-D approaches was confirmed in several benchmarks for first-principles modelling of cohesion of molecular crystals, such as X23,^37,39 Z20,^51,52 G60⁵³ and others.⁵⁴ However, sulfur-based heterocyclic molecules that feature widespread conjugation of π-electron systems have not been covered at all by any of those benchmarks, although such structural motifs are particularly suitable for designing OSC materials.^4,7,8 The dominance of the dispersion interactions and the presence of vast molecular π-electron systems in such materials impart substantial delocalization errors, especially for lower-tier DFT functionals.⁵⁵ These factors contribute to potential challenges associated with modeling crystals of large OSC materials with periodic DFT-D approaches.

This paper compares the performance of popular DFT-D3 and DFT-D4 models in predicting the thermodynamic properties of crystals of carbazole, dibenzothiophene, phenothiazine, and thianthrene molecules, which are depicted in Fig. 1. Information about their computationally addressed crystal structures is provided in Table 1. This molecular set has been selected based on the presence of similarly sized and shaped sulfur- or nitrogen-containing heterocycles in their molecules, on the availability of experimental crystal structures in the Cambridge Structure Database,⁵⁶ and on presumed sufficiently high vapor pressures enabling a direct experimental detection. On a more practical note, the materials targeted in the current study act as precursors for a wide range of optoelectronic applications, owing to their structures rich in π electrons that facilitate charge transport and shape persistence.^57–59 These molecules have been used to design low-band-gap materials,⁵⁷ high-performance redox cathodes,⁶⁰ solar cells,⁶¹ as well as lasers.⁶²

Fig. 1 Molecular structures of the target materials selected for this study. Top row – carbazole (CAZ) and dibenzothiophene (DBT); bottom row – phenothiazine (PTZ) and thianthrene (TTH).

Table 1 Summary of the target materials and parameters of their crystal structures considered in this work

Molecule	Formula	Acronym	CSD refcode	Space group	Z	Z′	Molecular size^a
a Largest interatomic distance within a molecule in its experimentally determined crystal phase structure.
Carbazole	C12H9N	CAZ	CRBZOL11^73	Pnma	4	0.5	8.63 Å
Dibenzothiophene	C12H8S	DBT	DBZTHP01^74	P21/n	4	1	8.87 Å
Phenothiazine	C12H9NS	PTZ	PHESAZ01^75	Pnma	4	0.5	9.22 Å
Thianthrene	C12H8S2	TTH	THIANT05^76	P21/c	4	1	8.80 Å

---

Similarities in molecular geometries of the target materials facilitate a straightforward comparison of performance of DFT-D3 and DFT-D4 models. Quasi-harmonic approximation (QHA)⁶³ is used in this work alongside the two DFT-D models to predict structural and thermodynamic properties of their crystals at variable finite temperatures. While the validity and reliability of QHA relying on DFT-D methods have already been demonstrated in benchmarks^32,37,64 and studies of compact^35,65,66 and flexible^67–70 molecules, their performance in our target OSC-precursor molecules is yet to be evaluated. In response to the sulfur-containing heterocyclic molecules being under-represented in existing benchmarks for molecular crystals, we address the question whether the established quasi-harmonic DFT-D protocols can be expected to retain their performance also for these materials.

Accuracy of first-principles calculations cannot be in general benchmarked without low-uncertainty reference data. Experimental data on heat capacities of the crystals and even on sublimation pressures of some of these materials are available in the literature, as discussed below. Still, novel measurements of the vapor pressures and of the phase behavior with an unified state-of-the-art experimental methodology, based on a simultaneous correlation of various thermodynamic properties,^71,72 were due to develop such low-uncertainty reference data in a consistent manner.

Methodology

DFT models of the crystals

Experimental crystal structures from the Cambridge Structural Database (CSD)⁵⁶ were used to define the initial atomic coordinates and unit-cell geometries of the target crystals. The geometric parameters of each crystal were then relaxed under a constrained space-group symmetry to determine its optimal volume V₀ based solely on its electronic energy. DFT calculations were performed using plane-wave basis sets using the projector-augmented wave (PAW) method,⁷⁷ imposing a kinetic energy cut-off 1000 eV.

The PBE functional⁷⁸ was used alongside the hard PAW potentials.⁷⁹ All the underlying PBE-TS⁸⁰ and PBE-D3(BJ)^20,47 calculations were performed in VASP 5.4.1, while PBE-D4^18,50 calculations were performed in VASP 6.4.2.^81–83 In addition, the more sophisticate r²SCAN functional,⁸⁴ as implemented in VASP 6.4.2,^81–83 was used for comparison along with the same dispersion-correction models. In particular, dispersion parameters were taken from ref. 85 for both r²SCAN-D3(BJ) and r²SCAN-D4 models, whereas our current r²SCAN-TS model relied on dispersion parameters originally developed for the SCAN functional.⁸⁰ The acronyms PBE-D3 and r²SCAN-D3 are hereafter used to refer to the PBE-D3(BJ) and r²SCAN-D3(BJ) levels of theory, respectively. The reciprocal k-space was sampled with a Monkhorst–Pack mesh⁸⁶ centered at the Γ-point such that approximately 30/a_i k-points were sampled in each direction, where a_i stand for the lengths of individual direct unit-cell vectors.⁸⁷

Quasi-harmonic approximation was used to mimic anharmonic effects on the total Helmholtz energy A_tot of each crystal structure.^32,63 Here, A_tot was formulated as the sum of electronic (E_el) and vibrational (A_vib) contributions, varying with temperature T and molar volume V_m:

A_tot(T, V_m) = E_el(V_m) + A_vib(T, V_m). (1)

Unit-cell energies at different volumes are necessary to get an expression for E_el(V_m). To sample these trends, the unit cell of each crystal structure was optimized at 14 constrained molar volumes ranging from 0.95 V₀ to 1.08 V₀. The resulting energies were then fitted to the Murnaghan equation of state,⁸⁸ separately for the compression and expansion branches to ensure a higher flexibility of the model.²⁹ In parallel, phonon calculations were performed for each optimized geometry to derive A_vib(T, V). The finite-displacement method⁸⁹ and the Phonopy code, version 2.20,⁹⁰ were used to calculate the density of phonon states and Γ-point vibrational frequencies within the QHA framework. Force constants of the normal phonon modes were evaluated from forces acting on atoms in perturbed supercells that exceeded 10 Å in size, where symmetry-unique atoms were displaced from their equilibrium positions by 0.01 Å.³² A_vib was modelled as a quadratic function of volume having temperature-dependent coefficients,⁸⁷ in effect representing A_vib as a function of both temperature and volume. Note that the full QHA workflow was executed only using the PBE-D3 and PBE-D4 theories due to the considerable computational costs of QHA performed with DFT functionals beyond the GGA level, although the respective QHA model fully performed at the r²SCAN-D levels of theory appears as promising in terms of its overall accuracy.⁹¹

Ab initio analysis of crystal cohesion

Fragment-based many-body expansion model^29,92 was used to refine the interaction energies of the closest two-molecule (dimers) and three-molecule (trimers) clusters extracted from crystal structures optimized at the PBE-D3 level of theory. A very short cut-off distance at 3.5 Å was imposed for the dimers as we focused only on the most important cohesive features within the first coordination shell, not aiming to thoroughly refine the overall cohesive energy of the crystals. For particular interactions found within that cut-off, we benchmarked the performance of both PBE-D3 and PBE-D4 models against a reference level of theory, namely the DLPNO formulation of the coupled clusters with iterative single and double excitations and a perturbative correction to triple excitations, CCSD(T).²⁶ Extrapolations towards the complete basis set⁹³ were performed from energy calculations with cc-pVTZ and cc-pVQZ basis sets for all dimers and with cc-pVDZ and cc-pVTZ basis sets for all trimers, as implemented in ORCA, version 5.0.3.^94,95 Within the DLPNO formalism, the VeryTightSCF and TightPNO convergence settings were applied within the self-consistent field iterative scheme and orbital localization, respectively. The AutoAux procedure was called in ORCA to set up all necessary auxiliary basis sets.⁹⁶

To localize important molecular contacts for the most intense attractive interactions, we analyzed non-covalent interactions within the framework of quantum theory or atoms in molecules.⁹⁷ For this purpose, reduced density gradients (RDG) within the closest dimers were analyzed with respect to the sign of the second eigenvalue of the electron-density Hessian (λ₂)⁹⁸ at the B3LYP/aug-cc-pVDZ level of theory. MultiWFN code⁹⁹ was employed for these RDG analyses.

DFT models of the vapor phase

Rigid rotor-harmonic oscillator (RRHO) model¹⁰⁰ was employed to model gas-phase thermodynamic properties. The gas phase was first modeled as an isolated molecule in a fixed cubic cell of size 20 Å.⁵¹ To ensure consistency with the crystal-phase description, optimization was done at the PBE-D3(BJ)/PAW and PBE-D4/PAW levels of theory and the same VASP versions as for the crystals.⁸⁷ Vibrational contributions to monomer entropy and heat capacity were identified through Phonopy without scaling the calculated frequencies. Contributions from molecular rotation modes were derived from the PBE-D/PAW optimized monomer geometries. The sublimation enthalpy was derived as a function of temperature from the sublimation enthalpy Δ_subH at 0 K and the isobaric heat capacities of ideal gas and crystalline phases C_p:In eqn (2), Δ_subH at 0 K was obtained as a sum of the cohesive energy E_coh (representing the difference between electronic energies of the crystal per unit molecule and of its corresponding monomer) and the difference in zero-point energies ΔE_ZP of the gas and solid phases:⁵¹

Δ_subH(0 K) = −E_coh + ΔE_ZP. (3)

The target molecules do not undergo conformational changes during transition from solid to gas, and hence the energy term for such a conformation relaxation upon the sublimation could be discarded. Saturated sublimation pressure p_sub was then calculated as a function of temperature from the equation:⁵²where p⁰ = 100 kPa, and Δ_subG⁰ and Δ_subS⁰ correspond to standard Gibbs energy and standard entropy changes due to sublimation to the ideal gas at p⁰.

In parallel, calculations for the gas phase were repeated using the state-of-the-art computational model to obtain low-uncertainty gas-phase properties suitable to be used along the experimental results in the development of reference sublimation properties. This time, molecular geometries were optimized at the B3LYP-D3/6-311+G(d,p) level of theory with the empirical D3 dispersion correction²⁰ in Gaussian 16 software, revision C01.¹⁰¹ Carbazole and dibenzothiophene were found to be stable in planar conformation, while for thianthrene and phenothiazine, C_S symmetry with a dihedral angle between the two benzene ring planes of 143° and 129°, respectively, was found to be the energy minimum. All degrees of freedom were treated with RRHO except for the deformation of this dihedral, where the potential energy during the deformation was scanned and represented using a symmetric double-well potential. The contributions to the thermodynamic functions were calculated as for a one-dimensional hindered rotor¹⁰² with a rotational axis passing through the two heteroatoms. The calculated fundamental vibrational frequencies were scaled by a double-linear scaling factor (0.9972–1.48 × 10⁻⁵ν cm⁻¹) and 0.960 for frequencies below and above 2000 cm⁻¹, respectively, developed on experimental vibrational frequencies of n-alkanes.¹⁰³

Material description

A description of the samples used in this study including their purities is given in Table 2. Purity of the commercial sample of CAZ was clearly insufficient both based on its color and our gas-chromatography (GC) analysis. At first, the sample was sublimed twice, but the mole fraction purity improved only to 0.993. Later, we succeeded to purify the sample by gram-scale zone refining. The CAZ sample slowly decomposed when melted and about a half of its volume turned black, but the refined part of the sample was clear white and of very high purity. Experiments with the sublimed CAZ are discussed throughout this paper for a comparison, but all tabulated results were obtained using the purer zone-refined sample. Purities of DBT and TTH determined by the supplier using GC seemed sufficient, but our analysis showed lower purity, around 0.999 and 0.993, respectively. The nonvolatile impurities contained were thus removed by a two-step sublimation, which also eliminated the greyish color of the samples. PTZ was commercially available at sufficient purity (>0.999) and no purification was performed. Sensitivity of PTZ to light and air was suggested by the supplier, therefore all the manipulation with the sample was performed in a glove box under an inert atmosphere.

Table 2 Description of samples of the target materials

Compound	Acronym	Supplier	CAS RN	M/g mol^−1	Original mole fraction purity	Purification method	Final mole fraction purity	Crystal structure^a
a Crystal structure determined at 293 ± 3 K with setup described in Section 2.6 is equivalent to the given CSD entry.b Purity determined by GLC using the Hewlett-Packard 6890A chromatograph equipped with a column HP-1, length 25 m, film thickness 0.52 μm, diameter 0.30 mm and FID detector in the temperature range of 353 to 513 K with heating rate of 20 K min^−1 and 60 minutes isotherm at the end of the heating ramp. Result is an average of two determinations. Purity of 1.0000 is stated if no impurities were detected.c Purity determined from the shape of fusion peaks (heating rate 0.5 K min^−1, for DBT 2 K min^−1) using the van’t Hoff method.^104d Mole fraction purity declared by supplier using GC.
Carbazole	CAZ	Urxovy závody	86-74-8	167.21	0.987^b	Zone refining	1.0000^b; 0.9995^c	CRBZOL03
Dibenzothiophene	DBT	Merck	132-65-0	184.26	1.000^d	Vacuum sublimation	0.9991^b; 0.9988^c	DBZTHP
Thianthrene	TTH	Merck	92-85-3	216.32	0.999^d	Vacuum sublimation	1.0000^b; 0.9985^c	THIANT03
Phenothiazine	PTZ	TCI	92-84-2	199.27	0.999^d	—	1.0000^b; 0.9987^c	PHESAZ01

---

Phase behavior measurements

The international temperature scale ITS90 was used for the calibration purposes and for all measurements. Molar masses of the compounds were calculated based on IUPAC recommendations.¹⁰⁵ The molar gas constant R = 8.314462618 J K⁻¹ mol⁻¹ was used for all calculations.¹⁰⁶

For the phase behavior study, a heat-flux calorimeter TA Discovery DSC 2500 (TA Instruments, New Castle, DE, United States) was used. Operating temperature range of the calorimeter is from 183 to 998 K. To avoid the contact of the sample (about 5 mg) with the atmosphere, aluminum pans were used and samples therein were hermetically sealed by cold welding. Experiments were carried out using heating rates of (20, 10, 5, 2 and 0.5) K min⁻¹ to examine the possibility of samples to crystallize to different forms. In a procedure similar to what we had implemented for previous model of calorimeter by TA Instruments,¹⁰⁷ the calorimeter was calibrated at the listed heating rates with seven standards. Based on the calibration measurements, the expanded uncertainties (k = 2, 0.95 level of confidence) of the measurement are U(T) = 0.3 K and U(ΔH) = 0.03 ΔH.

Heat capacity measurements

Heat capacities were studied using a commercial Tian-Calvet type calorimeter SETARAM Microcalvet (SETARAM, Caluire, France). The measurements were executed in continual heating mode with an isothermal step at the beginning and the end of each heating/cooling cycle. For this work, a temperature range from 235 to 355 K was used with a heating rate of 0.4 K min⁻¹. The device, measurement procedure, and calibration methodology were described in detail previously.¹⁰⁸ Based on the testing experiments, the combined expanded uncertainty of the measured data (k = 2, 0.95 level of confidence) is estimated to be U_c(C⁰_p,m) = 0.006 C⁰_p,m.

To extend the temperature range of our description of the heat capacity of crystalline phases to higher temperatures and to obtain the heat capacity of liquids, a power-compensated differential scanning calorimeter (PC-DSC) PE 8500 (PerkinElmer, Watham, MA, USA) was used. Heat capacities of CAZ, TTH, and PTZ were measured in a temperature range from 303 to 473 K using the temperature increment method (increments of 5 K at a heating rate of 5 K min⁻¹ separated by 1 min isotherms). Measurement and calibration procedures were described in detail previously.¹⁰⁹ Based on the calibration and testing experiments, the standard uncertainty of the temperature is u(T) = 0.1 K and the combined expanded uncertainty of the obtained heat capacities is 0.03C⁰_p,m (k = 2, 0.95 level of confidence). Due to the lower accuracy of the PC-DSC technique, the obtained data were scaled to match those from the Tian-Calvet type calorimeter in line with the common practice.¹¹⁰

X-ray powder diffraction

Crystal forms were identified by X-ray powder diffraction (XRPD), and their respective diffractograms were analyzed with the HighScore Plus software¹¹¹ joined with annually updated powder diffraction databases PDF4+ and PDF4/Organics.¹¹² The XRPD analysis was performed using a θ–θ powder diffractometer X’Pert3 Powder by PANalytical in Bragg–Brentano para-focusing geometry using wavelength CuKα radiation (λ = 1.5418 Å, U = 40 kV, I = 30 mA). The diffractograms were collected at 293.15 ± 3 K in the range of 2θ = (5° to 50°) with a step size of 2θ = 0.039° and a sampling period of 0.7 s.

Static vapor pressure measurements

To measure the vapor pressures, two custom-built static apparatuses (STAT8,¹¹³ with exchanged pressure gauge¹¹⁴ and STAT9¹¹⁵) covering a pressure range from 0.1 to 13330 Pa and a combined temperature range from 283 to 463 K were used. Experiments were carried out in stainless steel cells, which are vacuum-tight fastened to the apparatus. An optimal size of the sample was around one gram. Preliminary depletion of the sample might occur with low amounts, while overloading the cell would lead to sluggish degassing. Durations of the measurements were around one month except for DBT, where the measurement took only one week, as only short temperature range was covered by STAT8, while measurements over the dominant pressure range were executed by means of STAT9 and published elsewhere.¹¹⁵ Overall combined expanded uncertainties (k = 2, 0.95 level of confidence) of both static apparatuses are U_c(p/Pa) = 0.005 p/Pa + 0.05.

Simultaneous correlation of thermal properties

The method of simultaneous correlation (SimCor) enables to simultaneously correlate vapor pressure, fusion properties, and difference between heat capacities of ideal-gas phase and that of a condensed phase Δ^g_cdC⁰_p,m using a single suitable vapor pressure equation. Despite the computational difficulty, this method provides two significant advantages: (i) minimizing the error in description of the input properties over the combined temperature range, which is beneficial in regions where certain properties are immeasurable; and (ii) test of mutual consistency of the input data.

Thermodynamic relations are further described in ref. 71 and a detailed application on the case study of anthracene is presented in ref. 116. The pVT behavior of the gaseous phase was expressed by the means of second virial coefficient estimated by the method by Tsonopoulos¹¹⁷ using critical properties as input parameters. DBT and TTH were described with the term for sulfides, while the term for amides was used for CAZ and PTZ as suggested in Elliot et al.¹¹⁸ Next, Δ^g_cdC⁰_p,m was only correlated at p < 1000 Pa to avoid errors arising from uncertainty of the pVT correction. In this work, the Cox equation¹¹⁹ was used within the SimCor procedure, which can be written as follows:

Experimental results

This section presents a brief overview of reference thermodynamic data that were developed in this work on the basis of own measurements or consistent literature data. A detailed discussion of the development of reference data and their uncertainties is then provided in Section S4 or in the SI.

Temperatures of observed phase transitions and corresponding enthalpies of the materials studied in this work are listed in Table 3. All samples were stable upon melting except CAZ, which exhibited a gradual decomposition that was observable through broadening of the peak of fusion (depicted in Fig. S1). No signs of polymorphism were detected except for phase transitions in CAZ and PTZ observed already in the first runs. For CAZ, a step change of the heat flux was observed with an inflex temperature being highly dependent on the heating rate (documented in Table S1 and in Fig. S2). The reproducible solid–solid phase transition in PTZ exhibited the lambda shape (depicted in Fig. S3) with no significant thermal hysteresis.

Table 3 Calorimetrically determined phase-transition temperatures and enthalpies^a

Compound	i–j	Ti–j^b/K	Treatment	Δ^jiHm^c/kJ mol^−1	Number of repetitions
a Determined at 100 ± 10 kPa using a heating rate of 0.5 K min^−1 unless noted otherwise; uncertainties are expanded uncertainties (k = 2).b Ti–j is the temperature of transition from phase i to j, i.e. crI-l corresponds to a melting point.c Δ^jiHm is the molar enthalpy of transition from phases i to j.d Heating rate of 2 K min^−1 was used.
CAZ	crII-crI	410.0 ± 0.5^d	Temperature at inflex point	0	5
	crI-l	518.8 ± 0.3	Onsets extrapolated to absolute purity	26.5 ± 0.8	10
DBT^d	crI-l	371.7 ± 0.3	Onset corrected to impurity content	22.0 ± 0.7	9
TTH	cr-l	429.9 ± 0.3	Onset corrected to impurity content	27.3 ± 0.9	6
PTZ	crII-crI	249.4 ± 0.4	Peak temperature	0.13 ± 0.03	6
	crI-l	458.8 ± 0.3	Onset corrected to impurity content	25.8 ± 0.8	7

---

Calculated ideal-gas thermodynamic properties used within the SimCor treatment are listed in Table S2. Experimental isobaric heat capacities of condensed phases observed in this work by various methods are listed in Tables S3 and S4, while their trends are depicted in Fig. S4. Final reference data on condensed-phase heat capacities developed in this work can be reproduced with a linear equation:

with parameters a and b listed in Table 4.

Table 4 Parameters of eqn (6) for experimental isobaric heat capacity of condensed phases of the target materials measured in this work

Compound	Phase	Tmin–Tmax	a	b
CAZ	crII	240–410	−7.40905	67.8159
	crI	425–470	−53.7925	82.6645
TTH	l	375–470	143.773	45.2965
PTZ	crI	240–450	6.48684	70.7215
	l	445–470	156.981	41.6457

---

Experimental vapor pressures obtained in this work are presented in Tables S5 and S6. Mutually and internally consistent thermodynamic data sets (literature data printed in bold in Tables S7–S9 and S12, auxiliary data in Table S10, and our novel data listed in Tables S2–S4) were treated by the SimCor procedure to obtain the description of the sublimation equilibrium of the target materials. The resulting parameters of the Cox eqn (5) are summarized in Table 5. A graphical comparison of the trends of experimental sublimation pressures is given in Fig. S5.

Table 5 Parameters of the Cox eqn (5) describing the sublimation equilibrium of the target materials obtained by the SimCor procedure in this work^a

Compound	Phase	(Tmin–Tmax)/K	A0	A1 × 10^3	A2 × 10^6	T0/K	p0/Pa
a Recommended sublimation and/or vaporization enthalpies alongside their estimated uncertainties are listed in Table S12 in the SI.
CAZ	crII	240–410	3.231284	−0.197186	0	519.53	8340.57
	crI	410–519	3.281441	−0.347647	0	519.53	7764.98
	l	519–631	2.965491	−0.422732	0	519.53	7765.41
DBT	cr	240–372	3.456651	−0.179166	0	371.84	33.979
	l	359–662	3.353469	−0.788051	0.429511	371.84	33.979
TTH	cr	240–430	3.419090	−0.179873	0	429.49	198.37
	l	375–630	3.313317	−0.877775	0.470772	429.49	198.37
PTZ	crI	255–458	3.390750	−0.184194	0	458.22	341.45
	l	445–470	3.291155	−0.825985	0.452374	458.22	341.45

---

Discussion

Development of reference densities

Experimental data on crystal phase densities at various temperatures were searched in the CSD database.⁵⁶ Due to the very limited number of available experimental data points at sub-ambient temperatures, only a linear dependence of crystal density on the temperature could be evaluated. Solely mutually consistent data points corresponding to the given crystal form were included (namely for CAZ,^62,120,121 DBT,^{62,74,122,123} PTZ,^75,124,125 and TTH^76,126–128) in the development of reference densities, as illustrated in Fig. S6. Scatter of the crystal-phase densities corresponding to individual experimental ambient-temperature measurements indicates that the current reference densities at 298 K and at standard pressure are burdened with uncertainties ranging from 2 kg m⁻³ (for DBT) up to 11 kg m⁻³ (for CAZ).

Condensed-phase behavior

Table S7 lists an overview of available literature data on the melting temperatures and fusion enthalpies. For DBT and TTH, highly accurate data from adiabatic calorimetry are available, and thus no other literature sources are discussed. Both melting temperatures and fusion enthalpies observed in this work are in agreement within the expanded uncertainties with the values reported for DBT¹²⁹ and TTH.¹³⁰ Our observed data were also in very good agreement with earlier results for CAZ^131–134 and PTZ,^134–137 enabling us to develop reference values of the melting temperatures, 519.0 ± 0.5 K and 458.5 ± 0.3 K, respectively. More details on the uncertainty assessment and data processing are given in the SI in Section S2.

Earlier reports suggest a solid–solid phase transition in CAZ in the region between 420 and 460 K.^131,138 Existence of two orthorhombic polymorphs was indeed confirmed with XRPD measurements and no observed broadening of the NMR spectral lines excluded free rotation of the molecules.¹³¹ Our observation showed that the apparent phase-transition enthalpy changes with the heating rate as well as the phase transition temperature. At low heating rates, the anomaly resembled a glass transition without any relaxation enthalpy, while at 5 K min⁻¹ and higher, some (albeit small) relaxation enthalpy could be evaluated. This behavior satisfies the requirement that the enthalpy difference between crII and crI far from the phase transition must be independent of the heating rate. Heat capacity of crII is lower than that of crI (see Table S3) to the extent that the shift of the phase transition from 410 to 430 K observed with different heating rates corresponds to an enthalpic difference of 0.30 kJ mol⁻¹, which manifests through the ‘relaxation enthalpy’. Microscopic nature of the solid–solid transition in CAZ and its classification (first- or second-order) remains unclear. For further data treatment, we assume that the phases crI and crII have equal enthalpies and Gibbs energies at approximately 410 K.

Literature also contains reports of polymorphism of PTZ which are discussed in detail in the SI in Section S2. In short, a second-order phase transition between a monoclinic and an orthorhombic structure was originally reported,¹³⁹ but a recent study credibly reclassified it as a fist-order transition with only a minor energetic and structural change.⁷⁵ We observed the lambda-shaped peak at 249.4 ± 0.4 K (corresponding to the maximum of the DSC curves). Variation of the phase-transition temperature in literature (between 226 K and 251 K)^75,139–142 can be possibly attributed to shifts of the delicate equilibrium imparted by solid-soluble impurities. Interestingly, this work seems to be the first phase-behavior study of PTZ that reports the sample purity and treatment. The only remaining inconsistency seems to be the report of Bell et al.,¹³⁹ who described a mixture of plate and needle-shaped crystals at room temperature, while it seems that the monoclinic phase described by other authors^75,141,142 may not be superheated to room temperature at atmospheric pressure. The existence of a third PTZ polymorph within the P2₁ group that can be obtained by recrystallization according to Bell et al.¹³⁹ is therefore possible.

Development of reference heat capacities

An overview of literature data available for condensed-phase heat capacities is given in Table S8. For DBT and TTH, the agreement of heat capacities observed in this work with the reliable adiabatic heat capacity data^129,130 corroborates our experimental methodology. For crystalline CAZ, our results are in a reasonable agreement with the earlier results observed by PC-DSC that are burdened with higher uncertainty.¹⁴³ A detailed discussion about the agreement of individual data sets is given in the SI in Section S3. Fig. S7–S9 illustrate mutual agreement of the available heat capacity data sets for CAZ, DBT, and TTH, respectively.

Ideal-gas heat capacities were previously calculated for DBT¹⁴⁴ and PTZ¹³⁷ using the B3LYP/6-31G(d) level of theory and harmonic vibrational frequencies scaled by a single scale factor. Such data sets differ from our calculations by about 2% at 200 K, and the difference vanishes towards high temperatures, as depicted in Fig. S10. This behavior is typical for overcorrected harmonic frequencies that result from using a single scale factor value that is too low for at least the deformation modes.¹⁴⁵ For that purpose, we assume that our current calculations, which use a double-linear scaling of the frequencies, result in ideal-gas heat capacities that are superior to those in the literature in terms of their uncertainties. Note that using separate scaling of the low- and high-frequency proved as beneficial for the accuracy of thermodynamic modeling also for aromatic systems,¹⁴⁵ not only alkanes. Transferability of the double-linear scaling, developed originally for alkanes, also to other systems was verified earlier,¹⁴⁶ and the related computational model was extensively demonstrated to yield thermodynamic properties of the ideal-gas that are highly consistent with other properties relevant to vapor–liquid or sublimation equilibria.^{114,147–150} Our calculated heat capacities of CAZ and DBT are also in a reasonable agreement with an earlier study¹⁵¹ that reported heat capacities based on experimental molecular geometries in crystal and complete vibrational assignments. Those data sets are compared in Fig. S9 and discussed in detail in Section S3.

Development of reference sublimation data

Literature sources available for vapor pressure data for the compounds studied are listed in Table S9. Following the SimCor procedure, we are able to test the mutual consistency among individual data sets. That allows us to exclude obvious outliers as well as vapor pressure data that are inconsistent with the heat capacity difference between the solid and vapor phases or fusion enthalpies from further data correlations. More details about this procedure, including corrections for the non-ideal behavior of the vapor phase (parameters listed in Table S10), data assessment and a discussion on rejection of certain data sets from the data correlations are given in Section 2.8.

Available literature contains several calorimetrically determined sublimation and vaporization enthalpies listed in Table S12, which can be compared with our final reference sublimation data that are depicted in Fig. 2. For CAZ, Fig. S11 shows that our own sublimation-pressure measurements for the purified sample and a single literature data set on vapor pressures¹⁵² were found to be generally consistent with the available heat capacity data. For DBT, the data obtained in this work, by Štejfa et al.,¹¹⁵ Chirico et al.,¹²⁹ and Růžička et al.⁷² are in good agreement and were all included in the final fit as depicted in Fig. S12. A deviation of about 1% between our two static apparatuses (STAT8 and STAT9) in the overlapping temperature range is still covered by the combined uncertainties of the two apparatuses, or exceeds it negligibly. For TTH, the vapor pressure data included in the correlation (Monte et al.,¹⁵³ Steele et al.,¹³⁰ and this work) show excellent mutual agreement in Fig. S13. Notably, three different techniques and three independent laboratories involved yield a mutually perfectly compatible description of TTH sublimation. For PTZ, only the vapor pressure data measured in this work along with the data set reported by Freitas et al.¹³⁷ were found to be sufficiently consistent (Fig. S14).

Fig. 2 Temperature trends of reference sublimation and vaporization enthalpies of the considered compounds as obtained from the SimCor method. Figures correspond from left to right to CAZ, DBT, TTH, and PTZ.

For CAZ, differences between the recommended melting temperature and fusion enthalpy and the SimCor results reach about twice the experimental expanded uncertainty (see Table S11), pointing to a slight inconsistency between the fusion properties and the vapor pressures. The liquid vapor pressures by Senseman and Nelson¹⁵² would be most probably the critical segment. Table S11 further shows that fusion properties of other compounds are described correspondingly to their uncertainties.

For DBT, high-temperature vaporization enthalpies by Mraw and Keweshan,¹⁵⁴ and sublimation enthalpies reported originally by Freitas et al.¹⁴⁴ but recalculated using our reference ideal-gas heat capacities were found to be consistent with our description and were thus included in the final correlation. For TTH and PTZ, none of the literature data on sublimation enthalpies in Table S12 were found to be consistent with the SimCor description, as discussed in detail in Section S2 in the SI. Tabulated sublimation enthalpies described from the SimCor are given in Table S13 and their temperature trends are illustrated in Fig. S15–S18.

An ultimate verification of the mutual compatibility of the independent thermodynamic data sets included in the SimCor data treatment then represents the comparison of ideal-gas entropies from our B3LYP calculations and their experimental counterparts derived from the high-quality adiabatic heat capacity data and the sublimation properties derived from SimCor. These entropic data listed in Table 6 differ by no more than 1.0 J K⁻¹ mol⁻¹ which indicates a very good consistency of our SimCor description of the sublimation equilibrium. Note that any significant difference between those entropy values would either originate from residual entropy at 0 K (related to disorder in the crystal structure) or an error in one of the correlated data sets.

Table 6 Overview of thermodynamic properties of the crystalline phases: absolute molar entropies Δ^298.15K₀S⁰_m(cr), standard sublimation entropies Δ^g_crS⁰_m, and absolute molar entropies for the vapor phase. All entropy data are given in J mol⁻¹ K⁻¹ and are valid at 298.15 K and 100 kPa^a

Compound	Δ^298.15K0S^0m(cr)	Δ^gcrS^0m ^b	S^0m(g)^c	S^0m(g)^d
a Uncertainties given are the expanded uncertainties (k = 2).b SimCor results, values derived from the Cox equation with parameters in Table 5.c Experimental value, S^0m(g) = Δ^298.15K0S^0m(cr) + Δ^gcrS^0m(298.15 K).d Values calculated at the B3LYP-D3 level of theory, see Section 2.2.
CAZ (crII)	NA	186.5 ± 1.5	NA	383.4 ± 1.5
DBT	204.20 ± 0.46^129	186.8 ± 0.6	391.0 ± 0.8	390.4 ± 1.6
TTH	230.89 ± 0.52^130	194.6 ± 0.9	425.5 ± 1.1	424.5 ± 1.7
PTZ (crI)	NA	193.3 ± 1.7	NA	421.0 ± 1.7

---

Finally, for the purpose of benchmarking our DFT computational models, we derived low-uncertainty sublimation enthalpy values at 0 K (Table 7), which represent state-of-the-art reference description of sublimation of the target materials.

Table 7 Overview of thermal contributions to molar enthalpy of the crystalline and vapor phases Δ^298.15K₀H⁰_m(cr) and Δ^298.15K₀H⁰_m(g), and sublimation enthalpies Δ^g_crH⁰_m at 298.15 K and at 0 K. All enthalpy data are given in kJ mol^−1 a

Compound	Δ^298.15K0H^0m(cr)	Δ^gcrH^0m(298.15 K)^b	Δ^298.15K0H^0m(g)^c	Δ^gcrH^0m(0 K)
a Uncertainties given are the expanded uncertainties (k = 2).b SimCor results, values derived from the Cox equation with parameters in Table 5.c Values calculated at the B3LYP-D3 level of theory, see Section 2.2.
CAZ (crII)	NA	105.68 ± 0.61	26.45 ± 0.21	NA
DBT	30.44 ± 0.07^129	93.92 ± 0.23	26.91 ± 0.22	97.45 ± 0.32
TTH	34.36 ± 0.08^130	105.06 ± 0.34	30.77 ± 0.25	108.65 ± 0.43
PTZ (crI)	NA	109.12 ± 0.68	29.84 ± 0.24	NA

---

Microscopic quasi-harmonic modelling

This section presents computational results assembled for the quasi-harmonic modelling of the target crystal phases. Static electronic energies, densities of phonon states, and thence resulting vibrational Helmholtz energies, all calculated explicitly as functions of unit-cell volume using the PBE-D3 and PBE-D4 models, are compared at first to lay the foundations for further discussion of trends among macroscopic thermodynamic observables computed using the PBE-D3 and PBE-D4 models.

Trends in the calculated static electronic energies, E_el, are depicted in Fig. 3. The most important observation regarding the differences between E_el(V_m) curves obtained from PBE-D3 and PBE-D4 models is that PBE-D3 shifts the entire E_el(V_m) curves to smaller molar volumes. The respective difference of the minimum-energy volume coordinate V₀ is less than 1%, but it systematically affects the equilibrium densities of the crystals predicted by the two methods. Apart from this systematic shift, both models predict E_el(V_m) curves very similar in shape and curvature as depicted in Fig. S19. Only a detailed inspection of E_el(V_m) curves expressed in terms of reduced volumes (actual volume divided by the optimal volume) reveals that PBE-D4 exhibits a very weak tendency to yield less steep E_el(V_m) curves (see Fig. S20).

Fig. 3 Illustration of how static electronic energies E_el modelled by PBE-D3 and PBE-D4 for phenothiazine vary with respect to the molar volume V_m of the crystal (left); and a comparison of molar-volume coordinates V₀ of the minimum points of E_el(V_m) curves calculated for the target crystal forms using either PBE or r²SCAN functionals coupled with either D3, D4 or TS dispersion models (right).

Fig. 3 also contains additional results to extend this comparison of the predicted coordinates of E_el(V_m) curve minima, and to verify whether the observed hierarchy between D3 and D4 models holds also when a more sophisticated r²SCAN functional is applied, or a completely different TS dispersion model is adopted. These data confirm that the D4 model systematically yields larger optimized unit-cell volumes than D3 even when the r²SCAN functional is used. All dispersion-corrected r²SCAN results are somewhat lower (by 1–3%) than their respective dispersion-corrected PBE counterparts which seems promising as r²SCAN could yield more accurate quasi-harmonic equilibrium unit-cell volumes, as discussed below. Concerning the TS dispersion model, it is found to yield optimized unit-cell volumes that are by up to 2% higher than D3 when coupled with any of the two considered functionals (except CAZ), but very similar to D4 (differing no more than by 0.6%). A numerical comparison of these six data sets of optimized unit-cell volumes is listed in Table S14.

Dynamic degrees of freedom of the crystals are described in terms of the density of phonon states. Densities of phonon states calculated for the V₀ volumes are shown in Fig. S21. These characteristics of the individual crystal structures are closely similar for both PBE-D3 and PBE-D4 methods and the macroscopic view across the entire wavenumber range does not reveal any systematic phonon-frequency shifts between the models. Marginal variations in the densities of phonon states predicted by both models also lead to very small differences in zero-point vibrational energies of the crystals. These differences in individual zero-point energies amount to less than 0.1 kJ mol⁻¹ with individual zero-point energies listed in Table S15.

Fig. 4 depicts examples of calculated A_vib(V_m) trends for TTH using the PBE-D3 and PBE-D4 models (results for the remaining materials are depicted in Fig. S22). Both levels of theory yield A_vib(V_m) curves very similar in magnitude which is a consequence of the closely resembling phonon densities of state that rules out any large vertical A_vib(V_m) shifts. The actual differences of A_vib values obtained by both models depicted in Fig. 4 match the order of magnitude of the zero-point energy differences computed by the two models. Importantly, the zero-point energy is an increasing function of the phonon frequency, whereas the vibrational Helmholtz energy at finite temperatures is a decreasing function of the phonon frequency due to the behavior of its entropic component. Nonetheless, Fig. S23 reveals that the isothermal A_vib(V_m) functions from both models somewhat differ in their slopes, (∂A_vib/∂V_m)_T. At smaller molar volumes, and thus at elevated pressures or low-enough temperatures, PBE-D3 predicts slightly steeper A_vib(V_m) curves (except PTZ). However, at larger molar volumes (corresponding to low pressures or elevated temperatures), this hierarchy gets inverted with the PBE-D4 (∂A_vib/∂V_m)_T terms being higher in absolute value. Fig. 4 depicts that two of the four target materials follow the former regime (i.e. with steeper PBE-D3 A_vib(V_m) curves) at 298 K and at the standard pressure. Differences in the (∂A_vib/∂V_m)_T terms obtained from both models range from −6% (for CAZ) to 8% (for PTZ). At the standard pressure, such differences remain nearly constant with respect to temperature for DBT and TTH, whereas they steadily rise with temperature for PTZ, as depicted in Fig. S24. The differing A_vib(V_m) slopes observed for both models indicate a distinct performance of the models in describing variations of local curvatures of the potential energy surface, which correspond to individual normal phonon modes, with respect to the crystal volume, and thus also the length of intermolecular contacts. In addition, significant variations of the (∂A_vib/∂V_m)_T term over relevant volume ranges justify the use of non-linear A_vib(V_m) interpolation within our quasi-harmonic procedure.

Fig. 4 Illustration of the vibrational Helmholtz energy of crystalline thianthrene as a function of molar volume at selected temperatures (left) and a comparison of the slopes of A_vib(V_m) functions calculated using either the PBE-D3 or PBE-D4 models for all the target materials at 298 K and at predicted equilibrium volumes corresponding to the standard pressure. Value Δ represents the A_vib(V_m) slope from PBE-D3 minus that from PBE-D4, in J cm⁻³ (right).

Benchmarking structural and thermodynamic properties

Examples of total quasi-harmonic A_tot(T, V_m) functions, governing the thermodynamic behavior of the crystals, are given in Fig. S25. Fig. 5 compares calculated quasi-harmonic equilibrium crystal-phase densities at the standard pressure with their experimental counterparts. Complete results of variable-temperature density calculations are depicted in Fig. S26 and listed in Table S16.

Fig. 5 Illustration of trends of experimental reference densities ρ for dibenzothiophene and those calculated using the quasi-harmonic protocol and PBE-D3 and PBE-D4 models (left); and a comparison of experimental and calculated crystal-phase densities at 298.15 K and at the standard pressure for all the target materials (right).

While both computational models underestimate the densities, PBE-D3 is found to outperform PBE-D4, the former yielding densities systematically higher by up to 2%. This behavior is clearly traceable to the predicted hierarchy of the volume coordinates of the E_el(V) curve minima. Deviations in densities between both models change monotonically with temperature, exhibiting mutual divergence as the temperature increases (with the exception of PTZ, see Fig. S27). The PBE-D3 model yields densities that deviate from the experiment by −3% to −1% at the ambient temperature, as can be seen in Fig. S28. Since r²SCAN calculations yield these volume coordinates of the E_el(V) curve minimum at 1–3% lower volumes, it can be expected that such dispersion-corrected r²SCAN models would yield crystal-phase densities even closer to the experiment than the currently observed PBE data sets.

Coefficients of isobaric thermal expansion, α_p, derived from the temperature trends of calculated densities (Fig. S29 and S30), reveal that the PBE-D4 model predicts stronger thermal expansion (by up to 25%) for crystals for all the target materials except PTZ. These predicted thermal-expansion trends are an imprint of the interplay between the slightly less steep E_el(V_m) curves predicted by PBE-D4 and the hierarchies of the A_vib(V_m) slopes predicted by both models. In particular, the significantly steeper A_vib(V) slope yielded from PBE-D3 for PTZ translates to its more pronounced thermal expansion within that model when compared to PBE-D4.

Fig. 6 illustrates trends for calculated and experimental isobaric heat capacities. Results for CAZ are used there as a demonstrative example of the typical hierarchy among both computational models and the experimental values (see Fig. S31 for analogous comparisons for the remaining materials). In general, relative deviations between PBE-D3 and PBE-D4 predictions do not exceed 5%, and the latter model yields higher heat capacities for all target materials except for PTZ, as illustrated in Fig. S32. At 298 K, the largest deviation between the two methods is seen for CAZ (5.0 J mol⁻¹ K⁻¹ ≈ 2.4%) with the two data sets diverging as temperature increases, reflecting the overstated thermal expansion by the PBE-D4 model. Only for PTZ, the heat capacities predicted by PBE-D4 are lower than the PBE-D3 results, which agrees with its understated thermal expansion from the PBE-D4 model.

Fig. 6 Illustration of the temperature trends of experimental isobaric heat capacities C^cr_p,m for crystalline carbazole and those calculated using the quasi-harmonic protocol and PBE-D3 and PBE-D4 models (left); and a comparison of experimental and calculated crystal-phase heat capacities at 298.15 K and ambient pressure for all the target materials (right).

Heat capacities calculated through the PBE-D3 model are consistently closer to the experimental values at 298.15 K when compared to the PBE-D4 results. Over the temperature range covered by the current experiments, percentage deviations of computed heat capacities from experimental values range from −1.5% to 15% for PBE-D3 and from −2.0 to 20% for PBE-D4. The most outlying heat-capacity predictions belong to DBT. For the other materials, the respective computational errors are appreciably lower, corresponding to a very good computational accuracy in terms of the quasi-harmonic DFT-D heat capacities.^32,33,69,87 For PTZ and TTH, the relatively unimportant temperature trends of these relative deviations and small absolute deviations indicate that our predictions capture their experimental heat capacities very closely over broad temperature intervals, as depicted in Fig. S33. An adverse impact of overestimated thermal expansion for DBT can be seen in the steeply increasing computational errors of its heat capacities at elevated temperatures. Nonetheless, such results are still consistent with the typical accuracy of DFT-D quasi-harmonic models of heat capacities for molecular crystals.^32,33,68 Calculated crystal-phase heat capacities are listed in Table S17.

Concerning the gas-phase heat capacities, Fig. S34 shows the qualitative agreement of both PBE-D models with the reference values computed from scaled B3LYP-D3 vibrational frequencies. Impact of the dispersion correction on intramolecular vibration modes of relatively rigid molecules is generally small.

Benchmarking sublimation enthalpies

Absolute values of calculated cohesive energies |E_coh|, which naturally dominate predictions of the sublimation enthalpies, are listed in Table 8. In most cases, |E_coh| exceed 100 kJ mol⁻¹ which corresponds to the molecular size and the lack of proper hydrogen bonding within the crystal structures. Despite very similar structural and thermodynamic descriptions by the two PBE-D models, resulting E_coh values differ by nearly 4 kJ mol⁻¹ for all the target crystals, with PBE-D3 systematically predicting stronger cohesion. Fig. S35 shows that PBE-D3 gives higher (less negative) electronic energies for monomers than PBE-D4 does (by 4.5–6.1 kJ mol⁻¹), while the analogous differences for the crystal phase are smaller (0.1–2.1 kJ mol⁻¹). Note that these discrepancies in monomer energies between the two models arise indeed from the D3 and D4 dispersion energies. The underlying self-consistent PBE energies differ less than 0.03 kJ mol⁻¹ between the employed VASP5 and VASP6 versions.

Table 8 Calculated cohesive energies E_coh, zero-point contributions to sublimation enthalpies Δ_subE_ZP, sublimation enthalpies Δ_subH at 0 K, thermal contributions to Δ_subH at 298 K, and Δ_subH of target materials at 298 K. All data are computed using the PBE-D3/D4 models and values are given in kJ mol⁻¹

Material	−Ecoh		ΔsubEZP		ΔsubH (0 K)		ΔsubHtherm (298 K)		ΔsubH (298 K)
Dispersion model	D3	D4	D3	D4	D3	D4	D3	D4	D3	D4
CAZ	111.85	107.19	−4.90	−4.82	106.95	102.37	−3.46	−3.95	103.49	98.42
DBT	100.22	96.79	−3.77	−3.74	96.45	93.05	−4.67	−4.94	91.78	88.11
PTZ	120.31	116.16	−4.33	−4.29	115.98	111.87	−3.36	−3.26	112.63	108.61
TTH	112.85	108.85	−4.13	−4.10	108.64	104.75	−3.60	−3.79	105.04	100.96

---

Summation of |E_coh| with zero-point energy contributions Δ_subE_ZP leads to sublimation enthalpies at absolute zero. With differences in Δ_subE_ZP being negligible (<0.1 kJ mol⁻¹), Δ_subH at 0 K naturally retains the hierarchy and systematic offset between the two models observed for the same cohesive energies. Notably, the PBE-D3 model yielded Δ_subH at 0 K within 1 kJ mol⁻¹ or better for DBT and TTH, where the availability of experimental heat capacities of those crystals enabled to evaluate the reference values of Δ_subH at 0 K to be 97.45 kJ mol⁻¹ and 108.65 kJ mol⁻¹ for DBT and TTH, respectively.

While PBE-D4 generally predicted more negative thermal contributions to Δ_subH at 298 K than PBE-D3, differences were rather small (<0.5 kJ mol⁻¹), not altering in turn the trends in Δ_subH established at 0 K. This behavior was demonstrated for PTZ in Fig. 7 where temperature trends for calculated sublimation enthalpies are depicted. Nearly constant vertical offsets of the calculated Δ_subH curves confirm that both models differ predominantly in the cohesive energies, and that the impact of phonon-based thermal contributions is only minor. Fig. S36 depicts in detail that the thermal contributions to Δ_subH range from 2 to 6% at 298 K and always account for less than 11% of total Δ_subH near the upper-temperature bounds of thermodynamic stability of the materials.

Fig. 7 Illustration of the temperature trends of experimental sublimation enthalpy for crystalline phenothiazine and those calculated using the quasi-harmonic protocol and PBE-D3 and PBE-D4 models (left); and a comparison of experimental and calculated sublimation enthalpies at 298.15 K for all the target materials (right).

Fig. S37 shows the temperature trends of calculated Δ_subH in detail, while Table S18 lists Δ_subH calculated at selected temperatures. Fig. 7 illustrates the superior performance of PBE-D3 for predicting Δ_subH, occurring for all the target systems but PTZ. Surprisingly, values of Δ_subH predicted by PBE-D4 for PTZ are in a sub-kJ mol⁻¹ agreement with experiment, being in a stark contrast to the larger errors of 6–8 kJ mol⁻¹ observed for the remaining materials (see Fig. S38). Furthermore, the PBE-D4 model underestimates Δ_subH across all materials at (sub-)ambient temperatures and it further diverges from PBE-D3 at elevated temperatures (again except PTZ, see Fig. S39). These results highlight the challenges of enthalpic calculations that incorporate sophisticated dispersion corrections, such as the D4 model, particularly when coupled with relatively low-tier DFT functionals that are prone to delocalization errors, but allow for reasonable computational costs even for complex periodic systems. Δ_subH values at 298 K calculated by the PBE-D3 model for CAZ and DBT are underestimated, unlike the overestimated PBE-D3 results for PTZ that exhibit the largest error amounting to 3.5 kJ mol⁻¹ (≈3.2%) at 298 K. Considering the systematically lower Δ_subH obtained from the PBE-D4 model, their errors range from 0.5 to −7.3 kJ mol⁻¹.

In the context of other benchmarks of DFT-D-calculated sublimation enthalpies for molecular crystals,^37,51 these generally low computational errors observed for the PBE-D3 model fulfill the chemical accuracy criterion (errors below 4 kJ mol⁻¹). That demonstrates a fair performance of the current PBE-D3 computational model for the target heterocyclic materials. It is also useful to put the evaluated computational accuracy of the current Δ_subH results in the context of a typical scatter of individual literature experimental data entries on Δ_subH at 298 K originating from distinct laboratories. Mean average differences of available literature Δ_subH values from the current SimCor results amount to 8.3, 3.1, 4.2, and 3.6 kJ mol⁻¹ for CAZ, DBT, TTH, and PTZ, respectively. Comparing such an experimental scatter with the observed computational errors indicates that the accuracy of periodic DFT-D computational methods can approach the expected uncertainty of isolated experimental literature data entries for low-volatile materials that are similar to those targeted in this work. It needs to be emphasized, however, that the experimental uncertainties can be further lowered to the sub-chemical-accuracy region through the data-demanding SimCor approach and related critical data assessment.

Variation of Δ_subH with respect to temperature is described qualitatively well by both models. The slopes of reference and calculated Δ_subH(T), corresponding to the magnitude of a Δ_subC_p term, are mostly within 15 J K⁻¹ mol⁻¹ from each other, except for DBT at 350 K where computational errors for PBE-D3 and PBE-D4 reach 25 and 35 J K⁻¹ mol⁻¹ (≈79% and 111%), respectively.

The observed performance of periodic PBE-D3 and PBE-D4 predictions of sublimation enthalpies can be interpreted in terms of the accuracy of the individual interaction energies predicted by these models. In particular, one should focus on proximate molecular dimers and trimers that can be extracted from the optimized crystal structures. Given the dominance of short-range dispersion interactions in these materials, subsequent analyses will focus only on the first coordination shell around a reference molecule from each crystal structure, which spans molecular contacts roughly up to 3.5 Å.

Fig. 8 compares the pair-interaction energies E_int of proximate dimers calculated using the reference method DLPNO-CCSD(T)/CBS (which mimics the gold standard for modelling non-covalent interactions) with values predicted from appropriate combinations of the PBE and r²SCAN functionals and the three dispersion-correction models D3, D4, and TS. Within the narrow first coordination shell, individual values of E_int do not display a clear relationship with pair distances as mutual molecular orientation is very important at short distances. For close intermolecular contacts, corresponding to the first coordination shell, the full resolution of individual atoms to model the molecular interactions is very important as interaction energies depend in such cases in particular on which functional groups get into the closest contact. For instance, assuming a two-molecular cluster at a fixed short contact distance, its interaction energy will strongly depend on any of the following features: (i) proximity of multiple electronegative moieties; (ii) proximity of an electronegative and an electropositive moiety; (iii) presence of either a parallel or antiparallel alignment of electric dipoles; (iv) proximity of a strong permanent multipole with a highly polarizable moiety, and so forth. These orientational factors are extremely important as they contribute not only to electrostatic interactions, but also to induction, dispersion and exchange effects that are known for their short-range character.

Fig. 8 Comparison of the closest pair interactions extracted from the target crystal structures calculated using either PBE or r²SCAN functional and either D3, D4 or TS dispersion corrections with the reference results of DLPNO-CCSD(T)/CBS calculations. Geometries of the dimers corresponding to the strongest attractive interaction in individual crystal structures are depicted as insets.

For the most intense attractive pair interactions exceeding roughly −20 kJ mol⁻¹ in magnitude, PBE-D3 (more significantly) and PBE-D4 (slightly) underbinds the pair interactions with respect to the DLPNO-CCSD(T) reference results, whereas PBE-TS overbinds these dimers for all the considered materials. Interestingly, pair interactions modeled with r²SCAN are predominantly slightly less attractive than their PBE counterparts. The employed r²SCAN-TS model yields the most intense attractions among the considered r²SCAN data sets. However, the hierarchy between r²SCAN-D3 and r²SCAN-D4 pair interactions is predominantly inverted to that observed for PBE results, now with r²SCAN-D3 yielding the more bound dimers.

These trends are further underscored in Fig. S40, confirming that the r²SCAN-D4 and PBE-TS data sets represent predominantly the least and the most intense attractive interactions. Table 9 lists root-mean-squared deviations of all the considered DFT-D pair interactions energies from the reference DLPNO-CCSD(T) results. Based solely on this metric, the PBE-D4 theory could be assessed as the most accurate to describe the dimers, mimicking the reference values at the closest, whereas r²SCAN-D4 differs the most. Absolute pair interactions calculated using all DFT-D methods and the reference theory are once more compared graphically in a diagonal plot in Fig. S41. A complete list of the analyzed pair interaction energies and related mutual molecular orientations in all considered dimers is given in Table S19.

Table 9 Root-mean-squared deviations (in kJ mol⁻¹) of pair-interaction energies of the proximate dimers extracted from the first coordination shell of PBE-D3 optimized crystal structures of the target materials and calculated at individual DFT-D levels of theory from reference DLPNO-CCSD(T) results

Method	PBE-D3	PBE-D4	PBE-TS	r^2SCAN-D3	r^2SCAN-D4	r^2SCAN-TS
2-Body	1.56	1.02	2.06	2.15	2.62	1.31
3-Body	0.47	1.07	0.32	0.53	0.50	0.73

---

Focusing on dimers exhibiting the most intense attractive interactions, Fig. S42 depicts the so-called reduced density gradient (RDG) plots which enable the identification of the most important contact points within a cluster of interacting molecules. Comparison of the coordinates of the RDG critical points (sharp peaks where RDG approaches zero values) enables to characterize the most important cohesive features in the crystals. For DBT, such RDG critical points occur between −0.01 and 0.00 a.u. in terms of the signed electron density, which corresponds to weak dispersion attractions due to C–H⋯π contacts. Similar positions of the RDG critical points can be observed also for CAZ and TTH, being shifted slightly below −0.01 a.u., which rules out any proper hydrogen bonding in CAZ crystals. For PTZ, an RDG critical point is present at the lowest signed electron density, but still above −0.02 a.u., corresponding to the N–H⋯S hydrogen bond of a very weak intensity in its crystal.

The dimer analysis reveals that the PBE-D3 level of theory underbinds individual pair interactions (with respect to both PBE-D4 and coupled-clusters results), although PBE-D3 yields larger sublimation enthalpies than PBE-D4. That leaves space for additional interactions to cause the sublimation enthalpy to be underestimated by the PBE-D4 model. A similar analysis of three-body interaction energies of proximate molecules,^38,155 i.e. interaction terms contributing to the energy of a full trimer beyond the molecular-pairwise interactions, can identify systematic trends in treatment of molecular many-body interactions by individual PBE-D methods.

Fig. 9 illustrates that three-body interactions composed of the dimers from within the first-coordination shell exhibit a hierarchy distinct from what was observed for the pair interactions. In accordance with the Axilrod–Teller–Muto (ATM) model for three-body dispersion interactions of individual atoms,^48,49 also the currently observed three-body molecular interaction terms in the target crystals are predominantly positive, thus representing a repulsive contribution towards the overall crystal cohesion. Importantly, the repulsive character of molecular three-body terms effectively increases in a row: PBE-TS < PBE-D3 < PBE-D4, all such results being more repulsive than the reference DLPNO-CCSD(T)/CBS theory. On the other hand, r²SCAN results on the three-body terms are predominantly less repulsive than the reference theory with the dispersion models following the same hierarchy. The least and the most repulsive characteristics of the trimers from the PBE-D4 and r²SCAN-TS models are further emphasized in in Fig. S43 and S44. Note that only the D4 model contains the ATM terms by default, so that the observed hierarchies of DFT-D molecular three-body terms is also affected by three-body induction effects which can be, unlike the dispersion, handled by the uncorrected DFT functionals.

Fig. 9 Comparison of the closest three-body interactions extracted from the target crystal structures calculated using either PBE or r²SCAN functionals and either D3, D4 or TS dispersion model with reference DLPNO-CCSD(T)/CBS results. The three-body distance of an ABC trimer is defined using the closest-contact distances in its individual dimers: d = (d_AB³·d_AC³·d_BC³)^1/9.

RMSD values of DFT-calculated three-body terms from the reference data set, listed in Table 9, indicate that the PBE-TS exhibits the best performance in this case although that theory was not designed to explicitly treat three-body interactions. Such an observation can be thus regarded as a coincidence arising from fortuitous error cancellation. On the other hand, the considered levels of theory (including the D4 model) yield RMSD values larger by a factor of 2 to 3. The excessive three-body repulsion in the D4 model can be traced to its ATM component. Note that the ATM term is turned off in the PBE-D3 model by default. Clearly, these overly repulsive three-body terms in the D4-containing models contribute significantly to the underestimation of the sublimation enthalpies for all the four materials we considered.

Benchmarking sublimation pressures

Quasi-harmonic calculations for the crystal phase and RRHO calculations for the vapor phase also allow a temperature-dependent description of the standard sublimation entropy Δ_subS⁰. Fig. 10 shows a comparison of the calculated and the reference sublimation entropies at 298.15 K among the target materials, while Fig. S45 and S46 present the trends and corresponding errors over a broader temperature range. In all cases, the vapor phase is described as an ideal gas at the standard pressure of 100 kPa.

Fig. 10 Illustration of the temperature trends of experimental standard sublimation entropy for crystalline thianthrene and those calculated using the quasi-harmonic protocol and PBE-D3 and PBE-D4 models (left); and a comparison of experimental and calculated standard sublimation entropies at 298.15 K and 100 kPa for all the target materials (right).

The PBE-D3 and PBE-D4 models underestimate Δ_subS⁰ by 3–9% and 4–11%, respectively, with the computational errors predominantly increasing with temperature. PTZ is the only material where errors of calculated Δ_subS⁰ remain nearly constant around −7%, being almost independent of the choice of a particular PBE-D model. On the contrary, PBE-D4 yields lower Δ_subS⁰ values (thus less accurate) than PBE-D3 over broad temperature ranges for CAZ, DBT, and TTH, with the mutual differences reaching up to 4% (see Fig. S47). The larger errors observed in the Δ_subS⁰ obtained through PBE-D4 are a consequence of the overestimated quasi-harmonic heat capacities of the crystal phases. Results of both Δ_subS⁰ models then naturally diverge with temperature, although to a lesser extent than the sublimation enthalpy.

Finite-temperature entropy of a material depends on its phonon frequencies, with low-frequency modes contributing more importantly. As long as both computational models yield very similar phonon frequencies, also thence resulting sublimation entropies are very close. Computational uncertainties related to the low-frequency lattice and intramolecular modes of molecular crystals typically reach tens and low-units of percent, respectively,³² imparting the observed large errors of calculated sublimation entropies. In addition, differences between the treatments of ideal-gas entropy within the computational PBE data sets and the reference data set also contribute to this observed off-set.

Benchmarking sublimation pressures represents the most stringent test of accuracy for first-principles models of sublimation.⁵² Table S20 lists the calculated sublimation pressures at selected temperatures while Fig. S48 illustrates the comparison of trends of calculated and experimental sublimation pressures. Fig. 11 shows that PBE-D3 consistently underestimates p_sub across all materials at 298 K, whereas PBE-D4 does so only for PTZ. At the same time, PBE-D4 predicts p_sub values that are larger at least by a factor of three when compared to those from PBE-D3. The computational accuracy of both PBE-D3 and PBE-D4 p_sub results, exhibiting root-mean-squared deviations from the experiment amounting to 64% and 112%, respectively, can be then accepted as very good, especially in the context of the immense computational sensitivity of the sublimation pressures to any computational noise.⁵²

Fig. 11 Illustration of the temperature trends of experimental sublimation pressure for crystalline carbazole and those calculated using the quasi-harmonic protocol and PBE-D3 and PBE-D4 models (left); and a comparison of experimental and calculated sublimation pressures at 298.15 K for all the target materials (right).

The trends observed among calculated p_sub can be explained in terms of an interplay between the computational errors in the enthalpic and entropic contributions to the Gibbs energy of sublimation. The prevailing overestimation of PBE-D4 pressures is a consequence of underestimation of the related sublimation enthalpies. Such 4–7 kJ mol⁻¹ errors in the sublimation enthalpy can easily propagate to an order-of-magnitude error of the sublimation pressure.⁵² In this case, however, the concurrent underestimation of both the sublimation enthalpy and sublimation entropy represents a fortuitous coincidence, leading to large error cancellation between both terms.

As a consequence, the overall errors of p_sub values predicted with PBE-D4 reach factors of 2–3 rather than entire orders of magnitude. Fig. S49 presents an overview of the computational errors in the Gibbs energy of sublimation that are predominantly appreciably lower (by a factor of 2–3) than the corresponding sublimation enthalpy errors. Underestimation of calculated sublimation enthalpy, nevertheless, propagates to too mild slopes of the p_sub(T) functions, which is best illustrated in Fig. S48 by the PBE-D4 results for CAZ and DBT. Similarly, higher sublimation enthalpies from PBE-D3 translate to lower p_sub values than what PBE-D4 predicts.

Impact of the errors in sublimation entropy can be best illustrated for the case of PBE-D3 results for TTH. Its sublimation enthalpy is captured with a very small error (less than 0.1 kJ mol⁻¹), but its sublimation entropy is underestimated (6.6 J K⁻¹ mol⁻¹). Altogether, a nearly 2 kJ mol⁻¹ error in the standard Gibbs energy of sublimation at 298 K corresponds to a 56% error in the sublimation pressure. Similar observations on the entropic effect can be made for the PBE-D4 results for PTZ.

Given that the orders of magnitude of the sublimation pressures at 298 K range from 10⁻⁵ to 10⁻² Pa, the absolute errors of predicted p_sub amount to fractions of millipascals. Only for DBT exhibiting the highest volatility, such computational errors reach low tens of millipascals.

Finally, it remains to link the hierarchy of the modelled volatility of individual target materials with their cohesive interactions. Fig. S50 shows very strong correlations between the sublimation enthalpy (sublimation pressure likewise) of a bulk material and the energy of the (three) most intense attractive pair interactions in its crystal structure. The absence of the –NH– linker in DBT and TTH leaves out any potential for significant hydrogen bonding, which lowers their cohesion and augments their volatility.

In particular, PTZ is the least volatile material in our set (despite being a lighter molecule than TTH). PTZ also exhibits the strongest pair interaction dominated by the N–H⋯S hydrogen bond in its crystal structure, as illustrated in Fig. 8. The opposite holds for the most volatile material in our set, DBT. Its strongest dimer is governed only by substantially weaker C–H⋯π interactions. The large bend of the TTH molecule seems to rule out a simple herringbone packing of its crystal. Instead, a distinct packing motif with two adjacent molecules being interlocked, leaving a vast potential for both π–π and C–H⋯π interactions (depicted in Fig. 8), leads to a larger cohesion than what is observed in DBT. Comparing the interactions and properties within the CAZ–PTZ pair and within the DBT–TTH pair, where the latter molecules always contain one more sulfur atom, the most intense pair interactions within the latter (thus heavier) material are always stronger, and the volatility of the latter is always lower.

Occurrence of the N–H⋯S hydrogen bond in crystalline PTZ is an important feature that distinguishes this crystal from the remaining target materials. Concurrently, PTZ behaves as an outlier within our benchmark, exhibiting mostly opposite trends to the other materials. Importantly, it is also the only system where the PBE-D4 model captures all the considered properties better than PBE-D3. Our results thus indicate that the PBE-D4 model might better capture the interplay of dispersion and hydrogen bonding than the older, simpler PBE-D3 model. In the remaining target systems governed merely by dispersion, the imperfection of the low-tier PBE functional seems to fortuitously cancel out with the less sophisticated D3 correction, so that PBE-D3 is observed to perform better than PBE-D4.

For the reader's convenience, a summary of the main computational results of the PBE-D3 and PBE-D4 benchmark, root-mean-squared deviations of thermodynamic properties relevant for computational modeling of the sublimation equilibrium, is listed in Table 10.

Table 10 Root-mean-squared percentage deviations (RMSD) of quasi-harmonic properties of the target crystals at 298.15 K calculated in this work using either the PBE-D3 or PBE-D4 periodic quasi-harmonic models: density ρ, isobaric heat capacity C^cr_p, sublimation enthalpy Δ_subH, standard sublimation entropy Δ_subS⁰, and saturated sublimation pressure p_sub

Property	ρ		Cp^cr		ΔsubH		ΔsubS^0		psub
Model	D3	D4	D3	D4	D3	D4	D3	D4	D3	D4
a RMSD computed only for CAZ and DBT, containing a five-membered heterocycle in between two six-membered carbon rings.b RMSD computed only for PTZ and TTH, containing a six-membered heterocycle in between two six-membered carbon rings.
CAZ & DBT^a	2.6%	3.9%	4.6%	6.8%	2.2%	6.5%	6.5%	8.0%	44%	143%
PTZ & TTH^b	1.0%	1.6%	0.4%	1.0%	2.3%	2.8%	5.7%	6.1%	77%	78%
Global	2.0%	3.0%	3.3%	4.9%	2.2%	5.0%	6.1%	7.1%	63%	116%

---

Since the equilibrium volumes of the crystals are systematically burdened with an error ranging to low percentage units, the impact of this volume miscalculation on the computational accuracy of the remaining targeted thermodynamic properties should be discussed. A principal feature of the quasi-harmonic approximation is that the thermodynamic functions, such as entropy or heat capacity of a single phase depend largely on the slope of the A_vib(V) function, but not that importantly on its A_vib value itself. Results shown in Fig. S22 indicate that the vibrational Helmholtz energy behaves nearly linearly with respect to volume when one considers the width of a volume range corresponding to typical computational errors of quasi-harmonic crystal-phase volumes. This linearity over short volume ranges also leads to minimization of the entropy errors related to miscalculations of the volume. Computational errors of isobaric heat capacity, on the other hand, are affected significantly more by the quality of the description of thermal expansivity of a material. Finally, modelling phase equilibria depends on the actual Helmholtz energy values of the supposedly co-existing phases. However, in the case of polymorphism, one can assume that errors in the computed volumes are comparable across multiple polymorphs treated with the same theory, allowing for a massive error cancellation. In the case of sublimation, the actual lattice energy error due to the volume miscalculations can be estimated from the shapes of the electronic energy – volume curves. For the currently considered materials, Fig. S20 illustrates that this 3% error in volume can correspond roughly to 0.5–1.0 kJ mol⁻¹ error in the lattice energy.

Conclusions

In this work, we performed a thorough revision of thermodynamic data for four similar heterocyclic aromatic materials with potential for acting as precursors for organic semiconductors, namely carbazole, dibenzothiophene, phenothiazine, and thianthrene. Our efforts included extensive novel calorimetric and vapor pressure measurements, leading to the development of low-uncertainty reference sublimation data for these materials. Furthermore, we performed a benchmark analysis of the computational performance of multiple periodic dispersion-corrected DFT methods, relying either on PBE or r²SCAN functionals and either D3, D4 or TS dispersion corrections for modeling the structure and interactions within those heterocyclic molecular crystals. Although the accuracy of dispersion-corrected r²SCAN models in terms of equilibrium crystal densities was slightly better than for the PBE-D models, extensive benchmarking of quasi-harmonic models of thermodynamic properties of the crystals and their sublimation at finite temperatures was further performed only at the PBE-D3 and PBE-D4 levels of theory due to high computational cost of the quasi-harmonic r²SCAN models. Our results revealed that the presence of sulfur atoms in the heterocycles does not lead to sudden failures of the quasi-harmonic PBE-D models. The computational accuracy of these models for sulfur-heterocyclic materials was demonstrated to be fully comparable to what was common for other molecular materials benchmarked earlier.

The accuracy of the PBE-D4 model was, nevertheless, observed to be somewhat lower for materials composed of molecules with multiple differently sized condensed rings, i.e. carbazole and dibenzothiophene, which concurrently contain both five- and six-membered rings. On the other hand, PBE-D4 performed better for the sole system where an interplay of dispersion and weak hydrogen bonding occurred. Among the four target crystalline materials, the performance of PBE-D3 was observed to be better on average than PBE-D4 for modeling equilibrium densities, thermal expansion, heat capacities, sublimation enthalpies, and sublimation pressures. Although differences between both models are smaller than the actual computational errors, the PBE-D4 model yields computational errors larger by 65% (averaged over the benchmarked properties) than PBE-D3 does. Notably, the PBE-D4 model predicts excessively repulsive three-body interactions in the target crystal structure, translating to systematically underestimated cohesion of those molecular crystals. The better performance observed with PBE-D3 is, however, expected to be a consequence of a fortuitous cancellation of errors due to the imperfect DFT functional and less sophisticated dispersion-correction model.

Author contributions

Reynaldo Geronia II – investigation, formal analysis, writing – original draft, software, visualization. Štefan Kocian – investigation, formal analysis, writing – original draft. Vojtěch Štejfa – methodology, resources, supervision, software, validation, writing – review & editing, formal analysis, investigation, data curation, Ctirad Červinka – conceptualization, methodology, resources, supervision, software, formal analysis, data curation, funding acquisition, project administration, writing – review & editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

All data supporting the conclusions of this work are included within the SI.

Supplementary material 1 (.pdf file) contains: S1: details on experimental results; S2: discussion of reference thermodynamic data and phase behavior; S3: details on results of the quasi-harmonic model. Supplementary material 2 (.xml file) contains: machine-readable experimental and first-principles calculated data from Table 3 and Tables S1–S7 in the ThermoML format. See DOI: https://doi.org/10.1039/d5cp01485a

Acknowledgements

The authors acknowledge financial support from the Czech Science Foundation (GACR no. 23-05476M). Computational resources were provided by the e-INFRA CZ project (ID:90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic. R. G. and Š. K. acknowledge funding from the grant of specific university research A1_FCHI_2025_001 and A2_FCHI_2025_031.

References

1. S. L. Price, D. E. Braun and S. M. Reutzel-Edens, Chem. Commun., 2016, 52, 7065–7077.
2. C. Wang, H. Dong, L. Jiang and W. Hu, Chem. Soc. Rev., 2018, 47, 422–500.
3. J. Bernstein, Polymorphism in Molecular Crystals, Oxford University Press, 2020.
4. G. Schweicher, G. Garbay, R. Jouclas, F. Vibert, F. Devaux and Y. H. Geerts, Adv. Mater., 2020, 32, 1905909.
5. C. Červinka, CrystEngComm, 2025, 27, 2778–2794.
6. G. J. O. Beran, Chem. Rev., 2016, 116, 5567–5613.
7. H. Bronstein, C. B. Nielsen, B. C. Schroeder and I. McCulloch, Nat. Rev. Chem., 2020, 4, 66–77.
8. A. Y. Sosorev, M. V. Vener, O. G. Kharlanov, E. V. Feldman, O. V. Borshchev, N. I. Sorokina, T. V. Rybalova, S. A. Ponomarenko and D. Y. Paraschuk, J. Phys. Chem. C, 2023, 127, 17948–17957.
9. K. E. Riley, M. Pitoňák, P. Jurečka and P. Hobza, Chem. Rev., 2010, 110, 5023–5063.
10. M. O. Sinnokrot and C. D. Sherrill, J. Phys. Chem. A, 2006, 110, 10656–10668.
11. D. B. Ninković, J. P. Blagojević Filipović, M. B. Hall, E. N. Brothers and S. D. Zarić, ACS Cent. Sci., 2020, 6, 420–425.
12. E. A. Meyer, R. K. Castellano and F. Diederich, Angew. Chem., Int. Ed., 2003, 42, 1210–1250.
13. R. G. Huber, M. A. Margreiter, J. E. Fuchs, S. von Grafenstein, C. S. Tautermann, K. R. Liedl and T. Fox, J. Chem. Inf. Model., 2014, 54, 1371–1379.
14. C. D. Sherrill, Acc. Chem. Res., 2013, 46, 1020–1028.
15. M. O. Sinnokrot and C. D. Sherrill, J. Phys. Chem. A, 2003, 107, 8377–8379.
16. T. Togo, L. Tram, L. G. Denton, X. ElHilali-Pollard, J. Gu, J. Jiang, C. Liu, Y. Zhao, Y. Zhao, Y. Zheng, Y. Zheng, J. Yang, P. Fan, M. R. Arkin, H. Härmä, D. Sun, S. S. Canan, S. E. Wheeler and A. R. Renslo, J. Med. Chem., 2023, 66, 9784–9796.
17. M. Klajmon and C. Červinka, J. Chem. Theory Comput., 2021, 17, 6225–6239.
18. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
19. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473.
20. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
21. S. Grimme, A. Hansen, J. G. Brandenburg and C. Bannwarth, Chem. Rev., 2016, 116, 5105–5154.
22. Y. Huang, Y. Shao and G. J. O. Beran, J. Chem. Phys., 2013, 138, 224112.
23. M. Pitoňák and A. Heßelmann, J. Chem. Theory Comput., 2010, 6, 168–178.
24. J. Řezáč, C. Greenwell and G. J. O. Beran, J. Chem. Theory Comput., 2018, 14, 4711–4721.
25. J. Řezáč and P. Hobza, J. Chem. Theory Comput., 2013, 9, 2151–2155.
26. C. Riplinger, B. Sandhoefer, A. Hansen and F. Neese, J. Chem. Phys., 2013, 139, 134101.
27. J. Garcia, R. Podeszwa and K. Szalewicz, J. Chem. Phys., 2020, 152, 184109.
28. T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno and C. D. Sherrill, J. Chem. Phys., 2014, 140, 094106.
29. C. Červinka and G. J. O. Beran, Chem. Sci., 2018, 9, 4622–4629.
30. J. Yang, W. Hu, D. Usvyat, D. Matthews, M. Schütz and G. K.-L. Chan, Science, 2014, 345, 640–643.
31. J. Moellmann and S. Grimme, J. Phys. Chem. C, 2014, 118, 7615–7621.
32. C. Červinka, M. Fulem, R. P. Stoffel and R. Dronskowski, J. Phys. Chem. A, 2016, 120, 2022–2034.
33. J. Ludík, V. Kostková, Š. Kocian, P. Touš, V. Štejfa and C. Červinka, J. Chem. Theory Comput., 2024, 20, 2858–2870.
34. V. Pokorný, P. Touš, V. Štejfa, K. Růžička, J. Rohlíček, J. Czernek, J. Brus and C. Červinka, Phys. Chem. Chem. Phys., 2022, 24, 25904–25917.
35. P. Touš and C. Červinka, Cryst. Growth Des., 2023, 23, 4082–4097.
36. J. L. McKinley and G. J. O. Beran, Faraday Discuss., 2018, 211, 181–207.
37. G. A. Dolgonos, J. Hoja and A. D. Boese, Phys. Chem. Chem. Phys., 2019, 21, 24333–24344.
38. J. Hoja, A. List and A. D. Boese, J. Chem. Theory Comput., 2024, 20, 357–367.
39. A. M. Reilly and A. Tkatchenko, J. Chem. Phys., 2013, 139, 024705.
40. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
41. K. Lee, É. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
42. O. A. Vydrov and T. Van Voorhis, J. Chem. Phys., 2010, 133, 244103.
43. J. Řezáč, K. E. Riley and P. Hobza, J. Chem. Theory Comput., 2011, 7, 2427–2438.
44. D. J. Carter and A. L. Rohl, J. Chem. Theory Comput., 2014, 10, 3423–3437.
45. E. R. Johnson, in Non-Covalent Interactions in Quantum Chemistry and Physics, ed. A. Otero de la Roza and G. A. DiLabio, Elsevier, 2017, pp. 169–194 DOI:10.1016/B978-0-12-809835-6.00006-2.
46. A. Ambrosetti, A. M. Reilly, R. A. DiStasio, Jr. and A. Tkatchenko, J. Chem. Phys., 2014, 140, 18A508.
47. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
48. B. M. Axilrod and E. Teller, J. Chem. Phys., 1943, 11, 299–300.
49. Y. Muto, Proc. Phys.-Math. Soc. Jpn., 1943, 17, 629–631.
50. E. Caldeweyher, J.-M. Mewes, S. Ehlert and S. Grimme, Phys. Chem. Chem. Phys., 2020, 22, 8499–8512.
51. C. Červinka and M. Fulem, J. Chem. Theory Comput., 2017, 13, 2840–2850.
52. C. Červinka and M. Fulem, Cryst. Growth Des., 2019, 19, 808–820.
53. L. Maschio, B. Civalleri, P. Ugliengo and A. Gavezzotti, J. Phys. Chem. A, 2011, 115, 11179–11186.
54. S. P. Thomas, P. R. Spackman, D. Jayatilaka and M. A. Spackman, J. Chem. Theory Comput., 2018, 14, 1614–1623.
55. K. R. Bryenton, A. A. Adeleke, S. G. Dale and E. R. Johnson, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2023, 13, e1631.
56. C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 171–179.
57. A. Mishra, C.-Q. Ma and P. Bäuerle, Chem. Rev., 2009, 109, 1141–1276.
58. K. Takimiya, S. Shinamura, I. Osaka and E. Miyazaki, Adv. Mater., 2011, 23, 4347–4370.
59. S. Revoju, A. Matuhina, L. Canil, H. Salonen, A. Hiltunen, A. Abate and P. Vivo, J. Mater. Chem. C, 2020, 8, 15486–15506.
60. M. Fu, C. Zhang, Y. Chen, K. Fan, G. Zhang, J. Zou, Y. Gao, H. Dai, X. Wang and C. Wang, Chem. Commun., 2022, 58, 11993–11996.
61. V. Malytskyi, J.-J. Simon, L. Patrone and J.-M. Raimundo, RSC Adv., 2015, 5, 354–397.
62. C. Chen, Z. Chi, K. C. Chong, A. S. Batsanov, Z. Yang, Z. Mao, Z. Yang and B. Liu, Nat. Mater., 2021, 20, 175–180.
63. R. P. Stoffel, C. Wessel, M.-W. Lumey and R. Dronskowski, Angew. Chem., Int. Ed., 2010, 49, 5242–5266.
64. J. Nyman and G. M. Day, Phys. Chem. Chem. Phys., 2016, 18, 31132–31143.
65. J. George, R. Wang, U. Englert and R. Dronskowski, J. Chem. Phys., 2017, 147, 074112.
66. Y. N. Heit, K. D. Nanda and G. J. O. Beran, Chem. Sci., 2016, 7, 246–255.
67. J. G. Brandenburg, J. Potticary, H. A. Sparkes, S. L. Price and S. R. Hall, J. Phys. Chem. Lett., 2017, 8, 4319–4324.
68. C. Červinka, J. Comput. Chem., 2022, 43, 448–456.
69. C. Červinka and V. Štejfa, ChemPhysChem, 2020, 21, 1184–1194.
70. J. Hoja, H.-Y. Ko, M. A. Neumann, R. Car, R. A. DiStasio and A. Tkatchenko, Sci. Adv., 2019, 5, eaau3338.
71. K. Růžička and V. Majer, AIChE J., 1996, 42, 1723–1740.
72. K. Růžička, I. Mokbel, V. Majer, V. Růžička, J. Jose and M. Zábranský, Fluid Phase Equilib., 1998, 148, 107–137.
73. K. Gajda, B. Zarychta, K. Kopka, Z. Daszkiewicz and K. Ejsmont, Acta Crystallogr., Sect. C: Struct. Chem., 2014, 70, 987–991.
74. D. Yamazaki, T. Nishinaga and K. Komatsu, Org. Lett., 2004, 6, 4179–4182.
75. D. de Sa Pereira, D. R. Lee, N. A. Kukhta, K. H. Lee, C. L. Kim, A. S. Batsanov, J. Y. Lee and A. P. Monkman, J. Mater. Chem. C, 2019, 7, 10481–10490.
76. X. Cai, Z. Qiao, M. Li, X. Wu, Y. He, X. Jiang, Y. Cao and S.-J. Su, Angew. Chem., Int. Ed., 2019, 58, 13522–13531.
77. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
78. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
79. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
80. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
81. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
82. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
83. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
84. J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, J. Phys. Chem. Lett., 2020, 11, 8208–8215.
85. S. Ehlert, U. Huniar, J. Ning, J. W. Furness, J. Sun, A. D. Kaplan, J. P. Perdew and J. G. Brandenburg, J. Chem. Phys., 2021, 154, 061101.
86. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
87. C. Červinka, M. Klajmon and V. Štejfa, J. Chem. Theory Comput., 2019, 15, 5563–5578.
88. F. D. Murnaghan, Proc. Natl. Acad. Sci. U. S. A., 1944, 30, 244–247.
89. K. Parlinski, Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett., 1997, 78, 4063–4066.
90. A. Togo, J. Phys. Soc. Jpn., 2022, 92, 012001.
91. S. Grimme, A. Hansen, S. Ehlert and J.-M. Mewes, J. Chem. Phys., 2021, 154, 064103.
92. C. Červinka, M. Fulem and K. Růžička, J. Chem. Phys., 2016, 144, 064505.
93. F. Neese and E. F. Valeev, J. Chem. Theory Comput., 2011, 7, 33–43.
94. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108.
95. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1606.
96. G. L. Stoychev, A. A. Auer and F. Neese, J. Chem. Theory Comput., 2017, 13, 554–562.
97. R. F. W. Bader, Chem. Rev., 1991, 91, 893–928.
98. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
99. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
100. K. K. Irikura and D. J. Frurip, Computational Thermochemistry, American Chemical Society, Washington, DC, USA, 1998.
101. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. B. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A.Montgomery, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford, 2016.
102. J. Pfaendtner, X. Yu and L. J. Broadbelt, Theor. Chem. Acc., 2007, 118, 881–898.
103. V. Štejfa, M. Fulem and K. Růžička, J. Chem. Phys., 2019, 150, 224101.
104. L. Dmytro, N. Volovyk, O. Bevz, O. V. Vashchenko and O. Gryzodub, J. Pharm. Sci. Res., 2018, 10, 2709–2714.
105. J. Meija, T. B. Coplen, M. Berglund, W. A. Brand, P. De Bièvre, M. Gröning, N. E. Holden, J. Irrgeher, R. D. Loss, T. Walczyk and T. Prohaska, Pure Appl. Chem., 2016, 88, 265–291.
106. D. B. Newell, F. Cabiati, J. Fischer, K. Fujii, S. G. Karshenboim, H. S. Margolis, E. de Mirandés, P. J. Mohr, F. Nez, K. Pachucki, T. J. Quinn, B. N. Taylor, M. Wang, B. M. Wood and Z. Zhang, Metrologia, 2018, 55, L13–L16.
107. V. Štejfa, J. Rohlíček and C. Červinka, J. Chem. Thermodyn., 2021, 160, 106392.
108. V. Štejfa, O. Vojtíšková, V. Pokorný, J. Rohlíček, K. Růžička and M. Fulem, J. Therm. Anal. Calorim., 2024, 149, 6179–6193.
109. V. Pokorný, C. Červinka, V. Štejfa, J. Havlín, K. Růžička and M. Fulem, J. Chem. Eng. Data, 2020, 65, 1833–1849.
110. Y.-t Suzuki, Y. Yamamura, M. Sumita, S. Yasuzuka and K. Saito, J. Phys. Chem. B, 2009, 113, 10077–10080.
111. T. Degen, M. Sadki, E. Bron, U. König and G. Nénert, Powder Diffr., 2014, 29, S13–S18.
112. J. Faber, Crystallogr. Rev., 2004, 10, 97–107.
113. V. Štejfa, M. Fulem, K. Růžička and P. Morávek, J. Chem. Eng. Data, 2016, 61, 3627–3639.
114. V. Štejfa, S. Chun, V. Pokorný, M. Fulem and K. Růžička, J. Mol. Liq., 2020, 319, 114019.
115. V. Štejfa, M. Fulem and K. Růžička, J. Therm. Anal. Calorim., 2024, 149, 4709–4720.
116. T. Mahnel, V. Štejfa, M. Maryška, M. Fulem and K. Růžička, J. Chem. Thermodyn., 2019, 129, 61–72.
117. C. Tsonopoulos, AIChE J., 1974, 20, 263–272.
118. J. R. Elliott, V. Diky, T. A. I. V. Knotts and W. V. Wilding, Properties of Gases and Liquids, McGraw-Hill Education, New York, 2023.
119. E. R. Cox, J. Ind. Eng. Chem., 1936, 28, 613–616.
120. M. Kurahashi, M. Fukuyo, A. Shimada, A. Furusaki and I. Nitta, Bull. Chem. Soc. Jpn., 1969, 42, 2174–2179.
121. P. T. Clarke and J. M. Spink, Acta Crystallogr., Sect. B, 1969, 25, 162.
122. Y. Nambu, Y. Yoshitake, S. Yanagi, K. Mineyama, K. Tsurui, S. Kuwata, T. Takata, T. Nishikubo and K. Ishikawa, J. Mater. Chem. C, 2022, 10, 726–733.
123. S. Bhandary, R. Van Deun, A. M. Kaczmarek and K. Van Hecke, Chem. Sci., 2022, 13, 10308–10314.
124. J. McDowell, Acta Crystallogr., Sect. B, 1976, 32, 5–10.
125. A. Giunchi, L. Pandolfi, R. G. Della Valle, T. Salzillo, E. Venuti, N. Demitri, H. Riegler, C. Petschacher, J. Liu and O. Werzer, CrystEngComm, 2024, 26, 6573–6584.
126. H. Liu, Y. Gao, J. Cao, T. Li, Y. Wen, Y. Ge, L. Zhang, G. Pan, T. Zhou and B. Yang, Mater. Chem. Front., 2018, 2, 1853–1858.
127. S. B. Larson, S. H. Simonsen, G. E. Martin, K. Smith and S. Puig-Torres, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 1984, 40, 103–106.
128. H. Lynton and E. G. Cox, J. Chem. Soc., 1956, 0, 4886–4895.
129. R. D. Chirico, S. E. Knipmeyer, A. Nguyen and W. V. Steele, J. Chem. Thermodyn., 1991, 23, 431–450.
130. W. V. Steele, R. D. Chirico, S. E. Knipmeyer and A. Nguyen, J. Chem. Thermodyn., 1993, 25, 965–992.
131. P. M. Robinson and H. G. Scott, Mol. Cryst., 1969, 5, 387–404.
132. Z. Lisicki and M. E. Jamróz, J. Chem. Thermodyn., 2000, 32, 1335–1353.
133. K. Stark, P. Keil, S. Schug, K. Müller, P. Wasserscheid and W. Arlt, J. Chem. Eng. Data, 2016, 61, 1441–1448.
134. R. Sabbah and L. El Watik, Thermochim. Acta, 1989, 138, 241–247.
135. J. R. Donnelly, L. A. Drewes, R. L. Johnson, W. D. Munslow, K. K. Knapp and G. W. Sovocool, Thermochim. Acta, 1990, 167, 155–187.
136. R. K. Gupta, S. K. Singh and R. A. Singh, J. Cryst. Growth, 2007, 300, 415–420.
137. V. L. S. Freitas, J. R. B. Gomes and M. D. M. C. Ribeiro da Silva, J. Chem. Thermodyn., 2014, 73, 110–120.
138. A. Krajewska and K. Pigoń, Thermochim. Acta, 1980, 41, 187–197.
139. J. D. Bell, J. F. Blount, O. V. Briscoe and H. C. Freeman, Chem. Commun., 1968, 1656–1657.
140. D. Feil, M. H. Linck and J. J. H. McDowell, Nature, 1965, 207, 285–286.
141. H. Nakayama, K. Ishii, E. Chijiwa, M. Wada and A. Sawada, Solid State Commun., 1985, 55, 59–61.
142. B. W. van de Waal and D. Feil, Acta Crystallogr., Sect. B, 1977, 33, 314–315.
143. M. Radomska and R. Radomski, Thermochim. Acta, 1980, 40, 405–414.
144. V. L. S. Freitas, J. R. B. Gomes and M. D. M. C. Ribeiro da Silva, J. Chem. Thermodyn., 2009, 41, 1199–1205.
145. C. Červinka, M. Fulem and K. Růžička, J. Chem. Eng. Data, 2012, 57, 227–232.
146. V. Štejfa, M. Fulem and K. Růžička, J. Chem. Phys., 2019, 151, 144504.
147. V. Štejfa, M. Fulem and K. Růžička, J. Mol. Liq., 2023, 380, 121724.
148. V. Štejfa, T. Mahnel, E. Skořepová, J. Rohlíček, V. Eigner, B. Schröder, K. Růžička and M. Fulem, J. Chem. Thermodyn., 2020, 150, 106193.
149. V. Štejfa, V. Pokorný, J. Rohlíček, M. Fulem and K. Růžička, J. Chem. Thermodyn., 2021, 160, 106488.
150. K. Růžička, V. Štejfa, C. Červinka, M. Fulem and J. Šturala, Molecules, 2024, 29, 1110.
151. S. E. Stein and B. D. Barton, Thermochim. Acta, 1981, 44, 265–281.
152. C. E. Senseman and O. A. Nelson, Ind. Eng. Chem., 1923, 15, 382–383.
153. M. J. S. Monte, C. A. D. Sousa, J. M. S. Fonseca and L. M. N. B. F. Santos, J. Chem. Eng. Data, 2010, 55, 5264–5270.
154. S. C. Mraw and C. F. Keweshan, J. Chem. Thermodyn., 1984, 16, 873–883.
155. J. Řezáč, Y. Huang, P. Hobza and G. J. O. Beran, J. Chem. Theory Comput., 2015, 11, 3065–3079.

---

Footnote

† Both authors contributed equally to this work.

---