From 8fe531f515dca9c715e3fcfb44e60339ab101d2d Mon Sep 17 00:00:00 2001 From: kilian Date: Thu, 13 Aug 2026 14:29:55 +0200 Subject: [PATCH] Dateien nach "mod/butadien" hochladen --- mod/butadien/butadienCycleExclusion.py | 217 +++++++++++++++++++++++++ mod/butadien/graph.py | 136 ++++++++++++++++ 2 files changed, 353 insertions(+) create mode 100644 mod/butadien/butadienCycleExclusion.py create mode 100644 mod/butadien/graph.py diff --git a/mod/butadien/butadienCycleExclusion.py b/mod/butadien/butadienCycleExclusion.py new file mode 100644 index 0000000..bb9e580 --- /dev/null +++ b/mod/butadien/butadienCycleExclusion.py @@ -0,0 +1,217 @@ +#config.ilp.solver="CPLEX" +#import mod +#from mod import * +import networkx as nx + +# import Graph class from graph.py +from graph import GraphObj + +butadien = Graph.fromGMLString( +"""graph [ + node [ id 0 label "C" ] + node [ id 1 label "C" ] + node [ id 2 label "C" ] + node [ id 3 label "C" ] + node [ id 4 label "H" ] + node [ id 5 label "H" ] + node [ id 6 label "H" ] + node [ id 7 label "H" ] + node [ id 8 label "H" ] + node [ id 9 label "H" ] + edge [ source 0 target 1 label "=" ] + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 0 target 4 label "-" ] + edge [ source 0 target 5 label "-" ] + edge [ source 1 target 6 label "-" ] + edge [ source 2 target 7 label "-" ] + edge [ source 3 target 8 label "-" ] + edge [ source 3 target 9 label "-" ] +]""" +, name="Butadien") + + +#pentadien = Graph.fromGMLString( +"""graph +[ + node [ id 0 label "C" ] + node [ id 1 label "C" ] + node [ id 2 label "C" ] + node [ id 3 label "C" ] + node [ id 4 label "H" ] + node [ id 5 label "H" ] + node [ id 6 label "H" ] + node [ id 7 label "H" ] + node [ id 8 label "H" ] + node [ id 9 label "C" ] + node [ id 10 label "H" ] + node [ id 11 label "H" ] + node [ id 12 label "H" ] + edge [ source 0 target 1 label "=" ] + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 0 target 4 label "-" ] + edge [ source 0 target 5 label "-" ] + edge [ source 1 target 6 label "-" ] + edge [ source 2 target 7 label "-" ] + edge [ source 3 target 8 label "-" ] + edge [ source 3 target 9 label "-" ] + edge [ source 9 target 10 label "-" ] + edge [ source 9 target 11 label "-" ] + edge [ source 9 target 12 label "-" ] +] +""" +#, name="Pentadien") + +restswap = Rule.fromGMLString( +"""rule [ + left [ + edge [ source 1 target 2 label "=" ] + edge [ source 3 target 4 label "=" ] + ] + context [ + node [ id 1 label "C" ] + node [ id 2 label "C"] + node [ id 3 label "C"] + node [ id 4 label "C"] + ] + right [ + edge [ source 1 target 3 label "=" ] + edge [ source 2 target 4 label "=" ] + ] +]""" +) + +dielsalder = Rule.fromGMLString( +"""rule [ + left [ + edge [ source 1 target 2 label "=" ] + edge [ source 2 target 3 label "-" ] + edge [ source 3 target 4 label "=" ] + edge [ source 5 target 6 label "=" ] + + ] + context [ + node [ id 1 label "C" ] + node [ id 2 label "C"] + node [ id 3 label "C"] + node [ id 4 label "C"] + node [ id 5 label "C"] + node [ id 6 label "C"] + ] + right [ + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 3 target 4 label "-" ] + edge [ source 4 target 5 label "-" ] + edge [ source 5 target 6 label "-" ] + edge [ source 6 target 1 label "-" ] + ] +]""" +) + +#MCB? horton vs DePina, einfach Depth first ob nach vier cyclus +def cyclesizes(g): + nxGraph = GraphObj(g).nx_graph + cycles = sorted(list(nx.chordless_cycles(nxGraph))) + cycle_lengths = [len(x) for x in cycles] + #Find all cycles and list their sizes + #Only for C + #Chordless, elemenatary + #Is there a inbuild way to ignore all H? + # + if cycles == []: + return False + #One Atom can't be in four different cycles + for vertice in nxGraph.nodes: + counter = 0 + for cycle in cycles: + if vertice in cycle: + counter += 1 + if counter >= 4: + return True + #One Cycle can't overlapp with another on over 2 Connection points + for cycle in cycles: + for cycleref in cycles: + if cycle != cycleref and len(list(set(cycle) & set(cycleref))) >= 3: + return True + #Only cordless cycles of length 5,6 and 7 are acceptable + if min(cycle_lengths) >= 5 and max(cycle_lengths) <=7: + return False + return True + +#Disallows Allenes (Two Doublebonds on same carbon) +def doubledoublebond(g): + for vertice in g.vertices: + if (vertice.stringLabel == "C" and vertice.degree == 2): # and vertice.edges.label == ["=", "="] + for edge in vertice.incidentEdges: + print(type(edge.bondType)) + bonds = [type(edge.bondType) for edge in vertice.incidentEdges] + if (bonds[0] == bonds[1] and len(bonds) == 2): + return True + return False + + +def restriction(dg): + #Only rings 5 to 7 + #No neighbouring double bonds. Is this already in MOD? Can this even happen? + #No Carbon in 3 rings? Better 4? + for a in dg.right: + if a.vLabelCount("C") > 10: + return False + if doubledoublebond(a): + return False + if cyclesizes(a): + return False + return True + +flowPrinter = FlowPrinter() +flowPrinter.printUnfiltered = False + +postSection("Loaded Graphs") +for a in inputGraphs: + a.print() +postSection("Loaded Rules") +for a in inputRules: + a.print() + + +dg = DG(graphDatabase=inputGraphs) +dg.build().execute( + addSubset(inputGraphs) + >> rightPredicate[ + restriction + ]( + repeat(revive(inputRules)) #Revive not necessary + ) + +) +dg.print() +postSection("Product Graphs") +for a in dg.vertices: + a.graph.print() + +#flow = Flow(dg) +#flow.addSource(butadien) +#flow.findSolutions() +#flow.solutions.list() +#flow.solutions.print(flowPrinter) + + +sys.exit(0) + +rc = rcEvaluator(inputRules) +for dRef in dg.derivations: + der = dRef.derivation + educt = rcId(der.left[0]) + for i in range(1, len(der.left)): + educt = educt *rcParallel* rcId(der.left[i]) + product = rcId(der.right[0]) + for i in range(1, len(der.right)): + product = product *rcParallel* rcId(der.right[i]) + rcExp = educt *rcSuper(allowPartial=False)* der.rule *rcSuper(allowPartial=False)* product + res = rc.eval(rcExp) + dRef.print() + for a in res: + a.print() + a.printGML() \ No newline at end of file diff --git a/mod/butadien/graph.py b/mod/butadien/graph.py new file mode 100644 index 0000000..a787e41 --- /dev/null +++ b/mod/butadien/graph.py @@ -0,0 +1,136 @@ +import mod +import networkx as nx +import time + +# from equilibrator_assets.generate_compound import create_compounds + +# each vertex has attribute "label" which stores string label +# each edge has attribute "bond" which stores bond order (-, =, etc.) + +class GraphObj: + def __init__(self, graph, verbose=False): + self._verbose = verbose + + mod_lg_class = mod.libpymod.Rule.LeftGraph + mod_rg_class = mod.libpymod.Rule.RightGraph + mod_g_class = mod.libpymod.Graph + if isinstance(graph, mod_lg_class) or isinstance(graph, mod_rg_class) or isinstance(graph, mod_g_class): + self._nx_graph = self.mod_to_nx_graph(graph) + self._gml_string = self.nx_graph_to_GML_string(self._nx_graph) + elif isinstance(graph, nx.Graph): + self._nx_graph = graph + self._gml_string = self.nx_graph_to_GML_string(self._nx_graph) + elif isinstance(graph, str): + self._gml_string = graph + self._nx_graph = self.GML_to_nx_graph(self._gml_string) + else: + print(f"ERROR: Graph class cannot identify graph type in constructor: {type(graph)}") + + # up to the caller to check the number of components created + ccs = [self._nx_graph.subgraph(c).copy() for c in nx.connected_components(self._nx_graph)] + ccs = sorted(ccs, key=len, reverse=True) + self._num_components = len(ccs) + if self._num_components > 1: + self._components = [Graph(c) for c in ccs] + else: + self._components = [self] + + if self._num_components == 1: + self._mod_graph = self.nx_graph_to_mod(self._nx_graph) + + # self._equ_compounds = self.equ_compounds() + + def equ_compounds(self): + smiles = [c.mod_graph.smiles for c in self._components] + equ_comps = create_compounds(smiles, mol_format="smiles", bypass_chemaxon=True, save_empty_compounds=True) + return equ_comps + + def mod_to_nx_graph(self, modGraph: mod.Graph): + g = nx.Graph() + + for v in modGraph.vertices: + g.add_node(int(v.id), label=str(v.stringLabel), modID=int(v.id)) + + for e in modGraph.edges: + g.add_edge(int(e.source.id), int(e.target.id), bond=str(e.bondType)) + + return g + + def nx_graph_to_mod(self, nxGraph): + try: + return mod.graphGMLString(self.nx_graph_to_GML_string(nxGraph)) + except mod.libpymod.InputError: # graph is not connected probably + if self._verbose: + print("Error converting nxGraph to mod graph. Likely graph is not connected. This will not affect rule generation.") + return None + + def GML_to_nx_graph(self, gml_string): + g = nx.Graph() + lines = gml_string.split("\n") + for line in lines: + tokens = line.split() + if len(tokens) < 2: + continue + if tokens[0] == "node": + (_,_,_, mid, _, l, _) = tokens + g.add_node(int(mid), label=l[1:-1], modID=int(mid)) + elif tokens[0] == "edge": + (_, _, _, u, _, v, _, bondOrder, _ ) = tokens + g.add_edge(int(u), int(v), bond=bondOrder[1:-1]) + return g + + def nx_graph_to_GML_string(self, nxGraph): + out = [] + out.append("graph [") + + out.extend([f"\t\tnode [ id {nxGraph.nodes[node]['modID']} label \"{nxGraph.nodes[node]['label']}\" ]" for node in nxGraph.nodes]) + out.extend([f"\t\tedge [ source {u} target {v} label \"{nxGraph[u][v]['bond']}\" ]" + for (u,v) in nxGraph.edges]) + + out.append("]") + + return "\n".join(out) + + @property + def nx_graph(self): + return self._nx_graph + + @property + def mod_graph(self): + if self.num_components == 1: + return self._mod_graph + else: + return mod.graphGMLString(self.gml_string) + + @property + def gml_string(self): + return self._gml_string + + @property + def gml(self): + return self._gml_string + + @property + def edges(self): + return self._nx_graph.edges + + @property + def nodes(self): + return self._nx_graph.nodes + + @property + def num_components(self): + return self._num_components + + @property + def connected_components(self): + return self._components + + def mod_print(self): + self._mod_graph.print() + + def __str__(self) -> str: + return self._gml_string + + def __hash__(self) -> int: + return hash(self._gml_string)