diff --git a/environment-dev.yml b/environment-dev.yml index dd2ed3d7d..9763a26b5 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,5 @@ dependencies: - ipywidgets - ele >= 0.2.0 - pre-commit + - pip: + - "--editable=git+https://github.com/mosdef-hub/forcefield-utilities.git#egg=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/core/topology.py b/gmso/core/topology.py index 3c22d3898..4832bb800 100644 --- a/gmso/core/topology.py +++ b/gmso/core/topology.py @@ -844,7 +844,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() diff --git a/gmso/external/convert_parmed.py b/gmso/external/convert_parmed.py index 54377f128..a8982f40e 100644 --- a/gmso/external/convert_parmed.py +++ b/gmso/external/convert_parmed.py @@ -588,10 +588,12 @@ 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 diff --git a/gmso/parameterization/__init__.py b/gmso/parameterization/__init__.py new file mode 100644 index 000000000..303bb1ca7 --- /dev/null +++ b/gmso/parameterization/__init__.py @@ -0,0 +1,2 @@ +"""GMSO functions that generate parameterized Topologies.""" +from .parameterize import apply diff --git a/gmso/parameterization/foyer_utils.py b/gmso/parameterization/foyer_utils.py new file mode 100644 index 000000000..25bf56575 --- /dev/null +++ b/gmso/parameterization/foyer_utils.py @@ -0,0 +1,114 @@ +"""Utilities for atomtyping a gmso topology with foyer.""" +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 +from gmso.parameterization.subtopology_utils import subtop_bonds + + +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() + 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): + kwargs = ( + atomdata_populator(gmso_topology, atom) + if atomdata_populator + else {} + ) + if atom.name.startswith("_"): + top_graph.add_atom( + name=atom.name, + index=j, # Assumes order is preserved + atomic_number=None, + element=atom.name, + **kwargs, + ) + + else: + top_graph.add_atom( + name=atom.name, + index=j, # Assumes order is preserved + atomic_number=atom.element.atomic_number, + element=atom.element.symbol, + **kwargs, + ) + + for top_bond in gmso_topology.bonds: + atoms_indices = [ + atom_index_map[id(atom)] for atom in top_bond.connection_members + ] + top_graph.add_bond(atoms_indices[0], atoms_indices[1]) + + return top_graph + + +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)) + ) + + +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 + ) + + +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/isomorph.py b/gmso/parameterization/isomorph.py new file mode 100644 index 000000000..811cc429a --- /dev/null +++ b/gmso/parameterization/isomorph.py @@ -0,0 +1,63 @@ +"""TopologyGraph Functions that identify molecules from isomorphism.""" +from collections import deque + +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. + + 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 + """ + graph_queue = deque( + graph.subgraph(c) for c in nx.connected_components(graph) + ) + + graph_of_interest = graph_queue.popleft() + isomorphic_elements = { + graph_of_interest: [], + } + + 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, node_match=top_node_match + ) + 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/parameterization/parameterize.py b/gmso/parameterization/parameterize.py new file mode 100644 index 000000000..b525dde50 --- /dev/null +++ b/gmso/parameterization/parameterize.py @@ -0,0 +1,82 @@ +"""Functions used to atomtype a gmso.Topology.""" +from gmso.parameterization.topology_parameterizer import ( + TopologyParameterizationConfig, + TopologyParameterizer, +) + +__all__ = ["apply"] + + +def apply( + top, + forcefields, + identify_connections=False, + 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: 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. + 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 + 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. + + 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. 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. + """ + 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, + ) + ) + parameterizer = TopologyParameterizer( + topology=top, forcefields=forcefields, config=config + ) + + parameterizer.run_parameterization() + + return parameterizer.topology 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/topology_parameterizer.py b/gmso/parameterization/topology_parameterizer.py new file mode 100644 index 000000000..607ee9946 --- /dev/null +++ b/gmso/parameterization/topology_parameterizer.py @@ -0,0 +1,319 @@ +"""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 ParameterizationError(GMSOError): + """Raise when parameterization fails.""" + + +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.") + + 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", + ) + + 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): + return self.forcefields.get(key) + else: + return self.forcefields + + def _parameterize_sites(self, sites, typemap, ff): + """Parameterize sites with appropriate atom-types from the forcefield.""" + for j, site in enumerate(sites): + site.atom_type = ff.get_potential( + "atom_type", typemap[j]["atomtype"] + ).clone() # Always properly indexed or not? + + 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) + 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, 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 + ): + """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 ParameterizationError( + f"No parameters found for connection {connection}, group: {group}, " + f"identifiers: {connection_identifiers} in the Forcefield." + ) + 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 + ) + + 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 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 ParameterizationError( + "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 ParameterizationError( + "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): + 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( + 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, + self.config.identify_connected_components, + is_subtop=True, + ) + self._parameterize( + subtop, + typemap, + is_subtop=True, # This will be removed from the future iterations + ) + else: + typemap = self._get_atomtypes( + self.get_ff(), + self.topology, + self.config.identify_connected_components, + is_subtop=False, + ) + self._parameterize( + self.topology, + 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 + 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 the 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 + ) + typemap = {} + for graph, mirrors in isomorphic_substructures.items(): + typemap.update( + 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/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", +} diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py index 4c4702993..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 @@ -340,7 +338,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/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/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/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/__init__.py b/gmso/tests/parameterization/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/gmso/tests/parameterization/parameterization_base_test.py b/gmso/tests/parameterization/parameterization_base_test.py new file mode 100644 index 000000000..8650dbec3 --- /dev/null +++ b/gmso/tests/parameterization/parameterization_base_test.py @@ -0,0 +1,115 @@ +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 + + +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").to_gmso_ff() + + @pytest.fixture(scope="session") + def trappe_ua_gmso(self, xml_loader): + 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") + ).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") + ).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 + + @pytest.fixture + 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 + 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..922b2ebf0 --- /dev/null +++ b/gmso/tests/parameterization/test_impropers_parameterization.py @@ -0,0 +1,104 @@ +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 ParameterizationError +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(ParameterizationError): + 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_opls_gmso.py b/gmso/tests/parameterization/test_opls_gmso.py new file mode 100644 index 000000000..c5ed7639f --- /dev/null +++ b/gmso/tests/parameterization/test_opls_gmso.py @@ -0,0 +1,56 @@ +import glob +from pathlib import Path + +import parmed as pmd +import pytest +from pkg_resources import resource_filename + +from gmso.external.convert_parmed import from_parmed +from gmso.parameterization.parameterize import apply +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) + + +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(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, + 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 = 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) + + 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 new file mode 100644 index 000000000..db98145fb --- /dev/null +++ b/gmso/tests/parameterization/test_parameterization_options.py @@ -0,0 +1,130 @@ +import random + +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 ParameterizationError +from gmso.tests.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) + + +class TestParameterizationOptions(ParameterizationBaseTest): + 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(ParameterizationError): + 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(ParameterizationError): + 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_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( + 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 diff --git a/gmso/tests/parameterization/test_subtopology_utils.py b/gmso/tests/parameterization/test_subtopology_utils.py new file mode 100644 index 000000000..32c4b7690 --- /dev/null +++ b/gmso/tests/parameterization/test_subtopology_utils.py @@ -0,0 +1,112 @@ +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.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) +from gmso.utils.connectivity import identify_connections + + +class TestSubTopologyUtils(ParameterizationBaseTest): + @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_gmso.py b/gmso/tests/parameterization/test_trappe_gmso.py new file mode 100644 index 000000000..84e489e40 --- /dev/null +++ b/gmso/tests/parameterization/test_trappe_gmso.py @@ -0,0 +1,58 @@ +import glob +from pathlib import Path + +import pytest +from pkg_resources import resource_filename + +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.parameterization.parameterization_base_test import ( + ParameterizationBaseTest, +) + + +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(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_foyer, + trappe_ua_gmso, + 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_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) + + 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" + ) + assert_same_connection_params( + gmso_top, gmso_top_from_parmeterized_pmd, "dihedrals" + )