class ReactionTemplate:
"""
Data structure for reactions with R-groups
"""
def __init__(self, MG=[]):
"""
path identifies a folder containing the component .xyz and .txt files
MG - instance of MatchingGroupsData
"""
self.reactants = Molecule()
self.products = Molecule()
self.name = ""
self.conditions = Conditions()
# self.MG = MG # matching
self.rules = Rules()
self.major_products = Molecule()
self.major_reactants = Molecule()
self.minor_products = Molecule()
self.minor_reactants = Molecule()
# the spectators is a special beast.
# it is a molecule that overall is here to balance charges, but it has special properties.
# 1) it DOES NOT MATCH in the forward direction. It is a molecule that is a minor reactant, yes,
# but that is not strictly needed from the reaction to happen. Say you are adding HCl to a molecule
# to make the molecule react. The H+ is a minor reactant (it is needed by the reaction, even as catalyzer
# but when making the template, Cl- can be replaced by whatever counterion. It is a spectator.
# 2) When going in the back direction, it is a molecule that should be added in the graph as a minor product
# decided by the template maker (and it can match, or not, it does not matter in this direction).
self.spectators = Molecule()
self.changes_react_to_prod = None # Changes
self.changes_prod_to_react = None # Changes
self.side_reaction = False # To avoid making pots with a side_reaction template
self.react_with_itself = False # This is the logic for creating pots with double major reactants.
self.folder = ''
@property
def major_reactants_list(self):
return self.major_reactants.separate_graph_in_pieces()
def __str__(self):
string = self.name
return string
def _repr_html_(self):
return "<p>" + str(self).replace("\n", "<br>") + "</p>"
def rxn_smiles(self):
reactant_smiles = self.reactants.force_smiles()
product_smiles = self.products.force_smiles()
return '->'.join([reactant_smiles, product_smiles]).replace('[', '').replace(']', '')
@classmethod
def from_components(cls, name, reactants, changes_react_to_prod_dict, conditions, rules, MG=None, collapse_groups=None, spectators=[], side_reaction=False, react_with_itself=False, folder=''):
"""
Construct a new ReactionTemplate from individual components using a new edge-update-method
name - string
reactants - Molecule object, possibly with more than one connected component and with R groups
changes_react_to_prod - Dictionary of changes to update edges and create a Reactant
conditions - conditions Object
product_smiles - a smiles string for a non-template example of the reaction product
FG - MatchingGroupsData instance
spectators - this is a list of smiles
collapse_groups example:
[(6, 4, 'N-Me'), (22,23,'Me'), (28,29,'Me')]
"""
new_rxn = cls()
new_rxn.name = name
new_rxn.conditions = conditions
new_rxn.side_reaction = side_reaction
new_rxn.react_with_itself = react_with_itself
new_rxn.folder = folder
changes_react_to_prod = Changes.from_dict(changes_react_to_prod_dict)
products = reactants.update_edges_and_charges(changes_react_to_prod)
new_rxn.products = products
new_rxn.reactants = reactants
assert changes_react_to_prod.validate_changes(), "You need to doublecheck the changes dictionary!! Do not put duplicates!!"
unique_changes_react_to_prod = changes_react_to_prod
unique_changes_prod_to_react = changes_react_to_prod.reverse_changes(reactants)
if collapse_groups:
print(f"I am applying the collapse code using {collapse_groups}")
new_rxn.changes_prod_to_react = apply_permute(products.collapse_nodes(collapse_groups), unique_changes_prod_to_react)
new_rxn.changes_react_to_prod = apply_permute(reactants.collapse_nodes(collapse_groups), unique_changes_react_to_prod)
else:
new_rxn.changes_prod_to_react = apply_permute(products, unique_changes_prod_to_react)
new_rxn.changes_react_to_prod = apply_permute(reactants, unique_changes_react_to_prod)
new_rxn.MG = MG
new_rxn.rules = rules
spectators_graph = Molecule.from_list_of_smiles_to_single_graph(spectators)
major_products, minor_products = from_graph_to_pieces(products, spectators_graph=spectators_graph)
major_reactants, minor_reactants = from_graph_to_pieces(reactants, spectators_graph=spectators_graph)
# Above there are all the chemically relevant permutations of changes, but all of them are applied to the
# first set of "minor" and "major" reactant/product (the one that the template creator made).
# So the following attributes are still unique.
new_rxn.major_reactants = major_reactants
new_rxn.minor_reactants = minor_reactants
new_rxn.major_products = major_products
new_rxn.minor_products = minor_products
new_rxn.spectators = spectators_graph
return new_rxn
@staticmethod
def load(folder="."):
"""
loads everything from a folder
"""
reaction_name = os.path.basename(folder)
pickle_path = os.path.join(folder, f"{reaction_name}.p")
return pload(pickle_path)
def create_map_for_consecutive_rs(self):
i = 1
mapping = {}
for x in self.reactants.nodes:
if self.reactants.nodes[x]["element"][0] == "R":
mapping[self.reactants.nodes[x]["element"]] = f"R{i}"
i += 1
return mapping
def create_map_for_consecutive_rs_reverse(self):
i = 1
mapping = {}
for x in self.reactants.nodes:
if self.reactants.nodes[x]["element"][0] == "R":
mapping[f"R{i}"] = self.reactants.nodes[x]["element"]
i += 1
return mapping
def validate(self):
"""
this is a small helper function that makes sure that silly things are ok
"""
if not self.reactants.is_template:
print(f"\n\nWarning {self.name}!!!! This Template does not contain R groups!!! You forgetting something?\n\n")
assert self.products != self.reactants, f"Please control reaction {self.name}. Products and reactants are the same molecules."
assert self.names_in_template_are_good(), f"Please control {self.name}. It has strange names (ewg or base) in it."
assert self.is_there_rules_ismags_collision(), f"Ismags collision detected in rxn {self.name}. Two R groups have same symmetry but different rules."
def is_there_rules_ismags_collision(self):
rea = self.reactants
rul = self.rules
return rules_ismags_collision(rea, rul)
def generate_single_template_and_rules(self):
"""
This will return a list of TemplateAndRules from the template.
It is a list of major reactants with associated rules.
"""
list_of_tar = []
for piece in self.major_reactants.separate_graph_in_pieces():
new_rules = Rules.filter_entries_from_names(self.rules, piece.which_r_groups)
element = TemplateAndRules(template=piece, rules=new_rules)
list_of_tar.append(element)
return list_of_tar
def names_in_template_are_good(self):
"""
This functions makes a test to avoid strange names in templates.
Only R or elements are allowed.
"""
graph = self.reactants
al = [
"h",
"he",
"li",
"be",
"b",
"c",
"n",
"o",
"f",
"ne",
"na",
"mg",
"al",
"si",
"p",
"s",
"cl",
"ar",
"k",
"ca",
"sc",
"ti",
"v",
"cr",
"mn",
"fe",
"co",
"ni",
"cu",
"zn",
"ga",
"ge",
"as",
"se",
"br",
"kr",
"rb",
"sr",
"y",
"zr",
"nb",
"mo",
"tc",
"ru",
"rh",
"pd",
"ag",
"cd",
"in",
"sn",
"sb",
"te",
"i",
"xe",
"cs",
"ba",
"la",
"ce",
"pr",
"nd",
"pm",
"sm",
"eu",
"gd",
"tb",
"dy",
"ho",
"er",
"tm",
"yb",
"lu",
"hf",
"ta",
"w",
"re",
"os",
"ir",
"pt",
"au",
"hg",
"tl",
"pb",
"bi",
"po",
"at",
"rn",
"fr",
"ra",
"ac",
"th",
"pa",
"u",
"np",
"pu",
]
names = [graph.nodes[x]["element"] for x in graph.nodes()]
filtered = [x for x in names if x.lower() not in al]
filter_r = [x for x in filtered if x.lower()[0] != "r"]
return len(filter_r) == 0
def default_test(self, verbosity=0):
"""
This method just quickly test the apply method on the template, in order to save it.
"""
self.validate()
productsQ = self.products.qify(number=False)
reactantsQ = self.reactants.qify(number=False)
test_prods = self.apply_test(self.products, productsQ, self.changes_prod_to_react)
test_reacts = self.apply_test(self.reactants, reactantsQ, self.changes_react_to_prod)
bool1 = all([productsQ == x for x in test_reacts])
bool2 = all([reactantsQ == x for x in test_prods])
# Test charge is conserved
charges_react = self.reactants.charge
charges_prods = self.products.charge
bool3 = charges_prods == charges_react
final_test = bool1 and bool2 and bool3
return final_test
def rename_rs_from_map(self):
"""
Given a reaction template with non consecutive Rs, this will rename and remap Rs to be consecutive.
This is used before saving.
"""
mapping = self.create_map_for_consecutive_rs_reverse()
self.rules = self.rules.remap_rules_reverse(mapping)
self.reactants = self.reactants.remap_element_names_reverse(mapping)
self.products = self.products.remap_element_names_reverse(mapping)
self.major_products = self.major_products.remap_element_names_reverse(mapping)
self.major_reactants = self.major_reactants.remap_element_names_reverse(mapping)
return self
def save(self, folder=".", verbosity=0):
"""
Saves everything into folder
"""
self.folder = sys.argv[0]
self.rename_rs_from_map()
if self.default_test(verbosity=verbosity):
print(f"{self.name} - Basic test PASSED.")
file_name_no_spaces = self.name.replace(" ", "-")
full_path = os.path.join(folder, "reaction_templates", file_name_no_spaces)
if os.path.exists(full_path):
if verbosity > 0:
print(f"Folder {full_path} exists.")
else:
if verbosity > 0:
print(f"Creating folder {full_path}.")
os.makedirs(full_path, exist_ok=True)
psave(self, os.path.join(full_path, f"{file_name_no_spaces}.p"))
else:
print(f"Basic test failed. Template {self.name} has been NOT written to disc.")
def add_default_rules(self):
if self.rules is None:
rules = Rules()
print("self has no rules to enforce, creating default H and C rules")
groups_default = [("H", "Hydrogen", Molecule.from_smiles("[HH]").substitute_group(1, "L")), ("C", "Carbon", Molecule.from_smiles("[HH]").substitute_groups([(1, "L"), (0, "C")]))]
collection = {}
defaults = MatchingGroupsData.from_tuple_list(groups_default, collection)
# who are the Rs?
r_labels = [value[1]["element"] for i, value in enumerate(self.products.nodes(data=True)) if value[1]["element"][0] == "R"]
auto_enforce = {} # enforce rule dictionary
for single_r in r_labels:
auto_enforce[single_r] = ["H", "C"]
# I create the condition from the labels dictionary
r_condition = SingleCondition.from_dict(auto_enforce, defaults)
self.rules = rules.append_enforce(r_condition)
def changes_report(self):
print(f"template has {len(self.changes_react_to_prod)} changes forward and {len(self.changes_prod_to_react)} changes backwards")
print(self.changes_react_to_prod)
print(self.changes_prod_to_react)
def minor_reactant_break_list(self, mol_graph_list):
"""
from a list of molecule graph, we use filter_minor_reactants
this returns a list of list of molecules to break apart.
"""
return [self.filter_minor_reactants(x).separate_graph_in_pieces() for x in mol_graph_list]
def filter_minor_reactants(self, mol_graph):
"""
from a molecule graph, we filter out the minor reactant of the template.
this returns a molecule. This contains the spectator so it cannot be used on forward application anymore
"""
diff = mol_graph.graph_difference(self.minor_reactants + self.spectators)
return diff
def filter_minor_product(self, mol_graph):
"""
from a molecule graph, we filter out the minor reactant of the template.
this returns a molecule. This contains the spectator so it cannot be used on forward application anymore
"""
diff = mol_graph.graph_difference(self.minor_products + self.spectators)
return diff
def quick_apply(self, smiles, verbosity=0):
"""
this is a function to be used to quickly check application of a template with a molecule
"""
mol = Molecule.from_smiles(smiles)
fr = apply(self.reactants, mol + self.minor_reactants, self.changes_react_to_prod, self.rules, name="quick_apply", verbosity=verbosity)
final_results_unique = Molecule.delete_duplicates(fr)
return final_results_unique
def quick_apply_debug(self, mols, verbosity=0):
"""
a different quick_apply. This does not clear duplicates and can take molecules as input (this should really be fused with quick apply above one day)
"""
if isinstance(mols, str):
mol = Molecule.from_smiles(mols)
elif isinstance(mols, Molecule):
mol = mols
fr = apply(self.reactants, mol + self.minor_reactants, self.changes_react_to_prod, self.rules, "quick_apply", verbosity=verbosity)
return fr
def apply_forward(self, target_mol, filter_minor_products=True, skip_rules=False, verbosity=0):
"""
this is the one step forward for a single template to be used in the pot
"""
final_results = apply(self.reactants, target_mol, self.changes_react_to_prod, self.rules, name=self.name, skip_rules=skip_rules, verbosity=verbosity)
if filter_minor_products:
final_results = [self.filter_minor_product(x) for x in final_results]
final_results_unique = Molecule.delete_duplicates(final_results, verbosity=verbosity)
return final_results_unique
def apply_backward(self, target_mol, verbosity=0):
"""
this is the one step backward for a single template to be used to create a single pot
"""
# Backward should add minor things
target_molecule = target_mol + self.minor_products + self.spectators
final_results = apply(self.products, target_molecule, self.changes_prod_to_react, self.rules, name=self.name, verbosity=verbosity)
final_results_unique = Molecule.delete_duplicates(final_results)
return final_results_unique
def apply_backwards_convergence(self, mol, filtering=False):
"""
This will return every possible application of the same template backwards to one molecule graph.
This will exhaust every possible combination of them.
It returns a list of molecules.
"""
first_layer = self.apply_backward(mol)
if filtering:
filtered = [self.filter_minor_reactants(x) for x in first_layer]
else:
filtered = first_layer
final = self.permute_backwards(filtered, filtering=filtering)
unique = Molecule.delete_duplicates(final)
return unique
def permute_backwards(self, pots, filtering=False):
"""
This is the recursive part of the function, that goes through every combination.
"""
if pots == []:
return []
else:
d = []
for c in pots:
results = self.apply_backward(c)
if filtering:
filtered = [self.filter_minor_reactants(x) for x in results]
else:
filtered = results
d = d + filtered # I need to concatenate. Those are two lists.
# TODO test yield
return pots + self.permute_backwards(d, filtering=filtering)
def apply_test(self, templ_graph, target_mol, changes, verbosity=0):
"""
this is the test application logic
"""
final_results = apply(templ_graph, target_mol, changes, Rules(), name=self.name, verbosity=verbosity)
return final_results
def how_many_reactant_templates(rxn):
"""
this function will return the number of R molecules present in the reactants of this template
"""
number_of_components = nx.number_connected_components(rxn.major_reactants)
return number_of_components
def create_mapping_from_purchasable(self, new_moldict):
outer_list = []
print(f'Creating mapping for with {self.name}')
for single in self.generate_single_template_and_rules():
this_index_list = []
for smi, purch_index in tqdm(new_moldict.items()):
mol = Molecule.from_smiles(smi)
if single.is_subgraph_isomorphic_with_rules(mol):
this_index_list.append(purch_index)
outer_list.append(this_index_list)
recap = '\n'.join([f'{len(x)} molecules isomorphic to reactant n. {i}' for i, x in enumerate(outer_list)])
print(f'I found:\n{recap}')
n_pots = math.prod([len(x) for x in outer_list])
print(f'I can generate {n_pots} pots from this template, given the purchasables\n\n')
if n_pots == 0:
return None
else:
return outer_list
# DRAWING
def draw(self, size=(800, 800), table_width=60, string_mode=False, percentage=1, node_index=False, charges=False, mode="rdkit", draw_rules=True, fixed_bond_length=40, doi_titles=[], depth=2):
"""
Draw the reaction
"""
# I cannot draw in rdkit Nu Lg etc etc. Quick fix.
if any([self.reactants.nodes[x]["element"] in ["Lg", "Nu", "EWG", "EWG1", "EWG2", "EWG3", "CG"] for x in self.reactants.nodes]):
mode = "d3"
draw_rules = False
if doi_titles:
assert len(doi_titles) == len(self.conditions.doi)
unpacked_dois = pretty_print_dois(doi_titles, self.conditions.doi)
else:
unpacked_dois = ",".join([f'<a href="{x}">{a}</a>' for a, x in enumerate(self.conditions.doi)])
# reactantZ = self.products.update_edges_and_charges(self.changes_prod_to_react[0])
rwi_string = " REACT WITH ITSELF" if self.react_with_itself else ""
sr_string = " SIDE REACTION" if self.side_reaction else ""
specials_string = " ".join(filter(None, [rwi_string, sr_string]))
final_special_string = ""
if specials_string:
final_special_string = "<br>Special Conditions :" + specials_string + "<br>"
totalstring = f"""<h{depth}>{self.name}</h{depth}><p>References: {unpacked_dois} </p>
<div> {final_special_string} </div>
<div> {self.rules._repr_html_() if draw_rules else ""}</div>
<div style="width: {table_width}%; display: table;"> <div style="display: table-row;">
<div style="width: 40%; display: table-cell;">
<p style="text-align: center; font-weight: bold;">Reactants</p>
{self.reactants.draw(string_mode=True, percentage=percentage, size=size, node_index=node_index, charges=charges, mode=mode, fixed_bond_length=fixed_bond_length)}
</div>
<div style="width: 20%; display: table-cell; vertical-align: middle;">
{self.conditions.draw_reaction_arrow()}
</div>
<div style="width: 40%; display: table-cell;">
<p style="text-align: center; font-weight: bold;">Products</p>
{self.products.draw(string_mode=True, percentage=percentage, size=size, node_index=node_index, charges=charges, mode=mode, fixed_bond_length=fixed_bond_length)}
</div>
</div></div>
"""
if string_mode:
return totalstring
else:
return HTML(totalstring)
def get_header_html(self,
draw_rules=True,
doi_titles=[],
depth=2):
""" """
# I cannot draw in rdkit Nu Lg etc etc. Quick fix.
if any([self.reactants.nodes[x]["element"] in ["Lg", "Nu", "EWG", "EWG1", "EWG2", "EWG3", "CG"] for x in self.reactants.nodes]):
draw_rules = False
if doi_titles:
assert len(doi_titles) == len(self.conditions.doi)
unpacked_dois = pretty_print_dois(doi_titles, self.conditions.doi)
else:
unpacked_dois = ",".join([f'<a href="{x}">{a}</a>' for a, x in enumerate(self.conditions.doi)])
# reactantZ = self.products.update_edges_and_charges(self.changes_prod_to_react[0])
rwi_string = " REACT WITH ITSELF" if self.react_with_itself else ""
sr_string = " SIDE REACTION" if self.side_reaction else ""
specials_string = " ".join(filter(None, [rwi_string, sr_string]))
final_special_string = ""
if specials_string:
final_special_string = "Special Conditions :" + specials_string
totalstring = f"""<h{depth}>{self.name}</h{depth}><p>References: {unpacked_dois} </p>
<div> {final_special_string} </div>
<div> {self.rules._repr_html_() if draw_rules else ""}</div>
"""
return totalstring
def get_reaction_html(self,
size=(400, 400),
percentage=1,
node_index=False,
charges=False,
mode="rdkit",
draw_rules=True,
fixed_bond_length=40,
doi_titles=[],
depth=2,
):
""" """
total_string = f"""<div style="width: 60%; display: table;"> <div style="display: table-row;">
<div style="width: 40%; display: table-cell;">
<p style="text-align: center; font-weight: bold;">Reactants</p>
{self.reactants.draw(string_mode=True, percentage=percentage, size=size, node_index=node_index, charges=charges, mode=mode, fixed_bond_length=fixed_bond_length)}
</div>
<div style="width: 20%; display: table-cell; vertical-align: middle;">
{self.conditions.draw_reaction_arrow()}
</div>
<div style="width: 40%; display: table-cell;">
<p style="text-align: center; font-weight: bold;">Products</p>
{self.products.draw(string_mode=True, percentage=percentage, size=size, node_index=node_index, charges=charges, mode=mode, fixed_bond_length=fixed_bond_length)}
"""
return total_string
# OTHERS
def remove_empty(self, dict_list_mol):
"""
Remove empty molecule graphs from a dictionary of lists of molecules
"""
empty_mols = {k: [nx.is_empty(m) for m in v] for k, v in dict_list_mol.items()}
for k, v in dict_list_mol.items():
for i, mol in enumerate(v):
if empty_mols[k][i]:
del dict_list_mol[k][i]
return dict_list_mol
def simple_html(self, **kwargs):
list_of_mols = [self.reactants, self.products]
return Molecule.draw_list(list_of_mols, **kwargs)
def to_png(self, folder=None, file_name=None, arrow_size=12, fixed_bond_length=None):
if file_name is None:
file_name = f'{self.name}.png'
kitoptions = {
"transparent": "",
}
imgkit.from_string(self.simple_html(arrows=True,
columns=6,
size=(250, 250),
string_mode=True,
arrow_size=arrow_size,
fixed_bond_length=fixed_bond_length,
display_title=False),
f'{self.name}.png',
kitoptions)