We ran 14 different combinations over 4 hyperparameters to examine the performance as a function of SMARTS search space. The 4 parameters were bond length separation α, the bit depth b, the graph depth d, and the number of new branches n. Illustrated in Fig. , we found that α had the largest impact on the bond length RMSD. Nearly all runs with α = 0.01 Å resulted in an RMSD of around 0.01 Å, whereas the RMSD of the runs with α = 0.05 were around 0.02 Å. Increasing the search space for α = 0.01 Å did not impact the RMSD, although the hierarchies were indeed different since more SMARTS patterns were found with larger search spaces. The additional SMARTS patterns found from the increased search space resulted in reduced χ 2 scores that eventually no longer contained the ideal score of 1 inside their 95% confidence intervals. The increased search space for α = 0.01 therefore led to models that were significantly worse. The deviation in the χ 2 score away from the ideal score of 1 indicated that the hierarchy either underestimated the bond length variance, or the additional SMARTS patterns did not lead to any significant improvement in RMSD. In contrast, for α = 0.05 Å we found some trends in the RMSD when the search space was increased. Specifically, the bit depth b improved the RMSD when the graph depth d = 0. The reduced χ 2 for these runs were all near 1, indicating that the hierarchies were reasonable models of the bond lengths. These hierarchies contained fewer SMARTS than the hierarchies with α = 0.01 Å. The best performing hierarchy with α = 0.01 Å contained 406 SMARTS patterns with an RMSD of 0.0095 Å, whereas Sage contained 82 SMARTS patterns with an RMSD of 0.0129 Å. The reduced χ 2 score for these were within 95% confidence intervals of one another, which indicated that the increased number of SMARTS patterns and smaller RMSD lead to performance was statistically equivalent to the Sage SMARTS patterns. Thus, a hierarchy was found that was both data driven and automatic with good performance. See SI Sec. 3.3 for a listing of the SMARTS patterns of the best performing model corresponding to b = 3 d = 0 n = 0 model including various statistics of the bond lengths per each SMARTS pattern generated.In conclusion, we found that our approach is capable of generating a SMARTS hierarchy that was competitive with the expert-derived Sage force field and is thus a good starting point for further refinement and optimization of a general force field.