Repository: forlilab/Meeko (language: python)

Issue:
Allow template fetching 
This PR includes the utilities for residue template generation. **Currently it works for noncovalent cofactors only. ** 

Here's a summary of the changes: 

- Insert a check for unknown residues in the `__init__` of `LinkedRDKitChorizo`. 
If `raw_input_mols` or `set_templates` contain residues not in the current loaded residue template, a new residue template will be made. 

- Fetch from rcsb
The `input_resname` or specified resname will be used to fetch a definition CIF file from RCSB (https://files.rcsb.org/ligands/download/). The CIF files will be downloaded as temporary files. 

- Build process
The chemical component will be processed into the canonical protonation state (after deprotonating a default set of acidic protons). No further embedding, extension, capping, linker guessing will be attempted. The logging level of the build process is controlled by `ChemicalComponent_LoggingControler`. Other associated functions are imported from `chemtempgen.py`. 

- Added dependencies
**gemmi**, copy, urllib.request (and network connectivity), time, tempfile

- Example output
```
(mk_dev) amyhe@Amys-MBP Downloads % mk_prepare_receptor.py -i 1d3g_input/1d3g.pdb -o 1d3g_output -j --default_altloc A
/Users/amyhe/micromamba/envs/mk_dev/bin/mk_prepare_receptor.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('meeko==0.6.0a3')
@> 3144 atoms and 1 coordinate set(s) were parsed in 0.03s.
DEBUG:.prody:3144 atoms and 1 coordinate set(s) were parsed in 0.03s.
Input residue names {'DDQ', 'BRE', 'ORO', 'ACT', 'FMN', 'SO4'} not in residue_templates

WARNING:root:Molecule doesn't contain matching  with acidic_proton_loc {'[H][O][CX3](=O)': 0, '[H][O][SX4](=O)': 0, '[H][O][SX3](=O)': 0, '[H][O][PX4](=O)': 0, '[H][O][PX3](=O)': 0, '[H][O][CX3](=S)': 0, '[H][O][SX4](=S)': 0, '[H][SX2][a]': 0} -> deprotonate returning original mol... 
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "DDQ": {
            "smiles": "[H]C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])C([H])([H])[N+]([O-])(C([H])([H])[H])C([H])([H])[H]",
            "atom_name": ["H101", "C10", "H102", "H103", "C9", "H91", "H92", "C8", "H81", "H82", "C7", "H71", "H72", "C6", "H61", "H62", "C5", "H51", "H52", "C4", "H41", "H42", "C3", "H31", "H32", "C2", "H21", "H22", "C1", "H11", "H12", "N1", "O1", "CM1", "HM11", "HM12", "HM13", "CM2", "HM21", "HM22", "HM23"],
            "link_labels": {}
        }
    }
}
************************************************************
WARNING:root:Molecule doesn't contain matching  with acidic_proton_loc {'[H][O][CX3](=O)': 0, '[H][O][SX4](=O)': 0, '[H][O][SX3](=O)': 0, '[H][O][PX4](=O)': 0, '[H][O][PX3](=O)': 0, '[H][O][CX3](=S)': 0, '[H][O][SX4](=S)': 0, '[H][SX2][a]': 0} -> deprotonate returning original mol... 
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "ACT": {
            "smiles": "[H]C([H])([H])C(=O)[O-]",
            "atom_name": ["H1", "CH3", "H2", "H3", "C", "O", "OXT"],
            "link_labels": {}
        }
    }
}
************************************************************
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "BRE": {
            "smiles": "[H]C1=C([H])C([H])=C(C2=C([H])C([H])=C(C3=NC4=C(C([H])=C(F)C([H])=C4[H])C(C(=O)[O-])=C3C([H])([H])[H])C([H])=C2[H])C([H])=C1[H]",
            "atom_name": ["H211", "C21", "C20", "H201", "C18", "H181", "C17", "C16", "C12", "H121", "C8", "H81", "C23", "C3", "N1", "C7", "C6", "C10", "H101", "C14", "F1", "C15", "H151", "C11", "H111", "C4", "C5", "O1", "O2", "C2", "C1", "H11", "H12", "H13", "C9", "H91", "C13", "H131", "C19", "H191", "C22", "H221"],
            "link_labels": {}
        }
    }
}
************************************************************
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "ORO": {
            "smiles": "[H]C1=C(C(=O)[O-])N([H])C(=O)N([H])C1=O",
            "atom_name": ["H5", "C5", "C6", "C7", "O71", "O72", "N1", "HN1", "C2", "O2", "N3", "HN3", "C4", "O4"],
            "link_labels": {}
        }
    }
}
************************************************************
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "FMN": {
            "smiles": "[H]OC([H])(C([H])([H])OP(=O)([O-])[O-])C([H])(O[H])C([H])(O[H])C([H])([H])N1C2=NC(=O)N([H])C(=O)C2=NC2=C([H])C(C([H])([H])[H])=C(C([H])([H])[H])C([H])=C21",
            "atom_name": ["HO4'", "O4'", "C4'", "H4'", "C5'", "H5'1", "H5'2", "O5'", "P", "O1P", "O2P", "O3P", "C3'", "H3'", "O3'", "HO3'", "C2'", "H2'", "O2'", "HO2'", "C1'", "H1'1", "H1'2", "N10", "C10", "N1", "C2", "O2", "N3", "HN3", "C4", "O4", "C4A", "N5", "C5A", "C6", "H6", "C7", "C7M", "HM71", "HM72", "HM73", "C8", "C8M", "HM81", "HM82", "HM83", "C9", "H9", "C9A"],
            "link_labels": {}
        }
    }
}
************************************************************
WARNING:root:Molecule doesn't contain matching  with acidic_proton_loc {'[H][O][CX3](=O)': 0, '[H][O][SX4](=O)': 0, '[H][O][SX3](=O)': 0, '[H][O][PX4](=O)': 0, '[H][O][PX3](=O)': 0, '[H][O][CX3](=S)': 0, '[H][O][SX4](=S)': 0, '[H][SX2][a]': 0} -> deprotonate returning original mol... 
******************** New Template Built ********************
{
    "ambiguous": {},
    "residue_templates": {
        "SO4": {
            "smiles": "O=S(=O)([O-])[O-]",
            "atom_name": ["O1", "S", "O2", "O3", "O4"],
            "link_labels": {}
        }
    }
}
************************************************************

Files written:
1d3g_output.json <-- parameterized receptor

```

@diogomart Could you please take a quick look, and let me know if you like the template generation to be part of chorizo creating? If so, I have some more improvements to make and will continue working on this. If not, I will open another PR with just the check only. I can write some examples of how to make and edit chemical templates with the functions. 
Thank you for your time and kind advice in advance
Export sidechain full receptor



Discussion / hints:

Fix patch:
diff --git a/meeko/chemtempgen.py b/meeko/chemtempgen.py
new file mode 100644
index 00000000..c3421716
--- /dev/null
+++ b/meeko/chemtempgen.py
@@ -0,0 +1,769 @@
+import gemmi
+import json
+import pathlib
+import copy
+import urllib.request
+import time
+import tempfile
+import re
+import atexit
+import os
+
+from .utils.rdkitutils import covalent_radius
+
+from rdkit import Chem
+from rdkit.Chem import rdmolops
+
+from rdkit import RDLogger
+import sys, logging
+
+logger = logging.getLogger(__name__)
+
+
+list_of_AD_elements_as_AtomicNum = list(covalent_radius.keys())
+metal_AtomicNums = {12, 20, 25, 26, 30}  # Mg: 12, Ca: 20, Mn: 25, Fe: 26, Zn: 30
+
+# Utility Functions
+def mol_contains_unexpected_element(mol: Chem.Mol, allowed_elements: list[str] = list_of_AD_elements_as_AtomicNum) -> bool:
+    """Check if mol contains unexpected elements"""
+    for atom in mol.GetAtoms():
+        if atom.GetAtomicNum() not in allowed_elements:
+            return True
+    return False
+
+
+def get_atom_idx_by_names(mol: Chem.Mol, wanted_names: set[str] = set()) -> set[int]:
+    
+    if not wanted_names:
+        return set()
+    
+    wanted_atoms_by_names = {atom for atom in mol.GetAtoms() if atom.GetProp('atom_id') in wanted_names}
+    names_not_found = wanted_names.difference({atom.GetProp('atom_id') for atom in wanted_atoms_by_names})
+    if names_not_found: 
+        logger.warning(f"Molecule doesn't contain the requested atoms with names: {names_not_found} -> continue with found names... ")
+    return {atom.GetIdx() for atom in wanted_atoms_by_names}
+
+
+def get_atom_idx_by_patterns(mol: Chem.Mol, allowed_smarts: str, 
+                             wanted_smarts_loc: dict[str, set[int]] = None,
+                             allow_multiple: bool=False) -> set[int]:
+    
+    if wanted_smarts_loc is None:
+        return set()
+    
+    wanted_atoms_idx = set()
+    
+    tmol = Chem.MolFromSmarts(allowed_smarts)
+    match_allowed = mol.GetSubstructMatches(tmol)
+    if not match_allowed:
+        logger.warning(f"Molecule doesn't contain allowed_smarts: {allowed_smarts} -> no pattern-based action will be made. ")
+        return set()
+    if len(match_allowed) > 1 and not allow_multiple: 
+        logger.warning(f"Molecule contain multiple copies of allowed_smarts: {allowed_smarts} -> no pattern-based action will be made. ")
+        return set()
+    if len(match_allowed) > 1 and allow_multiple:
+        match_allowed = {item for sublist in match_allowed for item in sublist}
+    else:
+        match_allowed = match_allowed[0]
+    
+    atoms_in_mol = [atom for atom in mol.GetAtoms()]
+    for wanted_smarts in wanted_smarts_loc: 
+        lmol = Chem.MolFromSmarts(wanted_smarts)
+        match_wanted = mol.GetSubstructMatches(lmol)
+        if not match_wanted:
+            logger.warning(f"Molecule doesn't contain wanted_smarts: {wanted_smarts} -> continue with next pattern... ")
+            continue
+        for match_copy in match_wanted:
+            match_in_copy = (idx for idx in match_copy if match_copy.index(idx) in wanted_smarts_loc[wanted_smarts])
+            match_wanted_atoms = {atoms_in_mol[idx] for idx in match_in_copy if idx in match_allowed}
+            if match_wanted_atoms: 
+                wanted_atoms_idx.update(atom.GetIdx() for atom in match_wanted_atoms)
+    
+    return wanted_atoms_idx
+
+
+# Mol Editing Functions
+def embed(mol: Chem.Mol, allowed_smarts: str, 
+          leaving_names: set[str] = None, leaving_smarts_loc: dict[str, set[int]] = None, 
+          alsoHs: bool = True) -> Chem.Mol:
+    """
+    Remove atoms from the molecule based the union of
+    (a) leaving_names: list of atom IDs (names), and
+    (b) leaving_smarts_loc: dict to map substructure SMARTS patterns with 
+    tuple of 0-based indicies for atoms to delete (restricted by allowed_smarts)
+    """
+    if leaving_names is None and leaving_smarts_loc is None:
+        return mol
+
+    leaving_atoms_idx = set()
+
+    if leaving_names:
+        leaving_atoms_idx.update(get_atom_idx_by_names(mol, leaving_names))
+
+    if leaving_smarts_loc:
+        leaving_atoms_idx.update(get_atom_idx_by_patterns(mol, allowed_smarts, leaving_smarts_loc))
+
+    if leaving_atoms_idx and alsoHs:
+        atoms_in_mol = [atom for atom in mol.GetAtoms()]
+        leaving_Hs = (ne for atom_idx in leaving_atoms_idx.copy() for ne in atoms_in_mol[atom_idx].GetNeighbors() if ne.GetAtomicNum() == 1)
+        leaving_atoms_idx.update(atom.GetIdx() for atom in leaving_Hs)
+
+    if not leaving_atoms_idx:
+        logger.warning(f"No matched atoms to delete -> embed returning original mol...")
+        return mol
+    
+    rwmol = Chem.RWMol(mol)
+    for atom_idx in sorted(leaving_atoms_idx, reverse=True): 
+        rwmol.RemoveAtom(atom_idx)
+    rwmol.UpdatePropertyCache()
+    return rwmol.GetMol()
+
+
+def cap(mol: Chem.Mol, allowed_smarts: str, 
+        capping_names: set[str] = None, capping_smarts_loc: dict[str, set[int]] = None, 
+        protonate: bool = False) -> Chem.Mol:
+    """Add hydrogens to atoms with implicit hydrogens based on the union of
+    (a) capping_names: list of atom IDs (names), and
+    (b) capping_smarts_loc: dict to map substructure SMARTS patterns with 
+    tuple of 0-based indicies for atoms."""
+   
+    if capping_names is None and capping_smarts_loc is None:
+        return mol
+    
+    capping_atoms_idx = set()
+    
+    if capping_names:
+        capping_atoms_idx.update(get_atom_idx_by_names(mol, capping_names))
+
+    if capping_smarts_loc:
+        capping_atoms_idx.update(get_atom_idx_by_patterns(mol, allowed_smarts, capping_smarts_loc, allow_multiple = True))
+
+    if not capping_atoms_idx:
+        logger.warning(f"No matched atoms to cap -> cap returning original mol...")
+        return mol
+    
+    def get_max_Hid(mol: Chem.Mol) -> int:
+        all_Hids = (atom.GetProp('atom_id') for atom in mol.GetAtoms() if atom.GetAtomicNum()==1)
+        regular_ids = {Hid for Hid in all_Hids if Hid[0]=='H' and Hid[1:].isdigit()}
+        if len(regular_ids) > 0:
+            return max(int(x[1:]) for x in regular_ids)
+        else:
+            return 0
+    
+    rwmol = Chem.RWMol(mol)
+    new_Hid = get_max_Hid(mol) + 1
+    atoms_in_mol = [atom for atom in mol.GetAtoms()]
+    for atom_idx in capping_atoms_idx:
+        parent_atom  = atoms_in_mol[atom_idx]
+        if protonate: 
+            parent_atom.SetFormalCharge(parent_atom.GetFormalCharge() + 1)
+        needed_Hs = parent_atom.GetNumImplicitHs()
+        if needed_Hs == 0:
+            logger.warning(f"Atom # {atom_idx} ({parent_atom.GetProp('atom_id')}) in mol doesn't have implicit Hs -> continue with next atom... ")
+        else:
+            new_atom = Chem.Atom("H")
+            new_atom.SetProp('atom_id', f"H{new_Hid}")
+            new_Hid += 1
+            new_idx = rwmol.AddAtom(new_atom)
+            rwmol.AddBond(atom_idx, new_idx, Chem.BondType.SINGLE)
+    rwmol.UpdatePropertyCache()
+    return rwmol.GetMol()
+
+
+def deprotonate(mol: Chem.Mol, acidic_proton_loc: dict[str, int] = None) -> Chem.Mol:
+    """Remove acidic protons from the molecule based on acidic_proton_loc"""
+    # acidic_proton_loc is a mapping 
+    # keys: smarts pattern of a fragment
+    # value: the index (order in smarts) of the leaving proton
+
+    if acidic_proton_loc is None:
+        return mol
+
+    # deprotonate all matched protons
+    acidic_protons_idx = set()
+    for smarts_pattern, idx in acidic_proton_loc.items():
+        qmol = Chem.MolFromSmarts(smarts_pattern)
+        acidic_protons_idx.update(match[idx] for match in mol.GetSubstructMatches(qmol))
+    
+    if not acidic_protons_idx:
+        logger.warning(f"Molecule doesn't contain matching atoms for acidic_proton_loc:" + 
+                        f"{acidic_proton_loc}" + 
+                        f"-> deprotonate returning original mol... ")
+        return mol
+    
+    rwmol = Chem.RWMol(mol)
+    for atom_idx in sorted(acidic_protons_idx, reverse=True):
+        rwmol.RemoveAtom(atom_idx)
+        neighbors = mol.GetAtomWithIdx(atom_idx).GetNeighbors()
+        if neighbors:
+            neighbor_atom = rwmol.GetAtomWithIdx(neighbors[0].GetIdx())
+            neighbor_atom.SetFormalCharge(neighbor_atom.GetFormalCharge() - 1)
+    
+    rwmol.UpdatePropertyCache()
+    return rwmol.GetMol()
+
+
+def recharge(rwmol: Chem.RWMol) -> Chem.RWMol: 
+    # Recharging metal complexes
+    metal_atoms = set(atom for atom in rwmol.GetAtoms() if atom.GetAtomicNum() in metal_AtomicNums)
+    coord_valence = {v: k for k, v_set in {
+            1: {1, 9, 17, 35, 53}, # monovalent: H, halogens
+            2: {8, 16}, # divalent: O, S
+            3: {5, 7, 15}, # trivalent: B, N, P (phosphine)
+            4: {6, 14}, # tetravalent: C, Si
+            }.items() for v in v_set}
+    bond_order_mapping = {
+        Chem.BondType.SINGLE: 1,
+        Chem.BondType.DOUBLE: 2,
+        Chem.BondType.TRIPLE: 3,
+        Chem.BondType.AROMATIC: 1.5  
+    }
+
+    if not metal_atoms: 
+        return rwmol
+    else:
+        # Method a: If charge is ignored at metal center
+        # -> charge up nonmetal coordinated atoms 
+        # (metal ox state will be interpreted from valence of nonmetal coordinated atom)
+        metal_to_nonmetal_neighbors = {metal: {atom for atom in metal.GetNeighbors() 
+                                               if atom.GetAtomicNum() not in metal_AtomicNums} for metal in metal_atoms}
+        nonmetal_coord_atoms = set(atom for v_set in metal_to_nonmetal_neighbors.values() for atom in v_set)
+
+        if any(atom.GetFormalCharge() == 0 for atom in metal_atoms): 
+            logger.warning(f"Molecule contains metal with unspecified charge state -> charging nonmetal coordinated atoms...")
+            logger.warning(f"All metals will be neutralized and the nonmetal coordinated atoms will be charged according to their explicit valence. ")
+            for atom in metal_atoms: 
+                atom.SetFormalCharge(0)
+            for atom in nonmetal_coord_atoms: 
+                bond_sum = sum(bond_order_mapping.get(bond.GetBondType(), 0) for bond in atom.GetBonds())
+                atom.SetFormalCharge(bond_sum - coord_valence[atom.GetAtomicNum()])
+  
+        # Method b: If charge (ox) state is specified for all metal centers
+        # -> ionize (break bonds w/) and charge down nonmetal coordinated atoms
+        # (metal ox state will be from input charge state)
+        else: 
+            logger.warning(f"Molecule contains metal specified charge state -> ionizing metal-nonmetal coordinate bonds...")
+            logger.warning(f"All metal-nonmetal coordinate bonds will be polarized. ")
+            metal_bonds = {bond.GetIdx(): (bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()) 
+                               for metal in metal_atoms for bond in metal.GetBonds()}
+
+            for bond_idx in sorted(metal_bonds, reverse=True):
+                rwmol.RemoveBond(metal_bonds[bond_idx][0], metal_bonds[bond_idx][1])
+
+            for atom in nonmetal_coord_atoms: 
+                bond_sum = sum(bond_order_mapping.get(bond.GetBondType(), 0) for bond in atom.GetBonds())
+                atom.SetFormalCharge(bond_sum - coord_valence[atom.GetAtomicNum()])
+
+            # to avoid fragmentation, rebuild all metal-nonmetal bonds
+            for metal in metal_atoms: 
+                for nei in metal_to_nonmetal_neighbors[metal]:
+                    rwmol.AddBond(metal.GetIdx(), nei.GetIdx(), Chem.BondType.SINGLE)
+                    metal.SetFormalCharge(metal.GetFormalCharge() - 1)
+                    nei.SetFormalCharge(nei.GetFormalCharge() + 1)
+
+        logger.warning(f"Total charge of the molecule after recharging: {sum(atom.GetFormalCharge() for atom in rwmol.GetAtoms())}")
+
+    return rwmol
+
+
+# Attribute Formatters
+def get_smiles_with_atom_names(mol: Chem.Mol) -> tuple[str, list[str]]:
+    """Generate SMILES with atom names in the order of SMILES output."""
+    # allHsExplicit may expose the implicit Hs of linker atoms to Smiles; the implicit Hs don't have names
+    smiles_exh = Chem.MolToSmiles(mol, allHsExplicit=True)
+
+    smiles_atom_output_order = mol.GetProp('_smilesAtomOutputOrder')
+    delimiters = ('[', ']', ',')
+    for delimiter in delimiters:
+        smiles_atom_output_order = smiles_atom_output_order.replace(delimiter, ' ')
+    smiles_output_order = (int(x) for x in smiles_atom_output_order.split())
+
+    atom_id_dict = {atom.GetIdx(): atom.GetProp('atom_id') for atom in mol.GetAtoms()}
+    atom_name = [atom_id_dict[atom_i] for atom_i in smiles_output_order]
+
+    return smiles_exh, atom_name
+
+
+def get_pretty_smiles(smi: str) -> str: 
+    """Convert Smiles with allHsExplicit to pretty Smiles to be put on chem templates"""
+    # collect the inside square brackets contents
+    contents = set(re.findall(r'\[([^\]]+)\]', smi))
+
+    def is_nonmetal_element(symbol: str) -> bool:
+        """Check if a string represents a valid chemical element."""
+        try:
+            return Chem.GetPeriodicTable().GetAtomicNumber(symbol) not in metal_AtomicNums
+        # rdkit throws RuntimeError if invalid
+        except RuntimeError:
+            return False
+
+    for content in contents:
+        # keep [H] for explicit Hs
+        if content == 'H': 
+            continue
+        # drop H in the content to hide implicit Hs
+        H_stripped = content.split('H')[0]
+        # drop [ ] if the content is an uncharged nonmetal element symbol 
+        if is_nonmetal_element(content) or is_nonmetal_element(H_stripped):
+            smi = smi.replace(f"[{content}]", f"{H_stripped}" if 'H' in content else f"{content}")
+    return smi
+
+class ChemicalComponent_LoggingControler:
+
+    def __init__(self):
+        self.original_level = logger.level
+        self.rdkit_logger = RDLogger.logger()
+        self.default_rdkit_level = RDLogger.WARNING
+        self.handler = None
+
+    def __enter__(self):
+        self.rdkit_logger.setLevel(RDLogger.CRITICAL)
+        logger.setLevel(logging.WARNING)
+        self.handler = logging.StreamHandler(sys.stdout)
+        logger.addHandler(self.handler)
+        return self
+    
+    def __exit__(self, exc_type, exc_value, traceback):
+        self.rdkit_logger.setLevel(self.default_rdkit_level)
+        logger.setLevel(self.original_level)
+        logger.removeHandler(self.handler)
+
+
+class ChemicalComponent:
+
+    def __init__(self, rdkit_mol: Chem.Mol, resname: str, smiles_exh: str, atom_name: list[str]):
+        self.rdkit_mol = rdkit_mol
+        self.resname = resname
+        self.parent = resname # default parent to itself
+        self.smiles_exh = smiles_exh
+        self.atom_name = atom_name
+        self.link_labels = {} # default to empty dict (free molecular form)
+
+    def __eq__(self, other):
+        if isinstance(other, ChemicalComponent):
+            return self.smiles_exh == other.smiles_exh and self.atom_name == other.atom_name
+        return False
+
+    @classmethod
+    def from_cif(cls, source_cif: str, resname: str):
+        """Create ChemicalComponent from a chemical component dict file and a searchable residue name in file."""
+        
+        # Locate block by resname
+        doc = gemmi.cif.read(source_cif)
+        block = doc.find_block(resname)
+        
+        # Populate atom table
+        atom_category = '_chem_comp_atom.'
+        atom_attributes = ['atom_id', # atom names
+                           'type_symbol', # atom elements
+                           'charge', # (atomic) formal charges
+                           'pdbx_leaving_atom_flag', # flags for leaving atoms after polymerization
+                           ]
+        atom_table = block.find(atom_category, atom_attributes)
+        atom_cols = {attr: atom_table.find_column(f"{atom_category}{attr}") for attr in atom_attributes}
+
+        # Summon rdkit atoms into empty RWMol
+        rwmol = Chem.RWMol()
+        atom_elements = atom_cols['type_symbol']
+
+        for idx, element in enumerate(atom_elements):
+            if len(element)==2:
+                element = element[0] + element[1].lower()
+            rdkit_atom = Chem.Atom(element)
+            for attr in atom_attributes:
+                rdkit_atom.SetProp(attr, atom_cols[attr][idx])
+                # strip double quotes in names
+                raw_name = rdkit_atom.GetProp('atom_id')
+                rdkit_atom.SetProp('atom_id', raw_name.strip('"'))
+            target_charge = atom_cols['charge'][idx]
+            if target_charge!='0':
+                rdkit_atom.SetFormalCharge(int(target_charge)) # this needs to be int for rdkit
+            rwmol.AddAtom(rdkit_atom)
+
+        # Check if rwmol contains unexpected elements
+        if mol_contains_unexpected_element(rwmol):
+            logger.warning(f"Molecule contains unexpected elements -> template for {resname} will be None. ")
+            return None
+
+        # Map atom_id (atom names) with rdkit idx
+        name_to_idx_mapping = {atom.GetProp('atom_id'): idx for (idx, atom) in enumerate(rwmol.GetAtoms())}
+
+        # Populate bond table
+        bond_category = '_chem_comp_bond.'
+        bond_attributes = ['atom_id_1', # atom name 1
+                           'atom_id_2', # atom name 2
+                           'value_order', # bond order
+                           ]
+        bond_table = block.find(bond_category, bond_attributes)
+        bond_cols = {attr: bond_table.find_column(f"{bond_category}{attr}") for attr in bond_attributes}
+
+        # Connect atoms by bonds
+        bond_type_mapping = {
+            'SING': Chem.BondType.SINGLE,
+            'DOUB': Chem.BondType.DOUBLE,
+            'TRIP': Chem.BondType.TRIPLE,
+            'AROM': Chem.BondType.AROMATIC
+        }
+        bond_types = bond_cols['value_order']
+
+        for bond_i, bond_type in enumerate(bond_types):
+            rwmol.AddBond(name_to_idx_mapping[bond_cols['atom_id_1'][bond_i].strip('"')], 
+                          name_to_idx_mapping[bond_cols['atom_id_2'][bond_i].strip('"')], 
+                          bond_type_mapping.get(bond_type, Chem.BondType.UNSPECIFIED))
+
+        # Try recharging mol (for metals)
+        try:
+            rwmol = recharge(rwmol)
+        except Exception as e:
+            logger.error(f"Failed to recharge rdkitmol. Error: {e} -> template for {resname} will be None. ")
+            return None
+
+        # Finish eidting mol 
+        try:    
+            rwmol.UpdatePropertyCache()
+        except Exception as e:
+            logger.error(f"Failed to create rdkitmol from cif. Error: {e} -> template for {resname} will be None. ")
+            return None
+        
+        # Check implicit Hs
+        total_implicit_hydrogens = sum(atom.GetNumImplicitHs() for atom in rwmol.GetAtoms())
+        if total_implicit_hydrogens > 0:
+            logger.error(f"rdkitmol from cif has implicit hydrogens. -> template for {resname} will be None. ")
+            return None
+
+        rdkit_mol = rwmol.GetMol()
+            
+        # Get Smiles with explicit Hs and ordered atom names
+        smiles_exh, atom_name = get_smiles_with_atom_names(rdkit_mol)
+        
+        return cls(rdkit_mol, resname, smiles_exh, atom_name)
+
+
+    def make_canonical(self, acidic_proton_loc):
+        """Deprotonate acidic groups til the canonical (most deprotonated) state."""
+        self.rdkit_mol = deprotonate(self.rdkit_mol, acidic_proton_loc = acidic_proton_loc)
+        return self
+
+    def make_embedded(self, allowed_smarts, leaving_names = None, leaving_smarts_loc = None):
+        """Remove leaving atoms from the molecule by atom names and/or patterns."""
+        self.rdkit_mol = embed(self.rdkit_mol, allowed_smarts = allowed_smarts, 
+                               leaving_names = leaving_names, leaving_smarts_loc = leaving_smarts_loc)
+        return self
+        
+    def make_capped(self, allowed_smarts, capping_names = None, capping_smarts_loc = None, protonate = None):
+        """Build and name explicit hydrogens for atoms with implicit Hs by atom names and/or patterns."""
+        self.rdkit_mol = cap(self.rdkit_mol, allowed_smarts = allowed_smarts, 
+                             capping_names = capping_names, capping_smarts_loc = capping_smarts_loc,
+                             protonate = protonate)
+        return self
+        
+    def make_pretty_smiles(self):
+        """Build and name explicit hydrogens for atoms with implicit Hs by atom names and/or patterns."""
+        self.smiles_exh, self.atom_name = get_smiles_with_atom_names(self.rdkit_mol)
+        self.smiles_exh = get_pretty_smiles(self.smiles_exh)
+        return self
+
+    def make_link_labels_from_patterns(self, pattern_to_label_mapping = {}):
+        """Map patterns to link labels based on a given mapping."""
+        if not pattern_to_label_mapping:
+            return self
+
+        for pattern in pattern_to_label_mapping:
+            atom_idx = get_atom_idx_by_patterns(self.rdkit_mol, allowed_smarts = Chem.MolToSmarts(self.rdkit_mol), 
+                                                wanted_smarts_loc = {pattern: {0}})
+            atoms_in_mol = [atom for atom in self.rdkit_mol.GetAtoms()]
+            if not atom_idx:
+                logger.warning(f"Molecule doesn't contain pattern: {pattern} -> linker label for {pattern_to_label_mapping[pattern]} will not be made. ")
+            elif len(atom_idx) > 1:
+                logger.warning(f"Molecule contain multiple copies of pattern: {pattern} -> linker label for {pattern_to_label_mapping[pattern]} will not be made. ")
+            else:
+                atom_idx = next(iter(atom_idx))
+                name = atoms_in_mol[atom_idx].GetProp('atom_id')
+                self.link_labels.update({str(self.atom_name.index(name)): pattern_to_label_mapping[pattern]})
+
+        return self
+    
+    def make_link_labels_from_names(self, name_to_label_mapping = {}):
+        """Map atom names to link labels based on a given mapping."""
+        if not name_to_label_mapping:
+            return self
+
+        for atom in self.rdkit_mol.GetAtoms():
+            if atom.GetProp('atom_id') in name_to_label_mapping:
+                if atom.GetNumImplicitHs() > 0:
+                    name = atom.GetProp('atom_id')
+                    self.link_labels.update({str(self.atom_name.index(name)): name_to_label_mapping[name]})
+
+        return self
+    
+    def ResidueTemplate_check(self) -> bool:
+        # ResidueTemplate.check from linked_rdkit_chorizo
+        ps = Chem.SmilesParserParams()
+        ps.removeHs = False
+        mol = Chem.MolFromSmiles(self.smiles_exh, ps)
+        have_implicit_hs = {atom.GetIdx() for atom in mol.GetAtoms() if atom.GetTotalNumHs() > 0}
+        
+        if self.link_labels: 
+            int_labels = {int(atom_id) for atom_id in self.link_labels}
+        else:
+            int_labels = set()
+
+        if int_labels != have_implicit_hs:
+            raise ValueError(
+                f"expected any atom with non-real Hs ({have_implicit_hs}) to tagged in link_labels ({int_labels})"
+            )
+        
+        if not self.atom_name: 
+            return
+        
+        if len(self.atom_name) != mol.GetNumAtoms():
+            raise ValueError(f"{len(self.atom_name)=} differs from {mol.GetNumAtoms()=}")
+        return
+
+
+# Export/Writer Function
+def export_chem_templates_to_json(cc_list: list[ChemicalComponent], json_fname: str=""):
+    """Export list of chem templates to json"""
+
+    basenames = [cc.parent for cc in cc_list]
+    ambiguous_dict = {basename: [] for basename in basenames}
+    for cc in cc_list:
+        if cc.resname not in ambiguous_dict[cc.parent]:
+            ambiguous_dict[cc.parent].append(cc.resname)
+
+    data_to_export = {"ambiguous": {basename: basename+".ambiguous" for basename in basenames}}
+
+    residue_templates = {}
+    for cc in cc_list:
+        residue_templates[cc.resname] = {
+            "smiles": cc.smiles_exh,
+            "atom_name": cc.resname+".atom_names",
+        }
+        if cc.link_labels:
+            residue_templates[cc.resname]["link_labels"] = cc.resname+".link_labels"
+        else:
+            residue_templates[cc.resname]["link_labels"] = {}
+
+    data_to_export.update({"residue_templates": residue_templates})
+
+    json_str = json.dumps(data_to_export, indent = 4)
+
+    # format ambiguous resnames to one line
+    for basename in data_to_export["ambiguous"]:
+        single_line_resnames = json.dumps(ambiguous_dict[basename], separators=(', ', ': '))
+        json_str = json_str.replace(json.dumps(data_to_export["ambiguous"][basename], indent = 4), single_line_resnames)
+
+    # format link_labels and atom_name to one line
+    for cc in cc_list:
+        single_line_atom_name = json.dumps(cc.atom_name, separators=(', ', ': '))
+        json_str = json_str.replace(json.dumps(data_to_export["residue_templates"][cc.resname]["atom_name"], indent = 4), single_line_atom_name)
+        if cc.link_labels:
+            single_line_link_labels = json.dumps(cc.link_labels, separators=(', ', ': '))
+            json_str = json_str.replace(json.dumps(data_to_export["residue_templates"][cc.resname]["link_labels"], indent = 4), single_line_link_labels)
+
+    if json_fname:
+        with open(pathlib.Path(json_fname), 'w') as f:
+            f.write(json_str)
+        print(f"{json_fname} <-- Json File for New Chemical Templates")
+    else:
+        print(" New Template Built ".center(60, "*"))
+        print(json_str)
+        print("*"*60)
+        return json_str
+
+
+def fetch_from_pdb(resname: str, max_retries = 5, backoff_factor = 2) -> str: 
+    
+    url = f"https://files.rcsb.org/ligands/download/{resname}.cif"
+    file_path = f"{resname}.cif"
+    for retry in range(max_retries):
+        try:
+            with urllib.request.urlopen(url) as response:
+                content = response.read()
+            logger.info(f"File downloaded successfully: {file_path}")
+            with tempfile.NamedTemporaryFile(delete=False, suffix=".cif") as temp_file:
+                temp_file.write(content)
+            atexit.register(
+                lambda: os.remove(temp_file.name) 
+                if temp_file.name and os.path.exists(temp_file.name) 
+                else None
+            )
+            return temp_file.name 
+            
+        except Exception as e:
+            if isinstance(e, urllib.error.HTTPError) and e.getcode() == 404:
+                raise RuntimeError(f"Ligand {resname} not available from rcsb.org") from e
+            elif retry < max_retries - 1: 
+                wait_time = backoff_factor ** retry  
+                logger.info(f"Download failed for {resname}. Error: {e}. Retrying in {wait_time} seconds...")
+                time.sleep(wait_time)
+            else:
+                err = f"Max retries reached. Could not download CIF file for {resname}. Error: {e}"
+                raise RuntimeError(err) from e
+
+# Constants for deprotonate
+acidic_proton_loc_canonical = {
+        # any carboxylic acid, sulfuric/sulfonic acid/ester, phosphoric/phosphinic acid/ester
+        '[H][O]['+atom+'](=O)': 0 for atom in ('CX3', 'SX4', 'SX3', 'PX4', 'PX3')
+    } | {
+        # any thio carboxylic/sulfuric acid
+        '[H][O]['+atom+'](=S)': 0 for atom in ('CX3', 'SX4')
+    } | {
+        '[H][SX2][a]': 0, # thiophenol
+    }
+
+# Make free (noncovalent) CC
+def build_noncovalent_CC(basename: str) -> ChemicalComponent: 
+
+    with ChemicalComponent_LoggingControler(): 
+        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename), basename)
+        if cc_from_cif is None:
+            return None
+
+        cc = copy.deepcopy(cc_from_cif)
+        logger.info(f"*** using CCD ligand {basename} to construct residue {cc.resname} ***")
+
+        cc = cc.make_canonical(acidic_proton_loc = acidic_proton_loc_canonical)
+        if len(rdmolops.GetMolFrags(cc.rdkit_mol))>1:
+            err = f"Template Generation failed for {cc.resname}. Error: Molecule breaks into fragments during the deleterious editing. "
+            raise RuntimeError(err)
+
+        cc = cc.make_pretty_smiles()
+
+        # Check
+        try:
+            cc.ResidueTemplate_check()
+        except Exception as e:
+            err = f"Template {cc.resname} Failed to pass ResidueTemplate check. Error: {e}"
+            raise RuntimeError(err)
+            
+        logger.info(f"*** finish making {cc.resname} ***")
+    return cc
+
+
+def add_variants(cc_orig: ChemicalComponent, cc_list: list[ChemicalComponent] = [], 
+                 embed_allowed_smarts: str = None, 
+                 cap_allowed_smarts: str = None, cap_protonate: bool = False, 
+                 pattern_to_label_mapping_standard = dict[str, str], 
+                 variant_dict = dict[str, tuple]) -> list[ChemicalComponent]: 
+
+        for suffix in variant_dict:
+            cc = copy.deepcopy(cc_orig)
+            cc.resname += suffix
+            logger.info(f"*** using CCD residue {cc.parent} to construct {cc.resname} ***")
+
+            cc = (
+                cc
+                .make_canonical(acidic_proton_loc = acidic_proton_loc_canonical) 
+                .make_embedded(allowed_smarts = embed_allowed_smarts, 
+                            leaving_smarts_loc = variant_dict[suffix][0])
+                )
+            if len(rdmolops.GetMolFrags(cc.rdkit_mol))>1:
+                logger.warning(f"Molecule breaks into fragments during the deleterious editing of {cc.resname} -> skipping the vaiant... ")
+                continue
+
+            cc = (
+                cc
+                .make_capped(allowed_smarts = cap_allowed_smarts, 
+                            capping_smarts_loc = variant_dict[suffix][1],
+                            protonate = cap_protonate) 
+                .make_pretty_smiles()
+                .make_link_labels_from_patterns(pattern_to_label_mapping = pattern_to_label_mapping_standard)
+                )
+
+            try:
+                cc.ResidueTemplate_check()
+            except Exception as e:
+                err = f"Template {cc.resname} Failed to pass ResidueTemplate check. Error: {e}"
+                logger.error(err)
+                continue
+            
+            # Check redundancy
+            if any(cc == other_variant for other_variant in cc_list):
+                logger.error(f"Template Failed to pass redundancy check -> skipping the template... ")
+                continue
+
+            cc_list.append(cc)
+            logger.info(f"*** finish making {cc.resname} ***")
+
+        return cc_list
+
+
+class AA_recipe: 
+
+    embed_allowed_smarts = "[NX3]([H])([H])[CX4][CX3](=O)[O]"
+    cap_allowed_smarts = "[NX3][CX4][CX3](=O)"
+    cap_protonate = True
+    pattern_to_label_mapping_standard = {'[NX3h1]': 'N-term', '[CX3h1]': 'C-term'}
+
+    variant_dict = {
+            "":  ({"[NX3]([H])([H])[CX4][CX3](=O)[O]": {1, 6}}, None), # embedded amino acid
+            "_N": ({"[NX3]([H])([H])[CX4][CX3](=O)[O]": {6}}, {"[NX3][CX4][CX3](=O)": {0}}), # N-term amino acid
+            "_C": ({"[NX3]([H])([H])[CX4][CX3](=O)[O]": {1}}, None), # C-term amino acid
+        }
+    
+class NA_recipe: 
+
+    embed_allowed_smarts = "[O][PX4](=O)([O])[OX2][CX4][CX4]1[OX2][CX4][CX4][CX4]1[OX2][H]"
+    cap_allowed_smarts = "[OX2][CX4][CX4]1[OX2][CX4][CX4][CX4]1[OX2]"
+    cap_protonate = False
+    pattern_to_label_mapping_standard = {'[PX4h1]': '5-prime', '[O+0X2h1]': '3-prime'}
+    variant_dict = {
+            "_":  ({"[O][PX4](=O)([O])[OX2][CX4]": {0} ,"[CX4]1[OX2][CX4][CX4][CX4]1[OX2][H]": {6}}, None), # embedded nucleotide 
+            "_3": ({"[O][PX4](=O)([O])[OX2][CX4]": {0}}, None), # 3' end nucleotide 
+            "_5p": ({"[CX4]1[OX2][CX4][CX4][CX4]1[OX2][H]": {6}}, None), # 5' end nucleotide (extra phosphate than canonical X5)
+            "_5": ({"[O][PX4](=O)([O])[OX2][CX4]": {0,1,2,3}, "[CX4]1[OX2][CX4][CX4][CX4]1[OX2][H]": {6}}, {"[OX2][CX4][CX4]1[OX2][CX4][CX4][CX4]1[OX2]": {0}}), # 5' end nucleoside (canonical X5 in Amber)
+        }
+
+def build_linked_CCs(basename: str, embed_allowed_smarts: str = None, 
+                     cap_allowed_smarts: str = None, cap_protonate: bool = False, 
+                     pattern_to_label_mapping_standard = dict[str, str], 
+                     variant_dict = dict[str, tuple]) -> list[ChemicalComponent]: 
+
+    with ChemicalComponent_LoggingControler(): 
+        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename), basename)
+        if cc_from_cif is None:
+            return None
+        
+        # when no specific recipe is provided
+        if not embed_allowed_smarts:
+            # interpret possible embedding type
+            AA = cc_from_cif.rdkit_mol.GetSubstructMatch(Chem.MolFromSmarts(AA_recipe.embed_allowed_smarts))
+            NA = cc_from_cif.rdkit_mol.GetSubstructMatch(Chem.MolFromSmarts(NA_recipe.embed_allowed_smarts))
+
+            if not AA and not NA: 
+                logger.warning(f"Molecule doesn't contain the standard backbone of nucleotides or amino acids. -> no templates will be made. ")
+                return None
+            
+        else: 
+            # check if custom embedding is editable
+            editable = cc_from_cif.rdkit_mol.GetSubstructMatch(Chem.MolFromSmarts(embed_allowed_smarts))
+            if not editable:
+                logger.warning(f"Molecule doesn't contain embed_allowed_smarts: {embed_allowed_smarts} -> no templates will be made. ")
+                return None
+            
+            if not any(item is None for item in (cap_allowed_smarts, pattern_to_label_mapping_standard, variant_dict)):
+                return add_variants(cc_orig = cc_from_cif, cc_list = [], 
+                            embed_allowed_smarts = embed_allowed_smarts, 
+                            cap_allowed_smarts = cap_allowed_smarts, cap_protonate = cap_protonate, 
+                            pattern_to_label_mapping_standard = pattern_to_label_mapping_standard, 
+                            variant_dict = variant_dict)
+        
+        cc_variants = []
+        if AA: 
+            cc_variants = add_variants(cc_orig = cc_from_cif, cc_list = cc_variants, 
+                            embed_allowed_smarts = AA_recipe.embed_allowed_smarts, 
+                            cap_allowed_smarts = AA_recipe.cap_allowed_smarts, cap_protonate = AA_recipe.cap_protonate, 
+                            pattern_to_label_mapping_standard = AA_recipe.pattern_to_label_mapping_standard, 
+                            variant_dict = AA_recipe.variant_dict)
+        if NA:
+            cc_variants = add_variants(cc_orig = cc_from_cif, cc_list = cc_variants, 
+                embed_allowed_smarts = NA_recipe.embed_allowed_smarts, 
+                cap_allowed_smarts = NA_recipe.cap_allowed_smarts, cap_protonate = NA_recipe.cap_protonate, 
+                pattern_to_label_mapping_standard = NA_recipe.pattern_to_label_mapping_standard, 
+                variant_dict = NA_recipe.variant_dict)
+            
+    return cc_variants
+
+
+# XXX read from prepared? enumerate in stepwise? all protonation state variants, alter charge and update smiles/idx
diff --git a/meeko/linked_rdkit_chorizo.py b/meeko/linked_rdkit_chorizo.py
index f23a6294..172993e3 100644
--- a/meeko/linked_rdkit_chorizo.py
+++ b/meeko/linked_rdkit_chorizo.py
@@ -1,7 +1,9 @@
 import pathlib
 import json
 import logging
+import traceback
 from os import linesep as os_linesep
+from sys import exc_info
 from typing import Union
 from typing import Optional
 
@@ -20,7 +22,11 @@
 from .utils.rdkitutils import AtomField
 from .utils.rdkitutils import build_one_rdkit_mol_per_altloc
 from .utils.rdkitutils import _aux_altloc_mol_build
+from .utils.rdkitutils import covalent_radius
 from .utils.pdbutils import PDBAtomInfo
+from .chemtempgen import export_chem_templates_to_json
+from .chemtempgen import build_noncovalent_CC
+from .chemtempgen import build_linked_CCs
 
 import numpy as np
 
@@ -125,28 +131,7 @@ def find_inter_mols_bonds(mols_dict):
     -------
 
     """
-    covalent_radius = {  # from wikipedia
-        1: 0.31,
-        5: 0.84,
-        6: 0.76,
-        7: 0.71,
-        8: 0.66,
-        9: 0.57,
-        12: 0.00,  # hack to avoid bonds with metals
-        14: 1.11,
-        15: 1.07,
-        16: 1.05,
-        17: 1.02,
-        # 19: 2.03,
-        20: 0.00,
-        # 24: 1.39,
-        25: 0.00,  # hack to avoid bonds with metals
-        26: 0.00,
-        30: 0.00,  # hack to avoid bonds with metals
-        # 34: 1.20,
-        35: 1.20,
-        53: 1.39,
-    }
+
     allowance = 1.2  # vina uses 1.1 but covalent radii are shorter here
     max_possible_covalent_radius = (
         2 * allowance * max([r for k, r in covalent_radius.items()])
@@ -493,8 +478,36 @@ def _delete_residues(res_to_delete, raw_input_mols):
         raise ValueError(msg)
     return
 
+
 class ChorizoCreationError(RuntimeError):
-    pass
+
+    def __init__(self, error: str, recommendations: str = None): 
+        super().__init__(error) # main error message to pass to RuntimeError
+        self.error = error
+        self.recommendations = recommendations
+        exc_type, exc_value, exc_traceback = exc_info()
+        if exc_value is not None: 
+            self.traceback = ''.join(traceback.format_exception(exc_type, exc_value, exc_traceback))
+        else:
+            self.traceback = None
+
+    def __str__(self):
+        msg = "" + os_linesep
+        msg += "Error: Creation of data structure for receptor failed." + os_linesep
+        msg += "" + os_linesep
+        msg += "Details:" + os_linesep
+        msg += self.error + os_linesep
+
+        if self.traceback:
+            msg += self.traceback + os_linesep
+
+        if self.recommendations: 
+            msg += "Recommendations:" + os_linesep
+            msg += self.recommendations + os_linesep
+            msg += "" + os_linesep
+        
+        return msg
+
 
 def handle_parsing_situations(
     unmatched_res,
@@ -535,7 +548,13 @@ def handle_parsing_situations(
         err += msg
 
     if err:
-        raise ChorizoCreationError(err)
+        recs = "1. (for batch processing) Use -a/--allow_bad_res to automatically remove residues" + os_linesep
+        recs += "that do not match templates, and --default_altloc to set" + os_linesep
+        recs += "a default altloc variant. Use these at your own risk." + os_linesep
+        recs += "" + os_linesep
+        recs += "2. (processing individual structure) Inspecting and fixing the input structure is recommended." + os_linesep
+        recs += "Use --wanted_altloc to set variants for specific residues."
+        raise ChorizoCreationError(err, recs)
     return
 
 
@@ -706,8 +725,9 @@ def __init__(
         padders = residue_chem_templates.padders
         ambiguous = residue_chem_templates.ambiguous
 
-        # make sure all resiude_id in set_template exist
-        if set_template is not None:
+        if set_template is None:
+            set_template = {}
+        else:  # make sure all resiude_id in set_template exist
             missing = set(
                 [
                     residue_id
@@ -720,14 +740,92 @@ def __init__(
                     f"Residue IDs in set_template not found: {missing} {raw_input_mols.keys()}"
                 )
 
-        # currently allowing only one bond per residue pair
-        if any([len(v) > 1 for k, v in bonds.items()]):
-            msg = "got more than one bond between some residue pairs:"
-            for key, value in bonds.items():
-                if len(value) > 1:
-                    msg += f" {key}"
-            raise ValueError(msg)
-        bonds = {k: v[0] for k, v in bonds.items()}
+        # check if input assigned residue name in residue_templates
+        err = ""
+        supported_resnames = residue_templates.keys() | ambiguous.keys()
+        unknown_res_from_input = {res_id: raw_input_mols[res_id][1] 
+                                  for res_id in raw_input_mols 
+                                  if res_id not in set_template and raw_input_mols[res_id][1] not in supported_resnames
+                                  }
+        
+        if unknown_res_from_input:
+            unknown_valid_res_from_input = {k: v for k, v in unknown_res_from_input.items() if v != "UNL"}
+            if unknown_valid_res_from_input: 
+                err += f"Input residues {unknown_valid_res_from_input} not in residue_templates" + os_linesep
+            UNL_from_input = {k: v for k, v in unknown_res_from_input.items() if v == "UNL"}
+            if UNL_from_input: 
+                err += f"Input residues {UNL_from_input} do not have a concrete definition" + os_linesep
+        
+        unknown_res_from_assign = {}
+        if set_template:
+            unknown_res_from_assign = {res_id: resn for res_id, resn in set_template.items() if resn not in supported_resnames}
+            unknown_valid_res_from_assign = {k: v for k, v in unknown_res_from_assign.items() if v != "UNL"}
+            if unknown_valid_res_from_assign: 
+                err += f"Input residues {unknown_valid_res_from_assign} not in residue_templates" + os_linesep
+            UNL_from_assign = {k: v for k, v in unknown_res_from_assign.items() if v == "UNL"}
+            if UNL_from_assign: 
+                err += f"Input residues {UNL_from_assign} do not have a concrete definition" + os_linesep
+        
+        if err:
+            if "UNL" in err: 
+                err += "Resdiues that are named UNL can't be parameterized. " + os_linesep
+                rec = "1. (to parameterize the residues) Use --set_template to specify valid residue names, " + os_linesep
+                rec += "2. (to skip the residues) Use --delete_residues to ignore them. Residues will be deleted from the prepared receptor. "
+                raise ChorizoCreationError(err, rec)
+
+            print(err)
+            print("Trying to resolve unknown residues by building chemical templates... ")
+
+            all_unknown_res = unknown_res_from_input.copy()
+            all_unknown_res.update(unknown_res_from_assign)
+
+            bonded_unknown_res = {res_id: all_unknown_res[res_id] for res_id in all_unknown_res 
+                                  if any(res_id in respair for respair in bonds)}
+
+            unbound_unknown_res = all_unknown_res.copy()
+            for key in bonded_unknown_res:
+                unbound_unknown_res.pop(key, None) 
+
+            if unbound_unknown_res: 
+                for resname in set(unbound_unknown_res.values()): 
+                    try: 
+                        cc = build_noncovalent_CC(resname)
+                        fetch_template_dict = json.loads(export_chem_templates_to_json([cc]))['residue_templates'][cc.resname]
+                        residue_templates.update({resname: ResidueTemplate(
+                                                    smiles = fetch_template_dict['smiles'],
+                                                    atom_names = fetch_template_dict['atom_name'],
+                                                    link_labels = fetch_template_dict['link_labels'])})
+                        ambiguous[resname] = [cc.resname]
+                    except Exception as e: 
+                        print(f"Failed building template from CCD for {resname=}")
+                        raise ChorizoCreationError(str(e))
+
+            if bonded_unknown_res: 
+                failed_build = set()
+                try: 
+                    for resname in set(bonded_unknown_res.values()): 
+                        cc_list = build_linked_CCs(resname)
+                        if not cc_list: 
+                            failed_build.add(resname)
+                        else:
+                            for cc in cc_list:
+                                fetch_template_dict = json.loads(export_chem_templates_to_json([cc]))['residue_templates'][cc.resname]
+                                residue_templates.update({cc.resname: ResidueTemplate(
+                                                            smiles = fetch_template_dict['smiles'],
+                                                            atom_names = fetch_template_dict['atom_name'],
+                                                            link_labels = {int(key): value for key,value in fetch_template_dict['link_labels'].items()})})
+                                if resname in ambiguous: 
+                                    ambiguous[resname].append(cc.resname)
+                                else:
+                                    ambiguous[resname] = [cc.resname]
+                except Exception as e: 
+                    raise ChorizoCreationError(str(e))
+                            
+                if failed_build: 
+                    raise ChorizoCreationError(f"Template generation failed for unknown residues: {failed_build}, which appear to be linking fragments. " + os_linesep
+                                            + "Generation of chemical templates with modified backbones, which involves guessing of linker positions and types, are not currently supported. ", 
+                                            "1. (to parameterize the residues) Use --add_templates to pass the additional templates with valid linker_labels, " + os_linesep
+                                            + "2. (to skip the residues) Use --delete_residues to ignore them. Residues will be deleted from the prepared receptor. ")
 
         self.residues, self.log = self._get_residues(
             raw_input_mols,
@@ -739,14 +837,14 @@ def __init__(
         )
 
         _bonds = {}
-        for key, bond in bonds.items():
+        for key, bond_list in bonds.items():
             res1 = self.residues[key[0]]
             res2 = self.residues[key[1]]
             if res1.rdkit_mol is None or res2.rdkit_mol is None:
                 continue
             invmap1 = {j: i for i, j in res1.mapidx_to_raw.items()}
             invmap2 = {j: i for i, j in res2.mapidx_to_raw.items()}
-            _bonds[key] = (invmap1[bond[0]], invmap2[bond[1]])
+            _bonds[key] = [(invmap1[b[0]], invmap2[b[1]]) for b in bond_list]
         bonds = _bonds
 
         # padding may seem overkill but we had to run a reaction anyway for h_coord_from_dipep
@@ -1180,6 +1278,10 @@ def _get_residues(
             if set_template is not None and residue_key in set_template:
                 excess_H_ok = True  # e.g. allow set LYN (NH2) from LYS (NH3+)
                 template_key = set_template[residue_key]  # e.g. HID, NALA
+                if template_key not in residue_templates: 
+                    if template_key in ambiguous: 
+                        raise RuntimeError(f"Can't assign an ambiguous tamplate_key ({template_key}) to residue ({residue_key}). ")
+                    raise RuntimeError(f"Assigned tamplate_key ({template_key}) for residue ({residue_key}) is not in residue_templates. ")
                 template = residue_templates[template_key]
                 candidate_template_keys = [set_template[residue_key]]
                 candidate_templates = [template]
@@ -1206,11 +1308,12 @@ def _get_residues(
             if blunt_ends is None:
                 blunt_ends = []
             raw_atoms_with_bonds = []
-            for (r1, r2), (i, j) in bonds.items():
-                if r1 == residue_key:
-                    raw_atoms_with_bonds.append(i)
-                if r2 == residue_key:
-                    raw_atoms_with_bonds.append(j)
+            for (r1, r2), bond_list in bonds.items():
+                for i, j in bond_list:
+                    if r1 == residue_key:
+                        raw_atoms_with_bonds.append(i)
+                    if r2 == residue_key:
+                        raw_atoms_with_bonds.append(j)
 
             all_stats = {
                 "heavy_missing": [],
@@ -1372,10 +1475,16 @@ def _build_padded_mols(residues, bonds, padders):
                 atom.GetIdx(): atom.GetIdx() for atom in padded_mol.GetAtoms()
             }
             for atom_index, link_label in residue.link_labels.items():
-                adjacent_rid  = None
+                adjacent_rid = None
                 adjacent_mol = None
                 adjacent_atom_index = None
-                for (r1_id, r2_id), (i1, i2) in bonds.items():
+                for (r1_id, r2_id), bond_list in bonds.items():
+                    # TODO the second and subsequent bonds between a pair of
+                    # residues will not update the padding atoms with the
+                    # positions of the adjacent residues. This is OK, the same
+                    # happens for blunt residues, because the adjacent residue
+                    # is missing.
+                    i1, i2 = bond_list[0]
                     if r1_id == residue_id and i1 == atom_index:
                         adjacent_rid = r2_id
                         adjacent_atom_index = i2
@@ -2301,12 +2410,11 @@ def remove_atoms_with_mapping(product: Chem.Mol, mapping_numbers: set) -> Chem.M
     """Remove atoms with specific atom mapping numbers from a molecule."""
     editable_product = Chem.RWMol(product)
 
-    atoms_to_remove = []
-    for atom in editable_product.GetAtoms():
-        if atom.HasProp("molAtomMapNumber"):
-            map_num = int(atom.GetProp("molAtomMapNumber"))
-            if map_num in mapping_numbers:
-                atoms_to_remove.append(atom.GetIdx())
+    atoms_to_remove = [
+        atom.GetIdx() 
+        for atom in editable_product.GetAtoms() 
+        if atom.HasProp("molAtomMapNumber") and int(atom.GetProp("molAtomMapNumber")) in mapping_numbers
+    ]
     for idx in sorted(atoms_to_remove, reverse=True):
         editable_product.RemoveAtom(idx)
     
diff --git a/meeko/molsetup.py b/meeko/molsetup.py
index 2393c1ea..9d9216b4 100644
--- a/meeko/molsetup.py
+++ b/meeko/molsetup.py
@@ -1630,6 +1630,35 @@ def from_mol(
 
         return molsetup
 
+    @staticmethod
+    def remove_elements(mol, to_rm=(12, 20, 25, 26, 30)):
+        idx_to_rm = {}
+        neigh_idx_to_nr_h = {}
+        rm_to_neigh = {}
+        for atom in mol.GetAtoms():
+            if atom.GetAtomicNum() in to_rm:
+                idx_to_rm[atom.GetIdx()] = atom.GetFormalCharge()
+                rm_to_neigh[atom.GetIdx()] = set()
+                for neigh in atom.GetNeighbors():
+                    n = neigh.GetNumExplicitHs()
+                    neigh_idx_to_nr_h[neigh.GetIdx()] = n
+                    rm_to_neigh[atom.GetIdx()].add(neigh.GetIdx())
+        if not idx_to_rm:
+            return Chem.Mol(mol), idx_to_rm, rm_to_neigh
+        rwmol = Chem.EditableMol(mol)
+        for idx in sorted(idx_to_rm, reverse=True):
+            rwmol.RemoveAtom(idx)
+        mol = rwmol.GetMol()
+        for idx in neigh_idx_to_nr_h:
+            n = neigh_idx_to_nr_h[idx]
+            newidx = idx - sum([i < idx for i in idx_to_rm]) 
+            mol.GetAtomWithIdx(newidx).SetNumExplicitHs(n + 1)
+        mol.UpdatePropertyCache()
+        Chem.SanitizeMol(mol)
+        mol = Chem.AddHs(mol)
+        return mol, idx_to_rm, rm_to_neigh
+         
+
     def init_atom(self, assign_charges: bool, coords: list[np.ndarray]):
         """
         Generates information about the atoms in an RDKit Mol and adds them to an RDKitMoleculeSetup.
@@ -1647,12 +1676,46 @@ def init_atom(self, assign_charges: bool, coords: list[np.ndarray]):
         """
         # extract/generate charges
         if assign_charges:
-            copy_mol = Chem.Mol(self.mol)
+            things = self.remove_elements(self.mol)
+            copy_mol, idx_rm_to_formal_charge, rm_to_neigh = things
             for atom in copy_mol.GetAtoms():
                 if atom.GetAtomicNum() == 34:
                     atom.SetAtomicNum(16)
             rdPartialCharges.ComputeGasteigerCharges(copy_mol)
             charges = [a.GetDoubleProp("_GasteigerCharge") for a in copy_mol.GetAtoms()]
+            if idx_rm_to_formal_charge:
+                ok_charges = charges.copy()
+                for i in sorted(idx_rm_to_formal_charge, reverse=True):
+                    ok_charges.insert(i, 0.0)
+                nr_rm = len(idx_rm_to_formal_charge)
+                nr_added_h = copy_mol.GetNumAtoms() - self.mol.GetNumAtoms() + nr_rm
+                ok_charges = ok_charges[:len(ok_charges)-nr_added_h]
+                # print(f"{nr_added_h=}")
+                # print(f"{nr_rm=}")
+                # print(f"{idx_rm_to_formal_charge=}")
+                # print(f"{len(charges)=}")
+                # print(f"{len(ok_charges)=}")
+                # print(f"{copy_mol.GetNumAtoms()=}")
+                # print(f"{self.mol.GetNumAtoms()=}")
+                chrg_by_heavy_atom = {}
+                for i in range(nr_added_h):
+                    added_H_idx = self.mol.GetNumAtoms() + i - nr_rm
+                    # print(f"{added_H_idx=}")
+                    neighs = copy_mol.GetAtomWithIdx(added_H_idx).GetNeighbors()
+                    if len(neighs) != 1:
+                        raise RuntimeError("H should have 1 neighbor")
+                    if neighs[0].GetIdx() in chrg_by_heavy_atom:
+                        raise RuntimeError("expected only 1 added H per heavy atom, maybe deleted element had double bond to this heavy atom")
+                    chrg_by_heavy_atom[neighs[0].GetIdx()] = charges[added_H_idx]
+                # print(f"{chrg_by_heavy_atom=}")
+                for i, neighs in rm_to_neigh.items():
+                    # print(f"{i=}, {neighs=}")
+                    ok_charges[i] += idx_rm_to_formal_charge[i]
+                    for idx in neighs:
+                        newidx = idx - sum([i <= idx for i in idx_rm_to_formal_charge]) 
+                        # print(f"{idx=} {newidx=}")
+                        ok_charges[i] += chrg_by_heavy_atom[newidx]
+                charges = ok_charges
         else:
             charges = [0.0] * self.mol.GetNumAtoms()
         # register atom
diff --git a/meeko/utils/rdkitutils.py b/meeko/utils/rdkitutils.py
index a772c492..8c8851d1 100644
--- a/meeko/utils/rdkitutils.py
+++ b/meeko/utils/rdkitutils.py
@@ -249,3 +249,27 @@ def react_and_map(reactants: tuple[Chem.Mol], rxn: rdChemReactions.ChemicalReact
         outcomes.append((product, index_map))
 
     return outcomes
+
+
+covalent_radius = {  # from wikipedia
+    1: 0.31,
+    5: 0.84,
+    6: 0.76,
+    7: 0.71,
+    8: 0.66,
+    9: 0.57,
+    12: 0.00,  # hack to avoid bonds with metals
+    14: 1.11,
+    15: 1.07,
+    16: 1.05,
+    17: 1.02,
+    # 19: 2.03,
+    20: 0.00,
+    # 24: 1.39,
+    25: 0.00,  # hack to avoid bonds with metals
+    26: 0.00,
+    30: 0.00,  # hack to avoid bonds with metals
+    # 34: 1.20,
+    35: 1.20,
+    53: 1.39,
+}
diff --git a/scripts/mk_prepare_receptor.py b/scripts/mk_prepare_receptor.py
index d5fd23df..af9cde5e 100755
--- a/scripts/mk_prepare_receptor.py
+++ b/scripts/mk_prepare_receptor.py
@@ -520,20 +520,8 @@ def get_args():
                 wanted_altloc=wanted_altloc,
                 default_altloc=args.default_altloc,
             )
-        except ChorizoCreationError as err:
-            print()
-            print("Error: Creation of data structure for receptor failed.")
-            print()
-            print("Details: ")
-            print(err)
-            print()
-            print("Recommendations: ")
-            print("1. (for batch processing) Use -a/--allow_bad_res to automatically remove residues")
-            print("that do not match templates, and --default_altloc to set")
-            print("a default altloc variant. Use these at your own risk.")
-            print()
-            print("2. (processing individual structure) Inspecting and fixing the input structure is recommended.")
-            print("Use --wanted_altloc to set variants for specific residues.")
+        except ChorizoCreationError as e:
+            print(e)
             sys.exit(1)
 else:
     with open(args.read_pdb) as f:
@@ -550,17 +538,8 @@ def get_args():
             wanted_altloc=wanted_altloc,
             default_altloc=args.default_altloc,
         )
-    except ChorizoCreationError as err:
-        print()
-        print("Creation of data structure for receptor failed.")
-        print()
-        print(err)
-        print("Use -a/--allow_bad_res to automatically remove residues")
-        print("that do not match templates, and --default_altloc to set")
-        print("a default altloc variant. Use these at your own risk.")
-        print()
-        print("Inspecting and fixing the input structure is recommended.")
-        print("Use --wanted_altloc to set variants for specific residues.")
+    except ChorizoCreationError as e:
+        print(e)
         sys.exit(1)
 
 


Test patch:
diff --git a/test/linked_rdkit_chorizo_creation_test.py b/test/linked_rdkit_chorizo_creation_test.py
index 7987e01c..d8f4eae7 100644
--- a/test/linked_rdkit_chorizo_creation_test.py
+++ b/test/linked_rdkit_chorizo_creation_test.py
@@ -25,6 +25,7 @@
 insertion_code = pkgdir / "test/linked_rdkit_chorizo_data/1igy_B_82-83_has-icode.pdb"
 non_sequential_res = pkgdir / "test/linked_rdkit_chorizo_data/non-sequential-res.pdb"
 has_altloc = pkgdir / "test/linked_rdkit_chorizo_data/has-altloc.pdb"
+disulfide_adjacent = pkgdir / "test/linked_rdkit_chorizo_data/disulfide_bridge_in_adjacent_residues.pdb"
 
 
 # TODO: add checks for untested chorizo fields (e.g. input options not indicated here)
@@ -403,3 +404,15 @@ def test_altloc():
         if name == "OG":
             break
     assert abs(xyz[index][0] - 12.346) < 0.001
+
+def test_disulfide_adjacent():
+    """ disulfide bridge in adjacent residues broke a version of the code
+        that assumed only one bond between each pair of residues
+    """
+    with open(disulfide_adjacent, "r") as f:
+        pdb_text = f.read()
+    chorizo = LinkedRDKitChorizo.from_pdb_string(
+        pdb_text,
+        chem_templates,
+        mk_prep,
+    )
diff --git a/test/linked_rdkit_chorizo_data/disulfide_bridge_in_adjacent_residues.pdb b/test/linked_rdkit_chorizo_data/disulfide_bridge_in_adjacent_residues.pdb
new file mode 100644
index 00000000..6f600eaa
--- /dev/null
+++ b/test/linked_rdkit_chorizo_data/disulfide_bridge_in_adjacent_residues.pdb
@@ -0,0 +1,29 @@
+CRYST1  233.000  267.422   73.139  90.00  90.00  90.00 P 21 21 2    80
+ATOM      1  N   SER G 186      87.780 121.587  34.139  1.00 17.63      G    N  
+ATOM      2  CA  SER G 186      87.481 120.490  35.061  1.00 19.18      G    C  
+ATOM      3  C   SER G 186      86.239 120.784  35.901  1.00 19.64      G    C  
+ATOM      4  O   SER G 186      85.581 119.868  36.388  1.00 19.46      G    O  
+ATOM      5  CB  SER G 186      88.695 120.228  35.972  1.00 19.74      G    C  
+ATOM      6  OG  SER G 186      89.595 119.250  35.417  1.00 22.35      G    O  
+ATOM      7  N   CYS G 187      85.941 122.072  36.044  1.00 20.11      G    N  
+ATOM      8  CA  CYS G 187      84.819 122.610  36.791  1.00 20.95      G    C  
+ATOM      9  C   CYS G 187      83.465 122.202  36.250  1.00 21.24      G    C  
+ATOM     10  O   CYS G 187      82.504 122.084  37.011  1.00 21.20      G    O  
+ATOM     11  CB  CYS G 187      84.903 124.159  36.740  1.00 22.08      G    C  
+ATOM     12  SG  CYS G 187      84.272 124.967  35.209  1.00 25.14      G    S  
+ATOM     13  N   CYS G 188      83.378 122.046  34.927  1.00 21.63      G    N  
+ATOM     14  CA  CYS G 188      82.093 122.018  34.213  1.00 22.09      G    C  
+ATOM     15  C   CYS G 188      82.110 121.121  32.954  1.00 21.15      G    C  
+ATOM     16  O   CYS G 188      83.142 120.921  32.355  1.00 21.46      G    O  
+ATOM     17  CB  CYS G 188      81.690 123.456  33.849  1.00 22.89      G    C  
+ATOM     18  SG  CYS G 188      82.223 124.761  35.004  1.00 26.42      G    S  
+ATOM     19  N   PRO G 189      80.979 120.587  32.526  1.00 20.71      G    N  
+ATOM     20  CA  PRO G 189      81.000 119.599  31.427  1.00 20.34      G    C  
+ATOM     21  C   PRO G 189      81.070 120.183  29.988  1.00 19.58      G    C  
+ATOM     22  O   PRO G 189      81.602 119.530  29.095  1.00 18.94      G    O  
+ATOM     23  CB  PRO G 189      79.723 118.787  31.662  1.00 20.32      G    C  
+ATOM     24  CG  PRO G 189      78.748 119.767  32.390  1.00 20.43      G    C  
+ATOM     25  CD  PRO G 189      79.614 120.862  33.023  1.00 21.08      G    C  
+TER   
+HETATM   26  O   HOH G2056      88.268 122.391  38.090  1.00 63.88      UA   O  
+END