From 741039acd7c9d053c0f01fe77fd7aa0a8b715c89 Mon Sep 17 00:00:00 2001 From: JohannesKarwou Date: Fri, 8 Dec 2023 12:01:49 +0100 Subject: [PATCH 1/7] remove dependency on tf_routes --- devtools/conda-envs/fep_env.yaml | 1 - transformato/mutate.py | 253 +++++++++++++++++++++++++++---- 2 files changed, 220 insertions(+), 34 deletions(-) diff --git a/devtools/conda-envs/fep_env.yaml b/devtools/conda-envs/fep_env.yaml index 2c1a251e..b2fe8897 100644 --- a/devtools/conda-envs/fep_env.yaml +++ b/devtools/conda-envs/fep_env.yaml @@ -24,6 +24,5 @@ dependencies: - pytest-cov - codecov - pip: - - git+https://github.com/wiederm/tf_routes.git - git+https://github.com/ParmEd/ParmEd.git - git+https://github.com/wiederm/transformato_testsystems.git diff --git a/transformato/mutate.py b/transformato/mutate.py index 11bb9a54..3609cb11 100644 --- a/transformato/mutate.py +++ b/transformato/mutate.py @@ -526,49 +526,236 @@ def _match_terminal_dummy_atoms_between_common_cores( return (lj_default_cc1, lj_default_cc2) - @staticmethod - def _calculate_order_of_LJ_mutations( + def change_route_cycles(route, cycledict, degreedict, weightdict, G): + """ + preliminary mutation list is sorted using cycle and degree dictionary + currently used in _calculate_order_of_LJ_mutations_new + ---- + Args: + route: original mutation route + cycledict: dict of cycle participation of atoms + degreedict: dict of degree of atoms + weightdict: dict of weight of atoms + G: nx-graph of molecule + ---- + returns reordered array of the mutation route + """ + + for i in range(len(route) - 1): + routedict = route[i] + routeweight = weightdict.get(route[i]) + + routecycleval = cycledict.get(route[i]) + routedegreeval = degreedict.get(route[i]) + + for j in range(i, len(route)): + if routeweight == weightdict[route[j]]: + # if nodes have same weight (i.e. distance from root), the node participating in more cycles is removed later + + if routecycleval > cycledict[route[j]] or ( + routecycleval == cycledict[route[j]] + and routedegreeval > degreedict[route[j]] + ): + idx1 = route.index(route[i]) + idx2 = route.index(route[j]) + route[idx1], route[idx2] = route[idx2], route[idx1] + continue + + # if nodes have same weight (i.e. distance from root) and same cycle participation number, the node which has more neighbours already removed is removed earlier + + if routecycleval == cycledict[route[j]]: + edgesi = G.edges(routedict) + edgesj = G.edges(route[j]) + + iedgecounter = 0 + for edge in edgesi: + if edge[1] in route[0:i] or edge[0] in route[0:i]: + iedgecounter = iedgecounter + 1 + + jedgecounter = 0 + for edge in edgesj: + if edge[1] in route[0:i] or edge[0] in route[0:i]: + jedgecounter = jedgecounter + 1 + + if iedgecounter < jedgecounter: + idx1 = route.index(route[i]) + idx2 = route.index(route[j]) + route[idx1], route[idx2] = route[idx2], route[idx1] + + return route + + def cycle_checks_nx(G, use_actual_weight_for_mod=False): + """ + cycle processing, can be used in _calculate_order_of_LJ_mutations_new_iter and ..._new_iter_change (default is cycle_checks_nx_v2) + -------- + returns nx-graph-object with updated weights (according to cycle participation of the atom) + """ + + # search cycles using networkx + cycles = nx.cycle_basis(G) + + from collections import Counter + + cdict = Counter(x for xs in cycles for x in set(xs)) + + # modify weighted graph: nodes participating in many cycles get lower weight + for i in cdict: + edg = G.edges(i) + for el in edg: + if use_actual_weight_for_mod == True: + G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * ( + G[el[0]][el[1]]["weight"] ** (1 / 2) + ) + else: + G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * 5 + + return G + + def cycle_checks(G): + """ + cycle processing dictionary and degree dictionary for preferential removal (atoms which neighbours already have been removed are removed earlier), currently used in _calculate_order_of_LJ_mutations_new (via change_route_cycles) + ---- + returns dictionary containing number of cycle participation (cdict) and dict containing degree of atom (degreedict) + """ + + # search cycles using networkx + cycles = nx.cycle_basis(G) + + # alternatively, using rdkit + # ri = mol.GetRingInfo() + # cyclesrdkit = ri.AtomRings() + + import collections + from collections import Counter + + cdict = Counter(x for xs in cycles for x in set(xs)) + # cdictrdkit = Counter(x for xs in cyclesrdkit for x in set(xs)) + + # add atoms with no cycle participation + for key in G.nodes: + if key not in cdict: + cdict[key] = 0 + + degreedict = G.degree() + degreedict = {node: val for (node, val) in degreedict} + + return cdict, degreedict + + def exclude_Hs_from_mutations(connected_dummy_regions: list, G: nx.Graph): + """ + hydrogens are removed from the networkx-graph-representation and the list of connected dummy regions + ---- + Args: + connected_dummy_regions: list of connected dummy regions + G: nx-graph of molecule + ---- + returns list of connected dummy regions and networkx-graph without hydrogens + """ + + G_hydrogens = [x for x, y in G.nodes(data=True) if y["atom_type"] == "H"] + + G.remove_nodes_from(G_hydrogens) + connected_dummy_regions_copy = connected_dummy_regions + for hydroindex in G_hydrogens: + for indexregion, region in enumerate(connected_dummy_regions): + if hydroindex in region: + connected_dummy_regions_copy[indexregion].remove(hydroindex) + + return connected_dummy_regions_copy, G + + def _calculate_order_of_LJ_mutations_with_bfs( + self, connected_dummy_regions: list, match_terminal_atoms: dict, G: nx.Graph, + cyclecheck=True, + ordercycles=True, + exclude_Hs=True, ) -> list: - try: - from tf_routes.routes import ( - _calculate_order_of_LJ_mutations_new as _calculate_order_of_LJ_mutations_with_bfs, - ) + """ + bfs/djikstra-algorithm applied once for route (without iterations) + ----- + cyclecheck: updates weights according to cycle participation (should always be set to True) + ordercheck: if there is no possibility to decide between two nodes - i.e. the weight would be the exactly the same - weight updating according to preferential removal decides that the node in which neighbourhood nodes already have been removed is removed next + exclude_Hs: if True, hydrogens are removed before the mutation algorithm is applied - necessary for usual Transformato workflow + """ - return _calculate_order_of_LJ_mutations_with_bfs( - connected_dummy_regions, match_terminal_atoms, G + if exclude_Hs == True: + connected_dummy_regions, G = self.exclude_Hs_from_mutations( + connected_dummy_regions, G ) - except ModuleNotFoundError: - ordered_LJ_mutations = [] - for real_atom in match_terminal_atoms: - for dummy_atom in match_terminal_atoms[real_atom]: - for connected_dummy_region in connected_dummy_regions: - # stop at connected dummy region with specific dummy_atom in it - if dummy_atom not in connected_dummy_region: - continue - - G_dummy = G.copy() - # delete all nodes not in dummy region - remove_nodes = [ - node - for node in G.nodes() - if node not in connected_dummy_region - ] - for remove_node in remove_nodes: - G_dummy.remove_node(remove_node) + ordered_LJ_mutations = [] + + for real_atom in match_terminal_atoms: + for dummy_atom in match_terminal_atoms[real_atom]: + for connected_dummy_region in connected_dummy_regions: + # stop at connected dummy region with specific dummy_atom in it + if dummy_atom not in connected_dummy_region: + continue - # root is the dummy atom that connects the real region with the dummy region - root = dummy_atom + G_dummy = G.copy() + # delete all nodes not in dummy region + remove_nodes = [ + node for node in G.nodes() if node not in connected_dummy_region + ] + for remove_node in remove_nodes: + G_dummy.remove_node(remove_node) + + # root is the dummy atom that connects the real region with the dummy region + root = dummy_atom + + # process cycles + if cyclecheck == True and ordercycles == False: + G_dummy = self.cycle_checks_nx(G_dummy) - edges = list(nx.dfs_edges(G_dummy, source=root)) - nodes = [root] + [v for u, v in edges] - nodes.reverse() # NOTE: reverse the mutation - ordered_LJ_mutations.append(nodes) + # process cycles and correct order (according to 'preferential removal') + if cyclecheck == True and ordercycles == True: + cycledict, degreedict = self.cycle_checks(G_dummy) - return ordered_LJ_mutations + # dijkstra + ssource = nx.single_source_dijkstra( + G_dummy, source=root, weight="weight" + ) + # result of dijkstra algorithm is sorted + sortedssource = { + k: v + for k, v in sorted(ssource[0].items(), key=lambda item: item[1]) + } + + # get keys of sorted dict + sortedssource_edges = sortedssource.keys() + + sortedssource_edges_list = list(sortedssource_edges) + # sorted list contains the mutation route + nodes = sortedssource_edges_list + + # order has to be reversed - the most distant atom is the first to be removed + nodes.reverse() + + # sort nodes according to degree, cycle participation and removal order + if cyclecheck == True and ordercycles == True: + nodes = self.change_route_cycles( + nodes, cycledict, degreedict, sortedssource, G + ) + + logger.info("Final mutation route:") + logger.info(nodes) + ordered_LJ_mutations.append(nodes) + + return ordered_LJ_mutations + + @staticmethod + def _calculate_order_of_LJ_mutations( + self, + connected_dummy_regions: list, + match_terminal_atoms: dict, + G: nx.Graph, + ) -> list: + + return self._calculate_order_of_LJ_mutations_with_bfs( + connected_dummy_regions, match_terminal_atoms, G + ) def _check_for_lp( self, From 2da996a042bb7c99aa57d73e38be250e78435588 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Fri, 8 Dec 2023 11:03:20 +0000 Subject: [PATCH 2/7] Fix code style issues with Black --- transformato/mutate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/transformato/mutate.py b/transformato/mutate.py index 3609cb11..ea0dcba3 100644 --- a/transformato/mutate.py +++ b/transformato/mutate.py @@ -752,7 +752,6 @@ def _calculate_order_of_LJ_mutations( match_terminal_atoms: dict, G: nx.Graph, ) -> list: - return self._calculate_order_of_LJ_mutations_with_bfs( connected_dummy_regions, match_terminal_atoms, G ) From f4e71a53ca765160f2558eeacab4c01d47b0b974 Mon Sep 17 00:00:00 2001 From: JohannesKarwou Date: Mon, 11 Dec 2023 12:14:01 +0100 Subject: [PATCH 3/7] first clean up --- devtools/meta.yaml | 1 - transformato/analysis.py | 2 +- transformato/annihilation.py | 52 ---- transformato/charmm_factory.py | 2 +- transformato/constants.py | 60 ----- transformato/helper_functions.py | 232 ++++++++++++++++++ transformato/mutate.py | 169 ++----------- .../tests/test_alchemical_path_generation.py | 2 - transformato/tests/test_misc.py | 39 +-- transformato/tests/test_postprocessing.py | 9 +- transformato/tests/test_run_production.py | 3 +- transformato/tests/test_system_setup.py | 23 +- transformato/testsystems.py | 26 -- 13 files changed, 275 insertions(+), 345 deletions(-) delete mode 100644 transformato/annihilation.py delete mode 100644 transformato/constants.py create mode 100644 transformato/helper_functions.py delete mode 100644 transformato/testsystems.py diff --git a/devtools/meta.yaml b/devtools/meta.yaml index f141fda8..4def301c 100644 --- a/devtools/meta.yaml +++ b/devtools/meta.yaml @@ -42,7 +42,6 @@ test: - transformato - transformato.mutate - transformato.analysis - - transformato.constants - transformato.restraints - transformato.state - transformato.system diff --git a/transformato/analysis.py b/transformato/analysis.py index a154ff45..8063face 100644 --- a/transformato/analysis.py +++ b/transformato/analysis.py @@ -21,7 +21,7 @@ from openmm.app import CharmmPsfFile, Simulation from tqdm import tqdm -from transformato.constants import temperature +from transformato.helper_functions import temperature from transformato.utils import get_structure_name, isnotebook logger = logging.getLogger(__name__) diff --git a/transformato/annihilation.py b/transformato/annihilation.py deleted file mode 100644 index 5508ad95..00000000 --- a/transformato/annihilation.py +++ /dev/null @@ -1,52 +0,0 @@ -import logging -import parmed as pm - -import networkx as nx - - -logger = logging.getLogger(__name__) - - -def calculate_order_of_LJ_mutations_asfe(central_atoms: list, G: nx.Graph) -> list: - ordered_LJ_mutations = [] - root = central_atoms[0] - - final_order = [] - removearray = [] - removeG = nx.Graph() - G_dummy = G.copy() - - while len(G_dummy.nodes()) > 0: - G_origweights = G_dummy.copy() - - # dijkstra - ssource = nx.single_source_dijkstra(G_dummy, source=root, weight="weight") - - # find node with maximum distance - i.e next to be removed - max_node = max(ssource[0], key=ssource[0].get) - - # the currently most distant node is added to the final order - it is the next to be removed - final_order.extend([max_node]) - - # restore original weights - G_dummy = G_origweights - - # remove the most distant node (after added to the final_order-array); the next iteration will find the most distant node of the reduced graph - G_dummy.remove_node(max_node) - - # add to removeG - removeG.add_node(max_node) - removearray.append(max_node) - - # final_order contains mutation order for current dummy region - sortedssource_edges = final_order - - # sortedssource_edges_list already contains the nodes in right (reversed) order (including root) - nodes = sortedssource_edges - ordered_LJ_mutations.append(nodes) - - logger.info( - f"Turning off the LJ potential in the following order: {ordered_LJ_mutations}" - ) - - return ordered_LJ_mutations diff --git a/transformato/charmm_factory.py b/transformato/charmm_factory.py index d92a9cae..27355ac0 100644 --- a/transformato/charmm_factory.py +++ b/transformato/charmm_factory.py @@ -1,5 +1,5 @@ import datetime -from transformato.constants import temperature, charmm_gpu +from transformato.helper_functions import temperature, charmm_gpu from transformato.utils import check_switching_function from openmm import unit import os diff --git a/transformato/constants.py b/transformato/constants.py deleted file mode 100644 index 3ed58269..00000000 --- a/transformato/constants.py +++ /dev/null @@ -1,60 +0,0 @@ -import logging, sys -from openmm import unit - -logger = logging.getLogger(__name__) - -temperature = 303.15 * unit.kelvin -kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA -kT = kB * temperature -charmm_gpu = "" -# charmm_gpu = "domdec-gpu" # uncomment this if you want to use domdec-gpu -# charmm_gpu = "charmm-openmm" # uncomment this if you want to use charmm-openmm -this = sys.modules[__name__] -# we can explicitly make assignments on it -this.NUM_PROC = 1 - -##################################### -# config for tets -test_platform_openMM = "GPU" -test_platform_CHARMM = "CPU" -test_platform_override = True -##################################### - - -def initialize_NUM_PROC(n_proc): - if this.NUM_PROC == 1: - # also in local function scope. no scope specifier like global is needed - this.NUM_PROC = n_proc - else: - msg = "NUM_PROC is already initialized to {0}." - # raise RuntimeError(msg.format(this.NUM_PROC)) - print(msg) - - -def change_platform_to_test_platform(configuration: dict, engine: str): - if engine == "openMM": - change_to = test_platform_openMM - elif engine == "CHARMM": - change_to = test_platform_CHARMM - else: - raise RecursionError() - - if change_to.upper() == "GPU": - configuration["simulation"]["GPU"] = True - print("Setting platform to GPU") - elif change_to.upper() == "CPU": - configuration["simulation"]["GPU"] = False - print("Setting platform to CPU") - else: - raise RuntimeError("something went wrong") - - -def change_platform_to(configuration: dict, change_to: str): - if change_to.upper() == "GPU": - configuration["simulation"]["GPU"] = True - print("Setting platform to GPU") - elif change_to.upper() == "CPU": - configuration["simulation"]["GPU"] = False - print("Setting platform to CPU") - else: - raise RuntimeError("something went wrong") diff --git a/transformato/helper_functions.py b/transformato/helper_functions.py new file mode 100644 index 00000000..54f1f4c5 --- /dev/null +++ b/transformato/helper_functions.py @@ -0,0 +1,232 @@ +import logging +import parmed as pm +import networkx as nx +from openmm import unit +import sys + +logger = logging.getLogger(__name__) + + +#### Constants which might be needed sometimes! +temperature = 303.15 * unit.kelvin +kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA +kT = kB * temperature +charmm_gpu = "" +# charmm_gpu = "domdec-gpu" # uncomment this if you want to use domdec-gpu +# charmm_gpu = "charmm-openmm" # uncomment this if you want to use charmm-openmm +this = sys.modules[__name__] +# we can explicitly make assignments on it +this.NUM_PROC = 1 + + + +##################################### +# config for tets +test_platform_openMM = "GPU" +test_platform_CHARMM = "CPU" +test_platform_override = True +##################################### + + +def change_platform_to_test_platform(configuration: dict, engine: str): + if engine == "openMM": + change_to = test_platform_openMM + elif engine == "CHARMM": + change_to = test_platform_CHARMM + else: + raise RecursionError() + + if change_to.upper() == "GPU": + configuration["simulation"]["GPU"] = True + print("Setting platform to GPU") + elif change_to.upper() == "CPU": + configuration["simulation"]["GPU"] = False + print("Setting platform to CPU") + else: + raise RuntimeError("something went wrong") + + +def calculate_order_of_LJ_mutations_asfe(central_atoms: list, G: nx.Graph) -> list: + ordered_LJ_mutations = [] + root = central_atoms[0] + + final_order = [] + removearray = [] + removeG = nx.Graph() + G_dummy = G.copy() + + while len(G_dummy.nodes()) > 0: + G_origweights = G_dummy.copy() + + # dijkstra + ssource = nx.single_source_dijkstra(G_dummy, source=root, weight="weight") + + # find node with maximum distance - i.e next to be removed + max_node = max(ssource[0], key=ssource[0].get) + + # the currently most distant node is added to the final order - it is the next to be removed + final_order.extend([max_node]) + + # restore original weights + G_dummy = G_origweights + + # remove the most distant node (after added to the final_order-array); the next iteration will find the most distant node of the reduced graph + G_dummy.remove_node(max_node) + + # add to removeG + removeG.add_node(max_node) + removearray.append(max_node) + + # final_order contains mutation order for current dummy region + sortedssource_edges = final_order + + # sortedssource_edges_list already contains the nodes in right (reversed) order (including root) + nodes = sortedssource_edges + ordered_LJ_mutations.append(nodes) + + logger.info( + f"Turning off the LJ potential in the following order: {ordered_LJ_mutations}" + ) + + return ordered_LJ_mutations + + +def change_route_cycles(route, cycledict, degreedict, weightdict, G): + """ + preliminary mutation list is sorted using cycle and degree dictionary + currently used in _calculate_order_of_LJ_mutations_new + ---- + Args: + route: original mutation route + cycledict: dict of cycle participation of atoms + degreedict: dict of degree of atoms + weightdict: dict of weight of atoms + G: nx-graph of molecule + ---- + returns reordered array of the mutation route + """ + + for i in range(len(route) - 1): + routedict = route[i] + routeweight = weightdict.get(route[i]) + + routecycleval = cycledict.get(route[i]) + routedegreeval = degreedict.get(route[i]) + + for j in range(i, len(route)): + if routeweight == weightdict[route[j]]: + # if nodes have same weight (i.e. distance from root), the node participating in more cycles is removed later + + if routecycleval > cycledict[route[j]] or ( + routecycleval == cycledict[route[j]] + and routedegreeval > degreedict[route[j]] + ): + idx1 = route.index(route[i]) + idx2 = route.index(route[j]) + route[idx1], route[idx2] = route[idx2], route[idx1] + continue + + # if nodes have same weight (i.e. distance from root) and same cycle participation number, the node which has more neighbours already removed is removed earlier + + if routecycleval == cycledict[route[j]]: + edgesi = G.edges(routedict) + edgesj = G.edges(route[j]) + + iedgecounter = 0 + for edge in edgesi: + if edge[1] in route[0:i] or edge[0] in route[0:i]: + iedgecounter = iedgecounter + 1 + + jedgecounter = 0 + for edge in edgesj: + if edge[1] in route[0:i] or edge[0] in route[0:i]: + jedgecounter = jedgecounter + 1 + + if iedgecounter < jedgecounter: + idx1 = route.index(route[i]) + idx2 = route.index(route[j]) + route[idx1], route[idx2] = route[idx2], route[idx1] + + return route + + +def cycle_checks_nx(G, use_actual_weight_for_mod=False): + """ + cycle processing, can be used in _calculate_order_of_LJ_mutations_new_iter and ..._new_iter_change (default is cycle_checks_nx_v2) + -------- + returns nx-graph-object with updated weights (according to cycle participation of the atom) + """ + + # search cycles using networkx + cycles = nx.cycle_basis(G) + + from collections import Counter + + cdict = Counter(x for xs in cycles for x in set(xs)) + + # modify weighted graph: nodes participating in many cycles get lower weight + for i in cdict: + edg = G.edges(i) + for el in edg: + if use_actual_weight_for_mod == True: + G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * ( + G[el[0]][el[1]]["weight"] ** (1 / 2) + ) + else: + G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * 5 + + return G + + +def cycle_checks(G): + """ + cycle processing dictionary and degree dictionary for preferential removal (atoms which neighbours already have been removed are removed earlier), currently used in _calculate_order_of_LJ_mutations_new (via change_route_cycles) + ---- + returns dictionary containing number of cycle participation (cdict) and dict containing degree of atom (degreedict) + """ + + # search cycles using networkx + cycles = nx.cycle_basis(G) + + # alternatively, using rdkit + # ri = mol.GetRingInfo() + # cyclesrdkit = ri.AtomRings() + + import collections + from collections import Counter + + cdict = Counter(x for xs in cycles for x in set(xs)) + # cdictrdkit = Counter(x for xs in cyclesrdkit for x in set(xs)) + + # add atoms with no cycle participation + for key in G.nodes: + if key not in cdict: + cdict[key] = 0 + + degreedict = G.degree() + degreedict = {node: val for (node, val) in degreedict} + + return cdict, degreedict + + +def exclude_Hs_from_mutations(connected_dummy_regions: list, G: nx.Graph): + """ + hydrogens are removed from the networkx-graph-representation and the list of connected dummy regions + ---- + Args: + connected_dummy_regions: list of connected dummy regions + G: nx-graph of molecule + ---- + returns list of connected dummy regions and networkx-graph without hydrogens + """ + + G_hydrogens = [x for x, y in G.nodes(data=True) if y["atom_type"] == "H"] + + G.remove_nodes_from(G_hydrogens) + connected_dummy_regions_copy = connected_dummy_regions + for hydroindex in G_hydrogens: + for indexregion, region in enumerate(connected_dummy_regions): + if hydroindex in region: + connected_dummy_regions_copy[indexregion].remove(hydroindex) + + return connected_dummy_regions_copy, G diff --git a/transformato/mutate.py b/transformato/mutate.py index ea0dcba3..62351b97 100644 --- a/transformato/mutate.py +++ b/transformato/mutate.py @@ -17,7 +17,13 @@ IPythonConsole.ipython_useSVG = True # Change output to SVG from transformato.system import SystemStructure -from transformato.annihilation import calculate_order_of_LJ_mutations_asfe +from transformato.helper_functions import ( + calculate_order_of_LJ_mutations_asfe, + cycle_checks_nx, + cycle_checks, + exclude_Hs_from_mutations, + change_route_cycles, +) logger = logging.getLogger(__name__) @@ -526,144 +532,7 @@ def _match_terminal_dummy_atoms_between_common_cores( return (lj_default_cc1, lj_default_cc2) - def change_route_cycles(route, cycledict, degreedict, weightdict, G): - """ - preliminary mutation list is sorted using cycle and degree dictionary - currently used in _calculate_order_of_LJ_mutations_new - ---- - Args: - route: original mutation route - cycledict: dict of cycle participation of atoms - degreedict: dict of degree of atoms - weightdict: dict of weight of atoms - G: nx-graph of molecule - ---- - returns reordered array of the mutation route - """ - - for i in range(len(route) - 1): - routedict = route[i] - routeweight = weightdict.get(route[i]) - - routecycleval = cycledict.get(route[i]) - routedegreeval = degreedict.get(route[i]) - - for j in range(i, len(route)): - if routeweight == weightdict[route[j]]: - # if nodes have same weight (i.e. distance from root), the node participating in more cycles is removed later - - if routecycleval > cycledict[route[j]] or ( - routecycleval == cycledict[route[j]] - and routedegreeval > degreedict[route[j]] - ): - idx1 = route.index(route[i]) - idx2 = route.index(route[j]) - route[idx1], route[idx2] = route[idx2], route[idx1] - continue - - # if nodes have same weight (i.e. distance from root) and same cycle participation number, the node which has more neighbours already removed is removed earlier - - if routecycleval == cycledict[route[j]]: - edgesi = G.edges(routedict) - edgesj = G.edges(route[j]) - - iedgecounter = 0 - for edge in edgesi: - if edge[1] in route[0:i] or edge[0] in route[0:i]: - iedgecounter = iedgecounter + 1 - - jedgecounter = 0 - for edge in edgesj: - if edge[1] in route[0:i] or edge[0] in route[0:i]: - jedgecounter = jedgecounter + 1 - - if iedgecounter < jedgecounter: - idx1 = route.index(route[i]) - idx2 = route.index(route[j]) - route[idx1], route[idx2] = route[idx2], route[idx1] - - return route - - def cycle_checks_nx(G, use_actual_weight_for_mod=False): - """ - cycle processing, can be used in _calculate_order_of_LJ_mutations_new_iter and ..._new_iter_change (default is cycle_checks_nx_v2) - -------- - returns nx-graph-object with updated weights (according to cycle participation of the atom) - """ - - # search cycles using networkx - cycles = nx.cycle_basis(G) - - from collections import Counter - - cdict = Counter(x for xs in cycles for x in set(xs)) - - # modify weighted graph: nodes participating in many cycles get lower weight - for i in cdict: - edg = G.edges(i) - for el in edg: - if use_actual_weight_for_mod == True: - G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * ( - G[el[0]][el[1]]["weight"] ** (1 / 2) - ) - else: - G[el[0]][el[1]]["weight"] = G[el[0]][el[1]]["weight"] - cdict[i] * 5 - - return G - - def cycle_checks(G): - """ - cycle processing dictionary and degree dictionary for preferential removal (atoms which neighbours already have been removed are removed earlier), currently used in _calculate_order_of_LJ_mutations_new (via change_route_cycles) - ---- - returns dictionary containing number of cycle participation (cdict) and dict containing degree of atom (degreedict) - """ - - # search cycles using networkx - cycles = nx.cycle_basis(G) - - # alternatively, using rdkit - # ri = mol.GetRingInfo() - # cyclesrdkit = ri.AtomRings() - - import collections - from collections import Counter - - cdict = Counter(x for xs in cycles for x in set(xs)) - # cdictrdkit = Counter(x for xs in cyclesrdkit for x in set(xs)) - - # add atoms with no cycle participation - for key in G.nodes: - if key not in cdict: - cdict[key] = 0 - - degreedict = G.degree() - degreedict = {node: val for (node, val) in degreedict} - - return cdict, degreedict - - def exclude_Hs_from_mutations(connected_dummy_regions: list, G: nx.Graph): - """ - hydrogens are removed from the networkx-graph-representation and the list of connected dummy regions - ---- - Args: - connected_dummy_regions: list of connected dummy regions - G: nx-graph of molecule - ---- - returns list of connected dummy regions and networkx-graph without hydrogens - """ - - G_hydrogens = [x for x, y in G.nodes(data=True) if y["atom_type"] == "H"] - - G.remove_nodes_from(G_hydrogens) - connected_dummy_regions_copy = connected_dummy_regions - for hydroindex in G_hydrogens: - for indexregion, region in enumerate(connected_dummy_regions): - if hydroindex in region: - connected_dummy_regions_copy[indexregion].remove(hydroindex) - - return connected_dummy_regions_copy, G - - def _calculate_order_of_LJ_mutations_with_bfs( + def _calculate_order_of_LJ_mutations( self, connected_dummy_regions: list, match_terminal_atoms: dict, @@ -673,15 +542,16 @@ def _calculate_order_of_LJ_mutations_with_bfs( exclude_Hs=True, ) -> list: """ - bfs/djikstra-algorithm applied once for route (without iterations) + bfs/djikstra-algorithm applied to calculate the ordere for turning of the LJ interactions of the heavy atoms ----- + most functions for theses options are imported from the helper_functions file cyclecheck: updates weights according to cycle participation (should always be set to True) ordercheck: if there is no possibility to decide between two nodes - i.e. the weight would be the exactly the same - weight updating according to preferential removal decides that the node in which neighbourhood nodes already have been removed is removed next exclude_Hs: if True, hydrogens are removed before the mutation algorithm is applied - necessary for usual Transformato workflow """ if exclude_Hs == True: - connected_dummy_regions, G = self.exclude_Hs_from_mutations( + connected_dummy_regions, G = exclude_Hs_from_mutations( connected_dummy_regions, G ) @@ -707,11 +577,11 @@ def _calculate_order_of_LJ_mutations_with_bfs( # process cycles if cyclecheck == True and ordercycles == False: - G_dummy = self.cycle_checks_nx(G_dummy) + G_dummy = cycle_checks_nx(G_dummy) # process cycles and correct order (according to 'preferential removal') if cyclecheck == True and ordercycles == True: - cycledict, degreedict = self.cycle_checks(G_dummy) + cycledict, degreedict = cycle_checks(G_dummy) # dijkstra ssource = nx.single_source_dijkstra( @@ -735,7 +605,7 @@ def _calculate_order_of_LJ_mutations_with_bfs( # sort nodes according to degree, cycle participation and removal order if cyclecheck == True and ordercycles == True: - nodes = self.change_route_cycles( + nodes = change_route_cycles( nodes, cycledict, degreedict, sortedssource, G ) @@ -745,17 +615,6 @@ def _calculate_order_of_LJ_mutations_with_bfs( return ordered_LJ_mutations - @staticmethod - def _calculate_order_of_LJ_mutations( - self, - connected_dummy_regions: list, - match_terminal_atoms: dict, - G: nx.Graph, - ) -> list: - return self._calculate_order_of_LJ_mutations_with_bfs( - connected_dummy_regions, match_terminal_atoms, G - ) - def _check_for_lp( self, odered_connected_dummy_regions_cc_with_lp: list, diff --git a/transformato/tests/test_alchemical_path_generation.py b/transformato/tests/test_alchemical_path_generation.py index 55396161..24c86395 100644 --- a/transformato/tests/test_alchemical_path_generation.py +++ b/transformato/tests/test_alchemical_path_generation.py @@ -193,8 +193,6 @@ def test_generate_alchemical_path_for_neopentane_common_core(): @pytest.mark.rsfe def test_generate_alchemical_path_for_methanol_common_core(): - from transformato_testsystems.testsystems import perform_generic_mutation - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/test-methanol-methane-rsfe.yaml", input_dir=get_testsystems_dir(), diff --git a/transformato/tests/test_misc.py b/transformato/tests/test_misc.py index f1d612e2..80ca63d1 100644 --- a/transformato/tests/test_misc.py +++ b/transformato/tests/test_misc.py @@ -1,7 +1,7 @@ """ Unit and regression test for the transformato package. """ -from transformato.constants import temperature as T +from transformato.helper_functions import temperature from openmm import unit import numpy as np import mdtraj as md @@ -17,20 +17,20 @@ def test_reduced_energy(): # with openMM generated traj evaluated with openMM e = -41264.39524669979 * unit.kilojoule_per_mole - rE = return_reduced_potential(e, volume=0, temperature=T) + rE = return_reduced_potential(e, volume=0, temperature=temperature) print(rE) assert np.isclose(rE, -16371.30301422) # with openMM generated traj evaluated with CHARMM e = -1099.41855 * unit.kilocalorie_per_mole - rE = return_reduced_potential(e, volume=0, temperature=T) + rE = return_reduced_potential(e, volume=0, temperature=temperature) print(rE) assert np.isclose(rE, -1824.9982986086145) # energy term in CHARMM traj # DYNA> 0 0.00000 -7775.74490 1914.51007 -9690.25498 377.14828 e = -9690.25498 * unit.kilocalorie_per_mole # ENERGgy - rE = return_reduced_potential(e, volume=0, temperature=T) + rE = return_reduced_potential(e, volume=0, temperature=temperature) print(rE) assert np.isclose(rE, -16085.501605902184) @@ -43,37 +43,6 @@ def test_convert_to_kT(): rE = e * beta assert np.isclose(rE, -16371.30301422) - -def test_change_platform(): - from transformato.constants import ( - change_platform_to_test_platform, - test_platform_openMM, - test_platform_CHARMM, - ) - - configuration = load_config_yaml( - config=f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe.yaml", - input_dir=".", - output_dir=get_test_output_dir(), - ) - - change_platform_to_test_platform(configuration, engine="openMM") - print(configuration["simulation"]["GPU"]) - print(test_platform_openMM) - if test_platform_openMM.upper() == "CPU": - assert configuration["simulation"]["GPU"] == False - else: - assert configuration["simulation"]["GPU"] == True - - change_platform_to_test_platform(configuration, engine="CHARMM") - print(configuration["simulation"]["GPU"]) - print(test_platform_openMM) - if test_platform_CHARMM.upper() == "CPU": - assert configuration["simulation"]["GPU"] == False - else: - assert configuration["simulation"]["GPU"] == True - - def test_scaling(): for i in np.linspace(1, 0, 11): f = max((1 - ((1 - i) * 2)), 0.0) diff --git a/transformato/tests/test_postprocessing.py b/transformato/tests/test_postprocessing.py index d393f5d4..61bb54c5 100644 --- a/transformato/tests/test_postprocessing.py +++ b/transformato/tests/test_postprocessing.py @@ -6,9 +6,6 @@ import numpy as np import pytest import logging -from transformato.constants import ( - initialize_NUM_PROC, -) from transformato.utils import postprocessing # read in specific topology with parameters @@ -23,9 +20,6 @@ from transformato.tests.paths import get_test_output_dir from transformato_testsystems.testsystems import get_testsystems_dir -# rbfe_test_systemes_generated = os.path.isdir("data/2OJ9-original-2OJ9-tautomer-rbfe") - - ########################################### # 2OJ9-tautomer system ########################################### @@ -314,7 +308,6 @@ def test_2oj9_postprocessing_with_different_engines(): reason="Skipping tests that cannot pass in github actions", ) def test_2oj9_postprocessing_with_openMM(): - initialize_NUM_PROC(1) conf = f"{get_testsystems_dir()}/config/test-2oj9-tautomer-pair-rsfe.yaml" configuration = load_config_yaml( @@ -598,7 +591,7 @@ def test_acetylacetone_postprocessing_different_engines(): reason="Skipping tests that cannot pass in github actions", ) def test_acetylacetone_calculate_rsfe_with_different_engines_only_vacuum(): - from ..constants import kT + from transformato.helper_functions import kT from openmm import unit conf = f"{get_testsystems_dir()}/config/test-acetylacetone-tautomer-rsfe.yaml" diff --git a/transformato/tests/test_run_production.py b/transformato/tests/test_run_production.py index 8a683e72..2bc4ce5b 100644 --- a/transformato/tests/test_run_production.py +++ b/transformato/tests/test_run_production.py @@ -10,7 +10,7 @@ from transformato import load_config_yaml from transformato.tests.paths import get_test_output_dir from transformato.utils import run_simulation -from transformato.constants import change_platform_to_test_platform +from transformato.helper_functions import change_platform_to_test_platform from transformato_testsystems.testsystems import get_testsystems_dir @@ -225,7 +225,6 @@ def test_run_2OJ9_tautomer_pair_with_openMM(caplog): ) def test_run_2OJ9_tautomer_pair_with_CHARMM(caplog): caplog.set_level(logging.WARNING) - from transformato.constants import change_platform_to from .test_mutation import setup_2OJ9_tautomer_pair_rsfe workdir = get_test_output_dir() diff --git a/transformato/tests/test_system_setup.py b/transformato/tests/test_system_setup.py index 17b24a51..ce434395 100644 --- a/transformato/tests/test_system_setup.py +++ b/transformato/tests/test_system_setup.py @@ -26,6 +26,27 @@ warnings.filterwarnings("ignore", module="parmed") +def perform_generic_mutation(configuration: dict): + s1 = SystemStructure(configuration, "structure1") + s2 = SystemStructure(configuration, "structure2") + + s1_to_s2 = ProposeMutationRoute(s1, s2) + s1_to_s2.calculate_common_core() + + mutation_list = s1_to_s2.generate_mutations_to_common_core_for_mol2() + i = IntermediateStateFactory( + system=s2, + configuration=configuration, + ) + + perform_mutations( + configuration=configuration, + i=i, + mutation_list=mutation_list, + nr_of_mutation_steps_charge=1, + ) + + return i.output_files def test_read_yaml(): """Sample test, will check ability to read yaml files""" @@ -120,7 +141,6 @@ def test_initialize_systems(caplog): @pytest.mark.rsfe def test_setup_system_for_methane_common_core(): - from transformato_testsystems.testsystems import perform_generic_mutation print(get_testsystems_dir()) @@ -171,7 +191,6 @@ def test_setup_system_for_toluene_common_core_with_HMR(): @pytest.mark.rsfe def test_setup_system_for_methane_common_core_with_HMR(): - from transformato_testsystems.testsystems import perform_generic_mutation configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe-HMR.yaml", diff --git a/transformato/testsystems.py b/transformato/testsystems.py deleted file mode 100644 index c420bc07..00000000 --- a/transformato/testsystems.py +++ /dev/null @@ -1,26 +0,0 @@ -from transformato.mutate import ProposeMutationRoute, perform_mutations -from transformato.state import IntermediateStateFactory -from transformato.system import SystemStructure - - -def perform_generic_mutation(configuration: dict): - s1 = SystemStructure(configuration, "structure1") - s2 = SystemStructure(configuration, "structure2") - - s1_to_s2 = ProposeMutationRoute(s1, s2) - s1_to_s2.calculate_common_core() - - mutation_list = s1_to_s2.generate_mutations_to_common_core_for_mol2() - i = IntermediateStateFactory( - system=s2, - configuration=configuration, - ) - - perform_mutations( - configuration=configuration, - i=i, - mutation_list=mutation_list, - nr_of_mutation_steps_charge=1, - ) - - return i.output_files From fc4200db745c140e7396915288fe5e26fefb606a Mon Sep 17 00:00:00 2001 From: Lint Action Date: Mon, 11 Dec 2023 11:14:52 +0000 Subject: [PATCH 4/7] Fix code style issues with Black --- transformato/helper_functions.py | 1 - transformato/tests/test_misc.py | 1 + transformato/tests/test_postprocessing.py | 2 +- transformato/tests/test_system_setup.py | 4 ++-- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/transformato/helper_functions.py b/transformato/helper_functions.py index 54f1f4c5..29e9f607 100644 --- a/transformato/helper_functions.py +++ b/transformato/helper_functions.py @@ -19,7 +19,6 @@ this.NUM_PROC = 1 - ##################################### # config for tets test_platform_openMM = "GPU" diff --git a/transformato/tests/test_misc.py b/transformato/tests/test_misc.py index 80ca63d1..25d4a4ec 100644 --- a/transformato/tests/test_misc.py +++ b/transformato/tests/test_misc.py @@ -43,6 +43,7 @@ def test_convert_to_kT(): rE = e * beta assert np.isclose(rE, -16371.30301422) + def test_scaling(): for i in np.linspace(1, 0, 11): f = max((1 - ((1 - i) * 2)), 0.0) diff --git a/transformato/tests/test_postprocessing.py b/transformato/tests/test_postprocessing.py index 61bb54c5..aa148059 100644 --- a/transformato/tests/test_postprocessing.py +++ b/transformato/tests/test_postprocessing.py @@ -20,6 +20,7 @@ from transformato.tests.paths import get_test_output_dir from transformato_testsystems.testsystems import get_testsystems_dir + ########################################### # 2OJ9-tautomer system ########################################### @@ -308,7 +309,6 @@ def test_2oj9_postprocessing_with_different_engines(): reason="Skipping tests that cannot pass in github actions", ) def test_2oj9_postprocessing_with_openMM(): - conf = f"{get_testsystems_dir()}/config/test-2oj9-tautomer-pair-rsfe.yaml" configuration = load_config_yaml( config=conf, input_dir="data/", output_dir="data" diff --git a/transformato/tests/test_system_setup.py b/transformato/tests/test_system_setup.py index ce434395..c476d97b 100644 --- a/transformato/tests/test_system_setup.py +++ b/transformato/tests/test_system_setup.py @@ -26,6 +26,7 @@ warnings.filterwarnings("ignore", module="parmed") + def perform_generic_mutation(configuration: dict): s1 = SystemStructure(configuration, "structure1") s2 = SystemStructure(configuration, "structure2") @@ -48,6 +49,7 @@ def perform_generic_mutation(configuration: dict): return i.output_files + def test_read_yaml(): """Sample test, will check ability to read yaml files""" settingsMap = load_config_yaml( @@ -141,7 +143,6 @@ def test_initialize_systems(caplog): @pytest.mark.rsfe def test_setup_system_for_methane_common_core(): - print(get_testsystems_dir()) configuration = load_config_yaml( @@ -191,7 +192,6 @@ def test_setup_system_for_toluene_common_core_with_HMR(): @pytest.mark.rsfe def test_setup_system_for_methane_common_core_with_HMR(): - configuration = load_config_yaml( config=f"{get_testsystems_dir()}/config/test-toluene-methane-rsfe-HMR.yaml", input_dir=get_testsystems_dir(), From 09df630bfc79422d9b489d0c185d148809e8cbaa Mon Sep 17 00:00:00 2001 From: JohannesKarwou Date: Tue, 12 Dec 2023 14:04:03 +0100 Subject: [PATCH 5/7] set exclude Hs to False --- transformato/mutate.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/transformato/mutate.py b/transformato/mutate.py index 62351b97..584e5500 100644 --- a/transformato/mutate.py +++ b/transformato/mutate.py @@ -532,14 +532,14 @@ def _match_terminal_dummy_atoms_between_common_cores( return (lj_default_cc1, lj_default_cc2) + @staticmethod def _calculate_order_of_LJ_mutations( - self, connected_dummy_regions: list, match_terminal_atoms: dict, G: nx.Graph, cyclecheck=True, ordercycles=True, - exclude_Hs=True, + exclude_Hs=False, ) -> list: """ bfs/djikstra-algorithm applied to calculate the ordere for turning of the LJ interactions of the heavy atoms @@ -609,8 +609,8 @@ def _calculate_order_of_LJ_mutations( nodes, cycledict, degreedict, sortedssource, G ) - logger.info("Final mutation route:") - logger.info(nodes) + print("Final mutation route:") + print(nodes) ordered_LJ_mutations.append(nodes) return ordered_LJ_mutations From 2837e3c8f2699f410ab0fc8582d58dc493e24f0f Mon Sep 17 00:00:00 2001 From: JohannesKarwou Date: Tue, 12 Dec 2023 14:56:01 +0100 Subject: [PATCH 6/7] update test --- transformato/tests/test_misc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/transformato/tests/test_misc.py b/transformato/tests/test_misc.py index 25d4a4ec..73ce7a9c 100644 --- a/transformato/tests/test_misc.py +++ b/transformato/tests/test_misc.py @@ -37,7 +37,7 @@ def test_reduced_energy(): def test_convert_to_kT(): kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA - beta = 1.0 / (kB * T) + beta = 1.0 / (kB * temperature) e = -41264.39524669979 * unit.kilojoule_per_mole rE = e * beta From 0ae8a4315542b30045095a8e4e8b5e384bc2500a Mon Sep 17 00:00:00 2001 From: JohannesKarwou Date: Tue, 12 Dec 2023 15:03:50 +0100 Subject: [PATCH 7/7] update CI to newest versions --- .github/workflows/CI.yaml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index c6548204..0de6ef33 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -49,7 +49,7 @@ jobs: # miniforge-variant: Mambaforge-pypy3 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4.1.1 - name: Additional info about the build shell: bash @@ -63,7 +63,7 @@ jobs: - name: Make Cache (no worflow_dispatch) if: ${{ github.event_name != 'workflow_dispatch' }} - uses: actions/cache@v3 + uses: actions/cache@v3.3.2 with: path: ~/conda_pkgs_dir # ${{ matrix.prefix }} # Increase the last number (0) to reset cache manually @@ -72,7 +72,7 @@ jobs: - name: Make Cache (worflow_dispatch) if: ${{ github.event_name == 'workflow_dispatch' }} - uses: actions/cache@v2 + uses: actions/cache@v3.3.2 with: path: ~/conda_pkgs_dir # ${{ matrix.prefix }} key: ${{ matrix.label }}-conda-${{ hashFiles('devtools/conda-envs/fep_env.yaml') }}-${{ env.DATE }}-${{ inputs.CACHE_NUMBER }} @@ -80,7 +80,7 @@ jobs: # More info on options: https://github.com/conda-incubator/setup-miniconda - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3.0.1 with: python-version: ${{ matrix.python-version }} condarc-file: ${{ matrix.condarc-file }} @@ -124,7 +124,7 @@ jobs: pytest -v --cov=transformato --cov-report=xml --color=yes transformato/tests/ - name: CodeCov - uses: codecov/codecov-action@v1 + uses: codecov/codecov-action@v3.1.4 with: file: ./coverage.xml flags: unittests