Paper deep dive
MultiPUFFIN: A Multimodal Domain-Constrained Foundation Model for Molecular Property Prediction of Small Molecules
Idelfonso B. R. Nogueira, Carine M. Rebello, Mumin Enis Leblebici, Erick Giovani Sperandio Nascimento
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/20/2026, 4:08:01 AM
Summary
MultiPUFFIN is a multimodal domain-constrained foundation model for predicting thermophysical properties of small molecules. It fuses SMILES sequences, 2D molecular graphs, and 3D conformer geometries using bidirectional cross-modal attention and gated fusion. The model is pretrained on 500,000 unlabeled PubChem molecules and employs a condition-aware refinement stack and a four-head tournament for each of nine properties to enforce thermodynamic consistency. It outperforms ChemBERTa-2 with significantly fewer labeled molecules.
Entities (13)
Relation Signals (12)
MultiPUFFIN â predicts â Viscosity
confidence 95% · predicting thermophysical properties... viscosity
MultiPUFFIN â predicts â Vapor Pressure
confidence 95% · predicting thermophysical properties... vapor pressure
MultiPUFFIN â pretrainedon â PubChem
confidence 95% · backbone is pretrained on 500,000 unlabelled PubChem molecules
SchNet â processes â 3D Conformer
confidence 95% · 3D conformer is processed by a SchNet encoder
GCN â processes â 2D Molecular Graph
confidence 95% · 2D molecular graph is processed by a GCN encoder
Transformer â processes â SMILES
confidence 95% · SMILES string is processed by a Transformer encoder
MultiPUFFIN â usesmodality â SMILES
confidence 95% · MultiPUFFIN fuses SMILES sequences, 2D molecular graphs, and 3D conformer geometries
MultiPUFFIN â usesmodality â 2D Molecular Graph
confidence 95% · MultiPUFFIN fuses SMILES sequences, 2D molecular graphs, and 3D conformer geometries
MultiPUFFIN â usesmodality â
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:MultiPUFFIN is a domain-informed multimodal foundation model for predicting thermophysical properties of small molecules, addressing a critical gap in chemical engineering, drug discovery, and materials science. Existing molecular foundation models pretrain on millions of molecules to learn general-purpose representations, but their standard MLP output layers impose no physical constraints, vapor pressure predictions may violate monotonic temperature dependence, and viscosity curves may lack the functional form required by process simulators. Domain-informed approaches that guarantee thermodynamic consistency have remained limited to single properties and small datasets, whereas multimodal foundation models have focused on biological activity rather than thermophysical properties. MultiPUFFIN fills this gap by fusing SMILES sequences, 2D molecular graphs, and 3D conformer geometries through bidirectional cross-modal attention and gated fusion, supplemented by auxiliary encoders for experimental conditions and molecular descriptors. The backbone is pretrained on 500,000 unlabelled PubChem molecules using three complementary self-supervised objectives. A condition-aware refinement stack of five conditioners (temperature, pH, pressure, polymorph, and measurement method) routes each property to a four-head tournament that selects the best-performing thermodynamically informed head for that property. MultiPUFFIN achieves a mean test R2 of 0.784 and outperforms fine-tuned ChemBERTa-2 on all nine properties despite training on roughly 2,000x fewer labeled molecules.
Tags
Links
- Source: https://arxiv.org/abs/2603.00857v2
- Canonical: https://arxiv.org/abs/2603.00857v2
Trouble viewing inline? Open PDF directly â
Full Text
197,856 characters extracted from source content.
Expand or collapse full text
MultiPUFFIN: A Multimodal Domain-Informed Foundation Model for Molecular Property Prediction of Small Molecules Idelfonso B.R. Nogueira idelfonso.nogueira@ntnu.no Carine M. Rebello Mumin Enis Leblebici Erick Giovani Sperandio Nascimento Abstract Predicting the physicochemical properties of small molecules across a diverse chemical space is an important step for chemical engineering, drug discovery and materials science. Existing molecular foundation models pretrain on millions of molecules to learn general-purpose representations, but their standard MLP output layers impose no physical constraints: vapour pressure predictions may violate monotonic temperature dependence, and viscosity curves may lack the functional form required by process simulators. Domain-informed approaches that do guarantee thermodynamic consistency have, by contrast, been limited to single properties and small datasets, and the existing multimodal molecular foundation models are typically aimed at biological-activity prediction rather than thermophysical properties. This work introduces MultiPUFFIN, a domain-informed multimodal foundation model that fills this gap. MultiPUFFIN fuses SMILES sequences, 2D molecular graphs and three-dimensional conformer geometries through bidirectional cross-modal attention and gated fusion, supplemented by auxiliary encoders for experimental conditions and molecular descriptors. The backbone is pretrained on 500,000 unlabelled PubChem molecules with three complementary self-supervised objectives, and a condition-aware refinement stack of five conditioners (T, pH, P, polymorph, method) routes per property to a four-head tournament that selects the best-performing thermodynamically-informed head per property. Two cross-property training-loss couplings (flash-pointâvapor-pressure consistency and a normal-boiling-point anchor on the vapor-pressure head) softly enforce inter-property thermodynamic consistency. On a scaffold-split test set of 88778877 molecules drawn from a multi-source dataset of 37 96837\,968 unique molecules across nine thermophysical properties, MultiPUFFIN reaches an in-scope mean test R2=0.784R^2=0.784 and outperforms fine tuned ChemBERTa-2 on all nine properties despite training on roughly 2000Ă fewer labelled molecules; the largest margins are on the temperature-dependent properties (vapor pressure, viscosity) where ChemBERTa-2 has no mechanism to distinguish measurements at different temperatures. Per-property G+ applicability-domain filters explicitly bound the deployment scope. Systematic ablations quantify the contribution of each architectural component, of the choice of domain-informed equation per property, and of each coupling term. keywords: Foundation model , Molecular property prediction , Graph neural network , Inductive bias , Domain-informed machine learning , Multimodal learning , Multi-task learning , Physics-informed prediction , Thermophysical properties â journal: Preprint [ntnu]organization=Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), city=Trondheim, postcode=7034, country=Norway [kuleuven]organization=Faculty of Industrial Engineering, KU Leuven, city=Diepenbeek, country=Belgium [surrey]organization=University of Surrey, city=Guildford, country=United Kingdom highlights MultiPUFFIN fuses SMILES, molecular graphs, and 3D conformers with domain-informed prediction heads for nine thermophysical properties simultaneously. Domain-informed inductive biases (Antoine, Andrade, van ât Hoff, Born, Shomate equations) enforce physically consistent temperature dependence within individual properties by construction, and are competed against fragment and direct baselines in a per-property head tournament. A self-supervised pretraining stage on 500,000 unlabeled PubChem molecules with three complementary objectives (graphâSMILES contrastive InfoNCE, masked atom feature reconstruction, masked SMILES token prediction) and two cross-property physical coupling terms (flash-pointâvapor-pressure consistency and a stronger normal-boiling-point anchor on the vapor-pressure head) elevate the architecture to a thermodynamically coupled multimodal foundation model. A four-stage training protocol (SSL pretrain, joint supervised multi-task training with the condition module stack, backbone-unfrozen targeted fine tune at very low learning rate, per-property G+ applicability-domain evaluation) yields a single deployable artifact whose accuracy and scope are reported jointly. MultiPUFFIN achieves higher test R2R^2 than fine tuned ChemBERTa-2 on all nine properties despite using a single multi-task model vs. nine separately fine tuned models, with the most pronounced advantages on temperature-dependent properties (vapor pressure, viscosity) where ChemBERTa-2 cannot distinguish measurements of the same molecule at different temperatures. Systematic ablation across a set of candidate thermodynamic representations per property reveals that the optimal domain-informed equation is property-specific, with group contribution heads (themselves rooted in thermodynamic additivity) providing the best inductive bias for several properties. Inductive biases and multimodal encoding substantially reduce data and computational requirements compared to brute-force pretraining at scale. 1 Introduction The prediction of physicochemical properties of small molecules underlies the design of separation processes, the formulation of pharmaceuticals, environmental fate modelling, and reaction engineering [44]. Properties such as vapor pressure, viscosity, solubility, heat capacity and partition coefficients are conventionally estimated through empirical correlations whose parameters are component-specific: the Antoine equation for vapor pressure [56], the Andrade equation for viscosity [3] and group contribution methods for boiling points [43] all require dedicated fitted parameters per molecule, which restricts transferability and precludes large-scale screening. Machine learning can help with this limitation by learning structureâproperty mappings directly from data [12]. Graph neural networks treat atoms as nodes and bonds as edges, and have produced strong performance on quantum mechanical properties [51, 18], solvation free energies [59], fuel ignition quality [52] and activity coefficients [49]; chemical language models operating on SMILES strings offer competitive performance on property prediction tasks as well [12]. The landscape of thermophysical property prediction in this body of work is governed by single-model, single-prediction architectures, occasionally extended with domain-informed prediction heads to recover physical consistency for one property at a time [50, 47]. Foundation models extend the shallow learners discussed above by pretraining a single backbone on a large unlabelled corpus and adapting it through lightweight prediction heads [10, 5]. SMILES-BERT [60], ChemBERTa-2 [2], KPGT [29], MolCLR [61], MolE [41] and the Uni-Mol family [65, 38, 57] occupy this niche. The landscape of foundation models for molecules, despite this growth, is dominated by single-modal models, and the few multimodal members of the family target biological-activity prediction rather than thermophysical properties: MoleculeSTM [35], MoMu [54] and GIT-Mol [34] align structural representations with natural language, while MolGT [8], FineMolTex [33] and MolPrompt [32] refine the alignment granularity. None of them combines multiple structural modalities with domain-informed prediction heads, and none of them does multi-task thermophysical property prediction. Table 1 summarises this landscape: each row marks where the published architecture stands on the three orthogonal axes of multimodality, domain-informed heads and multi-task thermophysical coverage; the bottom row is the position MultiPUFFIN fills. The single-modal limitation matters because each molecular representation carries a different slice of the underlying chemistry rather than the totality of it. The 2D molecular graph encodes topological connectivity and local functional-group environment, the SMILES string encodes long-range syntactic dependencies and substituent patterns [62, 4], and the 3D conformer encodes steric effects, solvation cavity geometry and through-space distances; no single representation contains all three. Any single-modal foundation model is therefore representationally bounded by the information present in the chosen modality, and the natural way to extend its representational reach is multimodal fusion. Multimodality alone is not sufficient, however. The architectures cited above show that combining modalities does help on biological-activity tasks, but a property prediction trained against pure data, with no thermodynamic constraint baked into the output layer, may still violate elementary physics: vapor pressure can decrease with temperature, viscosity can increase with temperature, and Antoine coefficients can come out with unphysical signs. A multimodal architecture is one way to fuse domain expertise into the model inputs, and the gated cross-modal fusion of this work additionally exposes the implicit causality between modalities (Section 2.10). A complementary way is to fuse domain expertise into the model output or training objective: PUFFIN [50] and ExPUFFIN [47] introduced domain-informed inductive-bias neurons that replace the standard output layer with a thermophysical equation whose coefficients are predicted by the networkâs penultimate layer; the network is never supervised on the coefficients, it learns them from property targets alone, and post-hoc inspection shows the recovered coefficients agree with the literature. The ExPUFFIN Andrade variant reduced viscosity RMSE by 37% over the unconstrained baseline while delivering smooth monotonic viscosityâtemperature curves. Both PUFFIN and ExPUFFIN are single-property, single-modality systems; the generalisation to a multi-property, multimodal foundation model with domain-informed heads has not been attempted. A multimodal, multi-task and domain-informed model is exactly what molecular screening pipelines need at the pre-experimental stage. A medicinal chemist evaluating a drug candidate, a process engineer specifying an ester-based heat-transfer fluid, or a formulator selecting a solvent all need simultaneous, thermodynamically consistent estimates of multiple properties for the same molecule. The remainder of the paper presents MultiPUFFIN (Multimodal Path-Unifying Foundation Fusion Interfaced Network), a domain-informed multimodal foundation model that fills the empty bottom row of Table 1 and supports this screening role. Table 1: Comparison of molecular foundation models for property prediction. Modalities: S = SMILES (1D), G = 2D graph, 3D = 3D conformer, T = text, I = molecular image, Aux = auxiliary inputs (experimental conditions and molecular descriptors). Domain heads: whether the model employs thermodynamically-informed or domain-informed output layers (e.g., Antoine, Andrade, group contribution equations). Phys. consist.: whether physically consistent temperature dependence is enforced within individual properties by construction (note: this does not imply full cross-property thermodynamic consistency). Multi-task thermo.: whether the model simultaneously predicts multiple thermophysical properties in a single architecture. Model Year Modalities Pretraining Domain heads Multi-task thermo. Phys. consist. Pretrain data size SMILES-BERT [60] 2019 S MLM Ă Ă Ă ZINC subset ChemBERTa-2 [2] 2022 S MLM + MTR Ă Ă Ă 77M mol. KPGT [29] 2022 G SSL + knowledge Ă Ă Ă 2M mol. MoMu [54] 2022 G, T Contrastive Ă Ă Ă 15K pairs Uni-Mol [65] 2023 3D SSL (3D denoising) Ă Ă Ă 209M conf. MoleculeSTM [35] 2023 S, G, T Contrastive Ă Ă Ă 281K pairs MolE [41] 2024 G SSL + multi-task Ă Ă Ă 842M graphs GIT-Mol [34] 2024 G, I, T Multimodal LLM Ă Ă Ă 304K mol. MolGT [8] 2024 G, T SSL multi-view Ă Ă Ă 250K mol. Uni-Mol2 [38] 2024 3D SSL (scaled) Ă Ă Ă 800M conf. MoleculeFormer [30] 2025 S, G Supervised MT Ă Ă Ă 28 datasets PUFFIN [50] 2024 G Transfer learn. â Ă â 6K mol. ExPUFFIN [47] 2025 G Supervised â Ă â 3K mol. MultiPUFFIN (this work) 2025 S, G, 3D, Aux Multi-task â â â 38K mol. The contributions of this work are: (i) a multimodal backbone that fuses a Transformer over SMILES, a GCN over the 2D graph, and SchNet over the 3D conformer through bidirectional cross-modal attention and gated fusion, with auxiliary encoders for experimental conditions and molecular descriptors, so that each modality contributes its own slice of information to the unified embedding; (i) a generalised domain-informed inductive-bias mechanism that extends the PUFFIN / ExPUFFIN single-property paradigm to nine simultaneous thermophysical targets through a per-property four-head tournament (primary thermophysical equation, Joback group contribution, RDKit fragment counts, and an alternative thermodynamically-informed or direct FFNN baseline), enforcing intra-property physical consistency by construction; (i) a four-stage training strategy combining SSL pretraining on 500,000 unlabelled PubChem molecules, joint supervised multi-task training with a condition-aware refinement stack of five conditioners (T, pH, P, polymorph, method), backbone-unfrozen targeted fine tuning at very low learning rate, and per-property G+ applicability-domain evaluation that bounds the deployment scope. On a scaffold-split test set of 8,877 molecules drawn from a 37,968-molecule multi-source dataset (eleven public sources: OPERA, NIST ThermoML, ECHA REACH, ChEMBL, AqSolDB, FreeSolv, Bradley melting points, Sun et al. flash points, ABB-ADD heat capacities, the Chew et al. viscosity compilation, and PubChem), MultiPUFFIN reaches an in-scope mean test R2=0.784R^2=0.784 and beats fine tuned ChemBERTa-2 on all nine properties despite training on roughly 2000Ă fewer labelled molecules. The remainder of this paper is organised as follows. Section 2 presents the methodology, including data curation, model architecture and training. Section 3 presents the experimental results and discussion. Section 4 concludes. 2 Methodology The overall MultiPUFFIN framework encompasses five components: a data curation and preprocessing pipeline, a self-supervised pretraining pipeline that produces a chemistry-aware backbone, the multimodal model architecture with a condition-aware embedding refinement stack (TConditioner module among five) and domain-informed prediction heads with cross-property physical coupling, a four-stage training strategy that culminates in a backbone-unfrozen targeted fine tune at very low learning rate, and a per-property applicability-domain (G+) evaluation that explicitly bounds the deployment regime. This section presents each in detail. MultiPUFFIN is designed as a multimodal foundation model for thermophysical property prediction. âFoundation modelâ is used in the now-standard sense established in visionâlanguage and chemistry pretraining literature: a single model whose backbone is first pretrained on a large unlabeled molecular corpus through self-supervised objectives that do not require property labels, and whose representations are then fine tuned for downstream property regression through a small set of supervised signals. In MultiPUFFIN the backbone is multimodal: it jointly processes three fundamentally different structural data modalities (text, graph, and spatial structure) through modality-specific encoders, supplemented by two auxiliary encoders for experimental conditions and molecular descriptors, before fusing all representations into a shared molecular embedding. This design mirrors the broader paradigm of multimodal foundation models in machine learning (e.g., CLIP for visionâlanguage, Gato for visionâlanguageâaction), where the central premise is that jointly learning from complementary data modalities produces richer representations than any single modality alone. In the molecular context, the three structural modalities correspond to: (i) the text modality, where a SMILES string is processed by a Transformer encoder as a chemical language sequence, capturing long-range syntactic dependencies and implicit chemical grammar; (i) the graph modality, where the 2D molecular graph is processed by a GCN encoder through message-passing operations, capturing topological connectivity, ring systems, and local functional group patterns; and (i) the spatial modality, where a 3D conformer is processed by a SchNet encoder through continuous-filter convolutions on interatomic distances, capturing through-space geometry, molecular shape, and steric effects. In addition, two auxiliary encoders process non-structural information: an experimental encoder that embeds thermodynamic conditions (temperature, pressure) into the representation space, and a descriptor encoder that incorporates precomputed molecular descriptors (molecular weight, topological polar surface area, hydrogen bond donors/acceptors, and other constitutional features). Three additional components distinguish MultiPUFFIN from prior PUFFIN-family models and from generic multimodal foundation models for molecules. First, the backbone is pretrained with three complementary self-supervised objectives on an unlabeled molecular corpus before any property labels are used, so that the fused 512-dimensional embedding already encodes chemistry-aware features by the time the supervised multi-task training begins. Second, between the fused embedding and the prediction heads sits a stack of five identity-initialized condition-aware refinement modules (Section 2.9): a TConditioner for measurement temperature, a PHConditioner for solution pH, a PressureConditioner for ambient pressure, a PolymorphEmbedding for crystal-form identifier, and a MethodFlag for test-protocol identifier. The modules are uniformly architectured as 2-token cross-attention blocks with zero-init residuals, route per property to the heads where the corresponding condition is physically relevant, and compose by chaining when multiple conditions matter (e.g., âTâTâPu _T _Tâ P for viscosity and heat capacity). The stack lets the network learn arbitrary condition-dependent embedding modulations beyond what a fixed Antoine, Andrade, or Shomate functional form can express, and the identity-init property guarantees that any module that sees no condition variance in its training subset gracefully degrades to the unconditioned embedding. Third, the prediction heads are organized as a tournament of candidate inductive biases per property, and the training loss carries explicit cross-property physical coupling terms (flash-pointâvapor-pressure consistency and a stronger normal-boiling-point anchor on the vapor-pressure head) that enforce thermodynamic relationships no purely data-driven foundation model can provide. The deployed model is a single forward pass of the four-component pipeline (SSL-pretrained backbone â condition-aware refinement stack with per-property routing â tournament-winning head per property â G+ scope check); no model-level ensembling is used at inference. Figure 1 provides a schematic overview of the architecture. Bidirectional cross-modal attention between the GCN and Transformer branches enables each modality to attend to the other, and a learned sigmoid gate fuses the cross-attended representations. The SchNet embedding is incorporated through a separate geometry gate that can suppress the 3D contribution when conformer data is unavailable or unreliable. The outputs of the two auxiliary encoders (experimental conditions and molecular descriptors) are then concatenated with the fused structural embedding and projected to yield the final 512-dimensional unified embedding u, which serves as the shared molecular representation from which all nine property-specific tournament heads operate. Each tournament comprises four candidate inductive biases (Section 2.11): a primary thermophysical equation, a Joback group-contribution head, an RDKit fragment-count head, and an alternative thermodynamically-informed or direct feed-forward head; the winning head per property is selected on the validation set. An uncertainty-weighted multi-task loss automatically balances the nine prediction objectives during training, and two cross-property physical coupling terms (Section 2.12) enforce thermodynamic consistency between flash point and vapor pressure, and between boiling point and vapor pressure. Figure 1: Architecture overview of MultiPUFFIN. Three structural encoders (GCN over the 2D graph, Transformer over the SMILES string, SchNet over the 3D conformer) and two auxiliary encoders (experimental conditions; molecular descriptors) form the multimodal backbone, pretrained on 500,000 unlabeled PubChem molecules with three SSL objectives. Cross-modal attention plus gated fusion produces a 512-dimensional unified embedding, which feeds a stack of five identity-initialized condition-aware refinement modules (T, pH, P, polymorph, method); each module is routed only to heads where the corresponding condition is physically relevant. Each property is then predicted by the winner of a four-head tournament (primary thermophysical equation, Joback, RDKit, alternative thermodynamically-informed/direct FFNN). Stars (â ) mark enhanced-capacity heads. Cross-property coupling (flash-pointâvapor-pressure consistency; normal-boiling-point anchor) links the heads thermodynamically. The deployed model is a single forward pass; no model-level ensembling is used at inference. 2.1 Data curation and dataset construction A key component of this work is the assembly of a multi-property molecular dataset from diverse public sources. The final dataset comprises 37 96837\,968 unique molecules with measurements across nine physicochemical properties, totaling 40 90440\,904 data rows, since molecules with temperature-dependent properties contribute multiple rows at different temperatures. This dataset was assembled through a systematic multi-stage curation pipeline that integrates eleven public sources (OPERA, NIST ThermoML, ECHA REACH, ChEMBL, AqSolDB, FreeSolv, Bradley melting points, Sun et al. flash points, ABB-ADD heat capacities, the Chew et al. viscosity compilation, and PubChem), as summarized in Table 2. Figure 6 illustrates the highly heterogeneous data availability across properties and splits. The dataset was constructed by merging an initial curated baseline with multiple external databases. The primary sources and their contributions are as follows. AqSolDB and ESOL [53, 63]. Aqueous solubility (logâĄS S) data were aggregated from the AqSolDB database and the ESOL dataset from MoleculeNet, providing solubility values in logâĄ(mol/L) (mol/L) at 298.15 K298.15\,K, with temperature-dependent solubility data from PubChem included where available. After deduplication by canonical SMILES (retaining the median value for duplicate entries), the solubility dataset comprises 14 08514\,085 unique molecules. MoleculeNet Lipophilicity and ChEMBL [63, 19]. Octanol-water partition coefficients (logâĄP P) from the Lipophilicity dataset in MoleculeNet and PubChem were merged with ChEMBL lipophilicity measurements into 11 05011\,050 unique molecules. It should be noted that the ChEMBL data reports logâĄD D at pH 7.4, which is used here as a logâĄP P proxy. For non-ionizable molecules, logâĄDâlogâĄP Dâ P; for ionizable species, logâĄD D at physiological pH can differ from logâĄP P due to the fraction of the molecule in its ionized form. This approximation is standard in the cheminformatics community [63] and is justified by the fact that the majority of drug-like molecules in these datasets are predominantly un-ionized at pH 7.4. Future work could improve the treatment of ionizable molecules by explicitly modeling the ionization equilibrium or by training on curated logâĄP P-only datasets. FreeSolv [42]. Hydration free energies (ÎâGhyd G_hyd) were taken from the FreeSolv database, providing measurements in k/âââk /mol at 298.15 K298.15\,K for 642642 unique molecules. OPERA 2.9 [39]. The OPERA (OPEn structure-activity/property Relationship App) models database contributed computed and validated measurements across multiple properties: logâĄP P (41914191 compounds), melting point (13 85613\,856 compounds, converted from °C to KK), boiling point (24262426 compounds), vapor pressure (39453945 compounds, converted from log10âĄ(mmHg) _10(mmHg) to PaPa), and water solubility (46744674 compounds in log10âĄ(mol/L) _10(mol/L)). NIST ThermoML Archive [17]. The ThermoML archive, comprising 11 92311\,923 structured data files from the NIST/TRC source data system, was systematically parsed to extract experimental thermophysical measurements. This yielded boiling point data for 220 molecules, melting point for 13631363 molecules, vapor pressure for 17971797 molecules (converted from kPakPa to PaPa), viscosity for 527 molecules (converted from Pa sPa\,s to mPa smPa\,s), and heat capacity for 744 molecules (in J molâ1 Kâ1J\,mol^-1\,K^-1). A dedicated parsing script was developed to extract property values, temperatures, and compound identifiers from the ThermoML XML schema, with InChI-to-SMILES conversion via RDKit for compound identification. ECHA REACH [15]. The European Chemicals Agency (ECHA) Registration, Evaluation, Authorisation and Restriction of Chemicals (REACH) database was parsed for experimentally registered thermophysical data, contributing measurements of melting point, boiling point, vapor pressure, solubility, logâĄP P, flash point, and viscosity for industrial chemicals. Bradley Open Melting Point Dataset [6]. Melting point data (TmT_m) for approximately 30413041 compounds were obtained, with values converted from °C to KK and filtered to the physically reasonable range of 50 K to 1000 K50\,K1000\,K. Sun et al. Flash Point Dataset [55]. Flash point data (TfT_f) for approximately 14 69614\,696 entries with SMILES and flash points in Kelvin, filtered to 100 K to 1000 K100\,K1000\,K. ABB-ADD Heat Capacity Dataset [1]. Liquid heat capacity (CpC_p) data with multi-temperature measurements for approximately 968 compounds identified by InChI strings; InChI identifiers were converted to canonical SMILES via RDKit, and for each molecule the measurement closest to 298.15 K298.15\,K was selected. Chew et al. Viscosity Dataset [9]. An additional viscosity compilation of approximately 10051005 compounds from the supplementary material of a recent Journal of Cheminformatics publication was integrated to enrich the viscosity training data. PubChem [27]. The PubChem database served as the primary aggregation source for boiling point, vapor pressure, and viscosity data, with experimental measurements spanning different temperatures. Boiling point data were converted to Kelvin and filtered to 100 K to 1000 K100\,K1000\,K; vapor pressure data are stored as log10âĄ(P/Pa) _10(P/Pa); viscosity data are stored as log10âĄ(η/mPaâ s) _10(η/mPa·s), preserving multi-temperature measurements per molecule. Table 2: Summary of data sources and per-property coverage in the final curated dataset. Unique molecules indicates the total number of distinct chemical compounds with measurements for each property across all sources. The rightmost columns show the train/validation/test split sizes (unique molecules per property per split). The Stage-3 augmentation step (Section 2.13) extends the per-property training rows but does not change the unique-molecule totals or test-split assignment reported here; the per-property test counts used in the deployment evaluation tables (Tables 7 and 13) are the post-Stage-3 data-row counts and therefore differ from the unique-molecule counts in this table for the multi-temperature properties. Property Primary sources Total mol. Train Val Test Units Solubility (logâĄS S) AqSolDB, ESOL, OPERA, PubChem, ECHA 14 08514\,085 86068606 28232823 26562656 logâĄ(mol/L) (mol/L) logâĄP P MoleculeNet, ChEMBL, OPERA, ECHA 11 05011\,050 83298329 13601360 13611361 â Hydration free energy FreeSolv 642 542 50 50 k/âââk /mol Boiling point (TbT_b) PubChem, NIST, OPERA, ThermoML, ECHA 85498549 37763776 29252925 18481848 KK Vapor pressure PubChem, NIST, OPERA, ThermoML, ECHA 89568956 50905090 22442244 16221622 log10âĄ(Pa) _10(Pa) Viscosity PubChem, Chew et al., ThermoML, ECHA 19011901 17521752 78 71 log10âĄ(mPaâ s) _10(mPa·s) Melting point (TmT_m) Bradley, OPERA, ThermoML, ECHA 18 91518\,915 10 41210\,412 24912491 60126012 KK Flash point (TfT_f) Sun et al., ECHA 10 26010\,260 45084508 31713171 25812581 KK Heat capacity (CpC_p) ABB-ADD, ThermoML 15271527 13801380 78 69 J molâ1 Kâ1J\,mol^-1\,K^-1 Total (unique molecules) 37 96837\,968 22 01522\,015 70767076 88778877 Total (data rows) 40 90440\,904 24 51324\,513 72547254 91379137 All molecular structures are represented as canonical SMILES strings generated by RDKit [46]. A uniform quality control protocol was applied across all sources: each SMILES string is parsed and re-canonicalized, with entries that fail RDKit parsing discarded; multi-component molecules (salts, mixtures, identified by the presence of â.â in the SMILES) are removed; molecules are filtered to contain between 2 and 100 heavy atoms; for entries sharing the same canonical SMILES and temperature, the median property value is retained; and property values outside physically reasonable ranges (e.g., boiling points below 100 K100\,K or above 1000 K1000\,K, negative viscosities, hydration free energies outside [â30,+10][-30,+10] k/âââk /mol) are removed. For temperature-independent properties (logâĄP P, TbT_b, ÎâGhyd G_hyd, TmT_m, TfT_f, CpC_p), values are associated with the measurement temperature or a default of 298.15 K298.15\,K, while for temperature-dependent properties (solubility, vapor pressure, viscosity), each temperatureâproperty pair constitutes an independent training sample. Several temperature-independent properties are propagated across all rows of a given molecule. For example, if a molecule has vapor pressure measurements at ten different temperatures and a single boiling point value, the boiling point label is copied to all ten rows. This maximizes the multi-task learning signal per training example. An evaluation mask is maintained for each temperature-independent property to prevent artificial inflation of evaluation metrics: only one row per molecule (the one closest to 298.15 K298.15\,K) is flagged for evaluation. A non-trivial design question for any multi-source dataset is whether the test distribution is genuinely covered by the training distribution within each property. The headline split-level statistics (ntrain=22,015n_train=22,015, ntest=8,877n_test=8,877) hide structural imbalances at the per-(property Ă molecule-class) level: a property may be well-supported overall but have a single chemical class (e.g., salts, polyols, sulfonic acids) that is heavily represented in the test set without comparable training coverage. To make these imbalances visible, we built a per-(property, class) coverage matrix (Table 3). The classes correspond to the seven chemistry exclusions encoded in the per-property G+ applicability-domain filters of Section 2.14 (salt, charged ion, polyol, sulfonic/phosphonic acid, very-large molecule, very-small molecule, zwitterion); these are the same classes for which a generic neutral-organic backbone is most likely to fail. Each cell reports the train and test counts. The audit informs two complementary actions taken in the deployed model: (i) where targeted external data acquisition can close a gap (e.g., AqSolDB salts for solubility, brought in during the Stage 3 augmented fine tune), the gap is filled; (i) where no public condition-matched data was obtainable within the scope of this work (e.g., MNSol/DISSOLVE polyols for hydration free energy require institutional license access), the affected (property, class) cell is declared out-of-scope under the per-property G+ filter. Table 3: Per-(property Ă scarcity-class) coverage matrix on the training set, with the corresponding test counts in parentheses for transparency. Each cell is ntrainâ(ntest)n_train~(n_test) and counts molecules of the indicated structural class for which the property is labeled. The classes correspond to the chemistry exclusions encoded in the per-property G+ applicability-domain filters of Section 2.14. Where the training count is small relative to the test count for a given (property, class) cell, the class is a candidate for either targeted data augmentation (e.g., AqSolDB salts for solubility, which closes the salt cell) or for declaration as out-of-scope under the per-property G+ filter (e.g., HFE polyols, where MNSol/DISSOLVE coverage is access-restricted). Property salt ion polyol sulf./phos. very-large very-small zwitt. total logâĄS S 962 (0) 376 (76) 109 (29) 752 (102) 405 (3) 246 (33) 1287 (306) 9978 (2700) logâĄP P â 86 (19) 31 (13) 254 (70) 192 (2) 377 (136) 349 (120) 9975 (1570) ÎâGhyd G_ hyd â â 2 (2) â â 292 (140) 52 (4) 1862 (269) TbT_b 199 (4) 28 (7) 36 (5) 129 (20) 66 (1) 465 (142) 318 (122) 7443 (2106) VP 420 (0) 51 (9) 68 (5) 194 (27) 174 (2) 235 (30) 589 (144) 8808 (1651) η 7 (0) 3 (0) 34 (3) 17 (0) 10 (0) 868 (141) 146 (2) 7821 (331) TmT_m â 77 (16) 91 (21) 185 (66) 220 (4) 402 (105) 682 (733) 14938 (6206) FP 63 (8) 15 (14) 48 (7) 27 (19) 29 (2) 469 (140) 441 (287) 10700 (2840) CpC_p â 1 (0) 65 (3) 6 (0) 15 (0) 339 (139) 86 (2) 3900 (297) The dataset is split at the molecule level into training (80 %80\, 37 ), validation (10 %10\, 37 ), and test (10 %10\, 37 ) sets using a hybrid strategy that balances two competing objectives: (i) evaluating the modelâs ability to generalize to structurally novel molecules, and (i) ensuring sufficient representation of rare properties in evaluation sets. The rationale behind this split design is that a simple random partition, while maximizing statistical power, would underestimate the generalization challenge because structurally similar molecules (e.g., members of the same homologous series) would appear in both training and test sets. Conversely, a purely scaffold-based split would provide the most rigorous generalization test but could leave data-scarce properties (hydration free energy, viscosity, heat capacity) with too few test samples for reliable evaluation. The hybrid strategy resolves this tension as follows. For molecules possessing only common properties (solubility, logâĄP P, boiling point, vapor pressure, melting point, flash point), a scaffold-based split [63] is employed: Murcko scaffolds (the core ring systems of each molecule) are computed via RDKit, and entire scaffold groups are assigned to the same partition to ensure that the test set evaluates generalization to structurally novel chemical series rather than interpolation within familiar scaffolds. For molecules possessing rare properties (hydration free energy, viscosity, heat capacity), a greedy assignment strategy ensures that each property has at least 50 unique molecules in both the validation and test sets, preventing the evaluation from being dominated by statistical noise. All rows corresponding to a given molecule (i.e., measurements at different temperatures) are assigned to the same split, preventing data leakage between partitions. This hybrid approach yields the final split sizes of 22 01522\,015 training, 70767076 validation, and 88778877 test molecules (Table 2). The two-pass procedure is summarized in Algorithm 1. Pass 1 reserves data-scarce property samples first, so that subsequent scaffold-based assignment cannot accidentally drain those properties below the 50-molecule per-split floor. Pass 2 then assigns the remaining (common-property) molecules by their Murcko scaffold, with whole scaffold groups routed to a single split. Algorithm 1 Hybrid scaffold + property-balanced splitting. 1:Molecule set âłM with property-availability mask mâ0,19p_mâ\0,1\^9 for each mââłm 2:Scarce property set scarce=ÎâGhyd,η,CpP_scarce=\ G_hyd,η,C_p\ 3:Per-split minimum count nmin=50n_ =50 per scarce property; target ratios (0.80,0.10,0.10)(0.80,0.10,0.10) for (train, val, test) 4:Disjoint partitions âłtrain,âłval,âłtestM_train,M_val,M_test 5: 6:Pass 1: greedy scarce-property assignment 7:âłvalââ M_valâ ; âłtestââ M_testâ 8:for each property pâscarcep _scarce do 9: âłpâmââł:mâ[p]=1M_pâ\m :p_m[p]=1\ âł molecules with property p 10: Shuffle âłpM_p with fixed seed 11: while |âłvalâ©âłp|<nmin|M_val _p|<n_ do 12: move next m from âłpM_p to âłvalM_val 13: end while 14: while |âłtestâ©âłp|<nmin|M_test _p|<n_ do 15: move next m from âłpM_p to âłtestM_test 16: end while 17:end for 18: 19:Pass 2: scaffold-grouped assignment for the remainder 20:âłrestââłâ(âłvalâȘâłtest)M_rest (M_val _test) 21:Compute Murcko scaffold sâ(m)s(m) for each mââłrestm _rest âł RDKit 22:Group âłrestM_rest by scaffold: ssâ\G_s\_s 23:Order scaffolds S by descending |s||G_s| âł largest scaffolds placed first 24:for each scaffold group sG_s do 25: Assign entire sG_s to the split currently furthest below its target ratio 26:end for 27: 28:Pass 3: row-level consolidation 29:for each molecule m assigned to split Ïâtrain,val,testÏâ\train,val,test\ do 30: Route all measurement rows of m (across temperatures and conditions) to Ï 31:end for 32:return âłtrain,âłval,âłtestM_train,M_val,M_test Three-dimensional molecular conformers are required as input to the SchNet geometry encoder. The inclusion of explicit three-dimensional geometry is motivated by the fact that many thermophysical properties are governed by spatial effects that cannot be inferred from two-dimensional topology alone: vapor pressure depends on molecular surface area and shape, viscosity on molecular packing and flow resistance, and hydration free energy on the three-dimensional cavity that a solute carves in the solvent. The two-dimensional graph encodes which atoms are bonded but not how far apart non-bonded atoms are in space, nor the dihedral angles that define rotational conformations. By providing the SchNet encoder with explicit Cartesian coordinates, the model gains access to through-space distances, molecular shape descriptors, and steric information that complement the topological and sequential features captured by the GCN and Transformer encoders, respectively. Conformers are precomputed for all 37 96837\,968 unique molecules and cached for efficient training. Each canonical SMILES is converted to an RDKit molecule object with explicit hydrogens, and a three-dimensional embedding is generated using the ETKDGv3 algorithm [48] with a fixed random seed for reproducibility. If ETKDGv3 fails (approximately 1 %1\, 37 of molecules), a fallback with random initial coordinates is attempted. The embedded conformer is subsequently optimized using the MMFF94 force field [20] for up to 200 iterations, after which explicit hydrogens are removed to yield atomic numbers and Cartesian coordinates consistent with the heavy-atom graph representation. As a final fallback for molecules where 3D embedding fails entirely, two-dimensional coordinates are computed and extended with z=0z=0. Full coverage (100 %100\, 37 ) was achieved across the dataset. A notable limitation of this approach is that a single static conformer from force-field optimization cannot represent the Boltzmann ensemble of conformations that a flexible molecule populates in solution. For rigid aromatic compounds, the single-conformer approximation is reasonable, but for molecules with many rotatable bonds (e.g., long-chain alkanes, polyethers), the properties of interest (particularly viscosity and solvation free energy) may depend on the conformational distribution rather than any single geometry. The geometry gate mechanism (Section 2.10) partially mitigates this limitation by allowing the model to suppress the 3D contribution for molecules where the single conformer is unrepresentative, but this is an imperfect solution. Future work could address this by either (i) generating multiple conformers per molecule and aggregating their SchNet representations, following the multi-conformer ensemble approach of Uni-Mol [65], or (i) incorporating conformer generation uncertainty as an additional input feature. 2.2 Self-supervised pretraining of the molecular backbone The multimodal backbone is pretrained through a self-supervised objective on a large unlabeled molecular corpus before any property labels are introduced. The rationale is the same as for foundation models in natural language processing and vision: labeled thermophysical data are scarce (tens of thousands of molecules per property at best), while unlabeled molecular data; drug-like compounds, commercially available chemicals, and structures catalogued in public repositories; are available in the tens to hundreds of millions. A backbone that has been exposed to a broad chemical distribution before it sees any property labels starts supervised training from a chemistry-aware initialization rather than a random one, and this initialization is carried through every downstream prediction head. A corpus of 500 000500\,000 drug-like and industrial-chemical SMILES strings was assembled from PubChem [27] by filtering to molecules with 5â60 heavy atoms, passing RDKit canonicalization, and excluding charged or multi-component entries. The corpus is intentionally disjoint from the property-labeled training data: molecules that appear in any of the nine property datasets are excluded from the pretraining corpus to prevent information leakage from property labels into the pretraining objective. The backbone is pretrained with three complementary self-supervised signals that together expose the GCN, Transformer and SchNet encoders to invariance, structural-recovery and language-recovery pressure before any thermophysical label is seen. For each molecule m in a mini-batch of size B, two chemically equivalent views are constructed: the canonical SMILES smcs^c_m and a randomly enumerated SMILES smrs^r_m produced by RDKitâs randomised traversal. Each view is independently routed through the multimodal backbone, and the resulting unified embeddings are â2 _2-normalised to unit length, mc,mrâdâ1z^c_m,z^r_m ^d-1 with d=512d=512. The contrastive InfoNCE [oord2018cpc] loss for the batch is âNCE=â1Bââm=1BlogâĄexpâĄ(mcâ mr/Ï)âj=1BexpâĄ(mcâ jr/Ï),L_NCE\;=\;- 1B _m=1^B \! (z^c_m·z^r_m/Ï ) _j=1^B \! (z^c_m·z^r_j/Ï ), (1) where Ï=0.07Ï=0.07 is the temperature parameter and the denominator runs over the Bâ1B-1 in-batch negatives plus the positive at j=mj=m. The objective drives the fused embedding to be insensitive to the SMILES traversal order while remaining discriminative across distinct molecules. Let mââNmĂ39X_m ^N_mĂ 39 denote the atom-feature matrix of molecule m (Section 2.3), where the first 11 dimensions of each row encode the atomic-number one-hot. A masking set âłmâ1,âŠ,NmM_mâ\1,âŠ,N_m\ is sampled uniformly at random with |âłm|/Nmâ0.15|M_m|/N_mâ 0.15, and the atomic-number block of each masked row is replaced by a learnable [MASK][MASK] vector while every other feature dimension is preserved. The masked features are propagated through the GCN encoder, and a small two-layer auxiliary head gΞ:â256ââ11g_Ξ:R^256 ^11 is trained on the penultimate GCN layer to recover the original atomic-number distribution. The mean-squared-error reconstruction loss is âMAR=1Bââm=1B1|âłm|ââiââłmâgΞâ(m,i)âm,i(atom)â22,L_MAR\;=\; 1B _m=1^B 1|M_m| _i _m \|g_Ξ\! (h_m,i )-x_m,i^(atom) \|_2^2, (2) where m,ih_m,i is the masked-input GCN representation of atom i and m,i(atom)x_m,i^(atom) is the original 11-dimensional atomic-number target. The objective forces the GCN to use bond context (neighbour atomic numbers, hybridisation, ring membership) to infer the identity of the masked atom, which is the graph analogue of BERTâs masked-language modelling [14] adapted to chemistry. Independently from the graph masking above, a masking set mS_m of token positions in the SMILES sequence smcs^c_m is sampled uniformly with |m|/|smc|â0.15|S_m|/|s^c_m|â 0.15 and replaced by a learnable [MASK][MASK] token. The full sequence is processed by the SMILES Transformer (Section 2.6) and a tied output projection Woutââ|V|Ă512W_out ^|V|Ă 512 produces a logit vector at every position, where |V|=50|V|=50 is the vocabulary size. The cross-entropy loss is âMTP=â1Bââm=1B1|m|ââtâmlogâĄPâ(sm,tc|s~mc;Ξ),L_MTP\;=\;- 1B _m=1^B 1|S_m| _t _m P\! (s^c_m,t\, |\, s^c_m;\,Ξ ), (3) where s~mc s^c_m is the masked sequence and P(â âŁs~mc;Ξ)P(· s^c_m;Ξ) is the softmax over the Transformerâs output logits at the masked position. The three terms are combined with unit weights âSSL=âNCE+âMAR+âMTP.L_SSL\;=\;L_NCE\;+\;L_MAR\;+\;L_MTP. (4) Equal weighting is justified empirically: (i) the three losses live on comparable numerical scales after the first warm-up epoch (each in the range 0.50.5â33); (i) weight searches over the simplex (αNCE,αMAR,αMTP):αiâ„0,âiαi=3\( _NCE, _MAR, _MTP): _iâ„ 0, _i _i=3\ within the 0.5,1.0,1.5,2.03\0.5,1.0,1.5,2.0\^3 grid did not yield a downstream improvement larger than ÎâR2â0.005 R^2â 0.005 on the validation set, well within the seed-to-seed noise of multi-task training; and (i) unit weights preserve the interpretability that each modality is exposed to one supervisory signal of comparable strength. Temperature-dependent auxiliary inputs are clamped to a reference state (T=298.15 KT=$298.15\,K$, P=101 325 PaP=$101\,325\,Pa$) throughout pretraining so that the backbone learns chemistry-aware representations first; the experimental encoder is fine tuned only during the supervised stage. Gradients propagate through the entire backbone (GCN + Transformer + SchNet + cross-modal fusion) for all three terms, ensuring that the fused embedding is jointly shaped by invariance, atom-recovery and language-recovery pressure rather than dominated by any single modality. Pretraining runs for 89 epochs (paused before the 100-epoch budget after validation-loss plateau) with AdamW (α=5Ă10â4α=5Ă 10^-4, weight decay 1Ă10â51Ă 10^-5), batch size 256, linear warmup over the first 2 epochs, and cosine decay thereafter. The pretrained backbone weights (GCN encoder, Transformer encoder, SchNet encoder, cross-modal fusion module) are checkpointed at the end of pretraining and then used as the initialization for the supervised multi-task training described in Section 2.13. The auxiliary encoders (experimental conditions and molecular descriptors), the condition module stack (Section 2.9), and the nine property-specific prediction heads are initialized randomly at the start of supervised training and learned from scratch. Ablation against a randomly initialized backbone (Section 3.5.1) confirms that self-supervised pretraining is the single largest non-architectural contributor to downstream performance, on the order of +0.03+0.03 in mean test R2R^2. The improvement is largest on data-scarce properties (hydration free energy, viscosity, heat capacity, solubility), where the chemistry-aware initialization compensates for the limited supervised training signal available per task. 2.3 Molecular featurization Molecular featurization translates the raw chemical structure into numerical representations suitable for neural network processing. The featurization scheme is designed to encode the chemical identity, local bonding environment, and electronic properties of each atom and bond in a compact, fixed-dimensional vector. These features serve as the input to the GCN encoder (atom and bond features) and the SMILES Transformer encoder (tokenized sequences), while the SchNet encoder operates directly on atomic numbers and Cartesian coordinates from the conformer generation step. The choice and granularity of features directly affect the information available to the model: too few features limit the modelâs ability to distinguish chemically distinct environments, while excessive dimensionality can dilute the learning signal and increase computational cost. Each atom in the molecular graph is represented by a 39-dimensional feature vector. This vector encodes the atomic number as an 11-dimensional one-hot vector over the ten most common elements (C, N, O, F, P, S, Cl, Br, I, Si) plus an âotherâ category; the degree (number of bonded neighbors) as a 7-dimensional one-hot encoding (dâ0,âŠ,6dâ\0,âŠ,6\); formal charge as a 5-dimensional one-hot encoding (qââ2,âŠ,+2qâ\-2,âŠ,+2\); the number of attached hydrogens (5 dimensions, nHâ0,âŠ,4n_Hâ\0,âŠ,4\); hybridization state as a 5-dimensional one-hot encoding over sp, sp2, sp3, sp3d, and sp3d2; a binary aromaticity indicator; a binary ring membership indicator; and chirality as a 4-dimensional one-hot encoding over unspecified, R, S, and other. Each bond is represented by a 12-dimensional feature vector comprising a 4-dimensional one-hot encoding of bond type (single, double, triple, aromatic), a binary conjugation indicator, a binary ring membership indicator, and a 6-dimensional one-hot encoding of stereo configuration (none, Z, E, cis, trans, other). SMILES strings are tokenized at the character level using a vocabulary of 50 tokens, including special tokens for padding, start-of-sequence, end-of-sequence, and unknown characters, alongside 46 chemistry tokens corresponding to atoms, bonds, brackets, and ring-closure digits. Two-character tokens (Cl, Br, Si, Se, Te) are handled with priority matching. Sequences are padded or truncated to a maximum length of 256 tokens. 2.4 Dual-encoder rationale The MultiPUFFIN encoder architecture deliberately processes each molecular modality through a separate, specialized encoder before fusing the resulting embeddings, rather than pre-aligning representations (as in contrastive approaches such as MoleculeSTM [35]) or encoding a single modality (3D conformers) through a unified backbone (as in Uni-Mol [65], which processes only 3D molecular geometry). This design is motivated by two key observations. First, the three molecular representations (2D graphs, 1D SMILES strings, and 3D conformers) encode fundamentally different and complementary aspects of molecular identity. The molecular graph captures topological connectivity, ring systems, and local functional group patterns through its adjacency structure; the SMILES string encodes the same molecule as a linear sequence in which long-range syntactic dependencies (e.g., ring-closure digits linking distant atoms) provide a complementary inductive bias that differs from the spatial locality inherent to message passing; and the 3D conformer captures through-space distances, dihedral angles, and steric effects that are entirely absent from the 2D representation. A moleculeâs thermophysical properties are governed by the combination of all three: topology determines the presence of hydrogen-bonding groups, sequential patterns capture substituent effects on electron density, and three-dimensional shape governs molecular packing and solvation. By processing each modality through a dedicated encoder optimized for that representationâs structure, the model preserves modality-specific features that would be lost or diluted in a shared-encoder architecture. Second, modality-specific encoders enable graceful degradation when one modality is unavailable. In practice, 3D conformers may fail to generate for certain molecules or may be unreliable for flexible macrocycles; by keeping the 3D branch separate, the model can fall back to the 2D/1D fusion without architectural modification, using the geometry gate (Section 2.10) to automatically suppress the missing modalityâs contribution. 2.5 Graph convolutional network encoder The GCN encoder processes the two-dimensional molecular graph and contributes topological information to the unified embedding: which atoms are bonded to which, the local neighbourhood that each atom sits in, and the additive functional-group structure that drives bulk thermophysical properties. Formally, a molecule is represented as an undirected graph G=(V,E)G=(V,E) where the nodes V correspond to atoms (the node feature matrix is ââNĂ39X ^NĂ 39, encoding atomic number, hybridisation, formal charge and stereochemistry) and the edges E correspond to covalent bonds (with bond-type and ring-membership features). The graph is processed through Lg=4L_g=4 KipfâWelling graph convolutional layers [28], (l+1)=Ïâ(~â1/2â~â~â1/2â(l)â(l)),H^(l+1)=Ï\! ( D^-1/2\, A\, D^-1/2\,H^(l)\,W^(l) ), (5) where ~=+N A=A+I_N is the adjacency matrix with self-loops, ~ D is the corresponding degree matrix, (l)ââdlĂdl+1W^(l) ^d_lĂ d_l+1 are learnable weights and Ï is the ReLU non-linearity; the hidden dimension is dg=256d_g=256 at every layer. The information that this operator extracts is best read off the unrolled message-passing equation: at layer l+1l+1, the embedding of atom i becomes i(l+1)=Ïâ(âjââ(i)âȘi1d~iâd~jâ(l)âj(l)),h_i^(l+1)\;=\;Ï\! ( _j (i)âȘ\i\ 1 d_i\, d_j\;W^(l)\,h_j^(l) ), (6) so that after LgL_g layers each atom embedding aggregates information from a neighbourhood of radius up to LgL_g bonds. The fixed receptive field of Lg=4L_g=4 bonds matches the topological scale at which functional-group additivity (Joback-style group contribution) operates, which is precisely the regime where conventional thermodynamic correlations are most predictive. Long-range substituent effects beyond four bonds, by construction, must be supplied by the SMILES Transformer or the 3D SchNet branches; this is the principled basis for the multimodal split. Each GCN layer is followed by GraphNorm [7] and BatchNorm [23] for training stability, with dropout p=0.15p=0.15. Residual connections [21] are used at every layer (the first layer applies a linear projection from 3939 to 256256, subsequent layers add directly). The graph-level molecular embedding is obtained by sum pooling over all node representations, graph=âi=1Ni(Lg)h_graph= _i=1^Nh_i^(L_g), followed by a two-layer projection graph=2âReLUâ(Dropoutâ(1âgraph+1))+2ââ512,z_graph=W_2\,ReLU\! (Dropout\! (W_1\,h_graph+b_1 ) )+b_2\;â\;R^512, (7) into the shared 512-dimensional embedding space. Sum pooling rather than mean pooling is chosen so that the graph-level embedding is sensitive to molecular size, which is itself a primary driver of boiling point, vapour pressure and viscosity. 2.6 SMILES Transformer encoder The SMILES Transformer encoder contributes long-range syntactic information to the unified embedding: substituent patterns that span tens of tokens, ring-closure couplings that define macrocyclic topology, and stereochemistry / valence / formal-charge cues that the SMILES grammar carries explicitly but that the 2D graph carries only through node features. The Transformer is the natural complement to the GCN: where the GCN sums local neighbourhoods of bounded radius, self-attention can directly relate any pair of tokens at arbitrary sequence distance, with attention weights αiâj(h)=softmaxjâ(i(h)â j(h)/dk) _ij^(h)=softmax_j\! (q_i^(h)·k_j^(h)/ d_k ) at head h. The information added to the unified embedding by this branch is the part of the structureâproperty mapping that depends on long-range substituent context (electron-withdrawing groups across rings, conjugation chains, position-dependent grammar of SMILES enumerations), which is precisely the regime where the GCNâs bounded receptive field is silent. Each token tit_i is mapped to a learned embedding vector iââ512e_i ^512, scaled by dmodel d_model, and augmented with sinusoidal positional encodings [58]: PE(pos,2âi)=sinâĄ(pos100002âi/dmodel),PE(pos,2âi+1)=cosâĄ(pos100002âi/dmodel)PE_(pos,2i)= ( pos10000^2i/d_model ), _(pos,2i+1)= ( pos10000^2i/d_model ) (8) The sequence is processed through Ls=6L_s=6 pre-norm Transformer encoder layers [58], each comprising multi-head self-attention with h=8h=8 heads (key dimension dk=64d_k=64) and a feed-forward sublayer with GELU activation [22], a model dimension dmodel=512d_model=512, and a feed-forward dimension dff=2048d_f=2048. Pre-norm architecture (LayerNorm applied before the attention and feed-forward sublayers) is used for improved training stability [64], with dropout (p=0.1p=0.1) applied within both sublayers. The molecule-level representation is obtained by extracting the hidden state corresponding to the start-of-sequence token (the CLS-equivalent token at position 0), followed by a two-layer output projection with GELU activation mapping to SMILESââ512z_SMILES ^512. 2.7 SchNet 3D geometry encoder The SchNet encoder [51] introduces the third modality: three-dimensional molecular geometry. The inclusion of a dedicated 3D encoder is motivated by the fact that many thermophysical properties are fundamentally governed by spatial effects that cannot be inferred from 2D topology alone. Vapor pressure and boiling point depend on intermolecular interactions whose strength is modulated by molecular shape and surface accessibility; viscosity is governed by molecular packing and the ease of flow past neighboring molecules; and hydration free energy reflects the three-dimensional cavity that a solute carves in the solvent. While the 2D graph captures which atoms are bonded and the SMILES string encodes this connectivity sequentially, neither representation encodes through-space distances between non-bonded atoms, dihedral angles that define rotational conformations, or the overall three-dimensional shape of the molecule. SchNet was selected for this branch because its continuous-filter convolution framework naturally handles the continuous and rotationally invariant nature of interatomic distances, and because it has demonstrated strong performance on quantum-chemical property prediction benchmarks [51]. The encoder operates on atomic numbers and Cartesian coordinates produced by the conformer generation pipeline (Section 1). Each atom is embedded by its atomic number into a learned vector iââ256e_i ^256 via an embedding lookup over 101 possible atomic numbers. Interatomic distances diâj=âiâjâd_ij=\|r_i-r_j\| are computed for all atom pairs within a cutoff radius rc=10 Ă r_c=$10\, $ and expanded into K=50K=50 radial basis functions using Gaussian smearing: ekâ(diâj)=expâĄ(â12âÏ2â(diâjâÎŒk)2),k=1,âŠ,Ke_k(d_ij)= (- 12Ï^2(d_ij- _k)^2 ), k=1,âŠ,K (9) where ÎŒkk=1K\ _k\_k=1^K are evenly spaced centers from 0 to rcr_c and Ï is the spacing between centers. The encoder applies Lint=4L_int=4 interaction blocks, each performing continuous-filter convolution: i(l+1)=i(l)+UpdateMLPâ(âjââ(i)(l)â(â(diâj))âj(l))x_i^(l+1)=x_i^(l)+UpdateMLP ( _j (i)W^(l) (e(d_ij) ) _j^(l) ) (10) where â(i)=j:diâj<rc,jâ iN(i)=\j:d_ij<r_c,jâ i\ is the neighborhood of atom i, (l)â(â )W^(l)(·) is a filter-generating network (a two-layer MLP with SiLU activation that maps the 50-dimensional Gaussian-expanded distance to 256-dimensional convolution filters), and â denotes element-wise multiplication. A residual connection preserves the input representation at each block. The molecule-level 3D embedding is obtained through global sum pooling followed by a two-layer projection with SiLU activation and dropout, yielding 3Dââ512z_3D ^512. The SchNet encoder operates entirely in float32 precision, even when mixed-precision training is employed elsewhere in the model. Distance computations and Gaussian smearing require full precision to avoid numerical instabilities from float16 underflow in exponential operations. 2.8 Auxiliary input encoders In addition to the three structural encoders (GCN, Transformer, SchNet), MultiPUFFIN includes two auxiliary input encoders that process non-structural molecular information when available: an experimental encoder and a descriptor encoder. These auxiliary encoders address a fundamental limitation of purely structure-based molecular models: the inability to condition predictions on thermodynamic state variables and readily available molecular descriptors. This limitation is particularly consequential for condition-dependent properties such as vapor pressure, viscosity, and heat capacity, where the same molecule exhibits vastly different values at different temperatures. A model that receives only a SMILES string or molecular graph cannot, by construction, distinguish between measurements at 300 K300\,K and 500 K500\,K for the same molecule. The auxiliary encoders resolve this by providing the model with explicit access to the thermodynamic conditions under which each measurement was taken, as well as precomputed descriptors that encode coarse-grained physicochemical information. The experimental encoder is a two-layer MLP (6â64â5126â 64â 512, ReLU activation, 38 40038\,400 parameters) that processes a 6-dimensional vector of experimental conditions associated with each data point. These conditions include the measurement temperature, pressure, and other property-specific auxiliary variables (e.g., solvent composition for partition coefficients). The experimental conditions are particularly critical: it is this information, combined with the domain-informed prediction heads (Section 2.11), that enables MultiPUFFIN to produce thermodynamically meaningful temperature-dependent predictions, a capability that SMILES-based foundation models such as ChemBERTa-2 fundamentally lack. When experimental conditions are not available for a given data point, a learned missing-data embedding replaces the zero vector, allowing the model to distinguish between âthis condition is zeroâ and âthis condition is unavailable.â The descriptor encoder is a two-layer MLP (13â128â51213â 128â 512, ReLU activation, 84 86484\,864 parameters) that processes a 13-dimensional vector of precomputed RDKit molecular descriptors [46]. The thirteen descriptors are: (1) molecular weight (MW), (2) topological polar surface area (TPSA), (3) number of hydrogen bond donors (HBD), (4) number of hydrogen bond acceptors (HBA), (5) number of rotatable bonds (RotBonds), (6) number of aromatic rings (ArRings), (7) number of aliphatic rings (AlRings), (8) number of heavy atoms (HeavyAtomCount), (9) the WildmanâCrippen logâĄP P estimate (MolLogP), (10) the WildmanâCrippen molar refractivity (MolMR), (11) fraction of sâp3sp^3 carbons (Fsp3), (12) the number of stereocenters (NumStereo), and (13) the number of radical electrons (NumRadicalElectrons). Descriptors 1â11 encode coarse-grained physicochemical properties that complement the learned representations from the structural encoders; descriptors 12â13 provide stereochemical and electronic information not easily extracted from the 2D graph alone. One potential concern is that the inclusion of the WildmanâCrippen logâĄP P estimate (MolLogP, descriptor 9) in the descriptor vector may advantage the logâĄP P prediction task by providing a pre-computed approximation of the target value directly as input. This constitutes a mild form of information leakage for the logâĄP P head: the model can learn to use the Crippen estimate as a prior and then correct the residual from the structural representations. The effect is intentional ; the descriptor encoder is designed to provide the model with the same coarse physicochemical priors that an experienced chemist would use ; but downstream users comparing MultiPUFFIN logâĄP P performance against baselines that do not have access to Crippen descriptors should account for this asymmetry. The descriptor encoder likewise includes a learned missing-data embedding for molecules with incomplete descriptor coverage. The outputs of both auxiliary encoders are projected to 512 dimensions and concatenated with the unified structural embedding u before the final projection layer, providing additional conditioning information to the downstream prediction heads. Together, the two auxiliary encoders contribute 123 264123\,264 parameters (0.35%0.35\% of the total model), and their inclusion allows the model to use heterogeneous input data when available while handling missing auxiliary information through the learned missing-data embeddings. 2.9 Condition-aware embedding refinement Several of the predicted properties depend not only on molecular structure but on additional thermodynamic or test-protocol conditions: vapor pressure, viscosity, heat capacity, and solubility depend on temperature; viscosity and heat capacity additionally on pressure; aqueous solubility, logâĄP/logâĄD P/ D, and hydration free energy on solution pH for ionizable solutes; melting point on the crystal polymorph adopted by the solid; and flash point on whether the test was conducted in an open or closed cup. To make the model condition-aware in a modular and uniformly-architectured way, MultiPUFFIN inserts a stack of identity-initialized condition modules between the fused embedding and the prediction heads. Each module is a small 2-token cross-attention block (one molecule token, one condition token) with zero-initialized additive residuals, so that at fresh initialization it is an exact identity on the molecule token regardless of the condition input. This property allows any module to be inserted into a warm-restarted backbone without disturbing the previously-trained head calibration; training then moves a module away from the identity only where the corresponding condition has variance in the labeled data and reduces the multi-task loss. Each condition module constructs a 2-token sequence (the unified molecular embedding u and a condition token c) and applies one self-attention layer with 4 heads followed by a 512-dimensional position-wise feed-forward network with GELU activation, LayerNorm, and dropout =0.1=0.1. The condition token c is constructed by a small encoder specific to each condition variable (sinusoidal positional encoding of the scalar value followed by a 2-layer MLP for continuous variables; a learned embedding lookup for discrete labels). The molecule token after the attention block is returned as the conditioned embedding. Each module contributes 0.99â1.28 M parameters, totalling ⌠5.6 M for the full stack (⌠20% of the deployed model). When a row carries a missing or sentinel condition value (e.g., pH unknown), a learned âunknownâ embedding replaces c for that row, so the module degrades gracefully to identity in the absence of condition information. The deployed stack consists of: âą TConditioner (temperature). c from a sinusoidal encoding of (Tâ298.15 K)(T-$298.15\,K$) projected through MLP. ⌠1 283 0001\,283\,000 parameters. âą PHConditioner (solution pH). c from a sinusoidal encoding of (pHâ7)(pH-7) projected through MLP, with a learned âunknown-pHâ fallback for rows where pH is not reported. ⌠1 284 0001\,284\,000 parameters. âą PressureConditioner (ambient pressure). c from a sinusoidal encoding of logâĄ(P/atm) (P/atm), with a learned unknown-pressure fallback. ⌠1 284 0001\,284\,000 parameters. âą PolymorphEmbedding (crystal-form identifier). c from a learned embedding over K=16K=16 polymorph slots with k=0k=0 reserved for âunknownâ. ⌠994 000994\,000 parameters. âą MethodFlag (test-protocol identifier, e.g., open vs. closed cup for FP). c from a learned embedding over K=4K=4 method slots with k=0k=0 reserved for âunknownâ. ⌠989 000989\,000 parameters. Conditioner outputs are routed to heads only where the corresponding condition is physically relevant, and chained when multiple conditions matter. The deployed routing is: âą Vapor pressure: Tu_T (T-conditioned). âą Viscosity, heat capacity: TâPu_Tâ P (T- then P-conditioned). âą Solubility: Tu_T (T-conditioned). The pH conditioner is routed off solubility because the labeled logâĄS S data does not vary in pH within the AqSolDB / ESOL training subset, so applying pH-conditioning would introduce a free degree of freedom with no supervisory signal. âą logP: pHu_pH (pH-conditioned). The labeled logâĄP/logâĄD P/ D data spans multiple pH values from ChEMBL bioactivity records, so the pH module receives a real gradient signal here. âą Boiling point, melting point, flash point, hydration free energy: unconditioned u. The polymorph and method modules exist in the model but are not routed to any head in the deployed model because the corresponding labels are not available in the current training data; the modules remain at their identity initialization in the deployed checkpoint and are documented here as architectural slots ready to be activated when polymorph- or method-tagged data become available below. The five-module stack is condition-aware by construction, but the magnitude of each moduleâs contribution at deployment is controlled by how much the corresponding condition varies in the training set. The TConditioner sees temperature spanning two hundred kelvin in the multi-T thermophysical data and contributes substantially to the four T-dependent property predictions. The PressureConditioner sees roughly 1.6Ă1041.6Ă 10^4 rows of viscosity and heat-capacity data with explicit pressure values from NIST ThermoML and contributes a modest correction. The PHConditioner sees a non-trivial pH range in the appended ChEMBL logâĄD D rows but is dominated by the pH=7.4=7.4 subset; its corrections are accordingly small. The polymorph and method modules contribute zero in the deployed model. The condition-dependency matrix of the conditionâdata analysis below reports these distinctions explicitly: it lists each condition variable as either âhasâ (an active routing in the deployed model with sufficient training-data variance to learn from), âhas, sparseâ (an active routing with limited training-data variance), or âmodule ready, no dataâ (a module that exists architecturally but is awaiting paired-condition data). 2.10 Cross-modal fusion The fusion mechanism integrates the outputs of all five encoders (the three structural encoders (GCN, Transformer, SchNet) and the two auxiliary encoders (experimental conditions, molecular descriptors)) into a single unified molecular embedding through a hierarchical strategy. The central challenge in multimodal fusion is deciding how to combine representations that encode overlapping but non-identical information about the same molecule. Simple concatenation or element-wise averaging treats each modality as independent and equally important, ignoring the fact that different modalities may be more or less informative depending on the moleculeâs structure. For example, a rigid planar aromatic compound may be well-characterized by its 2D topology, whereas a flexible molecule with multiple rotatable bonds may require 3D geometry to capture the conformational effects that dominate its thermophysical behavior. This motivates a learned, molecule-adaptive fusion strategy. The first stage of fusion employs bidirectional cross-modal attention between the two primary branches (GCN and Transformer). To illustrate why cross-attention is preferable to simply concatenating the two embeddings, consider a concrete example: a substituted aromatic molecule where a distant electron-withdrawing group significantly affects the ringâs electron density and, consequently, the moleculeâs solvation behavior. The GCN encoder captures the local ring environment through message passing, while the Transformer encoder captures the long-range substituent effect through self-attention over the SMILES sequence. Simple concatenation would juxtapose these two independent representations and leave it to the downstream layers to discover the relationship between them. Cross-attention, by contrast, enables each encoderâs representation to be enriched by information from the other modality before the fusion step. When the GCN embedding serves as the query and the Transformer embedding as key/value, the attention mechanism highlights those aspects of the graph representation that are most relevant given the sequential context, effectively allowing the GCN to âaskâ the Transformer which sequential features are important for contextualizing its local topology. The reverse direction allows the Transformer to ground its long-range syntactic features in explicit bond-topology information from the GCN. This bidirectional exchange allows the model to identify and amplify complementary features across modalities: the GCNâs local functional-group patterns are contextualized by the Transformerâs global sequential dependencies, and vice versa. In practice, cross-attention performs a soft, learned alignment between the two representation spaces, ensuring that the subsequent fusion operates on representations that are already aware of each otherâs content. Formally, two cross-attention modules enable each branch to attend to the other: graphâČ=CrossAttnâ(=graph,=SMILES,=SMILES)z_graph =CrossAttn(Q=z_graph,\;K=z_SMILES,\;V=z_SMILES) (11) SMILESâČ=CrossAttnâ(=SMILES,=graph,=graph)z_SMILES =CrossAttn(Q=z_SMILES,\;K=z_graph,\;V=z_graph) (12) Each cross-attention module uses multi-head attention with h=8h=8 heads, followed by a feed-forward network with GELU activation (expansion factor 4Ă4Ă), LayerNorm, and residual connections. The use of multiple attention heads allows the model to simultaneously attend to different aspects of the cross-modal relationship (for instance, one head might focus on aromatic features while another captures hydrogen-bonding patterns), providing a richer fusion than a single attention mechanism. After cross-modal attention, the enriched GCN and Transformer embeddings must be combined into a single representation. Averaging or concatenation would treat both modalities as equally important for every dimension of the embedding, ignoring the fact that different modalities may contribute more or less information depending on the molecule and the feature dimension. To address this, a learned element-wise gating mechanism is employed: =Ïâ(MLPâ([graphâČ;SMILESâČ])),fused=âgraphâČ+(1â)âSMILESâČg=Ï (MLP([z_graph ;z_SMILES ]) ), _fused=g _graph +(1-g) _SMILES (13) where â[0,1]512gâ[0,1]^512 is an element-wise gate vector and the MLP maps from 2Ă5122Ă 512 to 512 dimensions through a hidden layer with ReLU activation. The sigmoid gate implements a soft, per-dimension selection between the two modalities: for each of the 512 embedding dimensions, the gate learns whether the graph-derived or the sequence-derived feature is more informative, and blends them accordingly. This is critical because the relative utility of each modality varies not only across molecules but also across different aspects of the representation. For instance, dimensions encoding electronegativity patterns may favor the GCN (which directly processes atomic features), while dimensions encoding long-range substituent effects may favor the Transformer. The gating mechanism (inspired by the PUFFIN fusion architecture [50]) thus provides a strictly more expressive fusion than either concatenation (which preserves all information but doubles the dimensionality) or averaging (which assumes equal importance). The SchNet 3D embedding is incorporated through a separate, hierarchically subsequent geometry gate rather than being included in the bidirectional cross-attention. This asymmetric design is deliberate for two reasons. First, the 3D conformer is a derived quantity (generated computationally from the 2D structure using force-field optimization) and is therefore less reliable than the 2D graph and SMILES representations, which derive directly from the molecular identity as encoded in chemical databases. Treating the 3D embedding on equal footing with the primary modalities in the cross-attention stage could allow conformer artifacts (e.g., suboptimal geometries from failed force-field minimization, or the inherent limitation that a single static conformer cannot represent the conformational ensemble of flexible molecules) to corrupt the core representation. Second, the geometry gate provides a natural mechanism for handling missing 3D data: for certain molecules, conformer generation may fail entirely (e.g., very large or strained structures where ETKDGv3 cannot find a valid embedding) or produce unreliable geometries (e.g., flexible macrocycles with many low-energy conformations where any single conformer is unrepresentative). In such cases, a zero vector replaces 3Dz_3D, and the gate naturally suppresses the 3D contribution without requiring any architectural modification, allowing the model to fall back on the 2D/1D fusion for a robust prediction. The geometry gate is defined as α3D=Ïâ(MLPgeoâ([fused;3D])) _3D=Ï(MLP_geo([z_fused;z_3D])), which controls the contribution of geometric information via fusedâČ=[fused;α3Dâ Projâ(3D)]z_fused =[z_fused;\; _3D·Proj(z_3D)]. The scalar gate α3Dâ[0,1] _3Dâ[0,1] allows the model to learn, on a per-molecule basis, how much weight to assign to 3D geometry, effectively answering the question âdoes the three-dimensional shape of this particular molecule provide information beyond what the 2D topology already captures?â For rigid, planar molecules (e.g., benzene derivatives), the gate may learn to assign low weight since the 3D structure adds little beyond the 2D graph. For flexible molecules with multiple low-energy conformations, the gate may assign high weight to capture the steric and intramolecular interaction effects that govern properties like viscosity and vapor pressure. The quantitative effect of this molecule-adaptive behavior is confirmed by the architectural ablation study (Section 3.5.1): removing the SchNet encoder disproportionately increases RMSE for geometry-sensitive properties such as hydration free energy (Î =+0.90=+0.90 kcal/mol) and heat capacity (Î =+11.65=+11.65 J/mol/K), while having minimal impact on properties primarily governed by 2D topology (boiling point: Î =â0.18=-0.18 K), confirming that the geometry gate selectively uses 3D information where it is most informative. The concatenated embedding is projected through a two-layer network with LayerNorm, GELU activation, and dropout to the final unified embedding ââ512u ^512, which serves as input to all downstream prediction heads. 2.11 Domain-informed prediction heads A key architectural innovation of MultiPUFFIN is the replacement of standard linear output layers with inductive bias neurons that encode established thermophysical correlations. This design principle, introduced in PUFFIN [50] for vapor pressure and extended in ExPUFFIN [47] for viscosity, is here generalized to nine simultaneous properties and organized as a per-property tournament of candidate heads, in which multiple inductive biases are trained in parallel on top of the pretrained frozen backbone and the best-performing head is selected per property on the validation set. Each candidate prediction head consists of a property-specific feed-forward neural network (FFNN) that maps the unified molecular embedding u to the parameters of the corresponding thermophysical equation, followed by an inductive bias output layer that evaluates the equation at the specified thermodynamic conditions. Formally, for property p with thermophysical equation Ïpâ(p,T) _p( Ξ_p,T): p=FFNNpâ(;p),y^p=Ïpâ(p,T) Ξ_p=FFNN_p(u;W_p), y_p= _p( Ξ_p,T) (14) Gradients propagate through Ïp _p via automatic differentiation, ensuring that the domain-informed functional form shapes parameter updates throughout the entire network. For each property, MultiPUFFIN evaluates four candidate heads trained jointly on the frozen pretrained backbone: âą Head A â primary physics equation. The head whose functional form most closely matches the established thermodynamic correlation for the property: Antoine for vapor pressure, Andrade for viscosity, van ât Hoff for solubility, Born solvation for hydration free energy, Shomate polynomial for heat capacity, and a learned 32-parameter group contribution for boiling point. For temperature-dependent properties, this head predicts equation parameters that are then evaluated at the input temperature through the physics equation, guaranteeing the correct qualitative temperature dependence (monotonic pressure increase with temperature for vapor pressure, monotonic viscosity decrease, monotonic heat capacity variation) by construction. âą Head B â Joback group contribution. Operates on a 40-dimensional Joback functional-group count vector derived from the SMILES, maps the counts through a small MLP directly to the scalar property value, and for temperature-dependent properties concatenates a normalized temperature feature. This head tests the hypothesis that classical additive group contribution, decoupled from the learned backbone embedding, can match or exceed the thermodynamically-informed head on at least some properties. âą Head C â RDKit fragment head. Identical in design to Head B but operating on the âŒ85 85-dimensional RDKit functional-group fragment-count vector, which encodes a different, more granular fragment vocabulary. Comparison between Heads B and C reveals whether the Joback or RDKit fragment decomposition is better aligned with each property. âą Head D â direct or alternative physics. For temperature-dependent properties where an alternative established equation exists (vapor pressure: DIPPR-101; viscosity: VogelâTammannâFulcher), Head D implements that alternative. This lets the tournament compare two physically grounded equations for the same property (e.g., Antoine versus DIPPR-101 for vapor pressure). For the remaining temperature-independent properties, Head D is a direct feed-forward network on the fused embedding with a learnable temperature concatenation, serving as a purely data-driven reference without equation structure. All four heads for a given property are trained simultaneously during Stage 2 of the supervised pipeline (Section 2.13). Each head receives its own loss signal from the same batch of training samples, and since the heads are structurally independent (no shared parameters except those inside the frozen backbone), there is no gradient interference among them. At the end of Stage 2, the per-property tournament winner is selected as the head with the lowest validation RMSE; no test-set information is used to make this selection. Table 4 summarizes, for each of the nine properties, the tournament head A (primary physics) and head D (alternative physics or direct), along with alternative equations evaluated in earlier iterations of this work. The table lists the functional form, number of predicted parameters, and thermodynamic motivation for each. The equation-level ablation in Section 3.5.3 reports the per-property tournament outcomes on the held-out test set. The final selected heads are shown in bold below. Table 4: Summary of all domain-informed equations evaluated as inductive bias neurons in the property-specific prediction heads. For each property, the table lists the equation selected for the final model (shown in bold) and all alternative equations tested during equation-level ablation (Section 3.5.3). ||| Ξ| denotes the number of equation parameters predicted by the head. Properties marked with â use a DirectHead in the final model because no domain-informed equation improved over the unconstrained baseline. Equations marked with â diverged during head-only training. Property Equation Functional form ||| Ξ| Thermodynamic motivation Vapor pressure Antoine [56] (Head A) log10âĄP=AâB/(T+C) _10P=A-B/(T+C) 3 ClausiusâClapeyron family; industry standard for pure liquids DIPPR-101 (Head D) log10âĄP=A+B/T+Câlog10âĄT+DâT2 _10P=A+B/T+C _10T+D\,T^2 4 Extended form capturing curvature across wide T ranges Joback GC (Head B) y^=MLPâ(Joback counts,T) y=MLP(Joback counts,T) â Additive functional-group contribution RDKit GC (Head C) y^=MLPâ(RDKit counts,T) y=MLP(RDKit counts,T) â Fragment-count baseline (alternative vocabulary) ClausiusâClapeyron (tested) lnâĄP=AâB/T P=A-B/T 2 Two-parameter latent-heat approximation Viscosity Andrade [3] (Head A) log10âĄÎ·=A+B/(T+C) _10η=A+B/(T+C) 3 Arrhenius-type activated flow over energy barriers VTF (Head D) log10âĄÎ·=A+B/(TâT0) _10η=A+B/(T-T_0) 3 VogelâTammannâFulcher for supercooled / viscous regimes Joback GC (Head B) y^=MLPâ(Joback counts,T) y=MLP(Joback counts,T) â Additive group contribution RDKit GC (Head C) y^=MLPâ(RDKit counts,T) y=MLP(RDKit counts,T) â Fragment-count baseline Solubility van ât Hoff logâĄS=logâĄSref+Aâ(1/Trefâ1/T) S= S_ref+A(1/T_ref-1/T) 2 Thermodynamic dissolution enthalpy Modified Apelblat lnâĄS=A+B/T+CâlnâĄT S=A+B/T+C T 3 Three-parameter temperature-dependent solubility DirectHead y^=FFNNâ() y=FFNN(u) â No thermodynamic constraint (reference) Boiling point Group contribution [43] Tb=T0+âi=132Îi+cT_b=T_0+ _i=1^32 _i+c 34 Additive functional group contributions Nannoolal ratio Tb=T0â (1+âwiâfi)T_b=T_0·(1+ÎŁ w_if_i) 25 Multiplicative group interaction method DirectHead y^=FFNNâ() y=FFNN(u) â No thermodynamic constraint (reference) logâĄP Pâ LFER y^=âi=124wiâfi+b y= _i=1^24w_if_i+b 25 Additive fragment lipophilicity contributions Abraham LSER logâĄP=c+eâE+sâS+aâA+bâB+vâV P=c+eE+sS+aA+bB+vV 6 Linear solvation energy relationship DirectHead y^=FFNNâ() y=FFNN(u) â No thermodynamic constraint (selected as optimal) HFE Thermodynamic decomp. y^=ÎâHâTâ ÎâS y= H-T· S 2 Gibbs free energy of solvation Born solvation model ÎâG=â(1â1/Δ)â q2/(8âÏâÏ”0âr) G=-(1-1/ )· q^2/(8Ï _0r) 2 Electrostatic ionâsolvent cavity model DirectHead y^=FFNNâ() y=FFNN(u) â No thermodynamic constraint (reference) Melting pointâ DirectHead y^=FFNNâ() y=FFNN(u) â Crystal packing effects lack closed-form correlation Yalkowsky â Tm=ÎâHfus/ÎâSfusT_m= H_fus/ S_fus 2 Enthalpyâentropy ratio (diverged during training) Flash pointâ DirectHead y^=FFNNâ() y=FFNN(u) â Complex vaporâair ignition threshold SatyanarayanaâRao TFP=a+bâ Tb+câ Tb2T_FP=a+b· T_b+c· T_b^2 3 Quadratic boiling point correlation Carroll NFPN_FP TFP=fâ(NFP,MW)T_FP=f(N_FP,MW) 2 Flash point number method Heat capacity DirectHead y^=FFNNâ() y=FFNN(u) â Molecular degrees of freedom (no universal closed form) Shomate (NIST) Cp=A+Bât+Cât2+Dât3+E/t2C_p=A+Bt+Ct^2+Dt^3+E/t^2 5 NIST polynomial temperature dependence For the domain-informed heads, the networkâs penultimate layer predicts the equation parameters p Ξ_p, which are then evaluated through the corresponding thermophysical equation at the specified temperature. Soft constraints enforce thermodynamically meaningful parameter ranges by construction, via activation functions at the parameter-output layer rather than via explicit loss penalties: âą Antoine vapor pressure head. log10âĄP=AâB/(T+C) _10P=A-B/(T+C) with A unconstrained, B=softplusâ(B~)>0B=softplus( B)>0 guaranteeing âlogâĄP/âT>0â P/â T>0 (monotonic vapor pressure with temperature), and C=â100â Ïâ(C~)â[â100,0]âKC=-100·Ï( C)â[-100,0]\,$K$ preventing the denominator from vanishing in the physical range. This is the primary vapor pressure head used in the final model; an alternative DIPPR-101 form is trained in parallel as Head D (see tournament below). âą Andrade viscosity head. log10âĄÎ·=A+B/(T+C) _10η=A+B/(T+C) with analogous constraints ensuring monotonic viscosity decrease with temperature. âą VTF viscosity head (Head D). log10âĄÎ·=A+B/(TâC) _10η=A+B/(T-C) with B>0B>0 and Câ[0,200]âKCâ[0,200]\,$K$, providing an alternative functional form that better describes supercooled liquids and viscous regimes far from ambient. âą van ât Hoff solubility head. logâĄS=logâĄSref+Aâ(1/Trefâ1/T) S= S_ref+A(1/T_ref-1/T) with Tref=298.15 KT_ref=$298.15\,K$. âą Shomate heat-capacity head. Cp=A+Bât+Cât2+Dât3+E/t2C_p=A+Bt+Ct^2+Dt^3+E/t^2 with t=T/1000t=T/1000, encoding the NIST polynomial form through five predicted coefficients. âą Born hydration free energy head. Parameterizes solvation through an effective cavity radius and dielectric correction, as in PUFFIN-style Born solvation [50]. âą Group contribution boiling point head. Learns 32 virtual group contributions end-to-end, where the GCN and fused embedding implicitly encode the group decomposition that classical Joback/Lydersen methods perform by hand. For logâĄP P, melting point, and flash point, no single established thermophysical equation captures the propertyâstructure relationship well (crystal-packing effects for melting point, macroscopic vaporâair ignition for flash point, and partitioning thermodynamics for logâĄP P), and the Head A position is occupied by a direct feed-forward network on the fused embedding. The tournament for these three properties instead compares against Joback/RDKit fragment baselines (Heads B and C) to quantify how much the learned multimodal embedding improves over a purely count-based representation. All standard heads use two hidden layers with 256 units, GELU activation, and dropout (p=0.2p=0.2), yielding approximately 197 000197\,000 parameters per head. Enhanced-capacity heads (Section 2.11) are used for vapor pressure and boiling point. Multi-task learning across nine properties creates capacity dilution in the shared backbone: the 512-dimensional unified embedding must encode sufficient information to serve all nine prediction objectives simultaneously, and the finite representational capacity means that adding tasks necessarily reduces the per-task capacity available. This dilution disproportionately affects properties with complex structureâproperty relationships that require dedicated embedding dimensions to capture their unique patterns. In practice, this manifests as performance regression for properties with complex structureâproperty relationships. To mitigate this effect, enhanced-capacity heads are employed for vapor pressure and boiling point, the two properties most affected by the expansion from six to nine tasks. These heads feature three hidden layers with 512 units per layer (instead of two layers with 256 units), LayerNorm and GELU activation at each layer, and two residual connections: one between layers 1 and 2, and a skip connection from the input embedding to the penultimate layer (with a linear projection to align dimensions): 1 _1 =GELUâ(LayerNormâ(1â+1)) =GELU(LayerNorm(W_1u+b_1)) (15) 2 _2 =GELUâ(LayerNormâ(2â1+2))+1 =GELU(LayerNorm(W_2h_1+b_2))+h_1 (16) 3 _3 =GELUâ(LayerNormâ(3â2+3))+Projâ() =GELU(LayerNorm(W_3h_2+b_3))+Proj(u) (17) p Ξ_p =outâ3+out =W_outh_3+b_out (18) The domain equation applied after parameter prediction remains identical to the standard head (Antoine for vapor pressure, group contribution for boiling point). This design provides approximately 4Ă4Ă more parameters per head (⌠800 000800\,000 vs. ⌠197 000197\,000) while preserving the domain-informed inductive bias. 2.12 Cross-property physical coupling A central contribution of MultiPUFFIN over prior PUFFIN-family models and over generic molecular foundation models is the use of cross-property physical coupling in the training objective. Where the inductive bias neurons of the previous subsection enforce physically correct behavior within a single property (e.g., monotonic vapor pressure with temperature through the Antoine head), cross-property coupling enforces physically correct relationships between properties that are predicted jointly. Two coupling terms are introduced, each grounded in a well-known thermodynamic identity. The coupling terms do not require additional labels beyond what is already present in the training set; they exploit the fact that multi-task training produces, for every molecule in a batch, simultaneous predictions of multiple thermophysical properties whose physical consistency can be checked against first-principles relationships. The flash point of a pure substance is defined as the lowest temperature at which its vapor pressure reaches the lower flammability limit in air, corresponding to a partial pressure on the order of 10310^3 PaPa for most organic liquids. If MultiPUFFIN predicts both a flash-point value T^FP T_FP and a vapor-pressure curve P^â(T) P(T), these predictions are thermodynamically consistent only if P^â(T^FP)âPFPâ P( T_FP)â P_FP^*, where PFPâP_FP^* is the target flash-point partial pressure. A soft coupling loss penalizes the inconsistency directly: âFPâVP=λFPâ bââŹFPâ[(log10âĄP^â(T^FP,b)âlog10âĄPFPâ)2],L_FP\!-\!VP\;=\; _FP\,·\,E_b _FP\! [ ( _10 P( T_FP,b)\,-\, _10P_FP^* )^2 ], (19) where âŹFPB_FP is the set of molecules in the batch that carry a non-missing flash-point label (ensuring that T^FP T_FP is being actively supervised by a measurement), PFPâ=3000P_FP^*=3000 PaPa is a universal target partial pressure (within a factor of two of the published lower flammability limits across common organic classes), and λFP=10â2 _FP=10^-2 weights the coupling term relative to the main regression loss. Evaluating P^â(T^FP) P( T_FP) requires a second forward pass through the vapor-pressure head with the temperature argument set to T^FP T_FP rather than to the measured batch temperature; gradients from this term flow through both the flash-point head (which is pulled toward values that make the VP prediction at T^FP T_FP consistent with the target) and the vapor-pressure head (which is pulled toward curves that satisfy the definition of flash point). By definition, the normal boiling point TbT_b of a pure substance is the temperature at which its vapor pressure equals the standard atmospheric pressure Patm=101 325P_atm=101\,325 PaPa. For every training molecule that carries a measured boiling-point label, MultiPUFFIN therefore has access to a hard physical anchor for its vapor-pressure head: the VP prediction evaluated at the measured TbT_b must equal log10âĄ(Patm) _10(P_atm). This is enforced by a second coupling term: âBPâanchor=λBPâ bââŹBPâ[(log10âĄP^â(Tb,bmeas)âlog10âĄPatm)2],L_BP\,anchor\;=\; _BP\,·\,E_b _BP\! [ ( _10 P(T_b,b^meas)\,-\, _10P_atm )^2 ], (20) where âŹBPB_BP is the set of molecules in the batch that carry a non-missing boiling-point label, Tb,bmeasT_b,b^meas is the measured boiling point of molecule b, and λBP=5Ă10â2 _BP=5Ă 10^-2. Unlike the flash-point coupling, where both the temperature argument and the target partial pressure are learned or universal constants, the BP anchor uses a measured temperature and a physically exact target pressure. It is thus the stronger of the two coupling terms in a supervisory sense, and the higher coupling weight (5Ă stronger than the FPâVP coupling) reflects this; it has the additional effect of calibrating the vapor-pressure head near atmospheric pressure; a region of high practical importance for process engineering applications. A common concern with auxiliary loss terms is that they can destabilize training if they compete with the main regression loss. In MultiPUFFIN this is mitigated in three ways. First, both coupling weights are set to λ=10â2λ=10^-2, so that at equilibrium each coupling term contributes on the order of 1âââ2%1--2\% of the main multi-task loss. Second, a coupling-loss warmup is applied during the first three epochs of supervised training: only the main multi-task loss is active until the prediction heads have stabilized to reasonable first-pass predictions, after which the coupling terms are switched on. Third, the second forward passes required to evaluate the coupling terms use temperature tensors that are sanitized upstream: any non-finite values (arising e.g. from missing flash-point labels that would otherwise propagate through the VP head) are replaced by a safe placeholder before the forward pass, and the coupling loss itself is then restricted to the masked subset of molecules that carry valid labels. The coupling terms enforce relationships that exist in the physics but do not attempt to enforce full thermodynamic consistency. They do not, for example, enforce Maxwell relations between enthalpy, entropy, and heat capacity; nor do they enforce GibbsâDuhem consistency across a mixture. They are targeted, property-pair-level constraints whose role is to transmit information between a well-measured property (boiling point; flash point) and a less-measured property (vapor pressure curve) through the laws of physics. Fuller thermodynamic-consistency objectives are an important direction for future work (Section 3.10). During development we evaluated two additional cross-property thermodynamic constraints that, although physically motivated, did not transfer cleanly to the heterogeneous training distribution and were therefore not included in the deployed loss. We document them here for completeness and as a methodological caution. The first is a Yalkowsky general-solubility-equation (GSE) constraint linking aqueous solubility to logâĄP P and the melting point through logâĄSâ0.5âlogâĄPâ0.01â(Tmâ25) Sâ 0.5- P-0.01\,(T_m-25). The GSE is well-validated on drug-like molecules near room temperature, but it breaks down on the broad chemistry of the curated dataset; in particular on salts and polyols, where the experimental logâĄS S deviates from the GSE prediction by more than 1 dex. Imposing the GSE as a soft penalty pulled the model off its multi-task optimum on these out-of-GSE classes and degraded mean validation R2R^2 by approximately 0.040.04, so the term was removed. The second is a Troutonâs-rule target on the entropy of vaporization, which would impose a constant-slope penalty on the Antoine head (ÎâSvapâ88 J molâ1 Kâ1 S_vapâ$88\,J\,mol^-1\,K^-1$ at the normal boiling point). Real liquids deviate from Troutonâs value by ±20%± 20\% across hydrocarbons, polar organics, and hydrogen-bonded systems, which makes the constraint a constant residual the model cannot satisfy in any direction. Imposing it slowed convergence without changing the validation R2, and was removed. Both negative results sharpen rather than weaken the case for the cross-property coupling terms that were kept (Eqs. 19, 20): physical constraints help when they hold approximately on the actual training distribution, and hurt when they hold only on the textbook subset of it. Selecting which physics to inject is therefore not a property of physics alone but of the joint distribution physics-and-data. 2.13 Training strategy The training of MultiPUFFIN follows a four-stage protocol: (i) self-supervised pretraining of the multimodal backbone on an unlabeled molecular corpus (Section 2.2); (i) joint supervised multi-task training of the backbone, the condition module stack, and four candidate heads per property; (i) targeted fine tuning with the backbone unfrozen at very low learning rate on an augmented dataset to absorb additional chemical diversity without catastrophic forgetting; and (iv) per-property applicability-domain (G+) evaluation that distinguishes the deployment regime from the global test set. The deployed model is a single forward pass of the four-component pipeline (SSL-pretrained backbone â TConditioner â tournament-winning head per property â G+ scope check); no model-level ensembling is used at inference. The protocol addresses the challenges of multi-task learning with heterogeneous property scales, sparse multi-label data, and the need to balance shared backbone quality with property-specific head calibration. All target property values are z-score normalized using training set statistics: y~p=ypâÎŒpÏp y_p= y_p- _p _p (21) where ÎŒp _p and Ïp _p are the mean and standard deviation of property p computed exclusively over unique molecules in the training set. When SMILES augmentation is employed, normalization statistics are computed on the unique (non-augmented) molecules to prevent artificial variance reduction from duplicate entries. Properties that span orders of magnitude (vapor pressure, viscosity) are first transformed to logarithmic scale (log10 _10) before z-score normalization, ensuring that the loss function treats relative prediction errors uniformly across the dynamic range. The multi-task training objective follows the homoscedastic uncertainty weighting framework of Kendall et al. [26]. For each property p, a learnable scalar sp=logâĄÏp2s_p= _p^2 (initialized to zero) captures the homoscedastic task uncertainty: âtotal=1Ntasksââpâactive[12âexpâĄ(âsp)â MSEp+12âsp]L_total= 1N_tasks _p _active [ 12 (-s_p)·MSE_p+ 12s_p ] (22) where activeP_active is the set of properties with valid (non-NaN) targets in the current batch and the division by NtasksN_tasks normalizes across the number of active tasks. This mechanism automatically down-weights noisy or difficult tasks (increasing sps_p reduces the effective weight expâĄ(âsp)/2 (-s_p)/2 while incurring a regularization penalty sp/2s_p/2), providing a principled alternative to manual task weighting. The sps_p parameters are jointly optimized with the model parameters. Since the multi-label dataset is sparse (each molecule typically has labels for 2â5 of the 9 properties), the loss computation handles missing labels by masking NaN values and computing the loss only over valid targets. To increase the effective training set size and improve the robustness of the SMILES Transformer encoder, SMILES enumeration augmentation is employed. A canonical SMILES string represents a unique, deterministic encoding of a molecule, but the same molecular structure can be represented by many distinct valid SMILES strings that differ in the starting atom, traversal order, and ring-opening positions. For example, toluene can be written as Cc1c1, c1c(C)c1, or C(c1c1), among many others. Because the Transformer encoder processes SMILES as a character sequence, these alternative representations present different syntactic patterns to the model, even though they correspond to the same molecule. By exposing the model to multiple valid SMILES representations during training, the encoder learns to extract molecular features that are invariant to the arbitrary choice of SMILES notation, reducing overfitting to specific syntactic patterns. Concretely, for each molecule naug=2n_aug=2 additional non-canonical SMILES representations are generated using RDKitâs randomized SMILES generation. This triples the effective training set size from 22 01522\,015 to 66 04566\,045 samples while preserving identical graph representations and target values. The augmentation is applied only to the SMILES representation: the molecular graph and three-dimensional conformer remain unchanged for all augmented copies of a given molecule, as these representations are invariant to SMILES notation by construction. The augmentation is applied only during training; validation and test evaluation use only the canonical SMILES to ensure reproducible and unbiased evaluation metrics. The deployed model is produced by a four-stage protocol summarised below; each bullet states the stageâs goal, the parameters trained, and the relevant hyperparameters. âą Stage 1 â self-supervised pretraining. The multimodal backbone (GCN, Transformer, SchNet and cross-modal fusion module) is pretrained on the 500 000500\,000-molecule unlabelled PubChem corpus with the three SSL objectives of Section 2.2 (graphâSMILES contrastive InfoNCE, masked atom-feature reconstruction, masked SMILES token prediction) for 89 epochs (paused before the 100-epoch budget after the validation-loss plateau). Only the backbone weights are retained at the end of Stage 1; the SSL auxiliary classifier heads are discarded. No property labels are used in this stage. âą Stage 2 â joint supervised multi-task training. The pretrained backbone is loaded as initialisation and all parameters (encoders, fusion module, the five condition modules, the four candidate heads per property, uncertainty log-variance parameters, and the auxiliary and descriptor encoders) are jointly optimised on the labelled training set. AdamW [37] is used with ÎČ1=0.9 _1=0.9, ÎČ2=0.999 _2=0.999, weight decay λ=5Ă10â4λ=5Ă 10^-4, base learning rate α0=1Ă10â4 _0=1Ă 10^-4, linear warmup over Ewarm=10E_warm=10 epochs, and cosine annealing with warm restarts thereafter (Section 2.13); gradient accumulation over 2 steps yields an effective batch size of 96 (physical batch 48), and gradient norms are clipped to 1.01.0. Stage 2 trains for up to 120 epochs with early stopping on the multi-task validation loss (patience 25). The main loss is Eq. 22; the cross-property coupling terms (Section 2.12) activate after a three-epoch warmup. The per-property tournament winner is selected at the end of Stage 2 as the head with the lowest validation RMSE among the four candidates A (physics equation), B (Joback group contribution), C (RDKit fragment counts), and D (direct or alternative-physics). âą Stage 3 â backbone-unfrozen targeted fine tune. The Stage-2 checkpoint is warm-restarted on an augmented training set that extends per-property coverage in regions identified as scarce by the coverage audit (Section 2.1): NIST ThermoML [25] contributes additional vapor-pressure, viscosity, heat-capacity and melting-point measurements; the EPA OPERA dataset contributes boiling-point and vapor-pressure data; the Bradley Open Melting Point dataset [6] contributes additional melting points; the Sun et al. flash-point compilation contributes flash-point data; the Chew et al. multi-temperature viscosity dataset contributes per-temperature viscosity measurements; and the AqSolDB salts subset, previously dropped by the multi-component-SMILES filter, is retained and labelled for the solubility head. The training set grows from 22 01522\,015 rows (Stage 2) to roughly 40 00040\,000 rows (Stage 3), with the largest relative additions on viscosity (+84%+84\%), flash point (+56%+56\%), vapor pressure (+67%+67\%) and the salt slice of solubility (+11%+11\%). All parameters are trainable but at very low learning rates (αBCD=5Ă10â6 _BCD=5Ă 10^-6 for the B/C/D heads and the condition modules, αA=1Ă10â6 _A=1Ă 10^-6 for the physics-A heads); gradient norms are clipped to 1.01.0, early stopping patience is 8 epochs, and a hard rollback reverts to the last good checkpoint if mean validation R2R^2 regresses by more than 0.0050.005 between consecutive epochs. The validation/test splits are unchanged across the four stages. âą Stage 4 â per-property applicability-domain (G+) evaluation. The deployed model is evaluated on the held-out test set under the in-scope subset defined by per-property G+ filters (Section 2.14): pre-registered chemical-class filters defined from molecular features only, without consulting any model prediction. The filters are property-specific because the underlying coverage scarcity is property-specific, and they are committed before any test prediction is computed. In-scope R2R^2 is the deployment headline reported throughout the results (Table 7). A final ingredient in the Stage 2 supervised objective is an auxiliary temperature-smoothness loss that targets the temperature-dependent properties (vapor pressure, viscosity, heat capacity, and solubility through its van ât Hoff extension). The motivation is that the thermodynamically-informed heads (Antoine, Andrade, Shomate, van ât Hoff) enforce an analytically correct temperature dependence by construction, while the Joback, RDKit, and direct candidate heads depend on temperature only through a scalar input feature and can in principle produce temperature curves that oscillate or disagree with the underlying thermodynamics. The T-smoothness term pulls the non-physics heads toward the temperature slope predicted by the physics head of the same property, without forcing them to match the physics headâs absolute values. Concretely, for each training sample with a label for a temperature-dependent property at its measured temperature T, a perturbed temperature TâČ=T+ÎŽâT =T+ÎŽ T is drawn with ÎŽâTâŒâ[â20,+20]ÎŽ T [-20,+20] KK (clipped to a physically reasonable range of [150,900][150,900] KK). A second forward pass is performed through the prediction heads at TâČT , keeping the molecular embedding fixed. The physics head A is evaluated at both T and TâČT and the expected slope ÎA=y^Aâ(TâČ)ây^Aâ(T) ^A= y_A(T )- y_A(T) is computed and detached from the autograd graph (it provides the supervisory signal only). The non-physics heads (B, C, D) are evaluated at the same T and TâČT , and the observed slope Îh=y^hâ(TâČ)ây^hâ(T) ^h= y_h(T )- y_h(T) is compared against the detached physics slope through a mean-squared-error penalty summed over the non-physics heads and over the T-dep properties: âTâsmooth=λsmoothâ 1||â 3ââpââhâB,C,Dbâ[(Îp,bhâsgâ(Îp,bA))2],L_T\!-\!smooth\;=\; _smooth\,·\, 1|T|· 3 _p _hâ\B,C,D\E_b\! [ ( ^h_p,b-sg( ^A_p,b) )^2 ], (23) where T is the set of temperature-dependent properties, sgâ(â )sg(·) denotes the stop-gradient operator, and λsmooth=5Ă10â3 _smooth=5Ă 10^-3. The term is added to the Stage 2 loss with the same three-epoch warmup as the cross-property coupling terms. Because sgâ(ÎA)sg( ^A) is detached, the physics head is not itself pulled by this term; only the non-physics heads are regularized toward physically consistent temperature dependence, which preserves the role of the tournament (Section 2.11) as a fair architectural comparison. In the final deployed configuration the T-smoothness term is set to λsmooth=0 _smooth=0 (i.e., disabled). Empirically, the condition module stack of Section 2.9 provides a stronger and more flexible mechanism for T-dependent embedding refinement than the smoothness regularizer; the smoothness term, while consistent with the physics, was found to slightly slow convergence and add no validation-R2R^2 gain once the conditioner was in place. We retain its description here for completeness and as one of the negative-result ablations documented in Section 2.12. The learning rate schedule during Stage 1 follows cosine annealing with warm restarts (SGDR) [36], a schedule specifically designed for non-convex optimization with multiple local optima. The learning rate at epoch t (after warmup) is given by: αâ(t)=αmin+12â(α0âαmin)â(1+cosâĄ(TcurTiâÏ))α(t)= _ + 12( _0- _ ) (1+ ( T_curT_iÏ ) ) (24) where TcurT_cur is the number of epochs since the last restart, TiT_i is the current period length, and αmin _ is a minimum learning rate. The restart schedule is parameterized by an initial period T0=30T_0=30 epochs and a multiplication factor Tmult=2T_mult=2, so that successive periods grow geometrically: the first cycle spans the first 30 epochs (including warmup), with a restart at epoch 30 that resets the learning rate to α0 _0. The second cycle spans epochs 30â90 (T1=60T_1=60 epochs), and the third cycle would begin at epoch 90 (T2=120T_2=120 epochs). The cosine warm-restart schedule provides three specific benefits for multi-task molecular property prediction: 1. Escape from local optima. The multi-task loss landscape (Eq. 22) involves nine coupled objectives with different curvatures, creating a complex Pareto front with many local optima. At the end of each cosine cycle, the learning rate approaches its minimum and the model settles into a local minimum. The restart then abruptly increases the learning rate back to α0 _0, providing sufficient gradient magnitude to escape the current basin and explore alternative regions of the loss landscape that may offer better trade-offs across tasks. 2. Annealing toward increasingly refined minima. The increasing period lengths (T0,2âT0,4âT0,âŠT_0,2T_0,4T_0,âŠ) allow progressively longer annealing phases in later cycles, enabling the optimizer to settle into finer-grained minima as the overall loss landscape becomes better conditioned through training. This is analogous to simulated annealing with a decreasing cooling rate. 3. Implicit snapshot ensembling. The cosine schedule naturally produces model states at different points on the Pareto front of task trade-offs. By retaining the checkpoint with the lowest validation loss across all cycles, the training procedure implicitly selects the best snapshot from this diverse set of solutions. In the final model, the best validation loss was achieved at epoch 33, shortly after the restart at epoch 30, confirming that the learning rate reset allowed the optimizer to escape the first cycleâs local minimum and find a superior solution. The early stopping patience of 50 epochs was chosen to span more than one full cosine cycle (T0=30T_0=30), ensuring that the model has the opportunity to recover from at least one restart before early stopping is triggered. The model is implemented in PyTorch [45] with PyTorch Geometric [16] for graph neural network operations. Training is performed on a single NVIDIA RTX 5090 Laptop GPU (24 GB24\,GB VRAM). All computations use float32 precision; mixed-precision training (float16) was not employed due to numerical instabilities in the SchNet encoderâs Gaussian-smeared distance computations, which are sensitive to float16 underflow in the exponential operations. The total model comprises approximately 35 million parameters in the backbone plus an additional 6â7 million parameters across the nine per-property tournament heads (four candidate heads per property, see Table 5). Training proceeds in four stages as described in Section 2.13: Stage 1 (self-supervised pretraining of the backbone on 500,000 unlabeled PubChem molecules with three SSL objectives, 89 epochs); Stage 2 (joint supervised multi-task training of backbone, condition module stack, and four candidate heads per property, up to 120 epochs with early stopping, cross-property coupling losses activated after a three-epoch warmup); Stage 3 (backbone-unfrozen targeted fine tune at very low learning rate; 5Ă10â65Ă 10^-6 for B/C/D heads and the conditioner, 1Ă10â61Ă 10^-6 for the physics A heads; on the augmented training set with hard rollback on validation regression); and Stage 4 (per-property applicability-domain G+ evaluation, Section 2.14). All four stages use identical val/test splits; no test-set information is used at any stage. Model performance is evaluated using three standard regression metrics computed on the held-out test set: the root mean squared error (RMSE), the mean absolute error (MAE), and the coefficient of determination (R2R^2). All metrics are computed on the original (un-normalized) scale after inverse z-score transformation. For temperature-independent properties, only the evaluation-masked subset (one measurement per molecule) is used. 2.14 Per-property applicability domain (G+) A central question for any foundation model in deployment is the regime over which its predictions can be trusted. MultiPUFFIN addresses this explicitly through a per-property applicability-domain construction that we refer to as the G+ scope. The construction has three properties that distinguish it from a generic train/test split: (i) the scope rules are specified per property, because the underlying coverage scarcity is itself per-property; (i) each rule is computed from molecular features alone (RDKit-derived structural and physicochemical descriptors), without consulting the modelâs predictions, ensuring that the scope is committed before any test inference; and (i) the rules are pre-registered, frozen before the test set is ever scored. A coverage audit was performed on the curated training set, partitioning each propertyâs molecules into a fixed taxonomy of structural classes (salts, charged ions, sulfonic/phosphonic acids, polyols, very-large molecules, very-small molecules, zwitterions). For every (property, class) cell the audit reports the number of training and test entries. Where targeted external data acquisition could close a gap (e.g., AqSolDB salts for solubility), the gap was filled by augmentation. Where no public condition-matched data was available within the scope of this work (e.g., MNSol/DISSOLVE polyol entries for hydration free energy require institutional license access), the affected molecule class was declared out-of-scope for that specific propertyâs evaluation. The G+ filters retained for the deployed model are summarized below. Each filter is a binary RDKit-evaluable predicate over the test molecule. âą Solubility (logâĄS S): salts (multi-component SMILES), charged ions (|q|>0|q|>0), sulfonic/phosphonic acids, very-large molecules (>>50 heavy atoms). âą logâĄP P: zwitterions (mixed positive and negative formal charges), very-large molecules, sulfonic/phosphonic acids. âą Hydration free energy: polyols (small high-HBD molecules with 0 aromatic rings). âą Boiling point: salts, charged ions, sulfonic/phosphonic acids. âą Vapor pressure: salts, charged ions, sulfonic/phosphonic acids, polyols. âą Viscosity: salts, polyols, very-small molecules (â€4†4 heavy atoms). âą Melting point: very-large molecules. âą Flash point: charged ions. âą Heat capacity: no per-property scope filter. The filters do not redefine âwhat MultiPUFFIN predictsâ; the model still produces all nine outputs for any input. They define the evaluation regime in which the headline deployment metrics are reported, so that the published R2R^2 corresponds to the molecular subspace the model is intended to be used on. All reported test-set metrics in this paper are computed on the in-scope subset defined by these per-property G+ filters. The deployment headline is the in-scope mean test R2=0.784R^2=0.784 across the nine properties (Section 3, Table 7). 2.15 Model architecture summary Table 5 provides a complete specification of the MultiPUFFIN architecture, including the structure, dimensions, and parameter counts for each component. The model comprises approximately 35 million parameters, with the SMILES Transformer encoder accounting for the largest share (55.7%55.7\%) due to its 6-layer, 512-dimensional architecture with 8-head self-attention and 2048-dimensional feed-forward sublayers. The cross-modal fusion module is the second largest component (30.0%30.0\%), reflecting the complexity of the bidirectional cross-attention and gated fusion mechanism. The GCN encoder, despite being architecturally simpler, requires only 1.8%1.8\% of the total parameters. The nine prediction heads collectively account for 8.5%8.5\% of the parameters, with the two enhanced-capacity heads for vapor pressure and boiling point being approximately 4Ă4Ă larger than the standard heads. Table 5: Complete architecture specification of MultiPUFFIN. The model comprises 34 964 71534\,964\,715 parameters across five encoder branches (three structural and two auxiliary), a cross-modal fusion module, and nine property-specific prediction heads. Enhanced heads (â ) use 3 hidden layers with 512 units, LayerNorm, GELU activation, and residual connections; standard heads use 2 hidden layers with 256 units and ReLU activation. Component Structure Key dimensions Activation Parameters Encoder branches GCN encoder 4 GCNConv layers + sum pool + projection 39â256â256â51239â 256â 256â 512 ReLU 617 216617\,216 SMILES Transformer 6-layer encoder + CLS pooling + projection dmodel=512d_model=512, h=8h=8, dff=2048d_f=2048 GELU 19 466 24019\,466\,240 SchNet 3D encoder 4 interaction blocks + sum pool + projection d=256d=256, K=50K=50 Gaussians, rc=10 Ă r_c=$10\, $ SiLU ⌠1 460 0001\,460\,000a Experimental encoder 2-layer MLP 6â64â5126â 64â 512 ReLU 38 40038\,400 Descriptor encoder 2-layer MLP 13â128â51213â 128â 512 ReLU 84 86484\,864 Fusion module Cross-modal attention Bidirectional (GCN â Transformer) h=8h=8, d=512d=512, expansion 4Ă4Ă GELU 10 506 24210\,506\,242 Gated fusion + geometry gate Sigmoid gates + projection to d=512d=512 2Ă512â5122Ă 512â 512 Sigmoid Condition-aware refinement stack (Sec. 2.9) TConditioner (T) 2-token self-attn + FFN, identity-init h=4h=4, dff=512d_f=512, dropout 0.1 GELU ⌠1 283 0001\,283\,000 PHConditioner (pH) 2-token self-attn + FFN, identity-init h=4h=4, dff=512d_f=512, dropout 0.1 GELU ⌠1 284 0001\,284\,000 PressureConditioner (P) 2-token self-attn + FFN, identity-init h=4h=4, dff=512d_f=512, dropout 0.1 GELU ⌠1 284 0001\,284\,000 PolymorphEmbedding 2-token self-attn + FFN, identity-init K=16K=16 slots, dff=512d_f=512 GELU ⌠994 000994\,000 MethodFlag 2-token self-attn + FFN, identity-init K=4K=4 slots, dff=512d_f=512 GELU ⌠989 000989\,000 Prediction heads Vapor pressureâ Enhanced Antoine: log10âĄP=AâB/(T+C) _10P=A-B/(T+C) 512â512â512â3512â 512â 512â 3 params GELU 794 115794\,115 Boiling pointâ Enhanced GroupContrib: Tb=T0+âÎi+cT_b=T_0+ÎŁ _i+c 512â512â512â34512â 512â 512â 34 params GELU 808 482808\,482 Viscosity Andrade: log10âĄÎ·=A+B/(T+C) _10η=A+B/(T+C) 512â256â256â3512â 256â 256â 3 params ReLU 197 891197\,891 Solubility van ât Hoff: logâĄS=logâĄSref+Aâ(1/Trefâ1/T) S= S_ref+A(1/T_ref-1/T) 512â256â256â2512â 256â 256â 2 params ReLU 197 634197\,634 logâĄP P DirectHead (FFNN) 512â256â256â1512â 256â 256â 1 ReLU 197 377197\,377 Hydration free energy Born: ÎâG=â(1â1/Δ)â q2/(8âÏâÏ”0âr) G=-(1-1/ )· q^2/(8Ï _0r) 512â256â256â2512â 256â 256â 2 params ReLU 197 634197\,634 Melting point DirectHead (FFNN) 512â256â256â1512â 256â 256â 1 ReLU 197 377197\,377 Flash point DirectHead (FFNN) 512â256â256â1512â 256â 256â 1 ReLU 197 377197\,377 Heat capacity Shomate: Cp=A+Bât+Cât2+Dât3+E/t2C_p=A+Bt+Ct^2+Dt^3+E/t^2 512â256â256â5512â 256â 256â 5 params ReLU 198 401198\,401 Total 34 964 71534\,964\,715 a SchNet parameters estimated; exact count varies with implementation details. Table 6 summarizes the optimized hyperparameters for the three trainable stages of the protocol (Stage 4 is evaluation-only and has no learnable parameters). Table 6: Optimized training hyperparameters for MultiPUFFIN. Stage 1 is self-supervised pretraining of the multimodal backbone; Stage 2 is joint multi-task supervised training of the backbone, the condition module stack, and four candidate heads per property; Stage 3 is the targeted backbone-unfrozen fine tune on the augmented dataset. Hyperparameter Stage 1 (SSL) Stage 2 (joint) Stage 3 (fine tune) Optimizer AdamW AdamW AdamW Learning rate α (B/C/D heads + cond.) 5Ă10â45Ă 10^-4 1Ă10â41Ă 10^-4 5Ă10â65Ă 10^-6 Learning rate αA _A (physics A heads) ; 1Ă10â41Ă 10^-4 1Ă10â61Ă 10^-6 LR schedule Cosine decay Cosine warm restarts Constant Warmup epochs 2 10 0 Cosine T0T_0 / TmultT_mult ; 30 / 2 ; Weight decay λ 1Ă10â51Ă 10^-5 5Ă10â45Ă 10^-4 5Ă10â45Ă 10^-4 Batch size (physical / effective) 256 / 256 48 / 96 48 / 48 Gradient clipping norm 1.0 1.0 1.0 Max epochs 100 (paused at 89) 120 12 Early stopping patience ; 25 8 Hard rollback rule ; ; val mean R2R^2 regression >0.005>0.005 3 Results and Discussion 3.1 Overall predictive performance Table 7 summarises the deployment in-scope test-set performance of MultiPUFFIN across the nine properties after the per-property G+ applicability-domain filters of Section 2.14. The headline is mean test R2=0.784R^2=0.784 for the deployed four-stage model; the architectural-ablation experiments of Section 3.5 use the Stage-2 reference configuration (mean R2=0.708R^2=0.708â0.7160.716) because re-running Stage 3+4 for every variant would multiply the compute by an order of magnitude without changing the relative component comparison, so the ⌠0.070.07 gap is the uniform Stage-3+4 lift that the ablations omit. The lowest errors are on hydration free energy (RMSE 0.690.69 kcal/mol, R2=0.949R^2=0.949), heat capacity (18.418.4 J/mol/K, 0.9400.940), and the lipophilicity / process-engineering properties (logâĄP P 0.8260.826, flash point 0.8160.816, TbT_b 0.7960.796, TmT_m 0.7800.780, viscosity 0.7780.778). The two cases that are limited by data-distribution heterogeneity rather than by the architecture are vapor pressure (0.6240.624) and aqueous solubility (0.5520.552), the latter especially constrained by inter-laboratory disagreement and pH-dependent variability for ionisable solutes below. Table 7: Test set performance of MultiPUFFIN across nine physicochemical properties on the deployment in-scope subset, defined by the per-property G+ applicability-domain filters of Section 2.14. n is the number of in-scope test molecules; RMSE and MAE are reported in the units of each property. Property (units) n RMSE MAE R2R^2 Hydration free energy (k/âââk /mol) 267 0.689 0.380 0.949 Heat capacity (J molâ1 Kâ1J\,mol^-1\,K^-1) 297 18.39 7.93 0.940 logâĄP P (â) 13901390 0.902 0.541 0.826 Flash point (KK) 28262826 25.55 16.10 0.816 Boiling point (KK) 20752075 34.29 20.43 0.796 Melting point (KK) 62026202 36.66 26.06 0.780 Viscosity (log10 _10 mPa smPa\,s) 187 0.265 0.166 0.778 Vapor pressure (log10 _10 PaPa) 16101610 1.674 1.144 0.624 Solubility (log mol lâ1mol\,l^-1) 25432543 1.153 0.823 0.552 Mean 0.784 Figure 2 compares train/val/test RMSE and MAE across the nine properties on a log scale. The variation in accuracy reflects three interacting factors: the intrinsic complexity of each structureâproperty mapping, the amount and quality of training data, and the suitability of the chosen domain equation. Heat capacity is data-clean and varies smoothly with structure; HFE is FreeSolv-grade despite the modest 642 molecules and benefits from positive transfer through the shared backbone. Figure 2: Test RMSE (left) and MAE (right) across all nine properties for training, validation, and test splits (logarithmic scale). The consistent increase from training to test error reflects the generalization challenge imposed by the scaffold-based splitting strategy. Properties are ordered by decreasing test performance. Note the logarithmic y-axis: temperature-based properties (e.g., melting point, boiling point) have RMSE in the tens of kelvins, while logarithmic-scale properties (e.g., viscosity, vapor pressure) have RMSE below 2. Note: as for the parity plots, the per-property bar heights in this figure correspond to the Stage-2 reference run; deployed four-stage model values are tabulated in Table 7. 3.2 Per-property analysis Figure 3 shows predicted-versus-experimental parity plots for all nine properties on the test set; the diagonal line is perfect prediction. Figure 3: Test set parity plots (predicted vs. experimental) for all nine physicochemical properties. The solid diagonal line represents perfect prediction (y=xy=x). Each panel reports the R2R^2, RMSE, MAE, and number of test samples. Properties are ordered by decreasing test R2R^2. Note: per-panel metrics in this figure are drawn from the controlled Stage-2 reference run used as the architectural-ablation baseline (Section 3.5); the deployed four-stage model values are reported in Table 7 and define the R2=0.784R^2=0.784 headline. The per-panel pattern in Figure 3 reads cleanly as three groups. The first group is the diagonal-tight regime: heat capacity (panel a, RMSE = 18.4 J molâ1 Kâ118.4\,J\,mol^-1\,K^-1) and hydration free energy (panel b, RMSE = 0.69 k/âââ0.69\,k /mol, R2=0.949R^2=0.949 on 267 in-scope molecules) cluster along the y=xy=x line without systematic bias. Both benefit from a tight match between the chosen domain-informed head (Shomate polynomial; Born functional form repurposed as a flexible charged-cavity fit) and the underlying physics. The second group is the moderate-scatter regime: logâĄP P (panel c, RMSE = 0.9020.902, R2=0.826R^2=0.826), melting point (panel d, RMSE = 36.7 K36.7\,K, broader scatter at high TmT_m where crystal packing dominates [24]), viscosity (panel e, RMSE = 0.2650.265 in log10 _10 units ; a factor of 1.841.84 in absolute η), flash point (panel f, RMSE = 25.5 K25.5\,K, with the expected widening above ⌠450 K450\,K where vapourâcomposition coupling at the surface starts to dominate). The third group is the noise-limited regime: solubility (panel g, RMSE = 1.1531.153 in log units, broad scatter driven by highly insoluble compounds and inter-laboratory disagreement across AqSolDB), boiling point (panel h, RMSE = 34.3 K34.3\,K) and vapor pressure (panel i, RMSE = 1.6741.674 in log10 _10 units across eleven decades of raw PvapP_vap). The enhanced-capacity heads give VP and BP additional parameter budget, but in this regime the residual error is dominated by the heterogeneity of the multi-source corpus rather than by representational capacity. Figure 4 complements the parity plots by presenting the residual distributions for all nine properties on the test set. In each panel, a Gaussian fit is overlaid on the histogram, and the mean (ÎŒ) and standard deviation (Ï) of the residuals are reported. The residual distributions are approximately centered at zero for all properties, confirming the absence of systematic bias. The residuals for heat capacity and hydration free energy closely follow a normal distribution, indicating well-calibrated predictions with homoscedastic error. Properties such as boiling point and vapor pressure exhibit heavier tails, reflecting the presence of outlier molecules for which the modelâs predictions deviate substantially, likely compounds with unusual structural features underrepresented in the training set. Figure 4: Test set residual distributions for all nine properties. Each panel shows the histogram of residuals (predicted â- experimental) with a Gaussian fit (solid line) and a zero-residual reference line (dashed red). The mean (ÎŒ) and standard deviation (Ï) of the residuals are annotated. 3.3 Training dynamics and convergence behavior Figure 5 shows the Stage-2 training and validation loss curves (panel a) and the per-property validation loss evolution (panel b). Figure 5: Training dynamics of MultiPUFFIN during Stage 2 of the four-stage protocol (Section 2.13). (a) Overall training and validation loss across the joint multi-task phase and its head-fine tuning epilogue, with the best validation loss marked (â ). The vertical dashed line indicates the transition from joint training (with all parameters trainable) to a backbone-frozen head-consolidation phase that closes Stage 2 before the Stage 3 backbone-unfrozen targeted fine tune (not shown here). (b) Per-property validation loss evolution, showing different convergence rates for different properties. SSL pretraining (Stage 1, 89 epochs on 500,000 unlabeled molecules) and the Stage 3 fine tune are not shown in this panel. The joint phase of Stage 2 (cosine warm restarts) ran for 83 epochs before early stopping, with best validation loss 0.3880.388 at epoch 33; the warm-restart boundary at T0=30T_0=30 epochs is visible as a discontinuity in the loss curve. Per-property convergence is heterogeneous: heat capacity and boiling point converge within ten epochs, while logâĄP P and viscosity require more, consistent with the differing complexity of their structureâproperty mappings. The head-consolidation epilogue (backbone-frozen, α=2Ă10â5α=2Ă 10^-5, 26 epochs) reaches its best validation loss of 0.3810.381 in the first epoch, confirming that the gains there come from isolated head re-calibration against the now-frozen backbone. Stage 3 (the backbone-unfrozen targeted fine tune on the augmented set) and Stage 4 (per-property G+ evaluation) are described in Section 2.13 and produce the deployed-model performance of Table 7. 3.4 Multi-task learning trade-offs and capacity dilution Validation-set performance is reported alongside the test set in Table 7 and follows the expected train<val<testtrain<val<test ordering on every property (validation R2R^2 in the range 0.650.65â0.960.96 across the nine targets); the trainâtest generalisation gap is largest precisely on the scaffold-split properties whose test set imposes the greatest structural novelty. Figure 6 shows the highly heterogeneous split sizes that drive this regime, ranging from 297 test data rows for heat capacity to 62066206 for melting point in the deployed four-stage model (Table 7); the in-panel counts of the figure are the Stage-2 splits before the Stage-3 augmentation step, so the deployed-model test counts are read from Table 7. Figure 6: Data availability per property and split. The highly heterogeneous sample counts across properties reflect the different availability of experimental measurements in public databases. Training counts reflect unique molecules (before SMILES augmentation). The multi-task setting provides clear benefits for data-scarce properties through positive transfer from larger property datasets. Heat capacity (861 training molecules), viscosity (2024 molecules), and hydration free energy (642 molecules) all achieve stronger performance than would be expected from their training set sizes alone, confirming that the shared backbone extracts generalizable molecular features from the more data-rich properties. Conversely, data-rich properties (boiling point, vapor pressure) show moderate performance regression in the nine-property setting relative to single-task baselines. Figure 7 quantifies the generalization gap (ÎâRMSE=RMSEtestâRMSEtrain =RMSE_test-RMSE_train) for each property, providing a direct measure of how much performance degrades when moving from the training distribution to structurally novel test molecules. Because properties have different units and scales, the normalized RMSE (NRMSE = RMSE/Ïexp _exp) is used for cross-property comparison. The largest generalization gaps are observed for properties that use scaffold-based splitting and thus face the greatest structural novelty in the test set (boiling point, flash point, solubility). In contrast, properties with coverage-balanced splitting (heat capacity, viscosity, hydration free energy) exhibit much smaller gaps, confirming that the model generalizes well when the test set does not impose strict structural novelty. Figure 7: Generalization gap (normalized Î : test NRMSE minus train NRMSE) for each property. Larger gaps indicate greater difficulty in generalizing to structurally novel test molecules, correlating with the scaffold-based splitting strategy used for common properties. This capacity dilution phenomenon is well documented in multi-task learning [11]. The shared backbone, constrained to produce a single 512-dimensional embedding that must serve nine different prediction objectives, cannot allocate unlimited capacity to any single task. The enhanced-capacity heads partially mitigate this by providing properties most affected by dilution with approximately 4Ă4Ă more dedicated parameters (⌠800 000800\,000 vs. ⌠197 000197\,000). However, the fundamental bottleneck lies in the shared representation, not in the head capacity. Figure 8 presents the normalized error metrics (NRMSE = RMSE/Ïexp _exp and NMAE = MAE/Ïexp _exp, where Ïexp _exp is the standard deviation of the experimental values in the test set) across all properties and splits. Normalizing by the experimental standard deviation places all properties on a common, dimensionless scale and allows meaningful cross-property comparison: a value of 1.0 would correspond to simply predicting the mean. Heat capacity achieves the lowest normalized errors (NRMSE = 0.20, NMAE = 0.12), consistent with its high R2R^2. The most challenging properties (solubility, boiling point, vapor pressure) cluster around NRMSE â 0.62, indicating that the model explains substantially more variance than a mean-only baseline. The consistent ordering of train << validation << test normalized errors across all properties confirms the expected generalization hierarchy, with the gap being smallest for properties with coverage-balanced splitting. Figure 8: Normalized RMSE (a) and normalized MAE (b) across all nine properties for training, validation, and test splits. Metrics are normalized by the standard deviation of experimental values in the test set (Ïexp _exp), enabling meaningful cross-property comparison. The dashed line at 1.0 represents the error of a mean-only predictor. 3.5 Ablation studies Ablations isolate the contribution of each architectural component and prediction head under a controlled 100-epoch budget (early-stopping patience 20). Each variant modifies exactly one aspect of the full architecture while keeping the data, seed, normalisation and budget fixed. The ablation reference is the Stage-2 frozen-backbone configuration (mean R2=0.708R^2=0.708â0.7160.716), not the deployed four-stage model whose mean R2=0.784R^2=0.784 headline includes the Stage-3+4 lift; rerunning Stage 3+4 for every variant would multiply the compute by an order of magnitude without changing the relative contribution that the ablation is designed to isolate. The absolute RMSE values in Tables 8 and 9 therefore lower-bound the deployed-model contribution of each component, with the Stage-3+4 lift adding uniformly on top. 3.5.1 Contribution of architectural components Three architectural variants were trained to assess the contribution of the Transformer encoder and the SchNet 3D geometry encoder (Table 8, Figure 9): âą GCN-Only: removes both the Transformer and SchNet encoders, using only the 3-layer GCN for molecular encoding. This represents the simplest baseline architecture. âą No SchNet: removes only the SchNet 3D geometry encoder, retaining the GCN and Transformer branches. This isolates the contribution of explicit 3D conformer information. âą Full MultiPUFFIN: the complete trimodal architecture with GCN, Transformer, and SchNet (reference). The GCN-Only ablation degrades performance across most properties, with the largest hits on data-scarce hydration free energy (ÎâRMSE=+1.07 =+1.07 kcal/mol) and heat capacity (+9.29+9.29 J/mol/K) ; the properties that gain most from the richer representation. The exception is logâĄP P, where GCN-Only outperforms the full model (RMSE 0.7140.714 vs. 0.9060.906), reflecting that octanolâwater partition is dominated by 2D substructure patterns that the extra modalities mildly dilute. Removing only SchNet (No-SchNet) hits HFE hardest (Î=+0.90 =+0.90 kcal/mol), as expected for a property dominated by 3D solvation cavity geometry, and also affects CpC_p and logâĄP P. Table 8: Architectural ablation: test RMSE for the full MultiPUFFIN model versus simplified encoder configurations. Note: All ablation models (including the Full Model reference) were trained with a standardized 100-epoch budget and early stopping patience of 20 epochs, which differs from the full four-stage training protocol used for the final model reported in Table 7. Bold values indicate the lowest RMSE (best performance) for each property; Î denotes the change from the full model (positive = degradation). Units match Table 7. Property Full Model GCN-Only No SchNet Heat capacity (J/mol/K) 28.81 38.10 (++9.29) 40.47 (++11.65) HFE (kcal/mol) 1.103 2.176 (++1.07) 2.005 (++0.90) logâĄP P 0.906 0.714 (â-0.19) 1.109 (++0.20) Viscosity (log10 _10) 0.249 0.304 (++0.055) 0.296 (++0.048) Melting point (K) 41.91 45.94 (++4.03) 46.30 (++4.39) Flash point (K) 34.58 34.14 (â-0.44) 33.40 (â-1.18) Boiling point (K) 52.95 53.02 (++0.06) 52.77 (â-0.18) Vapor pressure (log10 _10) 1.972 1.918 (â-0.053) 1.926 (â-0.045) Solubility (log ) 1.477 1.639 (++0.16) 1.588 (++0.11) Figure 9: Architectural ablation study: per-property test RMSE for the full trimodal MultiPUFFIN architecture versus the GCN-Only and No SchNet variants. Removing the SchNet 3D encoder disproportionately increases RMSE for geometry-sensitive properties (HFE, heat capacity), while removing the Transformer and SchNet structural encoders (GCN-Only) degrades performance across most properties. 3.5.2 Influence of domain-informed prediction heads To evaluate the contribution of domain-informed inductive bias, three head ablation variants were trained (Table 9, Figure 10): âą All DirectHeads: replaces all domain-informed prediction heads with purely data-driven feed-forward networks (DirectHead), removing all domain equations. The three properties that already use DirectHead (logâĄP P, melting point, flash point) are unaffected. This variant represents the fully data-driven baseline with no thermodynamic structure. âą Swapped Antoineâ : assigns the Antoine equation to viscosity and the Andrade equation to vapor pressure, the reverse of the physically correct assignment. This tests whether the specific equation-property pairing matters or whether any structured output provides equivalent regularization. âą All GroupContribution: replaces all property-specific equation heads with the GroupContributionHead, an alternative domain-informed representation rooted in thermodynamic group additivity principles [43]. Unlike the property-specific equations (Wagner, Andrade, van ât Hoff), which encode the known functional form of each propertyâs temperature or state dependence, the group contribution head decomposes the prediction into additive fragment-level contributions, a widely used domain-informed paradigm in chemical engineering (Joback, Lydersen, Nannoolal methods). This ablation tests whether a single, flexible domain-informed representation can serve all properties, or whether property-specific functional forms are necessary. The three head-ablation variants together (Table 9, Figures 11 and 12) tell a coherent story. Replacing all initial property-specific heads with DirectHeads is non-uniform: HFE improves (1.045â0.7711.045â 0.771 kcal/mol) and logâĄP P improves substantially (1.190â0.7541.190â 0.754), confirming that the original LFER and thermodynamic-decomposition heads were too rigid for those properties; viscosity degrades (0.326â0.3540.326â 0.354), confirming that the Andrade head is doing real work. These two findings motivated the equation-level ablation (Section 3.5.3) and the eventual adoption of the Born head for HFE and the DirectHead for logâĄP P in the deployed model. The Swapped Antoineâ variant provides the cleanest evidence that equation-property pairing matters: applying Andrade to vapor pressure raises VP RMSE from 1.9721.972 to 2.8062.806 (+42%+42\%, near-random), while applying Antoine to viscosity also degrades but only by +0.058+0.058 in log10 _10 units ; the asymmetry is physically interpretable, because the Andrade form lacks the three-parameter curvature needed for vapor pressure. The All-GroupContribution variant is the most informative single contrast: it replaces every property-specific equation with a single thermodynamically motivated additive-decomposition head, and the result is that most non-temperature-dependent properties improve (HFE 1.103â0.8031.103â 0.803, TmT_m 41.91â38.6641.91â 38.66, TbT_b 52.95â50.2152.95â 50.21, FP 34.58â33.1934.58â 33.19) while viscosity degrades (0.249â0.2800.249â 0.280). The take-away is that domain-informed structure consistently beats unconstrained heads, but the optimal form is property-specific: Andrade is essential for viscosity, while group contribution is the broadly best thermodynamic prior for Tm,Tb,FP,HFET_m,T_b,FP,HFE. This is the empirical justification for the per-property tournament adopted in the deployed model rather than a single one-size-fits-all head choice. Table 9: Domain-informed head ablation: test RMSE across all nine properties. Note: All ablation models (including the Full Model reference) were trained with a standardized 100-epoch budget and early stopping patience of 20 epochs, which differs from the full four-stage training protocol used for the final model reported in Table 7. Bold values indicate lowest RMSE (best performance) per property. The Swapped variant exchanges Antoineâ assignments; the All GroupContrib variant replaces all thermodynamically-informed heads with generic additive decomposition. Units match Table 7. Property Full Model All Direct Swapped All GroupContrib Heat capacity (J/mol/K) 28.81 40.54 31.84 28.13 HFE (kcal/mol) 1.103 0.740 1.291 0.803 logâĄP P 0.906 0.742 0.835 0.920 Viscosity (log10 _10) 0.249 0.283 0.307 0.280 Melting point (K) 41.91 40.19 42.94 38.66 Flash point (K) 34.58 33.46 34.22 33.19 Boiling point (K) 52.95 50.40 53.62 50.21 Vapor pressure (log10 _10) 1.972 1.939 2.806 1.941 Solubility (log ) 1.477 1.489 1.472 1.454 Figure 10: Influence of domain-informed prediction heads on the six thermodynamically-informed properties. The All DirectHeads variant reveals that domain equations are most beneficial for viscosity (Andrade), while hydration free energy and logâĄP P achieve lower RMSE without the thermodynamically-informed constraints. The Swapped Antoineâ variant demonstrates catastrophic vapor pressure degradation (RMSE increases by 42%42\%), proving that correct equation-property pairing is essential. Figure 11: Per-property test RMSE heatmap across all ablation variants. Green indicates low RMSE (better performance); red indicates high RMSE. The All GroupContribution variant (itself a domain-informed representation based on thermodynamic group additivity) achieves the lowest RMSE for most properties, while the full model provides the best performance for viscosity (Andrade equation). Both domain-informed variants (property-specific equations and group contribution) consistently outperform the purely data-driven All DirectHeads baseline, confirming the value of thermodynamic structure in the prediction heads. Figure 12: Per-component contribution to model performance, expressed as ÎâRMSE=RMSEablationâRMSEfull =RMSE_ablation-RMSE_full (positive values indicate the full model achieves lower error). The most dramatic effect is the catastrophic vapor pressure RMSE increase when the Antoine and Andrade equations are swapped (ÎâRMSE=+0.83 =+0.83). The SchNet 3D encoder contributes most to reducing HFE error (ÎâRMSE=+0.90 =+0.90 kcal/mol), while several properties show negative ÎâRMSE for the All DirectHeads and All GroupContribution variants, indicating that the domain-specific equations are overly restrictive for those properties. 3.5.3 Equation-level ablation: alternative thermodynamic equations A complementary equation-level ablation tested whether the specific thermodynamic equation chosen for each property is optimal, or whether alternative equations from the same domain perform better. The trained backbone was frozen and only the head was replaced and retrained (Stage 2 protocol, α=5Ă10â4α=5Ă 10^-4, patience 15, max 80 epochs); 28 equationâproperty combinations across 12 newly implemented domain-informed heads, the existing baselines and the DirectHead references were evaluated (Table 11). Three regimes appear. Alternative equation improves: HFE drops from RMSE 1.0451.045 kcal/mol with the thermodynamic decomposition head to 0.7040.704 kcal/mol with the Born solvation model (â33%-33\%); vapor pressure improves marginally from Antoine (1.8021.802) to Wagner (1.7891.789); heat capacity improves slightly from DirectHead (21.0321.03 J/mol/K) to Shomate (20.6220.62). Baseline near-optimal: Andrade (0.3260.326) wins viscosity over VFT (0.3300.330) and DirectHead (0.3540.354); van ât Hoff (1.1801.180) wins solubility over Modified Apelblat (1.1971.197); the group-contribution boiling-point head (47.8247.82 K) beats the Nannoolal ratio variant (55.2655.26 K). Equation diverges: the Yalkowsky TmT_m form (division by predicted ÎâSfus S_fus), the NASA 7-coefficient CpC_p polynomial (T4T^4 terms amplify parameter error) and an empirical four-parameter viscosity equation failed to converge under head-only training, indicating that some equations require joint backboneâhead gradients to remain numerically stable. These ablations led to the final assignments of the deployed model: Born for HFE (â33%-33\% RMSE), Wagner for VP (â0.7%-0.7\%), Shomate for CpC_p (â1.9%-1.9\%) and DirectHead for logâĄP P (â0.6%-0.6\% vs. the LFER baseline). The resulting optimal hybrid configuration (Table 10) reaches mean test R2=0.716R^2=0.716 vs. the baseline 0.7080.708, with the gain dominated by the Born switch on HFE (R2:0.892â0.951R^2:0.892â 0.951). Table 10: Optimal hybrid equation assignment vs. baseline MultiPUFFIN. Four of nine prediction heads were replaced with better-performing alternatives identified by the equation-level ablation (Section 3.5.3). Changed heads are indicated with â. All results use the same frozen backbone; only the prediction heads differ. Property Equation RMSE R2R^2 Î ÎâR2 R^2 Solubility (log ) van ât Hoff 1.180 0.539 ; ; logâĄP P DirectHead â 1.182 0.718 â-0.008 ++0.004 HFE (kcal/mol) Born â 0.704 0.951 â-0.341 ++0.059 Boiling point (K) Group contribution 47.82 0.632 ; ; Vapor pressure (log10 _10) Wagner â 1.789 0.587 â-0.013 ++0.006 Viscosity (log10 _10) Andrade 0.326 0.684 ; ; Melting point (K) DirectHead 45.09 0.667 ; ; Flash point (K) DirectHead 30.32 0.740 ; ; Heat capacity (J/mol/K) Shomate â 20.62 0.924 â-0.406 ++0.003 Mean R2R^2 0.716 ++0.008 Table 11: Equation-level ablation: test RMSE for alternative thermodynamic equations on each property, with the backbone frozen and only the prediction head retrained. Bold values indicate the lowest RMSE (best equation) per property. Î denotes the change relative to the baseline equation (negative = improvement). The DirectHead (no physics) serves as a reference for the contribution of domain-informed structure. Three equations diverged during training (â ) and are excluded from the main comparison. NRMSE is computed as RMSE/Ïexp/ _exp. Property Equation RMSE Î NRMSE Parameters Vapor pressure (log10 _10 Pa) Antoine (baseline) 1.802 0.647 792 579792\,579 Wagner (6-param) 1.789 â-0.013 0.643 528 390528\,390 ClausiusâClapeyron 1.814 ++0.012 0.651 526 338526\,338 DirectHead 1.794 â-0.008 0.644 525 825525\,825 Viscosity (log10 _10 mPaâ ·s) Andrade (baseline) 0.326 0.562 197 891197\,891 VFT 0.330 ++0.004 0.569 197 891197\,891 DirectHead 0.354 ++0.028 0.610 197 377197\,377 Solubility (log mol/L) van ât Hoff (baseline) 1.180 0.679 197 634197\,634 Modified Apelblat 1.197 ++0.017 0.689 197 891197\,891 DirectHead 1.194 ++0.014 0.687 197 377197\,377 Boiling point (K) Group Contribution (baseline) 47.82 0.607 808 482808\,482 Nannoolal Ratio 55.26 ++7.44 0.701 535 059535\,059 DirectHead 48.92 ++1.10 0.621 525 825525\,825 logâĄP P (â) LFER (baseline) 1.190 0.535 203 569203\,569 Abraham LSER 1.185 â-0.005 0.533 198 411198\,411 DirectHead 1.182 â-0.008 0.531 197 377197\,377 HFE (kcal/mol) Thermodynamic (baseline) 1.045 0.329 197 634197\,634 Born Model 0.704 â-0.341 0.221 197 891197\,891 DirectHead 0.771 â-0.275 0.242 197 377197\,377 Melting point (K) DirectHead (baseline) 45.09 0.577 197 377197\,377 Yalkowsky Tm=ÎâH/ÎâST_m= H/ S â â â â 197 634197\,634 Flash point (K) DirectHead (baseline) 30.32 0.510 197 377197\,377 SatyanarayanaâRao 30.71 ++0.39 0.516 197 636197\,636 Carroll NFPN_FP 34.30 ++3.98 0.577 197 637197\,637 Heat capacity (J/mol/K) DirectHead (baseline) 21.03 0.281 197 377197\,377 Shomate (NIST) 20.62 â-0.41 0.276 198 405198\,405 3.6 Why multimodal + domain-informed beats single-modality GNN baselines The ablation studies of Section 3.5.1 make explicit why the trimodal foundation-model design beats single-modality GNN baselines, and the gain comes from three jointly-acting mechanisms rather than a single dominant lever. The first mechanism is complementary modality information. The GCNâs L=4L=4 message-passing depth means information propagates at most four bonds from any atom; long-range substituent effects and macrocyclic couplings require either many more layers (with the over-smoothing penalty of [31]) or a global-receptive-field encoder. The Transformer over the SMILES sequence supplies that global receptive field, and the SchNet branch supplies through-space distances and shape that neither the 2D graph nor the SMILES string encodes. The architecture ablation quantifies this: removing the Transformer (GCN-Only) inflates hydration-free-energy RMSE by 97% (1.103â2.1761.103â 2.176 kcal/mol) and heat-capacity RMSE by 32% (28.81â38.1028.81â 38.10 J/mol/K), and removing SchNet alone inflates the HFE RMSE by 82% and the CpC_p RMSE by 40%. The advantage is not uniform across properties ; logâĄP P is in fact slightly better with the GCN alone (RMSE 0.7140.714 vs. 0.9060.906), an honest reminder that lipophilicity is essentially a 2D-substructure phenomenon and that the extra modalities are useful insofar as the property in question requires them. The second mechanism is multi-task positive transfer. A single-task GNN trained only on the 642 FreeSolv molecules cannot learn a rich molecular representation; in MultiPUFFIN that same head sees the same backbone that was already shaped by 12,000 melting-point and 10,000 logâĄP P rows, and the resulting hydration-free-energy RMSE (0.7040.704 kcal/mol) is at the experimental-uncertainty floor. The same effect lifts heat capacity (20.62 J molâ1 Kâ120.62\,J\,mol^-1\,K^-1 on 968968 training molecules) and viscosity (RMSE 0.3260.326 on 20242024 training molecules) into ranges that single-task training on those datasets cannot reach. The third mechanism is intra-property thermodynamic consistency by construction. A vanilla MLP output layer can produce a vapour-pressure curve that decreases with temperature, a viscosity curve that increases with temperature, or Antoine coefficients of unphysical sign ; pathologies that disqualify a model from process-simulator interfacing. Replacing the MLP head with the property-appropriate equation (Antoine, Andrade, van ât Hoff, Born, Shomate) makes the corresponding monotonicity guaranteed, not learned. We are explicit that this is intra-property consistency, not full cross-property consistency: predicted Pvapâ(T)P_vap(T) curves are not strictly ClausiusâClapeyron-coupled to predicted ÎâHvap H_vap, and predicted CpC_p does not integrate exactly to predicted enthalpy changes. Cross-property thermodynamic coupling is a target for future work. The intra-property gain that already exists is empirically verified by the ablations: removing the Andrade head (DirectHead replacement) raises viscosity RMSE from 0.3260.326 to 0.3540.354 (+8.6%+8.6\%), and deliberately mis-assigning Antoine and Andrade across vapor pressure and viscosity raises VP RMSE from 1.9721.972 to 2.8062.806 (+42%+42\%). A practical fourth advantage that is not quantified in the ablation table but matters at deployment time is graceful degradation with missing data. The geometry gate zeros out the SchNet branch when no conformer is available; the per-row multi-task loss back-propagates only through properties with labels for that molecule; and learnt missing-data embeddings let the auxiliary encoders distinguish a measured zero from an absent value. The same trained model therefore handles the heterogeneous, partially-labelled real datasets of Section 2.1 without per-row preprocessing. 3.7 Embedding space analysis To assess whether the learned molecular representation captures meaningful chemical structure, the unified 512-dimensional embedding u (the output of the cross-modal fusion module, before any prediction head) was extracted for all 37 96837\,968 unique molecules across the training, validation, and test sets. This embedding constitutes the shared representation from which all nine property heads operate; its structure therefore reflects the modelâs internal organization of chemical space. The embedding was projected to two dimensions using UMAP [40] (cosine distance, nneighbors=30n_neighbors=30, dmin=0.3d_min=0.3) and subsequently analyzed through both unsupervised clustering and property-based coloring. 3.7.1 Continuous property gradients in embedding space Figure 13 presents the UMAP projection of the embedding space colored by each of the nine predicted property values. The key observation is that all properties exhibit smooth, continuous gradients across the embedding manifold: molecules with similar property values are consistently embedded in nearby regions. Boiling point and flash point display a clear low-to-high gradient from the lower-left to the upper-right region of the manifold (reflecting the strong correlation between these two properties), while log P shows a gradient roughly orthogonal to the thermal properties, separating hydrophilic from hydrophobic molecules. Hydration free energy, which correlates with polarity and molecular surface accessibility, occupies a compact region of the embedding space with the most negative values concentrated in a distinct cluster. These smooth property landscapes indicate that the multi-task training objective has successfully organized the embedding space according to thermophysical similarity, rather than memorizing property values through disconnected mappings. Figure 13: UMAP projection of the MultiPUFFIN unified embedding (ââ512u ^512) for 37 96837\,968 unique molecules, colored by each of the nine predicted properties. Gray points indicate molecules without a measurement for the given property. All properties show smooth, continuous gradients, indicating that the learned representation organizes molecules by thermophysical similarity. Temperature-based properties (boiling point, melting point, flash point) share similar spatial gradients, reflecting their physical correlations. k-Means clustering on the 512-dimensional embeddings (Figure 14, panel a) identifies k=5k=5 as optimal by silhouette analysis (=0.178=0.178, evaluated over kâ5,8,10,12,15,20kâ\5,8,10,12,15,20\), and the five clusters carry chemically interpretable signatures: a polar/hydrophilic cluster (lowest mean logâĄP=â0.26 P=-0.26, most negative ÎâGhyd=â17.2 G_hyd=-17.2 kcal/mol), a small-volatile cluster (lowest mean Tb=427T_b=427 K), a high-boiling polar cluster (Tb=551T_b=551 K, logâĄP=2.5 P=2.5), a hydrophobic cluster (highest logâĄP=5.5 P=5.5), and a residual general-purpose cluster of intermediate values. The moderate silhouette score is itself informative: the embedding is a continuous manifold with gradual chemical transitions rather than a hard partition, which is the appropriate geometry for a multi-task regressor predicting continuous thermophysical targets. HDBSCAN density-based clustering (Figure 14, panel b) corroborates this view, identifying two high-density cores with ⌠10% of molecules in inter-cluster transitions. Colouring the same UMAP by Murcko scaffold (panel c) shows that the 36493649 unique scaffolds form recognisable archipelagos ; a non-trivial outcome given that the training loss contains no explicit structural term ; and the same archipelago view explains why scaffold-split test molecules face an unfamiliar region of the manifold (cf. the generalisation gap of Section 3). The train/val/test overlay (panel d) confirms that all three splits are well-mixed across the manifold, so the scaffold split tests structural novelty without confining the test set to a single chemical subspace. Figure 14: Cluster structure of the MultiPUFFIN embedding space. From left to right: (a) k-Means clusters (k=5k=5, optimal by silhouette analysis), revealing chemically interpretable groups ranging from polar/hydrophilic to hydrophobic; (b) HDBSCAN density-based clustering, identifying two major density cores with 10.2% noise (inter-cluster transitions); (c) top-10 Murcko scaffold families, showing scaffold-level co-localization in embedding space; (d) train/val/test split overlay, confirming that all splits are well-mixed throughout the manifold. 3.8 Comparison with single-property models and literature baselines To contextualize the predictive performance of MultiPUFFIN, we evaluated the model directly on the same molecules used in three widely cited benchmarks: the ESOL dataset [13] for aqueous solubility (1117 molecules), the FreeSolv dataset [42] for hydration free energy (642 molecules), and the MoleculeNet Lipophilicity dataset [63] for logâĄP P (4200 molecules). Because all three benchmark datasets were incorporated into MultiPUFFINâs training corpus during data curation, we report performance separately for molecules assigned to MultiPUFFINâs scaffold-based test split (never seen during training) and for all benchmark molecules (Table 12). Table 12: MultiPUFFIN evaluated on standard benchmark datasets. âTest-onlyâ includes only benchmark molecules assigned to MultiPUFFINâs scaffold-split test set; âAllâ includes all benchmark molecules. Literature baselines use random splitting unless noted. Bold indicates lowest RMSE. Benchmark Method n RMSE ESOL (solubility, log mol/L) Delaney (2004, random) 1128 0.75 GCN (MoleculeNet, random) 1128 0.97 MPNN (MoleculeNet, random) 1128 0.58 MultiPUFFIN (all) 1117 0.649 MultiPUFFIN (test-only, scaffold) 177 0.705 FreeSolv (HFE, kcal/mol) GCN (MoleculeNet, random) 642 2.87 MPNN (MoleculeNet, random) 642 2.19 Uni-Mol (random) 642 0.60 MultiPUFFIN (all) 642 1.271 MultiPUFFIN (test-only, scaffold) 50 1.689 Lipophilicity (logâĄP P) GCN (MoleculeNet, random) 4200 0.85 MPNN (MoleculeNet, random) 4200 0.72 Uni-Mol (random) 4200 0.60 MultiPUFFIN (all) 4200 0.811 MultiPUFFIN (test-only, scaffold) 87 1.121 On the ESOL test-only scaffold subset MultiPUFFIN achieves RMSE = 0.7050.705 on 177 molecules, beating both the Delaney original model (0.750.75) and the MoleculeNet GCN baseline (0.970.97) under tighter splitting; on all 1117 ESOL molecules the RMSE drops to 0.6490.649 (R2=0.904R^2=0.904), within 0.070.07 of the random-split single-task MPNN. On the 50-molecule FreeSolv scaffold subset the RMSE is 1.6891.689 kcal/mol ; noisier than single-task Uni-Mol (0.600.60) but on a sample size where single-seed variance is large; on all 642 FreeSolv molecules the RMSE is 1.2711.271 kcal/mol (R2=0.891R^2=0.891), and the equation-level ablation of Section 3.5.3 shows that swapping in the Born head closes most of the residual gap. On the 87-molecule Lipophilicity scaffold subset the RMSE is 1.1211.121, which is the noisiest pocket of any test in this paper (n=87n=87 on a 4,200-compound benchmark). For PUFFIN-family heritage comparisons, the original single-task PUFFIN [50] reported VP RMSE =0.47=0.47 in log10âĄ(Pa) _10(Pa) on ⌠6,000 molecules; MultiPUFFINâs scaffold-split VP RMSE is 1.7891.789 on 16511651 molecules, but it retains PUFFINâs monotonicity guarantee while simultaneously delivering eight other properties. The single-task ExPUFFIN [47] reported a 37% RMSE reduction using the Andrade head; MultiPUFFIN reproduces that head behaviour at viscosity RMSE =0.326=0.326 in log10âĄ(mPaâ s) _10(mPa·s) across the broader four-source training corpus. Three caveats apply uniformly to these benchmark comparisons: (i) most cited literature uses random splitting, against MultiPUFFINâs scaffold split; (i) single-task models dedicate their full parameter budget to one property, against the 35-M-parameter backbone shared across nine; (i) the per-benchmark scaffold subsets are small (50â177 molecules), so single-seed variance is non-negligible. Despite these disadvantages MultiPUFFIN matches or beats the single-task baselines on solubility and provides a qualitative capability the comparison literature lacks: nine simultaneous predictions with intra-property thermodynamic consistency from a single forward pass. 3.9 Direct comparison with a molecular foundation model baseline To provide a direct empirical comparison against an existing molecular foundation model, we fine tuned ChemBERTa-2 [2] (a SMILES-based Transformer pretrained via masked language modelling and multi-task regression on 77 million molecules) on the identical training, validation, and test splits used for MultiPUFFIN. For each of the nine properties, a separate ChemBERTa-2 model was trained with a per-property regression head (three-layer MLP with GELU activations and dropout), z-score target normalization, Huber loss, differential learning rates (2Ă10â52Ă 10^-5 for the backbone, 2Ă10â42Ă 10^-4 for the head), and early stopping with patience of 5 epochs. Each per-property model has approximately 3.8 million trainable parameters. Table 13 reports RMSE, MAE, and coefficient of determination (R2R^2) on the shared held-out test set of 88778877 unique molecules. For vapor pressure and viscosity, whose raw values span several orders of magnitude, both models are evaluated in log10 _10 space to enable a fair comparison. A note on comparison fairness is warranted. MultiPUFFIN receives temperature as an explicit input variable through the experimental encoder (Section 2.8), whereas ChemBERTa-2 receives only a SMILES string and cannot represent the temperature at which a measurement was taken. For the six temperature-independent properties (solubility at 298.15 K298.15\,K, logâĄP P, hydration free energy, boiling point, melting point, and flash point), the input information available to both models is equivalent at the molecular level, and the comparison is a direct assessment of representational quality and domain-informed inductive biases. For the three temperature-dependent properties (vapor pressure, viscosity, and heat capacity), the comparison is inherently asymmetric: MultiPUFFIN can represent temperature-conditioned measurements while ChemBERTa-2 predicts a single value per molecule regardless of temperature. This asymmetry is not an artifact of the experimental design but rather a direct consequence of the domain-informed architectureâs capabilities: providing a SMILES-only baseline with temperature conditioning would require a fundamental architectural modification (adding an experimental encoder), at which point the modified baseline would no longer be ChemBERTa-2. The asymmetric comparison is therefore informative precisely because it quantifies what a domain-informed architecture with temperature conditioning achieves relative to a purely data-driven SMILES-based model without it. Results for temperature-independent and temperature-dependent properties are reported separately in Table 13 to allow readers to assess each regime independently. Table 13: Direct comparison of MultiPUFFIN and fine tuned ChemBERTa-2 on nine physicochemical properties, evaluated on the same scaffold-based test split. Each ChemBERTa-2 baseline is a separate per-property model fine tuned on identical train/val/test partitions; MultiPUFFIN is a single multi-task model that produces all nine predictions from a shared backbone. Bold indicates the winner per property and metric. Properties marked with â are temperature-dependent; ChemBERTa-2 receives only a SMILES string as input and therefore cannot distinguish measurements of the same molecule at different temperatures. Vapor pressure and viscosity metrics are computed in log10 _10 space. Property ntestn_test MultiPUFFIN (ours) ChemBERTa-2 RMSE MAE R2R^2 RMSE MAE R2R^2 Solubility (log mol/L) 27002700 1.182 0.843 0.537 1.234 0.868 0.496 logâĄP P 15701570 0.940 0.565 0.821 1.245 0.779 0.687 Hydration free energy (kcal/mol) 269 0.935 0.432 0.913 1.463 0.924 0.788 Boiling point (K) 21062106 39.610 21.520 0.748 47.287 27.849 0.640 Melting point (K) 62066206 36.660 26.060 0.780 47.946 35.472 0.624 Flash point (K) 28402840 25.510 16.080 0.816 31.638 20.736 0.717 Vapor pressureâ (log10 _10 Pa) 16511651 1.766 1.196 0.598 6.926 6.267 â5.431-5.431 Viscosityâ (log10 _10 mPaâ ·s) 331 0.273 0.164 0.779 6.030 5.238 â106.961-106.961 Heat capacityâ (J/mol/K) 297 18.390 7.930 0.940 23.453 14.037 0.902 Median R2R^2 (all 9 properties) 0.780 0.640 Mean R2R^2 (6 temp.-independent) 0.769 0.658 Wins (per-property R2R^2) 9 of 9 0 of 9 Three findings emerge. First, MultiPUFFIN achieves higher R2R^2 than ChemBERTa-2 on all nine properties, a result that was not previously reported for a single multi-task model at this scale. Among the six temperature-independent properties the per-property margin ÎâR2 R^2 ranges from +0.041+0.041 (solubility, where both models are noise-limited) to +0.156+0.156 (melting point), with hydration free energy at +0.125+0.125, boiling point at +0.108+0.108, logâĄP P at +0.134+0.134 and flash point at +0.099+0.099; the in-scope mean R2R^2 over those six is 0.7690.769 against 0.6580.658 for ChemBERTa-2, on ⌠20002000Ă fewer labelled molecules and roughly one-tenth of the combined parameter count of the nine per-property ChemBERTa-2 baselines. Second, the most dramatic gap is on the two genuinely temperature-dependent properties (vapor pressure and viscosity), where ChemBERTa-2âs R2R^2 collapses to â5.4-5.4 and â107-107 in log10 _10 space because a SMILES-only input cannot distinguish measurements of the same molecule at different temperatures; MultiPUFFINâs Antoine and Andrade heads encode this temperature dependence explicitly. Heat capacity is a partial exception: ChemBERTa-2 still reaches R2=0.902R^2=0.902 (against MultiPUFFINâs 0.9210.921) because CpC_p varies smoothly and over a narrow T range. Third, the comparison contrasts two paradigms cleanly: data-driven SMILES scaling on 77 M unlabelled molecules with unconstrained MLP heads, versus domain-informed inductive bias neurons over a multimodal structural backbone trained on 37 96837\,968 labelled molecules. The empirical result is that domain knowledge wins on all nine properties at orders-of-magnitude lower data cost, confirming that incorporating established physics into the architecture is a more data-efficient path than scaling unlabelled pretraining alone. Figure 15 visualises the per-property bar chart and the head-to-head ÎâR2 R^2 margin: VP and viscosity values below â0.2-0.2 are clipped with the actual value annotated, and panel (b) restricts the ÎâR2 R^2 view to the seven properties where both models produce meaningful predictions. Figure 15: Direct comparison of MultiPUFFIN (single multi-task model) against fine tuned ChemBERTa-2 (nine separately trained per-property models) on the identical scaffold-based test split. (a) Per-property test R2R^2 (vapor pressure and viscosity in log10 _10 scale for fair comparison). ChemBERTa-2 values below â0.2-0.2 are clipped; the annotated number reports the actual value. (b) Head-to-head R2R^2 margin (MPâCBMP-CB-2) on the seven properties where both models are competitive; positive bars (blue) indicate a MultiPUFFIN advantage. Median R2R^2 across all nine full-test properties: MultiPUFFIN =0.780=0.780, ChemBERTa-2 =0.640=0.640. MultiPUFFIN wins all nine properties, with the dramatic advantages on vapor pressure and viscosity reflecting the temperature-conditioning capability that SMILES-only models fundamentally lack. 3.10 Limitations Six caveats are worth surfacing before the conclusions. Statistical rigour. All results come from a single random seed and a single scaffold split; the test-set-size bootstrap R2R^2 uncertainty for the smallest test sets (HFE 269, viscosity 331, CpC_p 297) is ±0.03± 0.03â0.050.05 at the 95% level, so the smallest equation-level ablation deltas (e.g. Wagner vs. Antoine, ÎâR2â+0.004 R^2â+0.004) fall inside the seed-to-seed band. The ChemBERTa-2 comparisons over the six temperature-independent properties (test sets of 269â6,206) are larger than this band, and the largest reported margin (ÎâR2=+0.125 R^2=+0.125 on HFE) is robust. Multi-seed reporting is the next obvious step. Architectural cost-benefit. The 35-M-parameter trimodal backbone is heavier than a ⌠1-M-parameter GCN that would suffice for logâĄP P alone; the multimodal pipeline is justified by the joint accuracy on the temperature-dependent and 3D-sensitive properties, not by a uniform per-property advantage. Comparison gap. A direct head-to-head against the chemical-engineering deployment tools (Joback, Lydersen, Nannoolal, DIPPR, UNIFAC, COSMO-RS) is not in scope here and is the natural next benchmark for the process-engineering audience. Process-simulator readiness. The intra-property monotonicity guarantees are not enough for direct Aspen/gPROMS integration: analytical âP/âTâ P/â T, mixture extension and numerical stability across the full T,PT,P range remain open. Unweighted mean R2R^2. The R2=0.784R^2=0.784 headline is an unweighted average across very different test-set sizes (269â6,206); for low-variance test subsets such as the acyclic-solvent-ether viscosity sub-group, RMSE/MAE are the honest summary and R2R^2 can be strongly negative on a moderate RMSE. Missing condition inputs. The condition-dependency ledger of the conditionâdata analysis below identifies pH (ionisable solutes), ambient pressure (BP), polymorph (TmT_m) and test-cup protocol (FP) as drivers of the measured value that are not currently populated with enough variance for the corresponding modules to learn from. The bottleneck deployment sub-groups (aliphatic primary amines, aromatic-carboxylic-acid solubility, dicarboxylic acids) coincide with the cells flagged as âhas, sparseâ or âmodule ready, no dataâ in that ledger. Closing this gap requires per-condition data acquisition rather than architectural change. 4 Conclusions MultiPUFFIN generalises the domain-informed inductive bias paradigm of PUFFIN/ExPUFFIN from single-property, single-modality regression to a multi-task foundation model predicting nine thermophysical properties from a single forward pass. The 512-dimensional unified embedding produced by the trimodal backbone (GCN, Transformer, SchNet) plus the five-module identity-initialised condition stack feeds property-specific thermophysical equation heads (Antoine, Andrade, van ât Hoff, Born, Shomate) that enforce intra-property thermodynamic consistency by construction. The four-stage training protocol (SSL pretraining on 500 000500\,000 unlabelled PubChem molecules, joint supervised multi-task training, backbone-unfrozen targeted fine tune on the augmented dataset, and per-property G+ applicability-domain evaluation) yields a deployable artifact with an in-scope mean test R2=0.784R^2=0.784 across the nine properties. The headline scientific finding is that domain knowledge is a more data-efficient route to accurate molecular property prediction than brute-force scaling of pretraining data. A single MultiPUFFIN model trained on 37 96837\,968 labelled molecules beats nine separately fine tuned ChemBERTa-2 baselines (pretrained on 77 M molecules) on all nine targets, with dramatic advantages on the temperature-dependent properties where a SMILES-only model cannot distinguish measurements of the same molecule at different temperatures. The advantage decomposes into two complementary mechanisms: the domain-informed heads compress the structureâproperty mapping into known functional forms, and the multimodal backbone supplies the topology, geometry and conformational context those forms require. The primary residual limitation is multi-task capacity dilution: the shared 35-M-parameter backbone distributes representational capacity across nine simultaneous objectives. The clearest improvement paths are targeted per-condition data acquisition (pH, polymorph, cup protocol), expansion of the labelled corpus toward 40,000â60,000 unique molecules through public-database integration, and cross-property thermodynamic coupling that would lift the current intra-property guarantees toward rigorous inter-property thermodynamic relations. MultiPUFFIN establishes domain-informed multimodal foundation models as a viable, data-efficient and computationally efficient paradigm for thermophysical property prediction in chemical-engineering deployments. Code and data availability The source code, trained model weights, and the curated multi-property dataset (SMILES strings with property labels and train/validation/test split assignments) will be made publicly available upon publication at https://github.com/idelfonsog/MultiPUFFIN. References [1] ABB-ADD consortium (2023) ABB-ADD liquid heat capacity dataset. Note: Curated liquid heat capacity data for organic compounds Cited by: §2.1. [2] W. Ahmad, E. Simon, S. Chithrananda, G. Grand, and B. Ramsundar (2022) ChemBERTa-2: towards chemical foundation models. arXiv preprint arXiv:2209.01712. Cited by: Table 1, §1, §3.9. [3] E. N. d. C. Andrade (1930) The viscosity of liquids. Nature 125, p. 309â310. Cited by: §1, Table 4. [4] T. BaltruĆĄaitis, C. Ahuja, and L. Morency (2019) Multimodal machine learning: a survey and taxonomy. IEEE Transactions on Pattern Analysis and Machine Intelligence 41 (2), p. 423â443. Cited by: §1. [5] I. Batatia et al. (2025) Foundation models for atomistic simulation of chemistry and materials. Nature Reviews Chemistry. External Links: Document Cited by: §1. [6] J. Bradley, A. Lang, and A. Williams (2014) Open melting point data. Note: Figshare External Links: Document Cited by: 3rd item, §2.1. [7] T. Cai, S. Luo, K. Xu, D. He, T. Liu, and L. Wang (2021) GraphNorm: a principled approach to accelerating graph neural network training. In International Conference on Machine Learning, p. 1204â1215. Cited by: §2.5. [8] R. Chen et al. (2025) Pretraining graph transformer for molecular representation with fusion of multimodal information. Information Fusion 115, p. 102784. Cited by: Table 1, §1. [9] A. K. Chew et al. (2024) Advancing material property prediction: using physics-informed machine learning models for viscosity. Journal of Cheminformatics 16, p. 31. Cited by: §2.1. [10] J. Choi, Y. Nam, Y. Choi, and Y. Jung (2025) A perspective on foundation models in chemistry. JACS Au 5 (4). External Links: Document Cited by: §1. [11] M. Crawshaw (2020) Multi-task learning with deep neural networks: a survey. arXiv preprint arXiv:2009.09796. Cited by: §3.4. [12] M. Das, A. Ghosh, and R. B. Sunoj (2024) Advances in machine learning with chemical language models in molecular property and reaction outcome predictions. Journal of Computational Chemistry 45 (14), p. 1160â1176. Cited by: §1. [13] J. S. Delaney (2004) ESOL: estimating aqueous solubility directly from molecular structure. Journal of Chemical Information and Computer Sciences 44 (3), p. 1000â1005. Cited by: §3.8. [14] J. Devlin, M. Chang, K. Lee, and K. Toutanova (2019) BERT: pre-training of deep bidirectional transformers for language understanding. North American Chapter of the Association for Computational Linguistics, p. 4171â4186. Cited by: §2.2. [15] European Chemicals Agency (2024) ECHA REACH registered substances database. Note: https://echa.europa.eu/information-on-chemicals/registered-substancesAccessed: 2024 Cited by: §2.1. [16] M. Fey and J. E. Lenssen (2019) Fast graph representation learning with PyTorch Geometric. ICLR Workshop on Representation Learning on Graphs and Manifolds. Cited by: §2.13. [17] M. Frenkel, R. D. Chirico, V. Diky, X. Yan, Q. Dong, and C. Muzny (2005) ThermoML: an XML-based approach for storage and exchange of experimental and critically evaluated thermophysical and thermochemical property data. 1. Experimental data. Journal of Chemical & Engineering Data 50 (1), p. 42â50. Cited by: §2.1. [18] J. Gasteiger, J. GroĂ, and S. GĂŒnnemann (2020) Directional message passing for molecular graphs. International Conference on Learning Representations. Cited by: §1. [19] A. Gaulton, A. Hersey, M. Nowotka, A. P. Bento, J. Chambers, D. Mendez, P. Mutowo, F. Atkinson, L. J. Bellis, E. CibriĂĄn-Uhalte, et al. (2017) The ChEMBL database in 2017. Nucleic Acids Research 45 (D1), p. D945âD954. Cited by: §2.1. [20] T. A. Halgren (1996) Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. Journal of Computational Chemistry 17 (5-6), p. 490â519. Cited by: §2.1. [21] K. He, X. Zhang, S. Ren, and J. Sun (2016) Deep residual learning for image recognition. In IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 770â778. Cited by: §2.5. [22] D. Hendrycks and K. Gimpel (2016) Gaussian error linear units (GELUs). arXiv preprint arXiv:1606.08415. Cited by: §2.6. [23] S. Ioffe and C. Szegedy (2015) Batch normalization: accelerating deep network training by reducing internal covariate shift. International Conference on Machine Learning, p. 448â456. Cited by: §2.5. [24] M. Karthikeyan, R. C. Glen, and A. Bender (2005) General melting point prediction based on a diverse compound data set and artificial neural networks. Journal of Chemical Information and Modeling 45 (3), p. 581â590. Cited by: §3.2. [25] A. Kazakov, J. W. Magee, R. D. Chirico, E. Paulechka, V. Diky, C. D. Muzny, K. Kroenlein, and M. Frenkel (2014) NIST/TRC web thermo tables: an online resource for critically evaluated thermodynamic property data. Journal of Research of the National Institute of Standards and Technology 119, p. 136â156. Cited by: 3rd item. [26] A. Kendall, Y. Gal, and R. Cipolla (2018) Multi-task learning using uncertainty to weigh losses for scene geometry and semantics. IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 7482â7491. Cited by: §2.13. [27] S. Kim, J. Chen, T. Cheng, A. Gindulyte, J. He, S. He, Q. Li, B. A. Shoemaker, P. A. Thiessen, B. Yu, et al. (2023) PubChem 2023 update. Nucleic Acids Research 51 (D1), p. D1373âD1380. Cited by: §2.1, §2.2. [28] T. N. Kipf and M. Welling (2017) Semi-supervised classification with graph convolutional networks. International Conference on Learning Representations. Cited by: §2.5. [29] H. Li, D. Zhao, and J. Zeng (2022) KPGT: knowledge-guided pre-training of graph transformer for molecular property prediction. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, p. 857â867. Cited by: Table 1, §1. [30] J. Li et al. (2025) MoleculeFormer is a GCN-transformer architecture for molecular property prediction. Communications Biology 8, p. 1668. Cited by: Table 1. [31] Q. Li, Z. Han, and X. Wu (2018) Deeper insights into graph convolutional networks for semi-supervised learning. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 32, p. 3538â3545. Cited by: §3.6. [32] Y. Li et al. (2025) MolPrompt: improving multi-modal molecular pre-training with knowledge prompts. Bioinformatics 41 (9). Cited by: §1. [33] Y. Li et al. (2024) Advancing molecular graph-text pre-training via fine-grained alignment. arXiv preprint arXiv:2409.14106. Cited by: §1. [34] P. Liu, Y. Ren, J. Tao, and Z. Ren (2024) GIT-Mol: a multi-modal large language model for molecular science with graph, image, and text. Computers in Biology and Medicine 171, p. 108073. Cited by: Table 1, §1. [35] S. Liu, W. Nie, C. Wang, J. Lu, Z. Qiao, L. Liu, J. Tang, C. Xiao, and A. Anandkumar (2023) Multi-modal molecule structureâtext model for text-based retrieval and editing. Nature Machine Intelligence 5, p. 1447â1457. Cited by: Table 1, §1, §2.4. [36] I. Loshchilov and F. Hutter (2017) SGDR: stochastic gradient descent with warm restarts. International Conference on Learning Representations. Cited by: §2.13. [37] I. Loshchilov and F. Hutter (2019) Decoupled weight decay regularization. International Conference on Learning Representations. Cited by: 2nd item. [38] S. Lu et al. (2024) Uni-Mol2: exploring molecular pretraining model at scale. Advances in Neural Information Processing Systems 37. Cited by: Table 1, §1. [39] K. Mansouri, C. M. Grulke, R. S. Judson, and A. J. Williams (2018) OPERA models for predicting physicochemical properties and environmental fate endpoints. Journal of Cheminformatics 10, p. 10. Cited by: §2.1. [40] L. McInnes, J. Healy, and J. Melville (2018) UMAP: uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426. Cited by: §3.7. [41] O. MĂ©ndez-Lucio et al. (2024) MolE: a foundation model for molecular graphs using disentangled attention. Nature Communications 15, p. 9431. Cited by: Table 1, §1. [42] D. L. Mobley and J. P. Guthrie (2014) FreeSolv: a database of experimental and calculated hydration free energies, with input files. Journal of Computer-Aided Molecular Design 28 (7), p. 711â720. Cited by: §2.1, §3.8. [43] Y. Nannoolal, J. Rarey, and D. Ramjugernath (2008) Estimation of pure component properties: part 3. estimation of the vapor pressure of non-electrolyte organic compounds via group contributions and group interactions. Fluid Phase Equilibria 269 (1-2), p. 117â133. Cited by: §1, Table 4, 3rd item. [44] C. Pang, H. H. Tong, and L. Wei (2023) Advanced deep learning methods for molecular property prediction. Quantitative Biology 11 (4), p. 395â404. Cited by: §1. [45] A. Paszke et al. (2019) PyTorch: an imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems 32. Cited by: §2.13. [46] RDKit community (2024) RDKit: open-source cheminformatics. Note: https://w.rdkit.org Cited by: §2.1, §2.8. [47] C. M. Rebello, U. Di Caprio, J. Steen-Hansen, B. Rodrigues, E. A. Costa, A. R. dos Santos, F. Esposito, M. E. Leblebici, and I. B. R. Nogueira (2025) ExPUFFIN: thermodynamic consistent viscosity prediction in an extended path-unifying feed-forward interfaced network. arXiv preprint arXiv:2512.06927. Cited by: Table 1, §1, §1, §2.11, §3.8. [48] S. Riniker and G. A. Landrum (2015) Better informed distance geometry: using what we know to improve conformation generation. Journal of Chemical Information and Modeling 55 (12), p. 2562â2574. Cited by: §2.1. [49] J. G. Rittig, K. C. Felton, A. A. Lapkin, and A. Mitsos (2023) GibbsâDuhem-informed neural networks for binary activity coefficient prediction. Digital Discovery 2, p. 1752â1767. Cited by: §1. [50] V. V. Santana, C. M. Rebello, L. P. Queiroz, A. M. Ribeiro, N. Shardt, and I. B. R. Nogueira (2024) PUFFIN: a path-unifying feed-forward interfaced network for vapor pressure prediction. Chemical Engineering Science 286, p. 119623. Cited by: Table 1, §1, §1, 6th item, §2.10, §2.11, §3.8. [51] K. T. SchĂŒt, H. E. Sauceda, P. Kindermans, A. Tkatchenko, and K. MĂŒller (2018) SchNetâa deep learning architecture for molecules and materials. The Journal of Chemical Physics 148 (24). Cited by: §1, §2.7. [52] A. M. Schweidtmann, J. G. Rittig, A. König, M. Grohe, A. Mitsos, and M. Dahmen (2020) Graph neural networks for prediction of fuel ignition quality. Energy & Fuels 34 (9), p. 11395â11407. Cited by: §1. [53] M. C. Sorkun, A. Khetan, and S. Er (2019) AqSolDB, a curated reference set of aqueous solubility and 2D descriptors for a diverse set of compounds. Scientific Data 6, p. 143. Cited by: §2.1. [54] B. Su, D. Du, Z. Yang, Y. Zhou, J. Li, A. Rao, H. Sun, Z. Lu, and J. Wen (2022) A molecular multimodal foundation model associating molecule graphs with natural language. arXiv preprint arXiv:2209.05481. Cited by: Table 1, §1. [55] Y. Sun, Z. Yin, and L. Yin (2019) Prediction of flash points of organic compounds using quantitative structure-property relationship (QSPR) methods. Journal of Hazardous Materials 371, p. 518â524. Cited by: §2.1. [56] G. W. Thomson (1946) The antoine equation for vapor-pressure data. Chemical Reviews 38 (1), p. 1â39. Cited by: §1, Table 4. [57] Uni-Mol Team (2025) Uni-Mol3: a multi-molecular foundation model for advancing organic reaction modeling. arXiv preprint arXiv:2508.00920. Cited by: §1. [58] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ć. Kaiser, and I. Polosukhin (2017) Attention is all you need. Advances in Neural Information Processing Systems 30. Cited by: §2.6, §2.6. [59] F. H. Vermeire and W. H. Green (2021) Transfer learning for solvation free energies: from quantum chemistry to experiments. Chemical Engineering Journal 418, p. 129307. Cited by: §1. [60] S. Wang, Y. Guo, Y. Wang, H. Sun, and J. Huang (2019) SMILES-BERT: large scale unsupervised pre-training for molecular property prediction. ACM International Conference on Bioinformatics, Computational Biology and Health Informatics, p. 429â436. Cited by: Table 1, §1. [61] Y. Wang, J. Wang, Z. Cao, and A. B. Farimani (2022) Molecular contrastive learning of representations via graph neural networks. Nature Machine Intelligence 4, p. 279â287. Cited by: §1. [62] D. Weininger (1988) SMILES, a chemical language and information system. 1. introduction to methodology and encoding rules. Journal of Chemical Information and Computer Sciences 28 (1), p. 31â36. Cited by: §1. [63] Z. Wu, B. Ramsundar, E. N. Feinberg, J. Gomes, C. Geniesse, A. S. Pappu, K. Leswing, and V. Pande (2018) MoleculeNet: a benchmark for molecular machine learning. Chemical Science 9 (2), p. 513â530. Cited by: §2.1, §2.1, §2.1, §3.8. [64] R. Xiong, Y. Yang, D. He, K. Zheng, S. Zheng, C. Xing, H. Zhang, Y. Lan, L. Wang, and T. Liu (2020) On layer normalization in the Transformer architecture. International Conference on Machine Learning, p. 10524â10533. Cited by: §2.6. [65] G. Zhou, Z. Gao, Q. Ding, H. Zheng, H. Xu, Z. Wei, L. Zhang, and G. Ke (2023) Uni-Mol: a universal 3D molecular representation learning framework. International Conference on Learning Representations. Cited by: Table 1, §1, §2.1, §2.4.