DDBST fragmentation database

In this notebook, we evaluate the performance of ugropy against the DDBST fragmentation database. The database file was originally published in Simon Müller’s article (https://doi.org/10.1186/s13321-019-0382-3) and is hosted in the official repository: https://github.com/simonmb/fragmentation_algorithm.

To run these tests, the database file has been replicated in the ugropy repository. All rights to the original data are retained by the authors. Both this notebook and the database file are available at: https://github.com/ipqa-research/ugropy/tree/main/docs/source/tests_book/against_ddbst

The database file in the ugropy repository has been slightly modified by adding custom comments for tracking purposes; however, this notebook can be executed using the original, unmodified database file with identical results.

These comments were added to the “code” column of the file to identify specific edge cases or behaviors for certain molecules of interest. The documentation for these characters can be found at the beginning of the database file, and their meanings are defined as follows:

• ?: Disagreement with the database solution. Indicates cases where ugropy yields a different fragmentation than the DDBST solution, and the database solution appears questionable. For the purpose of this test, these molecules are still accounted for as failed solutions.

• ??: Combinatorial explosion / Slow multi-solution retrieval. ugropy fails to obtain the correct solution on its first attempt. Reaching the correct solution requires exploring all possible alternative fragmentation pathways. Since these molecules exhibit a high number of symmetric multi-solutions, the computational time required to resolve them all is significantly high. Consequently, these molecules are skipped in this test and counted as failed solutions, even though the correct solution could eventually be reached with sufficient time.

• ?!: Aromaticity mismatches. The RDKit library detects aromatic rings that DDBST considers non-aromatic, which inevitably leads to a mismatch in the group assignment. In this test, these molecules are considered failed solutions.

from ugropy import unifac

import pandas as pd

# Skip the comment lines at the beginning of the file
df = pd.read_csv("ddbst_tests.csv", index_col=1, skiprows=3)

df

	code	number	solution
smiles
CC1=C(C(=CC(=C1Cl)Cl)Cl)Cl	ZOXPZFFPPKVNEA-UHFFFAOYSA-N	13870	9:1|11:1|53:4
CCC=CCC1=C(CCC1=O)C	XMLSXPIVAXONDL-UHFFFAOYSA-N	10261	1:2|2:3|19:1|70:1|6:1
C#CCCl	LJZPPWWHKPGCHS-UHFFFAOYSA-N	12221	65:1|44:1
C1=CN=CC=C1CCN	IDLHTECVNDEOIY-UHFFFAOYSA-N	83275	2:1|29:1|38:1
C[Si](C)(C)C[Si](C)(C)CCl	WBZTYASEAQTKBK-UHFFFAOYSA-N	87574	1:5|2:1|44:1|81:2
...	...	...	...
CN(C)C1=CC=CC2=C1C=CC=C2N(C)C	TYMBOHQMUQAOBA-UHFFFAOYSA-N	139088	1:2|10:4|34:2|9:6
C=CCCCCCCCCCC=C	BPHFKBMQSYYNGQ-UHFFFAOYSA-N	30874	2:9|5:2
C=CC#N.C1=CC(=CC(=C1)CN)CN	XMEXUJUMKRUUNG-UHFFFAOYSA-N	71311341	9:4|10:2|68:1|29:2
CCN(CC)CCOC(=O)C1=C(C=C(C=C1)N)O	GHSCYMOJHVOGDJ-UHFFFAOYSA-N	68097	1:2|2:3|35:1|36:1|9:3|10:1|77:1|17:1
CC1CC(C2=CC=CC=C12)C	IIJUYSSJMAITHJ-UHFFFAOYSA-N	20143	1:2|2:1|13:2|9:4

28675 rows × 3 columns

Auxiliary function

The following function converts the DDBST fragmentation strings into a dictionary format compatible with ugropy outputs. While the DDBST database provides fragmentation results as a serialized string of group IDs and counts, ugropy returns a Python dictionary. This function parses the DDBST string into a matching dictionary to facilitate a direct comparison.

Example:

For a given DDBST string representing a molecule’s functional groups: "9:1|11:1|53:4"

The function processes it as follows:

"9:1|11:1|53:4" -> get_sol_dict -> {"ACH": 1, "ACCH3": 1, "ACCL": 4}

# Function to convert the solution string into a dictionary of group
# occurrences. E.g. 9:1|11:1|53:4 -> {'ACH': 1, 'ACCH3': 1, 'ACCL': 4}
subgroup_lookup = (
    unifac.subgroups_info
    .reset_index()
    .set_index("subgroup_number")["group"]
    .to_dict()
)

def get_sol_dict(solution: str):
    return {
        subgroup_lookup[int(g)]: int(o)
        for g, o in (group.split(":") for group in solution.split("|"))
    }

First-attempt success rate

Due to the nature of the ugropy fragmentation algorithm, certain molecules may exhibit solution multiplicity. In these cases, the first generated solution might not match the expected one, even though the correct fragmentation can be obtained by exploring all the solutions. In this specific test, we evaluate only the success rate of the first obtained solution.

The following counters and lists track the benchmarking results:

• count_total: Total number of molecules evaluated (at the end is the same as the total rows processed from the dataframe).

• count_first_correct: Number of molecules where the first solution obtained by ugropy matches the DDBST database.

• count_skipped: Number of molecules skipped due to the ?? tag (tracked as failed solutions for this test).

• smiles_incorrects: SMILES strings of the molecules that failed on the first attempt. These are stored for later analysis to determine if they can be resolved via multi-solution search.

• sol_expected_incorrects: The expected DDBST solution dictionaries for the failed cases, stored to avoid parsing the raw strings again during post-analysis.

count_total = 0                 # Total number of solutions compared
count_first_correct = 0         # Number of times the first solution matches
count_skipped = 0               # Number of solutions skipped due to "??"
smiles_incorrects = []          # List of SMILES for which the first solution is incorrect
sol_expected_incorrects = []    # List of expected solutions for the incorrect cases

for smiles, code, solution in zip(df.index, df["code"], df["solution"]):

    count_total += 1

    if "??" in code:
        count_skipped += 1
        continue

    ddbst_sol = get_sol_dict(solution)
    ugropy_sol = unifac.get_groups(smiles, "smiles").subgroups

    if ddbst_sol == ugropy_sol:
        count_first_correct += 1
    else:
        smiles_incorrects.append(smiles)
        sol_expected_incorrects.append(ddbst_sol)

Following the initial evaluation (which takes approximately 2 minutes on a standard desktop computer), we can compute the primary performance metrics for ugropy:

• count_compared: The number of molecules evaluated in the comparison (total molecules minus those skipped).

• count_incorrect: The number of molecules that failed on the first attempt (the sum of the smiles_incorrects list plus the skipped molecules).

• pct_correct: The first-attempt success rate, calculated as: count_first_correct / count_total * 100.

• pct_incorrect: The first-attempt failure rate, calculated as: count_incorrect / count_total * 100.

count_compared = count_total - count_skipped
count_incorrect = len(smiles_incorrects) + count_skipped
pct_correct = count_first_correct / count_total * 100
pct_incorrect = count_incorrect / count_total * 100

print("=" * 45)
print(f"{'UGROPY VALIDATION RESULTS':^45}")
print("=" * 45)
print(f"  Total molecules:     {count_total:>6}")
print(f"  Skipped (??):        {count_skipped:>6}")
print(f"  Compared:            {count_compared:>6}")
print("-" * 45)
print(f"  ✅ Correct:          {count_first_correct:>6}  ({pct_correct:.3f}%)")
print(f"  ❌ Incorrect:        {count_incorrect:>6}  ({pct_incorrect:.3f}%)")
print("=" * 45)

=============================================
          UGROPY VALIDATION RESULTS
=============================================
  Total molecules:      28675
  Skipped (??):            29
  Compared:             28646
---------------------------------------------
  ✅ Correct:           26567  (92.649%)
  ❌ Incorrect:          2108  (7.351%)
=============================================

Solution multiplicity success rate

From the solutions that failed on the first attempt, we can analyze how many of them can be resolved by exploring all the alternative solutions generated by ugropy. This analysis provides insight into the algorithm’s ability to eventually find the correct fragmentation, even if it is not the first choice. In this analysis, a solution is considered correct if the expected DDBST solution is found within any of the alternative paths generated by ugropy. This process takes around 1 minute on a standard desktop computer.

The following counters and lists track the results of this analysis:

• count_multiplicity_correct: The number of molecules that were initially incorrect but could be resolved by exploring all alternative solutions (additional correct solutions found).

• store_multiple_solutions: A list that stores the alternative solutions generated by ugropy where one of them was correct. This is used later to test the different multi-solution filters that ugropy provides (avoiding the need to search for them again in subsequent tests).

• smiles_multiple_correct: The SMILES strings of the molecules successfully resolved by exploring all solutions. This is also used later to test the different multi-solution filters.

• expected_multiple_correct: The expected DDBST solution dictionaries for the molecules successfully resolved by exploring all solutions. This is also stored to test the different multi-solution filters later.

• incorrects_smiles: The SMILES strings of the molecules that remain incorrect even after exploring all alternative solutions.

• incorrects_expected: The expected DDBST solution dictionaries for the molecules that remain incorrect even after exploring all alternative solutions.

count_multiplicity_correct = 0      # Counts correct solutions for molecules with solutions multiplicity
store_multiple_solutions = []       # List multiple solutions for filters testing
smiles_multiple_correct = []        # Store SMILES of molecules with multiple correct solutions
expected_multiple_correct = []      # Store expected solutions of molecules with multiple correct solutions
incorrects_smiles = []              # Store SMILES of definitely incorrect solutions
incorrects_expected = []            # Store expected solutions of definitely incorrect solutions

for smiles, expected in zip(smiles_incorrects, sol_expected_incorrects):
    solutions = unifac.get_groups(smiles, "smiles", search_multiple_solutions=True)
    sols_groups = [sol.subgroups for sol in solutions]

    if expected in sols_groups:
        count_multiplicity_correct += 1
        smiles_multiple_correct.append(smiles)
        expected_multiple_correct.append(expected)
        store_multiple_solutions.append(solutions)
    else:
        incorrects_smiles.append(smiles)
        incorrects_expected.append(expected)

After this secondary evaluation, we can analyze the following cumulative metrics:

• count_total_correct: The total number of successfully fragmented molecules, either on the first attempt or through multi-solution search, calculated as: count_first_correct + count_multiplicity_correct.

• count_definitely_incorrect: The number of molecules that remain unresolved even after exploring all alternative paths, calculated as the length of incorrects_smiles. Note that skipped molecules are excluded here since they were not processed in the multi-solution search loop.

• count_total_incorrect: The total number of failed cases, combining both unresolved molecules and skipped ones, calculated as: count_definitely_incorrect + count_skipped.

• pct_total_correct: The overall success rate considering both first-attempt matches and multi-solution resolutions, calculated as: count_total_correct / count_total * 100.

• pct_multi: The percentage of total dataset molecules that were successfully resolved specifically via multi-solution search, calculated as: count_multiplicity_correct / count_total * 100. This metric reflects the fraction of the complete dataset that required exploring alternative pathways to achieve a correct fragmentation.

count_total_correct = count_first_correct + count_multiplicity_correct
count_definitely_incorrect = len(incorrects_smiles)
count_total_incorrect = count_definitely_incorrect + count_skipped
pct_total_correct = count_total_correct / count_total * 100
pct_multi = count_multiplicity_correct / count_total * 100
pct_not_correct = count_total_incorrect / count_total * 100

print("=" * 45)
print(f"{'UGROPY VALIDATION RESULTS':^45}")
print("=" * 45)
print(f"  Total molecules:     {count_total:>6}")
print(f"  Skipped (??):        {count_skipped:>6}")
print(f"  Compared:            {count_compared:>6}")
print("-" * 45)
print(f"  ✅ Correct (first):  {count_first_correct:>9}  ({pct_correct:.3f}%)")
print(f"  ✅ Correct (multiple):  {count_multiplicity_correct:>6}  ({pct_multi:.3f}%)")
print(f"                          {'─'*19}")
print(f"  ✅ Total correct:    {count_total_correct:>9}  ({pct_total_correct:.3f}%)")
print(f"  ❌ Incorrect:        {count_total_incorrect:>9}  ({pct_not_correct:.3f}%)")
print("=" * 45)

=============================================
          UGROPY VALIDATION RESULTS
=============================================
  Total molecules:      28675
  Skipped (??):            29
  Compared:             28646
---------------------------------------------
  ✅ Correct (first):      26567  (92.649%)
  ✅ Correct (multiple):    1817  (6.337%)
                          ───────────────────
  ✅ Total correct:        28384  (98.985%)
  ❌ Incorrect:              291  (1.015%)
=============================================

Testing ugropy multi-solution filters

Now we are going to evaluate the different multi-solution filters provided by ugropy. The objective is to determine which filter performs best at selecting the correct solution from the multiple alternatives generated by the algorithm. For a detailed explanation of how each filter operates, please refer to the tutorial documentation.

The variables smiles_multiple_correct and expected_multiple_correct contain the molecules that required a multi-solution search to reach the correct fragmentation (SMILES and expected solution). Additionally, store_multiple_solutions holds the alternative solutions generated by ugropy for these specific molecules.

The next variable count_multiple_needed represents the total number of molecules that required a multi-solution search. For each filter, we will count how many of these cases are correctly resolved and compute their respective success rates.

In some cases, a filter may fail to select a unique correct solution. Naturally, those ambiguous cases are not counted as successful resolutions, but the corresponding code block can be uncommented to print them for further analysis. The other way of obtaining a fail case is when filters retain an unique solution but, of course, is the incorrect one.

count_multiple_needed = len(smiles_multiple_correct)

print(f"Number of molecules that needed multiplicity to be correct: {count_multiple_needed}")

Number of molecules that needed multiplicity to be correct: 1817

First we will evaluate the performance of the base FragmentationModel filters which are not specifically designed for a GibbsModel (spoiler: they are very bad for UNIFAC).

filter_mostly_polyatomic

success_mostly_polyatomic = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_mostly_polyatomic(sols)

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions: {[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_mostly_polyatomic += 1

print(
    "Success rate of filter_mostly_polyatomic: "
    f"{success_mostly_polyatomic}/{count_multiple_needed}"
    f" ({success_mostly_polyatomic/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_mostly_polyatomic: 0/1817 (0.00%)

filter_mostly_polarity — polar

success_mostly_polarity__polar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_mostly_polarity(sols, polarity="polar")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions: {[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_mostly_polarity__polar += 1

print(
    "Success rate of filter_mostly_polarity (polar): "
    f"{success_mostly_polarity__polar}/{count_multiple_needed}"
    f" ({success_mostly_polarity__polar/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_mostly_polarity (polar): 78/1817 (4.29%)

filter_mostly_polarity — apolar

success_mostly_polarity__apolar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_mostly_polarity(sols, polarity="apolar")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions: {[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_mostly_polarity__apolar += 1

print(
    "Success rate of filter_mostly_polarity (apolar): "
    f"{success_mostly_polarity__apolar}/{count_multiple_needed}"
    f" ({success_mostly_polarity__apolar/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_mostly_polarity (apolar): 0/1817 (0.00%)

Now we are going to evaluate the performance of the multi-solution filter specifically designed for Gibbs models.

filter_bigger_polyatomics — Q

The \(Q\) parameters (surface area fractions) of each subgroup are highly relevant because they weight the residual term of the excess Gibbs free energy models.

success_biggers_polyatomics__Q = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_bigger_polyatomics(sols, criteria="Q")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions: {[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_biggers_polyatomics__Q += 1

print(
    "Success rate of filter_bigger_polyatomics: "
    f"{success_biggers_polyatomics__Q}/{count_multiple_needed}"
    f" ({success_biggers_polyatomics__Q/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_bigger_polyatomics: 1732/1817 (95.32%)

filter_bigger_polyatomics — R

success_biggers_polyatomics__R = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_bigger_polyatomics(sols, criteria="R")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions: {[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_biggers_polyatomics__R += 1

print(
    "Success rate of filter_bigger_polyatomics: "
    f"{success_biggers_polyatomics__R}/{count_multiple_needed}"
    f" ({success_biggers_polyatomics__R/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_bigger_polyatomics: 1631/1817 (89.76%)

filter_polarity_contribution — Q — polar

Definitely the most accurate heuristic for UNIFAC implemented in ugropy at the moment.

success_polarity_contribution__Q__polar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_polarity_contribution(sols, criteria="Q", polarity="polar")

    if len(filtered) > 1:
        print(f"SMILES: {smiles}")
        print(f"Expected solution: {expected}")
        print(f"All solutions:{[sol.subgroups for sol in sols]}")
        print(f"-"*100)
        continue

    if expected == filtered[0].subgroups:
        success_polarity_contribution__Q__polar += 1
    else:
        print(f"SMILES: {smiles}")
        print(f"Expected solution: {expected}")
        print(f"Filtered solution: {filtered[0].subgroups}")
        print(f"All solutions:{[sol.subgroups for sol in sols]}")
        print(f"-"*100)

print(
    "Success rate of filter_polarity_contribution: "
    f"{success_polarity_contribution__Q__polar}/{count_multiple_needed}"
    f" ({success_polarity_contribution__Q__polar/count_multiple_needed*100:.2f}%)"
)

SMILES: CCC(C)OCOC(C)CC
Expected solution: {'CH3': 4, 'CH2': 3, 'CHO': 2}
Filtered solution: {'CH3': 4, 'CH2': 2, 'CH': 1, 'CH2O': 1, 'CHO': 1}
All solutions:[{'CH3': 4, 'CH2': 2, 'CH': 1, 'CH2O': 1, 'CHO': 1}, {'CH3': 4, 'CH2': 3, 'CHO': 2}]
----------------------------------------------------------------------------------------------------
SMILES: CN(CC1COC2=CC=CC=C2O1)C(=O)N(CC=C)CC=C
Expected solution: {'CH2': 1, 'CH2N': 1, 'CH2=CH': 2, 'ACH': 4, 'AC': 2, 'AMCH3CH2': 1, 'CHO': 1, 'THF': 1}
All solutions:[{'CH3': 1, 'CH2=CH': 2, 'ACH': 4, 'AC': 2, 'CHO': 1, 'THF': 1, 'CH2N': 1, 'AM(CH2)2': 1}, {'CH2': 1, 'CH2=CH': 2, 'ACH': 4, 'AC': 2, 'CHO': 1, 'THF': 1, 'CH2N': 1, 'AMCH3CH2': 1}, {'CH2': 1, 'CH2=CH': 2, 'ACH': 4, 'AC': 2, 'CHO': 1, 'THF': 1, 'CH3N': 1, 'AM(CH2)2': 1}]
----------------------------------------------------------------------------------------------------
Success rate of filter_polarity_contribution: 1815/1817 (99.89%)

filter_polarity_contribution — Q — apolar

success_polarity_contribution__Q__apolar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_polarity_contribution(sols, criteria="Q", polarity="apolar")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions:{[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_polarity_contribution__Q__apolar += 1


print(
    "Success rate of filter_polarity_contribution: "
    f"{success_polarity_contribution__Q__apolar}/{count_multiple_needed}"
    f" ({success_polarity_contribution__Q__apolar/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_polarity_contribution: 10/1817 (0.55%)

filter_polarity_contribution — R — polar

success_polarity_contribution__R__polar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_polarity_contribution(sols, criteria="R", polarity="polar")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions:{[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_polarity_contribution__R__polar += 1


print(
    "Success rate of filter_polarity_contribution: "
    f"{success_polarity_contribution__R__polar}/{count_multiple_needed}"
    f" ({success_polarity_contribution__R__polar/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_polarity_contribution: 1715/1817 (94.39%)

filter_polarity_contribution — R — apolar

success_polarity_contribution__R__apolar = 0

for smiles, expected, sols in zip(smiles_multiple_correct, expected_multiple_correct, store_multiple_solutions):
    filtered = unifac.filter_polarity_contribution(sols, criteria="R", polarity="apolar")

    # if len(filtered) > 1:
    #     print(f"SMILES: {smiles}")
    #     print(f"Expected solution: {expected}")
    #     print(f"All solutions:{[sol.subgroups for sol in sols]}")
    #     print(f"-"*100)
    #     continue

    if expected == filtered[0].subgroups:
        success_polarity_contribution__R__apolar += 1


print(
    "Success rate of filter_polarity_contribution: "
    f"{success_polarity_contribution__R__apolar}/{count_multiple_needed}"
    f" ({success_polarity_contribution__R__apolar/count_multiple_needed*100:.2f}%)"
)

Success rate of filter_polarity_contribution: 143/1817 (7.87%)