Predicting solvation free energies with an implicit solvent machine learning potential

Solvation free energy, and notably hydration free energy, is generally recognized as a fundamental thermodynamic quantity of interest in computational chemistry. Defined as the work done when transferring a molecule from the gas phase to the solution (water in the case of hydration free energy), it enables the computation of several key physicochemical properties of molecules, such as solubility, partition coefficients, activity coefficients, and binding free energies in solutions.^1,2 These properties are of great importance to the pharmaceutical, environmental, and materials sciences,^3–9 prompting the organization of blind prediction SAMPL challenges^10–12 with hydration free energy as one of the main targets. In addition, Mobley and Guthrie compiled and curated a FreeSolv database of experimentally measured hydration free energies for small neutral molecules in water.^13,14

A wide spectrum of methods is available to calculate the solvation free energy, ranging from traditional approaches such as continuum solvation models^15,16 to recent machine learning (ML) algorithms^17–26 and their combinations.^27–29 The alchemical methods with Molecular Dynamics (MD) simulations^14,30,31 are typically assumed to be highly accurate but computationally expensive.^32,33 However, both the fidelity and the efficiency highly depend on the explicitly treated degrees of freedom and the employed potential energy model.

In implicit solvent models, the solvent molecules (e.g., water molecules) are not explicitly present in the system (as in explicit solvent models); instead, the interactions are modified to account for the solvent effects.³⁴ The number of degrees of freedom is thereby greatly reduced, resulting in large computational gains. Classical implicit models treat the solvent as a continuum medium with specific dielectric and interfacial properties. The typical approach decouples the electrostatic (polar) and nonpolar interactions. The former can be approximated by solving the Poisson–Boltzmann equation or further simplified with the popular generalized Born model, while the latter is most often estimated via the solvent-accessible surface area.^34,35 While recent advances, such as using ML to predict the generalized Born radii,³⁶ have increased the accuracy of these models, they are still in considerable disagreement with experimental data and explicit solvent models.^26,37,38 In particular, the solvation free energy root mean square error (RMSE) is ∼3.6 kcal/mol for both the Poisson–Boltzmann and generalized Born models. By specifically optimizing the nonpolar interactions using hydration free energy data, the test set RMSE can be decreased to 1.68 kcal/mol.³⁹ 

Higher accuracy can be achieved through ab initio methods, based on quantum chemical calculations, albeit at the cost of computational demand. Commonly utilized is the COSMO family,^40,41 with COSMO-RS (Conductor-like Screening Model for Real Solvents) being particularly notable in property predictions. In particular, extensions thereof have achieved high fidelity predictions of solvation free energies (MAE = 0.52 kcal/mol).⁴² 

An alternative to the computational hurdle imposed by ab initio methods, without significantly compromising predictive accuracy, is offered by the many-body, flexible potential energy surfaces that characterize ML potentials. During the past decade, ML potentials were deployed to derive atomistic⁴³ and coarse-grained⁴⁴ models based on training data provided by high-fidelity simulations,⁴⁵ experiments,⁴⁶ or both.⁴⁷ The derivation of empirical force fields from the experimental data has been extensively studied (see, for example, Refs. 48–51), and this approach naturally extends to ML potentials. In the context of implicit water models, the capacity of deep ML potentials is large enough to compensate for the removal of solvent degrees of freedom. The predicted solute properties, such as the conformational landscape of proteins, match those obtained from the reference explicit water model.^26,52–58 However, these studies used classical atomistic force fields as the data-generating model, inherently limiting the attainable accuracy of the resulting ML potential. As shown previously, and also in this work, classical atomistic models such as General Amber Force Field (GAFF) and CHARMM General Force Field (CGenFF) systematically overestimate hydration free energies.³⁰ Modifying the Lennard-Jones parameters yields an improved agreement with experiments but can, at the same time, negatively affect other properties.³⁰ 

Meanwhile, more accurate ab initio training data are prohibitively expensive for solutes in explicit bulk water. Consequently, frequently used Density Functional Theory (DFT) training databases such as QM9,⁵⁹ ANI,⁶⁰ and QM7-X⁶¹ contain small organic molecules in vacuum or at best contain samples with a few randomly placed water molecules around solutes.^62,63 Moreover, implicit solvent models are coarse-grained models. When trained with the common force-matching approach,⁶⁴ coarse-grained ML models require much more data than their atomistic counterparts.⁶⁵ The underlying reason is the surjective atomistic-to-coarse-grained mapping, resulting in noisy force labels. As an illustration, the parameterization of coarse-grained models for alanine dipeptide in implicit water required a dataset with 10⁶ configurations,^45,66 a number easily obtainable with classical force fields but vastly out of reach for ab initio calculations.

In this work, we present a Solvation Free Energy Path Reweighting (ReSolv) approach to parameterize an ML potential for small organic molecules in an implicit aqueous solvent. The above-mentioned difficulties are circumvented with a two-stage training procedure, utilizing first the DFT database of molecules in vacuum and then the experimental hydration free energy database. The second stage entails a non-trivial training since the hydration free energy is not a direct output of an ML model but instead involves molecular simulations. Constructing the free energy integration path along the ML model training process and utilizing the Zwanzig reweighting scheme enables us to perform efficient training that avoids differentiating through the molecular simulations. The ReSolv model predicts hydration free energies more accurately than the classical explicit solvent models despite being an implicit solvent (i.e., coarse-grained) model. In addition, ReSolv’s predictions are not systematically biased and are more robust for molecules with large negative hydration free energies. We also investigate error correlation between different modeling approaches and point out several potentially erroneous data points that should be reconsidered in future database curation.

The training methodology of ReSolv consists of two consecutive stages as shown in Fig. 1. In the first stage, we parameterize the ML potential for molecules in vacuum. The model takes as an input a configurational state of the molecule S and predicts the potential energy, i.e., U_vac = U(S; θ_vac). The forces on the atoms are computed as the negative derivative of the potential with respect to the atoms’ position vectors. In training, we aim to adjust the parameters of the model θ_vac such that the predictions match the corresponding energies U_DFT and forces F_DFT in the ab initio database, i.e., using a training approach via Eq. (5).

For the present task of learning the solvation free energy, DiffTRe could be directly employed as presented above. Nevertheless, the free energy difference between the potentials U_vac and U_sol is typically estimated by constructing a free energy integration path, e.g., with a new potential energy function U(λ) = (λ − 1)U_vac + λU_sol, and estimating the free energy differences between discrete steps of λ. This approach would require as many simulations as there are λ steps. In addition, intermediate steps would involve an MD simulation using a linear combination of two ML potentials, rendering the update step computationally expensive.

Here, we instead use the DiffTRe learning process itself as a free energy integration path. We initialize the parameters with the pre-trained θ_vac parameters and iteratively perturb them during training toward the desired θ_sol parameters (see Fig. 1). At each update step, we compute the free energy difference between the new and previous potentials using a hybrid of the free energy perturbation⁶⁷ and the Bennett Acceptance Ratio (BAR) methods.⁶⁹ The free energy differences are accumulated along the training, which, at the end of the training, yield the total free energy difference between U_vac and U_sol or ΔA.

Both vacuum and solvent ML potentials are based on the GNN architecture NequIP⁷⁰ as implemented in JAX MD.⁷¹ Computations were performed with double precision on NVIDIA A100 80GB GPUs. The architectural hyperparameters of NequIP are listed in Table S2 of the supplementary material.

In the second stage, we employ the RDKit⁷² to generate the 3D structure of a molecule based on the SMILES string provided by the FreeSolv database. Next, we perform the energy minimization of structures with the MMFF94 force field.^73–77 The obtained configurations are used to run a 300 ps initial equilibration simulation followed by a 200 ps production run for each molecule. All simulations are performed in the NVT ensemble and numerically integrated with the velocity Verlet scheme using a 1 fs time step. To match the conditions of the experimental data, the temperature is kept at 298.15 K with the Langevin thermostat with the damping factor set to 1/ps. We check the stability of our simulations by evaluating that the last configuration of the molecules can be represented by a single graph and that no atoms diverged from the molecule. In addition, to see whether our models yield physically reasonable trajectories, we compare the bonds, angles, and dihedrals of the U_sol trajectories for three randomly selected molecules to the bonds, angles, and dihedrals as obtained after an energy minimization with the UFF using RDKit⁷² and openFF-2.2.1 using openff-toolkit⁷⁸ and openMM;⁷⁹ see Figs. S9 and S10 of the supplementary material, respectively. In training, new simulations are initialized with the last configuration of the previous simulation and consist of a 50 ps equilibration and a 200 ps production run. We sample every 5 ps during the production run, i.e., 40 samples per trajectory. The reweighting effective sample size threshold is fixed to $N̄eff=0.9$⁠. The numerical optimization hyperparameters are reported in Table S3 of the supplementary material.

The reported hydration free energies are computed with the BAR method using only the end states, i.e., vacuum and water states. As previously reported, intermediate states are not necessary for implicit solvent models due to the sufficient overlap between the vacuum and solvated ensembles of the solute.³³ We first run 300 ps of equilibration and sample every 5 ps in the subsequent 200 ps production run. We perform convergence tests (Fig. S8 of the supplementary material) for 12 randomly selected molecules in the test set to ensure sufficient sampling.

In this work, we employ ReSolv to learn the hydration free energy of small organic molecules, which, by construction, also yields an ML potential for molecules in implicit water. The architecture of the ML potential is based on NequIP,⁷⁰ a data efficient E(3)-equivariant graph neural network. For stage 1 of training, we use the QM7-X dataset,⁶¹ providing the target DFT energies and forces. The dataset consists of 4.2 × 10⁶ small organic molecule samples, including equilibrium states, structural isomers, structural stereoisomers, and off-equilibrium configurations. The considered molecules are composed of the heavy atom set {C, N, O, S, Cl}. They have up to seven heavy atoms or four to twenty-three atoms in total, including hydrogen. We randomly split the dataset into train (89.8%), validation (10.0%), and test (0.2%) sets. The ReSolv vacuum model (U_vac) achieves a test set mean absolute error (MAE) within literature benchmarks^80,81 (Table I). The errors are close to the QM7-X dataset error, i.e., the DFT calculations were computed at the PBE0+MBD level, with a precision of 10⁻³ eV and 10⁻⁴ eV/Å for the energies and forces, respectively.

For stage 2 of training, we utilize the FreeSolv database.¹⁴ It contains hydration free energies of 643 small molecules with neutral charges, representing compounds relevant to drug-like molecules. Keeping consistent with the heavy atom types of the QM7-X dataset, i.e., {C, N, O, S, Cl}, we extract the corresponding subset with 559 molecules. We also excluded 17 molecules for which the ReSolv vacuum model yielded unstable simulations and an additional five molecules that were unstable during training. We split this subgroup into train and test sets of roughly 70% and 30% proportions, which amounts to 375 and 162 molecules, respectively. When dividing the data, we kept the occurrences of various heavy atom combinations consistent within the two sets, i.e., the training set contains roughly 70% of the molecules with heavy atom combinations N, N–C, N–C–O, etc. This choice ensures that the various heavy atom combinations are represented during training and testing. In addition, we ensured that the test set consists exclusively of molecules with chemical functional groups that appear in the training set, i.e., any molecule that is a lone representative of a functional group formed part of the training set. The functional group classification was taken from the FreeSolv database.¹³ The polyfunctional molecules were classified according to the first listed functional group as in Ref. 31. Apart from satisfying these two conditions, the molecules are split randomly.

The ReSolv’s test set errors are reported in Table I and compared with the explicit solvent classical atomistic models Amber (GAFF) and CHARMM (CGenFF). We compute the errors for the same test set molecules to enable a direct comparison. Data for Amber are provided in the FreeSolv database,¹³ while for CHARMM, the values are taken from Ref. 31. In the CHARMM error computation, two molecules (Mobley ID 6359135 and 2146331) are excluded because their data are not provided in Ref. 31 since the force field parameters could not be generated. The ReSolv’s MAE is close to the average uncertainty of the FreeSolv database (0.6 kcal/mol). Moreover, ReSolv outperforms the explicit solvent classical atomistic models by a large margin. The Amber and CHARMM models display a similar MAE/RMSE for the molecules in the training set (Fig. S1 of the supplementary material), confirming that the test set is representative of the entire considered database. In addition, we trained the ReSolv model using a random splitting of the data and achieved a similar performance. In particular, the obtained MAE and RMSE are 0.65 and 0.91 kcal/mol, respectively. The corresponding errors for the classical force fields remain substantially higher (Fig. S2 of the supplementary material). The comparison with the Amber and CHARMM models is somewhat unfair, given that the two models were not explicitly parameterized to reproduce the hydration free energy data. Nevertheless, reparameterizations and corrections in this direction were previously attempted, and the resulting models still scored lower than ReSolv. For example, Boulanger et al. rescaled the Amber (GAFF) molecule–water van der Waals dispersion interaction to better reproduce the hydration free energy and reported an MAE of 0.79 kcal/mol for the final optimized model.³⁰ In another example, Scheen et al. trained an ML model to correct the Amber (GAFF) force field predictions.²⁸ The hybrid Amber/ML approach achieved an MAE of 0.76 kcal/mol on the SAMPL4 test data containing 47 samples of the FreeSolv database.

Next, we perform an in-depth error analysis. The parity plot [Fig. 2(a)] reveals systematic errors for Amber and CHARMM predictions, with most points lying above the diagonal line. To further demonstrate this point, we plot the error distributions relative to the experimental values [Fig. 2(b)]. The distributions are skewed to the left for the classical force fields, indicating that the majority of predictions are overestimated. The Amber model, for example, overestimates the hydration free energies on average by more than 1 kcal/mol as already reported by Boulanger et al.³⁰ Conversely, no systematic errors are found for the ReSolv model, evidenced by the symmetric distribution about zero. These findings are reflected in the percentage of predictions within the experimental uncertainty [Fig. 2(c)]. The ReSolv model scores the highest, followed by CHARMM and Amber in last place, albeit Amber exhibiting a lower MAE and RMSE than the CHARMM model.

The parity plot [Fig. 2(a)] also highlights the FreeSolv’s uneven distribution with respect to the hydration free energy value. There are only a few molecules with a large and negative hydration free energy. The largest two (Mobley ID 9534740 and D-mannitol) are particularly problematic for the Amber and CHARMM models. These molecules, discussed further in Sec. III C, are part of a general trend. The MAE increases with decreasing hydration free energy [Fig. 2(d)] for all models. However, the ReSolv model is the most robust model in this respect. See Figs. S3 and S1 of the supplementary material for point-wise results and training dataset results, respectively. For Amber and CHARMM models, a possible explanation could be an inaccurate partial charge assignment. As shown by Jämbeck et al.,⁸² the hydration free energy is very sensitive to the choice of charge computation method and can differ by several kcal/mol depending on the method used. We found that the absolute sum of the Gasteiger partial charge per atom is highly correlated with the hydration free energy as is the polar surface area of the molecules (Fig. S4 of the supplementary material). Indeed, keeping in mind the net-neutrality of the molecules, the absolute Gasteiger partial charge per atom is sensitive to local charge gradients contributed by polar fragments within the molecule, whence its correlation with the polar surface area. The computed Pearson correlation coefficient is r = −0.71, which is comparable to the correlation with the related polar surface area (r = −0.74) but significantly larger than correlation with other molecular properties such as dipole moment (r = −0.42) or volume of a solute (r = −0.18). In particular, we find that the absolute sum of the charges increases with decreasing hydration free energy. An accurate partial charge assignment is, therefore, particularly important for molecules with large negative hydration free energies, potentially explaining why the predictions for classical force fields tend to worsen with decreasing hydration free energy. With this insight, we depict in Fig. S5 of the supplementary material the difference between the U_sol and U_vac atomic contributions to the potential energy for a fixed configuration and find a correspondence with the Gasteiger partial charges from the heavy atoms. The energy differences highlight that the changes are predominant around polar atoms, agreeing with what we would expect based on the above correlations, supporting that the learning of U_sol is physically informed.

In Fig. 3, the test set MAE is broken down into various molecular properties. First, we examine the error with respect to chemical functional groups. The chemical space covered by the FreeSolv database is quite extensive given the relatively small size of the database.^13,14,31 ReSolv demonstrates uniform error across the chemical functionalities, which is, to some extent, expected given the training set construction. However, low MAE is also found for the functional groups rarely appearing in the training set, i.e., for the “other” group where we merge functional groups with sparse occurrence in the train and test sets. In contrast, the Amber and CHARMM models exhibit larger fluctuations in MAE for some functional groups (e.g., primary alcohol for Amber or primary amine for CHARMM) displaying notably higher errors. Note that the same conclusion can be drawn for the molecules in the training set (Fig. S6 of the supplementary material). Next, we consider the models’ robustness with respect to the size of the molecule. From 14 heavy atoms onward, an increased MAE can be seen for all models, most notably for the Amber model. Finally, in terms of heavy atom types, all three predictive models display larger errors for molecules with four heavy atom types. However, we attribute this result mainly to the small sample size containing problematic molecules, i.e., the MAE is computed on only four molecules with the nitralin molecule (discussed in Sec. III C) contributing most to the error. The absence of a similar trend for the molecules in the training set (Fig. S6 of the supplementary material) further supports our claim.

The results thus far demonstrate ReSolv’s ability to accurately predict the hydration free energies of molecules with chemical functional groups seen during training. Given the typically limited availability of training samples, it is especially desirable for ML potentials to extrapolate effectively into unseen chemical space. To assess ReSolv’s chemical generalizability, we employ the same training procedure and dataset prescribed in Sec. III A but with all samples of a given functional group (including multifunctional molecules) removed from the training set and compare the performance to the cases when seen.

We conducted two such trainings, excluding samples with either primary amine or primary alcohol functionalities. These functional groups were selected due to their challenging nature for CHARMM and Amber, respectively [see Fig. 3(a)]. With the test sets identical, we found that the errors associated with the unseen functional groups were highly consistent with those observed when the same groups were seen in training; see Fig. 4. This consistency demonstrates ReSolv’s robustness and capacity to predict well even for unencountered regions in chemical space.

We noticed several outliers in the parity plot [Fig. 2(a)] for which the predictions are largely off for all three investigated models. Error correlation between different modeling approaches can be used to identify possible erroneous data points. Previous outlier analysis stimulated corrections of some experimental data points that are part of the FreeSolv database.⁸³ For example, the initially provided experimental value for D-mannitol was −27.79 kcal/mol, which was later corrected to −23.62 kcal/mol. Nevertheless, doubts about the validity of the experimental value remained as the COSMO-RS prediction deviated by 6 kcal/mol.⁸³ As a result, some recent studies excluded D-mannitol from the test set.^17,18 Interestingly, ReSolv’s prediction of D-mannitol’s hydration free energy is −27.42 kcal/mol, in excellent agreement with the original experimental value.

We compute the Pearson correlation coefficient of signed and absolute errors to investigate the error correlation between ReSolv, Amber, and CHARMM predictions (Fig. S7 of the supplementary material). Overall, we find positive correlations between models, with the largest correlation of 0.6 between the Amber and CHARMM models. In Fig. 5, we show point-wise error correlations and mark the molecules with large and correlated errors (exact values are given in Table S1 of the supplementary material). The nitralin and pirimor molecules particularly stand out, with all models’ predictions deviating from experimental values by more than 2.5 kcal/mol. These two molecules also have a large experimental uncertainty (1.93 kcal/mol), confirming that error correlation analysis is a useful approach for inaccurate data recognition.

Finally, we compare ReSolv’s performance to explicit solvent ML potentials. Using the same ML potential architecture, implementation, and common simulation setup in the literature, the measured computational cost differs by four orders of magnitude (see the supplementary material). The ReSolv’s speedup is due to two factors: (i) faster execution of MD steps due to the reduced number of particles and (ii) reduced number of required MD steps due to enhanced sampling.

Surprisingly, ReSolv’s substantial computational gains do not come at the cost of accuracy. To the contrary, a comparison with a recent study⁸⁴ on hydration free energy predictions for six selected molecules in the FreeSolv database suggests that the ReSolv’s predictions are in better agreement with the experimental results (Fig. 6). More precisely, the ReSolv model achieves better accuracy for five out of six molecules, including ethane, which is part of our test dataset, while the other five molecules were included in the training dataset.

This work presents the Solvation Free Energy Path Reweighting (ReSolv) to efficiently learn an ML-based potential energy surface in an implicit solvent by utilizing the experimental solvation free energy data and reweighting in training. Since the framework was showcased with QM7-X and FreeSolv datasets, the trained ReSolv model enables an accurate hydration free energy prediction for small organic compounds, which can be directly utilized in, for example, de novo drug design. However, the spanned chemical space of molecules could be enlarged with other datasets. In addition, the same methodology can be used to obtain an ML potential for implicit solvents other than water. For instance, extensively available octanol–water partition data could be used to derive a model in an implicit octanol environment. Predicting lipophilicity, correlated with oil–water partition coefficient, is also a critical parameter in the pharmaceutical industry as drug candidates must be sufficiently lipophilic to penetrate the lipid core of membranes but not too lipophilic that they remain there.⁸⁵ Transfer learning strategy would likely accelerate the learning of other solvents. It is reasonable to assume that similar solvents modify the interactions similarly. For example, the knowledge gained by learning an implicit water model could benefit the model for other polar solvents, such as ethanol. This setting could be particularly useful for solvents for which the solvation free energy data are scarce.

For the FreeSolv database, ReSolv’s MAE is close to experimental uncertainty. Compared to the classical force field models, ReSolv demonstrates better accuracy for molecules with large negative hydration free energies. In water, the molecular polar surface area strongly correlates with negative solvation free energy values. This improvement, therefore, may indicate that ReSolv implicitly captures local charge distributions in its energy predictions more effectively than the classical force fields. However, further work is necessary to clarify this point. The obtained accuracy is also similar to the previously reported structure–property models, employing ML algorithms to directly predict the hydration free energy from the physics-inspired fingerprints^8,17,18,29 or molecular structure.^{9,19,20,22–26} The latter employs graph-based neural network architectures and typically perform better than the former with the reported MAE in the range of 0.58–0.76 kcal/mol and RMSE in the range of 0.82–1.23 kcal/mol. These results indicate that further improvements will likely require an enlarged and improved experimental database rather than a better modeling approach or neural network architecture. In line with this conclusion is a recent study²¹ achieving a notably lower MAE of 0.42 kcal/mol where a graph neural network was first pre-trained on a large dataset and later fine-tuned on the FreeSolv database, i.e., by exploiting the transfer learning approach. ReSolv showcases good generalization across the functional space but a decreased generalization to larger molecules, which should be considered when planning future experiments. In addition, we found error correlations between our models and classical atomistic models that were not trained on the FreeSolv database. Molecules with large and correlated errors, such as nitralin and pirimor, call for a reevaluation of the experimental data. These points should be excluded when training future models to avoid detrimental effects.

In a broader context, ReSolv could be a first step toward a general implicit solvent ML potential with better accuracy and efficiency than classical atomistic models. Solvation free energy encodes solute–solvent interactions, which in turn govern many biomolecular processes, including folding, aggregation, and ligand binding. In vivo, these processes occur in an aqueous solution rendering the hydration free energy a general benchmark for classical force field validation,^1,10,86 fine-tuning,^2,87 comparison,^88,89 and calibration.^90,91 Nevertheless, before using the ReSolv model in general simulations, the model needs to be tested and potentially additionally trained on other available experimental data, which we leave for future studies. Similar as for other implicit solvent ML potentials, using a prior potential (i.e., a fixed simple potential) is likely also necessary to achieve stable simulations of large macromolecules.⁶⁵ Concerning computational efficiency, the ReSolv and ML-based implicit solvent models, in general, provide a speedup compared to the classical explicit solvent force fields due to reduced degrees of freedom and associated accelerated dynamics,^92,93 even though evaluating an ML potential requires much more operations than an empirical potential.⁹⁴ Further gains could be obtained with constrained dynamics or a coarse-grained representation of solute molecules, enabling larger integration time steps.

See the supplementary material for hyperparameters and additional results including performance and error analysis on the train dataset, absolute error vs experimental hydration free energy, correlation of molecular properties with experimental hydration free energy, potential energy difference vs Gasteiger charges, Pearson correlation of errors, and convergence of the hydration free energy computation.

This work was funded by the ERC StG under Grant No. 101077842—SupraModel. The authors thank Paul Fuchs for implementing the BAR method into the chemtrain library https://github.com/tummfm/chemtrain and for insightful discussions on theory and implementation of the ReSolv method.

The authors have no conflicts to disclose.

S.R. and A.F.B. contributed equally to this work. S.R. and A.F.B. implemented and applied the ReSolv method and conducted simulations and postprocessing. All authors planned the study, analyzed and interpreted the results, and wrote the paper.

Sebastien Röcken: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Anton F. Burnet: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Julija Zavadlav: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and its supplementary material and will be made openly available at https://github.com/tummfm/ReSolv upon acceptance of this paper. The QM7-X dataset by Hoja et al.⁶¹ is available at https://doi.org/10.5281/zenodo.4288677. The FreeSolv dataset by Mobley and Guthrie¹³ is available at https://escholarship.org/uc/item/6sd403pz. The code will be made publicly available at https://github.com/tummfm/ReSolv.git upon acceptance of this paper.