From 8e69d6e9ea54907458ea6efb5ef5efae8385f58a Mon Sep 17 00:00:00 2001 From: CalCraven Date: Thu, 14 Apr 2022 15:27:03 -0500 Subject: [PATCH 01/32] Atomtype parameterization methods for parameterizing a topology in GMSO. This PR will add the atomtyping module to GMSO for passing a GMSO.Forcefield and a GMSO.Topology and match atomtype using foyer as the backend. Then, the corresponding connection types will be found in the Forcefield and applied to the connections in the topology. * Create parameterize.py, which has the apply function which can take a topology, and a gmso forcefield to apply to it. This can use subgraph isomorphism to identify molecular structures in the topology through the bondgraph and bin those into unique molecules that match a specified forcefield. This apply function can also do the standard atomtyping of the entire topology in one step. * Create isomorph.py which uses networkx graphs to identify disconnected components and isomorphism to identify repeated structures. --- gmso/__init__.py | 1 + gmso/atomtyping/__init__.py | 1 + gmso/atomtyping/isomorph.py | 59 +++++++ gmso/atomtyping/parameterize.py | 275 ++++++++++++++++++++++++++++++++ 4 files changed, 336 insertions(+) create mode 100644 gmso/atomtyping/__init__.py create mode 100644 gmso/atomtyping/isomorph.py create mode 100644 gmso/atomtyping/parameterize.py diff --git a/gmso/__init__.py b/gmso/__init__.py index 9a5985491..2e9112396 100644 --- a/gmso/__init__.py +++ b/gmso/__init__.py @@ -1,4 +1,5 @@ """GMSO: General Molecular Simulation Object.""" +from .atomtyping.parameterize import apply from .core.angle import Angle from .core.angle_type import AngleType from .core.atom import Atom diff --git a/gmso/atomtyping/__init__.py b/gmso/atomtyping/__init__.py new file mode 100644 index 000000000..729a6c586 --- /dev/null +++ b/gmso/atomtyping/__init__.py @@ -0,0 +1 @@ +"""GMSO functions that generate parameterized Topologies.""" diff --git a/gmso/atomtyping/isomorph.py b/gmso/atomtyping/isomorph.py new file mode 100644 index 000000000..27ba74d8d --- /dev/null +++ b/gmso/atomtyping/isomorph.py @@ -0,0 +1,59 @@ +"""TopologyGraph Functions that identify molecules from isomorphism.""" +from collections import deque + +import networkx as nx + + +def partition_isomorphic_topology_graphs(graph): + """Return a dictionary of subgraphs that are partitioned by isomorphism. + + Parameters + ---------- + graph: foyer.topology_graph.TopologyGraph + The networkx subclassed TopologyGraph with data identifying the nodes + and edges that make up a topology atom and bonds structure + + Returns + ------- + isomorphic_elements: dict + The keys are unique disconnected graphs, and the values are all + identical subgraphs in the graph + + Notes + ----- + See https://github.com/networkx/networkx/blob/main/networkx/algorithms/isomorphism/isomorphvf2.py + from the networkx documentation about identifying isomorphic components + """ + + subgraphs_gen = (graph.subgraph(c) for c in nx.connected_components(graph)) + subgraphs_list = list(subgraphs_gen) + + graph_queue = deque(subgraphs_list) + graph_of_interest = graph_queue.popleft() + isomorphic_elements = { + graph_of_interest: [], + } + + last_element_popped = None + count = 0 + first_mismatch = None + while graph_queue: + if graph_queue[0] == first_mismatch: + count = 0 + graph_of_interest = graph_queue.popleft() + isomorphic_elements[graph_of_interest] = [] + if graph_queue: + graph = graph_queue.popleft() + matcher = nx.algorithms.isomorphism.GraphMatcher( + graph, graph_of_interest + ) + if matcher.is_isomorphic(): + isomorphic_elements[graph_of_interest].append( + (graph, matcher.mapping) + ) + else: + if count == 0: + first_mismatch = graph + graph_queue.append(graph) + count += 1 + return isomorphic_elements diff --git a/gmso/atomtyping/parameterize.py b/gmso/atomtyping/parameterize.py new file mode 100644 index 000000000..0dace83bd --- /dev/null +++ b/gmso/atomtyping/parameterize.py @@ -0,0 +1,275 @@ +"""Functions used to atomtype a gmso.Topology.""" +from operator import attrgetter + +import networkx as nx +from foyer.atomtyper import AtomTypingRulesProvider, find_atomtypes +from foyer.smarts import SMARTS +from foyer.topology_graph import TopologyGraph + +from gmso import Atom +from gmso.atomtyping.isomorph import ( # check refs + partition_isomorphic_topology_graphs, +) + +__all__ = ["apply"] + + +def apply( + top, + forcefields, + identify_connected_components=True, + use_residue_info=False, + assert_bond_params=True, + assert_angle_params=True, + assert_dihedral_params=True, + assert_improper_params=False, +): + """Set Topology parameter types from GMSO Forcefields. + + Parameters + ---------- + top: gmso.core.topology.Topology, required + The GMSO topology on which to apply forcefields + forcefields: dict, required + The keys are labels that match the subtopology label or site residue_name, and the + values are gmso Forcefield objects that gets applied to the specified molecule + identify_connected_components: bool, optional, default=True + A flag to determine whether or not to search the topology for repeated disconnected + structures, otherwise known as molecules and type each molecule only once. + use_residue_info: bool, optional, default=False + A flag to determine whether or not to look at site.residue_name to look parameterize + each molecule only once. Will only be used if identify_connected_components=False + assert_bond_params : bool, optional, default=True + If True, Foyer will exit if parameters are not found for all system + bonds. + assert_angle_params : bool, optional, default=True + If True, Foyer will exit if parameters are not found for all system + angles. + assert_dihedral_params : bool, optional, default=True + If True, Foyer will exit if parameters are not found for all system + proper dihedrals. + assert_improper_params : bool, optional, default=False + If True, Foyer will exit if parameters are not found for all system + improper dihedrals. + """ + top.identify_connections() + top_graph = get_modified_topology_graph(top) + if identify_connected_components: # True, (True,False) + # populate molecule information + isomorphic_top_graphs = partition_isomorphic_topology_graphs(top_graph) + apply_atomtypes_from_isomorphic( + top, isomorphic_top_graphs, forcefields + ) # what if only one forcefield?, traversing done in this function + elif not use_residue_info: # False, False + # get entire topgraph + topgraph = get_modified_topology_graph(top) + # get rules provider from gmso_ff + # Assumes forcefields is just one forcefield + at_rules_provider = get_atomtyping_rules_provider(gmso_ff=forcefields) + # create a typemap of that topgraph + typemap = find_atomtypes(topgraph, at_rules_provider) + # apply typemap + traverse_typemap(top, typemap, forcefields) + elif use_residue_map: # False, True + """Not Implemented + maps = {} + for key, val in forcefields.items(): + # generate topgraph of items with "key" residue name + # NOTE CAL: this iter_sites_by method won't work currently, can only iter one at a time. + res_iter = top.iter_sites_by(["residue_name", "residue_number"], [key, 1]) #Need to note that resnumber index 1 + # create a top/graph from this set of sites + #subgraph = FUNCTION(res_iter) + # append that typemap to a dict with residue name as key + #maps[key] = find_atomtypes(subgraph, val) + # generate a full typemap by iterating through all sites and copying over map info from `maps` + #typemap = FUNCTION(top, maps) + # apply typemap + traverse_typemap(top, typemap, forcefields) + """ + raise ( + GMSOError, + "Using site residue matching to match substructures to a given forcefield is not implemented yet", + ) + + if assert_bond_params: + apply_connection_types(top, forcefields, "bond") + if assert_angle_params: + apply_connection_types(top, forcefields, "angle") + if assert_dihedral_params: + apply_connection_types(top, forcefields, "dihedral") + if assert_improper_params: + apply_connection_types(top, forcefields, "improper") + + return top + + +def get_modified_topology_graph(gmso_topology): + """Return a TopologyGraph with relevant attributes from an GMSO topology. + + Parameters + ---------- + gmso_topology: gmso.Topology + The GMSO Topology + + Returns + ------- + TopologyGraph + The equivalent TopologyGraph of the openFF Topology `openff_topology` + """ + top_graph = TopologyGraph() + for atom in gmso_topology.sites: + if isinstance(atom, Atom): + if atom.name.startswith("_"): + top_graph.add_atom( + name=atom.name, + index=gmso_topology.get_index(atom), + atomic_number=None, + element=atom.name, + subtopology_label=atom.label, + ) + + else: + top_graph.add_atom( + name=atom.name, + index=gmso_topology.get_index(atom), + atomic_number=atom.element.atomic_number, + element=atom.element.symbol, + subtopology_label=atom.label, + ) + + for top_bond in gmso_topology.bonds: + atoms_indices = [ + gmso_topology.get_index(atom) + for atom in top_bond.connection_members + ] + top_graph.add_bond(atoms_indices[0], atoms_indices[1]) + + return top_graph + + +def get_atomtyping_rules_provider(gmso_ff): + """Return a foyer AtomTypingRulesProvider from a GMSO forcefield. + + Parameters + ---------- + gmso_ff: gmso.core.forcefield.Forcefield + The GMSO forcefield object to extract the rules from + Returns + ------- + AtomTypingRulesProvider + The foyer.atomtyper.AtomTypingRulesProvider object used to parse atomtype definitions. + Typically, SMARTS is the ruleset of choice. See https://github.com/mosdef-hub/foyer/issues/63 + for curently supported features in Foyer. + """ + atom_type_defs = {} + atom_type_overrides = {} + parser = SMARTS({}) + for atom_type_name, atom_type in gmso_ff.atom_types.items(): + if atom_type.definition: + atom_type_defs[atom_type_name] = atom_type.definition + if atom_type.overrides: + atom_type_overrides[atom_type_name] = atom_type.overrides + + return AtomTypingRulesProvider( + atom_type_defs, atom_type_overrides, {}, parser + ) + + +def apply_atomtypes_from_isomorphic(top, isomorphic_top_graphs, forcefields): + """Set Topology atomtypes from isomorphic structures matching supplied forcefields. + + Parameters + ---------- + top: gmso.core.topology.Topology, required + The GMSO topology on which to apply atomtypes + isomorphic_top_graphs: dict, required + The dictionary mapping which breaks up molecules into one graph structure, and a list of + the corresponding subgraphs in the TopologyGraph that are isomorphic + forcefields: dict, required + The keys are labels that match the subtopology labels for each isomorphic graph, and the + values are gmso Forcefield objects that gets applied to the specified molecule + """ + for top_graph, identical in isomorphic_top_graphs.items(): + sub_top_label = next( + iter( + nx.classes.function.get_node_attributes( + top_graph, "atom_data" + ).values() + ) + ).subtopology_label + ff_of_interest = forcefields[sub_top_label] + at_rules_provider = get_atomtyping_rules_provider( + gmso_ff=ff_of_interest + ) + sub_type_map = find_atomtypes(top_graph, at_rules_provider) + traverse_typemap(top, sub_type_map, ff_of_interest) + for graph, mapping in identical: + for node in graph: + mapped = sub_type_map[mapping[node]]["atom_type"] + top._sites[node].atom_type = mapped.clone() + + +def traverse_typemap(top, type_map, forcefield): + """Take a typemap and applies those atomtypes to the Topology. + + Parameters + ---------- + top: gmso.core.topology.Topology, required + The GMSO topology on which to apply atomtypes + typemap: dict, required + The dictionary of atom_index with iterated matching to get the atomtype of that site + forcefield: gmso.forcefield.Forcefield, required + The GMSO forcefield object which has information associated with each unique atomtype + """ + for index in type_map: + site_of_interest = top._sites[index] + site_of_interest.atom_type = forcefield.atom_types[ + type_map[index]["atomtype"] + ].clone() + type_map[index]["atom_type"] = site_of_interest.atom_type + + +def apply_residue_names(top, more): # TODO + """Take a topology and applies residue info to all sites. + + Parameters + ---------- + top: gmso.core.topology.Topology, required + The GMSO topology on which to apply atomtypes + more: TBD, required + The information necessary to map residue name and number for each site + """ + return + + +def apply_connection_types(top, forcefields, parameter): + """Set Topology connection types from matching atomclass information in forcefield. + + Parameters + ---------- + top: gmso.core.topology.Topology, required + The GMSO topology on which to apply paramter_types + forcefields: dict, required + The keys are labels that match the subtopology labels for each molecule, and the + values are gmso Forcefield objects that gets applied to the specified molecule + parameter: str, required + The connection type to parameterize. Can any of "bond", "angle", "dihedral", "improper" + """ + connect_members = { + "bond": [0, 1], + "angle": [0, 1, 2], + "dihedral": [0, 1, 2, 3], + "improper": [0, 1, 2, 3], + } # TODO validate improper order + for connect in getattr(top, "_" + parameter + "s"): + sub_top_label = connect.connection_members[0].label_ + ff_of_interest = forcefields[sub_top_label] + type_getter = attrgetter("atom_type.atomclass") + member_types = [ + type_getter(getattr(connect, "connection_members")[i]) + for i in connect_members[parameter] + ] + ctype = ff_of_interest.get_potential( + group=parameter + "_type", key=member_types + ).clone() + setattr(connect, parameter + "_type", ctype) From 9074909c4df6dffec53ddbfcd2d6ac8b2aadb56e Mon Sep 17 00:00:00 2001 From: CalCraven Date: Thu, 14 Apr 2022 16:52:32 -0500 Subject: [PATCH 02/32] Move module imports for apply into atomtyping to prevent circular imports --- gmso/__init__.py | 1 - gmso/atomtyping/__init__.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/gmso/__init__.py b/gmso/__init__.py index 2e9112396..9a5985491 100644 --- a/gmso/__init__.py +++ b/gmso/__init__.py @@ -1,5 +1,4 @@ """GMSO: General Molecular Simulation Object.""" -from .atomtyping.parameterize import apply from .core.angle import Angle from .core.angle_type import AngleType from .core.atom import Atom diff --git a/gmso/atomtyping/__init__.py b/gmso/atomtyping/__init__.py index 729a6c586..303bb1ca7 100644 --- a/gmso/atomtyping/__init__.py +++ b/gmso/atomtyping/__init__.py @@ -1 +1,2 @@ """GMSO functions that generate parameterized Topologies.""" +from .parameterize import apply From 80c1850dfb16b992953fd97a525b5d67659bda36 Mon Sep 17 00:00:00 2001 From: CalCraven Date: Sun, 24 Apr 2022 23:25:40 -0500 Subject: [PATCH 03/32] Add a quick fix which will do atomtyping if no residue flag is passed --- gmso/atomtyping/parameterize.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/gmso/atomtyping/parameterize.py b/gmso/atomtyping/parameterize.py index 0dace83bd..fd722961e 100644 --- a/gmso/atomtyping/parameterize.py +++ b/gmso/atomtyping/parameterize.py @@ -65,11 +65,12 @@ def apply( topgraph = get_modified_topology_graph(top) # get rules provider from gmso_ff # Assumes forcefields is just one forcefield - at_rules_provider = get_atomtyping_rules_provider(gmso_ff=forcefields) + ff_of_interest = next(iter(forcefields.values())) + at_rules_provider = get_atomtyping_rules_provider(gmso_ff=ff_of_interest) # create a typemap of that topgraph typemap = find_atomtypes(topgraph, at_rules_provider) # apply typemap - traverse_typemap(top, typemap, forcefields) + traverse_typemap(top, typemap, ff_of_interest) elif use_residue_map: # False, True """Not Implemented maps = {} From 6376da79525264ceb2e5aeae5546baf18bf0ea32 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 25 Apr 2022 04:26:06 +0000 Subject: [PATCH 04/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gmso/atomtyping/parameterize.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gmso/atomtyping/parameterize.py b/gmso/atomtyping/parameterize.py index fd722961e..4b9c3ba6a 100644 --- a/gmso/atomtyping/parameterize.py +++ b/gmso/atomtyping/parameterize.py @@ -66,7 +66,9 @@ def apply( # get rules provider from gmso_ff # Assumes forcefields is just one forcefield ff_of_interest = next(iter(forcefields.values())) - at_rules_provider = get_atomtyping_rules_provider(gmso_ff=ff_of_interest) + at_rules_provider = get_atomtyping_rules_provider( + gmso_ff=ff_of_interest + ) # create a typemap of that topgraph typemap = find_atomtypes(topgraph, at_rules_provider) # apply typemap From 99c0793b07a3896659414eabd7a381d151e54751 Mon Sep 17 00:00:00 2001 From: Co Quach Date: Wed, 27 Apr 2022 00:38:07 +0700 Subject: [PATCH 05/32] correctly update sites when adding subtop --- gmso/core/topology.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gmso/core/topology.py b/gmso/core/topology.py index 12fec0f81..ff2f806bd 100644 --- a/gmso/core/topology.py +++ b/gmso/core/topology.py @@ -624,7 +624,7 @@ def add_subtopology(self, subtop, update=True): """ self._subtops.add(subtop) subtop.parent = self - self._sites.union(subtop.sites) + self._sites = self._sites.union(subtop.sites) if update: self.update_topology() From e1f3d4e0ae94b17eddf92395f7e6e2920f17c9df Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 2 May 2022 10:33:29 -0500 Subject: [PATCH 06/32] Changes to doc strings for clarity. Add a subtop_label variable to generalize where the molecule definition is pulled from --- gmso/atomtyping/isomorph.py | 2 +- gmso/atomtyping/parameterize.py | 20 +++++++++++++------- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/gmso/atomtyping/isomorph.py b/gmso/atomtyping/isomorph.py index 27ba74d8d..49b0ce736 100644 --- a/gmso/atomtyping/isomorph.py +++ b/gmso/atomtyping/isomorph.py @@ -5,7 +5,7 @@ def partition_isomorphic_topology_graphs(graph): - """Return a dictionary of subgraphs that are partitioned by isomorphism. + """Return a collection of isomorphic sets of the subgraphs of the Topology Graph.. Parameters ---------- diff --git a/gmso/atomtyping/parameterize.py b/gmso/atomtyping/parameterize.py index 4b9c3ba6a..2785d0fb7 100644 --- a/gmso/atomtyping/parameterize.py +++ b/gmso/atomtyping/parameterize.py @@ -40,16 +40,16 @@ def apply( A flag to determine whether or not to look at site.residue_name to look parameterize each molecule only once. Will only be used if identify_connected_components=False assert_bond_params : bool, optional, default=True - If True, Foyer will exit if parameters are not found for all system + If True, an error is raised if parameters are not found for all system bonds. assert_angle_params : bool, optional, default=True - If True, Foyer will exit if parameters are not found for all system + If True, an error is raised if parameters are not found for all system angles. assert_dihedral_params : bool, optional, default=True - If True, Foyer will exit if parameters are not found for all system + If True, an error is raised if parameters are not found for all system proper dihedrals. assert_improper_params : bool, optional, default=False - If True, Foyer will exit if parameters are not found for all system + If True, an error is raised if parameters are not found for all system improper dihedrals. """ top.identify_connections() @@ -128,7 +128,7 @@ def get_modified_topology_graph(gmso_topology): index=gmso_topology.get_index(atom), atomic_number=None, element=atom.name, - subtopology_label=atom.label, + subtopology_label=atom.label, #TODO: Change based on PR #638 conventions ) else: @@ -137,7 +137,7 @@ def get_modified_topology_graph(gmso_topology): index=gmso_topology.get_index(atom), atomic_number=atom.element.atomic_number, element=atom.element.symbol, - subtopology_label=atom.label, + subtopology_label=atom.label, #TODO: Change based on PR #638 conventions ) for top_bond in gmso_topology.bonds: @@ -264,8 +264,14 @@ def apply_connection_types(top, forcefields, parameter): "dihedral": [0, 1, 2, 3], "improper": [0, 1, 2, 3], } # TODO validate improper order + molecule_label = "label_" # place to grab info referencing molecule name for connect in getattr(top, "_" + parameter + "s"): - sub_top_label = connect.connection_members[0].label_ + sub_top_label = connect.connection_members[0].get(molecule_label) + if sub_top_label is None: + msg = f"Site {connect.connection_members[0]}in this connection has no molecule_label {molecule_label}.\ + This string should correspond to one of the keys associated with the forcefield that has been \ + passed. Currently provided name:GMSOForcefields pairs are {forcefields.items()}" + raise ParameterizationError(msg) ff_of_interest = forcefields[sub_top_label] type_getter = attrgetter("atom_type.atomclass") member_types = [ From 15035781e1dd7ee73e7187de1acd6ae6353b8f9c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 2 May 2022 15:33:54 +0000 Subject: [PATCH 07/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gmso/atomtyping/parameterize.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/gmso/atomtyping/parameterize.py b/gmso/atomtyping/parameterize.py index 2785d0fb7..a0adf0864 100644 --- a/gmso/atomtyping/parameterize.py +++ b/gmso/atomtyping/parameterize.py @@ -128,7 +128,7 @@ def get_modified_topology_graph(gmso_topology): index=gmso_topology.get_index(atom), atomic_number=None, element=atom.name, - subtopology_label=atom.label, #TODO: Change based on PR #638 conventions + subtopology_label=atom.label, # TODO: Change based on PR #638 conventions ) else: @@ -137,7 +137,7 @@ def get_modified_topology_graph(gmso_topology): index=gmso_topology.get_index(atom), atomic_number=atom.element.atomic_number, element=atom.element.symbol, - subtopology_label=atom.label, #TODO: Change based on PR #638 conventions + subtopology_label=atom.label, # TODO: Change based on PR #638 conventions ) for top_bond in gmso_topology.bonds: @@ -264,7 +264,7 @@ def apply_connection_types(top, forcefields, parameter): "dihedral": [0, 1, 2, 3], "improper": [0, 1, 2, 3], } # TODO validate improper order - molecule_label = "label_" # place to grab info referencing molecule name + molecule_label = "label_" # place to grab info referencing molecule name for connect in getattr(top, "_" + parameter + "s"): sub_top_label = connect.connection_members[0].get(molecule_label) if sub_top_label is None: From 7be5f5739e740197887eecdaf235e016c63c494f Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 6 May 2022 14:11:18 -0500 Subject: [PATCH 08/32] More modular architecture --- .../__init__.py | 0 gmso/parameterization/foyer_utils.py | 111 ++++++++ .../isomorph.py | 38 +-- .../parameterize.py | 2 +- gmso/parameterization/parameterizer.py | 241 ++++++++++++++++++ gmso/parameterization/subtopology_utils.py | 50 ++++ gmso/parameterization/utils.py | 15 ++ 7 files changed, 439 insertions(+), 18 deletions(-) rename gmso/{atomtyping => parameterization}/__init__.py (100%) create mode 100644 gmso/parameterization/foyer_utils.py rename gmso/{atomtyping => parameterization}/isomorph.py (63%) rename gmso/{atomtyping => parameterization}/parameterize.py (99%) create mode 100644 gmso/parameterization/parameterizer.py create mode 100644 gmso/parameterization/subtopology_utils.py create mode 100644 gmso/parameterization/utils.py diff --git a/gmso/atomtyping/__init__.py b/gmso/parameterization/__init__.py similarity index 100% rename from gmso/atomtyping/__init__.py rename to gmso/parameterization/__init__.py diff --git a/gmso/parameterization/foyer_utils.py b/gmso/parameterization/foyer_utils.py new file mode 100644 index 000000000..bc1b695bb --- /dev/null +++ b/gmso/parameterization/foyer_utils.py @@ -0,0 +1,111 @@ +"""Utilities for atomtyping a gmso topology with foyer.""" +from collections import namedtuple + +from foyer.atomtyper import AtomTypingRulesProvider, find_atomtypes +from foyer.topology_graph import TopologyGraph + +from gmso.core.atom import Atom +from gmso.parameterization.subtopology_utils import subtop_bonds +from gmso.utils.decorators import experimental + + +@experimental +def get_topology_graph(gmso_topology, atomdata_populator=None): + """Return a TopologyGraph with relevant attributes from an GMSO topology. + + Parameters + ---------- + gmso_topology: gmso.Topology-like + The GMSO Topology + + atomdata_populator: callable, default=None + A function that will be called with the following arguments `gmso_topology` as well as `atom` to pass extra + arguments to the foyer.topology_graph.AtomData object + + Notes + ----- + The gmso topology here is a duck type. + + Returns + ------- + foyer.topology_graph.TopologyGraph + A light networkx representation of the topology + """ + top_graph = TopologyGraph() + for atom in gmso_topology.sites: + if isinstance(atom, Atom): + kwargs = ( + atomdata_populator(gmso_topology, atom) + if atomdata_populator + else {} + ) + if atom.name.startswith("_"): + top_graph.add_atom( + name=atom.name, + index=gmso_topology.get_index(atom), + atomic_number=None, + element=atom.name, + **kwargs, + ) + + else: + top_graph.add_atom( + name=atom.name, + index=gmso_topology.get_index(atom), + atomic_number=atom.element.atomic_number, + element=atom.element.symbol, + **kwargs, + ) + + for top_bond in gmso_topology.bonds: + atoms_indices = [ + gmso_topology.get_index(atom) + for atom in top_bond.connection_members + ] + top_graph.add_bond(atoms_indices[0], atoms_indices[1]) + + return top_graph + + +@experimental +def get_topology_graph_from_subtop(subtopology): + """Get an equivalent topology graph for a sub-topology.""" + subtop_named_tuple = namedtuple("subtopology", ("sites", "bonds")) + return get_topology_graph( + subtop_named_tuple(subtopology.sites, subtop_bonds(subtopology)) + ) + + +@experimental +def get_atomtyping_rules_provider(gmso_ff): + """Return a foyer AtomTypingRulesProvider from a GMSO forcefield. + + Parameters + ---------- + gmso_ff: gmso.core.forcefield.Forcefield + The GMSO forcefield object to extract the rules from + + Returns + ------- + AtomTypingRulesProvider + The foyer.atomtyper.AtomTypingRulesProvider object used to parse atomtype definitions. + Typically, SMARTS is the ruleset of choice. See https://github.com/mosdef-hub/foyer/issues/63 + for curently supported features in Foyer. + """ + atom_type_defs = {} + atom_type_overrides = {} + for atom_type_name, atom_type in gmso_ff.atom_types.items(): + if atom_type.definition: + atom_type_defs[atom_type_name] = atom_type.definition + if atom_type.overrides: + atom_type_overrides[atom_type_name] = atom_type.overrides + + return AtomTypingRulesProvider( + atom_type_defs, atom_type_overrides, gmso_ff.non_element_types + ) + + +@experimental +def typemap_dict(topology_graph, atomtyping_rules_provider, max_iter=10): + """Return a dictionary of typemap, by finding atomtypes in foyer.""" + return find_atomtypes(topology_graph, atomtyping_rules_provider, max_iter) diff --git a/gmso/atomtyping/isomorph.py b/gmso/parameterization/isomorph.py similarity index 63% rename from gmso/atomtyping/isomorph.py rename to gmso/parameterization/isomorph.py index 49b0ce736..811cc429a 100644 --- a/gmso/atomtyping/isomorph.py +++ b/gmso/parameterization/isomorph.py @@ -4,39 +4,43 @@ import networkx as nx +def top_node_match(n1, n2): + """Match two nodes in the topology graph based on their elements.""" + return n1["atom_data"].element == n2["atom_data"].element + + def partition_isomorphic_topology_graphs(graph): - """Return a collection of isomorphic sets of the subgraphs of the Topology Graph.. + """Return a collection of isomorphic sets of the subgraphs of the Topology Graph. - Parameters - ---------- - graph: foyer.topology_graph.TopologyGraph - The networkx subclassed TopologyGraph with data identifying the nodes - and edges that make up a topology atom and bonds structure + Parameters + ---------- + graph: foyer.topology_graph.TopologyGraph + The networkx subclassed TopologyGraph with data identifying the nodes + and edges that make up a topology atom and bonds structure - Returns - ------- - isomorphic_elements: dict - The keys are unique disconnected graphs, and the values are all - identical subgraphs in the graph + Returns + ------- + isomorphic_elements: dict + The keys are unique disconnected graphs, and the values are all + identical subgraphs in the graph Notes ----- See https://github.com/networkx/networkx/blob/main/networkx/algorithms/isomorphism/isomorphvf2.py from the networkx documentation about identifying isomorphic components """ + graph_queue = deque( + graph.subgraph(c) for c in nx.connected_components(graph) + ) - subgraphs_gen = (graph.subgraph(c) for c in nx.connected_components(graph)) - subgraphs_list = list(subgraphs_gen) - - graph_queue = deque(subgraphs_list) graph_of_interest = graph_queue.popleft() isomorphic_elements = { graph_of_interest: [], } - last_element_popped = None count = 0 first_mismatch = None + while graph_queue: if graph_queue[0] == first_mismatch: count = 0 @@ -45,7 +49,7 @@ def partition_isomorphic_topology_graphs(graph): if graph_queue: graph = graph_queue.popleft() matcher = nx.algorithms.isomorphism.GraphMatcher( - graph, graph_of_interest + graph, graph_of_interest, node_match=top_node_match ) if matcher.is_isomorphic(): isomorphic_elements[graph_of_interest].append( diff --git a/gmso/atomtyping/parameterize.py b/gmso/parameterization/parameterize.py similarity index 99% rename from gmso/atomtyping/parameterize.py rename to gmso/parameterization/parameterize.py index a0adf0864..08e0258e4 100644 --- a/gmso/atomtyping/parameterize.py +++ b/gmso/parameterization/parameterize.py @@ -7,7 +7,7 @@ from foyer.topology_graph import TopologyGraph from gmso import Atom -from gmso.atomtyping.isomorph import ( # check refs +from gmso.parameterization.isomorph import ( # check refs partition_isomorphic_topology_graphs, ) diff --git a/gmso/parameterization/parameterizer.py b/gmso/parameterization/parameterizer.py new file mode 100644 index 000000000..30d09d7ea --- /dev/null +++ b/gmso/parameterization/parameterizer.py @@ -0,0 +1,241 @@ +"""The parameterizer module for a gmso Topology.""" + +import warnings +from typing import Dict, Union + +from pydantic import Field + +from gmso.abc.gmso_base import GMSOBase +from gmso.core.forcefield import ForceField +from gmso.core.topology import Topology +from gmso.exceptions import GMSOError +from gmso.parameterization.foyer_utils import ( + get_atomtyping_rules_provider, + get_topology_graph, + get_topology_graph_from_subtop, + typemap_dict, +) +from gmso.parameterization.isomorph import partition_isomorphic_topology_graphs +from gmso.parameterization.subtopology_utils import ( + assert_no_boundary_bonds, + subtop_angles, + subtop_bonds, + subtop_dihedrals, + subtop_impropers, +) +from gmso.parameterization.utils import POTENTIAL_GROUPS + + +class GMSOParameterizationError(GMSOError): + """Raise when parameterization fails.""" + + +class Parameterizer(GMSOBase): + """Utility class to parameterize a topology with gmso Forcefield.""" + + topology: Topology = Field(..., description="The gmso topology.") + + forcefields: Union[ForceField, Dict[str, ForceField]] = Field( + ..., + description="The gmso forcefield/ a dictionary of gmso " + "forcefields per sub-topology, where the keys " + "should match the subtopology names", + ) + + def get_ff(self, key=None): + """Return the forcefield of choice by looking up the forcefield dictionary.""" + if isinstance(self.forcefields, Dict): + return self.forcefields.get(key) + else: + return self.forcefields + + def _parameterize_sites(self, top_or_subtop, typemap, ff): + """Parameterize sites with appropriate atom-types from the forcefield.""" + for j, site in enumerate(top_or_subtop.sites): + site.atom_type = ff.get_potential( + "atom_type", typemap[j]["atom_type"] + ).clone() # Always properly indexed or not? + + def _parameterize_connections( + self, + top_or_subtop, + ff, + is_subtop=False, + assert_bond_params=True, + assert_angle_params=True, + assert_dihedral_params=True, + assert_improper_params=False, + ): + """Parameterize connections with appropriate potentials from the forcefield.""" + if is_subtop: + bonds = subtop_bonds(top_or_subtop) + angles = subtop_angles(top_or_subtop) + dihedrals = subtop_dihedrals(top_or_subtop) + impropers = subtop_impropers(top_or_subtop) + else: + bonds = top_or_subtop.bonds + angles = top_or_subtop.angles + dihedrals = top_or_subtop.dihedrals + impropers = top_or_subtop.impropers + + self._apply_connection_parameters(bonds, ff, assert_bond_params) + self._apply_connection_parameters(angles, ff, assert_angle_params) + self._apply_connection_parameters(dihedrals, ff, assert_dihedral_params) + self._apply_connection_parameters(impropers, ff, assert_improper_params) + + def _apply_connection_parameters( + self, connections, ff, error_on_missing=True + ): + """Find and assign potentials from the forcefield for the provided connections.""" + visited = dict() + + for connection in connections: + group, connection_identifiers = self.connection_identifier( + connection + ) + match = None + for identifier_key in connection_identifiers: + if tuple(identifier_key) in visited: + match = visited[tuple(identifier_key)] + break + + match = ff.get_potential( + group=group, key=identifier_key, warn=True + ) + if match: + visited[tuple[identifier_key]] = match + break + + if not match and error_on_missing: + raise GMSOParameterizationError( + f"No parameters found for connection {connection} in the Forcefield." + ) + setattr(connection, group, match.clone()) + + def _parameterize( + self, + subtop_or_top, + typemap, + is_subtop=False, + assert_bond_params=True, + assert_angle_params=True, + assert_dihedral_params=True, + assert_improper_params=False, + ): + """Parameterize a topology/subtopology based on an atomtype map.""" + forcefield = self.get_ff(subtop_or_top.name) + self._parameterize_sites(subtop_or_top.sites, typemap, forcefield) + self._parameterize_connections( + subtop_or_top, + forcefield, + is_subtop=is_subtop, + assert_bond_params=assert_bond_params, + assert_angle_params=assert_angle_params, + assert_dihedral_params=assert_dihedral_params, + assert_improper_params=assert_improper_params, + ) + + def apply( + self, + identify_connected_components=True, + use_residue_info=False, + assert_bond_params=True, + assert_angle_params=True, + assert_dihedral_params=True, + assert_improper_params=False, + ): + """Apply the current forcefield(s) to the topology/ various subtopologies.""" + parameterize_kwargs = dict( + assert_bond_params=assert_bond_params, + assert_angle_params=assert_angle_params, + assert_dihedral_params=assert_dihedral_params, + assert_improper_params=assert_improper_params, + ) + if self.topology.is_typed(): + raise GMSOParameterizationError( + "Cannot parameterize a typed topology. Please provide a topology without any types" + ) + + if isinstance(self.forcefields, Dict): + for subtop in self.topology.subtops: + if subtop.name not in self.forcefields: + warnings.warn( + f"Subtopology {subtop.name} will not be parameterized, as the forcefield to parameterize it " + f"is missing." + ) # FixMe: Will warning be enough? + else: + assert_no_boundary_bonds(subtop) + typemap = self._get_atomtypes( + self.get_ff(subtop.name), + subtop, + identify_connected_components, + is_subtop=True, + ) + self._parameterize( + typemap, + subtop, + is_subtop=True, # This will be removed from the future iterations + **parameterize_kwargs, + ) + else: + typemap = self._get_atomtypes( + self.get_ff(), + self.topology, + identify_connected_components, + is_subtop=False, + ) + self._parameterize( + self.topology, + typemap, + is_subtop=False, # This will be removed from the future iterations + **parameterize_kwargs, + ) + + @staticmethod + def connection_identifier( + connection, + ): # This can extended to incorporate a pluggable object from the forcefield. + """Return the group and list of identifiers for a connection to query the forcefield for its potential.""" + group = POTENTIAL_GROUPS[type(connection)] + return group, [ + list( + member.atom_type.atomclass + for member in connection.connection_members + ), + list( + member.atom_type.name + for member in connection.connection_members + ), + ] + + @staticmethod + def _get_atomtypes( + forcefield, topology, use_isomprohic_checks=False, is_subtop=False + ): + """Run atom-typing in foyer and return a typemap.""" + atom_typing_rules_provider = get_atomtyping_rules_provider(forcefield) + + if is_subtop: + foyer_topology_graph = get_topology_graph_from_subtop(topology) + else: + foyer_topology_graph = get_topology_graph(topology) + + if use_isomprohic_checks: + isomorphic_substructures = partition_isomorphic_topology_graphs( + foyer_topology_graph + ) + for graph, mirrors in isomorphic_substructures.items(): + typemap = typemap_dict( + atomtyping_rules_provider=atom_typing_rules_provider, + topology_graph=graph, + ) + for mirror, mapping in mirrors: + for node in mirror: + typemap[node] = typemap[mapping[node]] + return typemap + + else: + return typemap_dict( + topology_graph=foyer_topology_graph, + atomtyping_rules_provider=atom_typing_rules_provider, + ) diff --git a/gmso/parameterization/subtopology_utils.py b/gmso/parameterization/subtopology_utils.py new file mode 100644 index 000000000..de32564b9 --- /dev/null +++ b/gmso/parameterization/subtopology_utils.py @@ -0,0 +1,50 @@ +"""Utilities for application of a particular forcefield to a subtopology.""" + + +def _members_in_subtop(connection, subtop): + """Check if all the members in a connection belong to a subtopology.""" + return all(site in subtop.sites for site in connection.connection_members) + + +def _subtop_connections(subtop, attr): + """Return all the connections belonging to a subtopology.""" + return filter( + lambda conn: _members_in_subtop(conn, subtop), + getattr(subtop._parent, attr), + ) + + +def subtop_bonds(subtop): + """Given a subtopology, return its bonds.""" + return _subtop_connections(subtop, "bonds") + + +def subtop_angles(subtop): + """Given a subtopology, return its angles.""" + return _subtop_connections(subtop, "angles") + + +def subtop_dihedrals(subtop): + """Given a subtopology, return its dihedrals.""" + return _subtop_connections(subtop, "dihedrals") + + +def subtop_impropers(subtop): + """Given a subtopology, return its impropers.""" + return _subtop_connections(subtop, "impropers") + + +def assert_no_boundary_bonds(subtop): + """Given a subtopology, assert that no bonds exist between its sites and external sites.""" + for bond in subtop._parent.bonds: + site_pairs = bond.connection_members + assertion_msg = "Site {} is in the subtopology {}, but its bonded partner {} is not." + + if site_pairs[0] in subtop.sites: + assert site_pairs[1] in subtop.sites, assertion_msg.format( + site_pairs[0].name, subtop.name, site_pairs[1].name + ) + elif site_pairs[1] in subtop.sites: + assert site_pairs[0] in subtop.sites, assertion_msg.format( + site_pairs[1].name, subtop.name, site_pairs[0].name + ) diff --git a/gmso/parameterization/utils.py b/gmso/parameterization/utils.py new file mode 100644 index 000000000..3fb94404e --- /dev/null +++ b/gmso/parameterization/utils.py @@ -0,0 +1,15 @@ +"""Generic utilities for parameterizing a gmso Topology.""" + +from gmso.core.angle import Angle +from gmso.core.atom import Atom +from gmso.core.bond import Bond +from gmso.core.dihedral import Dihedral +from gmso.core.improper import Improper + +POTENTIAL_GROUPS = { + Bond: "bond_type", + Angle: "angle_type", + Dihedral: "dihedral_type", + Improper: "improper_type", + Atom: "atom_type", +} From 73664522832a73bb7ac754cfe76f8c337a636835 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 6 May 2022 16:55:38 -0500 Subject: [PATCH 09/32] WIP- Add testing, minor refactoring of the APIs. --- gmso/parameterization/foyer_utils.py | 5 - gmso/parameterization/parameterize.py | 260 ++---------------- ...meterizer.py => topology_parameterizer.py} | 152 ++++++---- gmso/parameterization/utils.py | 1 + gmso/tests/parameterization/__init__.py | 0 gmso/tests/parameterization/test_opls_gmso.py | 54 ++++ .../test_subtopology_utils.py | 110 ++++++++ .../parameterization/test_trappe_foyer.py | 0 8 files changed, 282 insertions(+), 300 deletions(-) rename gmso/parameterization/{parameterizer.py => topology_parameterizer.py} (65%) create mode 100644 gmso/tests/parameterization/__init__.py create mode 100644 gmso/tests/parameterization/test_opls_gmso.py create mode 100644 gmso/tests/parameterization/test_subtopology_utils.py create mode 100644 gmso/tests/parameterization/test_trappe_foyer.py diff --git a/gmso/parameterization/foyer_utils.py b/gmso/parameterization/foyer_utils.py index bc1b695bb..6b86c957c 100644 --- a/gmso/parameterization/foyer_utils.py +++ b/gmso/parameterization/foyer_utils.py @@ -6,10 +6,8 @@ from gmso.core.atom import Atom from gmso.parameterization.subtopology_utils import subtop_bonds -from gmso.utils.decorators import experimental -@experimental def get_topology_graph(gmso_topology, atomdata_populator=None): """Return a TopologyGraph with relevant attributes from an GMSO topology. @@ -67,7 +65,6 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): return top_graph -@experimental def get_topology_graph_from_subtop(subtopology): """Get an equivalent topology graph for a sub-topology.""" subtop_named_tuple = namedtuple("subtopology", ("sites", "bonds")) @@ -76,7 +73,6 @@ def get_topology_graph_from_subtop(subtopology): ) -@experimental def get_atomtyping_rules_provider(gmso_ff): """Return a foyer AtomTypingRulesProvider from a GMSO forcefield. @@ -105,7 +101,6 @@ def get_atomtyping_rules_provider(gmso_ff): ) -@experimental def typemap_dict(topology_graph, atomtyping_rules_provider, max_iter=10): """Return a dictionary of typemap, by finding atomtypes in foyer.""" return find_atomtypes(topology_graph, atomtyping_rules_provider, max_iter) diff --git a/gmso/parameterization/parameterize.py b/gmso/parameterization/parameterize.py index 08e0258e4..f6b9a673d 100644 --- a/gmso/parameterization/parameterize.py +++ b/gmso/parameterization/parameterize.py @@ -1,14 +1,7 @@ """Functions used to atomtype a gmso.Topology.""" -from operator import attrgetter - -import networkx as nx -from foyer.atomtyper import AtomTypingRulesProvider, find_atomtypes -from foyer.smarts import SMARTS -from foyer.topology_graph import TopologyGraph - -from gmso import Atom -from gmso.parameterization.isomorph import ( # check refs - partition_isomorphic_topology_graphs, +from gmso.parameterization.topology_parameterizer import ( + TopologyParameterizationConfig, + TopologyParameterizer, ) __all__ = ["apply"] @@ -17,6 +10,7 @@ def apply( top, forcefields, + identify_connections=False, identify_connected_components=True, use_residue_info=False, assert_bond_params=True, @@ -33,6 +27,12 @@ def apply( forcefields: dict, required The keys are labels that match the subtopology label or site residue_name, and the values are gmso Forcefield objects that gets applied to the specified molecule + + identify_connections: bool, optional, default=False + If true, add connections identified using networkx graph matching to match + the topology's bonding graph to smaller sub-graphs that correspond to an angle, + dihedral, improper etc + identify_connected_components: bool, optional, default=True A flag to determine whether or not to search the topology for repeated disconnected structures, otherwise known as molecules and type each molecule only once. @@ -52,233 +52,21 @@ def apply( If True, an error is raised if parameters are not found for all system improper dihedrals. """ - top.identify_connections() - top_graph = get_modified_topology_graph(top) - if identify_connected_components: # True, (True,False) - # populate molecule information - isomorphic_top_graphs = partition_isomorphic_topology_graphs(top_graph) - apply_atomtypes_from_isomorphic( - top, isomorphic_top_graphs, forcefields - ) # what if only one forcefield?, traversing done in this function - elif not use_residue_info: # False, False - # get entire topgraph - topgraph = get_modified_topology_graph(top) - # get rules provider from gmso_ff - # Assumes forcefields is just one forcefield - ff_of_interest = next(iter(forcefields.values())) - at_rules_provider = get_atomtyping_rules_provider( - gmso_ff=ff_of_interest + config = TopologyParameterizationConfig.parse_obj( + dict( + identify_connections=identify_connections, + identify_connected_components=identify_connected_components, + use_residue_info=use_residue_info, + assert_bond_params=assert_bond_params, + assert_angle_params=assert_angle_params, + assert_dihedral_params=assert_dihedral_params, + assert_improper_params=assert_improper_params, ) - # create a typemap of that topgraph - typemap = find_atomtypes(topgraph, at_rules_provider) - # apply typemap - traverse_typemap(top, typemap, ff_of_interest) - elif use_residue_map: # False, True - """Not Implemented - maps = {} - for key, val in forcefields.items(): - # generate topgraph of items with "key" residue name - # NOTE CAL: this iter_sites_by method won't work currently, can only iter one at a time. - res_iter = top.iter_sites_by(["residue_name", "residue_number"], [key, 1]) #Need to note that resnumber index 1 - # create a top/graph from this set of sites - #subgraph = FUNCTION(res_iter) - # append that typemap to a dict with residue name as key - #maps[key] = find_atomtypes(subgraph, val) - # generate a full typemap by iterating through all sites and copying over map info from `maps` - #typemap = FUNCTION(top, maps) - # apply typemap - traverse_typemap(top, typemap, forcefields) - """ - raise ( - GMSOError, - "Using site residue matching to match substructures to a given forcefield is not implemented yet", - ) - - if assert_bond_params: - apply_connection_types(top, forcefields, "bond") - if assert_angle_params: - apply_connection_types(top, forcefields, "angle") - if assert_dihedral_params: - apply_connection_types(top, forcefields, "dihedral") - if assert_improper_params: - apply_connection_types(top, forcefields, "improper") - - return top - - -def get_modified_topology_graph(gmso_topology): - """Return a TopologyGraph with relevant attributes from an GMSO topology. - - Parameters - ---------- - gmso_topology: gmso.Topology - The GMSO Topology - - Returns - ------- - TopologyGraph - The equivalent TopologyGraph of the openFF Topology `openff_topology` - """ - top_graph = TopologyGraph() - for atom in gmso_topology.sites: - if isinstance(atom, Atom): - if atom.name.startswith("_"): - top_graph.add_atom( - name=atom.name, - index=gmso_topology.get_index(atom), - atomic_number=None, - element=atom.name, - subtopology_label=atom.label, # TODO: Change based on PR #638 conventions - ) - - else: - top_graph.add_atom( - name=atom.name, - index=gmso_topology.get_index(atom), - atomic_number=atom.element.atomic_number, - element=atom.element.symbol, - subtopology_label=atom.label, # TODO: Change based on PR #638 conventions - ) - - for top_bond in gmso_topology.bonds: - atoms_indices = [ - gmso_topology.get_index(atom) - for atom in top_bond.connection_members - ] - top_graph.add_bond(atoms_indices[0], atoms_indices[1]) - - return top_graph - - -def get_atomtyping_rules_provider(gmso_ff): - """Return a foyer AtomTypingRulesProvider from a GMSO forcefield. - - Parameters - ---------- - gmso_ff: gmso.core.forcefield.Forcefield - The GMSO forcefield object to extract the rules from - Returns - ------- - AtomTypingRulesProvider - The foyer.atomtyper.AtomTypingRulesProvider object used to parse atomtype definitions. - Typically, SMARTS is the ruleset of choice. See https://github.com/mosdef-hub/foyer/issues/63 - for curently supported features in Foyer. - """ - atom_type_defs = {} - atom_type_overrides = {} - parser = SMARTS({}) - for atom_type_name, atom_type in gmso_ff.atom_types.items(): - if atom_type.definition: - atom_type_defs[atom_type_name] = atom_type.definition - if atom_type.overrides: - atom_type_overrides[atom_type_name] = atom_type.overrides - - return AtomTypingRulesProvider( - atom_type_defs, atom_type_overrides, {}, parser + ) + parameterizer = TopologyParameterizer( + topology=top, forcefields=forcefields, config=config ) + parameterizer.run_parameterization() -def apply_atomtypes_from_isomorphic(top, isomorphic_top_graphs, forcefields): - """Set Topology atomtypes from isomorphic structures matching supplied forcefields. - - Parameters - ---------- - top: gmso.core.topology.Topology, required - The GMSO topology on which to apply atomtypes - isomorphic_top_graphs: dict, required - The dictionary mapping which breaks up molecules into one graph structure, and a list of - the corresponding subgraphs in the TopologyGraph that are isomorphic - forcefields: dict, required - The keys are labels that match the subtopology labels for each isomorphic graph, and the - values are gmso Forcefield objects that gets applied to the specified molecule - """ - for top_graph, identical in isomorphic_top_graphs.items(): - sub_top_label = next( - iter( - nx.classes.function.get_node_attributes( - top_graph, "atom_data" - ).values() - ) - ).subtopology_label - ff_of_interest = forcefields[sub_top_label] - at_rules_provider = get_atomtyping_rules_provider( - gmso_ff=ff_of_interest - ) - sub_type_map = find_atomtypes(top_graph, at_rules_provider) - traverse_typemap(top, sub_type_map, ff_of_interest) - for graph, mapping in identical: - for node in graph: - mapped = sub_type_map[mapping[node]]["atom_type"] - top._sites[node].atom_type = mapped.clone() - - -def traverse_typemap(top, type_map, forcefield): - """Take a typemap and applies those atomtypes to the Topology. - - Parameters - ---------- - top: gmso.core.topology.Topology, required - The GMSO topology on which to apply atomtypes - typemap: dict, required - The dictionary of atom_index with iterated matching to get the atomtype of that site - forcefield: gmso.forcefield.Forcefield, required - The GMSO forcefield object which has information associated with each unique atomtype - """ - for index in type_map: - site_of_interest = top._sites[index] - site_of_interest.atom_type = forcefield.atom_types[ - type_map[index]["atomtype"] - ].clone() - type_map[index]["atom_type"] = site_of_interest.atom_type - - -def apply_residue_names(top, more): # TODO - """Take a topology and applies residue info to all sites. - - Parameters - ---------- - top: gmso.core.topology.Topology, required - The GMSO topology on which to apply atomtypes - more: TBD, required - The information necessary to map residue name and number for each site - """ - return - - -def apply_connection_types(top, forcefields, parameter): - """Set Topology connection types from matching atomclass information in forcefield. - - Parameters - ---------- - top: gmso.core.topology.Topology, required - The GMSO topology on which to apply paramter_types - forcefields: dict, required - The keys are labels that match the subtopology labels for each molecule, and the - values are gmso Forcefield objects that gets applied to the specified molecule - parameter: str, required - The connection type to parameterize. Can any of "bond", "angle", "dihedral", "improper" - """ - connect_members = { - "bond": [0, 1], - "angle": [0, 1, 2], - "dihedral": [0, 1, 2, 3], - "improper": [0, 1, 2, 3], - } # TODO validate improper order - molecule_label = "label_" # place to grab info referencing molecule name - for connect in getattr(top, "_" + parameter + "s"): - sub_top_label = connect.connection_members[0].get(molecule_label) - if sub_top_label is None: - msg = f"Site {connect.connection_members[0]}in this connection has no molecule_label {molecule_label}.\ - This string should correspond to one of the keys associated with the forcefield that has been \ - passed. Currently provided name:GMSOForcefields pairs are {forcefields.items()}" - raise ParameterizationError(msg) - ff_of_interest = forcefields[sub_top_label] - type_getter = attrgetter("atom_type.atomclass") - member_types = [ - type_getter(getattr(connect, "connection_members")[i]) - for i in connect_members[parameter] - ] - ctype = ff_of_interest.get_potential( - group=parameter + "_type", key=member_types - ).clone() - setattr(connect, parameter + "_type", ctype) + return parameterizer.topology diff --git a/gmso/parameterization/parameterizer.py b/gmso/parameterization/topology_parameterizer.py similarity index 65% rename from gmso/parameterization/parameterizer.py rename to gmso/parameterization/topology_parameterizer.py index 30d09d7ea..a34c152e0 100644 --- a/gmso/parameterization/parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -30,7 +30,63 @@ class GMSOParameterizationError(GMSOError): """Raise when parameterization fails.""" -class Parameterizer(GMSOBase): +class TopologyParameterizationConfig(GMSOBase): + """Configuration options for parameterizing a topology.""" + + clone_topology: bool = Field( + default=False, + description="If true, clone the topology and apply parameters to the cloned one.", + ) # Unused + + identify_connections: bool = Field( + default=False, + description="If true, add connections identified using networkx graph matching to match" + "the topology's bonding graph to smaller sub-graphs that correspond to an " + "angle, dihedral, improper etc", + ) + + identify_connected_components: bool = Field( + default=False, + description="A flag to determine whether or not to search the topology" + " for repeated disconnected structures, otherwise known as " + "molecules and type each molecule only once.", + ) + + use_residue_info: bool = Field( + default=False, + description="A flag to determine whether or not to look at site.residue_name " + "to look parameterize each molecule only once. Will only be used if " + "identify_connected_components=False", + ) # Unused + + assert_bond_params: bool = Field( + default=True, + description="If True, an error is raised if parameters are not found for " + "all system bonds.", + ) + + assert_angle_params: bool = Field( + default=True, + description="If True, an error is raised if parameters are not found for " + "all system angles", + ) + + assert_dihedral_params: bool = ( + Field( + default=True, + description="If True, an error is raised if parameters are not found for " + "all system dihedrals.", + ), + ) + + assert_improper_params: bool = Field( + default=False, + description="If True, an error is raised if parameters are not found for " + "all system impropers.", + ) + + +class TopologyParameterizer(GMSOBase): """Utility class to parameterize a topology with gmso Forcefield.""" topology: Topology = Field(..., description="The gmso topology.") @@ -42,6 +98,10 @@ class Parameterizer(GMSOBase): "should match the subtopology names", ) + config: TopologyParameterizationConfig = Field( + ..., description="The configuration options for the parameterizer." + ) + def get_ff(self, key=None): """Return the forcefield of choice by looking up the forcefield dictionary.""" if isinstance(self.forcefields, Dict): @@ -49,23 +109,14 @@ def get_ff(self, key=None): else: return self.forcefields - def _parameterize_sites(self, top_or_subtop, typemap, ff): + def _parameterize_sites(self, sites, typemap, ff): """Parameterize sites with appropriate atom-types from the forcefield.""" - for j, site in enumerate(top_or_subtop.sites): + for j, site in enumerate(sites): site.atom_type = ff.get_potential( - "atom_type", typemap[j]["atom_type"] + "atom_type", typemap[j]["atomtype"] ).clone() # Always properly indexed or not? - def _parameterize_connections( - self, - top_or_subtop, - ff, - is_subtop=False, - assert_bond_params=True, - assert_angle_params=True, - assert_dihedral_params=True, - assert_improper_params=False, - ): + def _parameterize_connections(self, top_or_subtop, ff, is_subtop=False): """Parameterize connections with appropriate potentials from the forcefield.""" if is_subtop: bonds = subtop_bonds(top_or_subtop) @@ -78,10 +129,18 @@ def _parameterize_connections( dihedrals = top_or_subtop.dihedrals impropers = top_or_subtop.impropers - self._apply_connection_parameters(bonds, ff, assert_bond_params) - self._apply_connection_parameters(angles, ff, assert_angle_params) - self._apply_connection_parameters(dihedrals, ff, assert_dihedral_params) - self._apply_connection_parameters(impropers, ff, assert_improper_params) + self._apply_connection_parameters( + bonds, ff, self.config.assert_bond_params + ) + self._apply_connection_parameters( + angles, ff, self.config.assert_angle_params + ) + self._apply_connection_parameters( + dihedrals, ff, self.config.assert_dihedral_params + ) + self._apply_connection_parameters( + impropers, ff, self.config.assert_improper_params + ) def _apply_connection_parameters( self, connections, ff, error_on_missing=True @@ -103,59 +162,36 @@ def _apply_connection_parameters( group=group, key=identifier_key, warn=True ) if match: - visited[tuple[identifier_key]] = match + visited[tuple(identifier_key)] = match break if not match and error_on_missing: raise GMSOParameterizationError( f"No parameters found for connection {connection} in the Forcefield." ) - setattr(connection, group, match.clone()) - - def _parameterize( - self, - subtop_or_top, - typemap, - is_subtop=False, - assert_bond_params=True, - assert_angle_params=True, - assert_dihedral_params=True, - assert_improper_params=False, - ): + elif match: + setattr(connection, group, match.clone()) + + def _parameterize(self, subtop_or_top, typemap, is_subtop=False): """Parameterize a topology/subtopology based on an atomtype map.""" forcefield = self.get_ff(subtop_or_top.name) self._parameterize_sites(subtop_or_top.sites, typemap, forcefield) self._parameterize_connections( - subtop_or_top, - forcefield, - is_subtop=is_subtop, - assert_bond_params=assert_bond_params, - assert_angle_params=assert_angle_params, - assert_dihedral_params=assert_dihedral_params, - assert_improper_params=assert_improper_params, + subtop_or_top, forcefield, is_subtop=is_subtop ) - def apply( - self, - identify_connected_components=True, - use_residue_info=False, - assert_bond_params=True, - assert_angle_params=True, - assert_dihedral_params=True, - assert_improper_params=False, - ): - """Apply the current forcefield(s) to the topology/ various subtopologies.""" - parameterize_kwargs = dict( - assert_bond_params=assert_bond_params, - assert_angle_params=assert_angle_params, - assert_dihedral_params=assert_dihedral_params, - assert_improper_params=assert_improper_params, - ) + def run_parameterization(self): + """Run parameterization of the topology with give forcefield(s) and configuration.""" if self.topology.is_typed(): raise GMSOParameterizationError( "Cannot parameterize a typed topology. Please provide a topology without any types" ) + if self.config.identify_connections: + """ToDo: This mutates the topology and is agnostic to downstream + errors. So, here we should use index only option""" + self.topology.identify_connections() + if isinstance(self.forcefields, Dict): for subtop in self.topology.subtops: if subtop.name not in self.forcefields: @@ -168,27 +204,25 @@ def apply( typemap = self._get_atomtypes( self.get_ff(subtop.name), subtop, - identify_connected_components, + self.config.identify_connected_components, is_subtop=True, ) self._parameterize( typemap, subtop, is_subtop=True, # This will be removed from the future iterations - **parameterize_kwargs, ) else: typemap = self._get_atomtypes( self.get_ff(), self.topology, - identify_connected_components, + self.config.identify_connected_components, is_subtop=False, ) self._parameterize( self.topology, typemap, is_subtop=False, # This will be removed from the future iterations - **parameterize_kwargs, ) @staticmethod @@ -212,7 +246,7 @@ def connection_identifier( def _get_atomtypes( forcefield, topology, use_isomprohic_checks=False, is_subtop=False ): - """Run atom-typing in foyer and return a typemap.""" + """Run atom-typing in foyer and return the typemap.""" atom_typing_rules_provider = get_atomtyping_rules_provider(forcefield) if is_subtop: diff --git a/gmso/parameterization/utils.py b/gmso/parameterization/utils.py index 3fb94404e..d04d10659 100644 --- a/gmso/parameterization/utils.py +++ b/gmso/parameterization/utils.py @@ -1,5 +1,6 @@ """Generic utilities for parameterizing a gmso Topology.""" +from gmso.abc.gmso_base import GMSOBase from gmso.core.angle import Angle from gmso.core.atom import Atom from gmso.core.bond import Bond diff --git a/gmso/tests/parameterization/__init__.py b/gmso/tests/parameterization/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py new file mode 100644 index 000000000..17d386d02 --- /dev/null +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -0,0 +1,54 @@ +import glob +from pathlib import Path + +import parmed as pmd +import pytest +from forcefield_utilities.xml_loader import FoyerFFs +from pkg_resources import resource_filename + +from gmso.external.convert_parmed import from_parmed +from gmso.parameterization.parameterize import apply +from gmso.tests.base_test import BaseTest + + +def get_foyer_opls_test_dirs(): + all_dirs = glob.glob(resource_filename("foyer", "opls_validation") + "/*") + with open( + resource_filename("foyer", "tests/implemented_opls_tests.txt") + ) as impl_file: + correctly_implemented = set(impl_file.read().strip().split("\n")) + + parent_dirs = map(Path, all_dirs) + parent_dirs = list( + filter( + lambda p: p.name in correctly_implemented + and (p / f"{p.name}.top").exists(), + parent_dirs, + ) + ) + return parent_dirs + + +class TestOPLSGMSO(BaseTest): + @pytest.fixture(scope="session") + def xml_loader(self): + return FoyerFFs() + + @pytest.fixture(scope="session") + def oplsaa(self, xml_loader): + return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() + + @pytest.mark.parametrize( + "system_dir", get_foyer_opls_test_dirs(), ids=lambda p: p.name + ) + def test_foyer_oplsaa_files( + self, system_dir, oplsaa, are_equivalent_topologies + ): + top_file = str(system_dir / f"{system_dir.name}.top") + gro_file = str(system_dir / f"{system_dir.name}.gro") + struct = pmd.load_file(top_file, xyz=gro_file) + gmso_top_unparam = from_parmed(struct, refer_type=False) + apply(gmso_top_unparam, oplsaa, identify_connected_components=False) + assert are_equivalent_topologies( + gmso_top_unparam, from_parmed(struct, refer_type=True) + ), f"{system_dir} Failed" diff --git a/gmso/tests/parameterization/test_subtopology_utils.py b/gmso/tests/parameterization/test_subtopology_utils.py new file mode 100644 index 000000000..b6aba7049 --- /dev/null +++ b/gmso/tests/parameterization/test_subtopology_utils.py @@ -0,0 +1,110 @@ +import mbuild as mb +import pytest + +from gmso.external.convert_mbuild import from_mbuild +from gmso.parameterization.subtopology_utils import ( + _members_in_subtop, + assert_no_boundary_bonds, + subtop_angles, + subtop_bonds, + subtop_dihedrals, + subtop_impropers, +) +from gmso.tests.base_test import BaseTest +from gmso.utils.connectivity import identify_connections + + +class TestSubTopologyUtils(BaseTest): + @pytest.fixture(scope="session") + def ethane_box_gmso(self): + ethane_box = mb.fill_box( + mb.lib.molecules.Ethane(), n_compounds=20, density=2 + ) + ethane_box_gmso = from_mbuild(ethane_box) + identify_connections(ethane_box_gmso) + return ethane_box_gmso + + def test_no_boundary_bonds_ethane(self, ethane): + with pytest.raises(AssertionError): + assert_no_boundary_bonds(ethane.subtops[0]) + + def test_no_boundary_bonds_ethane_box(self, ethane_box_gmso): + for subtop in ethane_box_gmso.subtops: + assert_no_boundary_bonds(subtop) + + def test_subtopology_bonds(self, ethane_box_gmso): + for subtop in ethane_box_gmso.subtops: + bonds = list(subtop_bonds(subtop)) + assert len(bonds) == 7 + for bond in bonds: + assert _members_in_subtop(bond, subtop) + + bond_members = map( + lambda b: tuple(map(lambda s: s.name, b.connection_members)), + bonds, + ) + expected_members = {("C", "H"), ("C", "C"), ("H", "C")} + assert all( + b_member in expected_members for b_member in bond_members + ) + + def test_subtopology_angles(self, ethane_box_gmso): + for subtop in ethane_box_gmso.subtops: + angles = list(subtop_angles(subtop)) + assert len(list(angles)) == 12 + for angle in angles: + assert _members_in_subtop(angle, subtop) + + angle_members = map( + lambda a: tuple(map(lambda s: s.name, a.connection_members)), + angles, + ) + expected_members = { + ("H", "C", "H"), + ("H", "C", "C"), + ("C", "C", "H"), + } + assert all( + a_member in expected_members for a_member in angle_members + ) + + def test_subtopology_dihedrals(self, ethane_box_gmso): + for subtop in ethane_box_gmso.subtops: + dihedrals = list(subtop_dihedrals(subtop)) + assert len(dihedrals) == 9 + for dihedral in dihedrals: + assert _members_in_subtop(dihedral, subtop) + + dihedral_members = map( + lambda d: tuple(map(lambda s: s.name, d.connection_members)), + dihedrals, + ) + expected_members = {("H", "C", "C", "H")} + assert all( + a_member in expected_members for a_member in dihedral_members + ) + + def test_subtopology_impropers(self, ethane_box_gmso): + for subtop in ethane_box_gmso.subtops: + impropers = list(subtop_impropers(subtop)) + assert len(impropers) == 8 + for improper in impropers: + assert _members_in_subtop(improper, subtop) + + improper_members = list( + map( + lambda i: tuple( + map(lambda s: s.name, i.connection_members) + ), + impropers, + ) + ) + expected_members = { + ("C", "C", "H", "H"), + ("C", "H", "H", "C"), + ("C", "H", "H", "H"), + ("C", "H", "C", "H"), + } + assert all( + a_member in expected_members for a_member in improper_members + ) diff --git a/gmso/tests/parameterization/test_trappe_foyer.py b/gmso/tests/parameterization/test_trappe_foyer.py new file mode 100644 index 000000000..e69de29bb From cc46491e8cf8478bb4c9d26e2310128a4bd8a38d Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Mon, 9 May 2022 10:06:29 -0500 Subject: [PATCH 10/32] WIP- Better error handling --- gmso/parameterization/topology_parameterizer.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index a34c152e0..bd4bf9cfa 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -167,7 +167,8 @@ def _apply_connection_parameters( if not match and error_on_missing: raise GMSOParameterizationError( - f"No parameters found for connection {connection} in the Forcefield." + f"No parameters found for connection {connection}, group: {group}, " + f"identifiers: {connection_identifiers} in the Forcefield." ) elif match: setattr(connection, group, match.clone()) From f5633e7ccbfd44f25f0cfc647c9aca2c98caec49 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 1 Jun 2022 13:17:03 -0500 Subject: [PATCH 11/32] Misc Changes 1. Update Topology after parametrization 2. Add dependency for forcefield_utilities in env files 3. Add tests for trappe forcefield 4. Patch parmed (should be moved to a cal's PR #658 --- environment-dev.yml | 3 +- environment.yml | 1 + gmso/external/convert_parmed.py | 18 +++-- .../topology_parameterizer.py | 1 + .../parameterization/test_trappe_foyer.py | 0 .../parameterization/test_trappe_gmso.py | 74 +++++++++++++++++++ 6 files changed, 90 insertions(+), 7 deletions(-) delete mode 100644 gmso/tests/parameterization/test_trappe_foyer.py create mode 100644 gmso/tests/parameterization/test_trappe_gmso.py diff --git a/environment-dev.yml b/environment-dev.yml index dd2ed3d7d..d1c5b8a4b 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -12,7 +12,7 @@ dependencies: - pytest - mbuild >= 0.11.0 - openbabel >= 3.0.0 - - foyer >= 0.9.4 + - foyer >= 0.11.1 - gsd >= 2.0 - parmed >= 3.4.3 - pytest-cov @@ -22,3 +22,4 @@ dependencies: - ipywidgets - ele >= 0.2.0 - pre-commit + - forcefield-utilities diff --git a/environment.yml b/environment.yml index df1145fcf..869017908 100644 --- a/environment.yml +++ b/environment.yml @@ -10,3 +10,4 @@ dependencies: - pydantic < 1.9.0 - networkx - ele >= 0.2.0 + - forcefield-utilities diff --git a/gmso/external/convert_parmed.py b/gmso/external/convert_parmed.py index 403fa4ace..0f5ef4ae2 100644 --- a/gmso/external/convert_parmed.py +++ b/gmso/external/convert_parmed.py @@ -455,11 +455,15 @@ def to_parmed(top, refer_type=True): # Set up Parmed structure and define general properties structure = pmd.Structure() structure.title = top.name - structure.box = np.concatenate( - ( - top.box.lengths.to("angstrom").value, - top.box.angles.to("degree").value, + structure.box = ( + np.concatenate( + ( + top.box.lengths.to("angstrom").value, + top.box.angles.to("degree").value, + ) ) + if top.box + else None ) # Maps @@ -479,8 +483,10 @@ def to_parmed(top, refer_type=True): pmd_atom = pmd.Atom( atomic_number=atomic_number, name=site.name, - mass=site.mass.to(u.amu).value, - charge=site.charge.to(u.elementary_charge).value, + mass=site.mass.to(u.amu).value if site.mass else None, + charge=site.charge.to(u.elementary_charge).value + if site.charge + else None, ) pmd_atom.xx, pmd_atom.xy, pmd_atom.xz = site.position.to( "angstrom" diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index bd4bf9cfa..085af1ba3 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -225,6 +225,7 @@ def run_parameterization(self): typemap, is_subtop=False, # This will be removed from the future iterations ) + self.topology.update_topology() @staticmethod def connection_identifier( diff --git a/gmso/tests/parameterization/test_trappe_foyer.py b/gmso/tests/parameterization/test_trappe_foyer.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py new file mode 100644 index 000000000..47e0e5d86 --- /dev/null +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -0,0 +1,74 @@ +import glob +from pathlib import Path + +import foyer +import parmed as pmd +import pytest +import unyt as u +from forcefield_utilities.xml_loader import FoyerFFs +from pkg_resources import resource_filename + +import gmso +from gmso.external.convert_parmed import from_parmed, to_parmed +from gmso.parameterization.parameterize import apply +from gmso.tests.base_test import BaseTest + + +def get_foyer_trappe_test_dirs(): + all_dirs = glob.glob(resource_filename("foyer", "trappe_validation") + "/*") + with open( + resource_filename("foyer", "tests/implemented_trappe_tests.txt") + ) as impl_file: + correctly_implemented = set(impl_file.read().strip().split("\n")) + + parent_dirs = map(Path, all_dirs) + parent_dirs = list( + filter( + lambda p: p.name in correctly_implemented + and (p / f"{p.name}.mol2").exists(), + parent_dirs, + ) + ) + return parent_dirs + + +class TestTrappeGMSO(BaseTest): + @pytest.fixture(scope="session") + def xml_loader(self): + return FoyerFFs() + + @pytest.fixture(scope="session") + def trappe_ua_gmso(self, xml_loader): + return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() + + @pytest.fixture(scope="session") + def trappe_ua_openmm_foyer(self): + return foyer.forcefields.load_TRAPPE_UA() + + @pytest.mark.parametrize( + "system_dir", get_foyer_trappe_test_dirs(), ids=lambda p: p.name + ) + def test_foyer_trappe_files( + self, + system_dir, + trappe_ua_openmm_foyer, + trappe_ua_gmso, + are_equivalent_connections, + ): + mol2_file = system_dir / f"{system_dir.name}.mol2" + top_gmso = gmso.Topology.load(mol2_file) + struct_pmd = trappe_ua_openmm_foyer.apply(to_parmed(top_gmso)) + apply(top_gmso, trappe_ua_gmso, identify_connected_components=False) + gmso_top_from_parmeterized_pmd = from_parmed(struct_pmd) + for atom, mirror in zip( + top_gmso.sites, gmso_top_from_parmeterized_pmd.sites + ): + assert atom.name == mirror.name + assert u.allclose_units(atom.mass, mirror.mass, 1e-4) + + assert u.allclose_units(atom.atom_type.charge, mirror.charge, 1e-4) + + atom_params = atom.atom_type.get_parameters() + mirror_params = mirror.atom_type.get_parameters() + for k in atom_params: + assert u.allclose_units(atom_params[k], mirror_params[k]) From 72137ee02ef8094dc4b8462f853cd46fc3ec7662 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 1 Jun 2022 16:32:52 -0500 Subject: [PATCH 12/32] WIP- Skip opls for now, full TrappE tests --- gmso/tests/base_test.py | 1 - gmso/tests/parameterization/test_opls_gmso.py | 1 + .../parameterization/test_trappe_gmso.py | 49 +++++++++++++++++-- 3 files changed, 45 insertions(+), 6 deletions(-) diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index 4c4702993..68faa73e6 100644 --- a/gmso/tests/base_test.py +++ b/gmso/tests/base_test.py @@ -340,7 +340,6 @@ def test_atom_equality(atom1, atom2): ) for prop in atom1.dict(by_alias=True): if not equal(atom2.dict().get(prop), atom1.dict().get(prop)): - print(atom2.dict().get(prop), atom1.dict().get(prop), prop) return False return True diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py index 17d386d02..84e52dae4 100644 --- a/gmso/tests/parameterization/test_opls_gmso.py +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -29,6 +29,7 @@ def get_foyer_opls_test_dirs(): return parent_dirs +@pytest.mark.skip class TestOPLSGMSO(BaseTest): @pytest.fixture(scope="session") def xml_loader(self): diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py index 47e0e5d86..cb4a9287b 100644 --- a/gmso/tests/parameterization/test_trappe_gmso.py +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -2,7 +2,6 @@ from pathlib import Path import foyer -import parmed as pmd import pytest import unyt as u from forcefield_utilities.xml_loader import FoyerFFs @@ -32,6 +31,38 @@ def get_foyer_trappe_test_dirs(): return parent_dirs +def _match_connection_parameters( + from_gmso_top, from_parmed_top, connection_type="bonds" +): + connection_types_original = {} + connection_types_mirror = {} + for connection in getattr(from_parmed_top, connection_type): + connection_types_mirror[ + tuple( + from_parmed_top.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for connection in getattr(from_gmso_top, connection_type): + connection_types_original[ + tuple( + from_gmso_top.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for key in connection_types_original: + conn = connection_types_original[key] + print(key) + conn_mirror = connection_types_mirror[key] + for param in conn_mirror.bond_type.parameters: + assert u.allclose_units( + conn_mirror.bond_type.parameters[param], + conn.bond_type.parameters[param], + ) + + class TestTrappeGMSO(BaseTest): @pytest.fixture(scope="session") def xml_loader(self): @@ -56,12 +87,12 @@ def test_foyer_trappe_files( are_equivalent_connections, ): mol2_file = system_dir / f"{system_dir.name}.mol2" - top_gmso = gmso.Topology.load(mol2_file) - struct_pmd = trappe_ua_openmm_foyer.apply(to_parmed(top_gmso)) - apply(top_gmso, trappe_ua_gmso, identify_connected_components=False) + gmso_top = gmso.Topology.load(mol2_file) + struct_pmd = trappe_ua_openmm_foyer.apply(to_parmed(gmso_top)) + apply(gmso_top, trappe_ua_gmso, identify_connected_components=False) gmso_top_from_parmeterized_pmd = from_parmed(struct_pmd) for atom, mirror in zip( - top_gmso.sites, gmso_top_from_parmeterized_pmd.sites + gmso_top.sites, gmso_top_from_parmeterized_pmd.sites ): assert atom.name == mirror.name assert u.allclose_units(atom.mass, mirror.mass, 1e-4) @@ -72,3 +103,11 @@ def test_foyer_trappe_files( mirror_params = mirror.atom_type.get_parameters() for k in atom_params: assert u.allclose_units(atom_params[k], mirror_params[k]) + + _match_connection_parameters(gmso_top, gmso_top_from_parmeterized_pmd) + _match_connection_parameters( + gmso_top, gmso_top_from_parmeterized_pmd, "angles" + ) + _match_connection_parameters( + gmso_top, gmso_top_from_parmeterized_pmd, "dihedrals" + ) From 263ef780f3ccad3aabd6d111dcced6e0d1186f74 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 3 Jun 2022 09:40:53 -0500 Subject: [PATCH 13/32] Avoid accidental overwriting of typemap when using isomorphism --- gmso/parameterization/topology_parameterizer.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index 085af1ba3..742baa6e8 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -260,10 +260,13 @@ def _get_atomtypes( isomorphic_substructures = partition_isomorphic_topology_graphs( foyer_topology_graph ) + typemap = {} for graph, mirrors in isomorphic_substructures.items(): - typemap = typemap_dict( - atomtyping_rules_provider=atom_typing_rules_provider, - topology_graph=graph, + typemap.update( + typemap_dict( + atomtyping_rules_provider=atom_typing_rules_provider, + topology_graph=graph, + ) ) for mirror, mapping in mirrors: for node in mirror: From c7ea34632ab6fdc1d1a977354cdc77acf2c17e87 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 3 Jun 2022 09:55:07 -0500 Subject: [PATCH 14/32] WIP- Remove unused import --- gmso/parameterization/utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gmso/parameterization/utils.py b/gmso/parameterization/utils.py index d04d10659..3fb94404e 100644 --- a/gmso/parameterization/utils.py +++ b/gmso/parameterization/utils.py @@ -1,6 +1,5 @@ """Generic utilities for parameterizing a gmso Topology.""" -from gmso.abc.gmso_base import GMSOBase from gmso.core.angle import Angle from gmso.core.atom import Atom from gmso.core.bond import Bond From d267d2fc7da2bccea3b07ac44e462909034b7d16 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 3 Jun 2022 15:48:12 -0500 Subject: [PATCH 15/32] Use enumerate for atom index while converting to TopologyGraph --- gmso/parameterization/foyer_utils.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/gmso/parameterization/foyer_utils.py b/gmso/parameterization/foyer_utils.py index 6b86c957c..3f25b3610 100644 --- a/gmso/parameterization/foyer_utils.py +++ b/gmso/parameterization/foyer_utils.py @@ -30,7 +30,9 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): A light networkx representation of the topology """ top_graph = TopologyGraph() - for atom in gmso_topology.sites: + atom_index_map = {} + for j, atom in enumerate(gmso_topology.sites): + atom_index_map[id(atom)] = j if isinstance(atom, Atom): kwargs = ( atomdata_populator(gmso_topology, atom) @@ -40,7 +42,7 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): if atom.name.startswith("_"): top_graph.add_atom( name=atom.name, - index=gmso_topology.get_index(atom), + index=j, # Assumes order is preserved atomic_number=None, element=atom.name, **kwargs, @@ -49,7 +51,7 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): else: top_graph.add_atom( name=atom.name, - index=gmso_topology.get_index(atom), + index=j, # Assumes order is preserved atomic_number=atom.element.atomic_number, element=atom.element.symbol, **kwargs, @@ -57,8 +59,7 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): for top_bond in gmso_topology.bonds: atoms_indices = [ - gmso_topology.get_index(atom) - for atom in top_bond.connection_members + atom_index_map[id(atom)] for atom in top_bond.connection_members ] top_graph.add_bond(atoms_indices[0], atoms_indices[1]) From 9aab52841b8324e07e1971dd4edd4dea476824f8 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 3 Jun 2022 16:06:18 -0500 Subject: [PATCH 16/32] Fix argument order --- gmso/parameterization/topology_parameterizer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index 742baa6e8..cb0298d71 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -209,8 +209,8 @@ def run_parameterization(self): is_subtop=True, ) self._parameterize( - typemap, subtop, + typemap, is_subtop=True, # This will be removed from the future iterations ) else: From 51b10815361eb3f4e488617355bd5abf4725bce9 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Fri, 3 Jun 2022 16:34:18 -0500 Subject: [PATCH 17/32] WIP- Add test for subtopology parameterization --- .../test_subtopology_parameterization.py | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 gmso/tests/parameterization/test_subtopology_parameterization.py diff --git a/gmso/tests/parameterization/test_subtopology_parameterization.py b/gmso/tests/parameterization/test_subtopology_parameterization.py new file mode 100644 index 000000000..4443490ca --- /dev/null +++ b/gmso/tests/parameterization/test_subtopology_parameterization.py @@ -0,0 +1,23 @@ +import forcefield_utilities as ffutils +import mbuild as mb +import pytest +from mbuild.lib.molecules import Ethane, Methane + +from gmso.external.convert_mbuild import from_mbuild +from gmso.parameterization.parameterize import apply +from gmso.tests.base_test import BaseTest + + +class TestSubTopologyParameterization(BaseTest): + @pytest.fixture(scope="session") + def ethane_methane_top(self): + cmpd = mb.Compound() + cmpd.add(Ethane()) + cmpd.add(Methane()) + gmso_top = from_mbuild(cmpd) + gmso_top.identify_connections() + return gmso_top + + def test_different_ffs_apply(self, ethane_methane_top): + opls = ffutils.FoyerFFs().load(ffname="oplsaa").to_gmso_ff() + apply(ethane_methane_top, {"Ethane": opls, "Methane": opls}) From 06df2ba2a5089310b91b93eb2be936096012ee3f Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Mon, 6 Jun 2022 12:22:17 -0500 Subject: [PATCH 18/32] Make opls/trappe global fixtures, Add tests for isomorphism --- gmso/tests/base_test.py | 13 +++++++++++++ gmso/tests/parameterization/test_opls_gmso.py | 15 ++++----------- ...ation.py => test_parameterization_options.py} | 16 +++++++++++++++- gmso/tests/parameterization/test_trappe_gmso.py | 8 -------- 4 files changed, 32 insertions(+), 20 deletions(-) rename gmso/tests/parameterization/{test_subtopology_parameterization.py => test_parameterization_options.py} (57%) diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index 68faa73e6..6e297a06c 100644 --- a/gmso/tests/base_test.py +++ b/gmso/tests/base_test.py @@ -4,6 +4,7 @@ import numpy as np import pytest import unyt as u +from forcefield_utilities.xml_loader import FoyerFFs from gmso.core.angle import Angle from gmso.core.atom import Atom @@ -545,3 +546,15 @@ def pentane_ua_parmed(self, pentane_ua_mbuild): @pytest.fixture(scope="session") def pentane_ua_gmso(self, pentane_ua_mbuild): return from_mbuild(pentane_ua_mbuild) + + @pytest.fixture(scope="session") + def xml_loader(self): + return FoyerFFs() + + @pytest.fixture(scope="session") + def oplsaa_gmso(self, xml_loader): + return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() + + @pytest.fixture(scope="session") + def trappe_ua_gmso(self, xml_loader): + return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py index 84e52dae4..897028e79 100644 --- a/gmso/tests/parameterization/test_opls_gmso.py +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -3,7 +3,6 @@ import parmed as pmd import pytest -from forcefield_utilities.xml_loader import FoyerFFs from pkg_resources import resource_filename from gmso.external.convert_parmed import from_parmed @@ -31,25 +30,19 @@ def get_foyer_opls_test_dirs(): @pytest.mark.skip class TestOPLSGMSO(BaseTest): - @pytest.fixture(scope="session") - def xml_loader(self): - return FoyerFFs() - - @pytest.fixture(scope="session") - def oplsaa(self, xml_loader): - return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() - @pytest.mark.parametrize( "system_dir", get_foyer_opls_test_dirs(), ids=lambda p: p.name ) def test_foyer_oplsaa_files( - self, system_dir, oplsaa, are_equivalent_topologies + self, system_dir, oplsaa_gmso, are_equivalent_topologies ): top_file = str(system_dir / f"{system_dir.name}.top") gro_file = str(system_dir / f"{system_dir.name}.gro") struct = pmd.load_file(top_file, xyz=gro_file) gmso_top_unparam = from_parmed(struct, refer_type=False) - apply(gmso_top_unparam, oplsaa, identify_connected_components=False) + apply( + gmso_top_unparam, oplsaa_gmso, identify_connected_components=False + ) assert are_equivalent_topologies( gmso_top_unparam, from_parmed(struct, refer_type=True) ), f"{system_dir} Failed" diff --git a/gmso/tests/parameterization/test_subtopology_parameterization.py b/gmso/tests/parameterization/test_parameterization_options.py similarity index 57% rename from gmso/tests/parameterization/test_subtopology_parameterization.py rename to gmso/tests/parameterization/test_parameterization_options.py index 4443490ca..e4d9db831 100644 --- a/gmso/tests/parameterization/test_subtopology_parameterization.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -8,7 +8,7 @@ from gmso.tests.base_test import BaseTest -class TestSubTopologyParameterization(BaseTest): +class TestParameterizationOptions(BaseTest): @pytest.fixture(scope="session") def ethane_methane_top(self): cmpd = mb.Compound() @@ -18,6 +18,20 @@ def ethane_methane_top(self): gmso_top.identify_connections() return gmso_top + @pytest.fixture(scope="session") + def ethane_box_with_methane(self): + cmpd_box = mb.fill_box([Ethane(), Methane()], [50, 50], density=1.0) + return from_mbuild(cmpd_box) + def test_different_ffs_apply(self, ethane_methane_top): opls = ffutils.FoyerFFs().load(ffname="oplsaa").to_gmso_ff() + ethane_methane_top.identify_connections() apply(ethane_methane_top, {"Ethane": opls, "Methane": opls}) + + def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): + apply( + ethane_box_with_methane, + oplsaa_gmso, + identify_connections=False, + identify_connected_components=True, + ) diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py index cb4a9287b..9c5b0bf19 100644 --- a/gmso/tests/parameterization/test_trappe_gmso.py +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -64,14 +64,6 @@ def _match_connection_parameters( class TestTrappeGMSO(BaseTest): - @pytest.fixture(scope="session") - def xml_loader(self): - return FoyerFFs() - - @pytest.fixture(scope="session") - def trappe_ua_gmso(self, xml_loader): - return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() - @pytest.fixture(scope="session") def trappe_ua_openmm_foyer(self): return foyer.forcefields.load_TRAPPE_UA() From d7b08d6df3d28b9a379d85af0de3dd565aac3e98 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Mon, 6 Jun 2022 12:35:59 -0500 Subject: [PATCH 19/32] Further testing isomorphism --- .../test_parameterization_options.py | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index e4d9db831..6b67bb994 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -1,3 +1,5 @@ +import random + import forcefield_utilities as ffutils import mbuild as mb import pytest @@ -29,9 +31,34 @@ def test_different_ffs_apply(self, ethane_methane_top): apply(ethane_methane_top, {"Ethane": opls, "Methane": opls}) def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): + ethane_box_with_methane.identify_connections() apply( ethane_box_with_methane, oplsaa_gmso, identify_connections=False, identify_connected_components=True, ) + + ethane_subtops = list( + filter( + lambda subtop: subtop.name == "Ethane", + ethane_box_with_methane.subtops, + ) + ) + methane_subtops = list( + filter( + lambda subtop: subtop.name == "Methane", + ethane_box_with_methane.subtops, + ) + ) + ethane_a = random.choice(ethane_subtops) + ethane_b = random.choice(ethane_subtops) + for atom_a, atom_b in zip(ethane_a.sites, ethane_b.sites): + assert atom_a.atom_type == atom_b.atom_type + assert atom_a.atom_type is not None + + methane_a = random.choice(methane_subtops) + methane_b = random.choice(methane_subtops) + for atom_a, atom_b in zip(methane_a.sites, methane_b.sites): + assert atom_a.atom_type == atom_b.atom_type + assert atom_a.atom_type is not None From 4a523e62831a14e45b28effeac33e0a334c46ca9 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Tue, 7 Jun 2022 10:19:58 -0500 Subject: [PATCH 20/32] REVERT - skip OPLS tests --- gmso/tests/parameterization/test_opls_gmso.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py index 897028e79..a2adbba93 100644 --- a/gmso/tests/parameterization/test_opls_gmso.py +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -28,7 +28,6 @@ def get_foyer_opls_test_dirs(): return parent_dirs -@pytest.mark.skip class TestOPLSGMSO(BaseTest): @pytest.mark.parametrize( "system_dir", get_foyer_opls_test_dirs(), ids=lambda p: p.name From 341ecc5cb8cec4f51c9aa17b73d18637dd787f4f Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Tue, 7 Jun 2022 13:06:03 -0500 Subject: [PATCH 21/32] Copy scaling factors and combining rules after parametrization --- .../topology_parameterizer.py | 31 ++++++++++++ .../test_parameterization_options.py | 48 +++++++++++++++++++ 2 files changed, 79 insertions(+) diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index cb0298d71..f32b07a28 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -181,8 +181,36 @@ def _parameterize(self, subtop_or_top, typemap, is_subtop=False): subtop_or_top, forcefield, is_subtop=is_subtop ) + def _verify_forcefields_metadata(self): + """Verify all the provided forcefields have the same scaling factors and combining rule.""" + if isinstance(self.forcefields, dict): + ffs = list(self.forcefields.values()) + init_scaling_factors = ffs[0].scaling_factors + init_combining_rule = ffs[0].combining_rule + for ff in ffs[1:]: + if ff.scaling_factors != init_scaling_factors: + raise GMSOParameterizationError( + "Scaling factors of the provided forcefields do not" + "match, please provide forcefields with same scaling" + "factors that apply to a Topology" + ) + + if ff.combining_rule != init_combining_rule: + raise GMSOParameterizationError( + "Combining rules of the provided forcefields do not" + "match, please provide forcefields with same scaling" + "factors that apply to a Topology" + ) + return init_scaling_factors, init_combining_rule + else: + return ( + self.forcefields.scaling_factors, + self.forcefields.combining_rule, + ) + def run_parameterization(self): """Run parameterization of the topology with give forcefield(s) and configuration.""" + scaling_factors, combining_rule = self._verify_forcefields_metadata() if self.topology.is_typed(): raise GMSOParameterizationError( "Cannot parameterize a typed topology. Please provide a topology without any types" @@ -225,6 +253,9 @@ def run_parameterization(self): typemap, is_subtop=False, # This will be removed from the future iterations ) + + self.topology.scaling_factors.update(scaling_factors) + self.topology.combining_rule = combining_rule self.topology.update_topology() @staticmethod diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index 6b67bb994..d87028f01 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -5,8 +5,12 @@ import pytest from mbuild.lib.molecules import Ethane, Methane +from gmso.core.forcefield import ForceField from gmso.external.convert_mbuild import from_mbuild from gmso.parameterization.parameterize import apply +from gmso.parameterization.topology_parameterizer import ( + GMSOParameterizationError, +) from gmso.tests.base_test import BaseTest @@ -25,10 +29,54 @@ def ethane_box_with_methane(self): cmpd_box = mb.fill_box([Ethane(), Methane()], [50, 50], density=1.0) return from_mbuild(cmpd_box) + def test_parameterization_error_different_scaling_factors( + self, ethane_methane_top + ): + ff1 = ForceField() + ff1.name = "FF1" + ff1.scaling_factors = { + "electrostatic14Scale": 1.0, + "columbic14Scale": 2.0, + } + ff2 = ForceField() + ff2.name = "FF2" + ff2.scaling_factors = { + "electrostatic14Scale": 3.0, + "columbic14Scale": 2.0, + } + + with pytest.raises(GMSOParameterizationError): + apply(ethane_methane_top, {"Ethane": ff1, "Methane": ff2}) + + def test_parameterization_different_combining_rule( + self, ethane_methane_top + ): + ff1 = ForceField() + ff1.name = "FF1" + ff1.scaling_factors = { + "electrostatic14Scale": 1.0, + "columbic14Scale": 1.0, + } + ff1.combining_rule = "lorrentz" + ff2 = ForceField() + ff2.name = "FF2" + ff2.scaling_factors = { + "electrostatic14Scale": 1.0, + "columbic14Scale": 1.0, + } + + ff2.combining_rule = "geometric" + + with pytest.raises(GMSOParameterizationError): + apply(ethane_methane_top, {"Ethane": ff1, "Methane": ff2}) + def test_different_ffs_apply(self, ethane_methane_top): opls = ffutils.FoyerFFs().load(ffname="oplsaa").to_gmso_ff() ethane_methane_top.identify_connections() apply(ethane_methane_top, {"Ethane": opls, "Methane": opls}) + assert ethane_methane_top.combining_rule == "geometric" + for key, v in opls.scaling_factors.items(): + assert ethane_methane_top.scaling_factors[key] == v def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): ethane_box_with_methane.identify_connections() From fe352fde1632dce5b4494db54fe1a42cd43a2bd3 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Tue, 7 Jun 2022 14:15:34 -0500 Subject: [PATCH 22/32] Proper OPLS tests --- gmso/tests/base_test.py | 4 ++ gmso/tests/parameterization/test_opls_gmso.py | 34 ++++++++++---- .../parameterization/test_trappe_gmso.py | 44 +++---------------- gmso/tests/utils.py | 35 +++++++++++++++ 4 files changed, 70 insertions(+), 47 deletions(-) diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index 6e297a06c..e0dce7930 100644 --- a/gmso/tests/base_test.py +++ b/gmso/tests/base_test.py @@ -558,3 +558,7 @@ def oplsaa_gmso(self, xml_loader): @pytest.fixture(scope="session") def trappe_ua_gmso(self, xml_loader): return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() + + @pytest.fixture(scope="session") + def oplsaa_foyer(self): + return foyer.forcefields.load_OPLSAA() diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py index a2adbba93..1d44c1525 100644 --- a/gmso/tests/parameterization/test_opls_gmso.py +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -3,11 +3,13 @@ import parmed as pmd import pytest +import unyt as u from pkg_resources import resource_filename from gmso.external.convert_parmed import from_parmed from gmso.parameterization.parameterize import apply from gmso.tests.base_test import BaseTest +from gmso.tests.utils import match_connection_parameters def get_foyer_opls_test_dirs(): @@ -25,7 +27,7 @@ def get_foyer_opls_test_dirs(): parent_dirs, ) ) - return parent_dirs + return parent_dirs[0:5] class TestOPLSGMSO(BaseTest): @@ -33,15 +35,29 @@ class TestOPLSGMSO(BaseTest): "system_dir", get_foyer_opls_test_dirs(), ids=lambda p: p.name ) def test_foyer_oplsaa_files( - self, system_dir, oplsaa_gmso, are_equivalent_topologies + self, system_dir, oplsaa_gmso, oplsaa_foyer, are_equivalent_topologies ): top_file = str(system_dir / f"{system_dir.name}.top") gro_file = str(system_dir / f"{system_dir.name}.gro") struct = pmd.load_file(top_file, xyz=gro_file) - gmso_top_unparam = from_parmed(struct, refer_type=False) - apply( - gmso_top_unparam, oplsaa_gmso, identify_connected_components=False - ) - assert are_equivalent_topologies( - gmso_top_unparam, from_parmed(struct, refer_type=True) - ), f"{system_dir} Failed" + gmso_top = from_parmed(struct, refer_type=False) + apply(gmso_top, oplsaa_gmso, identify_connected_components=False) + + pmd_struct_param = oplsaa_foyer.apply(struct) + gmso_top_from_pmd = from_parmed(pmd_struct_param, refer_type=True) + + for atom, mirror in zip(gmso_top.sites, gmso_top_from_pmd.sites): + assert atom.name == mirror.name + assert u.allclose_units(atom.mass, mirror.mass, 1e-3) + + assert u.allclose_units(atom.atom_type.charge, mirror.charge, 1e-3) + + atom_params = atom.atom_type.get_parameters() + mirror_params = mirror.atom_type.get_parameters() + + for k in atom_params: + assert u.allclose_units(atom_params[k], mirror_params[k]) + + match_connection_parameters(gmso_top, gmso_top_from_pmd) + match_connection_parameters(gmso_top, gmso_top_from_pmd, "angles") + match_connection_parameters(gmso_top, gmso_top_from_pmd, "dihedrals") diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py index 9c5b0bf19..a2197f28a 100644 --- a/gmso/tests/parameterization/test_trappe_gmso.py +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -4,13 +4,13 @@ import foyer import pytest import unyt as u -from forcefield_utilities.xml_loader import FoyerFFs from pkg_resources import resource_filename -import gmso +from gmso.core.topology import Topology from gmso.external.convert_parmed import from_parmed, to_parmed from gmso.parameterization.parameterize import apply from gmso.tests.base_test import BaseTest +from gmso.tests.utils import match_connection_parameters def get_foyer_trappe_test_dirs(): @@ -31,38 +31,6 @@ def get_foyer_trappe_test_dirs(): return parent_dirs -def _match_connection_parameters( - from_gmso_top, from_parmed_top, connection_type="bonds" -): - connection_types_original = {} - connection_types_mirror = {} - for connection in getattr(from_parmed_top, connection_type): - connection_types_mirror[ - tuple( - from_parmed_top.get_index(member) - for member in connection.connection_members - ) - ] = connection - - for connection in getattr(from_gmso_top, connection_type): - connection_types_original[ - tuple( - from_gmso_top.get_index(member) - for member in connection.connection_members - ) - ] = connection - - for key in connection_types_original: - conn = connection_types_original[key] - print(key) - conn_mirror = connection_types_mirror[key] - for param in conn_mirror.bond_type.parameters: - assert u.allclose_units( - conn_mirror.bond_type.parameters[param], - conn.bond_type.parameters[param], - ) - - class TestTrappeGMSO(BaseTest): @pytest.fixture(scope="session") def trappe_ua_openmm_foyer(self): @@ -79,7 +47,7 @@ def test_foyer_trappe_files( are_equivalent_connections, ): mol2_file = system_dir / f"{system_dir.name}.mol2" - gmso_top = gmso.Topology.load(mol2_file) + gmso_top = Topology.load(mol2_file) struct_pmd = trappe_ua_openmm_foyer.apply(to_parmed(gmso_top)) apply(gmso_top, trappe_ua_gmso, identify_connected_components=False) gmso_top_from_parmeterized_pmd = from_parmed(struct_pmd) @@ -96,10 +64,10 @@ def test_foyer_trappe_files( for k in atom_params: assert u.allclose_units(atom_params[k], mirror_params[k]) - _match_connection_parameters(gmso_top, gmso_top_from_parmeterized_pmd) - _match_connection_parameters( + match_connection_parameters(gmso_top, gmso_top_from_parmeterized_pmd) + match_connection_parameters( gmso_top, gmso_top_from_parmeterized_pmd, "angles" ) - _match_connection_parameters( + match_connection_parameters( gmso_top, gmso_top_from_parmeterized_pmd, "dihedrals" ) diff --git a/gmso/tests/utils.py b/gmso/tests/utils.py index 6a79e9adb..6d587308a 100644 --- a/gmso/tests/utils.py +++ b/gmso/tests/utils.py @@ -33,3 +33,38 @@ def allclose_units_mixed(u_iter1, u_iter2): if not u.allclose_units(q1, q2): return False return True + + +def match_connection_parameters( + from_gmso_top, from_parmed_top, connection_type="bonds" +): + """Match connection parameters between two gmso topologies.""" + connection_types_original = {} + connection_types_mirror = {} + for connection in getattr(from_parmed_top, connection_type): + connection_types_mirror[ + tuple( + from_parmed_top.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for connection in getattr(from_gmso_top, connection_type): + connection_types_original[ + tuple( + from_gmso_top.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for key in connection_types_original: + conn = connection_types_original[key] + conn_mirror = connection_types_mirror[key] + conn_type_attr = connection_type[:-1] + "_type" + conn_type_mirror = getattr(conn_mirror, conn_type_attr) + conn_type = getattr(conn, conn_type_attr) + for param in conn_type.parameters: + assert u.allclose_units( + conn_type_mirror.parameters[param], + conn_type.parameters[param], + ) From b3e045f086c6ba28916fc2464edc8318fa1cb732 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Tue, 7 Jun 2022 18:08:32 -0500 Subject: [PATCH 23/32] WIP- Refactor the test module --- gmso/tests/base_test.py | 17 ---- .../parameterization_base_test.py | 85 +++++++++++++++++++ gmso/tests/parameterization/test_opls_gmso.py | 43 ++++------ .../test_parameterization_options.py | 6 +- .../test_subtopology_utils.py | 6 +- .../parameterization/test_trappe_gmso.py | 37 +++----- gmso/tests/utils.py | 35 -------- 7 files changed, 123 insertions(+), 106 deletions(-) create mode 100644 gmso/tests/parameterization/parameterization_base_test.py diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index e0dce7930..68faa73e6 100644 --- a/gmso/tests/base_test.py +++ b/gmso/tests/base_test.py @@ -4,7 +4,6 @@ import numpy as np import pytest import unyt as u -from forcefield_utilities.xml_loader import FoyerFFs from gmso.core.angle import Angle from gmso.core.atom import Atom @@ -546,19 +545,3 @@ def pentane_ua_parmed(self, pentane_ua_mbuild): @pytest.fixture(scope="session") def pentane_ua_gmso(self, pentane_ua_mbuild): return from_mbuild(pentane_ua_mbuild) - - @pytest.fixture(scope="session") - def xml_loader(self): - return FoyerFFs() - - @pytest.fixture(scope="session") - def oplsaa_gmso(self, xml_loader): - return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() - - @pytest.fixture(scope="session") - def trappe_ua_gmso(self, xml_loader): - return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() - - @pytest.fixture(scope="session") - def oplsaa_foyer(self): - return foyer.forcefields.load_OPLSAA() diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py new file mode 100644 index 000000000..b0fd43b5a --- /dev/null +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -0,0 +1,85 @@ +import foyer +import pytest +import unyt as u +from forcefield_utilities.xml_loader import FoyerFFs + +from gmso.tests.base_test import BaseTest + + +class ParameterizationBaseTest(BaseTest): + @pytest.fixture(scope="session") + def xml_loader(self): + return FoyerFFs() + + @pytest.fixture(scope="session") + def oplsaa_gmso(self, xml_loader): + return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() + + @pytest.fixture(scope="session") + def trappe_ua_gmso(self, xml_loader): + return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() + + @pytest.fixture(scope="session") + def oplsaa_foyer(self): + return foyer.forcefields.load_OPLSAA() + + @pytest.fixture(scope="session") + def trappe_ua_foyer(self): + return foyer.forcefields.load_TRAPPE_UA() + + @pytest.fixture(scope="session") + def assert_same_connection_params(self): + def _assert_same_connection_params(top1, top2, connection_type="bonds"): + """Match connection parameters between two gmso topologies.""" + connection_types_original = {} + connection_types_mirror = {} + for connection in getattr(top2, connection_type): + connection_types_mirror[ + tuple( + top2.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for connection in getattr(top1, connection_type): + connection_types_original[ + tuple( + top1.get_index(member) + for member in connection.connection_members + ) + ] = connection + + for key in connection_types_original: + conn = connection_types_original[key] + conn_mirror = connection_types_mirror[key] + conn_type_attr = connection_type[:-1] + "_type" + conn_type_mirror = getattr(conn_mirror, conn_type_attr) + conn_type = getattr(conn, conn_type_attr) + for param in conn_type.parameters: + assert u.allclose_units( + conn_type_mirror.parameters[param], + conn_type.parameters[param], + ) + + return _assert_same_connection_params + + @pytest.fixture(scope="session") + def assert_same_atom_params(self): + def _assert_same_atom_params(top1, top2): + """Match atom parameters between two gmso topologies. + + Notes + ----- + This is specific + """ + for atom, mirror in zip(top1.sites, top2.sites): + assert atom.name == mirror.name + assert u.allclose_units(atom.mass, mirror.mass, 1e-3) + + atom_params = atom.atom_type.get_parameters() + mirror_params = mirror.atom_type.get_parameters() + + for k in atom_params: + assert u.allclose_units(atom_params[k], mirror_params[k]) + + return _assert_same_atom_params diff --git a/gmso/tests/parameterization/test_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py index 1d44c1525..c5ed7639f 100644 --- a/gmso/tests/parameterization/test_opls_gmso.py +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -3,13 +3,13 @@ import parmed as pmd import pytest -import unyt as u from pkg_resources import resource_filename from gmso.external.convert_parmed import from_parmed from gmso.parameterization.parameterize import apply -from gmso.tests.base_test import BaseTest -from gmso.tests.utils import match_connection_parameters +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) def get_foyer_opls_test_dirs(): @@ -27,37 +27,30 @@ def get_foyer_opls_test_dirs(): parent_dirs, ) ) - return parent_dirs[0:5] + return parent_dirs -class TestOPLSGMSO(BaseTest): +class TestOPLSGMSO(ParameterizationBaseTest): @pytest.mark.parametrize( "system_dir", get_foyer_opls_test_dirs(), ids=lambda p: p.name ) def test_foyer_oplsaa_files( - self, system_dir, oplsaa_gmso, oplsaa_foyer, are_equivalent_topologies + self, + system_dir, + oplsaa_gmso, + oplsaa_foyer, + assert_same_connection_params, + assert_same_atom_params, ): top_file = str(system_dir / f"{system_dir.name}.top") gro_file = str(system_dir / f"{system_dir.name}.gro") - struct = pmd.load_file(top_file, xyz=gro_file) + struct = oplsaa_foyer.apply(pmd.load_file(top_file, xyz=gro_file)) + + gmso_top_from_pmd = from_parmed(struct, refer_type=True) gmso_top = from_parmed(struct, refer_type=False) apply(gmso_top, oplsaa_gmso, identify_connected_components=False) - pmd_struct_param = oplsaa_foyer.apply(struct) - gmso_top_from_pmd = from_parmed(pmd_struct_param, refer_type=True) - - for atom, mirror in zip(gmso_top.sites, gmso_top_from_pmd.sites): - assert atom.name == mirror.name - assert u.allclose_units(atom.mass, mirror.mass, 1e-3) - - assert u.allclose_units(atom.atom_type.charge, mirror.charge, 1e-3) - - atom_params = atom.atom_type.get_parameters() - mirror_params = mirror.atom_type.get_parameters() - - for k in atom_params: - assert u.allclose_units(atom_params[k], mirror_params[k]) - - match_connection_parameters(gmso_top, gmso_top_from_pmd) - match_connection_parameters(gmso_top, gmso_top_from_pmd, "angles") - match_connection_parameters(gmso_top, gmso_top_from_pmd, "dihedrals") + assert_same_atom_params(gmso_top, gmso_top_from_pmd) + assert_same_connection_params(gmso_top, gmso_top_from_pmd) + assert_same_connection_params(gmso_top, gmso_top_from_pmd, "angles") + assert_same_connection_params(gmso_top, gmso_top_from_pmd, "dihedrals") diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index d87028f01..8add09374 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -11,10 +11,12 @@ from gmso.parameterization.topology_parameterizer import ( GMSOParameterizationError, ) -from gmso.tests.base_test import BaseTest +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) -class TestParameterizationOptions(BaseTest): +class TestParameterizationOptions(ParameterizationBaseTest): @pytest.fixture(scope="session") def ethane_methane_top(self): cmpd = mb.Compound() diff --git a/gmso/tests/parameterization/test_subtopology_utils.py b/gmso/tests/parameterization/test_subtopology_utils.py index b6aba7049..32c4b7690 100644 --- a/gmso/tests/parameterization/test_subtopology_utils.py +++ b/gmso/tests/parameterization/test_subtopology_utils.py @@ -10,11 +10,13 @@ subtop_dihedrals, subtop_impropers, ) -from gmso.tests.base_test import BaseTest +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) from gmso.utils.connectivity import identify_connections -class TestSubTopologyUtils(BaseTest): +class TestSubTopologyUtils(ParameterizationBaseTest): @pytest.fixture(scope="session") def ethane_box_gmso(self): ethane_box = mb.fill_box( diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py index a2197f28a..02dd958ed 100644 --- a/gmso/tests/parameterization/test_trappe_gmso.py +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -9,8 +9,9 @@ from gmso.core.topology import Topology from gmso.external.convert_parmed import from_parmed, to_parmed from gmso.parameterization.parameterize import apply -from gmso.tests.base_test import BaseTest -from gmso.tests.utils import match_connection_parameters +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) def get_foyer_trappe_test_dirs(): @@ -31,43 +32,29 @@ def get_foyer_trappe_test_dirs(): return parent_dirs -class TestTrappeGMSO(BaseTest): - @pytest.fixture(scope="session") - def trappe_ua_openmm_foyer(self): - return foyer.forcefields.load_TRAPPE_UA() - +class TestTrappeGMSO(ParameterizationBaseTest): @pytest.mark.parametrize( "system_dir", get_foyer_trappe_test_dirs(), ids=lambda p: p.name ) def test_foyer_trappe_files( self, system_dir, - trappe_ua_openmm_foyer, + trappe_ua_foyer, trappe_ua_gmso, - are_equivalent_connections, + assert_same_connection_params, + assert_same_atom_params, ): mol2_file = system_dir / f"{system_dir.name}.mol2" gmso_top = Topology.load(mol2_file) - struct_pmd = trappe_ua_openmm_foyer.apply(to_parmed(gmso_top)) + struct_pmd = trappe_ua_foyer.apply(to_parmed(gmso_top)) apply(gmso_top, trappe_ua_gmso, identify_connected_components=False) gmso_top_from_parmeterized_pmd = from_parmed(struct_pmd) - for atom, mirror in zip( - gmso_top.sites, gmso_top_from_parmeterized_pmd.sites - ): - assert atom.name == mirror.name - assert u.allclose_units(atom.mass, mirror.mass, 1e-4) - - assert u.allclose_units(atom.atom_type.charge, mirror.charge, 1e-4) - - atom_params = atom.atom_type.get_parameters() - mirror_params = mirror.atom_type.get_parameters() - for k in atom_params: - assert u.allclose_units(atom_params[k], mirror_params[k]) - match_connection_parameters(gmso_top, gmso_top_from_parmeterized_pmd) - match_connection_parameters( + assert_same_atom_params(gmso_top_from_parmeterized_pmd, gmso_top) + assert_same_connection_params(gmso_top, gmso_top_from_parmeterized_pmd) + assert_same_connection_params( gmso_top, gmso_top_from_parmeterized_pmd, "angles" ) - match_connection_parameters( + assert_same_connection_params( gmso_top, gmso_top_from_parmeterized_pmd, "dihedrals" ) diff --git a/gmso/tests/utils.py b/gmso/tests/utils.py index 6d587308a..6a79e9adb 100644 --- a/gmso/tests/utils.py +++ b/gmso/tests/utils.py @@ -33,38 +33,3 @@ def allclose_units_mixed(u_iter1, u_iter2): if not u.allclose_units(q1, q2): return False return True - - -def match_connection_parameters( - from_gmso_top, from_parmed_top, connection_type="bonds" -): - """Match connection parameters between two gmso topologies.""" - connection_types_original = {} - connection_types_mirror = {} - for connection in getattr(from_parmed_top, connection_type): - connection_types_mirror[ - tuple( - from_parmed_top.get_index(member) - for member in connection.connection_members - ) - ] = connection - - for connection in getattr(from_gmso_top, connection_type): - connection_types_original[ - tuple( - from_gmso_top.get_index(member) - for member in connection.connection_members - ) - ] = connection - - for key in connection_types_original: - conn = connection_types_original[key] - conn_mirror = connection_types_mirror[key] - conn_type_attr = connection_type[:-1] + "_type" - conn_type_mirror = getattr(conn_mirror, conn_type_attr) - conn_type = getattr(conn, conn_type_attr) - for param in conn_type.parameters: - assert u.allclose_units( - conn_type_mirror.parameters[param], - conn_type.parameters[param], - ) From 4834b5f1a4ccc60249522a9939a62a1c7fcda2b0 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Tue, 7 Jun 2022 18:09:24 -0500 Subject: [PATCH 24/32] WIP- Remove unused import --- gmso/tests/base_test.py | 2 -- gmso/tests/parameterization/test_trappe_gmso.py | 2 -- 2 files changed, 4 deletions(-) diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index 68faa73e6..bc6bd074d 100644 --- a/gmso/tests/base_test.py +++ b/gmso/tests/base_test.py @@ -1,6 +1,5 @@ import foyer import mbuild as mb -import mbuild.recipes import numpy as np import pytest import unyt as u @@ -11,7 +10,6 @@ from gmso.core.bond import Bond from gmso.core.box import Box from gmso.core.dihedral import Dihedral -from gmso.core.element import Hydrogen, Oxygen from gmso.core.forcefield import ForceField from gmso.core.improper import Improper from gmso.core.pairpotential_type import PairPotentialType diff --git a/gmso/tests/parameterization/test_trappe_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py index 02dd958ed..84e489e40 100644 --- a/gmso/tests/parameterization/test_trappe_gmso.py +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -1,9 +1,7 @@ import glob from pathlib import Path -import foyer import pytest -import unyt as u from pkg_resources import resource_filename from gmso.core.topology import Topology From 0a4c78042f73325d3039658bbb8d1c5035486253 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Thu, 9 Jun 2022 13:38:54 -0500 Subject: [PATCH 25/32] WIP- Add test for parameterization with impropers --- gmso/tests/files/fake_ethane_impropers.xml | 28 ++++++++++ .../parameterization_base_test.py | 7 +++ .../test_parameterization_options.py | 51 +++++++++++++++++++ 3 files changed, 86 insertions(+) create mode 100644 gmso/tests/files/fake_ethane_impropers.xml diff --git a/gmso/tests/files/fake_ethane_impropers.xml b/gmso/tests/files/fake_ethane_impropers.xml new file mode 100644 index 000000000..98d063ebb --- /dev/null +++ b/gmso/tests/files/fake_ethane_impropers.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py index b0fd43b5a..d31de0d71 100644 --- a/gmso/tests/parameterization/parameterization_base_test.py +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -4,6 +4,7 @@ from forcefield_utilities.xml_loader import FoyerFFs from gmso.tests.base_test import BaseTest +from gmso.tests.utils import get_path class ParameterizationBaseTest(BaseTest): @@ -19,6 +20,12 @@ def oplsaa_gmso(self, xml_loader): def trappe_ua_gmso(self, xml_loader): return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() + @pytest.fixture(scope="session") + def fake_improper_ff_gmso(self, xml_loader): + return xml_loader.load( + get_path("fake_ethane_impropers.xml"), rel_to_module=True + ).to_gmso_ff() + @pytest.fixture(scope="session") def oplsaa_foyer(self): return foyer.forcefields.load_OPLSAA() diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index 8add09374..3ea5c7942 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -3,10 +3,13 @@ import forcefield_utilities as ffutils import mbuild as mb import pytest +import unyt as u from mbuild.lib.molecules import Ethane, Methane from gmso.core.forcefield import ForceField +from gmso.core.views import PotentialFilters from gmso.external.convert_mbuild import from_mbuild +from gmso.lib.potential_templates import PotentialTemplateLibrary from gmso.parameterization.parameterize import apply from gmso.parameterization.topology_parameterizer import ( GMSOParameterizationError, @@ -112,3 +115,51 @@ def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): for atom_a, atom_b in zip(methane_a.sites, methane_b.sites): assert atom_a.atom_type == atom_b.atom_type assert atom_a.atom_type is not None + + def test_improper_parameterization(self, fake_improper_ff_gmso, ethane): + ethane.identify_connections() + apply(ethane, fake_improper_ff_gmso, assert_improper_params=True) + + lib = PotentialTemplateLibrary() + template_improper_type = lib["PeriodicImproperPotential"] + + assert ( + len( + ethane.improper_types( + filter_by=PotentialFilters.UNIQUE_NAME_CLASS + ) + ) + == 2 + ) + for improper_type in ethane.improper_types: + assert improper_type.expression == template_improper_type.expression + assert improper_type.member_classes in { + ("CT", "CT", "HC", "HC"), + ("CT", "HC", "HC", "HC"), + } + + if improper_type.member_classes == ("CT", "CT", "HC", "HC"): + assert u.allclose_units( + improper_type.parameters["phi_eq"], [180.0] * u.degree + ) + assert u.allclose_units( + improper_type.parameters["k"], [4.6024] * u.kJ / u.mol + ) + assert u.allclose_units( + improper_type.parameters["n"], [2] * u.dimensionless + ) + + elif improper_type.member_classes == ("CT", "CT", "HC", "HC"): + assert u.allclose_units( + improper_type.parameters["phi_eq"], [180.0] * u.degree + ) + assert u.allclose_units( + improper_type.parameters["k"], [2.5560] * u.kJ / u.mol + ) + assert u.allclose_units( + improper_type.parameters["n"], [2] * u.dimensionless + ) + + def test_improper_assertion_error(self, ethane_methane_top, oplsaa_gmso): + with pytest.raises(GMSOParameterizationError): + apply(ethane_methane_top, oplsaa_gmso) From 3711172c6b05ba3e1b06037870c46309a5c94185 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 15 Jun 2022 13:19:10 -0500 Subject: [PATCH 26/32] WIP- Additional impropers test; Separate module for testing impropers --- gmso/tests/files/benzene_aa.mol2 | 38 +++++++ ...benzene_and_alkane_branched_benzene_aa.xml | 59 ++++++++++ gmso/tests/files/ethyl_benzene_aa.mol2 | 50 +++++++++ gmso/tests/files/methyl_benzene_aa.mol2 | 44 ++++++++ .../parameterization_base_test.py | 24 ++++ .../test_impropers_parameterization.py | 106 ++++++++++++++++++ .../test_parameterization_options.py | 68 ----------- 7 files changed, 321 insertions(+), 68 deletions(-) create mode 100644 gmso/tests/files/benzene_aa.mol2 create mode 100644 gmso/tests/files/benzene_and_alkane_branched_benzene_aa.xml create mode 100644 gmso/tests/files/ethyl_benzene_aa.mol2 create mode 100644 gmso/tests/files/methyl_benzene_aa.mol2 create mode 100644 gmso/tests/parameterization/test_impropers_parameterization.py diff --git a/gmso/tests/files/benzene_aa.mol2 b/gmso/tests/files/benzene_aa.mol2 new file mode 100644 index 000000000..26c3fe27b --- /dev/null +++ b/gmso/tests/files/benzene_aa.mol2 @@ -0,0 +1,38 @@ +@MOLECULE +BEN + 12 12 1 0 0 +SMALL +NO_CHARGES +**** +Energy = 0 + +@ATOM + 1 C 0.0000 0.0000 0.0000 C 1 BEN 0.000000 + 2 C 1.4000 0.0000 0.0000 C 1 BEN 0.000000 + 3 C 2.1000 1.2124 0.0000 C 1 BEN 0.000000 + 4 C 1.4000 2.4249 0.0000 C 1 BEN 0.000000 + 5 C 0.0000 2.4249 0.0000 C 1 BEN 0.000000 + 6 C -0.7000 1.2124 0.0000 C 1 BEN 0.000000 + 7 H -0.5000 -0.8660 0.0000 H 1 BEN 0.000000 + 8 H 1.9000 -0.8660 0.0000 H 1 BEN 0.000000 + 9 H 3.1000 1.2124 0.0000 H 1 BEN 0.000000 + 10 H 1.9000 3.2909 0.0000 H 1 BEN 0.000000 + 11 H -0.5000 3.2909 0.0000 H 1 BEN 0.000000 + 12 H -1.7000 1.2124 0.0000 H 1 BEN 0.000000 +@BOND + 1 1 2 2 + 2 1 6 1 + 3 1 7 1 + 4 2 3 1 + 5 2 8 1 + 6 3 4 2 + 7 3 9 1 + 8 4 5 1 + 9 4 10 1 + 10 5 6 2 + 11 5 11 1 + 12 6 12 1 +@SUBSTRUCTURE +1 **** 1 TEMP 0 **** **** 0 ROOT + +#generated by VMD diff --git a/gmso/tests/files/benzene_and_alkane_branched_benzene_aa.xml b/gmso/tests/files/benzene_and_alkane_branched_benzene_aa.xml new file mode 100644 index 000000000..7c3b6496f --- /dev/null +++ b/gmso/tests/files/benzene_and_alkane_branched_benzene_aa.xml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/gmso/tests/files/ethyl_benzene_aa.mol2 b/gmso/tests/files/ethyl_benzene_aa.mol2 new file mode 100644 index 000000000..b46f2354f --- /dev/null +++ b/gmso/tests/files/ethyl_benzene_aa.mol2 @@ -0,0 +1,50 @@ +@MOLECULE +EBN + 18 18 1 0 0 +SMALL +NO_CHARGES +**** +Energy = 0 + +@ATOM + 1 C 0.0000 0.0000 0.0000 C 1 EBN 0.000000 + 2 C 1.4000 0.0000 0.0000 C 1 EBN 0.000000 + 3 C 2.1000 1.2124 0.0000 C 1 EBN 0.000000 + 4 C 1.4000 2.4249 0.0000 C 1 EBN 0.000000 + 5 C 0.0000 2.4249 0.0000 C 1 EBN 0.000000 + 6 C -0.7000 1.2124 0.0000 C 1 EBN 0.000000 + 7 C -2.2000 1.2123 0.0000 C 1 EBN 0.000000 + 8 C -2.9501 2.5113 0.0000 C 1 EBN 0.000000 + 9 H -0.5000 -0.8660 0.0000 H 1 EBN 0.000000 + 10 H 1.9000 -0.8660 0.0000 H 1 EBN 0.000000 + 11 H 3.1000 1.2123 0.0000 H 1 EBN 0.000000 + 12 H 1.9000 3.2909 0.0000 H 1 EBN 0.000000 + 13 H -0.5000 3.2909 0.0000 H 1 EBN 0.000000 + 14 H -2.4903 0.7094 -0.8141 H 1 EBN 0.000000 + 15 H -2.4903 0.7094 0.8141 H 1 EBN 0.000000 + 16 H -2.2941 3.2660 -0.0000 H 1 EBN 0.000000 + 17 H -3.5223 2.5568 -0.8188 H 1 EBN 0.000000 + 18 H -3.5223 2.5568 0.8188 H 1 EBN 0.000000 +@BOND + 1 1 2 1 + 2 1 6 2 + 3 1 9 1 + 4 2 3 2 + 5 2 10 1 + 6 3 4 1 + 7 3 11 1 + 8 4 5 2 + 9 4 12 1 + 10 5 6 1 + 11 5 13 1 + 12 6 7 1 + 13 7 8 1 + 14 7 14 1 + 15 7 15 1 + 16 8 16 1 + 17 8 17 1 + 18 8 18 1 +@SUBSTRUCTURE +1 **** 1 TEMP 0 **** **** 0 ROOT + +#generated by VMD diff --git a/gmso/tests/files/methyl_benzene_aa.mol2 b/gmso/tests/files/methyl_benzene_aa.mol2 new file mode 100644 index 000000000..670c1fbcd --- /dev/null +++ b/gmso/tests/files/methyl_benzene_aa.mol2 @@ -0,0 +1,44 @@ +@MOLECULE +MBN + 15 15 1 0 0 +SMALL +NO_CHARGES +**** +Energy = 0 + +@ATOM + 1 C 0.0000 0.0000 0.0000 C 1 MBN 0.000000 + 2 C 1.4000 0.0000 0.0000 C 1 MBN 0.000000 + 3 C 2.1000 1.2124 0.0000 C 1 MBN 0.000000 + 4 C 1.4000 2.4249 0.0000 C 1 MBN 0.000000 + 5 C 0.0000 2.4249 0.0000 C 1 MBN 0.000000 + 6 C -0.7000 1.2124 0.0000 C 1 MBN 0.000000 + 7 C -2.2000 1.2123 0.0000 C 1 MBN 0.000000 + 8 H -0.5000 -0.8660 0.0000 H 1 MBN 0.000000 + 9 H 1.9000 -0.8660 0.0000 H 1 MBN 0.000000 + 10 H 3.1000 1.2123 0.0000 H 1 MBN 0.000000 + 11 H 1.9000 3.2909 0.0000 H 1 MBN 0.000000 + 12 H -0.5000 3.2909 0.0000 H 1 MBN 0.000000 + 13 H -2.5255 0.2668 0.0000 H 1 MBN 0.000000 + 14 H -2.5256 1.6850 0.8188 H 1 MBN 0.000000 + 15 H -2.5256 1.6850 -0.8188 H 1 MBN 0.000000 +@BOND + 1 1 2 1 + 2 1 6 2 + 3 1 8 1 + 4 2 3 2 + 5 2 9 1 + 6 3 4 1 + 7 3 10 1 + 8 4 5 2 + 9 4 11 1 + 10 5 6 1 + 11 5 12 1 + 12 6 7 1 + 13 7 13 1 + 14 7 14 1 + 15 7 15 1 +@SUBSTRUCTURE +1 **** 1 TEMP 0 **** **** 0 ROOT + +#generated by VMD diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py index d31de0d71..dea30e41c 100644 --- a/gmso/tests/parameterization/parameterization_base_test.py +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -1,8 +1,11 @@ import foyer +import mbuild as mb import pytest import unyt as u from forcefield_utilities.xml_loader import FoyerFFs +from mbuild.lib.molecules import Ethane, Methane +from gmso.external.convert_mbuild import from_mbuild from gmso.tests.base_test import BaseTest from gmso.tests.utils import get_path @@ -26,6 +29,13 @@ def fake_improper_ff_gmso(self, xml_loader): get_path("fake_ethane_impropers.xml"), rel_to_module=True ).to_gmso_ff() + @pytest.fixture(scope="session") + def benzene_alkane_aa_ff_gmso(self, xml_loader): + return xml_loader.load( + get_path("benzene_and_alkane_branched_benzene_aa.xml"), + rel_to_module=True, + ).to_gmso_ff() + @pytest.fixture(scope="session") def oplsaa_foyer(self): return foyer.forcefields.load_OPLSAA() @@ -90,3 +100,17 @@ def _assert_same_atom_params(top1, top2): assert u.allclose_units(atom_params[k], mirror_params[k]) return _assert_same_atom_params + + @pytest.fixture(scope="session") + def ethane_methane_top(self): + cmpd = mb.Compound() + cmpd.add(Ethane()) + cmpd.add(Methane()) + gmso_top = from_mbuild(cmpd) + gmso_top.identify_connections() + return gmso_top + + @pytest.fixture(scope="session") + def ethane_box_with_methane(self): + cmpd_box = mb.fill_box([Ethane(), Methane()], [50, 50], density=1.0) + return from_mbuild(cmpd_box) diff --git a/gmso/tests/parameterization/test_impropers_parameterization.py b/gmso/tests/parameterization/test_impropers_parameterization.py new file mode 100644 index 000000000..ef786a7e9 --- /dev/null +++ b/gmso/tests/parameterization/test_impropers_parameterization.py @@ -0,0 +1,106 @@ +from pathlib import Path + +import pytest +import unyt as u + +from gmso.core.topology import Topology +from gmso.core.views import PotentialFilters +from gmso.lib.potential_templates import PotentialTemplateLibrary +from gmso.parameterization.parameterize import apply +from gmso.parameterization.topology_parameterizer import ( + GMSOParameterizationError, +) +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) +from gmso.tests.utils import get_path + + +class TestImpropersParameterization(ParameterizationBaseTest): + def test_improper_parameterization(self, fake_improper_ff_gmso, ethane): + ethane.identify_connections() + apply(ethane, fake_improper_ff_gmso, assert_improper_params=True) + + lib = PotentialTemplateLibrary() + template_improper_type = lib["PeriodicImproperPotential"] + + assert ( + len( + ethane.improper_types( + filter_by=PotentialFilters.UNIQUE_NAME_CLASS + ) + ) + == 2 + ) + for improper_type in ethane.improper_types: + assert improper_type.expression == template_improper_type.expression + assert improper_type.member_classes in { + ("CT", "CT", "HC", "HC"), + ("CT", "HC", "HC", "HC"), + } + + if improper_type.member_classes == ("CT", "CT", "HC", "HC"): + assert u.allclose_units( + improper_type.parameters["phi_eq"], [180.0] * u.degree + ) + assert u.allclose_units( + improper_type.parameters["k"], [4.6024] * u.kJ / u.mol + ) + assert u.allclose_units( + improper_type.parameters["n"], [2] * u.dimensionless + ) + + elif improper_type.member_classes == ("CT", "CT", "HC", "HC"): + assert u.allclose_units( + improper_type.parameters["phi_eq"], [180.0] * u.degree + ) + assert u.allclose_units( + improper_type.parameters["k"], [2.5560] * u.kJ / u.mol + ) + assert u.allclose_units( + improper_type.parameters["n"], [2] * u.dimensionless + ) + + def test_improper_assertion_error(self, ethane_methane_top, oplsaa_gmso): + with pytest.raises(GMSOParameterizationError): + apply(ethane_methane_top, oplsaa_gmso, assert_improper_params=True) + + @pytest.mark.parametrize( + "mol2_loc", + [ + get_path("methyl_benzene_aa.mol2"), + get_path("benzene_aa.mol2"), + get_path("ethyl_benzene_aa.mol2"), + ], + ids=lambda p: Path(p).stem, + ) + def test_benzene_aa_ff(self, mol2_loc, benzene_alkane_aa_ff_gmso): + gmso_top = Topology.load(filename=mol2_loc) + apply(gmso_top, benzene_alkane_aa_ff_gmso, identify_connections=True) + for improper in gmso_top.impropers: + if improper.improper_type: + if [ + improper.member_types[0], + set(improper.member_types[1:]), + ] == ["CH_sp2", {"HCE", "CH_sp2", "C_sp2"}]: + params = improper.improper_type.get_parameters() + assert u.allclose_units(params["k"], 4.0 * u.kJ / u.mol) + assert u.allclose_units(params["n"], 2.0 * u.dimensionless) + assert u.allclose_units(params["phi_eq"], 0 * u.radian) + + elif [ + improper.member_types[0], + set(improper.member_types[1:]), + ] == ["C_sp2", {"CH3_sp3", "CH_sp2", "CH_sp2"}]: + params = improper.improper_type.get_parameters() + assert u.allclose_units(params["k"], 4.0 * u.kJ / u.mol) + assert u.allclose_units(params["n"], 1.0 * u.dimensionless) + assert u.allclose_units( + params["phi_eq"], 180 * u.degree, atol=10e-4 + ) + else: + assert len(improper.member_classes) == 4 + assert set(improper.member_classes) not in [ + {"CE", "HCE", "CE", "CE"}, + {"CE", "CT", "CE", "CE"}, + ] diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index 3ea5c7942..1d529db04 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -1,15 +1,9 @@ import random import forcefield_utilities as ffutils -import mbuild as mb import pytest -import unyt as u -from mbuild.lib.molecules import Ethane, Methane from gmso.core.forcefield import ForceField -from gmso.core.views import PotentialFilters -from gmso.external.convert_mbuild import from_mbuild -from gmso.lib.potential_templates import PotentialTemplateLibrary from gmso.parameterization.parameterize import apply from gmso.parameterization.topology_parameterizer import ( GMSOParameterizationError, @@ -20,20 +14,6 @@ class TestParameterizationOptions(ParameterizationBaseTest): - @pytest.fixture(scope="session") - def ethane_methane_top(self): - cmpd = mb.Compound() - cmpd.add(Ethane()) - cmpd.add(Methane()) - gmso_top = from_mbuild(cmpd) - gmso_top.identify_connections() - return gmso_top - - @pytest.fixture(scope="session") - def ethane_box_with_methane(self): - cmpd_box = mb.fill_box([Ethane(), Methane()], [50, 50], density=1.0) - return from_mbuild(cmpd_box) - def test_parameterization_error_different_scaling_factors( self, ethane_methane_top ): @@ -115,51 +95,3 @@ def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): for atom_a, atom_b in zip(methane_a.sites, methane_b.sites): assert atom_a.atom_type == atom_b.atom_type assert atom_a.atom_type is not None - - def test_improper_parameterization(self, fake_improper_ff_gmso, ethane): - ethane.identify_connections() - apply(ethane, fake_improper_ff_gmso, assert_improper_params=True) - - lib = PotentialTemplateLibrary() - template_improper_type = lib["PeriodicImproperPotential"] - - assert ( - len( - ethane.improper_types( - filter_by=PotentialFilters.UNIQUE_NAME_CLASS - ) - ) - == 2 - ) - for improper_type in ethane.improper_types: - assert improper_type.expression == template_improper_type.expression - assert improper_type.member_classes in { - ("CT", "CT", "HC", "HC"), - ("CT", "HC", "HC", "HC"), - } - - if improper_type.member_classes == ("CT", "CT", "HC", "HC"): - assert u.allclose_units( - improper_type.parameters["phi_eq"], [180.0] * u.degree - ) - assert u.allclose_units( - improper_type.parameters["k"], [4.6024] * u.kJ / u.mol - ) - assert u.allclose_units( - improper_type.parameters["n"], [2] * u.dimensionless - ) - - elif improper_type.member_classes == ("CT", "CT", "HC", "HC"): - assert u.allclose_units( - improper_type.parameters["phi_eq"], [180.0] * u.degree - ) - assert u.allclose_units( - improper_type.parameters["k"], [2.5560] * u.kJ / u.mol - ) - assert u.allclose_units( - improper_type.parameters["n"], [2] * u.dimensionless - ) - - def test_improper_assertion_error(self, ethane_methane_top, oplsaa_gmso): - with pytest.raises(GMSOParameterizationError): - apply(ethane_methane_top, oplsaa_gmso) From 956977c5249587f6039ade1196449ab049a2d6cc Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 15 Jun 2022 14:33:25 -0500 Subject: [PATCH 27/32] Minor refacotors; additional edge cases coverage/tests --- gmso/parameterization/foyer_utils.py | 7 +++ .../topology_parameterizer.py | 18 +++++--- .../parameterization_base_test.py | 4 +- .../test_impropers_parameterization.py | 6 +-- .../test_parameterization_options.py | 43 ++++++++++++++++--- 5 files changed, 62 insertions(+), 16 deletions(-) diff --git a/gmso/parameterization/foyer_utils.py b/gmso/parameterization/foyer_utils.py index 3f25b3610..25bf56575 100644 --- a/gmso/parameterization/foyer_utils.py +++ b/gmso/parameterization/foyer_utils.py @@ -2,6 +2,7 @@ from collections import namedtuple from foyer.atomtyper import AtomTypingRulesProvider, find_atomtypes +from foyer.exceptions import FoyerError from foyer.topology_graph import TopologyGraph from gmso.core.atom import Atom @@ -31,6 +32,12 @@ def get_topology_graph(gmso_topology, atomdata_populator=None): """ top_graph = TopologyGraph() atom_index_map = {} + + if len(gmso_topology.sites) == 0: + raise FoyerError( + "Cannot create a topology graph from a topology with no sites." + ) + for j, atom in enumerate(gmso_topology.sites): atom_index_map[id(atom)] = j if isinstance(atom, Atom): diff --git a/gmso/parameterization/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py index f32b07a28..607ee9946 100644 --- a/gmso/parameterization/topology_parameterizer.py +++ b/gmso/parameterization/topology_parameterizer.py @@ -26,7 +26,7 @@ from gmso.parameterization.utils import POTENTIAL_GROUPS -class GMSOParameterizationError(GMSOError): +class ParameterizationError(GMSOError): """Raise when parameterization fails.""" @@ -166,7 +166,7 @@ def _apply_connection_parameters( break if not match and error_on_missing: - raise GMSOParameterizationError( + raise ParameterizationError( f"No parameters found for connection {connection}, group: {group}, " f"identifiers: {connection_identifiers} in the Forcefield." ) @@ -189,14 +189,14 @@ def _verify_forcefields_metadata(self): init_combining_rule = ffs[0].combining_rule for ff in ffs[1:]: if ff.scaling_factors != init_scaling_factors: - raise GMSOParameterizationError( + raise ParameterizationError( "Scaling factors of the provided forcefields do not" "match, please provide forcefields with same scaling" "factors that apply to a Topology" ) if ff.combining_rule != init_combining_rule: - raise GMSOParameterizationError( + raise ParameterizationError( "Combining rules of the provided forcefields do not" "match, please provide forcefields with same scaling" "factors that apply to a Topology" @@ -212,7 +212,7 @@ def run_parameterization(self): """Run parameterization of the topology with give forcefield(s) and configuration.""" scaling_factors, combining_rule = self._verify_forcefields_metadata() if self.topology.is_typed(): - raise GMSOParameterizationError( + raise ParameterizationError( "Cannot parameterize a typed topology. Please provide a topology without any types" ) @@ -222,6 +222,14 @@ def run_parameterization(self): self.topology.identify_connections() if isinstance(self.forcefields, Dict): + if self.topology.n_subtops == 0: + raise ParameterizationError( + f"The provided gmso topology doesn't have any subtopologies." + f"Either use a single forcefield to apply to to whole topology " + f"or provide an appropriate topology whose sub-topology names are " + f"the keys of the `forcefields` dictionary. Provided Forcefields: " + f"{self.forcefields}, Topology: {self.topology}" + ) for subtop in self.topology.subtops: if subtop.name not in self.forcefields: warnings.warn( diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py index dea30e41c..ea0cc0bfb 100644 --- a/gmso/tests/parameterization/parameterization_base_test.py +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -101,7 +101,7 @@ def _assert_same_atom_params(top1, top2): return _assert_same_atom_params - @pytest.fixture(scope="session") + @pytest.fixture def ethane_methane_top(self): cmpd = mb.Compound() cmpd.add(Ethane()) @@ -110,7 +110,7 @@ def ethane_methane_top(self): gmso_top.identify_connections() return gmso_top - @pytest.fixture(scope="session") + @pytest.fixture def ethane_box_with_methane(self): cmpd_box = mb.fill_box([Ethane(), Methane()], [50, 50], density=1.0) return from_mbuild(cmpd_box) diff --git a/gmso/tests/parameterization/test_impropers_parameterization.py b/gmso/tests/parameterization/test_impropers_parameterization.py index ef786a7e9..922b2ebf0 100644 --- a/gmso/tests/parameterization/test_impropers_parameterization.py +++ b/gmso/tests/parameterization/test_impropers_parameterization.py @@ -7,9 +7,7 @@ from gmso.core.views import PotentialFilters from gmso.lib.potential_templates import PotentialTemplateLibrary from gmso.parameterization.parameterize import apply -from gmso.parameterization.topology_parameterizer import ( - GMSOParameterizationError, -) +from gmso.parameterization.topology_parameterizer import ParameterizationError from gmso.tests.parameterization.parameterization_base_test import ( ParameterizationBaseTest, ) @@ -62,7 +60,7 @@ def test_improper_parameterization(self, fake_improper_ff_gmso, ethane): ) def test_improper_assertion_error(self, ethane_methane_top, oplsaa_gmso): - with pytest.raises(GMSOParameterizationError): + with pytest.raises(ParameterizationError): apply(ethane_methane_top, oplsaa_gmso, assert_improper_params=True) @pytest.mark.parametrize( diff --git a/gmso/tests/parameterization/test_parameterization_options.py b/gmso/tests/parameterization/test_parameterization_options.py index 1d529db04..db98145fb 100644 --- a/gmso/tests/parameterization/test_parameterization_options.py +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -2,12 +2,13 @@ import forcefield_utilities as ffutils import pytest +from foyer.exceptions import FoyerError from gmso.core.forcefield import ForceField +from gmso.core.subtopology import SubTopology +from gmso.core.topology import Topology from gmso.parameterization.parameterize import apply -from gmso.parameterization.topology_parameterizer import ( - GMSOParameterizationError, -) +from gmso.parameterization.topology_parameterizer import ParameterizationError from gmso.tests.parameterization.parameterization_base_test import ( ParameterizationBaseTest, ) @@ -30,7 +31,7 @@ def test_parameterization_error_different_scaling_factors( "columbic14Scale": 2.0, } - with pytest.raises(GMSOParameterizationError): + with pytest.raises(ParameterizationError): apply(ethane_methane_top, {"Ethane": ff1, "Methane": ff2}) def test_parameterization_different_combining_rule( @@ -52,7 +53,7 @@ def test_parameterization_different_combining_rule( ff2.combining_rule = "geometric" - with pytest.raises(GMSOParameterizationError): + with pytest.raises(ParameterizationError): apply(ethane_methane_top, {"Ethane": ff1, "Methane": ff2}) def test_different_ffs_apply(self, ethane_methane_top): @@ -63,6 +64,38 @@ def test_different_ffs_apply(self, ethane_methane_top): for key, v in opls.scaling_factors.items(): assert ethane_methane_top.scaling_factors[key] == v + def test_no_subtops_dict_ff(self, oplsaa_gmso): + top = Topology(name="topWithNoSubTops") + with pytest.raises(ParameterizationError): + apply(top, {"subtopA": oplsaa_gmso}) + + def test_missing_subtop_name_ff(self, oplsaa_gmso): + top = Topology(name="top1") + for j in range(0, 10, 2): + top.add_subtopology(SubTopology(name=f"subtop{j+1}")) + with pytest.warns( + UserWarning, + match=r"Subtopology subtop\d will not be parameterized," + r" as the forcefield to parameterize it is missing.", + ): + apply(top, {"subtopA": oplsaa_gmso}) + + def test_diff_combining_rules_error(self, ethane_methane_top): + ff1 = ForceField() + ff1.combining_rule = "lorrentz" + ff2 = ForceField() + ff2.combining_rule = "geometric" + with pytest.raises(ParameterizationError, match=""): + apply(ethane_methane_top, {"Ethane": ff1, "Methane": ff2}) + + def test_empty_ff_foyer_error(self, ethane_methane_top): + with pytest.raises(FoyerError): + apply(ethane_methane_top, ForceField()) + + def test_empty_top_parameterization(self, oplsaa_gmso): + with pytest.raises(FoyerError): + apply(top=Topology(), forcefields=oplsaa_gmso) + def test_isomporhic_speedups(self, ethane_box_with_methane, oplsaa_gmso): ethane_box_with_methane.identify_connections() apply( From f0fd1e506fe8efddbda6d7a1c0b3a4b4c6395caf Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 15 Jun 2022 14:38:55 -0500 Subject: [PATCH 28/32] Docstring minor fix --- gmso/parameterization/parameterize.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/gmso/parameterization/parameterize.py b/gmso/parameterization/parameterize.py index f6b9a673d..70c72f6c1 100644 --- a/gmso/parameterization/parameterize.py +++ b/gmso/parameterization/parameterize.py @@ -24,9 +24,11 @@ def apply( ---------- top: gmso.core.topology.Topology, required The GMSO topology on which to apply forcefields - forcefields: dict, required - The keys are labels that match the subtopology label or site residue_name, and the - values are gmso Forcefield objects that gets applied to the specified molecule + + forcefields: Forcefield or dict, required + The forcefield to apply. If a dictionary is used the keys are labels that match + the subtopology name, and the values are gmso Forcefield objects that gets applied + to the specified subtopology identify_connections: bool, optional, default=False If true, add connections identified using networkx graph matching to match @@ -36,18 +38,23 @@ def apply( identify_connected_components: bool, optional, default=True A flag to determine whether or not to search the topology for repeated disconnected structures, otherwise known as molecules and type each molecule only once. + use_residue_info: bool, optional, default=False A flag to determine whether or not to look at site.residue_name to look parameterize - each molecule only once. Will only be used if identify_connected_components=False + each molecule only once. Currently unused. + assert_bond_params : bool, optional, default=True If True, an error is raised if parameters are not found for all system bonds. + assert_angle_params : bool, optional, default=True If True, an error is raised if parameters are not found for all system angles. + assert_dihedral_params : bool, optional, default=True If True, an error is raised if parameters are not found for all system proper dihedrals. + assert_improper_params : bool, optional, default=False If True, an error is raised if parameters are not found for all system improper dihedrals. From 1e866d7370e479e508bc5df2ede6006fb7027aa8 Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 15 Jun 2022 14:57:42 -0500 Subject: [PATCH 29/32] Remove rel_to_module as is obsolete in forcefield_utilities --- .../tests/parameterization/parameterization_base_test.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py index ea0cc0bfb..a2eac0135 100644 --- a/gmso/tests/parameterization/parameterization_base_test.py +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -17,23 +17,22 @@ def xml_loader(self): @pytest.fixture(scope="session") def oplsaa_gmso(self, xml_loader): - return xml_loader.load("oplsaa", rel_to_module=True).to_gmso_ff() + return xml_loader.load("oplsaa").to_gmso_ff() @pytest.fixture(scope="session") def trappe_ua_gmso(self, xml_loader): - return xml_loader.load("trappe_ua", rel_to_module=True).to_gmso_ff() + return xml_loader.load("trappe_ua").to_gmso_ff() @pytest.fixture(scope="session") def fake_improper_ff_gmso(self, xml_loader): return xml_loader.load( - get_path("fake_ethane_impropers.xml"), rel_to_module=True + get_path("fake_ethane_impropers.xml") ).to_gmso_ff() @pytest.fixture(scope="session") def benzene_alkane_aa_ff_gmso(self, xml_loader): return xml_loader.load( - get_path("benzene_and_alkane_branched_benzene_aa.xml"), - rel_to_module=True, + get_path("benzene_and_alkane_branched_benzene_aa.xml") ).to_gmso_ff() @pytest.fixture(scope="session") From da3a23e66c73a0aacaaa4558b437b05db9e01f2b Mon Sep 17 00:00:00 2001 From: Umesh Timalsina Date: Wed, 15 Jun 2022 15:04:51 -0500 Subject: [PATCH 30/32] Change trappe_ua to trappe-ua for correct loading --- gmso/tests/parameterization/parameterization_base_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py index a2eac0135..8650dbec3 100644 --- a/gmso/tests/parameterization/parameterization_base_test.py +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -21,7 +21,7 @@ def oplsaa_gmso(self, xml_loader): @pytest.fixture(scope="session") def trappe_ua_gmso(self, xml_loader): - return xml_loader.load("trappe_ua").to_gmso_ff() + return xml_loader.load("trappe-ua").to_gmso_ff() @pytest.fixture(scope="session") def fake_improper_ff_gmso(self, xml_loader): From b018a6c6d417e20c56549e36b9a94d5b1befce2f Mon Sep 17 00:00:00 2001 From: Co Quach Date: Wed, 15 Jun 2022 15:08:22 -0500 Subject: [PATCH 31/32] fix typo, add note about specific use case --- gmso/parameterization/parameterize.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/gmso/parameterization/parameterize.py b/gmso/parameterization/parameterize.py index 70c72f6c1..b525dde50 100644 --- a/gmso/parameterization/parameterize.py +++ b/gmso/parameterization/parameterize.py @@ -18,17 +18,20 @@ def apply( assert_dihedral_params=True, assert_improper_params=False, ): - """Set Topology parameter types from GMSO Forcefields. + """Set Topology parameter types from GMSO ForceFields. Parameters ---------- top: gmso.core.topology.Topology, required The GMSO topology on which to apply forcefields - forcefields: Forcefield or dict, required + forcefields: ForceField or dict, required The forcefield to apply. If a dictionary is used the keys are labels that match - the subtopology name, and the values are gmso Forcefield objects that gets applied - to the specified subtopology + the subtopology name, and the values are gmso ForceField objects that gets applied + to the specified subtopology. + Note: if a Topology with no subtopologies is provided, this option will only take + a ForceField object. If a dictionary of ForceFields is provided, this method will + fail. identify_connections: bool, optional, default=False If true, add connections identified using networkx graph matching to match From 318c5ada26e802dea6c6df5b3e2429ecee6d2ae7 Mon Sep 17 00:00:00 2001 From: Co Quach Date: Wed, 15 Jun 2022 18:19:29 -0500 Subject: [PATCH 32/32] pip install forcefield-utilites until new release --- environment-dev.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/environment-dev.yml b/environment-dev.yml index d1c5b8a4b..9763a26b5 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -22,4 +22,5 @@ dependencies: - ipywidgets - ele >= 0.2.0 - pre-commit - - forcefield-utilities + - pip: + - "--editable=git+https://github.com/mosdef-hub/forcefield-utilities.git#egg=forcefield-utilities"