As mentioned in the earlier section, a transfer learning approach is used to first train the model on the pre-training set and subsequently fine-tune the model on the fine-tuning set with fewer epochs. Within each of the 5 folds of the pre-training set, an ensemble of 5 different models are generated by using different random initialization of model parameters. This results in a total of 25 individual models from the 5-fold pre-training set. The optimized parameters from the 25 pre-trained models are then used to initialize the 25 fine-tuned models that are trained on the fine-tuning set. When evaluating the reaction-and solvent-split errors on the pre-training test set, average predictions from the ensemble of 5 models are employed to compute the error in each fold. When evaluating the final error on the experimental set, an ensemble of all 25 fine-tuned models are used to make the average predictions. We ensured that no models are trained or validated on the tested reactions or solvents. Five additional sets of features that are listed in Table are explored to improve the model performance. The RP-solv features represent the solvation free energies and solvation enthalpies of reactants and products that are calculated in this work. We also tested the Abraham parameters of solvents (S-abr) that are obtained from the ML model by Chung et al. The Abraham parameters consist of five descriptors that can describe various solute/solvent interactions. The RDKit-mol represents the 2D molecular features generated from RDKit. There are 200 2D molecular features available within RDKit, and 20 and 15 features were selected for reactants/products and solvents, respectively, based on the variance threshold and random forest methods imported from scikit-learn. Six additional RDKit atomic features (Addit-atom) that are adopted from Vermeire and Green 62 are also tested as they have shown to improve solvation free energy predictions. These include the number of radical electrons, ring size, number of lone electron pairs, H-bond donating and accepting characters, and electronegativity. Lastly, the QM atomic and bond descriptors (QM-desc) obtained from the ML model by Guan et el. are explored. The QM-desc contains 4 atomic descriptors (Hirshfeld partial charge, two Fukui indices, NMR shielding constants) and 2 bond descriptors (bond lengths, bond orders) that were shown to improve reaction property predictions. The performances of the additional features are compared using the pre-training test set. Within the ML model, additional molecular features are concatenated with the reaction and solvent embeddings and fed into the FFN layer to make the predictions. Additional atom and bond features are concatenated with the initial atom and bond features prior to the CGR/D-MPNN layers. Note that nearly all features can be calculated instantly or predicted by existing ML models. The only exception is the RP-solv features which are computed with the COSMO-RS method. Yet, several ML models are available for predicting solvation energy and enthalpy of closed-shell compounds, and the RP-solv features can be therefore estimated with the ML models if fast approximations are needed. We did not consider the 3D structures of the reactants and products as additional inputs in our study as they are usually not readily available and prone to calculation noise and error. Furthermore, Spiekermann et al. showed that the 2D D-MPNN model outperformed the 3D ML model for gas phase barrier height predictions on the Grambow et al.'s reactions. Since the same data set and similar model architecture are used in our study, we expect the result to be similar and hence do not consider the 3D ML model in this work.