feat: account for radical cofactors

docs/source/rec_cli_options.rst

@@ -107,6 +107,10 @@ Receptor Perception Options
    Default is gasteiger, 'zero' sets all zeros, 'read' requires --read_pqr
 
 
+.. option:: --radical_cofactor
+   Meeko fails if the receptor contains a radical cofactor. This option relaxes some of the 
+   checks so that radical cofactors are handled properly. 
+
 Grid Box Options
 ~~~~~~~~~~~~~~~~
 

meeko/chemtempgen.py

@@ -372,22 +372,27 @@ def recharge(rwmol: Chem.RWMol) -> Chem.RWMol:
 
 
 # Attribute Formatters
-def get_smiles_with_atom_names(mol: Chem.Mol) -> tuple[str, list[str]]:
+def get_smiles_with_atom_names(mol: Chem.Mol, radical_cofactor) -> tuple[str, list[str]]:
     """
     Generate SMILES with atom names in the order of SMILES output
 
     Parameters
     ----------
     mol : Chem.Mol
         Input molecule
+    radical_cofactor:
+        true of the input is meant to be a radical
 
     Returns
     -------
     tuple[str, list[str]]
         Tuple containing the SMILES string and a list of 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)
+    if radical_cofactor:
+        smiles_exh = Chem.MolToSmiles(mol)
+    else:
+        smiles_exh = Chem.MolToSmiles(mol, allHsExplicit=True)
 
     smiles_atom_output_order = mol.GetProp('_smilesAtomOutputOrder')
     delimiters = ('[', ']', ',')
@@ -515,7 +520,7 @@ def __eq__(self, other):
         return False
 
     @classmethod
-    def from_cif(cls, source_cif: str, resname: str):
+    def from_cif(cls, source_cif: str, resname: str, radical_cofactor):
         """
         Create ChemicalComponent from a chemical component dict file and a searchable residue name in file.
 
@@ -553,6 +558,8 @@ def from_cif(cls, source_cif: str, resname: str):
             if len(element)==2:
                 element = element[0] + element[1].lower()
             rdkit_atom = Chem.Atom(element)
+            if radical_cofactor:
+                rdkit_atom.SetNoImplicit(True)
             for attr in atom_attributes:
                 rdkit_atom.SetProp(attr, atom_cols[attr][idx])
                 # strip double quotes in names
@@ -607,6 +614,8 @@ def from_cif(cls, source_cif: str, resname: str):
         # Finish eidting mol 
         try:    
             rwmol.UpdatePropertyCache()
+            if radical_cofactor:
+                Chem.SanitizeMol(rwmol)
         except Exception as e:
             logger.error(f"Failed to create rdkitmol from cif. Error: {e} -> template for {resname} will be None. ")
             return None
@@ -620,7 +629,7 @@ def from_cif(cls, source_cif: str, resname: str):
         rdkit_mol = rwmol.GetMol()
 
         # Get Smiles with explicit Hs and ordered atom names
-        smiles_exh, atom_name = get_smiles_with_atom_names(rdkit_mol)
+        smiles_exh, atom_name = get_smiles_with_atom_names(rdkit_mol, radical_cofactor)
 
         return cls(rdkit_mol, resname, smiles_exh, atom_name)
 
@@ -693,7 +702,7 @@ def make_capped(self, allowed_smarts, capping_names = None, capping_smarts_loc =
 
     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, self.atom_name = get_smiles_with_atom_names(self.rdkit_mol, False)
         self.smiles_exh = get_pretty_smiles(self.smiles_exh)
         return self
 
@@ -912,7 +921,9 @@ def fetch_from_pdb(resname: str, max_retries = 5, backoff_factor = 2, ignore_htt
     }
 
 # Standard pipelines
-def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalComponent: 
+def build_noncovalent_CC(basename: str, 
+                         ignore_https_cert = False, 
+                         radical_cofactor=False) -> ChemicalComponent: 
     """
     Build a noncovalent chemical component from a CIF file.
 
@@ -927,7 +938,7 @@ def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalCo
         The constructed ChemicalComponent instance.
     """
     with ChemicalComponent_LoggingControler(): 
-        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename)
+        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename, radical_cofactor)
         if cc_from_cif is None:
             return None
 
@@ -939,7 +950,8 @@ def build_noncovalent_CC(basename: str, ignore_https_cert = False) -> ChemicalCo
             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()
+        if not radical_cofactor:
+            cc = cc.make_pretty_smiles()
 
         # Check
         try:
@@ -1057,7 +1069,8 @@ class NA_recipe:
 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], ignore_https_cert = False) -> list[ChemicalComponent]: 
+                     variant_dict = dict[str, tuple], ignore_https_cert = False, 
+                     radical_cofactor=False) -> list[ChemicalComponent]: 
     """
     Build a linked chemical component from a CIF file.
 
@@ -1084,7 +1097,7 @@ def build_linked_CCs(basename: str, embed_allowed_smarts: str = None,
         List of ChemicalComponent instances with the added variants.
     """
     with ChemicalComponent_LoggingControler(): 
-        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename)
+        cc_from_cif = ChemicalComponent.from_cif(fetch_from_pdb(basename, ignore_https_cert=ignore_https_cert), basename, radical_cofactor)
         if cc_from_cif is None:
             return None
 

meeko/cli/mk_prepare_receptor.py

@@ -255,7 +255,6 @@ def get_args():
         help='specify the residues for which to make terminal functional group rotatable by the chain ID and residue number, e.g. -t ":42,B:23" is equivalent to -t ":42" -t "B:23" (leave chain ID empty if omitted in input PDB or mmCIF)',
     )
 
-
     config_group.add_argument(
         "--compute_charges",
         help="compute charges from scratch with the given charge model instead of reading from template (note: this option is slower)",
@@ -268,6 +267,11 @@ def get_args():
         help="default is gasteiger, 'zero' sets all zeros, 'read' requires --read_pqr",
         default=None,
     )
+    config_group.add_argument(
+        "--radical_cofactor",
+        help="if the receptor has a radical cofactor and meeko faisl, try this option",
+        action="store_true",
+    )
 
     box_group = parser.add_argument_group("Size and center of grid box")
     box_group.add_argument(
@@ -590,6 +594,8 @@ def main():
 
     mk_config["compute_charges"] = args.compute_charges
 
+    mk_config["radical_cofactor"] = args.radical_cofactor
+
     # update config by inputs from arguments
 
     if args.charge_model is not None: 

meeko/polymer.py

@@ -1243,7 +1243,9 @@ def __init__(
             if unbound_unknown_res: 
                 for resname in set(unbound_unknown_res.values()): 
                     try: 
-                        cc = build_noncovalent_CC(resname, ignore_https_cert=ignore_https_cert)
+                        cc = build_noncovalent_CC(resname, 
+                                                  ignore_https_cert=ignore_https_cert, 
+                                                  radical_cofactor=mk_prep.radical_cofactor)
                         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'],

meeko/preparation.py

@@ -115,6 +115,7 @@ def __init__(
         crippen_as_solpar=False,
         override_ad4sol_par_including_q=False,
         override_ad4sol_par_including_q_qasp=0.0,
+        radical_cofactor=False,
     ):
         """
 
@@ -246,6 +247,8 @@ def __init__(
         self.add_index_map = add_index_map
         self.remove_smiles = remove_smiles
 
+        self.radical_cofactor=radical_cofactor
+
         self._bond_typer = BondTyperLegacy()
         self._macrocycle_typer = FlexMacrocycle(
             self.min_ring_size,

test/chemtempgen_test.py

@@ -55,6 +55,16 @@ def test_build_noncovalent_CC():
 
     assert template_equality_check(default_template_file, basename, "_fl-ccd", cc)
 
+
+def test_radical_cofactor_NO():
+    basename = "NO"  # free ligand from CCD
+    cc = build_noncovalent_CC(basename, radical_cofactor=True)
+
+    assert cc is not None
+    assert isinstance(cc, ChemicalComponent)
+
+    assert cc.smiles_exh == '[N]=O'
+
 def test_add_variants():
     basename = "AMP" 
     cc_list = build_linked_CCs(basename)