diff --git a/cinnabar/femap.py b/cinnabar/femap.py index 221c147..93f079d 100644 --- a/cinnabar/femap.py +++ b/cinnabar/femap.py @@ -492,69 +492,3 @@ def draw_graph(self, title: str = "", filename: Union[str, None] = None): else: plt.savefig(filename, bbox_inches="tight") - def get_cycle_closure(self): - """Calculate the sum of DDG along all ligand cycles as a measure of convergence.""" - network = self.to_legacy_graph() - y = [x[2]["calc_DDG"] for x in network.edges(data=True)] - - # Find all ligand cycles - cycles = sorted(nx.simple_cycles(network.to_undirected())) - edges = network.edges - - # Loop over cycles, calculate sum of DG along cycle - dict = {} - for cycle in cycles: - - # Store DDG values along the cycle - sum_ddgs = 0 - for inx, ligand in enumerate(cycle): - if inx < len(cycle) - 1: - ligA = ligand - ligB = cycle[inx + 1] - # Last ligand is connected to first ligand - else: - ligA = ligand - ligB = cycle[0] - # depending on the direction the edge was calculated, - # the sign of the DDG has to change - if (ligA, ligB) in list(edges): - ddg = y[list(edges).index((ligA, ligB))] - elif (ligB, ligA) in edges: - ddg = -y[list(edges).index((ligB, ligA))] - # sum up DDGs along cycle - sum_ddgs += ddg - - # divide by sqrt of number of ligands in the cycle - # to get cycle closure error PER EDGE - cc = abs(sum_ddgs / math.sqrt(len(cycle))) - # Store cycle and cycle closure in dict - dict[','.join(cycle)] = round(cc, 2) - - # Sort cycle closure from high to low - sorted_list = sorted(dict.items(), key=lambda x: x[1], reverse=True) - - return sorted_list - - def store_cycle_closure_to_csv(self, file='cycle_closure.csv'): - """Save cycle closure, sorted from highest cycle closure to lowest - in a csv file.""" - sorted_list = get_cycle_closure(self) - - # CSV file to store results - f = open(file, 'w') - writer = csv.writer(f, lineterminator='\n') - header = '# ligands in cycle, sum(DDGs) / ' \ - 'sqrt(number of ligands in cycle)\n' - f.write(header) - writer.writerows(sorted_list) - f.close() - - return - - def plot_hist_cycle_closure(self, file='cycle_closure_hist.png'): - - sorted_list = get_cycle_closure(self) - - plt.hist([s[1] for s in sorted_list]) - plt.xlabel('Cycle closure in kcal/mol') - plt.savefig(file, bbox_inches="tight")