Glass transition temperature prediction of disordered molecular solids

Amorphous molecular solids—disordered organic molecular materials—are central to electronic and optoelectronic applications (e.g., OLEDs and perovskite solar cells) due to their processability and tunable properties. The glass transition temperature (T_g) governs morphological stability and device lifetimes, motivating the search for high-T_g compounds and reliable predictive tools. Previous approaches include QSPR and machine learning models, as well as direct extraction of T_g from molecular dynamics (MD) simulations. While tailor-made forcefields have yielded accurate T_g for specific systems, a deterministic, transferable workflow for diverse small conjugated molecules has been lacking, in part due to the high computational and manual cost of parameterization and MD, and difficulties in creating a forcefield protocol that transfers across diverse chemistries. The authors aim to develop a computational methodology that combines an objective fitting protocol for density–temperature data with an automated non-bonded parameterization (via atoms-in-molecule partitioning) to accurately and reproducibly predict T_g for a set of OLED host and hole-transport materials.

Prior work has pursued T_g prediction via statistical models (QSPR, machine learning) and direct MD simulations. QSPR and ML efforts (e.g., for conjugated polymers and OLED materials) can be effective but may lack transferability across chemistries or depend on training data quality. Direct MD extraction of T_g is conceptually appealing but often depends on tailor-made forcefields and suffers from analysis subjectivity (choosing fitting ranges in bilinear density–temperature fits). The authors’ earlier multiscale protocols accurately predicted electronic properties (ionization energies, DOS, mobility) for OLED hosts using OPLS parameters and MK charges; however, T_g predictions were inaccurate, suggesting shortcomings in non-bonded parameterization. Literature indicates that non-bonded interactions and dihedral correlations critically influence glass transition behavior, and that uncertainty in analysis (e.g., bilinear fitting choices) can be substantial. These gaps motivate a transferable non-bonded parameterization (DDEC6-based) and an objective fitting strategy that reduces human bias.

The workflow comprises three main components: (i) forcefield parameterization, (ii) classical MD simulations to generate density–temperature data, and (iii) an objective fitting protocol for T_g extraction. (i) Forcefield parameterization: Bonded terms (except proper/improper dihedrals) are taken from OPLS-AA and prior work. Non-bonded parameters (atomic partial charges and Lennard-Jones coefficients) are derived via an atoms-in-molecule, electron density-based protocol: electronic densities are partitioned using DDEC6 to obtain atomic charges, and Lennard-Jones A and B parameters are obtained via the Tkatchenko–Scheffler scheme using free-atom radii. Electronic structure calculations use Gaussian 16 at the wB97X-D3/6-311G(d,p) level; DDEC6 partitioning uses Chargemol (2017-09-26). Molecules are decomposed into rigid fragments. Missing torsional terms connecting rigid fragments are parameterized by constrained dihedral scans at the wB97X-D3/def2-TZVP level using ORCA 4.2.1. (ii) Classical MD simulations: Simulations are performed with GROMACS 2020.3. Electrostatics use PME with 0.12 nm Fourier spacing; non-bonded cutoff 13 Å. Thermostat: velocity rescaling with stochastic term (τ_T = 0.5 ps). Barostat: isotropic Berendsen (P0 = 1 bar; τ_P = 0.5 ps; compressibility parameterization as specified). Systems contain 3000 molecules initially packed at low density (50–150 kg m^-3) using Packmol. Thermal protocol: heat from 100 K to 300 K at 0.67 K ps^-1, equilibrate at 300 K to stabilize density; then heat to 800 K at 0.5 K ps^-1 and equilibrate 10 ns at 800 K; finally, linearly cool from 800 K to 0 K at 100 K ns^-1 while recording density and temperature for T_g analysis. Stepwise vs continuous cooling and other protocol variations (system size, initial configurations, barostats, cooling rates, density fluctuations) were tested; main results use continuous cooling. (iii) Fitting protocol: Traditional T_g extraction uses a bilinear fit to the density–temperature (ρ–T) curve, but non-ideal behavior and slope variability introduce subjective choice of fitting ranges. The authors propose plotting R^2 versus temperature for sliding-window linear fits over [T, T+α), where α is the fitting window (tested α = 50–300 K). The R^2–T plot typically exhibits a valley separating two linear regimes (glassy and liquid-like). The two fitting ranges are chosen automatically as the two local maxima (“hill tops”) of R^2 on either side of the valley, eliminating manual bias. Sensitivity to α is assessed; α = 200 K is selected for final comparisons due to overall best agreement with experiments across the dataset. Experimental T_g measurement (for validation): Differential scanning calorimetry (DSC) at Merck KGaA using 10–15 mg samples in a Netsch DSC 204/1/G Phönix: heat 5 °C min^-1 to 370 °C, cool 20 °C min^-1 to 0 °C, reheat 20 °C min^-1 to 370 °C; T_g taken at half-step in heat flow. For TMBT (and BCP where standard protocol lacked a clear kink), alternative protocol on TA Instruments DSC Discovery in N2 with 5 mg sample, heat 20 °C min^-1 to 320 °C, quench with liquid N2, reheat 20 °C min^-1 to 320 °C to observe T_g.

- Bilinear fit ambiguity and human bias: For BCP, varying fitting ranges yields T_g from 51.2 °C to 147.7 °C, evidencing large analysis-induced variability. - R^2–T fitting strategy: The R^2 valley robustly identifies transition region across compounds; selecting the R^2 maxima on either side reduces bias. Dependence on fitting window α is modest overall; α = 150–200 K performs best, with α = 200 K chosen. - Forcefield impact: Using OPLS with MK charges (OPLS-MK) yields poor agreement with experiment (MAE = 59.1 °C; R^2 = 0.61; slope ≈ 1.01). Reparameterizing non-bonded terms via DDEC6 + TS dramatically improves predictions (MAE = 20.5 °C; R^2 = 0.87; slope ≈ 1.25) over 24 compounds (excluding outliers MTDATA and 2-TNATA). Experimental T_g span ≈ 50–230 °C. - Comparison to Patrone et al. hyperbolic fit: On the same MD data (DDEC6), Patrone’s method gives MAE = 64.7 °C; R^2 = 0.75; slope ≈ 0.97, whereas the proposed R^2-based bilinear procedure achieves MAE = 20.5 °C and R^2 = 0.87, indicating superior accuracy and correlation. - Outliers and mechanistic insight: MTDATA and 2-TNATA (starburst TPA-based dendrimers) are outliers for both OPLS-MK and DDEC6, likely due to missing non-local dihedral correlations in TPA moieties; improved or ML-augmented forcefields are suggested. - Computational efficiency: Total ∼846,029 CPU hours across compounds; distribution: ~50.1% MD compressing/heating, ~29.9% MD cooling, ~19.9% QM constrained dihedral scans, ~0.1% DDEC charge/LJ parameter derivation. Human effort ≈ 6 h per compound, ~70% spent on topology and parameter generation. - Protocol robustness: Cooling rate reduction to 10 K ns^-1 did not eliminate non-idealities; chosen 100 K ns^-1 protocol with objective fitting provides practical accuracy for screening.

The study addresses the core challenge of reliably predicting T_g of disordered molecular solids by tackling two sources of error: (1) analysis bias in extracting T_g from non-ideal ρ–T data and (2) insufficient transferability of non-bonded forcefield parameters. The R^2–T-based fitting protocol objectively determines linear regimes, substantially reducing human-induced variability. Reparameterizing non-bonded interactions via DDEC6/TS provides physically grounded partial charges and van der Waals parameters, improving agreement with experiments across diverse OLED hosts and transport materials. Together, these advances reduce the MAE to ~20 °C and raise correlation to R^2 ≈ 0.87 over a wide T_g range, enabling practical in silico prescreening. Remaining discrepancies in TPA-based starburst dendrimers highlight the importance of capturing correlated dihedral behavior; integrating such correlations (e.g., via enhanced functional forms or ML) is a logical next step. Despite significant computational cost, the workflow is automatable and the open-access OPLS-DDEC parameters can facilitate broader multiscale studies (e.g., guest–host and interfacial systems).

The authors present a reliable and automated methodology for T_g prediction in amorphous organic semiconductors by combining an objective R^2–T fitting protocol with DDEC6-based non-bonded forcefield parameterization. Across 26 OLED-relevant compounds (24 used for metrics), the approach achieves MAE ≈ 20.5 °C and R^2 ≈ 0.87 versus experiment, a substantial improvement over OPLS-MK parameters and alternative fitting schemes. The method reduces analysis bias, offers transferable non-bonded parameters, and is suitable for prescreening in materials design workflows. Future work should focus on improving forcefields for systems with strongly correlated torsions (e.g., TPA-based starbursts), optimizing MD protocols to lower computational cost (e.g., smaller systems, denser initial morphologies), and extending applicability to more complex device-relevant morphologies and chemistries.

- Two clear outliers (MTDATA and 2-TNATA) indicate limitations in capturing non-local dihedral correlations in TPA-based starburst structures with current forcefields. - Non-ideal ρ–T curves and density fluctuations can complicate fitting; while the R^2-driven approach reduces bias, dependence on fitting window size remains. - High cooling rates (100 K ns^-1) used for practicality may not fully represent experimental conditions; slower rates are computationally prohibitive for large systems. - Computational cost is substantial (≈8.5×10^5 CPU hours total) and human effort in parameterization remains nontrivial, suggesting a need for further automation and efficiency improvements. - Transferability, while improved, may still be challenged by chemistries outside the tested fragment set.