From 6c4c80bc9420906dfb6d870e3cbca0be5d122f8c Mon Sep 17 00:00:00 2001 From: shashankNREL Date: Thu, 20 Aug 2026 16:14:29 -0600 Subject: [PATCH 1/8] Script to decompose molecules in to structural groups defined by Gani --- tutorials/decompose_cg.py | 670 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 670 insertions(+) create mode 100644 tutorials/decompose_cg.py diff --git a/tutorials/decompose_cg.py b/tutorials/decompose_cg.py new file mode 100644 index 0000000..861db11 --- /dev/null +++ b/tutorials/decompose_cg.py @@ -0,0 +1,670 @@ +""" +Constantinou-Gani (CG) group-contribution method decomposition for SAF-relevant +hydrocarbons. + +Reference: Constantinou & Gani, AIChE J. 40(10), 1994. + "New group contribution method for estimating properties of pure compounds" + +This tool decomposes a SMILES string into first-order and second-order CG groups. +The output is a count vector matching the 121-group column ordering used by +FuelLib/gcmTableData/gcmTable.csv. + +Scope: aliphatic + aromatic hydrocarbons (n-paraffins, iso-paraffins, terminal +alpha-olefins, alkylbenzenes, alkylnaphthalenes, monocycloparaffins, +dicycloparaffins, cycloaromatics). + +Validation: compared against existing hand-decomposed data in +FuelLib/fuelData/groupDecompositionData/refCompounds.csv +""" + +import csv +import os +import sys + +try: + from rdkit import Chem + from rdkit.Chem import rdmolops +except ImportError as e: + raise ImportError( + "RDKit is required for the CG decomposition tool. " + "Install with: pip install rdkit" + ) from e + + +class UnsupportedGroupError(ValueError): + """Raised when a molecule contains atoms or groups outside the SAF subset.""" + + +# ============================================================================= +# First-order group definitions +# ============================================================================= + +# (C, H) atom counts per first-order subgroup. +# ACCH3/ACCH2/ACCH are 2-atom groups (aromatic C + aliphatic C bundled). +FIRST_ORDER_CH = { + "CH3": (1, 3), + "CH2": (1, 2), + "CH": (1, 1), + "C": (1, 0), + "CH2=CH": (2, 3), + "CH=CH": (2, 2), + "CH2=C": (2, 2), + "CH=C": (2, 1), + "C=C": (2, 0), + "CH2=C=CH": (3, 3), + "ACH": (1, 1), + "AC": (1, 0), + "ACCH3": (2, 3), + "ACCH2": (2, 2), + "ACCH": (2, 1), +} + +# Canonical column ordering for the 121 groups (first-order + second-order) +# matching gcmTable.csv columns 2..122 (0-indexed). +# Second-order groups start at index 78 in this list. +CG_GROUP_NAMES = [ + # -- First-order groups (indices 0-77) -- + "CH3", "CH2", "CH", "C", + "CH2=CH", "CH=CH", "CH2=C", "CH=C", "C=C", "CH2=C=CH", + "ACH", "AC", "ACCH3", "ACCH2", "ACCH", + "OH", "ACOH", "CH3CO", "CH2CO", "CHO", + "CH3COO", "CH2COO", "HCOO", "CH3O", "CH2O", "CH-O", "FCH2O", + "CH2NH2", "CHNH2", "CH3NH", "CH2NH", "CHNH", "CH3N", "CH2N", + "ACNH2", "C5H4N", "C5H3N", "CH2CN", "COOH", + "CH2CL", "CHCL", "CCL", "CHCL2", "CCL2", "CCL3", "ACCL", + "CH2NO2", "CHNO2", "ACNO2", "CH2SH", + "I", "Br", "CH≡C", "C≡C", "CL—(C=C)", "ACF", + "HCON(CH2)2", "CF3", "CF2", "CF", "COO", "CCL2F", "HCCLF", "CCLF2", + "Fspecial", "CONH2", "CONHCH3", "CONHCH2", "CON(CH3)2", "CONCH3CH2", + "CON(CH2)2", "C2H5O2", "C2H4O2", "CH3S", "CH2S", "CHS", "C4H3S", "C4H2S", + # -- Second-order groups (indices 78-120) -- + "(CH3)2CH", "(CH3)3C", "CH(CH3)CH(CH3)", "CH(CH3)C(CH3)2", "C(CH3)2C(CH3)2", + "3 membered ring", "4 membered ring", "5 membered ring", + "6 membered ring", "7 membered ring", + "CHn=CHm—CHp=CHk k,n,m,p in (0,2)", + "CH3-CHm=CH, m in (0,1), n in (0,2)", + "CH2-CHm=CHn, m, n in (0,2)", + "CH-CHm=CHn or C-CHm=CHn, m,n m in (0,2)", + "Alicyclic side-chain CcyclicCm m > 1", + "CH3CH3", + "CHCHO or CCHO", "CH3COCH2", "CH3COCH or CH3COC", "Ccyclic(=0)", + "ACCHO", "CHCOOH or CCOOH", "ACCOOH", + "CH3COOCH or CH3COOC", "COCH2COO or COCHCOO or COCCOO", + " CO-O-CO", "ACCOO", + "CHOH", "COH", "CHm(OH)CHn(OH), m,n in (0,2)", + "CHm cyclic-OH, m in (0,1)", "CHm(OH)CHn(NHp), m,n,p in (0,3)", + "CHm(NH2)CHn(NH2)", "CHm cyclic-NHp-CHn cyclic, m,n,p in (0,2)", + "Chm=Chn-F, m,n in (0,2)", "AC-O-CHm", + "CHm cyclic-S-CHn cyclic, m,n in (0,2)", + "CHm=CHn—F, m,n in (0,2)", "CHm=CHn—Br, m,n in (0,2)", + "CHm=CHn—I, m,n in (0,2)", "ACBr", "ACI", + "CHm(NH2)-COOH, m,n in (0,2)", +] + +assert len(CG_GROUP_NAMES) == 121, f"Expected 121 groups, got {len(CG_GROUP_NAMES)}" + + +# ============================================================================= +# First-order decomposition +# ============================================================================= + +def _find_terminal_vinyls(mol): + """ + Locate terminal alpha-olefin vinyl groups (CH2=CH-R). + Returns (count, set_of_covered_atom_indices). + """ + covered = set() + n_vinyl = 0 + for bond in mol.GetBonds(): + if bond.GetBondTypeAsDouble() != 2.0: + continue + a, b = bond.GetBeginAtom(), bond.GetEndAtom() + if a.GetIsAromatic() or b.GetIsAromatic(): + continue + ah, bh = a.GetTotalNumHs(), b.GetTotalNumHs() + an, bn = a.GetDegree(), b.GetDegree() + # CH2=CH pattern: one end has 2H and degree 1, other has 1H and degree 2 + if (ah == 2 and an == 1) and (bh == 1 and bn == 2): + tail, head = a, b + elif (bh == 2 and bn == 1) and (ah == 1 and an == 2): + tail, head = b, a + else: + raise UnsupportedGroupError( + "C=C bond is not a terminal alpha-olefin " + f"(atoms {a.GetIdx()} H={ah} deg={an}, " + f"{b.GetIdx()} H={bh} deg={bn})." + ) + covered.add(tail.GetIdx()) + covered.add(head.GetIdx()) + n_vinyl += 1 + return n_vinyl, covered + + +def _find_aromatic_substituents(mol, excluded): + """ + Assign aromatic-substituent carbons to ACCH3/ACCH2/ACCH. + The matching aromatic ring carbon is consumed (not counted as AC later). + Returns (subgroups_dict, consumed_aromatic_set). + """ + subgroups = {} + consumed_aromatic = set() + for atom in mol.GetAtoms(): + if atom.GetIdx() in excluded: + continue + if atom.GetSymbol() != "C" or atom.GetIsAromatic(): + continue + arom_C_neighbors = [ + n for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and n.GetIsAromatic() + ] + if not arom_C_neighbors: + continue + if len(arom_C_neighbors) > 1: + raise UnsupportedGroupError( + f"Aliphatic carbon with {len(arom_C_neighbors)} aromatic neighbors " + f"(atom idx={atom.GetIdx()})." + ) + arom_C = arom_C_neighbors[0] + if arom_C.GetIdx() in consumed_aromatic: + raise UnsupportedGroupError( + f"Aromatic ring carbon (idx={arom_C.GetIdx()}) has more than one " + "aliphatic neighbor." + ) + consumed_aromatic.add(arom_C.GetIdx()) + + n_alC = sum( + 1 for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and not n.GetIsAromatic() + ) + if n_alC == 0: + subgroups[atom.GetIdx()] = "ACCH3" + elif n_alC == 1: + subgroups[atom.GetIdx()] = "ACCH2" + elif n_alC == 2: + subgroups[atom.GetIdx()] = "ACCH" + else: + raise UnsupportedGroupError( + f"Aromatic-substituent carbon with {n_alC} aliphatic neighbors " + f"(atom idx={atom.GetIdx()})." + ) + return subgroups, consumed_aromatic + + +def _classify_aliphatic_atom(atom): + """ + Assign one first-order group name (CH3/CH2/CH/C) to an aliphatic carbon. + Only for atoms not already assigned to vinyl or aromatic-substituent groups. + """ + if atom.GetSymbol() != "C": + raise UnsupportedGroupError( + f"Non-carbon atom (symbol={atom.GetSymbol()})." + ) + for bond in atom.GetBonds(): + bt = bond.GetBondTypeAsDouble() + if bt not in (1.0, 1.5): + raise UnsupportedGroupError( + f"Non-single, non-aromatic bond (order={bt}) on atom " + f"idx={atom.GetIdx()}." + ) + h_count = atom.GetTotalNumHs() + if atom.GetIsAromatic(): + if h_count == 1: + return "ACH" + if h_count == 0: + return "AC" + raise UnsupportedGroupError( + f"Aromatic carbon with H count={h_count} (atom idx={atom.GetIdx()})." + ) + # Aliphatic: classify by number of aliphatic C neighbors + aliphatic_C_neighbors = sum( + 1 for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and not n.GetIsAromatic() + ) + if aliphatic_C_neighbors == 1: + return "CH3" + if aliphatic_C_neighbors == 2: + return "CH2" + if aliphatic_C_neighbors == 3: + return "CH" + if aliphatic_C_neighbors == 4: + return "C" + raise UnsupportedGroupError( + f"Aliphatic carbon with {aliphatic_C_neighbors} C neighbors " + f"(atom idx={atom.GetIdx()})." + ) + + +def _first_order_decomposition(mol): + """ + Decompose molecule into first-order CG groups. + Returns dict mapping group name → count. + """ + counts = {} + + # 1. Terminal alpha-olefins + n_vinyl, vinyl_idxs = _find_terminal_vinyls(mol) + if n_vinyl: + counts["CH2=CH"] = n_vinyl + + # 2. Aromatic substituents (ACCH3/ACCH2/ACCH) + acch_subgroups, consumed_aromatic = _find_aromatic_substituents(mol, vinyl_idxs) + for sg_name in acch_subgroups.values(): + counts[sg_name] = counts.get(sg_name, 0) + 1 + + # 3. Classify remaining atoms + for atom in mol.GetAtoms(): + idx = atom.GetIdx() + if idx in vinyl_idxs or idx in acch_subgroups or idx in consumed_aromatic: + continue + sg = _classify_aliphatic_atom(atom) + counts[sg] = counts.get(sg, 0) + 1 + + return counts + + +# ============================================================================= +# Second-order decomposition +# ============================================================================= + +def _detect_branching_groups(mol): + """ + Detect second-order branching groups: + - (CH3)2CH: aliphatic CH with exactly 2 CH3 neighbors + - (CH3)3C: quaternary C with exactly 3 CH3 neighbors + - CH(CH3)CH(CH3): adjacent pair of CH's, each with at least 1 CH3 + - CH(CH3)C(CH3)2: adjacent CH (1 CH3) and C (2 CH3) + - C(CH3)2C(CH3)2: adjacent quaternary C's, each with 2 CH3 + """ + counts = {} + + # Classify each atom: count how many terminal CH3 neighbors it has + def _ch3_neighbor_count(atom): + """Count terminal-CH3 neighbors of an aliphatic atom.""" + n = 0 + for nbr in atom.GetNeighbors(): + if (nbr.GetSymbol() == "C" and not nbr.GetIsAromatic() + and nbr.GetTotalNumHs() == 3 and nbr.GetDegree() == 1): + n += 1 + return n + + def _is_aliphatic_C(atom): + return atom.GetSymbol() == "C" and not atom.GetIsAromatic() + + def _aliphatic_C_degree(atom): + """Number of C-C bonds (aliphatic neighbors).""" + return sum(1 for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and not n.GetIsAromatic()) + + # (CH3)2CH: CH with 2 CH3 neighbors + ch3_2_ch_atoms = set() + for atom in mol.GetAtoms(): + if not _is_aliphatic_C(atom): + continue + if atom.GetTotalNumHs() == 1 and _aliphatic_C_degree(atom) == 3: + # This is a CH (3 aliphatic C neighbors, 1 H) + if _ch3_neighbor_count(atom) >= 2: + ch3_2_ch_atoms.add(atom.GetIdx()) + if ch3_2_ch_atoms: + counts["(CH3)2CH"] = len(ch3_2_ch_atoms) + + # (CH3)3C: quaternary C with 3 CH3 neighbors + ch3_3_c_atoms = set() + for atom in mol.GetAtoms(): + if not _is_aliphatic_C(atom): + continue + if atom.GetTotalNumHs() == 0 and _aliphatic_C_degree(atom) == 4: + if _ch3_neighbor_count(atom) >= 3: + ch3_3_c_atoms.add(atom.GetIdx()) + if ch3_3_c_atoms: + counts["(CH3)3C"] = len(ch3_3_c_atoms) + + # CH(CH3)CH(CH3): adjacent CH-CH pair, each with at least 1 CH3 + ch_ch_pairs = set() + for bond in mol.GetBonds(): + a, b = bond.GetBeginAtom(), bond.GetEndAtom() + if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): + continue + # Both must be CH (1H, degree 3 aliphatic neighbors) + a_is_ch = (a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3) + b_is_ch = (b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3) + if a_is_ch and b_is_ch: + if _ch3_neighbor_count(a) >= 1 and _ch3_neighbor_count(b) >= 1: + pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) + ch_ch_pairs.add(pair) + if ch_ch_pairs: + counts["CH(CH3)CH(CH3)"] = len(ch_ch_pairs) + + # CH(CH3)C(CH3)2: adjacent CH (1 CH3) and quaternary C (2 CH3) + ch_c_pairs = set() + for bond in mol.GetBonds(): + a, b = bond.GetBeginAtom(), bond.GetEndAtom() + if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): + continue + # Check a=CH with CH3, b=C with 2 CH3 + a_is_ch = (a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3) + b_is_quat = (b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4) + if a_is_ch and b_is_quat: + if _ch3_neighbor_count(a) >= 1 and _ch3_neighbor_count(b) >= 2: + pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) + ch_c_pairs.add(pair) + # Symmetric check + b_is_ch = (b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3) + a_is_quat = (a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4) + if b_is_ch and a_is_quat: + if _ch3_neighbor_count(b) >= 1 and _ch3_neighbor_count(a) >= 2: + pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) + ch_c_pairs.add(pair) + if ch_c_pairs: + counts["CH(CH3)C(CH3)2"] = len(ch_c_pairs) + + # C(CH3)2C(CH3)2: adjacent quaternary C's, each with 2 CH3 + c_c_pairs = set() + for bond in mol.GetBonds(): + a, b = bond.GetBeginAtom(), bond.GetEndAtom() + if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): + continue + a_quat = (a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4) + b_quat = (b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4) + if a_quat and b_quat: + if _ch3_neighbor_count(a) >= 2 and _ch3_neighbor_count(b) >= 2: + pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) + c_c_pairs.add(pair) + if c_c_pairs: + counts["C(CH3)2C(CH3)2"] = len(c_c_pairs) + + return counts + + +def _detect_rings(mol): + """ + Count non-aromatic rings by size (3-7 membered). + Uses the Smallest Set of Smallest Rings (SSSR). + + A ring is counted if it is NOT fully aromatic. This handles fused + aromatic-alicyclic systems (e.g., tetralin has one fully aromatic ring + and one ring with 4 non-aromatic + 2 aromatic atoms at the junction; + only the latter is counted as a "6 membered ring"). + """ + counts = {} + ring_info = mol.GetRingInfo() + for ring in ring_info.AtomRings(): + # Skip fully aromatic rings (e.g., benzene ring in tetralin) + all_aromatic = all( + mol.GetAtomWithIdx(idx).GetIsAromatic() for idx in ring + ) + if all_aromatic: + continue + size = len(ring) + if 3 <= size <= 7: + name = f"{size} membered ring" + counts[name] = counts.get(name, 0) + 1 + return counts + + +def _detect_alicyclic_sidechain(mol): + """ + Detect alicyclic side-chain CcyclicCm (m > 1). + + NOTE: Based on empirical evidence from FuelLib refCompounds.csv, ALL + monocycloparaffins (including ethylcyclohexane, propylcyclohexane, etc.) + have this group = 0. This suggests the "alicyclic side-chain" correction + does NOT apply to simple alkyl substituents on cycloparaffin rings in the + SAF context. The exact structural requirement for this group is unclear + from the paper alone. + + Current implementation: disabled (always returns 0) to match FuelLib data. + TODO: revisit if non-SAF compounds need this correction. + """ + return 0 + + +def _detect_ch3ch3(mol): + """Detect CH3CH3 (ethane) second-order group. Only for ethane itself.""" + if mol.GetNumAtoms() == 2: + a, b = mol.GetAtomWithIdx(0), mol.GetAtomWithIdx(1) + if (a.GetSymbol() == "C" and b.GetSymbol() == "C" + and a.GetTotalNumHs() == 3 and b.GetTotalNumHs() == 3): + return 1 + return 0 + + +def _second_order_decomposition(mol): + """ + Decompose molecule into second-order CG groups. + Returns dict mapping group name → count. + """ + counts = {} + + # Branching groups + branching = _detect_branching_groups(mol) + counts.update(branching) + + # Ring corrections + rings = _detect_rings(mol) + counts.update(rings) + + # Alicyclic side-chain + n_sidechain = _detect_alicyclic_sidechain(mol) + if n_sidechain: + counts["Alicyclic side-chain CcyclicCm m > 1"] = n_sidechain + + # CH3CH3 (ethane) + n_ethane = _detect_ch3ch3(mol) + if n_ethane: + counts["CH3CH3"] = n_ethane + + return counts + + +# ============================================================================= +# Public API +# ============================================================================= + +def decompose(smiles): + """ + Decompose a hydrocarbon SMILES into CG first-order and second-order group counts. + + :param smiles: SMILES string. + :return: dict mapping group name → count (only non-zero entries). + :raises UnsupportedGroupError: If the molecule is outside the SAF subset. + :raises ValueError: If the SMILES cannot be parsed. + """ + mol = Chem.MolFromSmiles(smiles) + if mol is None: + raise ValueError(f"RDKit could not parse SMILES: {smiles!r}") + + counts = {} + + # First-order + fo = _first_order_decomposition(mol) + counts.update(fo) + + # Second-order + so = _second_order_decomposition(mol) + counts.update(so) + + return counts + + +def to_vector(counts): + """ + Convert a group-count dict to a 121-element list in canonical order. + """ + return [counts.get(name, 0) for name in CG_GROUP_NAMES] + + +def verify_formula(smiles, counts): + """ + Cross-check that first-order subgroup counts reproduce the molecular formula. + Only checks C and H from the first 15 groups (hydrocarbon groups). + """ + mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) + expected_C = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == "C") + expected_H = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == "H") + + got_C = 0 + got_H = 0 + for gname in list(FIRST_ORDER_CH.keys()): + n = counts.get(gname, 0) + c, h = FIRST_ORDER_CH[gname] + got_C += c * n + got_H += h * n + + if got_C == expected_C and got_H == expected_H: + return True, f"C{got_C}H{got_H} (matches)" + return False, f"got C{got_C}H{got_H}, expected C{expected_C}H{expected_H}" + + +# ============================================================================= +# Validation against FuelLib +# ============================================================================= + +def _load_refcompounds(): + """ + Load refCompounds.csv from FuelLib + Returns dict: compound_name → list of 121 int counts. + """ + ref_path = os.path.join( + os.path.dirname(os.path.dirname(os.path.abspath(__file__))), + "fuelData", "groupDecompositionData", "refCompounds.csv" + ) + if not os.path.exists(ref_path): + return None + + data = {} + with open(ref_path) as f: + reader = csv.reader(f) + header = next(reader) + for row in reader: + name = row[0] + # Columns 1-78 are first-order, 79-121 are second-order + # Total data columns: 121 + values = [int(row[i]) for i in range(1, min(122, len(row)))] + # Pad to 121 if needed + while len(values) < 121: + values.append(0) + data[name] = values + return data + + +def compare_with_fuellib(name, computed_vector, ref_data): + """ + Compare computed decomposition vector against FuelLib refCompounds. + Returns (match, diff_report). + """ + if ref_data is None or name not in ref_data: + return None, f"No reference data for '{name}'" + + ref_vector = ref_data[name] + diffs = [] + for i, (comp, ref) in enumerate(zip(computed_vector, ref_vector)): + if comp != ref: + diffs.append(f" {CG_GROUP_NAMES[i]}: computed={comp}, ref={ref}") + + if not diffs: + return True, "MATCH" + return False, "MISMATCH:\n" + "\n".join(diffs) + + +# ============================================================================= +# Test harness +# ============================================================================= + +if __name__ == "__main__": + # Test cases with SMILES and expected name + # Maps test name → (smiles, fuellib_name_or_None) + test_cases = [ + # n-paraffins + ("CCCCCCC", "n-heptane", "n-C07"), + ("CCCCCCCCCC", "n-decane", "n-C10"), + # iso-paraffins (2-methylalkanes) + ("CC(C)CCCC", "2-methylhexane", "C07-Isoparaffin"), + ("CC(C)CCCCC", "2-methylheptane", "C08-Isoparaffin"), + # More branched + ("CC(C)C(C)C", "2,3-dimethylbutane", None), + ("CC(C)(C)C", "neopentane", None), + # Monocycloparaffins + ("CC1CCCCC1", "methylcyclohexane", "C07-Monocycloparaffin"), + ("CCC1CCCCC1", "ethylcyclohexane", "C08-Monocycloparaffin"), + # Dicycloparaffins + ("C1CCC2CCCCC2C1", "decalin (trans)", "C10-Dicycloparaffin"), + ("C1CC[C@H]2CCCC[C@H]2C1", "cis-decalin", "C10-Dicycloparaffin"), + # Aromatics + ("Cc1ccccc1", "toluene", "Toluene"), + ("CCc1ccccc1", "ethylbenzene", "C2-Benzene"), + ("CCCc1ccccc1", "propylbenzene", "C3-Benzene"), + # Naphthalenes + ("c1ccc2ccccc2c1", "naphthalene", "Diaromatic-C10"), + ("Cc1cccc2ccccc12", "1-methylnaphthalene", "Diaromatic-C11"), + # Cycloaromatics + ("C1Cc2ccccc2C1", "indane", "Cycloaromatic-C09"), + ("C1CCc2ccccc2C1", "tetralin", "Cycloaromatic-C10"), + ("CC1CCc2ccccc2C1", "2-methyltetralin", "Cycloaromatic-C11"), + # Alkenes + ("C=CCCCCCCCCCC", "1-dodecene", "C12-Alkene"), + ] + + # Load FuelLib reference data + ref_data = _load_refcompounds() + if ref_data: + print(f"Loaded {len(ref_data)} compounds from FuelLib refCompounds.csv") + else: + print("WARNING: Could not load FuelLib refCompounds.csv for comparison") + print() + + print(f"{'Name':<22} {'SMILES':<28} {'Formula':12} {'FO Groups':<40} {'SO Groups'}") + print("=" * 130) + + n_pass = 0 + n_fail = 0 + n_skip = 0 + issues = [] + + for smi, name, ref_name in test_cases: + try: + d = decompose(smi) + ok, msg = verify_formula(smi, d) + status = "✓" if ok else "✗" + + # Separate first-order and second-order for display + fo_parts = {k: v for k, v in d.items() if k in FIRST_ORDER_CH} + so_parts = {k: v for k, v in d.items() if k not in FIRST_ORDER_CH} + + fo_str = str(fo_parts) if fo_parts else "{}" + so_str = str(so_parts) if so_parts else "{}" + + print(f"{status} {name:<20} {smi:<28} {msg:12} {fo_str:<40} {so_str}") + + if not ok: + n_fail += 1 + issues.append(f" FORMULA MISMATCH: {name} — {msg}") + continue + + # Compare with FuelLib + if ref_name and ref_data: + vec = to_vector(d) + match, report = compare_with_fuellib(ref_name, vec, ref_data) + if match is True: + n_pass += 1 + print(f" → FuelLib comparison: PASS") + elif match is False: + n_fail += 1 + print(f" → FuelLib comparison: FAIL") + print(f" {report}") + issues.append(f" FUELLIB MISMATCH: {name} vs {ref_name}") + else: + n_skip += 1 + else: + n_skip += 1 + n_pass += 1 # formula check passed + + except UnsupportedGroupError as e: + print(f"✗ {name:<20} {smi:<28} → UNSUPPORTED: {e}") + n_fail += 1 + + print() + print(f"Results: {n_pass} pass, {n_fail} fail, {n_skip} skip (no ref)") + if issues: + print("\nIssues found:") + for issue in issues: + print(issue) From 79eeb92adddf051e40c6421be6f1ab20eac4e22f Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 10:41:49 -0600 Subject: [PATCH 2/8] Format w/black --- tutorials/decompose_cg.py | 211 ++++++++++++++++++++++++++++---------- 1 file changed, 155 insertions(+), 56 deletions(-) diff --git a/tutorials/decompose_cg.py b/tutorials/decompose_cg.py index 861db11..23a5f3e 100644 --- a/tutorials/decompose_cg.py +++ b/tutorials/decompose_cg.py @@ -14,7 +14,7 @@ dicycloparaffins, cycloaromatics). Validation: compared against existing hand-decomposed data in -FuelLib/fuelData/groupDecompositionData/refCompounds.csv +FuelLib/fuelData/groupDecompositionData/refCompounds.csv """ import csv @@ -64,40 +64,127 @@ class UnsupportedGroupError(ValueError): # Second-order groups start at index 78 in this list. CG_GROUP_NAMES = [ # -- First-order groups (indices 0-77) -- - "CH3", "CH2", "CH", "C", - "CH2=CH", "CH=CH", "CH2=C", "CH=C", "C=C", "CH2=C=CH", - "ACH", "AC", "ACCH3", "ACCH2", "ACCH", - "OH", "ACOH", "CH3CO", "CH2CO", "CHO", - "CH3COO", "CH2COO", "HCOO", "CH3O", "CH2O", "CH-O", "FCH2O", - "CH2NH2", "CHNH2", "CH3NH", "CH2NH", "CHNH", "CH3N", "CH2N", - "ACNH2", "C5H4N", "C5H3N", "CH2CN", "COOH", - "CH2CL", "CHCL", "CCL", "CHCL2", "CCL2", "CCL3", "ACCL", - "CH2NO2", "CHNO2", "ACNO2", "CH2SH", - "I", "Br", "CH≡C", "C≡C", "CL—(C=C)", "ACF", - "HCON(CH2)2", "CF3", "CF2", "CF", "COO", "CCL2F", "HCCLF", "CCLF2", - "Fspecial", "CONH2", "CONHCH3", "CONHCH2", "CON(CH3)2", "CONCH3CH2", - "CON(CH2)2", "C2H5O2", "C2H4O2", "CH3S", "CH2S", "CHS", "C4H3S", "C4H2S", + "CH3", + "CH2", + "CH", + "C", + "CH2=CH", + "CH=CH", + "CH2=C", + "CH=C", + "C=C", + "CH2=C=CH", + "ACH", + "AC", + "ACCH3", + "ACCH2", + "ACCH", + "OH", + "ACOH", + "CH3CO", + "CH2CO", + "CHO", + "CH3COO", + "CH2COO", + "HCOO", + "CH3O", + "CH2O", + "CH-O", + "FCH2O", + "CH2NH2", + "CHNH2", + "CH3NH", + "CH2NH", + "CHNH", + "CH3N", + "CH2N", + "ACNH2", + "C5H4N", + "C5H3N", + "CH2CN", + "COOH", + "CH2CL", + "CHCL", + "CCL", + "CHCL2", + "CCL2", + "CCL3", + "ACCL", + "CH2NO2", + "CHNO2", + "ACNO2", + "CH2SH", + "I", + "Br", + "CH≡C", + "C≡C", + "CL—(C=C)", + "ACF", + "HCON(CH2)2", + "CF3", + "CF2", + "CF", + "COO", + "CCL2F", + "HCCLF", + "CCLF2", + "Fspecial", + "CONH2", + "CONHCH3", + "CONHCH2", + "CON(CH3)2", + "CONCH3CH2", + "CON(CH2)2", + "C2H5O2", + "C2H4O2", + "CH3S", + "CH2S", + "CHS", + "C4H3S", + "C4H2S", # -- Second-order groups (indices 78-120) -- - "(CH3)2CH", "(CH3)3C", "CH(CH3)CH(CH3)", "CH(CH3)C(CH3)2", "C(CH3)2C(CH3)2", - "3 membered ring", "4 membered ring", "5 membered ring", - "6 membered ring", "7 membered ring", + "(CH3)2CH", + "(CH3)3C", + "CH(CH3)CH(CH3)", + "CH(CH3)C(CH3)2", + "C(CH3)2C(CH3)2", + "3 membered ring", + "4 membered ring", + "5 membered ring", + "6 membered ring", + "7 membered ring", "CHn=CHm—CHp=CHk k,n,m,p in (0,2)", "CH3-CHm=CH, m in (0,1), n in (0,2)", "CH2-CHm=CHn, m, n in (0,2)", "CH-CHm=CHn or C-CHm=CHn, m,n m in (0,2)", "Alicyclic side-chain CcyclicCm m > 1", "CH3CH3", - "CHCHO or CCHO", "CH3COCH2", "CH3COCH or CH3COC", "Ccyclic(=0)", - "ACCHO", "CHCOOH or CCOOH", "ACCOOH", - "CH3COOCH or CH3COOC", "COCH2COO or COCHCOO or COCCOO", - " CO-O-CO", "ACCOO", - "CHOH", "COH", "CHm(OH)CHn(OH), m,n in (0,2)", - "CHm cyclic-OH, m in (0,1)", "CHm(OH)CHn(NHp), m,n,p in (0,3)", - "CHm(NH2)CHn(NH2)", "CHm cyclic-NHp-CHn cyclic, m,n,p in (0,2)", - "Chm=Chn-F, m,n in (0,2)", "AC-O-CHm", + "CHCHO or CCHO", + "CH3COCH2", + "CH3COCH or CH3COC", + "Ccyclic(=0)", + "ACCHO", + "CHCOOH or CCOOH", + "ACCOOH", + "CH3COOCH or CH3COOC", + "COCH2COO or COCHCOO or COCCOO", + " CO-O-CO", + "ACCOO", + "CHOH", + "COH", + "CHm(OH)CHn(OH), m,n in (0,2)", + "CHm cyclic-OH, m in (0,1)", + "CHm(OH)CHn(NHp), m,n,p in (0,3)", + "CHm(NH2)CHn(NH2)", + "CHm cyclic-NHp-CHn cyclic, m,n,p in (0,2)", + "Chm=Chn-F, m,n in (0,2)", + "AC-O-CHm", "CHm cyclic-S-CHn cyclic, m,n in (0,2)", - "CHm=CHn—F, m,n in (0,2)", "CHm=CHn—Br, m,n in (0,2)", - "CHm=CHn—I, m,n in (0,2)", "ACBr", "ACI", + "CHm=CHn—F, m,n in (0,2)", + "CHm=CHn—Br, m,n in (0,2)", + "CHm=CHn—I, m,n in (0,2)", + "ACBr", + "ACI", "CHm(NH2)-COOH, m,n in (0,2)", ] @@ -108,6 +195,7 @@ class UnsupportedGroupError(ValueError): # First-order decomposition # ============================================================================= + def _find_terminal_vinyls(mol): """ Locate terminal alpha-olefin vinyl groups (CH2=CH-R). @@ -154,8 +242,7 @@ def _find_aromatic_substituents(mol, excluded): if atom.GetSymbol() != "C" or atom.GetIsAromatic(): continue arom_C_neighbors = [ - n for n in atom.GetNeighbors() - if n.GetSymbol() == "C" and n.GetIsAromatic() + n for n in atom.GetNeighbors() if n.GetSymbol() == "C" and n.GetIsAromatic() ] if not arom_C_neighbors: continue @@ -173,7 +260,8 @@ def _find_aromatic_substituents(mol, excluded): consumed_aromatic.add(arom_C.GetIdx()) n_alC = sum( - 1 for n in atom.GetNeighbors() + 1 + for n in atom.GetNeighbors() if n.GetSymbol() == "C" and not n.GetIsAromatic() ) if n_alC == 0: @@ -196,9 +284,7 @@ def _classify_aliphatic_atom(atom): Only for atoms not already assigned to vinyl or aromatic-substituent groups. """ if atom.GetSymbol() != "C": - raise UnsupportedGroupError( - f"Non-carbon atom (symbol={atom.GetSymbol()})." - ) + raise UnsupportedGroupError(f"Non-carbon atom (symbol={atom.GetSymbol()}).") for bond in atom.GetBonds(): bt = bond.GetBondTypeAsDouble() if bt not in (1.0, 1.5): @@ -217,8 +303,7 @@ def _classify_aliphatic_atom(atom): ) # Aliphatic: classify by number of aliphatic C neighbors aliphatic_C_neighbors = sum( - 1 for n in atom.GetNeighbors() - if n.GetSymbol() == "C" and not n.GetIsAromatic() + 1 for n in atom.GetNeighbors() if n.GetSymbol() == "C" and not n.GetIsAromatic() ) if aliphatic_C_neighbors == 1: return "CH3" @@ -266,6 +351,7 @@ def _first_order_decomposition(mol): # Second-order decomposition # ============================================================================= + def _detect_branching_groups(mol): """ Detect second-order branching groups: @@ -282,8 +368,12 @@ def _ch3_neighbor_count(atom): """Count terminal-CH3 neighbors of an aliphatic atom.""" n = 0 for nbr in atom.GetNeighbors(): - if (nbr.GetSymbol() == "C" and not nbr.GetIsAromatic() - and nbr.GetTotalNumHs() == 3 and nbr.GetDegree() == 1): + if ( + nbr.GetSymbol() == "C" + and not nbr.GetIsAromatic() + and nbr.GetTotalNumHs() == 3 + and nbr.GetDegree() == 1 + ): n += 1 return n @@ -292,8 +382,11 @@ def _is_aliphatic_C(atom): def _aliphatic_C_degree(atom): """Number of C-C bonds (aliphatic neighbors).""" - return sum(1 for n in atom.GetNeighbors() - if n.GetSymbol() == "C" and not n.GetIsAromatic()) + return sum( + 1 + for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and not n.GetIsAromatic() + ) # (CH3)2CH: CH with 2 CH3 neighbors ch3_2_ch_atoms = set() @@ -325,8 +418,8 @@ def _aliphatic_C_degree(atom): if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): continue # Both must be CH (1H, degree 3 aliphatic neighbors) - a_is_ch = (a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3) - b_is_ch = (b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3) + a_is_ch = a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3 + b_is_ch = b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3 if a_is_ch and b_is_ch: if _ch3_neighbor_count(a) >= 1 and _ch3_neighbor_count(b) >= 1: pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) @@ -341,15 +434,15 @@ def _aliphatic_C_degree(atom): if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): continue # Check a=CH with CH3, b=C with 2 CH3 - a_is_ch = (a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3) - b_is_quat = (b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4) + a_is_ch = a.GetTotalNumHs() == 1 and _aliphatic_C_degree(a) == 3 + b_is_quat = b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4 if a_is_ch and b_is_quat: if _ch3_neighbor_count(a) >= 1 and _ch3_neighbor_count(b) >= 2: pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) ch_c_pairs.add(pair) # Symmetric check - b_is_ch = (b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3) - a_is_quat = (a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4) + b_is_ch = b.GetTotalNumHs() == 1 and _aliphatic_C_degree(b) == 3 + a_is_quat = a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4 if b_is_ch and a_is_quat: if _ch3_neighbor_count(b) >= 1 and _ch3_neighbor_count(a) >= 2: pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) @@ -363,8 +456,8 @@ def _aliphatic_C_degree(atom): a, b = bond.GetBeginAtom(), bond.GetEndAtom() if not (_is_aliphatic_C(a) and _is_aliphatic_C(b)): continue - a_quat = (a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4) - b_quat = (b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4) + a_quat = a.GetTotalNumHs() == 0 and _aliphatic_C_degree(a) == 4 + b_quat = b.GetTotalNumHs() == 0 and _aliphatic_C_degree(b) == 4 if a_quat and b_quat: if _ch3_neighbor_count(a) >= 2 and _ch3_neighbor_count(b) >= 2: pair = tuple(sorted([a.GetIdx(), b.GetIdx()])) @@ -389,9 +482,7 @@ def _detect_rings(mol): ring_info = mol.GetRingInfo() for ring in ring_info.AtomRings(): # Skip fully aromatic rings (e.g., benzene ring in tetralin) - all_aromatic = all( - mol.GetAtomWithIdx(idx).GetIsAromatic() for idx in ring - ) + all_aromatic = all(mol.GetAtomWithIdx(idx).GetIsAromatic() for idx in ring) if all_aromatic: continue size = len(ring) @@ -422,8 +513,12 @@ def _detect_ch3ch3(mol): """Detect CH3CH3 (ethane) second-order group. Only for ethane itself.""" if mol.GetNumAtoms() == 2: a, b = mol.GetAtomWithIdx(0), mol.GetAtomWithIdx(1) - if (a.GetSymbol() == "C" and b.GetSymbol() == "C" - and a.GetTotalNumHs() == 3 and b.GetTotalNumHs() == 3): + if ( + a.GetSymbol() == "C" + and b.GetSymbol() == "C" + and a.GetTotalNumHs() == 3 + and b.GetTotalNumHs() == 3 + ): return 1 return 0 @@ -460,6 +555,7 @@ def _second_order_decomposition(mol): # Public API # ============================================================================= + def decompose(smiles): """ Decompose a hydrocarbon SMILES into CG first-order and second-order group counts. @@ -516,17 +612,20 @@ def verify_formula(smiles, counts): # ============================================================================= -# Validation against FuelLib +# Validation against FuelLib # ============================================================================= + def _load_refcompounds(): """ - Load refCompounds.csv from FuelLib + Load refCompounds.csv from FuelLib Returns dict: compound_name → list of 121 int counts. """ ref_path = os.path.join( os.path.dirname(os.path.dirname(os.path.abspath(__file__))), - "fuelData", "groupDecompositionData", "refCompounds.csv" + "fuelData", + "groupDecompositionData", + "refCompounds.csv", ) if not os.path.exists(ref_path): return None @@ -604,7 +703,7 @@ def compare_with_fuellib(name, computed_vector, ref_data): ("C=CCCCCCCCCCC", "1-dodecene", "C12-Alkene"), ] - # Load FuelLib reference data + # Load FuelLib reference data ref_data = _load_refcompounds() if ref_data: print(f"Loaded {len(ref_data)} compounds from FuelLib refCompounds.csv") From ad2b092fc329c8db6f5cb80fdf823f25c403bd13 Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 10:56:01 -0600 Subject: [PATCH 3/8] Add rdkit and update path for refCompounds --- pyproject.toml | 1 + tutorials/decompose_cg.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index eaa67fd..53b7f26 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ dependencies = [ "scipy>=1.5.0", "pyyaml>=5.0", "matplotlib>=3.0", + "rdkit>=2022.3.1", "importlib-resources>=5.0; python_version < '3.9'", ] diff --git a/tutorials/decompose_cg.py b/tutorials/decompose_cg.py index 23a5f3e..6cad013 100644 --- a/tutorials/decompose_cg.py +++ b/tutorials/decompose_cg.py @@ -623,6 +623,8 @@ def _load_refcompounds(): """ ref_path = os.path.join( os.path.dirname(os.path.dirname(os.path.abspath(__file__))), + "fuellib", + "data", "fuelData", "groupDecompositionData", "refCompounds.csv", From b4bea089df8434b47a8aa7f1a48ffb8f5647e49a Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 10:58:16 -0600 Subject: [PATCH 4/8] Add rdkit to list of dependencies in docs --- docs/index.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/index.rst b/docs/index.rst index 132f040..1cdcba5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -56,6 +56,7 @@ FuelLib requires: - pandas ≥1.0.0 - scipy ≥1.5.0 - matplotlib ≥3.0.0 +- rdkit ≥2022.3.1 Development tools (Sphinx, Black, pytest) are available for developers installing from source; see the installation instructions in the `Contributing `_ section. From b82d2551849f30514bb91e03e9aa9ab0aea4d09f Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 11:25:59 -0600 Subject: [PATCH 5/8] Resolve codespell failure --- .codespellrc | 6 ++++++ .github/workflows/ci.yml | 2 +- pyproject.toml | 1 + tutorials/decompose_cg.py | 30 +++++++++++------------------- 4 files changed, 19 insertions(+), 20 deletions(-) create mode 100644 .codespellrc diff --git a/.codespellrc b/.codespellrc new file mode 100644 index 0000000..0393183 --- /dev/null +++ b/.codespellrc @@ -0,0 +1,6 @@ +[codespell] +# Codespell configuration for FuelLib + +skip = *.bib,*.csv,*.pdf + +ignore-words-list = mape,coo,aci \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 68ff807..974f755 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -39,7 +39,7 @@ jobs: python -m pip install --upgrade pip pip install codespell - name: Run codespell - run: codespell --skip="*.bib,*.csv,*.pdf" --ignore-words-list="mape" + run: codespell Accuracy: runs-on: ubuntu-latest diff --git a/pyproject.toml b/pyproject.toml index 53b7f26..872cd55 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,7 @@ dev = [ "sphinx>=4.0", "sphinx-rtd-theme>=1.0", "sphinxcontrib-bibtex>=2.0", + "codespell>=2.2.0", ] [project.scripts] diff --git a/tutorials/decompose_cg.py b/tutorials/decompose_cg.py index 6cad013..cf5eb9a 100644 --- a/tutorials/decompose_cg.py +++ b/tutorials/decompose_cg.py @@ -19,16 +19,8 @@ import csv import os -import sys - -try: - from rdkit import Chem - from rdkit.Chem import rdmolops -except ImportError as e: - raise ImportError( - "RDKit is required for the CG decomposition tool. " - "Install with: pip install rdkit" - ) from e +from rdkit import Chem +from rdkit.Chem import rdmolops class UnsupportedGroupError(ValueError): @@ -572,12 +564,12 @@ def decompose(smiles): counts = {} # First-order - fo = _first_order_decomposition(mol) - counts.update(fo) + first_order = _first_order_decomposition(mol) + counts.update(first_order) # Second-order - so = _second_order_decomposition(mol) - counts.update(so) + second_order = _second_order_decomposition(mol) + counts.update(second_order) return counts @@ -728,13 +720,13 @@ def compare_with_fuellib(name, computed_vector, ref_data): status = "✓" if ok else "✗" # Separate first-order and second-order for display - fo_parts = {k: v for k, v in d.items() if k in FIRST_ORDER_CH} - so_parts = {k: v for k, v in d.items() if k not in FIRST_ORDER_CH} + first_order_parts = {k: v for k, v in d.items() if k in FIRST_ORDER_CH} + second_order_parts = {k: v for k, v in d.items() if k not in FIRST_ORDER_CH} - fo_str = str(fo_parts) if fo_parts else "{}" - so_str = str(so_parts) if so_parts else "{}" + first_order_str = str(first_order_parts) if first_order_parts else "{}" + second_order_str = str(second_order_parts) if second_order_parts else "{}" - print(f"{status} {name:<20} {smi:<28} {msg:12} {fo_str:<40} {so_str}") + print(f"{status} {name:<20} {smi:<28} {msg:12} {first_order_str:<40} {second_order_str}") if not ok: n_fail += 1 From 88a9fd74445843549eb36c67664d34c3318eb403 Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 11:45:06 -0600 Subject: [PATCH 6/8] Format... --- tutorials/decompose_cg.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tutorials/decompose_cg.py b/tutorials/decompose_cg.py index cf5eb9a..9b21bb7 100644 --- a/tutorials/decompose_cg.py +++ b/tutorials/decompose_cg.py @@ -726,7 +726,9 @@ def compare_with_fuellib(name, computed_vector, ref_data): first_order_str = str(first_order_parts) if first_order_parts else "{}" second_order_str = str(second_order_parts) if second_order_parts else "{}" - print(f"{status} {name:<20} {smi:<28} {msg:12} {first_order_str:<40} {second_order_str}") + print( + f"{status} {name:<20} {smi:<28} {msg:12} {first_order_str:<40} {second_order_str}" + ) if not ok: n_fail += 1 From e13dc9363699b776f51ce4763faf363bda26f228 Mon Sep 17 00:00:00 2001 From: d-montgomery Date: Fri, 21 Aug 2026 12:14:01 -0600 Subject: [PATCH 7/8] Update version in pyproject.toml for rdkit dependency update --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 872cd55..0ede447 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "fuellib" -version = "3.0.0" +version = "3.0.1" description = "FuelLib: A Python library for Group Contribution Method (GCM) calculations of fuel properties" readme = "README.md" license = {text = "BSD-3-Clause"} From 75c674eacfb216c55ca2d4be413e6cef8ddea313 Mon Sep 17 00:00:00 2001 From: shashankNREL Date: Mon, 31 Aug 2026 09:59:24 -0600 Subject: [PATCH 8/8] Helper routines to build unifac 2.0 decompositions --- tutorials/build_unifac_decompositions.py | 255 +++++++++++++++++ tutorials/decompose_unifac.py | 331 +++++++++++++++++++++++ 2 files changed, 586 insertions(+) create mode 100644 tutorials/build_unifac_decompositions.py create mode 100644 tutorials/decompose_unifac.py diff --git a/tutorials/build_unifac_decompositions.py b/tutorials/build_unifac_decompositions.py new file mode 100644 index 0000000..bcb3e0f --- /dev/null +++ b/tutorials/build_unifac_decompositions.py @@ -0,0 +1,255 @@ +""" +Build fuelData/unifacDecomposition/.csv for every fuel in fuelData/gcData/. + +Run after editing the SMILES table or adding fuels. + +Usage: + python tools/build_unifac_decompositions.py + +The script resolves each compound's name (Reference Compound column, falling back +to Compound column for single-component fuels) to a SMILES string via NAME_TO_SMILES. +Tricyclic compounds whose PelePhysics Key column already holds a SMILES string are +used directly. Every resulting decomposition is cross-checked against the molecular +formula derived from the SMILES; mismatches abort the build. +""" + +import csv +import os +import sys + +# Ensure the FuelLib repo root is on sys.path so we can import paths.py. +THIS_FILE = os.path.abspath(__file__) +REPO_ROOT = os.path.dirname(os.path.dirname(THIS_FILE)) +if REPO_ROOT not in sys.path: + sys.path.insert(0, REPO_ROOT) +sys.path.insert(0, os.path.join(REPO_ROOT, "tools")) + +import pandas as pd +from rdkit import Chem + +import paths +from decompose_unifac import decompose, verify_formula, UnsupportedGroupError + +# --- Compound name → SMILES map --- +# +# Names match the "Reference Compound" column (or "Compound" column for single- +# component fuels). Whitespace is trimmed and case is preserved as written in the +# gcData CSVs. Where the gcData typo'd a name ("2-ethly tetralin", "2-pentyldacalin") +# the typo is preserved as the key so the lookup succeeds. +NAME_TO_SMILES = { + # --- n-paraffins (linear alkanes) --- + "n-heptane": "CCCCCCC", + "n-octane": "CCCCCCCC", + "n-nonane": "CCCCCCCCC", + "n-decane": "CCCCCCCCCC", + "n-undecane": "CCCCCCCCCCC", + "n-dodecane": "CCCCCCCCCCCC", + "n-tridecane": "CCCCCCCCCCCCC", + "n-tetradecane": "CCCCCCCCCCCCCC", + "n-pentadecane": "CCCCCCCCCCCCCCC", + "n-hexadecane (cetane)": "CCCCCCCCCCCCCCCC", + "n-heptadecane": "CCCCCCCCCCCCCCCCC", + "n-octadecane": "CCCCCCCCCCCCCCCCCC", + # Single-component fuels use a different naming convention. + "NC7H16": "CCCCCCC", + "NC10H22": "CCCCCCCCCC", + "NC12H26": "CCCCCCCCCCCC", + # --- iso-paraffins (2-methyl-alkanes) --- + "2-methyl hexane": "CC(C)CCCC", + "2-methyl heptane": "CC(C)CCCCC", + "2-methyl octane": "CC(C)CCCCCC", + "2-methyl nonane": "CC(C)CCCCCCC", + "2-methyl decane": "CC(C)CCCCCCCC", + "2-methyl undecane": "CC(C)CCCCCCCCC", + "2-methyl dodecane": "CC(C)CCCCCCCCCC", + "2-methyl tridecane": "CC(C)CCCCCCCCCCC", + "2-methyl tetradecane": "CC(C)CCCCCCCCCCCC", + "2-methyl pentadecane": "CC(C)CCCCCCCCCCCCC", + "2-methyl hexadecane": "CC(C)CCCCCCCCCCCCCC", + "2-methyl heptadecane": "CC(C)CCCCCCCCCCCCCCC", + "2-methyl octadecane": "CC(C)CCCCCCCCCCCCCCCC", + "2-methyl nonadecane": "CC(C)CCCCCCCCCCCCCCCCC", + "2-methyltricosane": "CC(C)CCCCCCCCCCCCCCCCCCCCC", # C24 + # --- alpha-olefins (terminal alkenes) --- + "1-dodecene": "C=CCCCCCCCCCC", + "1-hexadecene": "C=CCCCCCCCCCCCCCC", + # --- alkyl benzenes (C6H5-R) --- + "toluene": "Cc1ccccc1", + "ethyl benzene": "CCc1ccccc1", + "propyl benzene": "CCCc1ccccc1", + "butyl benzene": "CCCCc1ccccc1", + "pentyl benzene": "CCCCCc1ccccc1", + "hexyl benzene": "CCCCCCc1ccccc1", + "heptyl benzene": "CCCCCCCc1ccccc1", + "octyl benzene": "CCCCCCCCc1ccccc1", + "nonyl benzene": "CCCCCCCCCc1ccccc1", + # --- naphthalenes --- + "naphthalene": "c1ccc2ccccc2c1", + "1-methyl naphthalene": "Cc1cccc2ccccc12", + "1-ethyl naphthalene": "CCc1cccc2ccccc12", + "1-propyl naphthalene": "CCCc1cccc2ccccc12", + # --- alicyclic-aromatic fused (indane, tetralins) --- + "indane": "C1Cc2ccccc2C1", + "tetralin": "C1CCc2ccccc2C1", + "2-methyl tetralin": "CC1CCc2ccccc2C1", + "2-ethly tetralin": "CCC1CCc2ccccc2C1", # gcData typo: "ethly" + "2-propyl tetralin": "CCCC1CCc2ccccc2C1", + "2-butyl tetralin": "CCCCC1CCc2ccccc2C1", + # --- alkyl cyclohexanes --- + "methyl cyclohexane": "CC1CCCCC1", + "ethyl cyclohexane": "CCC1CCCCC1", + "propyl cyclohexane": "CCCC1CCCCC1", + "butyl cyclohexane": "CCCCC1CCCCC1", + "pentyl cyclohexane": "CCCCCC1CCCCC1", + "hexyl cyclohexane": "CCCCCCC1CCCCC1", + "heptyl cyclohexane": "CCCCCCCC1CCCCC1", + "octyl cyclohexane": "CCCCCCCCC1CCCCC1", + "nonyl cyclohexane": "CCCCCCCCCC1CCCCC1", + "decyl cyclohexane": "CCCCCCCCCCC1CCCCC1", + "undecyl cyclohexane": "CCCCCCCCCCCC1CCCCC1", + # --- bicyclic naphthenes (decalin family + hydrindane + octahydropentalene) --- + "Octahydropentalene": "C1CCC2CCCC12", # cis/trans-bicyclo[3.3.0]octane + "Hydrindane": "C1CCC2CCCC2C1", # cis/trans-bicyclo[4.3.0]nonane + "Decalin": "C1CCC2CCCCC2C1", # cis/trans-bicyclo[4.4.0]decane + "2-methyldecalin": "CC1CCC2CCCCC2C1", + "2-ethyldecalin": "CCC1CCC2CCCCC2C1", + "2-propyldecalin": "CCCC1CCC2CCCCC2C1", + "2-butyldecalin": "CCCCC1CCC2CCCCC2C1", + "2-pentyldacalin": "CCCCCC1CCC2CCCCC2C1", # gcData typo: "dacalin" + # --- tricyclic naphthenes (SMILES taken directly from PelePhysics Key column) --- + "C1CC2C(C1)C1CCCC21": "C1CC2C(C1)C1CCCC21", + "C1CC2CC3CCCC3C2C1": "C1CC2CC3CCCC3C2C1", + "C1CC2CC3CCCC3CC2C1": "C1CC2CC3CCCC3CC2C1", +} + + +def get_compound_name(row): + """ + Pick the most informative compound identifier from a gcData row. + + :param row: One row of the gcData CSV as a dict. + :type row: dict + :return: The lookup name (Reference Compound if present, else Compound). + :rtype: str + """ + name = row.get("Reference Compound") or row.get("Compound") or "" + return name.strip() + + +def load_subgroup_columns(): + """ + Return ordered subgroup metadata (113 entries) used to build CSV headers. + + Column headers in ``fuelData/unifacDecomposition/.csv`` are the integer + ``Subgroup_No`` values (always unique), not the human-readable + ``Subgroup_Name`` (two subgroups share the name ``CHO``, which would collide + on a ``pd.read_csv`` roundtrip). + + :return: Tuple ``(subgroup_numbers, subgroup_names)`` where ``subgroup_numbers`` + is a list of strings (header labels) and ``subgroup_names`` is the + matching list of human-readable names in the same order. + :rtype: tuple[list[str], list[str]] + """ + df = pd.read_csv(paths.UNIFAC_SUBGROUP_FILE) + numbers = [str(int(n)) for n in df["Subgroup_No"].tolist()] + names = df["Subgroup_Name"].tolist() + return numbers, names + + +def build_one_fuel(fuel_name, gc_path, out_path, subgroup_numbers, subgroup_names): + """ + Build one fuelData/unifacDecomposition/.csv from the gcData rows. + + :param fuel_name: Fuel identifier (e.g., 'heptane', 'posf10325'). + :type fuel_name: str + :param gc_path: Path to fuelData/gcData/_init.csv. + :type gc_path: str + :param out_path: Path to fuelData/unifacDecomposition/.csv to write. + :type out_path: str + :param subgroup_numbers: Ordered list of 113 ``Subgroup_No`` strings, used as + column headers (always unique). + :type subgroup_numbers: list[str] + :param subgroup_names: Ordered list of 113 subgroup names parallel to + ``subgroup_numbers`` — used to map decompose() + output (keyed by name) into the right column. + :type subgroup_names: list[str] + :return: Number of compounds decomposed. + :rtype: int + :raises KeyError: If a compound name is not in NAME_TO_SMILES. + :raises AssertionError: If formula balance fails for any compound. + """ + with open(gc_path, encoding="utf-8-sig") as fh: + rows = list(csv.DictReader(fh)) + + # Map from subgroup name → column header (subgroup number as string). + # When two subgroups share the same name (e.g., "CHO"), only the first + # occurrence's column is used by the decomposer; the second's column stays + # at 0. The classifier in decompose_unifac.py emits only the first-occurrence + # names by construction (it never generates "CHO" for ether-like groups). + name_to_header = {} + for header, name in zip(subgroup_numbers, subgroup_names): + name_to_header.setdefault(name, header) + + out_rows = [] + for r in rows: + name = get_compound_name(r) + if name not in NAME_TO_SMILES: + raise KeyError( + f"{fuel_name}: no SMILES for compound {name!r}. " + "Add it to NAME_TO_SMILES in this script." + ) + smiles = NAME_TO_SMILES[name] + counts = decompose(smiles) + ok, msg = verify_formula(smiles, counts) + assert ok, f"{fuel_name}/{name}: formula balance failed — {msg}" + compound_col = r.get("Compound", name).strip() + row_out = {header: 0 for header in subgroup_numbers} + for sg_name, n in counts.items(): + header = name_to_header[sg_name] + row_out[header] = n + row_out["Compound"] = compound_col + out_rows.append(row_out) + + out_df = pd.DataFrame(out_rows, columns=["Compound"] + subgroup_numbers) + out_df.to_csv(out_path, index=False) + return len(out_rows) + + +def main(): + """ + Build all 13 fuelData/unifacDecomposition/.csv files. + + :return: None. + :rtype: NoneType + """ + subgroup_numbers, subgroup_names = load_subgroup_columns() + assert len(subgroup_numbers) == 113, "Expected 113 UNIFAC subgroups." + + gc_dir = paths.FUELDATA_GC_DIR + out_dir = paths.FUELDATA_UNIFAC_DIR + os.makedirs(out_dir, exist_ok=True) + + fuel_files = sorted(f for f in os.listdir(gc_dir) if f.endswith("_init.csv")) + summary = [] + for fname in fuel_files: + fuel_name = fname.replace("_init.csv", "") + gc_path = os.path.join(gc_dir, fname) + out_path = os.path.join(out_dir, f"{fuel_name}.csv") + try: + n = build_one_fuel( + fuel_name, gc_path, out_path, subgroup_numbers, subgroup_names + ) + print(f" ✓ {fuel_name:<18} {n:>3} compounds → {out_path}") + summary.append((fuel_name, n, "ok")) + except (KeyError, UnsupportedGroupError, AssertionError) as e: + print(f" ✗ {fuel_name:<18} FAILED: {e}") + summary.append((fuel_name, 0, f"failed: {e}")) + + n_ok = sum(1 for _, _, s in summary if s == "ok") + print(f"\n{n_ok}/{len(summary)} fuels built successfully.") + if n_ok < len(summary): + sys.exit(1) + + +if __name__ == "__main__": + main() diff --git a/tutorials/decompose_unifac.py b/tutorials/decompose_unifac.py new file mode 100644 index 0000000..84a1395 --- /dev/null +++ b/tutorials/decompose_unifac.py @@ -0,0 +1,331 @@ +""" +UNIFAC subgroup decomposition for SAF-relevant hydrocarbons. + +helper for generating ``fuelData/unifacDecomposition/.csv`` +files. Requires RDKit (install with ``pip install rdkit``). + +Scope: aliphatic + aromatic hydrocarbons whose subgroups lie in UNIFAC main groups +1 (CH2), 2 (C=C, terminal alpha-olefins only), 3 (ACH/AC), 4 (ACCH2). These cover +all SAF jet-fuel surrogate molecules: n-paraffins, iso-paraffins, terminal +alpha-olefins, alkylbenzenes, alkylnaphthalenes, monocycloparaffins, +dicycloparaffins (decalin family), and aromatic-alicyclic fused systems (indane, +tetralin family). + +UNIFAC subgroup conventions (Magnussen 1981; same as the existing FuelLib GCM +decomposition files): +- ``ACCH3``, ``ACCH2``, ``ACCH`` are **2-atom subgroups** containing one + aromatic ring carbon AND its directly-bonded aliphatic substituent carbon + together. Their R/Q values equal AC + (CH3/CH2/CH). +- ``AC`` is **only** used for aromatic ring carbons that have no H AND no + aliphatic substituent (i.e., ring-fusion junctions like the 4a,8a positions + of naphthalene). +- ``ACH`` is each aromatic ring carbon with one H. + +A compound containing atoms or groups outside this subset raises +``UnsupportedGroupError`` so the build script flags it instead of silently +producing wrong counts. +""" + +import os +import sys + +try: + from rdkit import Chem +except ImportError as e: + raise ImportError( + "RDKit is required for the UNIFAC decomposition tool. " + "Install with: pip install rdkit" + ) from e + + +class UnsupportedGroupError(ValueError): + """Raised when a molecule contains atoms or groups outside the SAF subset.""" + + +# Subgroup number per Magnussen 1981 / UNIFAC 2.0 ``unifac_subgroups.csv``. +# Names match the ``Subgroup_Name`` column verbatim. +SUBGROUP_NUMBERS = { + "CH3": 1, + "CH2": 2, + "CH": 3, + "C": 4, + "CH2=CH": 5, # terminal alpha-olefin vinyl group + "ACH": 9, + "AC": 10, + "ACCH3": 11, + "ACCH2": 12, + "ACCH": 13, +} + + +# (C, H) atom counts per subgroup, encoding the convention above. +SUBGROUP_CH = { + "CH3": (1, 3), + "CH2": (1, 2), + "CH": (1, 1), + "C": (1, 0), + "CH2=CH": (2, 3), # =CH2 + =CH-, 2 carbons total + "ACH": (1, 1), + "AC": (1, 0), + "ACCH3": (2, 3), # AC + CH3 bundled + "ACCH2": (2, 2), # AC + CH2 bundled + "ACCH": (2, 1), # AC + CH bundled +} + + +def _find_terminal_vinyls(mol): + """ + Locate terminal alpha-olefin vinyl groups (CH2=CH-R) in a molecule. + + :param mol: RDKit molecule (implicit-H form). + :type mol: rdkit.Chem.rdchem.Mol + :return: Tuple ``(vinyl_count, covered_atom_idxs)`` — number of CH2=CH groups + found and the set of atom indices they cover. + :rtype: tuple[int, set[int]] + :raises UnsupportedGroupError: If a C=C bond is internal, di-substituted, or + otherwise not a clean terminal vinyl. + """ + covered = set() + n_vinyl = 0 + for bond in mol.GetBonds(): + if bond.GetBondTypeAsDouble() != 2.0: + continue + a, b = bond.GetBeginAtom(), bond.GetEndAtom() + if a.GetIsAromatic() or b.GetIsAromatic(): + continue # aromatic ring bond, not olefinic + ah, bh = a.GetTotalNumHs(), b.GetTotalNumHs() + an, bn = a.GetDegree(), b.GetDegree() + if (ah == 2 and an == 1) and (bh == 1 and bn == 2): + tail, head = a, b + elif (bh == 2 and bn == 1) and (ah == 1 and an == 2): + tail, head = b, a + else: + raise UnsupportedGroupError( + "C=C bond is not a terminal alpha-olefin " + f"(atoms {a.GetIdx()} H={ah} deg={an}, " + f"{b.GetIdx()} H={bh} deg={bn})." + ) + covered.add(tail.GetIdx()) + covered.add(head.GetIdx()) + n_vinyl += 1 + return n_vinyl, covered + + +def _find_aromatic_substituents(mol, excluded): + """ + Assign each aliphatic-carbon aromatic substituent its ACCH3/ACCH2/ACCH subgroup. + + The matching aromatic ring carbon is marked as "consumed" so it is not + later counted as an AC. + + :param mol: RDKit molecule (implicit-H form). + :type mol: rdkit.Chem.rdchem.Mol + :param excluded: Atom indices already assigned to other subgroups (e.g., + vinyl atoms). Atoms in this set are skipped. + :type excluded: set[int] + :return: Tuple ``(subgroups, consumed_aromatic)`` where ``subgroups`` is a + dict mapping substituent atom index → subgroup name, and + ``consumed_aromatic`` is the set of aromatic ring atom indices + that have been bundled into an ACCH_x. + :rtype: tuple[dict[int, str], set[int]] + :raises UnsupportedGroupError: If an aliphatic carbon bonds to more than one + aromatic carbon, or if the substituent has + too many aliphatic neighbors (quaternary + aromatic substituent is outside the SAF subset). + """ + subgroups = {} + consumed_aromatic = set() + for atom in mol.GetAtoms(): + if atom.GetIdx() in excluded: + continue + if atom.GetSymbol() != "C" or atom.GetIsAromatic(): + continue + arom_C_neighbors = [ + n for n in atom.GetNeighbors() if n.GetSymbol() == "C" and n.GetIsAromatic() + ] + if not arom_C_neighbors: + continue + if len(arom_C_neighbors) > 1: + raise UnsupportedGroupError( + f"Aliphatic carbon with {len(arom_C_neighbors)} aromatic neighbors " + f"(atom idx={atom.GetIdx()}). Bridge carbon between aromatic rings " + "is not in the SAF subset." + ) + arom_C = arom_C_neighbors[0] + if arom_C.GetIdx() in consumed_aromatic: + raise UnsupportedGroupError( + f"Aromatic ring carbon (idx={arom_C.GetIdx()}) has more than one " + "aliphatic neighbor — not in the SAF subset." + ) + consumed_aromatic.add(arom_C.GetIdx()) + + n_alC = sum( + 1 + for n in atom.GetNeighbors() + if n.GetSymbol() == "C" and not n.GetIsAromatic() + ) + if n_alC == 0: + subgroups[atom.GetIdx()] = "ACCH3" + elif n_alC == 1: + subgroups[atom.GetIdx()] = "ACCH2" + elif n_alC == 2: + subgroups[atom.GetIdx()] = "ACCH" + else: + raise UnsupportedGroupError( + f"Aromatic-substituent carbon with {n_alC} aliphatic neighbors " + f"(atom idx={atom.GetIdx()}). Quaternary aromatic substituent is " + "not in the SAF subset." + ) + return subgroups, consumed_aromatic + + +def _classify_atom(atom): + """ + Assign one UNIFAC subgroup name to a single atom whose subgroup is one of + ACH, AC, CH3, CH2, CH, C. Aromatic-substituent carbons (ACCH3/2/H) and vinyl + atoms are handled separately by :func:`_find_aromatic_substituents` and + :func:`_find_terminal_vinyls`; this function must not be called on them. + + :param atom: RDKit atom (must be carbon). + :type atom: rdkit.Chem.rdchem.Atom + :return: Subgroup name. + :rtype: str + :raises UnsupportedGroupError: If the atom or its bonds cannot be mapped. + """ + if atom.GetSymbol() != "C": + raise UnsupportedGroupError( + f"Non-carbon atom encountered (symbol={atom.GetSymbol()}). " + "Only hydrocarbons are supported." + ) + + for bond in atom.GetBonds(): + bt = bond.GetBondTypeAsDouble() + if bt not in (1.0, 1.5): + raise UnsupportedGroupError( + f"Non-single, non-aromatic bond (order={bt}) on atom " + f"idx={atom.GetIdx()}. Internal alkenes and alkynes are not in " + "the SAF subset." + ) + + h_count = atom.GetTotalNumHs() + + if atom.GetIsAromatic(): + if h_count == 1: + return "ACH" + if h_count == 0: + return "AC" + raise UnsupportedGroupError( + f"Aromatic carbon with H count={h_count} (atom idx={atom.GetIdx()})." + ) + + aliphatic_C_neighbors = sum( + 1 for n in atom.GetNeighbors() if n.GetSymbol() == "C" and not n.GetIsAromatic() + ) + if aliphatic_C_neighbors == 1: + return "CH3" + if aliphatic_C_neighbors == 2: + return "CH2" + if aliphatic_C_neighbors == 3: + return "CH" + if aliphatic_C_neighbors == 4: + return "C" + raise UnsupportedGroupError( + f"Aliphatic carbon with {aliphatic_C_neighbors} C neighbors " + f"(atom idx={atom.GetIdx()})." + ) + + +def decompose(smiles): + """ + Decompose a hydrocarbon SMILES into UNIFAC subgroup counts. + + :param smiles: SMILES string of the molecule. + :type smiles: str + :return: Mapping from subgroup name to integer count + (e.g., ``{"CH3": 2, "CH2": 5}``). + :rtype: dict[str, int] + :raises UnsupportedGroupError: If the molecule contains atoms or groups + outside the SAF subset (see module docstring). + :raises ValueError: If the SMILES cannot be parsed. + """ + mol = Chem.MolFromSmiles(smiles) + if mol is None: + raise ValueError(f"RDKit could not parse SMILES: {smiles!r}") + + counts = {} + + # 1. Peel off terminal alpha-olefins as CH2=CH groups. + n_vinyl, vinyl_idxs = _find_terminal_vinyls(mol) + if n_vinyl: + counts["CH2=CH"] = n_vinyl + + # 2. Assign aromatic-substituent carbons to ACCH3/ACCH2/ACCH and consume + # their aromatic ring partners. + accH_subgroups, consumed_aromatic = _find_aromatic_substituents(mol, vinyl_idxs) + for sg_name in accH_subgroups.values(): + counts[sg_name] = counts.get(sg_name, 0) + 1 + + # 3. Classify remaining atoms (ACH, AC, CH3/CH2/CH/C). + for atom in mol.GetAtoms(): + idx = atom.GetIdx() + if idx in vinyl_idxs or idx in accH_subgroups or idx in consumed_aromatic: + continue + sg = _classify_atom(atom) + counts[sg] = counts.get(sg, 0) + 1 + + return counts + + +def verify_formula(smiles, subgroup_counts): + """ + Cross-check that subgroup counts reproduce the molecular formula (C, H only). + + :param smiles: SMILES string of the molecule. + :type smiles: str + :param subgroup_counts: Mapping from subgroup name to count + (as returned by :func:`decompose`). + :type subgroup_counts: dict[str, int] + :return: Tuple ``(ok, message)`` where ``ok`` is ``True`` iff the subgroup + counts match the molecular formula derived from the SMILES. + :rtype: tuple[bool, str] + """ + mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) + expected_C = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == "C") + expected_H = sum(1 for a in mol.GetAtoms() if a.GetSymbol() == "H") + + got_C = sum(SUBGROUP_CH[sg][0] * n for sg, n in subgroup_counts.items()) + got_H = sum(SUBGROUP_CH[sg][1] * n for sg, n in subgroup_counts.items()) + + if got_C == expected_C and got_H == expected_H: + return True, f"C{got_C}H{got_H} (matches)" + return False, f"got C{got_C}H{got_H}, expected C{expected_C}H{expected_H}" + + +if __name__ == "__main__": + test_cases = [ + ("CCCCCCC", "n-heptane"), + ("CCCCCCCCCC", "n-decane"), + ("CC(C)CCCC", "2-methylhexane"), + ("Cc1ccccc1", "toluene"), + ("CCc1ccccc1", "ethylbenzene"), + ("CCCCCc1ccccc1", "pentylbenzene"), + ("CC1CCCCC1", "methylcyclohexane"), + ("c1ccc2ccccc2c1", "naphthalene"), + ("Cc1cccc2ccccc12", "1-methylnaphthalene"), + ("C1CCc2ccccc2C1", "tetralin"), + ("C1Cc2ccccc2C1", "indane"), + ("CC1CCc2ccccc2C1", "2-methyltetralin"), + ("C1CCC2CCCCC2C1", "decalin"), + ("C1CC[C@H]2CCCC[C@H]2C1", "cis-decalin"), + ("CC1CCC2CCCCC2C1", "2-methyldecalin"), + ("C=CCCCCCCCCCC", "1-dodecene"), + ("C=CCCCCCCCCCCCCCC", "1-hexadecene"), + ] + for smi, name in test_cases: + try: + d = decompose(smi) + ok, msg = verify_formula(smi, d) + status = "✓" if ok else "✗" + print(f"{status} {name:<22} {smi:<25} → {d} [{msg}]") + except UnsupportedGroupError as e: + print(f"✗ {name:<22} {smi:<25} → UNSUPPORTED: {e}")